Climbing Image Nudged Elastic Band in VASP
Vanilla VASP can only do Nudged Elastic Band (NEB) without the climbing image option. There exists a package from Henkelmann's group at http://theory.cm.utexas.edu/vtsttools/ that enables climbing images and provides other optimizers for your NEB calculations.
To turn on the climbing image option, one simply sets the tag LCLIMB = .TRUE. in one's INCAR file. Alternatively ASE does also understand it and in the calculator one sets lclimb = True as follows:
calculator = VASP(..., lclimb=True)
N.B. when the climbing image option is turned on, the relaxation algorithm used in the simulation must not expect conservative forces. The different optimizers you can use from Henkelmann's group and how to load them is shown on http://theory.cm.utexas.edu/vtsttools/optimizers.html
Using ASE to interpolate between initial and final states
In ASE one has the option to interpolate linearly between the final and initial states, also one can interpolate using the image dependent pair potential(IDPP) scheme published by Smidstrup et al. http://dx.doi.org/10.1063/1.4878664. The IDPP interpolation is based on keeping the pairwise distances as constant as possible, and this can help errors where the atoms come too close and also save some calculation time. The idpp interpolation is a fairly recent addition to ASE(>= version 3.9) and might not be available in the version you have loaded.
A full example on how to generate the initial path with linear and idpp interpolation is shown below. The example is based on the user having the initial and final configuration in a .traj file; here the initial state is methane dissociation over a Pd(111) surface. The example is run by the command python neb_init.py and then calling VASP from a submit script afterwards.
# neb_init.py from ase import Atoms, Atom from ase.calculators.vasp import Vasp from ase.io.vasp import write_vasp import numpy as np from ase.constraints import FixAtoms from ase.neb import NEB from ase.io import read,write import os # Constants a0 = 3.95612 k0 = 6 Nim = 5 #Number of intermediate images ini_fn = "initial.traj" fin_fn = "final.traj" # Read in structures ini = read(ini_fn) fin = read(fin_fn) # Fix the Pd atoms exept two layers: del ini.constraints pos_fix = float(ini.positions[8][2]) cs = FixAtoms(mask=[atom.position[2]<pos_fix-1. for atom in ini ] ) ini.set_constraint(cs) # Gather images for NEB images = [ini] images += [ini.copy() for i in range(Nim)] images += [fin] calc = Vasp(prec='Normal',xc='PBE', encut = 450.0, #PW cut-off ispin=1, #spin-polarization ibrion=3, #Relaxation RM-DIIS nsw = 500, #Max. no of relaxation steps. kpts = [k0,k0,1], sigma=0.05, potim=0, iopt = 7, # Use FIRE to optimize it maxmove = 0.05, # Max move atoms 0.05 Ang per iteration isif=0, ediffg=-1E-1, # Modify after climbing lreal = "Auto", npar = 4, kpar = 4, isym=False, #Do not use symmetries, can help NEB problems lwave=False, lcharg=False, lclimb=False, # Do not climb yet. images=Nim, spring=-5.0 ) # Interpolate. neb = NEB(images) neb.interpolate("idpp") #intepolate by IDPP #neb.interpolate() # - remove comment for linearinterpolation neb.images[0].set_calculator(calc) print "Initializing..." calc.initialize(neb.images[0]) print "Done!.." calc.write_incar(neb.images[0]) calc.write_potcar() calc.write_kpoints() #Now for each image, create its directory, cd into that #and create the POSCAR: for i in range(len(images)): os.system("mkdir %02i" % i) #makes folder 00, 01, 02 write("%02i/POSCAR" % i,neb.images[i],format="vasp") # Write the poscar
The submit script must call your custom VASP installation compiled with the VTST tools as the following example shows:
#!/usr/bin/env bash #SBATCH -N 5 ... export NCORES=$(($SLURM_JOB_NUM_NODES*16)) module purge module load mkl/11.1/100 module load openmpi/1.5.4 mpirun /c3se/users/mikjorge/Glenn/software/custom_vasp/vasp
Interpolating selectively and selectively freeing different atoms in the NEB
Sometimes it is feasible to only interpolate the some atoms, but despite that let some of the non-interpolated be free in the NEB calculations. This is done easily in ASE by first fixing all atoms that you would not like to interpolate. In other words: ASE does not interpolate atoms that are fixed, and we can just free the ones afterwards. The following example shows how to achieve the goal:
# Fix all Pd atoms to not interpolate the surface atoms : del ini.constraints cs = FixAtoms(mask=[atom.symbol=="Pd" for atom in ini ] ) ini.set_constraint(cs) images = [ini.copy for i in range(Nim)] ... neb.interpolate("idpp")
Then after interpolation and before writing the INCAR and POSCARS free the atoms that are free in the NEB calculations:
# Fix atoms that are lower down than atom number 8 in z direction pos_fix = float(ini.positions[8][2]) cs2 = FixAtoms(mask=[atom.position[2]<pos_fix-1. for atom in ini ] ) for i in range(len(images)): del images[i].constraints images[i].set_constraint(cs2) # Now write the VASP files...
General considerations
A good way to begin a TST search with NEB is to set up the atoms in their prefered positions for the initial and final state. Then converging the simulation a little with with sufficiently many images to sample the path. Then you might note where the interesting part of the barrier is located and zoom in on that, so to say. Here one can typically decrease the number of images and perhaps climb after a few iterations.
The optimizers in vanilla VASP are not all force based and using climbing images introduces a dissipative force. Thus it is good to use the Henkelmann group's optimizers that are turned on setting POTIM = 0 and IBRION=3 in your INCAR file. Then the specific Henkelmann VTST optimizer is chosen with the IOPT = X keyword. You can read about the optimizers at http://theory.cm.utexas.edu/vtsttools/optimizers.html . Also remember that you might need to lower the MAXMOVE parameter since it allows for the atoms to move 0.2 Å between each iteration, which is too high sometimes.