Finding transition states with NEB-TS#

Suppose you wanted to predict the transition state (TS) for the hydrolysis of methyl-acetate into acetic acid and methanol:

../_images/react_2d.png — Figure: Hydrolysis reaction of methyl-acetate.#

Maybe you need to investigate the mechanism, or even to compute the \(\Delta G^\ddagger\) and try to predict reaction rates. How could one do that using theory?

In ORCA there is a black box method called NEB-TS (from Nudged Elastic Band with TS optimization), that can find the TS structure only from the geometries of the reactants and products!

Defining your reactants and products#

So the first step is to optimize the "reactants" and "products".

Important

Please, have in mind that reactants and products means everything that is left and right of the arrow. Accordingly, the encounter complex of the reacting molecules has to be provided in one single structure file. The same is required for the products.

Hint

Reasonably pre-oriented encounter complexes significantly increase the probability of a successful TS search.

In this example we will use Grimme's efficient PBEh-3c composite method. The corresponding input looks like:

!PBEh-3c OPT

../_images/reac-prod.png — Figure: Optimized reactants and products.#

Important

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 reactant, should have the same order in the xyz table in the product. The best way to guarantee this is to start from one geometry, make a copy and move the atoms to draw the next structure.

Note

If using Avogadro, you have to disable the Adjust Hydrogen feature, otherwise atoms will be automatically placed when a bond is broken, and not necessarily where you want them to be.

The NEB-TS input#

With the optimized reactants and products at hand, we can create the NEB-TS input which is quite simple. NEB-TS is requested via the NEB-TS keyword and all other settings are defined in the %NEB block. To verify that our probably found TS is reasonable we also request a final vibrational frequencies calculation.

!PBEh-3c NEB-TS FREQ

%NEB 
 NEB_END_XYZFILE "product.xyz" 
END

*XYZFILE 0 1 reactant.xyz

In case you have a guess structure for a potential TS as well, that can be included as:

!PBEh-3c NEB-TS FREQ

%NEB 
 NEB_END_XYZFILE "product.xyz"
 NEB_TS_XYZFILE "guessTS.xyz"
END

* XYZFILE 0 1 reactant.xyz

Note

The reactant and products will not be optimized during the default run, so make sure you are using the same method you chose for the previous optimization. If you want to reoptimize them automatically, set PREOPT_ENDS TRUE in the %NEB block.

Running the NEB-TS#

With our input being set up, we can start the NEB-TS calculation. The first output contains some details of the method.

----------------------
NEB settings
----------------------
Method type                             ....  climbing image
Threshold for climbing image            ....  2.00e-02 Eh/Bohr
Free endpoints                          ....  off
Tangent type                            ....  improved
Number of intermediate images           ....  8
Number of images free to move           ....  8
Spring type for image distribution      ....  distance between adjacent images
Spring constant                         ....  energy weighted (0.0100 -to- 0.1000) Eh/Bohr^2
Spring force perp. to the path          ....  none
Generation of initial path              ....  image dependent pair potential
Initial path via TS guess               ....  off
[...]

Then, comes the construction of a initial path:

Generation of  the initial path:
Two fragments detected in initial structure. Preparing initial structure... Done
Minimize RMSD between reactant and product configurations     .... done
       RMSD before minimization         ....  2.1833 Bohr
       RMSD after minimization          ....  1.7369 Bohr
Performing linear interpolation         ....  done
Fragment preparation                           ....  true
Bond cutoff factor                             ....  1.200e+00
Max. fragment distance                         ....  3.500e+00 Ang.
Interpolation using                            ....  IDPP (Image Dependent Pair Potential)
IDPP-Settings:
       Initial path for IDPP optimization      ....  Linear interpolation
       Pairwise distance interpolation         ....  Linear
       Remove global transl. and rot.          ....  true
       Convergence tolerance                   ....  0.0100 1/Ang.^3
       Max. number of iterations               ....    7000
       Spring constant                         ....  1.0000 1/Ang.^4
       Time step                               ....  0.0100 fs
       Max. movement per iteration             ....  0.0500 Ang.
       Full print                              ....  false
idpp initial path generation successfully converged in 130 iterations
Displacement along initial path:             12.1839 Bohr
Writing initial trajectory to file      ....  nebts_initial_path_trj.xyz

The initial path is saved in basename_initial_path_trj.xyz, and it is recommended to check if that makes sense. In this case, we get the following sequence of eight "images":

../_images/init-path.gif — Animation: Initial path images from NEB-TS.#

which looks reasonable, with the proton being transferred from the water to the methanol as it leaves.

Note

Trajectory files are simply multixyz files containing several structures one after the other. Chemcraft, Molden and some other software can open those directly. Avogadro requires going to "Extensions" -> "Animation..." and the file can be opened and even animated.

Then, the iterative process starts:

Starting iterations:

Optim.  Iteration  HEI  E(HEI)-E(0)  max(|Fp|)   RMS(Fp)    dS 
Switch-on CI threshold               0.020000
   LBFGS     0      5    0.372340    0.129615   0.029150  12.1839
   LBFGS     1      5    0.352980    0.104320   0.024586  12.1349       
   LBFGS     2      5    0.340743    0.140351   0.026619  12.1315
   LBFGS     3      5    0.326472    0.097524   0.019838  12.1426
[...]

and after some convergence, a "climbing image" (CI) is found:

Image  6 will be converted to a climbing image in the next iteration (max(|Fp|) < 0.0200)

Optim.  Iteration  CI   E(CI)-E(0)   max(|Fp|)   RMS(Fp)    dS     max(|FCI|)   RMS(FCI)
Convergence thresholds               0.020000   0.010000            0.002000    0.001000
   LBFGS    49      6    0.130964    0.016940   0.003526  15.3056    0.017426    0.005167
   LBFGS    50      6    0.129316    0.045416   0.007419  15.4226    0.040132    0.010358
   LBFGS    51      6    0.126481    0.067751   0.011014  15.6647    0.050468    0.013683
   LBFGS    52      6    0.124766    0.022067   0.004391  15.6905    0.022698    0.006567
[...]

This is our first guess to the TS final structure.

Optimization of the TS#

After some iterations, the NEB-CI converges:

                                .--------------------.
          ----------------------| CI-NEB convergence |-------------------------
          Item                value                   Tolerance       Converged 
          ---------------------------------------------------------------------
          RMS(Fp)             0.0013467524            0.0100000000      YES
          MAX(|Fp|)           0.0142691584            0.0200000000      YES
          RMS(FCI)            0.0006034194            0.0010000000      YES
          MAX(|FCI|)          0.0016044259            0.0020000000      YES
          ---------------------------------------------------------------------  

The elastic band and climbing image have converged successfully to a MEP in 102 iterations!


                    *********************H U R R A Y*********************
                    ***        THE NEB OPTIMIZATION HAS CONVERGED     ***
                    *****************************************************

and a transition state optimization starts. First, a guess Hessian matrix is built using information from the previous NEB and a saddle point search is initiated.

If the TS is fully optimized, a summary for the whole NEB-TS is printed:

                    ***********************HURRAY************************
                    ***        THE TS OPTIMIZATION HAS CONVERGED      ***
                    *****************************************************

          
          
---------------------------------------------------------------
                      PATH SUMMARY FOR NEB-TS             
---------------------------------------------------------------
All forces in Eh/Bohr. Global forces for TS.

Image     E(Eh)   dE(kcal/mol)  max(|Fp|)  RMS(Fp)
  0    -344.08225     0.00       0.00014   0.00003
  1    -344.07675     3.45       0.00232   0.00071 
  2    -344.07387     5.26       0.01427   0.00363
  3    -344.07067     7.27       0.00282   0.00105
  4    -344.04841    21.23       0.00130   0.00058
  5    -344.01506    42.16       0.00124   0.00056
  6    -343.99728    53.32       0.00150   0.00060 <= CI
 TS    -343.99864    52.47       0.00013   0.00005 <= TS
  7    -344.03091    32.22       0.00151   0.00061 
  8    -344.06788     9.02       0.00179   0.00073
  9    -344.07676     3.45       0.00026   0.00009

and the frequencies of the TS are computed, to verify if that has only one negative frequency. The result for this case is:

-----------------------
VIBRATIONAL FREQUENCIES
-----------------------

Scaling factor for frequencies =  1.000000000  (already applied!)

     0:       0.00 cm**-1
     1:       0.00 cm**-1
     2:       0.00 cm**-1
     3:       0.00 cm**-1 
     4:       0.00 cm**-1
     5:       0.00 cm**-1
     6:   -1192.35 cm**-1 ***imaginary mode***   
     7:      40.41 cm**-1
     8:     136.38 cm**-1
     9:     209.70 cm**-1
[...]

proving that the optimized TS structure is indeed a saddle point.

The final trajectory from the reactants to products can be read from basename_MEP_trj.xyz, and the final transition state geometry, saved in basename_NEB-TS_converged.xyz is:

../_images/ts1.png — Figure: Optimized transition state structure.#

One can see that the TS has a tetrahedral carbon atom, as expected from the experimental mechanism, and there is a concerted transfer of proton to the leaving methoxy group.

Note

This calculation was run in gas phase to simplify the explanation. If one includes the solvent, the mechanism will be more complex and the proton transfer can occur, possibly, through an external water molecule.

Possible intermediates#

During the NEB-CI process, the following message could appear:

Possible intermediate minimum found at image(s): 7
  => have a look at the .interp file.
  => it is often a good idea to divide the path calculation into smaller fragments.

This message indicates that a local minimum was detected along the reaction path. Actually, there could be even more than one. If this is the case as demonstrated in the exemplary reaction diagram below, one should consider to further optimize these intermediates and rerun NEB-TS to localize their connecting transition states to obtain the full picture of the reaction mechanism.

../_images/interp-2.png — Figure: Exemplary reaction path with two intermediates.#

NEB for excited states#

The NEB and NEB-TS algorithm can also be used for excited states in combination with TD-DFT. To do so, simply add the respective %tddftblock to your input:

! NEB-TS

%NEB
 PRODUCT "product.xyz"
END
%TDDFT
 NROOTS 1
 IROOT 1
END

Finding transition states with NEB-TS#

Defining your reactants and products#

The NEB-TS input#

Running the NEB-TS#

Optimization of the TS#

Possible intermediates#

NEB for excited states#

Structures#