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.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

r2.xyz

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

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.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

a7.xyz

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:

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
    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 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.

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

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:

Transition State Search

Example: NEB and QST2/Dimer for Amine-Catalyzed Aldol Reaction

Example: dMECP for F-+CH3Cl = CH3F+Cl-

Example: NEB and QST2/Dimer for Amine-Catalyzed Aldol Reaction 

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