Theory

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.

General

Experiment

Theory

PmWiki

pmwiki.org

edit SideBar

Blix theme adapted by David Gilbert, powered by PmWiki