Calculation of elastic constants

From EVOCD
Jump to: navigation, search

< Homework 1

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.

E-e relation.png

(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 by
chmod u+x <script name>.py
Also 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

Personal tools
Namespaces

Variants
Actions
home
Materials
Material Models
Design
Resources
Projects
Education
Toolbox