Transition State search using CI-NEB by VASP-VTST

Xijun Wang

12/3/2021

The transition state (TS) of a chemical reaction is a particular configuration along the reaction coordinate. It is defined as the state corresponding to the highest potential energy along this reaction coordinate. Left figure is an illustration of the potential energy surface (PES) of a system, where the two local minimas (reactant and product) are connected by a red-marked reaction coordinate. The maximum along such a path is called a TS, which is a first order saddle point, a maximum in one coordinate and minima in all others.

Climbing Image Nudged Elastic Band (CI-NEB)

The nudged elastic band (NEB) is a widely-used method for searching TS on the PES. The method works by optimizing a number of intermediate images along the reaction path connecting reactants and products. Each image finds the lowest energy possible while maintaining equal spacing to neighboring images. This constrained optimization is done by adding spring forces along the band between images and by projecting out the component of the force due to the potential perpendicular to the band.

The climbing image is an advanced version of the NEB method in which the highest energy image is driven up to the saddle point. In stead of adding spring forces, the inverted true forces along the tangent is added to this image. Such that the image's energy can be maximized along the band, and minimized in all other directions. With this approach, the exact saddle point (TS) can be found when this image converges.

Because the highest image is moved to the saddle point and it does not feel the spring forces, the spacing of images on either side of this image will be different. It can be important to do some minimization with the regular NEB method before this flag is turned on, both to have a good estimate of the reaction co-ordinate around the saddle point, and so that the highest image is close to the saddle point. If the maximum image is initially very far from the saddle point, and the climbing image was used from the outset, the path would develop very different spacing on either side of the saddle point.

Installation of VTST codes into VASP

I would like to skip this step because:

1. VTST is so popular that it can be easily loaded on most clusters by some simple module commands like: module load vasp/5.3.3.vtst.

2. Very detialed step-by-step tutorials has been provied by the official site of VTST:

https://theory.cm.utexas.edu/vtsttools/installation.html

Example: NH3 activation on the surface of V2O5.

Using CI-NEB implemented by VASP-VTST, the target of this example is to search the involved TS and calculate the energy required for the NH3 activation (*NH3 *NH2 + *H), in which a proton transfers from the N site on NH3 to an oxygen site on the surface of V2O5.

Step 1. Prepartion

Download all the wide-used scripts from the office-site of VTST:

https://theory.cm.utexas.edu/vtsttools/scripts.html

tar -xvf vtstscripts.tgz

cp * ~/bin

chmod +x ~/bin/*

Step 2. IS and FS Optimization

Optimize the initial structure (IS, *NH3) and final structure (FS, *NH2 + *H) using VASP respectively prior to TS search. Copy the generated CONTCAR of IS and FS to a new directory as IS.vasp and FS.vasp, respectively.

cp IS/CONTCAR /your_directory/IS.vasp

cp FS/CONTCAR /your_directory/FS.vasp

Step 3. Check the serial number of atoms in IS.vasp and FS.vasp

It is required that the serial number of atoms in IS.vasp and FS.vasp are the same. A script ~/bin/dist.pl can be used to check how similar the IS.vasp and the FS.vasp are, which generates the average distance of the same atom in IS.vasp and FS.vasp. In most cases the output value is smaller than 5 Å .

In addition, it is required to replace the first line of IS.vasp with the atomic type in the order of POTCAR as shown below.

[xwang99@login04 test]$ dist.pl IS.vasp FS.vasp

2.24955063206987

[xwang99@login04 test]$ head IS.vasp

V O N H # The first line of IS.vasp must be the atomic type in the order of POTCAR

1.00000000000000

11.5241003035999992 0.0000000000000000 0.0000000000000000

0.0000000000000000 12.6374998092999995 0.0000000000000000

0.0000000000000000 0.0000000000000000 20.0000000000000000

V O N H

18 45 1 3

Selective dynamics

Direct

0.1495900019999965 0.0364399999999989 0.0001299999999986 F F F


Step 4. Interpolate images between IS and FS

This step can be easily achieved by using ~/bin/nebmake.pl IS.vasp FS.vasp N, where N is the number of intermediate images as shown below. The generated sub-directories 00 contains the POSCAR of the initial geometry and 0(N+1) contains the POSCAR of the final geometry. The sub-directories 01 - 0N contains the POSCAR of intermediate images. After generating the sub-directories 01/, 02/, 03/, ... check the validity of these intermediate images (POSCAR) by visualizing them in VESTA.

[xwang99@login04 test]$ nebmake.pl IS.vasp FS.vasp 3 # Here 3 is the number of images, which can be replaced with any number you want.

filetype1: vasp5

filetype2: vasp5


OK, ALL SETUP HERE

FOR LATER ANALYSIS, PUT OUTCARs IN FOLDERS 00 and 04 !!!

Step 5. Prepare INCAR, POTCAR, KPOINTS

In INCAR, it is required to add/modify the related keywords of CI-NEB, including but not limitted to:

-------------------------------------------------------------------------------------------------------------------------

#NEB

LCLIMB=.True.

ICHAIN=0 # Indicates which method to run. NEB (ICHAIN=0) is the default

IOPT=0 # = 1, 2, 3, 4, 7 ... for different optimizers

IBRION=1 # If you want to use IOPT = 1, 2, 3 ..., set IBRION=3

POTIM=0.2 # Change the step size based on your systems

IMAGES = 3 # The number of images

SPRING = -5 # The spring constant, in eV/Ang^2 between the images; negative value turns on nudging

EDIFFG=-0.02 # Convergence criterion of force

NSW=500 # The maximum steps allowed

ISIF=2 # Fix the lattice constants or not

-------------------------------------------------------------------------------------------------------------------------

POTCAR and KPOINTS are the same as structural optimization of IS and FS. Also, a script is required to submit the job.

After preparing all the following files, submit the job.

[xwang99@login04 test]$ ls

00 01 02 03 04 FS.vasp INCAR IS.vasp KPOINTS POTCAR script-vasp5.4-NEB

Step 6. Check the output files

After the job ends normally, copy the IS/OUTCAR and FS/OUTCAR into /your_directory/00/ and 04/, respectively, which could make it easier to obtain the energy of IS, TS, and FS.

Subsequently, double-check the validity of the obtained images (CONTCAR) in 01/, 02/, and 03/ by visualizing them in VESTA.

If everything looks good, use nebresults.pl to generate all the information about this job. As an example, the highlighted information are the energies of the IS.vasp, interpolated images. and FS.vasp, in which the one with the highest energy is the TS.

[xwang99@login04 test]$ nebresults.pl


Unziping the OUTCARs ... done

Do nebbarrier.pl ; nebspline.pl

Do nebef.pl

Do nebmovie.pl

Do nebjmovie.pl

Do nebconverge.pl

Warning: empty x range [0:0], adjusting to [-1:1]

Warning: empty y range [0.035996:0.035996], adjusting to [0.035636:0.036356]

Warning: empty y2 range [-463.462:-463.462], adjusting to [-468.097:-458.828]

Warning: empty x range [0:0], adjusting to [-1:1]

Warning: empty y range [0.049215:0.049215], adjusting to [0.0487229:0.0497071]

Warning: empty y2 range [-462.495:-462.495], adjusting to [-467.12:-457.87]

Warning: empty x range [0:0], adjusting to [-1:1]

Warning: empty y range [0.040008:0.040008], adjusting to [0.0396079:0.0404081]

Warning: empty y2 range [-462.893:-462.893], adjusting to [-467.522:-458.264]

Zipping the OUTCARs again ... done


Forces and Energy:

0 0.000000 -463.541800 0.000000

1 0.015996 -463.462200 0.079600

2 0.019215 -462.495000 1.046800

3 0.010008 -462.893300 0.648500

4 0.000000 -462.999100 0.542700


Extremum 1 found at image 1.999395 with energy: 1.046749

Extremum 2 found at image 3.000000 with energy: 0.648451

Extremum 3 found at image 3.666667 with energy: 159015.437564


Some Empirical Tips for TS searching

  1. In most cases, the automatically interpolated images are not as good as the ones you built based on your chemical intuition. Feel free to adjust the interpolated images.

  2. I would suggest to try IOPT=0 IBRION=1 first. If it doesn't work properly, try IOPT=3 IBRION=3 and then try IOPT=2, 4, ...

  3. In some cases, it would be efficient to use IOPT=3 IBRION=3 to find a rough TS configuration, and then switch to IOPT=0 IBRION=1 to find the exact TS. Remenber to cp CONTCAR POSCAR in each interpolated image directory 01/, 02/, ...

  4. Keep monitoring the status of the jobs.