Pylada High-throughput Workflow for SS-NEB Calculations
This tutorial is based on VASP+VTST, vtstscripts, and Pylada
Written by Cheng-Wei Lee (clee2 [at] mines [dot] edu)
This tutorial assumes that users already have the relaxed endpoints for (SS-)NEB calculations, which can be also done in a high-throughput fashion using Pylada scripts like HT_relax.py or any other Python workflow that can automate the VASP calculations. The tutorials cover three topics: (1) modification to existing Pylada package, (2) customized Pylada workflows, and (3) user defined inputs.
Modification to exiting Pylada
Since nudged elastic band (NEB) calculations are performed differently from the typical VASP calculations, i.e., there are multiple working subfolders within a working folder, there are two required modifications to the published Pylada package.
modification 1 applies to "pylada/vasp/functional.py"
At line 865-868 (see blow) of the python script, there is a IF condition block to ensure that the typical VASP runs are executed properly. Since the file structure are different for NEB calculations from typical VASP calculations, we need to comment these lines out as a temporary solution.
865 ## Last yield should be an extraction object.
866 #if not program.success:
867 # raise RuntimeError("Vasp failed to execute correctly.")
868 #return program
modification 2 applies to "pylada/vasp/keywords.py"
For the two blocks that determine ICHARG and ISTART for continuing VASP calculations, the original Pylada package only consider the main working folder and we need to add one of the working subfolders to the checklist:
original ones at line 702 and 796
directories = [outdir]
modified one for both
directories = [outdir, outdir+'/01/']
Customized Pylada workflow
Here's the customized Pylada workflow that takes the number of NEB iterations and the number of intermediate NEB images as input.
The user needs to modify the path to vtstscripts
User-defined inputs
Here's an exemplar Pylada submission script to utilize the customized Pylada workflow above. Pylada will prepare the INCAR, KPOINTS, POTCAR for each provided pair of endpoints and interact with HPC system to submit VASP jobs. User still need to create the initial intermediate images.
Here an exmeplar script based on Pylada and vtstscript to generate initial neb intermediate images
######################################
class SSNEB_Chain(object):
# Defining the folder tree structure
def __init__(self, no_iterations=10, no_images=24):
self.no_iterations = no_iterations
self.fnames = ['iter_%s'%(x+1) for x in range(no_iterations)]
self.no_images = no_images
# Defining the Extraction object
def Extract(self, jobdir):
'''
The success condition for the ss-neb calculation is based on all images
for each image, success = True when a VASP run is finished properly (with the timing information shown at the end of OUTCAR)
'''
from pylada.vasp import MassExtract, Extract
extract = MassExtract(jobdir)
#extract = Extract(jobdir)
success={}
from glob import iglob
if self.no_images > 99: digit = '/???'
else: digit = '/??'
for image in iglob(jobdir+digit):
success[image]=Extract(image).success
#for name in self.fnames:
# success[name]=Extract(jobdir+'/'+name).success
success=all(success.values())
#success = all(extract.success.values())
extract.success=success
return extract
def __call__(self, structure, outdir=None, vasp=None, **kwargs ):
from copy import deepcopy
from os import getcwd
import os
from os.path import exists
from os.path import join
from pylada.misc import RelativePath
from pylada.error import ExternalRunFailed
from pylada.vasp.extract import Extract, MassExtract
from pylada.vasp import Vasp
from pylada.vasp.relax import Relax
# make this function stateless.
structure_ = structure.copy()
outdir = getcwd() if outdir is None else RelativePath(outdir).path
############## NEB Loop ########################
for name in self.fnames:
#print("Stage of "+name)
iternum = int(name.split('_')[-1])
fulldir = join(outdir, name)
# create the folder if it does not exist
if not os.path.exists(fulldir):
os.makedirs(fulldir)
# check if the given iteration finished
if exists(fulldir+'/progress.dat'): continue
# restart from unfinished calculations;
# make sure it is not a empty CONTCAR before mv CONTCAR POSCAR
if exists(outdir+'/01/CONTCAR') and os.path.getsize(outdir+'/01/CONTCAR') > 0:
for i in range(1,self.no_images+1):
num = '%02d'%(i)
os.system('mv %s %s'%(outdir+'/'+num+'/CONTCAR',outdir+'/'+num+'/POSCAR'))
## functional
neb = Vasp(copy=vasp)
neb.add_keyword('IMAGES', self.no_images)
if iternum > 1:
neb.istart = 1 # use previous step's WF; but get new FFT grid
neb.icharg = 0
params = deepcopy(kwargs)
## if this calculation has not been done, run it
output = neb(structure_, outdir=outdir, **params)
#print(name+" finished")
# run nebefs.pl and record volume, energy, force, and stress
os.chdir(outdir)
os.system('/home/clee2/software/vtstscripts-967/nebefs.pl > %s'%(fulldir+'/progress.dat'))
os.system('cp stdout %s'%(fulldir+'/stdout'))
#os.chdir('-')
# move the CONTCAR to POSCAR and cp the new POSCAR to fulldir,i.e., the progress folders
for i in range(1,self.no_images+1):
num = '%02d'%(i)
os.system('mv %s %s'%(outdir+'/'+num+'/CONTCAR',outdir+'/'+num+'/POSCAR'))
os.system('cp %s %s'%(outdir+'/'+num+'/POSCAR',fulldir+'/POSCAR_'+num))
if os.system('grep "reached required accuracy" %s'%(outdir+'/stdout')) == 0:
break
else:
# remove the OUTCAR in 01/ such that extract.success remains False before SS-NEB converges
os.system('/home/clee2/software/vtstscripts-967/nebbarrier.pl')
os.system('rm %s'%(outdir+'/01/OUTCAR'))
#if not output.success:
# raise ExternalRunFailed(name+" calculation did not succeed.")
#########################
return self.Extract(outdir)
"custom_chain_SSNEB.py" is the name of the workflow above
import numpy as np
import os
from pylada.crystal import read
from glob import iglob
from pylada.vasp import Extract
# loop through all compounds
# choose an integer number for the number of intermediate images,
numI = 39 #
enforce = False # determine whether to create new POSCARs and overwrite existing poscars
# get current working directory
pathcwd = os.getcwd()
for names in iglob(pathcwd+'/end_images/AlHfN_p083_random*/'):
strc_p = names+'POSCAR_p/CONTCAR'
strc_a = names+'POSCAR_a/CONTCAR'
if Extract(names+'POSCAR_p/').success == True and Extract(names+'POSCAR_a/').success == True:
print('Relaxation of '+names+' finished')
strc = read.poscar(strc_p)
name = names.replace('/end_images/','/SS-NEB_alloy/')
if not os.path.exists(name):
os.makedirs(name)
print("Create new folder "+name)
os.chdir(name)
os.system('/home/clee2/software/vtstscripts-967/nebmake.pl %s %s %s'%(strc_p,strc_a,numI)) # 0.5 is that half of the atoms is cation ;
os.chdir(pathcwd)
elif enforce == True:
os.chdir(name)
os.system('/home/clee2/software/vtstscripts-967/nebmake.pl %s %s %s'%(strc_p,strc_a,numI))
os.chdir(pathcwd)
# copy OUTCAR to first and last folder
first = name+'00/'
last = name+ "%02d/" % (int(numI+1))
os.system('cp %s %s'%(names+'POSCAR_p/OUTCAR',first))
os.system('cp %s %s'%(names+'POSCAR_a/OUTCAR',last))
