Part 1: Basic transition states ¶

Content¶

 1 Reaction path with nudged elastic band
 2 Improved dimer method
 3 H₂/Cu(111) elbow plot
 4 Intrinsic reaction coordinates: S_N2

1 Reaction path with nudged elastic band

$\uparrow$

¶

By the end of this tutorial, you will be able to:

• prepare and perform a nudged elastic band (NEB) calculation
• optimize the image structures and the reaction pathway
• visualize the image structures and energetic pathway

1.1 Task¶

Calculate the reaction pathway using nudged elastic band with 4 images for the self-diffusion of a Si atom to a vacancy site in a 15-atom Si primitive supercell.

All chemical reactions go from reactant to product via a transition state. One approach for modeling transition states is the nudged elastic band (NEB) method [J. Chem. Phys. 113, 9978 (2000)]. NEB takes the reactant and product structures and interpolates points between them to model the minimum energy path (MEP), i.e. the reaction pathway. The points consist of images, typically 4-20, gradually progressing from the initial (reactant) to the final (product) structures. In this exercise, we will look at a sparse set of only 4 images. The more images used, the more expensive the calculation and the greater the difficulty with convergence. The images are connected to each other, as well as to the initial and final structures, by spring interactions. Connecting these images creates a continuous "elastic" band.

The force of the band is described by the conventional ("true force") perpendicular to the band $F_{i, \perp }$ and the spring force parallel to the band $F_{i, spring}$, referred to as "nudging". The band is relaxed according to its forces $F_{i}$ and the reaction pathway is thereby obtained.

In this exercise, we will apply NEB to a 15-atom Si vacancy supercell shown below. The blue and cyan atoms are Si in the cell, while the red shows the position of the vacancy site.

The Si atom in blue will move to the position of the red vacancy site during the NEB calculation, as shown by the arrow. The diffusion between sites in diamond cubic silicon is known as self-diffusion [Phys. Rev. B 59, 3969 (1999)] and is simple to model. In this exercise, we will investigate the reaction pathway associated with the self-diffusion. There is additional information on NEB and other transition state methods on the VASP Wiki.

1.2 Input¶

The input files to run this example are prepared at $TUTORIALS/transition_states/e01_NEB.

In this exercise, the POSCAR files are provided. The structure of the vacancies (V1 and V2) and images have been pre-relaxed and so the relaxation will only take a few NEB steps.

For reference, a build and interpolation script are provided below:

from ase.build import bulk
from ase.build import make_supercell
from ase.visualize import view
from ase.io.vasp import write_vasp
import numpy as np

# Build the Si diamond primitive unit cell
Si = bulk('Si', 'diamond', a=5.47)
N = 2
M = [[N, 0, 0], [0, N, 0], [0, 0, N]]

# Create a 16-atom supercell
Si_Bulk = make_supercell(Si, M)
atoms = Si_Bulk.get_positions()

# Defect migration in crystalline silicon, Lindsey J. Munro and David J. Wales
# https://journals.aps.org/prb/abstract/10.1103/PhysRevB.59.3969

# Remove one atom to make an initial vacancy
Si_V1 = make_supercell(Si, M) # Rebuild the supercell - a copy points to the original in the memory, so a new supercell must be created
del Si_V1[-1] # Remove the last atom
write_vasp("./e01_NEB_coarse/POSCAR.V1", Si_V1, direct = True) # Write the output in VASP format

# Remove an alternative, adjacent atom to make a second vacancy
Si_V2 = make_supercell(Si, M) # Rebuild the supercell
del Si_V2[-2] # Remove the second to last atom
write_vasp("./e01_NEB_coarse/POSCAR.V2", Si_V2, direct = True) # Write the output in VASP format

# Uncomment these lines to view the structures!
#view(Si_Bulk, viewer='x3d')
#view(Si_V1, viewer='x3d')
#view(Si_V2, viewer='x3d')

System = Si
Electronic minimisation
ENCUT  = 250
PREC = Normal                # precision tag
EDIFFG = -0.01
EDIFF = 1e-6

DOS related values
ISMEAR = 0                   # Gaussian smearing
SIGMA = 0.05                 # Smearing width

Ionic relaxation
NSW = 500                    # Max number of ionic steps
IBRION=2                     # Conjugate gradient algorithm
ISIF=0                       # Cell shape and volume fixed, ions free

from ase.io.vasp import read_vasp
from ase.io.vasp import write_vasp
from ase.mep.neb import NEB

# Read final and initial vacancy structures from their respective POSCAR files V1 and V2
initial = read_vasp("./e01_NEB/POSCAR.V1")
final = read_vasp("./e01_NEB/POSCAR.V2")

# Make a band image of 4 files and concatenate them into a list
no_images = 4 
images = [initial]
images += [initial.copy() for i in range(no_images)]
images += [final]
neb = NEB(images) # convert into ASE format

# Interpolate the four middle images from the initial and final
neb.interpolate(method="idpp", apply_constraint=True) # idpp ensures that the image atoms do not lie too close to other atoms; apply_constraint turns reads for selective dynamics, if you have set them
a=0
# Write out the POSCAR for each image
for image in images:
    a += 1
    write_vasp(("./e01_NEB/0"+str(a-1)+"/POSCAR"), image)
    if a == no_images+2: break

POSCAR

Si
 1.0000000000000000
     5.4699999999999998    5.4699999999999998    0.0000000000000000
     0.0000000000000000    5.4699999999999998    5.4699999999999998
     5.4699999999999998    0.0000000000000000    5.4699999999999998
 Si
  15
Cartesian
  0.0000000000000000  0.0000000000000000  0.0000000000000000
  1.3671586518987400  1.3671586518987400  1.3671586518987400
  2.7349454006694627  0.0000000000000000  2.7349454006694627
  4.1029484152820723  1.3648589317780140  4.1029484152820723
  0.0000000000000000  2.7349454006694627  2.7349454006694627
  1.3648589317780140  4.1029484152820723  4.1029484152820723
  2.7382543911184447  2.7382543911184447  5.4765087822368894
  4.1112502755029183  4.1112502755029183  6.8426350428209535
  2.7349454006694627  2.7349454006694627  0.0000000000000000
  4.1029484152820723  4.1029484152820723  1.3648589317780140
  5.4765087822368894  2.7382543911184447  2.7382543911184447
  6.8426350428209535  4.1112502755029174  4.1112502755029183
  2.7382543911184447  5.4765087822368894  2.7382543911184447
  4.1112502755029174  6.8426350428209535  4.1112502755029174
  5.4967101394738123  5.4967101394738123  5.4967101394738123

Si
 1.0000000000000000
     5.4699999999999998    5.4699999999999998    0.0000000000000000
     0.0000000000000000    5.4699999999999998    5.4699999999999998
     5.4699999999999998    0.0000000000000000    5.4699999999999998
 Si
  15
Selective Dynamics
Cartesian
  0.0000000000000000  0.0000000000000000  0.0000000000000000 F F F
  1.3670300457234879  1.3670300457234879  1.3670300457234879 T T T
  2.7344532895378659  0.0000000000000000  2.7344532895378659 T F T
  4.1023111536722139  1.3645333876656973  4.1023111536722139 T T T
  0.0000000000000000  2.7344532895378659  2.7344532895378659 F T T
  1.3645333876656973  4.1023111536722139  4.1023111536722139 T T T
  2.7368303796572948  2.7368303796572948  5.4736607593145896 T T T
  4.1105054591614119  4.1105054591614119  6.8390331629742418 T T T
  2.7344532895378659  2.7344532895378659  0.0000000000000000 T T F
  4.1023111536722139  4.1023111536722139  1.3645333876656973 T T T
  5.4736607593145896  2.7368303796572948  2.7368303796572948 T T T
  6.8390331629742418  4.1105054591614110  4.1105054591614119 T T T
  2.7368303796572948  5.4736607593145896  2.7368303796572948 T T T
  4.1105054591614110  6.8390331629742418  4.1105054591614110 T T T
  5.7629786857685836  5.7629786857685836  5.7629786857685836 T T T

Si
 1.0000000000000000
     5.4699999999999998    5.4699999999999998    0.0000000000000000
     0.0000000000000000    5.4699999999999998    5.4699999999999998
     5.4699999999999998    0.0000000000000000    5.4699999999999998
 Si
  15
Selective Dynamics
Cartesian
  0.0000000000000000  0.0000000000000000  0.0000000000000000 F F F
  1.3669014395482357  1.3669014395482357  1.3669014395482357 T T T
  2.7339611784062692  0.0000000000000000  2.7339611784062692 T F T
  4.1016738920623554  1.3642078435533804  4.1016738920623554 T T T
  0.0000000000000000  2.7339611784062692  2.7339611784062692 F T T
  1.3642078435533804  4.1016738920623554  4.1016738920623554 T T T
  2.7354063681961454  2.7354063681961454  5.4708127363922907 T T T
  4.1097606428199054  4.1097606428199054  6.8354312831275301 T T T
  2.7339611784062692  2.7339611784062692  0.0000000000000000 T T F
  4.1016738920623554  4.1016738920623554  1.3642078435533804 T T T
  5.4708127363922907  2.7354063681961454  2.7354063681961454 T T T
  6.8354312831275301  4.1097606428199045  4.1097606428199054 T T T
  2.7354063681961454  5.4708127363922907  2.7354063681961454 T T T
  4.1097606428199045  6.8354312831275301  4.1097606428199045 T T T
  6.0292472320633550  6.0292472320633550  6.0292472320633550 T T T

Si
 1.0000000000000000
     5.4699999999999998    5.4699999999999998    0.0000000000000000
     0.0000000000000000    5.4699999999999998    5.4699999999999998
     5.4699999999999998    0.0000000000000000    5.4699999999999998
 Si
  15
Selective Dynamics
Cartesian
  0.0000000000000000  0.0000000000000000  0.0000000000000000 F F F
  1.3667728333729838  1.3667728333729838  1.3667728333729838 T T T
  2.7334690672746720  0.0000000000000000  2.7334690672746720 T F T
  4.1010366304524979  1.3638822994410638  4.1010366304524979 T T T
  0.0000000000000000  2.7334690672746720  2.7334690672746720 F T T
  1.3638822994410638  4.1010366304524979  4.1010366304524979 T T T
  2.7339823567349955  2.7339823567349955  5.4679647134699909 T T T
  4.1090158264783989  4.1090158264783989  6.8318294032808193 T T T
  2.7334690672746720  2.7334690672746720  0.0000000000000000 T T F
  4.1010366304524979  4.1010366304524979  1.3638822994410638 T T T
  5.4679647134699909  2.7339823567349955  2.7339823567349955 T T T
  6.8318294032808193  4.1090158264783989  4.1090158264783989 T T T
  2.7339823567349955  5.4679647134699909  2.7339823567349955 T T T
  4.1090158264783989  6.8318294032808193  4.1090158264783989 T T T
  6.2955157783581264  6.2955157783581264  6.2955157783581264 T T T

Si
 1.0000000000000000
     5.4699999999999998    5.4699999999999998    0.0000000000000000
     0.0000000000000000    5.4699999999999998    5.4699999999999998
     5.4699999999999998    0.0000000000000000    5.4699999999999998
 Si
  15
Selective Dynamics
Cartesian
  0.0000000000000000  0.0000000000000000  0.0000000000000000 F F F
  1.3666442271977317  1.3666442271977317  1.3666442271977317 T T T
  2.7329769561430752  0.0000000000000000  2.7329769561430752 T F T
  4.1003993688426394  1.3635567553287469  4.1003993688426394 T T T
  0.0000000000000000  2.7329769561430752  2.7329769561430752 F T T
  1.3635567553287469  4.1003993688426394  4.1003993688426394 T T T
  2.7325583452738460  2.7325583452738460  5.4651166905476920 T T T
  4.1082710101368924  4.1082710101368924  6.8282275234341077 T T T
  2.7329769561430752  2.7329769561430752  0.0000000000000000 T T F
  4.1003993688426394  4.1003993688426394  1.3635567553287469 T T T
  5.4651166905476920  2.7325583452738460  2.7325583452738460 T T T
  6.8282275234341077  4.1082710101368924  4.1082710101368924 T T T
  2.7325583452738460  5.4651166905476920  2.7325583452738460 T T T
  4.1082710101368924  6.8282275234341077  4.1082710101368924 T T T
  6.5617843246528968  6.5617843246528968  6.5617843246528968 T T T

Si
 1.0000000000000000
     5.4699999999999998    5.4699999999999998    0.0000000000000000
     0.0000000000000000    5.4699999999999998    5.4699999999999998
     5.4699999999999998    0.0000000000000000    5.4699999999999998
 Si
  15
Cartesian
  0.0000000000000000  0.0000000000000000  0.0000000000000000
  1.3665156210224796  1.3665156210224796  1.3665156210224796
  2.7324848450114785  0.0000000000000000  2.7324848450114785
  4.0997621072327810  1.3632312112164302  4.0997621072327810
  0.0000000000000000  2.7324848450114785  2.7324848450114785
  1.3632312112164302  4.0997621072327810  4.0997621072327810
  2.7311343338126961  2.7311343338126961  5.4622686676253922
  4.1075261937953860  4.1075261937953860  6.8246256435873960
  2.7324848450114785  2.7324848450114785  0.0000000000000000
  4.0997621072327810  4.0997621072327810  1.3632312112164302
  5.4622686676253922  2.7311343338126961  2.7311343338126961
  6.8246256435873960  4.1075261937953860  4.1075261937953860
  2.7311343338126961  5.4622686676253922  2.7311343338126961
  4.1075261937953860  6.8246256435873960  4.1075261937953860
  6.8280528709476682  6.8280528709476682  6.8280528709476682

INCAR.sp

System = Si
! Electronic minimization
ENCUT = 250.0                # plane wave energy cutoff in eV
EDIFF = 1e-6                 # electronic convergence criteria

! DOS related values
ISMEAR = 0                   # Gaussian smearing
SIGMA = 0.05                 # Smearing width

! single point
NSW = 0                      # number of maximum ionic steps

INCAR.neb

System = Si
! Electronic minimization
ENCUT  = 250
PREC = Normal                # precision tag
EDIFFG = -0.04
EDIFF = 1e-6

! DOS related values
ISMEAR = 0                   # Gaussian smearing
SIGMA = 0.05                 # Smearing width

! Ionic relaxation
NSW = 100                    # Max number of ionic steps
IBRION = 1                   # Quasi-Newton algorithm
ISIF = 2                     # Cell shape and volume fixed, ions free

! NEB
IMAGES = 4                   # 4 intermediate geometries for the NEB
SPRING = -5                  # Spring constant

KPOINTS

K-Points
 0
Gamma
 2  2  2
 0  0  0

POTCAR

 Pseudopotentials of Si.

Check the tags set in the INCAR.sp and INCAR.neb files!

The INCAR.sp file is for single-point calculations, while the INCAR.neb file is for performing the NEB calculation. Since they both use identical tags, only INCAR.neb will be gone through here.

The first four tags define the electronic minimization: ENCUT defines the plane wave basis, PREC defines how precise the calculations are, EDIFFG is the break condition for ionic relaxation, and EDIFF defines the break condition for the global, electronic self-consistency step. EDIFFG is the break condition for ionic relaxation. Negative values indicate a force threshold in eV Å^-1; positive values lead to a break when the energy difference between two ionic steps is less than it. We recommend using force thresholds. The smearing is defined by ISMEAR and SIGMA.

The remaining tags specify the settings for the ionic relaxation. NSW is the maximum number of ionic steps taken. For a single-point calculation, this is set to 0. A quasi-Newton algorithm is set up for ionic relaxation using IBRION = 1. POTIM defines the ionic step width during the relaxation. ISIF = 2 allows the ionic positions to relax while keeping the cell shape and size fixed. You might see "BRIONS problems: POTIM should be increased" in the output. This indicates difficulties with the convergence. It may be fixed by increasing POTIM or by changing the algorithm used by IBRION.

The final two tags define the NEB. IMAGES sets the number of images (intermediate structures) used in the NEB calculation, in this case 4 (note that the MPI rank must be a product of the number of images). SPRING = -5 sets the spring constant for NEB.

The KPOINTS file specifies a Γ-centered $\mathbf{k}$-mesh.

The settings used in this tutorial are sufficient to run calculations quickly. For higher-quality results, increase the k-point mesh and the supercell size until results, e.g. formation energies, converge. EDIFFG should also be tightened, e.g. to -0.01 eV Å^-1

1.3 Calculation¶

Open a terminal and navigate to this example's directory, then prepare the following:

cd $TUTORIALS/transition_states/e01_*

We will perform NEB using four images. Alongside the initial and final structure, there are six directories in total: initial (00), sequential images (01, 02, 03, 04), and final (05). Copies of the corresponding POSCAR files are provided with the other input files. The number of images is defined by IMAGES in the INCAR file.

An NEB calculation must be run from the parent ./e01_NEB directory. The calculations for each image are then performed in the sub-directories by VASP. Now that you have finished preparing, you can run the NEB calculation by entering the ./e01_NEB directory and running VASP:

cd $TUTORIALS/transition_states/e01_*
mpirun -np 4 vasp_std

Note that the stdout printed on the screen is from image 1. The stdout for the remaining images are found in their respective directories.

You have completed your first NEB calculation. Well done!

You have the energy of the NEB images as a Si atom diffuses to a neighboring vacancy site. The energies of the initial and final structures will also be required to investigate the energy profile, as these structures are fixed during the NEB calculation. Let us obtain these by doing single-point calculations for each of these. First, the initial structure:

cd $TUTORIALS/transition_states/e01_*
cp INCAR.sp 00/INCAR
cp POTCAR KPOINTS 00
cd ./00
mpirun -np 4 vasp_std

Then for the final structure:

cd $TUTORIALS/transition_states/e01_*
cp INCAR.sp 05/INCAR
cp POTCAR KPOINTS 05
cd ./05
mpirun -np 4 vasp_std

Due to the cheap settings used in this example, the final energies will not match those seen experimentally.

! ab initio
ISMEAR= -5          ! tetrahedron method with Bloechl corrections

ENCUT = 250.0       ! plane wave energy cutoff in eV
EDIFF = 1e-6        ! electronic convergence relaxation

With all the calculations finished, it is helpful to look at how the structure changes. Begin by taking a look at the 3D structure using py4vasp:

The diffusion of the Si atom is largely along a straight line. By taking a cross-section of the Si supercell and projecting the positions onto it, the progression of the diffusion Si atom can be readily visualized from its initial position to the position of the original defect. The following script plots the projected structure of each image onto a pre-defined 2D-plane:

Most of the Si atoms remain in the same position for all the images. The diffusing Si atom is very distinct as it progresses from its initial position to occupy the vacancy site. The change of structure is clear throughout the process. Each image seems to be a similar distance apart. Let us see how the energy changes between images across the reaction profile using py4vasp:

The energetic profile is roughly symmetric with a peak ~0.7 eV. If we were to further investigate this transition, we would add additional images in between the four we currently have, to produce a fuller picture of the transition. Note that better convergence in a single calculation is achieved using a small number of images, e.g. up to four. If too many are used, however, then the optimization can become very computationally expensive.

With the example of Si diffusion across vacancies, you have now seen that it can be a complex process investigating transition states. You could take an intermediate image as a trial structure for finding the transition state. In the next section, you will learn how to optimize a trial transition state.

1.4 Questions¶

1. Which tags do you need to specify an NEB calculation?
2. What is the pathway that NEB finds?
3. Why should you not optimize too many images at once?

2 Improved dimer method

$\uparrow$

¶

By the end of this tutorial, you will be able to:

• calculate imaginary modes using a vibrational calculation
• optimize the transition state for H₂ dissociation using the improved dimer method (IDM)
• check for imaginary frequencies for the transition state

2.1 Task¶

Optimize the transition state with the improved dimer method for H₂ dissociation on Cu(111) using a 3-layered (2x2) supercell.

NEB allows you to find the reaction pathway. Other methods must be used to find the transition state (TS) itself and optimize along the reaction pathway to find the energetic maximum at a saddle point where there is only one imaginary frequency, i.e. the TS. One such method is the improved dimer method (IDM) [J. Chem. Phys. 111 (1999) 7010]. IDM is a local saddle-point search algorithm using the forces, as opposed to using the Hessian (i.e. the first derivative with respect to the energy and position rather than the second). The IDM proceeds in four steps illustrated and described below:

  a) An initial direction $u_{\xi}$ is taken from the most negative vibrational mode of the trial structure. The trial structure is the first point $q$.
  b) An additional point on the potential energy surface (PES) forward along the trial direction is defined $q + \delta u_{\xi}$. The first and second points together define the dimer.
  c) The dimer is then rotated on the PES about $q$ by angle $\phi_1$ to $\tilde{q}$.
  d) A new direction is defined by rotating $u_{\xi}$ by $\phi_{min}$ to minimize the negative curvature of the PES $\lambda$. A search direction $\bar{N}$ is defined using a minimization algorithm, then a translation is taken in the direction of $\bar{N}$ to $q + \epsilon \bar{N}$, where $\epsilon$ is a step distance.

The new structure $q + \epsilon \bar{N}$ is used to define a new dimer and steps a-d are repeated, rotation followed by translation, until the saddle point is reached, yielding a saddle point with only one negative vibrational mode, i.e. the transition state.

A simple example of reactive adsorption is the dissociative adsorption of H₂ on Cu(111). The dissociation has been well-studied, with good experimental and theoretical data [J. Chem. Theory Comput. 13 (2017) 3208–3219]. In this exercise, you will use IDM to relax the transition state.

2.2 Preparing structures¶

The input files to run this example are prepared in $TUTORIALS/transition_states/e02_IDM. Note that the additional lines in the POSCAR file after the positions are, generally, velocities.

POSCAR.idm

Cu H
   1.00000000000000
     5.1053109601668734    0.0000000000000000    0.0000000000000000
     2.5526554800834367    4.4213289857236360    0.0000000000000000
     0.0000000000000000    0.0000000000000000   17.0220689435490975
  Cu              H
    12     2
Selective dynamics
Direct
      0.16666667       0.16666667       0.37755692 F F F
      0.66666667       0.16666667       0.37755692 F F F
      0.16666667       0.66666667       0.37755692 F F F
      0.66666667       0.66666667       0.37755692 F F F
      0.83333333       0.33333333       0.50000000 F F F
      0.33333333       0.33333333       0.50000000 F F F
      0.83333333       0.83333333       0.50000000 F F F
      0.33333333       0.83333333       0.50000000 F F F
      0.99238939       0.00352484       0.62042734 T T T
      0.48015918       0.00298666       0.62122341 T T T
      0.98637412       0.50378217       0.62136799 T T T
      0.49056648       0.50136365       0.62111959 T T T
      0.68831531       0.88592943       0.69883863 T T T
      0.62492060       0.13310462       0.68602475 T T T
! Trial unstable direction
 0 0 0
 0 0 0
 0 0 0
 0 0 0
 0 0 0
 0 0 0
 0 0 0
 0 0 0
 0.027098 -0.003023 -0.015080
-0.022306 -0.004165 -0.029077
 0.004290 -0.027089 -0.023693
 0.003362  0.018196 -0.023398
 0.061633  0.599759  0.325308
-0.179319 -0.518661  0.474306

POSCAR.freq1

Cu H
   1.00000000000000
     5.1053109601668734    0.0000000000000000    0.0000000000000000
     2.5526554800834367    4.4213289857236360    0.0000000000000000
     0.0000000000000000    0.0000000000000000   17.0220689435490975
  Cu              H
    12     2
Selective dynamics
Direct
      0.16666667       0.16666667       0.37755692 F F F
      0.66666667       0.16666667       0.37755692 F F F
      0.16666667       0.66666667       0.37755692 F F F
      0.66666667       0.66666667       0.37755692 F F F
      0.83333333       0.33333333       0.50000000 F F F
      0.33333333       0.33333333       0.50000000 F F F
      0.83333333       0.83333333       0.50000000 F F F
      0.33333333       0.83333333       0.50000000 F F F
      0.99238939       0.00352484       0.62042734 T T T
      0.48015918       0.00298666       0.62122341 T T T
      0.98637412       0.50378217       0.62136799 T T T
      0.49056648       0.50136365       0.62111959 T T T
      0.68831531       0.88592943       0.69883863 T T T
      0.62492060       0.13310462       0.68602475 T T T

INCAR.freq1

System: fcc Cu (111), 3 layers
! Ab initio
PREC = Normal
EDIFF = 1e-6               # electronic convergence

! DOS
ISMEAR = 1                 # First-order Methfessel-Paxton smearing
SIGMA = 0.2                # 0.2 eV smearing width

! Frequency
IBRION = 5                 # Finite differences
NFREE = 2                  # No. displacements for finite differences
NSW = 1                    # Positive number of ionic steps

! Machine learned force field
ML_LMLFF = T               # Turns on Machine learned force field (MLFF)
ML_MODE = run              # Sets MLFF mode to running, rather than training

INCAR.idm

System: fcc Cu (111), 3 layers
! Ab initio
PREC = Normal
EDIFF = 1e-6              # electronic convergence
NELMIN = 5                # at least 5 el. scf steps for each ionic step

! DOS
ISMEAR = 1                # First-order Methfessel-Paxton smearing
SIGMA = 0.2               # 0.2 eV smearing width

! IDM
NSW = 20
IBRION=44                 # use the dimer method as optimization engine
EDIFFG=-0.01

KPOINTS

K-Points
 0
Gamma
 3  3  1
 0  0  0

POTCAR

 Pseudopotentials of Cu and H.

Check the tags set in the INCAR file! These have already been described in Example 1, with the following exceptions:

The frequency calculation is defined in INCAR.freq1. IBRION = 5 sets up a finite-difference calculation and NFREE the number of displacements. ML_LMLFF = T and L_MODE = run enables the use of a pre-trained machine-learned force field (MLFF). The IDM calculation is defined in INCAR.idm and set up with IBRION = 44.

Two remaining tags are used: NELMIN defines the minimum number of elementary electronic steps.

The KPOINTS file specifies a 2D $\mathbf{k}$-mesh. There is only the Γ-point in the direction perpendicular to the surface. Since this is a surface, we do not want interaction between the periodic images, therefore the minimum number of $\mathbf{k}$-points, one, and a large vacuum are used.

The settings used in this tutorial are sufficient to run quickly. For higher-quality results, converge the results, e.g. adsorption energies, with respect to the k-point mesh, the size of the vacuum, and the number of layers used (four or more is standard in quantum chemical calculations). Note that for the IDM, we have set the number of ionic steps NSW = 20. This was selected for speed, such that there is only one imaginary frequency remaining. True convergence would take many more steps, on the order of hundreds or thousands of steps.

2.3 Calculation¶

Open a terminal, navigate to this example's directory, and navigate to the directory:

cd $TUTORIALS/transition_states/e02_*

Once the calculation is finished, search for the frequency modes using the following:

grep "cm-1" $TUTORIALS/transition_states/e02_*/freq1/OUTCAR

The vibrational modes with their frequencies in several different units will be printed out:

   1 f  =   45.115633 THz   283.469886 2PiTHz 1504.895547 cm-1   186.583184 meV
   2 f  =   41.417254 THz   260.232281 2PiTHz 1381.530879 cm-1   171.287921 meV
   3 f  =   14.024817 THz    88.120523 2PiTHz  467.817533 cm-1    58.001956 meV
...
  15 f  =    1.895701 THz    11.911042 2PiTHz   63.233786 cm-1     7.839987 meV
  16 f  =    1.644956 THz    10.335564 2PiTHz   54.869830 cm-1     6.802989 meV
  17 f/i=    2.504934 THz    15.738964 2PiTHz   83.555604 cm-1    10.359570 meV
  18 f/i=   28.911252 THz   181.654755 2PiTHz  964.375566 cm-1   119.567278 meV

You can see that there are two imaginary modes. For a transition state, there should only be one. The IDM will relax along the larger mode, removing the smaller one.

These modes may be visualized using the following script:

python3 ../../scripts/vibFreq/vibFreq.py 0.0 0.2 0.0 > freq.dat

vibFreq.py generates the hesseModes_shifted.molden file. Download and visualize it in JMol or molden. The image below shows an image of the H₂ molecule on the Cu(111) surface with the strongest vibrational mode, the one that we are interested in, shown with a double-headed orange arrow.

The vibration is the stretch of the H-H bond. During dissociation, this bond breaks. Therefore, we will perform an IDM following this unstable vibrational mode to find the transition state for the dissociation on the surface.

2.3.2 Improved Dimer Method¶

For IDM we want to relax along the most unstable imaginary mode (largest frequency) toward the transition state. We take this direction from the unstable phonon modes calculated in Section 2.3.1. Here, we have only a few unstable modes because the structure is close to relaxed. Take the mode with the largest frequency using the following script. You can find the expected result in POSCAR.idm.

cd $TUTORIALS/e02_*/freq1/
grep -A 15 "18 f/i=" OUTCAR | awk '{print $4 " " $5 " " $6}' | tail -n 14 > temp
echo "! Unstable trial direction" >> POSCAR
cat temp >> POSCAR
rm temp
cp POSCAR ../idm

With the POSCAR.idm file containing the trial unstable direction, you now have everything required to run the IDM. Copy it to POSCAR and run the IDM relaxation with the following:

cd $TUTORIALS/transition_states/e02_*/idm
cp ../POSCAR.idm ./POSCAR
mpirun -np 4 vasp_std

This relaxation will last only 20 steps before terminating (defined by NSW), so will not be fully converged. However, you will perform a second frequency calculation in the next section to confirm that there is only one imaginary frequency in the relaxed structure.

Take a look at the structures with py4vasp:

Once this vibrational analysis is finished, you may search for frequency modes again:

cd $TUTORIALS/transition_states/e02_*/freq2/
grep "cm-1" OUTCAR

Resulting in the following output:

   1 f  =   50.817114 THz   319.293345 2PiTHz 1695.076468 cm-1   210.162602 meV
   2 f  =   38.497867 THz   241.889230 2PiTHz 1284.150605 cm-1   159.214312 meV
   3 f  =   15.383430 THz    96.656942 2PiTHz  513.135994 cm-1    63.620726 meV
...
  15 f  =    3.793943 THz    23.838049 2PiTHz  126.552329 cm-1    15.690482 meV
  16 f  =    1.905720 THz    11.973993 2PiTHz   63.567980 cm-1     7.881421 meV
  17 f  =    1.586471 THz     9.968091 2PiTHz   52.918973 cm-1     6.561114 meV
  18 f/i=   15.906183 THz    99.941492 2PiTHz  530.573138 cm-1    65.782656 meV

You can see that there is only one imaginary mode. This mode may be visualized using the following script:

python3 ../../scripts/vibFreq/vibFreq.py 0.0 0.2 0.0 > freq.dat

Congratulations, you have successfully found a transition state for the dissociation of H₂ on Cu(111)! The remaining imaginary frequency may be visualized using the following vibFreq.py, which generates the hesseModes_shifted.molden file containing vibrational data in molden format. Visualize it in JMol or molden once downloaded.

Let us take a look at the optimization:

Even from a small number of steps, it is clear that there are occasionally rapid drops in energy. For comparison, let us take a look at a longer IDM calculation for the same system:

POSCAR.idm file that you have calculated was taken from ~600 steps into this longer calculation. You can immediately see that there is a gradual decay before this, then a rapid drop is seen around 600 before plateauing. It is this drop where the additional imaginary frequency is removed.

You have successfully optimized the transition state! From here, you could perform an IRC calculation to obtain the product and reactant systems, which would enable comparison between barrier heights and experiment. The intrinsic reaction coordinate (IRC) will be gone through in Example 4.

2.4 Questions¶

1. What optimization approach does IBRION = 44 set?
2. How many imaginary frequencies should there be in the transition state?
3. What sort of calculation should you perform to obtain a trial imaginary frequency?

3 H₂/Cu(111) elbow plots

$\uparrow$

¶

By the end of this tutorial, you will be able to:

• prepare structures for an elbow plot
• visualize a 2D cross-section of the PES

3.1 Task¶

Create an elbow plot for H₂ dissociation on Cu(111) using a 3-layered (2x2) supercell.

IDM gives an optimized structure for the transition state. This is the midpoint of the reaction pathway between reactants and products, which may be calculated using NEB. Another technique may be used for diatomic molecules. They have six degrees of freedom and so a meaningful subset can be taken. With the surface fixed in place, a 2D potential energy surface (PES) for the diatomic can be produced with the height of the diatomic center of mass z and the bond length r as the two degrees of freedom. This 2D PES is called an elbow plot and shows a slice of the PES as a contour plot [Phys. Rev. Lett. 73, 1400 (1994), Phys. Rev. B 62, 8295].

Following on from the previous exercise, you will create an elbow plot for the dissociative adsorption of H₂ on Cu(111) in this one.

3.2 Preparing structures¶

The input files to run this example are prepared in $TUTORIALS/transition_states/e03_Elbow_plot.

POSCAR.surf

Cu 
 1.0000000000000000
     5.1053109601668734    0.0000000000000000    0.0000000000000000
     2.5526554800834367    4.4213289857236360    0.0000000000000000
     0.0000000000000000    0.0000000000000000   17.0220689435490975
 Cu  
  12   
Selective dynamics
Cartesian
  1.2763277400417001  0.7368881642872622  6.4268000000000303   F   F   F
  3.8289832201251368  0.7368881642872622  6.4268000000000303   F   F   F
  2.5526554800834185  2.9475526571490800  6.4268000000000303   F   F   F
  5.1053109601668556  2.9475526571490800  6.4268000000000303   F   F   F
  5.1053109601668920  1.4737763285745558  8.5110344717745487   F   F   F
  2.5526554800834549  1.4737763285745558  8.5110344717745487   F   F   F
  6.3816387002086099  3.6844408214363740  8.5110344717745487   F   F   F
  3.8289832201251732  3.6844408214363740  8.5110344717745487   F   F   F
  0.0000000000000000  0.0000000000000000 10.5987747510960109   T   T   T
  2.5526554800834367  0.0000000000000000 10.5987747510960109   T   T   T
  1.2763277400417183  2.2106644928618180 10.5987747510960109   T   T   T
  3.8289832201251550  2.2106644928618180 10.5987747510960109   T   T   T

INCAR

System: fcc Cu (111), 3 layers
! Ab initio
NSW = 0                    # Single ionic step, i.e. single point

! DOS
ISMEAR = 1                 # First-order Methfessel-Paxton smearing
SIGMA = 0.2                # 0.2 eV smearing width

LCHARG=.FALSE.             # Does not save the CHG or CHGCAR files to disk
LWAVE=.FALSE.              # Does not save the WAVECAR files to disk

KPOINTS

K-Points
 0
Gamma
 1  1  1
 0  0  0

POTCAR

 Pseudopotentials of Cu and H.

Check the tags set in the INCAR file! Two new tags are LCHARG and LWAVE. These determine whether or not the CHGCAR and CHG, and WAVECAR files, respectively, are written to disk. Since many calculations are being run in this exercise, these are set to ".FALSE." to conserve disk space. Otherwise, all INCAR tags are have been previously defined.

Note that NSW = 0 so that single-point calculations are performed.

A Γ-point only $\mathbf{k}$-mesh is used. We stress that the $\mathbf{k}$-mesh here is very poor and is only used for illustrative purposes in this tutorial.

3.3 Calculation¶

The structures for the elbow plot are created by taking a bare Cu(111) surface POSCAR.surf and placing H₂ on the surface, for a set of heights above the surface z and bond lengths r. To begin with, open a terminal and navigate to this example's directory:

cd $TUTORIALS/transition_states/e03_*/elbow/

In the following script, using ASE, the POSCAR files are generated for all H₂ heights above the surface z and bond lengths r. Use the following:

This generates a lot of structures. Prepare the directories for each structure using the following:

cd $TUTORIALS/transition_states/e03_*/elbow/
prefix=POSCAR. 
for a in POS*; do 
    temp=${a#"$prefix"}
    mkdir $temp
    cp $a $temp/POSCAR
    cp ../INCAR ../POTCAR ../KPOINTS $temp
done

They may then be run using the provided run script:

cd $TUTORIALS/transition_states/e03_*/elbow
bash run

#!/bin/bash
VASP="mpirun -np 4 vasp_gam"

for a in */; do
        cd $a
            echo $a
            $VASP
        cd ..
done

You have successfully generated all of the necessary data. Now let us produce a contour plot using it! The following script extracts the energies for each structure with the corresponding z and r coordinates. It first loops over the range of z- and r-positions, then uses py4vasp to extract the energy for each calculation. The energy is then presented as a contour plot on a mesh z-r grid.

Run the following script to create the elbow plot:

Congratulations! You've produced your first elbow plot! Below is an example using tighter settings. Run the following script to visualize it. What differences can you see?

The two plots look visually similar. However, the second plot is an 11² point grid, so there are 121 points in total. More points mean that the contours are less jagged and smoother as they offer a better description of the PES. Also, energetically the latter is far more negative, due to a denser $\mathbf{k}$-point mesh being used, 3x3x1, rather than only the Γ-point. For fully converged results, much finer meshes would be required, but the difference is apparent here.

Note that there are two minima, which correspond to the H₂ molecule on the bridge site and two H atoms adsorbed in fcc and hcp sites, i.e. the dissociated molecule adsorbed on the surface.

3.4 Questions¶

1. Why do the two images differ so much?
2. What systems can elbow plots be used to describe?
3. What adsorption structures do the two different minima correspond to?

4 Instrinsic reaction coordinates: S_N2 substitution reaction

$\uparrow$

¶

By the end of this tutorial, you will be able to:

• follow the intrinsic reaction coordinate (IRC) forward to the product and backward to the reactant
• visualize the reaction profile and the reactants/ products

4.1 Task¶

Follow the intrinsic reaction coordinate for S_N2 Cl^- substitution of chloromethane from the transition state in a 12×12×12 ų supercell using a machine-learned force field.

A final transition state approach is to calculate the intrinsic reaction coordinate (IRC) [J. Phys. Chem. A 106 (2002) 165]. The steepest descent path from the transition state in mass-weighted Cartesian coordinates defines the IRC, producing the reaction profile along this coordinate. Beginning from the transition state TS, the structure is propagated along the TS' imaginary vibrational mode $u_{\xi}$ forward toward the product P and backward to the reactant R using a damped-velocity Verlet (DVV) algorithm. The structure of the TS must be already well-optimized, e.g. using IDM.

The S_N2 substitution reaction of CH₃Cl by Cl^- will be followed using IRC. A machine-learned force field (MLFF) will be used to speed up the calculation so that it finishes within a few minutes.

4.2 Input¶

The input files to run this example are prepared at $TUTORIALS/transition_states/e04_IRC.

POSCAR

transition state for the SN2 reaction of CH3Cl with Cl-
  1.00000000000000
   12.0000000000000000    0.0000000000000000    0.0000000000000000
    0.0000000000000000   12.0000000000000000    0.0000000000000000
    0.0000000000000000    0.0000000000000000   12.0000000000000000
  C    H    Cl
    1     3     2
 Direct
 0.5850734246784209  0.5520914206000113  0.6939676081570985 ! coordinates for atom 1
 0.6198799799863561  0.5166828822917992  0.6192173893914401
 0.5120438703668218  0.5155341457726120  0.7311383922253840
 0.6232988535554732  0.6240585575728448  0.7315296645899942
 0.4727762112937258  0.6675904758190063  0.5870243322899322
 0.6973145053782144  0.4366510550400121  0.8010186348593505 ! coordinates for atom N
! unstable direction optimized by the dimer method
 0.50309310E+00 -0.52116977E+00  0.48103627E+00             ! components for atom 1
-0.29990134E-01  0.31072124E-01 -0.24160721E-01
-0.37734828E-01  0.43431417E-01 -0.39422281E-01
-0.41647034E-01  0.41752749E-01 -0.38921427E-01
-0.20205452E+00  0.20762993E+00 -0.17619553E+00
-0.19166655E+00  0.19728353E+00 -0.20233630E+00             ! components for atom N

INCAR

! IRC settings
IBRION = 40                  # invokes the IRC calculation
IRC_DIRECTION = -1           # negative initial direction of movement
IRC_STOP = 20                # terminate when IRC_STOP energies in row increase
IRC_VNORM0 = 0.0005          # affects accuracy, the smaller the better but the more
                             # ionic steps needed
NSW = 5000                   # maximal number of steps

! Ab initio settings
ISIF = 2
NELECT = 22                  # the overall charge of this system is -1 hence the number
                             # of electrons is increased (not really needed when
                             # the calculation is done at the MLFF level)

! ML settings
ML_LMLFF = .TRUE.            # invoke MLFF
ML_MODE = run                # sets the mode of MLFF to production

KPOINTS

Automatic
0
Gamma
1  1  1
0. 0. 0.

POTCAR

 Pseudopotentials of C, H, and Cl.

The POSCAR file is a little different to what you have seen previously. The coordinates are not given in cartesian coordinates (xyz) but in direct (in terms of the lattice vectors). Additionally, there are more lines after the ionic coordinates, which define the unstable direction for the IRC to follow obtained from an IDM calculation.

Check the tags set in the INCAR file! Most of these have been defined in the previous exercises. We will go through the ones specific to here. There are two ML tags, ML_LMLFF and ML_MODE, which use a pre-generated machine-learned force field (MLFF).

The IRC calculation is defined using IBRION = 40. The direction of the initial displacement of the transition state structure along the unstable direction is set using IRC_DIRECTION; 1 sets it parallel to the unstable direction toward the product, while -1 sets it antiparallel to the unstable direction and toward the reactant. IRC_STOP = 20 sets the number of steps of relatively constant energy before the IRC calculation terminates and IRC_VNORM0 sets the constant velocity vector.

The only remaining tag is NELECT = 22, which adds an extra electron to the system (the chloride ion holds a negative charge).

A Γ-point only $\mathbf{k}$-mesh is used, sufficient for a gas-phase reaction.

4.3 Calculation¶

Open a terminal and navigate to this example's directory by entering the following:

cd $TUTORIALS/transition_states/e04_*

Congratulations, you've completed your first IRC calculation! This IRC calculation took only a few seconds to do ~450 steps. The machine-learning force field is the reason that this was affordable. It would have taken several hours if this had been done fully ab initio.

Now let us analyze the data. To extract it, simply use the ircShift.sh script with the following commands:

ircShift.sh extracts the energy and the IRC coordinate at each step:

 DAMPED VELOCITY VERLET ALGORITHM:
  -----------------------------------------------------------
    IRC (A):  -0.77141163 E(eV): -.28118857E+02

What do you expect the shape of the reaction pathway to be?

Let us visualize the data to find out:

The graph shows a neatly symmetric curve. The ion attacks the carbon directly opposite the chlorine atom, then inverts the methyl group, to eject the original Cl as Cl^-. Overall, this reaction conserves the reactants, so the input equals the output, thus creating a symmetric curve.

To get a better sense of the reaction, you may visualize the reactant and product structures using the following script:

You can view the animated the reaction using the anim.xyz file in JMol. A recording of anim.xyz is given below:

You can see that the entire process progresses smoothly from reactant to product. This reaction is a simple, symmetric one. Later, the techniques learned here will be applied to a more complex and interesting reaction.

4.4 Questions¶

1. What tag must you change to get the forward and backward transitions?
2. How long would an equivalent ab initio IRC calculation take?
3. Why is the IRC curve symmetric in this reaction?

Good job! You have finished Part 1!¶

Go to Top $\uparrow$