讓VASP實現定壓計算的python腳本

02-01

最近在研究指定的應力張量下晶體缺陷的行為。本來想用VASP來進行相應的DFT模擬。但查了半天，發現VASP並不支持指定應力的計算（除非是isotropic的應力）。具體來說，你可以輸入一個應變，VASP可以算出對應的應力。但是，它不能自動調整應變以獲得你需要的應力。當然，你可以手動的調整模型來實現這個繁瑣的步驟，我有個師兄就是這麼乾的（心疼他一秒...）。不過，秉承著能自動化就就絕不手動乾的精神，我還是寫了個python腳本來實現這個過程（請無視我的渣渣python水平）。該腳本會根據廣義胡克定律自動修改應變，然後迭代地調用VASP直到應力張量收斂到你所期望的值。

我已經用這個代碼對W，Mo，Cr，alpha-Fe體系做了一些模擬，到目前為止表現還不錯。儘管如此，這始終是個beta版本，你可能會遇到一些我沒有考慮到的bug。若然，你可以email至 jie.hou2@mail.mcgill.ca向我報告bug，或者自己動手修復bug。任何人都可以自由的使用，修改和傳播這段代碼。

要使用這個代碼，你需要安裝一個python解釋器（我用的是python 2.7.3）。詳細的使用說明寫在代碼頭部的注釋中。代碼不長，約100行左右，用起來應該不難。

2018/01/11

#!/usr/bin/pythonn#PBS -l nodes=6:ppn=12n#PBS -l walltime=3:00:00n#PBS -Vn#PBS -A emm-484-adn#PBS -N W-120H2nnInstructions:n1. Input the pressure (stress) tensor in the script(see below)n2. Input estimated elastic modulus in the scriptn3. Set ISIF=2 in the INCARn4. Copy POSCAR to poscar.0n5. Configurate mpiexec commands in the script according to your VASP directoryn6. cd to job directory and use qsub fixpressure.py to submit your jobnnSuggestions for accurate force calculation:nUse a high ECUT. See VASP manual for details.nUse PREC=High or Accurate, and LREAL=AutonnSuggestions for better convergency:nUse LWAVE=.T. and default ISTART and ICHARGnnGeneral algorithm:n1. run vasp with ISIF=2 to get current pressure tensorn2. using generalized Hookes law and elastic modulus you input, estimate the POSCAR modification required to get target pressuren3. repeat 1-2 until the current pressure - target pressure is within convergency criterianShould you find any bug, please contact jie.hou2@mail.mcgill.cannnimport sys,osnimport subprocessnimport shutilnimport stringnimport numpy as npnimport linecachenn#######################You need to configurate the following variables##############nSetpress= np.array([120,0,0,0,0,0])n#Set the pressure tensor (xx yy zz xy yz zx) in kB. Note pressure=-stress in VASPnpresscirt= 0.1 n#convergency criteria for pressure, unit in kBnE=2790 n#Youngs modulus in kBnv=0.21 n#Possion rationG=E/(2-2*v)n#Shear modulus in kBnimax=30n#maximum iteration cyclesnmpiexec=mpiexec -n `cat $PBS_NODEFILE | wc -l` /gs/project/emm-484-aa/jiehou/software/vasp54/vasp_std > results.txtn#set your VASP executable directory heren##################################################################################nnnnn########do not alter the following codes unless you know what you are doing#######nos.chdir(os.getenv(PBS_O_WORKDIR, ))nshutil.copy (poscar.0,POSCAR)nsubprocess.call("echo starting calculation > pressure.all", shell=True)nsubprocess.call(mpiexec, shell=True)n#initial run to calculate pressure tensornsubprocess.call("cat OSZICAR > oszicar.all", shell=True)nsubprocess.call("cat OUTCAR > outcar.all", shell=True)nsubprocess.call("grep Total CPU time used (sec): OUTCAR | tail -n 1 > end.txt", shell=True)nniteration=0nwhile iteration<=imax:n iteration=iteration+1n P_coeff=1-iteration*1.0/imaxn #graduately reduce P_coeff to avoid flucations in convergencyn print iteration=,iterationn subprocess.call("grep Total CPU time used (sec): OUTCAR | tail -n 1 > end.txt", shell=True)n notempty=os.path.getsize(end.txt)n os.remove(end.txt)n #scan if the previous calculation is completedn if notempty:n subprocess.call(" grep in kB OUTCAR | tail -n 1 > pressure.txt", shell=True)n subprocess.call(" grep in kB OUTCAR | tail -n 1 >> pressure.all", shell=True)n for line in open("pressure.txt"): n press=np.array([ float(x) for x in line.split()[2:]])n print pressure=,pressn os.remove(pressure.txt)n #read pressure from OUTCARn addpress=Setpress-press[0:6]n print adding pressure:,addpressn #calculate the additional pressure needed to achieve Setpressn abs_P=max([abs(y) for y in addpress])n if (abs_P<presscirt):n subprocess.call("echo aborting calculation for the convergence is reached >> pressure.all", shell=True)n breakn #stop if the additional pressure is too smalln else:n #copy CONTCAR to POSCAR, n shutil.copy (CONTCAR,POSCAR)n POS=linecache.getlines(POSCAR)n M=np.zeros((3,3))n addM=np.zeros((3,3))n M[0,:]=np.array([float(x) for x in POS[2].split()])n M[1,:]=np.array([float(x) for x in POS[3].split()])n M[2,:]=np.array([float(x) for x in POS[4].split()])n #read current POSCAR matrixn addM[0,0]=(addpress[0]-v*(addpress[1]+addpress[2]))/(-E)n addM[1,1]=(addpress[1]-v*(addpress[0]+addpress[2]))/(-E)n addM[2,2]=(addpress[2]-v*(addpress[1]+addpress[0]))/(-E)n addM[0,1]=addM[1,0]=addpress[3]/(-G)/2n addM[1,2]=addM[2,1]=addpress[4]/(-G)/2n addM[0,2]=addM[2,0]=addpress[5]/(-G)/2n #calculate how much strain we need to add according to addpressn #beware that VASP defines compressive pressure as positive, hence /E becomes /(-E) and /G becomes /(-G)n addM=np.diag([1,1,1])+addM*P_coeffn M=np.dot(M,addM)n print adjusting matrix to:n print Mn np.savetxt("M.txt",M)n linecache.clearcache()n addPOS=linecache.getlines(M.txt)n POS[2:5]=addPOSn subprocess.call("rm M.txt", shell=True)n fo = open("POSCAR", "w+")n fo.writelines(POS)n fo.close()n posname=poscar.+str(iteration)n #shutil.copy (POSCAR,posname)nn #run vasp to calculate pressure for this POSCARn subprocess.call(mpiexec, shell=True)n subprocess.call("cat OSZICAR >> oszicar.all", shell=True)n #run VASP program, configure the shell command according to your linux environmentn else:n print calculation abortedn #previous vasp simulation is not completed, calculation abortnsubprocess.call("rm WAVECAR", shell=True)nsubprocess.call("rm outcar.all", shell=True)nsubprocess.call("rm PROCAR", shell=True)nsubprocess.call("rm DOSCAR", shell=True)nsubprocess.call("rm EIGENVAL", shell=True)n