Tutorial on using nudged elastic band method

This tutorial is based on VASP 5.4.4 and the VTST package developed and maintained by the Henkelman research group at UT Austin

Installing VASP with VTST package

The general information of how to install VASP with VTST can be found at the VTST website. Here is a demonstration of installation VASP 5.4.4. + VTST on NREL EAGLE HPC system:

1. Load modules needed for VASP

1) gcc/8.4.0  2) comp-intel/2020.1.217 3) intel-mpi/2020.1.217              
4) fftw/3.3.8/intel-impi 5) mkl/2020.1.217

2. Modify main.F file within the src/ folder of VASP 5.4.4:

Find and Replace the following lines

CALL CHAIN_FORCE(T_INFO%NIONS,DYN%POSION,TOTEN,TIFOR, &
    LATT_CUR%A,LATT_CUR%B,IO%IU6)

with

CALL CHAIN_FORCE(T_INFO%NIONS,DYN%POSION,TOTEN,TIFOR, &
    TSIF,LATT_CUR%A,LAT  T_CUR%B,IO%IU6)

3. Apply patches (Need to make sure the patch file is located at the correct folder by checking the first few lines of the patch file). The patch is included in the VTST package

patch -p0 < vasp-5.4.4-mpmd.patch

4. Copy source file from vtstcode5 folders. For example, the version 182 has the following source files:

bbm.F, bfgs.F, chain.F, dynamic.F, fire.F, lanczos.F, neb.F, qm.F, 
cg.F, dimer.F, dynmat.F, instanton.F, lbfgs.F, opt.F, sd.F.

5. Modify the compiling order within the hidden .objects file within the folder of src/. Note that the added .o file should be placed before chain.o

 ....
 hamil_rot.o \
 bfgs.o dynmat.o instanton.o lbfgs.o sd.o cg.o dimer.o bbm.o \
 fire.o lanczos.o neb.o qm.o opt.o \
 chain.o \
 dyna.o \ 
 ....

6. compile VASP as it is done typically (e.g. make std)

There are significant differences in the installing procedure for VASP 6, see details on the VTST website.

High-level idea of nudged elastic band method

Only very high-level picture is provided here, further details can be found at the original papers by Henkelman and Jónsson[1][2].

Schematic minimum energy path (red dashed curve) between two local minima (A and B sites) on a 2D potential energy surface

The general goal of NEB calculation is to find the minimal energy path (MEP) between two user-defined local minima on the potential energy surface (PES) of a given system (see figure below). The potential energy surface is a function of configuration coordinate and the relation can be described by semi-empirical force field or electronic-structure method like DFT under the Born-Oppenheimer approximation. Using the 2 dimensional PES as an example, one naive way to find the MEP is to sample a set of evenly distributed 2D grid points and get the PES for each point. In principle, we can find the MEP when the grid is dense enough.

To find the MEP more efficiently, nudged elastic band method (NEB) is commonly used. The method start with an initial discrete pathway constructed by linear interpolation between two end points (A and B sites in the figure). Such discrete linear pathway is generally a good initial guess but by no means has the MEP. In order to find MEP, ionic relaxation is needed but these discrete points will relax into either site A or site B, depending on its initial position on the PES. To make these points "hung" along the path between two local minima, springs are introduced to connect adjacent discrete points and thus the method has the name of "elastic band". The concept of "nudging" is to separate forces along and perpendicular to the tangence of the elastic band. The calculation only take into account spring forces along the tangence of the band and the true force (e.g. forces from DFT) perpendicular to the tangence of the band. As a result, the distance between discrete points will remain roughly the same while the path is nudged toward minimum energy path.

Building on top of the NEB method, the climbing image NEB (CI-NEB) method was introduced to find the saddle point with higher energy resolution. The saddle point is the maximum along the minimum energy path while the minimum along other orthogonal directions. The idea of climbing image is to make the discrete point or the intermediate image that has the higher energy climb to the maximum. In order to achieve this numerically, for this specific image, the sign of the forces along the elastic band is flipped such that the saddle point becomes a local minimum.

Furthermore, in addition to atomic degree of freedom (DOF), and unit cell DOF can be also included and such approach is referred to as the generalized solid-state NEB (SS-NEB). This is a powerful method to find solid-solid phase transformation, ferroelectric switching, etc. Generally speaking, much higher number of images may be needed for SS-NEB since the minimum energy path could be more complicated than phenomena that can already captured by NEB. (SS-NEB generally involves more DOF)

VASP input for NEB calculations

Here's an overview of the steps to perform NEB calculations and detail

1. Prepare end point/images. Using vacancy-mediated diffusion as an example. The end images are supercell with vacancy on two different but adjacent sites. The end images should be relaxed.

2. Prepare INCAR and initial intermediate images specific for NEB calculations.

3. Run NEB first to have a clear intermediate image with highest energy difference from end points. Then run CI-NEB to refine the energy barrier

4. Extract and analyze the data

If multiple local minima are found along the MEP, the MEP can be broken into smaller pieces and calculated separately.

Preparing intermediate images

The script (nebmake.pl) provided within the VTST package is handy in terms of making linear interpolation between end images.

nebmake.pl POSCAR_i POSCAR_f 5

where i and f indicates initial and final images and 5 is the user-defined number of intermediate images to create. Odd number is generally preferred for systems that are expected to have single peak along the minimum energy path. Once this is done, 5+2 folders (+2 are for end images) will be created with the associated linear interpolated images inside. Note: intermediate images are the configuration coordinate, i.e. positions of atoms, along the path and are embodied by the POSCAR files for VASP program.

To use the script properly, the order of atoms within the VASP POSCAR should be the same for end images. Otherwise, weird initial path will be created and lead to weird results.

VASP+VTST INCAR specific for NEB calculations

Only flags different from normal VASP relaxation are highlighted below

ISIF = 3  (need for SS-NEB; ISF=2 for NEB)
NELMIN = 5 (avoid ionic stpes with only two electronic scf steps)
ICHAIN =0  (for VTST package, this is the flag for NEB)
IMAGES = 5 (the number of intermediate images, end images not included)
SPRINGS = -5  (Spring constant for the elastic band. -5 generally works well)
LCLIMB = .FALSE. (FALSE is default value. Only turn on this flag when there is a clear maximum)
IOPT = 3  (3 for SS-NEB; 3→ 2 → 1 for CI-NEB)
IBRION  = 3 (required for VTST)
POTIM = 0 (required for VTST)
LNEBCELL = .TRUE. (True for SS-NEB; otherwise normal NEB)

Generally, higher the number of images, better sampling of the minimum energy path. Simple MEPs only require 3 images while complicated ones need 10+. Unable to reach force convergence criteria is usually an indication to increase the number of intermediate images.

Restart the simulation

The neb calculation generally requires multiple restart to converge the results. To restart the simulation, like other VASP relaxation calculations, we need to make CONTCAR the new POSCAR. Below is an example bash script to do so, assuming there are 13 intermediate images:

for i in `seq -w 1 13`; do mv $i/CONTCAR $i/POSCAR;done

The scripts provided by the VTST package are very useful! The following are few examples

Data extraction and analysis

nebbarrier.pl   (create neb.dat file which has the information on relative energy, distance between image and forces along the bands)
nebspline.pl (forced based cubic spline interpolation for the MEP)

There is a different set of scripts for solid-state NEB (e.g nebefs.pl and nebsplines.pl)