ICME-QM1
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.
Overview
The generalized stacking fault curve represents the energy dependency of a crystal which is rigidly sheared along a plane. Simply it is the Energy-displacement curve as a perfect crystal is sheared.It gives information about a crystals shear properties such as shear strength.
Using the following scripts the generalized stacking fault energy (GSFE) curve for Al (111) surface can be calculated using Density Functional Theory. Test for k-point convergence and compare to values calculated by C. Brandl, P. M. Derlet, and H. Van Swygenhoven in PRB.
Fcc (111) surface generation script
Save the following script as <name>.py and give permissions as executable(chmod u+x <name>.py). It can be executed as
<name>.py <a> <nx> <ny> <nz> <t> <p/d>
The arguments (input parameters) are equilibrium lattice constant (a), periodicity in x,y and z (nx,ny,nz), displacement of top half of the simulation box (t), and output file type (p/d).
#!/usr/bin/python
# Purpose: generate fcc(111) surface
usage="""
Usage:fcc_111_gsfe_curve.py a nx ny nz t p/d
a - equilibrium lattice parameter
nx,ny,nz - periodicity in x,y and z
t - displacement of top block of atoms
p/d - output filetype p(POSCAR)/d(LAMMPS datafile)
"""
import os
import re
import sys
import math
# Default setting
a = 4.0
nx = 1
ny = 4
nz = 1
vaccum = False
def get_fcc(a,t,nx=1,ny=1,nz=1):
""" fcc burger vector = 1/3*[112] slip plane = (111): [112]:x, [111]:y, [110]:z """
xa = []; ya = []; za = [];ty = [];ly=[]
ax = a*math.sqrt(6)/2.
ay = a*math.sqrt(3)
az = a*math.sqrt(2)/2.
bx,by,bz = ax*nx,ay*ny,az*nz
# ux = 0
ux = ax*t
layer=0
for j in range(ny/2):
for i in range(nx):
for k in range(nz):
xa.append((0/6.+i)*ax); ya.append((0/3.+j)*ay); za.append((0/2.+k)*az);ty.append(1);ly.append(layer)
xa.append((3/6.+i)*ax); ya.append((0/3.+j)*ay); za.append((1/2.+k)*az);ty.append(1);ly.append(layer);layer+=1
xa.append((2/6.+i)*ax); ya.append((1/3.+j)*ay); za.append((0/2.+k)*az);ty.append(1);ly.append(layer)
xa.append((5/6.+i)*ax); ya.append((1/3.+j)*ay); za.append((1/2.+k)*az);ty.append(1);ly.append(layer);layer+=1
xa.append((4/6.+i)*ax); ya.append((2/3.+j)*ay); za.append((0/2.+k)*az);ty.append(1);ly.append(layer)
xa.append((1/6.+i)*ax); ya.append((2/3.+j)*ay); za.append((1/2.+k)*az);ty.append(1);ly.append(layer);layer+=1
for j in range(ny/2,ny):
for i in range(nx):
for k in range(nz):
xa.append((0/6.+i)*ax + ux ); ya.append((0/3.+j)*ay); za.append((0/2.+k)*az);ty.append(1);ly.append(layer)
xa.append((3/6.+i)*ax + ux ); ya.append((0/3.+j)*ay); za.append((1/2.+k)*az);ty.append(1);ly.append(layer);layer+=1
xa.append((2/6.+i)*ax + ux ); ya.append((1/3.+j)*ay); za.append((0/2.+k)*az);ty.append(1);ly.append(layer)
xa.append((5/6.+i)*ax + ux ); ya.append((1/3.+j)*ay); za.append((1/2.+k)*az);ty.append(1);ly.append(layer);layer+=1
xa.append((4/6.+i)*ax + ux ); ya.append((2/3.+j)*ay); za.append((0/2.+k)*az);ty.append(1);ly.append(layer)
xa.append((1/6.+i)*ax + ux ); ya.append((2/3.+j)*ay); za.append((1/2.+k)*az);ty.append(1);ly.append(layer);layer+=1
for i in range(len(xa)):
xa[i] = ( (xa[i] + bx)/bx - int((xa[i] + bx)/bx) )* bx
if vaccum == True: ya[i] += 4.0
za[i] = ( (za[i] + bz)/bz - int((za[i] + bz)/bz) )* bz
if vaccum == True: by += 8.0
print layer
return ty,xa,ya,za,bx,by,bz,ly
def shift(a,ly,xa,step=1.0):
ux = a*math.sqrt(6)/2.*step
for i in range(len(ly)):
if ly[i]>ly[len(ly)-1]/2:
xa[i] += ux
return xa
def gen_datafile(ty,xa,ya,za,box):
fout = open("datafile","w")
fout.write("Atomic locations\n\n")
fout.write("%d atoms\n"%len(xa))
fout.write("2 atom types\n")
fout.write(" 0.0 %22.16f xlo xhi\n"%box[0])
fout.write(" 0.0 %22.16f ylo yhi\n"%box[1])
fout.write(" 0.0 %22.16f zlo zhi\n"%box[2])
fout.write("\nAtoms\n\n")
for i in range(len(xa)):
fout.write("%4d %3d %22.16f %22.16f %22.16f\n"%(i+1,ty[i],xa[i],ya[i],za[i]))
fout.close()
return len(xa)
def gen_poscar(xa,ya,za,box):
fout = open("POSCAR","w")
fout.write("Al fcc (111)\n")
fout.write("1.0\n")
fout.write(" %22.16f %22.16f %22.16f\n"%(box[0],0,0))
fout.write(" %22.16f %22.16f %22.16f\n"%(0,box[1],0))
fout.write(" %22.16f %22.16f %22.16f\n"%(0,0,box[2]))
fout.write("%d\n"%len(xa))
fout.write("Selective Dynamics\n")
fout.write("Cart\n")
for i in range(len(xa)):
fout.write("%22.16f %22.16f %22.16f F F T\n"%(xa[i],ya[i],za[i]))
fout.close()
return len(xa)
# Main Program
if len(sys.argv) > 6:
a=float(sys.argv[1])
nx=int(sys.argv[2])
ny=int(sys.argv[3])
nz=int(sys.argv[4])
t=float(sys.argv[5])
option=sys.argv[6]
ty,xa,ya,za,bx,by,bz,ly = get_fcc(a,t,nx,ny,nz)
#xa = shift(a,ly,xa,t)
if option == "d":
gen_datafile(ty,xa,ya,za,[bx,by,bz])
elif option == "p":
gen_poscar(xa,ya,za,[bx,by,bz])
else:
print "Error: wrong number of arguments!!!"
print usage
Running the calculation
Once you generate a POSCAR file from the above script, make an INCAR and a KPOINTS file and copy the POTCAR file relevant to your group to the same directory. After that make follow the directions to submit job to cluster. Example files are given below.
INCAR
LWAVE = .FALSE. LCHARG = .FALSE. LREAL = Auto ISMEAR = -5 ENCUT = 240.3 EDIFF = 1e-6 NSW=100 ISIF=2 IBRION=2
KPOINTS
Automatic mesh 0 ! number of k-points = 0 ->automatic generation scheme Monkhorst-Pack ! Use Monkhorst-Pack scheme 4 1 4 ! only need one k-point in the direction perpendicular to stacking fault 0. 0. 0. ! optional shift of the mesh
GSFE curve generation script
The following script invokes the former fcc (111) surface generation script at different displacement for the top half of the simulation box ranging from 0 to 3.0. Due to the computational cost of this calculation this should be executed using a parallel version of VASP and utilizing a cluster at HPCC using a pbs command script
Save the following script as <name> and give permissions as executable(chmod u+x <name>). It requires two arguments, number of processors and the path of the VASP executable.
#!/bin/sh
rm EvsA
nproc=$1
vasp_exe=$2
totT=0
for a in `seq -w 0.0 0.1 3.0`
do
echo "a= $a"
./<name>.py 4.05 1 4 1 $a p
mpiexec -np $nproc ./job.sh
E=`tail -n1 OSZICAR | awk '{ print $5}'`
echo $a $E >> EvsA
cleanvaspfiles;rm C* W*
done
Submitting job to cluster
Two files are needed to submit to the cluster. One is the pbs command script and the other is a job.sh shell script to invoke ulimit command on all allocated processors. They are presented below.Both files should be in your work directory with the rest of the VASP input files.
Pbs command script
#PBS -N <name of output files> #PBS -l nodes=4:ppn=4 #PBS -l walltime=48:00:00 #PBS -q q64p48h@raptor #PBS -mea #PBS -r n #PBS -V cd $PBS_O_WORKDIR ./<name of GSFE generation script> <no. of processors>
Job.sh script
#!/bin/bash ulimit -s unlimited <Absolute path of VASP exe>Make the file and make it executable by
chmod u+x job.sh
Extracting Data
Energy v. displacement data is written to EvsA file. You could use any plotting program to generate the GSFE plot.
Convergence
- The curve should be converged in terms of k-point grid and number of layers in the direction parallel to the stacking falut . Therefore try different number of layers. Remember that the ratio of k-points should be inversely proportional to the lengths of the lattice vectors (in this case the edges of the simulation box). Since the supercell is elongated only one k-point is needed in the elongated direction.
References
- ↑ C. Brandl, P. M. Derlet, and H. Van Swygenhoven, General-stacking-fault energies in highly strained metallic environments: Ab initio calculations, PHYSICAL REVIEW B 76, 054124 (2007)
