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:

Maybe you need to investigate the mechanism, or even to compute the \(\Delta G^\dagger\) 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". Please, have in mind that here it means everything that is left and right of the arrow, the collection of molecules. First, for the reactant, after drawing the initial guess and doing a Geometry optimization using:

!B3LYP DEF2-SVP OPT D4

one gets, for instance, the following geometry:

with the nucleophilic water making a hydrogen bond with the carbonyl group. Rearranging the atoms to make a guess product structure, and after optimizing with the same method, one gets:

and we are ready to find the transition state.

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.

Note

The idea is not to optimize all molecules separately and just join then in a single file, but rather already build your "reactant" and "product" adducts in a certain orientation that makes sense with respect to the expected reaction, like done in the example.

Making the NEB-TS input#

The NEB-TS input is quite simple, one has to add only NEB-TS in the main input, and give the name of the product xyz file on the %NEB options:

!B3LYP DEF2-SVP D4 NEB-TS FREQ
%NEB NEB_END_XYZFILE "product.xyz" END
* XYZFILE 0 1 reactant.xyz

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

!B3LYP DEF2-SVP D4 NEB-TS FREQ
%NEB NEB_END_XYZFILE "product.xyz"
     NEB_TS_XYZFILE "guessTS.xyz"
END
* XYZFILE 0 1 reactant.xyz

and you can already run it.

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, set PREOPT_ENDS TRUE under %NEB.

Running the NEB#

After the calculation starts, the first output is:

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

which contains some details of the method. Then, comes the construction of a initial path:

Generation of  the initial path:
Minimize RMSD between reactant and product configurations     .... done
       RMSD before minimization         ....  2.1153 Bohr
       RMSD after minimization          ....  1.6611 Bohr
Performing linear interpolation         ....  done
Interpolation using                     ....  IDPP (Image Dependent Pair Potential)
IDPP-Settings:
       Remove global transl. and rot.   ...   true
       Convergence tolerance            ....  0.0100 1/Ang.^3
       Max. numer 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 117 iterations
Displacement along initial path:             11.6511 Bohr
Writing initial trajectory to file      ....  NEBTS_initial_path_trj.xyz

The idea here is to use a method called IDPP to create a series of images, from the reactant to products, that will be optimized together using the NEB formalism.

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":

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.297110    0.116137   0.026058  11.6511
   LBFGS     1      5    0.277441    0.092448   0.021661  11.5993
   LBFGS     2      6    0.262453    0.111904   0.022265  11.5925
   LBFGS     3      6    0.252519    0.082144   0.017474  11.6060

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    23      6    0.160709    0.029259   0.004374  12.8932    0.030101    0.007764
   LBFGS    24      6    0.159481    0.051391   0.006905  13.0230    0.054092    0.012911
   LBFGS    25      6    0.156310    0.068675   0.009114  13.3077    0.072840    0.017253
   LBFGS    26      6    0.153131    0.031865   0.004433  13.3046    0.031225    0.007908

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.0006879134            0.0100000000      YES
          MAX(|Fp|)           0.0032907070            0.0200000000      YES
          RMS(FCI)            0.0006591030            0.0010000000      YES
          MAX(|FCI|)          0.0016597640            0.0020000000      YES
          ---------------------------------------------------------------------

The elastic band and climbing image have converged successfully to a MEP in 128 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     0.000    -344.38971      0.00       0.00014   0.00003
  1     2.024    -344.38507      2.91       0.00144   0.00053
  2     3.395    -344.38062      5.70       0.00150   0.00055
  3     4.863    -344.38237      4.61       0.00167   0.00064
  4     5.663    -344.36742     13.98       0.00329   0.00102
  5     6.043    -344.33717     32.97       0.00274   0.00097
  6     6.336    -344.32049     43.44       0.00178   0.00058 <= CI
  7     6.754    -344.34916     25.45       0.00260   0.00098
  8     7.426    -344.38096      5.49       0.00143   0.00058
  9     9.661    -344.38666      1.91       0.00023   0.00007

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.00 cm**-1
        0.00 cm**-1
        0.00 cm**-1
        0.00 cm**-1
        0.00 cm**-1
        0.00 cm**-1
    -1149.57 cm**-1 ***imaginary mode***
      106.27 cm**-1
      135.34 cm**-1
      188.41 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:

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.

That was because a local minimum was detected along the reaction path. Actually, there could be even more than one.

The energies following the reaction minimum energy path are printed in the "interp" files and can be plotted to have a better look at it:

Here one clearly sees that there are two intermediates and one should optimize this structures as well in order to predict reaction energies correctly, and maybe rerun the NEB-TS with these new rectant and product guesses.

Finding transition states with NEB-TS#

Defining your reactants and products#

Making the NEB-TS input#

Running the NEB#

Optimization of the TS#

Possible intermediates#

Starting structures#