Calculation of elastic constants
From EVOCD
Author(s): Laalitha Liyanage
Contents |
License
By using the codes provided here you accept the the Mississippi State University's license agreement. Please read the agreement carefully before usage.
Abstract
This page gives the tools necessary to calculate the elastic constants of fcc crystal structures using VASP density functional theory package. Chapter 9 of "Intermetallic compounds principles and practice vol.1: principles" by Mehl et al. [1] explains the method to calculate elastic properties of metallic compounds through first principles calculation. The python script presented can calculate elastic constants of any cubic crystal with modification to crystal structure section of the code.
Elastic constant calculation script
The following script can be used to calculate the elastic constants of fcc crystal structure. The script uses the Energy-strain relationship [2],[3] as given by the following equation to calculate the elastic constants.
(Here E0 is the equilibrium energy, V is the volume of the undistorted lattice, P(V) is the pressure at that volume, ΔV is the change in the volume of the lattice due to strain, C[i,j] is the stiffness tensor, e[i], e[j] are strain tensors written in vector form, and O[e^3] are neglected higher order terms.).
It requires the equilibrium lattice parameter, type of calculation (relaxed/unrelaxed), percentage strain and the number of steps to include. VASP will be executed with in the script, therefore the following string in the script which is under "Default setting" has to be updated.executable = "mpirun -np 8 ~/bin/vasp_5212_talon"The script requires 3 types of INCAR files named INCAR.isif3, INCAR.isif2, INCAR.static, a KPOINTS file, a POTCAR file, and a cleanvaspfiles file. The converged k-point grid from the Energy-Volume curve calculation can be used here. The 3 INCAR files are given below after the script.
Running the script
To run the script have all three INCAR files and the script saved in the same directory. Give executable permissions to the python script bychmod u+x <script name>.pyAlso have KPOINTS file with converged k-point grid obtained from Energy-Volume curve results and the POTCAR file (LDA, GGA, or PBE) in the same directory. You will have to submit the script to a cluster (Raptor), therefore use a pbs command script. The command to execute the script is
./<script name>.py 4.07 relaxed 1 20
Output of script
The script will write the elastic constants to a file called responses.txt. It will also generate three plots relating to the VASP calculated Energy-Strain curve and its fitted plot. These are in gif files named *_fit.gif where * relates to c11,c12 and c44. The output of VASP calculations are written to *report.vasp files. If any error is encountered please refer to the *report.vasp files.
The gif files can be visualized using the following command
display name.gif
The Script
#!/usr/local/bin/python
#Purpose : Calculate C11,C12 and C44 of fcc crystal structure
#Author : Laalitha Liyanage laalitha@cavs.msstate.edu
import os
import re
import sys
import math
import numpy as np
usage = """
Usage: fcc_elconst.py a type perc nt
a - fcc lattice parameter
type - relaxed/unrelaxed
perc - maximum percentage strain
nt - number of points
E.g. : fcc_elconst.py 4.05 unrelaxed 2 20
"""
def gen_gpfile_for_parabola_fit(a,b,c,prefix):
#Generate input script for GNUplot for fitting
os.system("rm -f fit.log")
fn = "%s_fit.gp"%prefix
fout = open(fn,"w")
fout.write('set term gif\n')
fout.write('set output "%s_fit.gif"\n'%prefix)
fout.write("f(x) = a*x**2 + b*x+ c\n")
b = 0.0
fout.write("a = %f; b = %f; c = %f\n"%(a,b,c))
fout.write("FIT_LIMIT = 1e-9\n")
if b == 0.0:
fout.write("fit f(x) '%s_summary' via a,c\n"%prefix)
else:
fout.write("fit f(x) '%s_summary' via a,b,c\n"%prefix)
fout.write('fx_title = sprintf("(%f)*x**2+(%f)*x+(%f)",a,b,c)\n')
fout.write("plot '%s_summary' t 'MEAM' with linespoints 1 5,f(x) t fx_title"%prefix)
fout.close()
return fn
def get_fitted_param(a,b,c,fn):
#Extract fitted parameters from GNUplot output
fin = file(fn,"r")
while True:
line = fin.readline()
if line == "": sys.exit(123)
if re.search("^Final set of parameters",line): break
for i in range(2): line = fin.readline()
for i in range(2):
line = fin.readline().split()
if len(line)>0 and line[0] == "a": a = float(line[2])
elif len(line)>0 and line[0] == "b": b = float(line[2])
elif len(line)>0 and line[0] == "c": c = float(line[2])
fin.close()
return a,b,c
def gen_poscar(box,strain_tensor):
#Generate crystal structure in VASP format.
ax,ay,az = box[0],box[1],box[2]
e1,e2,e3,e4,e5,e6 = strain_tensor
ax_p = [(1.0+e1)*ax,0.5*e6*ax,0.5*e5*ax]
ay_p = [0.5*e6*ay,(1.0+e2)*ay,0.5*e4*ay]
az_p = [0.5*e5*az,0.5*e4*az,(1.0+e3)*az]
fout = open("POSCAR","w")
fout.write("Al\n")
fout.write("1.0\n")
fout.write(" %22.16f %22.16f %22.16f\n"%(ax_p[0],ax_p[1],ax_p[2]))
fout.write(" %22.16f %22.16f %22.16f\n"%(ay_p[0],ay_p[1],ay_p[2]))
fout.write(" %22.16f %22.16f %22.16f\n"%(az_p[0],az_p[1],az_p[2]))
fout.write("4\n")
fout.write("Direct\n")
fout.write("0.0 0.0 0.0\n")
fout.write("0.0 0.5 0.5\n")
fout.write("0.5 0.0 0.5\n")
fout.write("0.5 0.5 0.0\n")
def get_E0(a):
# Get equilibrium energy, volume and lattice parameters through geometric optimization
gen_poscar([a,a,a],strain_tensor=[0,0,0,0,0,0])
os.system("cp INCAR.isif3 INCAR")
os.system("%s > E0.report.vasp"%(executable))
os.system("cp CONTCAR POSCAR")
os.system("./cleanvaspfiles")
os.system("cp INCAR.isif3 INCAR")
os.system("%s >> E0.report.vasp"%(executable))
os.system("cp CONTCAR POSCAR")
os.system("./cleanvaspfiles")
os.system("cp INCAR.static INCAR")
os.system("%s >> E0.report.vasp"%(executable))
E0 = get_data("tail -n1 OSZICAR |awk '{print $5}'")
V0 = get_data("grep 'volume of cell' OUTCAR |tail -1|awk '{print $5}'")
A0 = get_data("head -n3 CONTCAR |tail -n1|awk '{print $1}'")
os.system("mv OSZICAR OSZICAR.E0")
os.system("mv CONTCAR CONTCAR.E0")
os.system("mv OUTCAR OUTCAR.E0")
return(E0, V0, A0)
def get_data(expr):
#Grab data from file
os.system("%s>eout"%(expr))
fin = file("eout",'r')
line = fin.readline().split()
return (float(line[0]))
def get_unrelaxed_eng(t):
#execute VASP
os.system("echo 't=%f' >> report.vasp "%(t))
os.system("cp INCAR.static INCAR")
os.system("%s >> report.vasp"%(executable))
return(get_data("tail -n1 OSZICAR |awk '{print $5}'"))
def get_relaxed_eng(t):
#execute VASP
os.system("echo 't=%f' >> report.vasp "%(t))
os.system("cp INCAR.isif2 INCAR")
os.system("%s >> report.vasp"%(executable))
os.system("cp CONTCAR POSCAR")
os.system("./cleanvaspfiles")
os.system("cp INCAR.static INCAR")
os.system("%s >> report.vasp"%(executable))
return(get_data("tail -n1 OSZICAR |awk '{print $5}'"))
def get_cxx(cname,a0,a,b,c):
#Calculate elastic constants
energy = []
os.system("echo %s > report.vasp"%(cname))
for i in range(nt+1):
t = t0+i*dt
#Generate lattice according to strain tensor
ax,ay,az = a0,a0,a0
if cname == 'c11':
gen_poscar([ax,ay,az],strain_tensor=[t,0,0,0,0,0])
elif cname == 'c12':
gen_poscar([ax,ay,az],strain_tensor=[t,-t,0,0,0,0])
elif cname == 'c44':
gen_poscar([ax,ay,az],strain_tensor=[0,0,0,2*t,0,0])
#Run VASP and get energy
if type == 'unrelaxed':
energy.append(get_unrelaxed_eng(t))
elif type == 'relaxed':
energy.append(get_relaxed_eng(t))
fout = open("%s_summary"%cname,"w")
for i in range(nt+1):
fout.write("%22.16f %22.16f\n"%(t0+i*dt,energy[i]))
fout.close()
gpfile = gen_gpfile_for_parabola_fit(a*V0,1.0,E0,cname)
os.system("gnuplot %s > gnuplot.report 2>&1"%gpfile)
a,b,c = get_fitted_param(a*V0,b,c,"fit.log")
os.system("mv fit.log %s_fit.log"%cname)
os.system("mv report.vasp %s_report.vasp"%cname)
return a/V0
# Main Program
#==========================================================================
# Default setting
#=====================================================================================
unit_conversion1 = 160.217646 # 1 eV/A^3 = 1.60217646E-19 / 1E-30 Pa = 160.217646 GPa
#PLEASE MODIFY the following command to suit your system
executable = "mpirun -np 8 ~/bin/vasp_5212_talon"
if len(sys.argv) > 4:
a = float(sys.argv[1])
type = sys.argv[2]
perc = float(sys.argv[3])
nt = int(sys.argv[4])
t0,t1 = -perc/100,perc/10
dt = (t1 - t0)/nt
E0,V0,a0 = get_E0(a)
c11 = get_cxx("c11",a0,0.5,1.0,E0)*2*unit_conversion1
k12 = get_cxx("c12",a0,0.5,1.0,E0)*2*unit_conversion1
c22 = c11
c12 = (c11+c22-k12)/2
c44 = get_cxx("c44",a0,2.0,1.0,E0)/2*unit_conversion1
fout = file("responses.txt",'w')
S = "C11\t%f\t(107)\nC12\t%f\t(61)\nc44\t%f\t(28)\n"%(c11,c12,c44)
fout.write(S)
fout.close()
else:
print "\n\tERROR: Wrong number of arguments !!!"
print usage
INCAR files
Copy and save these files with the same file name. These files should be in the same directory as the python script.
INCAR.isif3
LWAVE = .FALSE. LCHARG = .FALSE. PREC = Accurate ISMEAR = 1 ENCUT = 318 LREAL = .FALSE. ISIF = 3 NSW = 100 IBRION =2
INCAR.isif2
LWAVE = .FALSE. LCHARG = .FALSE. PREC = Accurate ISMEAR = 1 ENCUT = 318 LREAL = .FALSE. ISIF = 2 NSW = 100 IBRION =2
INCAR.static
LWAVE = .FALSE. LCHARG = .FALSE. PREC = Accurate ISMEAR = -5 ENCUT = 318 LREAL = .FALSE. #ISIF = 3 #NSW = 100 IBRION =2
