ICME-QM2
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 is a tutorial to calculate intrinsic and extrinsic stacking fault energies of fcc (111) surface using VASP .
Overview
Using the following scripts intrinsic stacking fault energy and extrinsic stacking fault energy of fcc (111) surface can be calculated. The figure illustrates the different stacking fault configurations. To calculate the stacking fault energy the following equation can be used.
Where ESFE is the stacking fault energy, Etot is the energy of the system with the stacking fault, N is the number of atoms, 
is energy per atom in bulk structure, and A is the area of the surface perpendicular to the stacking fault.
_Stacking_Faults.jpg.html)
Crystal structure generation
Save the following script as<name>.pyand give permissions as executable
chmod u+x <name>.py
#!/usr/local/bin/python
usage = """
Usage: gen_datafile_fcc.py a nx ny type p/d
a - Lattice parameter
nx,ny - periodicity in x and y
type - Stacking fault type I(intrinsic)/E(extrinsic)
p/d - output file format type p(POSCAR)/d (LAMMPS)
"""
import os
import re
import sys
import math
# Default setting
a=4.05
def gen_data_for_stacking_fault(a,nx=1,ny=1,type="I"):
""" Generate datafile of FCC surface: 110:x, 112:y, 111:z """
xa=[]; ya=[]; za=[]
ax = a*math.sqrt(2)/2
ay = a*math.sqrt(6)/2
az = a*math.sqrt(3)
x0 = 0.0
x2 = math.sqrt(2)/4 * a
y2 = math.sqrt(6)/4 * a
y3 = math.sqrt(6)/6 * a
y4 = math.sqrt(6)*5/12 * a
y5 = math.sqrt(6)*2/6 * a
y6 = math.sqrt(6)/12 * a
if type == "I":
nlayer = 11
nfaults = 1
elif type == "E":
nlayer = 13
nfaults = 1
bx,by,bz = ax*nx, ay*ny, az/3.0*nlayer
for i in range(nx):
for j in range(ny):
layer = 0
for k in range(2):
xa.append(x0+i*ax); ya.append(x0+j*ay); za.append(layer/3.0*az)
xa.append(x2+i*ax); ya.append(y2+j*ay); za.append(layer/3.0*az); layer += 1
xa.append(x0+i*ax); ya.append(y3+j*ay); za.append(layer/3.0*az)
xa.append(x2+i*ax); ya.append(y4+j*ay); za.append(layer/3.0*az); layer += 1
xa.append(x0+i*ax); ya.append(y5+j*ay); za.append(layer/3.0*az)
xa.append(x2+i*ax); ya.append(y6+j*ay); za.append(layer/3.0*az); layer += 1
if type == "I":
xa.append(x0+i*ax); ya.append(y3+j*ay); za.append(layer/3.0*az)
xa.append(x2+i*ax); ya.append(y4+j*ay); za.append(layer/3.0*az); layer += 1
xa.append(x0+i*ax); ya.append(y5+j*ay); za.append(layer/3.0*az)
xa.append(x2+i*ax); ya.append(y6+j*ay); za.append(layer/3.0*az); layer += 1
xa.append(x0+i*ax); ya.append(x0+j*ay); za.append(layer/3.0*az)
xa.append(x2+i*ax); ya.append(y2+j*ay); za.append(layer/3.0*az); layer += 1
xa.append(x0+i*ax); ya.append(y3+j*ay); za.append(layer/3.0*az)
xa.append(x2+i*ax); ya.append(y4+j*ay); za.append(layer/3.0*az); layer += 1
xa.append(x0+i*ax); ya.append(y5+j*ay); za.append(layer/3.0*az)
xa.append(x2+i*ax); ya.append(y6+j*ay); za.append(layer/3.0*az); layer += 1
elif type == "E":
xa.append(x0+i*ax); ya.append(y3+j*ay); za.append(layer/3.0*az)
xa.append(x2+i*ax); ya.append(y4+j*ay); za.append(layer/3.0*az); layer += 1
for k in range(2):
xa.append(x0+i*ax); ya.append(x0+j*ay); za.append(layer/3.0*az)
xa.append(x2+i*ax); ya.append(y2+j*ay); za.append(layer/3.0*az); layer += 1
xa.append(x0+i*ax); ya.append(y3+j*ay); za.append(layer/3.0*az)
xa.append(x2+i*ax); ya.append(y4+j*ay); za.append(layer/3.0*az); layer += 1
xa.append(x0+i*ax); ya.append(y5+j*ay); za.append(layer/3.0*az)
xa.append(x2+i*ax); ya.append(y6+j*ay); za.append(layer/3.0*az); layer += 1
return xa,ya,za,bx,by,bz,nfaults
def gen_datafile(xa,ya,za,box):
fout = open("datafile","w")
fout.write("FCC 111 surface in Z\n\n")
fout.write("%d atoms\n"%len(xa))
fout.write("1 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 1 %3d %22.16f %22.16f %22.16f\n"%(i+1,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) > 5:
a=float(sys.argv[1])
nx=int(sys.argv[2])
ny=int(sys.argv[3])
type=sys.argv[4]
option=sys.argv[5]
xa,ya,za,bx,by,bz,nfaults = gen_data_for_stacking_fault(a,nx,ny,type)
if option == "d":
gen_datafile(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 a pbs command script and submit the calculation to a cluster (Raptor). 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
Convergence of k-point grid should be tested. A single k-point should be sufficient in the direction perpendicular to the stacking fault.
Automatic mesh 0 ! number of k-points = 0 ->automatic generation scheme Monkhorst-Pack ! Use Monkhorst-Pack scheme 4 4 1 ! only need one k-point in the direction perpendicular to the surface 0. 0. 0. ! optional shift of the mesh
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 mpiexec -np 16 ./job.sh
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
Once the calculation is complete you will see some output files. CONTCAR gives the final geometry of the system. You could visualize CONTCAR using OVITO to make sure the structure and surface are stable. OSZICAR gives data of each electronic and ionic step. The final energy of the system (Etot) is in the OSZICAR file as well. To get the final energy execute the following command
tail -n1 OSZICAR|awk '{print $5}'
Convergence
- The curve should be converged in terms of k-point grid. Therefore try different k-point grids. 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
- ↑ M. Muzyk, Z. Pakiela, K.J. Kurzydlowski, Ab initio calculations of the generalized stacking fault energy in aluminium alloys, Scripta Materialia, Volume 64, Issue 9, May
