Tip

All input files can be downloaded: Files.

Tip

For more information of this section, please refer to these pages:

• opt

• mecp

Transition State Search

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:

The first step is to prepare the input file for the reactant and product geometries. They are given in here:

r1.xyzr2.xyz

r1.xyz

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

r1.xyzr2.xyz

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:

anion-neb.inp

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:

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:

a5.xyza7.xyz

a5.xyz

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

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:

anion-dft-qst2.inp

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:

Example: dMECP for F^-+CH₃Cl = CH₃F+Cl^-

For AB+C=A+BC type reactions, we can use the dMECP method to search the transition state. For example, for the reaction F^-+CH₃Cl = CH₃F+Cl^-, we can use the following input file:

fch3cl.inp

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
    dmecp 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 dmecp 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 eda and 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:

The frag1 and frag2 definitions are shown below:

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. dMECP only needs reactant complex struture. This is very convenient for the transition state search.

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:

In fch3cl.out, we can find following content:

fch3cl.out

Adiabatic state energy:        -599.69936976 Hartree
Diabatic state 1 energy:       -599.62911003 Hartree
Diabatic state 2 energy:       -599.62485338 Hartree
|E(2)-E(1)|:                      0.00425665 Hartree
delta|E(2)-E(1)|:                 0.00000520 Hartree
RMS gradient:                     0.00008376 Hartree/Bohr
RMS displacement:                 0.00098426 Angstrom

Here, Adiabatic state energy is the energy of standard DFT. Diabatic state 1 energy and Diabatic state 2 energy is the diabatic energy of frag1 and frag2.