.. _reactivity: Reactivity workflow =================== The tendency of substances to undergo a chemical reaction is determined by the activation and reaction energies of the reaction. The activation energy is the energy difference between the reactants and the transition state, while the reaction energy is the energy difference between the reactants and the products. To conveniently study reactions, you can use the Reactivity Workflow, which automatically calculates the energies of the reactants, transition states, and products, and the corresponding activation and reaction energies. The workflow uses the `Nudged Elastic Band (NEB) `_ technique, which optimizes intermediate structures along a path connecting the reactants and products of a chemical reaction. NEB calculations are commonly used to determine the activation energies and transition states of chemical reactions. In the following, some examples are given where the Reactivity Workflow has been used. The basics of the workflow are explained using the example of the :ref:`rotational barrier of ethane`, followed by the :ref:`Diels-Alder mechanism` and the study of :ref:`rotational groups`. .. _rotbarrier: Rotational barrier ------------------ We have a look at the rotational barrier of ethane, from the staggered configuration via the eclipsed to the staggered configuration. .. figure:: reactivity/reaction.png :align: center :width: 500 Rotation around the central bond of ethane. How to run the calculation ........................... In principle, the xyz coordinates of ethane introduced in :ref:`Scans` can be used. However, the reactivity workflow requires both reactant as well as product structure as input. The XYZ files of the :download:`reactant ` (``ethane_reac.xyz``) and the :download:`product ` (``ethane_prod.xyz``) ethane can be downloaded to the directory in which you want to run the calculation. .. important:: In case you create the XYZ files yourself: In order for the NEB-TS to work, the atom numbering in both reactant and product geometries must be the same. Which means that a given carbon on the first, should have the same order on the xyz table in the last. The best way to guarantee this is to start from one geometry and move the atoms to draw the other. The workflow can be run with the following command, specifying the product structure with the ``-product`` keyword. .. prompt:: bash $ weasel ethane_reac.xyz -W Reactivity -product ethane_prod.xyz Steps of the reactivity workflow ................................. The workflow optimizes the reactant and the product geometry and it finds the TS that connects both. The following calculation steps are the basis of the workflow. .. figure:: reactivity/workflow_reactivity.jpg :align: center :width: 650 The frequency calculation is required to * check for the (minimum and TS) nature of the optimized structures and * the thermochemical correction using the RRHO approximation. The final SP energy calculation recomputes all three structures at a more accurate DFT level, here based on the basic workflow with wB97X-V/def2-TZVP. .. admonition:: Recommendation :class: tip Highly accurate activation and reaction free energies can be critical for reactivity calculations. To obtain a more accurate result, it is recommended to recalculate the final SP energy with a SP wavefunction based method, e.g. DLPNO-CCSD(T) by adding the keyword ``-spwf``. .. note:: In addition to the ``Reactivity`` workflow, there is also the ``Reactivity-Explore`` workflow. They differ in that ``Reactivity`` uses the standard methods for pre-optimization (XTB), optimization and frequency calculation (r2Scan-3c), while ``Reactivity-Explore`` uses XTB by default for all three steps. This makes the latter faster but less accurate. The methods can also be changed via the respective keyword input. .. note:: The solvent in the ``Reactivity`` standard workflow is water and needs to be adapted if necessary. Structure of output files .......................... The calculation produces the following files:: . ├── ethane_reac.xyz ├── ethane_prod.xyz └── ethane ethane_reac_ethane_prod_exploreReactivity ├── ethane_reac_ethane_prod_Reactivity.input.xyz ├── ethane_reac_ethane_prod_Reactivity_MEP_trj.xyz ├── ethane_reac_ethane_prod_Reactivity_prod_Opt.xyz ├── ethane_reac_ethane_prod_Reactivity_Product.input.xyz ├── ethane_reac_ethane_prod_Reactivity_reac_Opt.xyz ├── ethane_reac_ethane_prod_Reactivity.report ├── ethane_reac_ethane_prod_Reactivity.restart.tsv.bz2 ├── ethane_reac_ethane_prod_Reactivity.summary ├── ethane_reac_ethane_prod_Reactivity_TS_Opt_Freq.hess.v006.xyz ├── ethane_reac_ethane_prod_Reactivity_TS_Opt.xyz ├── BuildTopo │   └── ORCA job files ├── PreOpt │   └── ORCA job files ├── Opt │   └── ORCA job files ├── NEB │   └── ORCA job files ├── SP_DFT    └── ORCA job files Below is a table that outlines the most important files created and their respective descriptions. +--------------------------------+----------------------------------+ | File | Description | +================================+==================================+ | stem_reac_Opt.xyz | optimized reactant structure | +--------------------------------+----------------------------------+ | stem_prod_Opt.xyz | optimized product structure | +--------------------------------+----------------------------------+ | stem_TS_Opt.xyz | optimized TS structure | +--------------------------------+----------------------------------+ | stem_TS_Opt_Freq.hess.v006.xyz | Trajectory file with the TS mode.| | | Can be visualized with avogadro | +--------------------------------+----------------------------------+ | stem_NEBTS_MEP_trj.xyz | Minimum energy path of the | | | converged NEB calculation | +--------------------------------+----------------------------------+ | stem_exploreReactivity.report | Output file of the WEASEL run | +--------------------------------+----------------------------------+ | stem_exploreReactivity.summary | Summary file for the WEASEL run | +--------------------------------+----------------------------------+ .. figure:: reactivity/reactionReTsPr.png :align: center :width: 500 Reactant (left), TS (middle) and product (right) structure. .. note:: Further information on the minimum-energy path is found in the NEB task directory. .. _reactivitydEdG: Reaction and activation free energies ..................................... The summary file stores different sets of data, but we'll focus on the most important ones here. At its very end, the activation and reaction energy and free energy at DFT level:: -0.000233 Reaction Energy dE TD RIJCOSX-B3LYP def2-TZVP(-F) CPCM(Water) 0 1 NEB-TS Summary DFT 2.571387 Activation Energy dE TS RIJCOSX-B3LYP def2-TZVP(-F) CPCM(Water) 0 1 NEB-TS Summary DFT -0.007393 Reaction Free Energy dG TD RIJCOSX-B3LYP def2-TZVP(-F) CPCM(Water) 0 1 NEB-TS Summary DFT 2.537263 Activation Free Energy dG TS RIJCOSX-B3LYP def2-TZVP(-F) CPCM(Water) 0 1 NEB-TS Summary DFT as well as at wavefunction level:: 0.000257 Reaction Energy dE TD DLPNO-CCSD(T) def2-TZVP CPCM(Water) 0 1 NEB-TS Summary WF 2.818516 Activation Energy dE TS DLPNO-CCSD(T) def2-TZVP CPCM(Water) 0 1 NEB-TS Summary WF -0.006903 Reaction Free Energy dG TD DLPNO-CCSD(T) def2-TZVP CPCM(Water) 0 1 NEB-TS Summary WF 2.784392 Activation Free Energy dG TS DLPNO-CCSD(T) def2-TZVP CPCM(Water) 0 1 NEB-TS Summary WF are provided as key results of the simulation. .. _dielsalder: Diels-Alder reaction -------------------- More interesting is the formation of new bonds, as in the following example of the Diels-Alder reaction of cyclopentadiene. .. figure:: reactivity/DA_reaction.png :align: center :width: 450 Diels-Alder reaction of cyclopentadiene. In the following, we will compare the Diels-Alder reaction between the two cyclopentadienes with a reaction between a cyclopentadiene and a methylcyclopentadiene. The purpose of this comparison is to demonstrate that the reactivity workflow in WEASEL can be used for reactions involving bond formation, and to illustrate the potential for transferring results from one WEASEL calculation to another, thereby streamlining and accelerating the process. Let's start with the example reaction between the two cyclopentadienes. Reaction between two cyclopentadienes ........................................ If you want to do the calculation between the two cyclopentadienes yourself, you can download the XYZ file ``DA_reac.xyz`` of the :download:`reactant ` complex and the file ``DA_prod.xyz`` for the :download:`product ` and save them in the directory where you want to do the calculation. If you create the structures yourself, make sure that the atoms have the same order in both the reactant and product structures. The file ``DA_reac.xyz`` already contains the two molecules that will react, positioned close together to form the reactant complex. Instead of manually searching for a suitable reactant complex with multiple interacting molecules, as in this Diels-Alder example, you can use the Docker in WEASEL to generate the required complex. The Docker tutorial can be found :ref:`here`. To perform the calculation, you can use the reactivity workflow, as already applied to the rotational barrier of ethane, by running the following command: .. prompt:: bash $ weasel DA_reac.xyz -W Reactivity -product DA_prod.xyz .. admonition:: Recommendation :class: tip In this scenario, the reactant complex has many degrees of freedom. Consequently, during structure optimization, there is a chance that points on the potential energy surface with multiple negative frequencies will be mistaken for a local minimum. To reduce the likelihood of this happening in your calculation, you can use the ``-correctimagmodes`` keyword. This will automatically repeat the optimization if any (unwanted) imaginary modes are detected during the frequency analysis. The calculation produces the minimum energy path shown below: .. figure:: reactivity/DA_trj.gif :align: center :width: 450 Minimum energy path for the Diels-Alder reaction from the NEB-TS calculation. The resulting reaction and activation (free) energies are printed again in the summary file as following:: -24.523926 Reaction Energy dE TD wB97X-V def2-TZVP CPCM(Water) 0 1 NEB-TS Summary DFT 21.780152 Activation Energy dE TS wB97X-V def2-TZVP CPCM(Water) 0 1 NEB-TS Summary DFT -17.070249 Reaction Free Energy dG TD wB97X-V def2-TZVP CPCM(Water) 0 1 NEB-TS Summary DFT 24.855281 Activation Free Energy dG TS wB97X-V def2-TZVP CPCM(Water) 0 1 NEB-TS Summary DFT Reaction between cyclopentadiene and methylcyclopentadiene - Using a template Hessian ....................................................................................... Now that we have the results of the Diels-Alder reaction between the two cyclopentadienes, we can proceed to see how the reactivity changes when we add a methyl group to one of the cyclopentadienes. Instead of starting the whole process from scratch, we can use the optimized reactant complex, product, and transition state from the previous calculation and modify them accordingly. The modification involves replacing one of the hydrogen atoms attached to one of the carbon atoms with a methyl group. This can be done in any molecular builder tool of your choice. If you are not familiar with it, you can refer to our tutorial on how to use Avogadro. When adding the methyl group, be sure to replace the exact hydrogen in the reactant complex, product, and transition state. This step will keep the numbering consistent in all three systems. .. figure:: reactivity/DA2_systems.png :align: center :width: 600 Input structures for Diels-Alder reaction between cyclopentadiene and methylcyclopentadiene: reactants (left), transition state (center), and product (right). You can save the updated structures in the files ``DA_reactant2.xyz``, ``DA_product2.xyz`` and ``DA_ts2.xyz``. From here, there are two approaches to calculating the activation and reaction (free) energies for this second reaction. One is to follow the same steps as before, using the reactivity workflow. In this case, however, you can use the previously prepared transition state structure as a guess for the transition state using the ``-tsguess`` keyword. The WEASEL command would then read:: $ weasel DA_reactant2.xyz -W Reactivity -product DA_product2.xyz -tsguess DA_ts2.xyz Running this calculation would produce the following results in the summary file:: -24.502095 Reaction Energy dE TD wB97X-V def2-TZVP CPCM(Water) 0 1 NEB-TS Summary DFT 24.679614 Activation Energy dE TS wB97X-V def2-TZVP CPCM(Water) 0 1 NEB-TS Summary DFT -17.130277 Reaction Free Energy dG TD wB97X-V def2-TZVP CPCM(Water) 0 1 NEB-TS Summary DFT 27.645757 Activation Free Energy dG TS wB97X-V def2-TZVP CPCM(Water) 0 1 NEB-TS Summary DFT Comparing the results of the Diels-Alder reaction of two cyclopentadienes and a cyclopentadiene with a methylcyclopentadiene, it can be seen that both have a similar reaction (free) energy. However, the activation (free) energy is higher in the case of the cyclopentadiene with methylcyclopentadiene, showing that the methyl group inhibits the reaction. While this approach is convenient in that WEASEL would provide you directly with the activation and reaction (free) energies, it is not the most computationally efficient way. The reactivity workflow, despite the transition state guess, would perform a full NEB calculation, which is the most time-consuming step, and the transition state guess only makes the search more accurate, not necessarily significantly faster. However, due to the previous calculation of the two cyclopentadienes, you already know your transition state quite accurately and therefore do not really need the NEB calculation. This brings us to the second approach: To skip the NEB calculation when calculating your activation and reaction (free) energies, you can do three separate calculations. You need to optimize the reactants and product structures using the following WEASEL commands:: $ weasel DA_reactant2.xyz -freq and:: $ weasel DA_product2.xyz -freq The ``-freq`` is needed to calculate the free energies. .. _templateHess: To optimize the transition state you must use the ``-optts`` keyword. To further increase the efficiency of the calculation you can use the Hessian from the previous calculation (``DA_ts.hess``) as a template for the current run using the keyword ``-optts-inputhessian``. This template Hessian strategy can make the process much faster and more accurate. The complete WEASEL command to compute the transition state would then be:: $ weasel DA_ts2.xyz -optts -optts-inputhessian DA_ts.hess -freq From the separate results of the reactant, product, and transition state optimizations, you can then manually calculate the reaction and activation (free) energies. The reaction (free) energy must be calculated as the difference between the total (free) energy of the products and the reactants, and the activation (free) energy must be calculated as the difference between the total (free) energy of the transition state and the reactants. .. note:: To further simplify the transition state search, you can also use the keywords ``-optts-activeatoms A [A ...]`` or ``-optts-activeatoms-template A [A ...]``. In the first case you have to specify the atom numbers of the active atoms of the transition state, i.e. the atoms involved in the reaction. In the second case you have to specify the active atoms of the transition state of the template Hessian (``-optts-inputhessian``). .. _rotgroup: Rotating groups --------------- For reactions, in which larger chemical groups are rotating, or for the comparison of multiple paths (e.g. rotation via two different directions), providing a guess for the TS structure is recommended. .. figure:: reactivity/IrRotation.png :align: center :width: 500 Rotation of a large group around the Ir-C bond. This TS guess can be a simple "educated guess", and will be used for the generation of the initial path. It can be provided with the keyword ``-tsguess``" .. prompt:: bash $ weasel Ir_HCl_II.xyz -W exploreReactivity -product Ir_HCl_III.xyz -tsguess Ir_TSGuessI.xyz The exploreReactivity workflow produces the reaction's MEP: .. figure:: reactivity/IrComplexMEP.gif :align: center :width: 450 Minimum energy path from the NEB calculation. TS structure, and activation and reaction free energies:: -10.635936 Reaction Energy dE TD RIJCOSX-B3LYP def2-TZVP(-F) Gas phase 0 1 NEB-TS Summary DFT 16.795442 Activation Energy dE TS RIJCOSX-B3LYP def2-TZVP(-F) Gas phase 0 1 NEB-TS Summary DFT -9.908759 Reaction Free Energy dG TD RIJCOSX-B3LYP def2-TZVP(-F) Gas phase 0 1 NEB-TS Summary DFT 19.399550 Activation Free Energy dG TS RIJCOSX-B3LYP def2-TZVP(-F) Gas phase 0 1 NEB-TS Summary DFT Keywords and Remarks --------------------- **Keywords for reactivity calculations** .. list-table:: :widths: 30 70 :header-rows: 1 :class: fixed-width-table * - Keyword - Description * - ``-W Reactivity`` - | Calculates the reaction and activation free energy for a reaction between | provided reactant and product. * - ``-W Reactivity-Explore`` - | Calculates the reaction and activation free energy for a reaction between | provided reactant and product on a lower level of theory. * - ``-product PRODUCT`` - | File with product structure. * - ``-reactivity-nebts`` - | Find TS with NEB-TS Feature. * - ``-neb-usemepguess`` - | Use guess for minimum energy path in NEB calculation. Guess provided in | ORCA's allxyz file format as input structure file. * - ``-neb-fixends`` - | Reactant and product structure should be kept as they are. * - ``-neb-sp-fullreaction`` - | If requested, compute SP and frequency for reactant, TS, and product. * - ``-neb-boost`` - | Use Fast-NEB-TS to boost the convergence. Comes with slightly lower | robustness. * - ``-neb-nimages NEB_NIMAGES`` - | Use this number of images for the NEB calculation. * - ``-correctimagmodes`` - | Rerun optimization in case (unwanted) imaginary modes are found with | the frequency analysis. * - ``-no-correctimagmodes`` - | Do not rerun optimization in case (unwanted) imaginary modes are found | with the frequency analysis. * - ``-reactivity-en`` - | Compute electro- and nucleophilicity. **Keywords for transition state calculations** .. list-table:: :widths: 30 70 :header-rows: 1 :class: fixed-width-table * - Keyword - Description * - ``-optts`` - | Do TS optimization * - ``-optts-inputhessian INPUTHESSIAN`` - | Input hessian of (previously optimized) TS structure. Needs to be an | ORCA Hessian file. * - ``-optts-activeatoms A [A ...]`` - | Active atoms of transition state, i.e. atoms that are involved in reaction. * - ``-optts-activeatoms-template A [A ...]`` - | Active atoms of transition state of template Hessian (specified via | ``-optts-inputhessian``), i.e. atoms that are involved in reaction. * - ``-scants`` - | Do relaxed surface scan and subsequent TS optimization. * - ``-reactivity-nebts`` - | Find TS with NEB-TS Feature. * - ``-neb-sp-onlyts`` - | If requested, compute SP and frequency only for TS structure. * - ``-tsguess TSGUESS`` - | File with TS guess structure for finding TS (optional). * - ``-freq`` - | Compute vibrational frequencies, e.g. at the end of a (TS) optimization. | Free energies are automatically calculated. * - ``-no-freq`` - | Do not compute vibrational frequencies, e.g. at the end of a (TS) | optimization.