.. tip:: All input files can be downloaded: :download:`Files `. .. tip:: For more information of this section, please refer to these pages: - :doc:`../keywords/opt` - :doc:`../keywords/mecp` Transition State Search ========================================= .. contents:: :local: This tutorial will lead you step by step to search transition states using Qbics. There are two kinds of methods to search transition states in Qbics: - **The QST2/dimer method and the nudged elastic band (NEB) method.** For both methods, a reactant and a product geometry have to be provided. - **The diabatic minimum energy crossing point (dMECP) method.** This method is highly useful for AB+C=A+BC type reactions. Here, only the reactant geometry is needed. .. hint:: In most cases, searching transition states is not a trivial task. You may have to try several starting structures and combine NEB and dimer methods to find the transition state. The following examples are just for demonstration purposes, and may not be suitable for your own calculations. Example: NEB and QST2/Dimer for Amine-Catalyzed Aldol Reaction ----------------------------------------------------------------------------------- In this section, we will search the transition state of the following amine-catalyzed aldol reaction: .. figure:: figs/a20.jpg The first step is to prepare the input file for the reactant and product geometries. They are given in here: .. tabs:: .. tab:: r1.xyz .. code-block:: :caption: r1.xyz :linenos: 18 Reactant N -0.44812625 -0.05503381 -0.29939653 C 0.28940454 -1.28430107 0.02594994 C 0.42365103 0.85305301 -1.05852939 H 1.14858130 -1.03947388 0.61482595 H -0.34515443 -1.94528943 0.57851603 H 0.60162891 -1.76291266 -0.87867534 C -0.11218661 1.70819672 -1.96305927 H 1.48128134 0.83886150 -0.89692269 H 0.52237224 2.36918599 -2.51562444 H -1.16981687 1.72238788 -2.12466634 C -1.82051925 1.96886205 -0.83651451 H -1.13776031 2.16858444 -0.03723284 C -2.45646906 3.13519149 -1.61551609 O -2.10383295 0.78091580 -1.13997394 H -1.25109519 -0.28384426 -0.84974794 H -3.37275441 3.42197565 -1.14320685 H -1.78384110 3.96729827 -1.62387764 H -2.65467293 2.82667229 -2.62071813 .. tab:: r2.xyz .. code-block:: :caption: r2.xyz :linenos: 18 Product N -0.22474732025333 -0.31730283410046 -0.51491194374209 C 0.54138841877455 -1.42626714378869 0.01659619702235 C 0.37581142979957 0.57244190213102 -1.18506998859001 H 1.59956963538446 -1.37223656022845 -0.26069633867705 H 0.44285772700755 -1.41562576183705 1.10277773062867 H 0.10612456734704 -2.35145036624173 -0.36248430677077 C -0.37357663580959 1.73548805223194 -1.76480799540763 H 1.45524939783346 0.53941227186085 -1.38723095744134 H 0.27481589890769 2.61572199377069 -1.78645396480960 H -0.63716949392230 1.49356209011140 -2.80064855509307 C -1.66735730523003 2.04995476069694 -0.98212699921514 H -1.39146640788106 2.28056626272631 0.06585987026323 C -2.39173812629474 3.24996769168947 -1.58779580073698 O -2.54492763771693 0.95029342074102 -1.01960271506481 H -2.03040745192793 0.17070983855352 -0.71586070424615 H -3.30862690757118 3.43506393231313 -1.02954259805264 H -1.76246901303224 4.13842282177014 -1.54692923821636 H -2.65264077541507 3.04160762759991 -2.62544169185061 .. tabs:: .. tab:: r1.xyz .. figure:: figs/a21.jpg .. tab:: r2.xyz .. figure:: figs/a22.jpg Both structures can be constructed by chemical ituition and have been optimized. However, even with these two semmingly reasonable structures, it is not easy to search the transition state. So, we will combine NEB and dimer, semiempirical and DFT methods to search the transition state. **Step 1: NEB Search.** NEB is a method to search the path by connecting the reactant and product geometries. It can be a useful starting point for the dimer method. To rapidly get this, we will use xTB method: .. code-block:: :caption: anion-neb.inp :linenos: mol r1.xyz end mol2 r2.xyz end basis def2-svp end opt type neb num_images 15 neb_k 0.01 end task opt xtb end Let's examine the input file: - The reactant and product geometries are defined by ``mol`` and ``mol2``. - The NEB method is activated by ``opt...end`` with ``type neb``. The default option of ``type`` is ``min`` which is for geometry optimization, so we have to set ``type neb`` explicitly. - ``num_images`` is the number of images in the NEB path, which can be adjusted according to your needs. A good choice is between ``10`` and ``20``. - ``neb_k`` is the spring constant for the NEB path, which can also be adjusted according to your needs. A good choice is between ``0.01``. After running the input file, we can get the following output file: - ``anion-neb-opt.xyz`` The "transition state" structure. But in many cases, we suggest you to visualize it and check whether it is reasonable. - ``anion-neb-opt-traj.xyz`` The NEB trajectory, which can be visualized by `Qbics-MolStar `_: .. figure:: figs/a23.gif Indeed, this trajectory is our target reaction. We can choose the two closest images on the path that the trasition state may lie between them. In this case, we can choose the 5th and 7th images in ``anion-neb-opt-traj.xyz``, which are saved in ``a5.xyz`` and ``a7.xyz``: .. tabs:: .. tab:: a5.xyz .. code-block:: :caption: a5.xyz :linenos: 18 element x y z N -0.03436 -0.42780 -1.21081 C 0.19877 -1.31325 -0.09641 C 0.54963 0.79676 -1.27792 H 1.21009 -1.16425 0.27978 H -0.51219 -1.13914 0.71792 H 0.09714 -2.34661 -0.43350 C 0.12282 1.79024 -2.06831 H 1.40495 0.93556 -0.62766 H 0.65351 2.72274 -2.12058 H -0.66378 1.63364 -2.78419 C -1.86900 2.10882 -0.51273 H -1.16547 2.64884 0.14563 C -2.61217 2.99713 -1.46905 O -2.13799 0.94532 -0.31815 H -0.98415 -0.45631 -1.55950 H -3.41739 3.47081 -0.90569 H -1.97325 3.77562 -1.87311 H -3.05645 2.41221 -2.27009 .. tab:: a7.xyz .. code-block:: :caption: a7.xyz :linenos: 18 element x y z N -0.04991 -0.40765 -1.08610 C 0.35194 -1.39200 -0.13587 C 0.50298 0.72490 -1.16520 H 1.36153 -1.23389 0.25614 H -0.36195 -1.32388 0.69157 H 0.25447 -2.38320 -0.57375 C -0.26333 1.81516 -1.81912 H 1.42373 0.97712 -0.62468 H 0.32790 2.72274 -1.93602 H -0.65750 1.50523 -2.78663 C -1.50488 2.07052 -0.89479 H -1.14590 2.52877 0.04533 C -2.47962 3.02762 -1.58061 O -2.17108 0.87844 -0.60177 H -1.60693 0.12971 -0.93634 H -3.25074 3.33649 -0.87990 H -1.95994 3.90581 -1.95044 H -2.96008 2.50844 -2.40618 Now we will use the dimer/QST2 algorithm to search the transition state between these two images. Below, we show how to do this on B3LYP/def2-SVP levels of theory: .. code-block:: :caption: anion-dft-qst2.inp :linenos: mol a5.xyz end mol2 a7.xyz end basis def2-svp end scf charge 0 spin2p1 1 end opt type qst2 # You can also use dimer. end task opt b3lyp end Fortunately, we can get the transition state structure, which is saved in ``anion-dft-qst2-opt.xyz`` and is shown below: .. figure:: figs/a24.jpg Example: dMECP for F\ :sup:`-`\ +CH\ :sub:`3`\ Cl = CH\ :sub:`3`\ F+Cl\ :sup:`-` ----------------------------------------------------------------------------------- For AB+C=A+BC type reactions, we can use the dMECP method to search the transition state. For example, for the reaction F\ :sup:`-`\ +CH\ :sub:`3`\ Cl = CH\ :sub:`3`\ F+Cl\ :sup:`-`, we can use the following input file: .. code-block:: :caption: fch3cl.inp :linenos: basis def2-svp end mol C -0.43654823 1.13197968 0.00000000 H -0.07987539 1.63637787 0.87365150 H -0.07987539 1.63637787 -0.87365150 H -1.50654823 1.13199286 0.00000000 Cl 0.15009830 -0.52737135 0.00000000 F -1.63124586 2.75454539 -0.44024554 end scf type u # This must be set. charge -1 spin2p1 1 end mecp num_steps 200 energy_cov 1.E-5 # Reactant. frag1 0 1 1-5 frag1 -1 1 6 # Product. frag2 0 1 1-4 6 frag2 -1 1 5 end task mecp b3lyp end We can see that we want to search the transition state at B3LYP/def2-SVP level of theory. The dMECP method is activated by ``mecp`` in ``task...end``, and is controlled by ``mecp...end`` option. Below are some important points: - dMECP is implemented with the TSO method, so in ``scf...end``, we have to set ``type u``. - ``num_steps`` is the number of steps for the dMECP search, which can be adjusted according to your needs. - The reactant and product connectivity are defined by ``frag1`` and ``frag2``, respectively. Their grammar is similar to ``frag`` in :doc:`../keywords/eda` and :doc:`../keywords/scfguess`, but here we have to define the connectivity of the reactant and product separately by ``frag1`` and ``frag2``. The geometry structure in ``fch3cl.inp`` is shown below: .. figure:: figs/a16.jpg The ``frag1`` and ``frag2`` definitions are shown below: .. figure:: figs/a17.jpg Therefore, we just need to **define** the "connectivity" and "charge" of the reactant and product, **no matter what the geometry is.** The dMECP method will automatically optimize the geometry to find the transition state. After running the input file, we can get the following output file: - ``fch3cl-mecp.xyz`` The final transition state. - ``fch3cl-mecp-traj.xyz``. The dMECP trajectory. Both files can be visualized by `Qbics-MolStar `_: .. figure:: figs/a18.jpg .. figure:: figs/a19.gif In ``fch3cl.out``, we can find following content: .. code-block:: :caption: fch3cl.out :linenos: MECP search step: 1 Reference state energy: -599.67127343 Hartree State 1 energy: -599.62533338 Hartree State 2 energy: -599.49026107 Hartree |E(2)-E(1)|: 0.13507232 Hartree RMS gradient: 2.14932461 Hartree/Bohr Here, ``Reference state energy`` is the energy of standard DFT (**adiabatic energy**), ``State 1 energy`` and ``State 2 energy`` is the **diabatic** energy of ``frag1`` and ``frag2``. For the transition state, the two diabatic energies should be identical: .. figure:: figs/a19.jpg