DFT Assignment
Contents |
Overview
The assignment was broken into five parts, a - e. The base scripts for the KPOINTS, INCAR, POSCAR, and POTCAR files were provided through the ICME website. This group used the PAW-GGA POTCAR file provided. Part a) was to plot the energy volume curves for aluminum in fcc, bcc, and hcp crystal structures and determine the equilibrium lattice constant and the bulk modulus. Part b) was to calculate vacancy and interstitial formation energies. Part c) was to calculate surface formation energy for the (111), (110), and (001) surfaces. Part d) was to calculate the extrinsic and intrinsic stacking fault energies. Part e) was to plot the generalized stacking fault curve.
Part a) Energy Volume Curves
The energy volume curve was plotted for aluminum in the fcc, bcc, and hcp crystal structures. The scripts for each structure are provided below.
Scripts
FCC
rm EvsV EvsA SUMMARY evfit.*
for a in `seq -w 3.75 0.02 4.25`
do
echo "a= $a"
cat >POSCAR <<!
Al_FCC
1
$a 0.00 0.00
0.00 $a 0.00
0.00 0.00 $a
4
Direct
0.00000000 0.00000000 0.0000000
0.50000000 0.00000000 0.5000000
0.50000000 0.50000000 0.0000000
0.00000000 0.50000000 0.5000000
!
mpirun -np 4 /home/ben28/ICME/vasp_5212_talon
V=`grep volume/ion OUTCAR|awk '{print $5}'`
E=`tail -n1 OSZICAR | awk '{ print $5/4}'`
totT=`echo $totT $T|awk '{print $1+$2}'`
echo $a $E >> EvsA
echo $V $E >> EvsV
./cleanvaspfiles
done
cat <<! | evfit
FCC
4
EvsA
evfit.4
!
Emin=`grep Emin evfit.4|awk '{print $2}'`
A0=`grep A0 evfit.4|awk '{print $2}'`
K0=`grep A0 evfit.4|awk '{print $4}'`
tail +5 evfit.4|awk '{print $1,"\t",$2,"\t",$3}' > grace
echo "Data for Al in FCC" >> SUMMARY
echo "==========================================" >> SUMMARY
echo "Equilibrium Energy per Atom = $Emin">>SUMMARY
echo "Equilibrium Lattice constant = $A0">>SUMMARY
echo "Bulk Modulus = $K0">>SUMMARY
echo "Total CPU time (sec) = $totT">>SUMMARY
BCC
rm EvsA EvsV SUMMARY evfit.*
for a in `seq -w 3.0 0.02 3.5`
do
echo "a= $a"
cat >POSCAR <<!
Al_BCC
1
$a 0.00 0.00
0.00 $a 0.00
0.00 0.00 $a
2
Direct
0.00000000 0.00000000 0.0000000
0.50000000 0.50000000 0.5000000
!
mpirun -np 4 /home/ben28/ICME/vasp_5212_talon
V=`grep volume/ion OUTCAR|awk '{print $5}'`
E=`tail -n1 OSZICAR | awk '{ print $5/4}'`
totT=`echo $totT $T|awk '{print $1+$2}'`
echo $a $E >> EvsA
echo $V $E >> EvsV
./cleanvaspfiles
done
cat <<! | evfit
BCC
4
EvsA
evfit.4
!
Emin=`grep Emin evfit.4|awk '{print $2}'`
A0=`grep A0 evfit.4|awk '{print $2}'`
K0=`grep A0 evfit.4|awk '{print $4}'`
tail +5 evfit.4|awk '{print $1,"\t",$2,"\t",$3}' > grace
echo "Data for Al in FCC" >> SUMMARY
echo "==========================================" >> SUMMARY
echo "Equilibrium Energy per Atom = $Emin">>SUMMARY
echo "Equilibrium Lattice constant = $A0">>SUMMARY
echo "Bulk Modulus = $K0">>SUMMARY
echo "Total CPU time (sec) = $totT">>SUMMARY
HCP
rm EvsV EvsA SUMMARY evfit.*
for a in `seq -w 2.75 0.02 3.25`
do
echo "a= $a"
b=`echo $a|awk '{print $1*3**.5}'`
c=`echo $a|awk '{print $1*1.633}'`
cat >POSCAR <<!
Al_HCP
1.0
$a 0.00000000000000 0.00000
0.000000 $b 0.00000
0.000000 0.00000000000000 $c
4
Direct
0.25 0.333333 0.25
0.75 0.833333 0.25
0.25 0.00 0.75
0.75 0.50 0.75
!
mpirun -np 4 /home/ben28/ICME/vasp_5212_talon
V=`grep volume/ion OUTCAR|awk '{print $5}'`
E=`tail -n1 OSZICAR | awk '{ print $5/4}'`
totT=`echo $totT $T|awk '{print $1+$2}'`
echo $a $E >> EvsA
echo $V $E >> EvsV
./cleanvaspfiles
done
cat <<! | evfit
GEN
1
4
EvsA
evfit.4
!
Emin=`grep Emin evfit.4|awk '{print $2}'`
A0=`grep A0 evfit.4|awk '{print $2}'`
K0=`grep A0 evfit.4|awk '{print $4}'`
tail +5 evfit.4|awk '{print $1,"\t",$2,"\t",$3}' > grace
echo "Data for Al in HCP" >> SUMMARY
echo "==========================================" >> SUMMARY
echo "Equilibrium Energy per Atom = $Emin">>SUMMARY
echo "Equilibrium Lattice constant = $A0">>SUMMARY
echo "Bulk Modulus = $K0">>SUMMARY
echo "Total CPU time (sec) = $totT">>SUMMARY
Results
These scripts resulted in the following energy volume curves:
With these results gathered the evfit program extrapolated the minimum energy, equilibrium lattice parameter, and the bulk constant. These results were tabulated in the table below:
In order to obtain the correct bulk modulus a k-points convergence test had to be preformed. The scripts above were placed in a loop to preform the convergence test. This loop ran the code for k-points of 3, 5, 7, 9, 11, 13, 15, and 17. This test resulted in a converged k-points value of 13. The results for the 13 k-points simulation are shown in the table below.
Part b) Vacancy and Interstitial Formation Energies
The methods for calculating the vacancy and interstitial energies are very similar. Both utilize the ground state energy of a 108 atom system. The vacancy formation energy is calculated by removing the center most atom from the system and calculating the new ground state energy. The difference between these two energies is the vacancy formation energy. In order to calculate the interstitial energy an atom must be added to the system. The location of this atom can be determined using one of two methods: octahedral or tetrahedral. The ground state energy of the system with the added atom is then compared to the ground state of the 108 atom system. This difference is the interstitial formation energy. A k-point convergence test was performed to resulting in a converged system at 8 k-points. The interstitial formation simulations were not completed in time.
Vacancy Formation
In order to ensure proper results the system had to be relaxed. This was done using the INCAR file below.
LWAVE = .FALSE. LCHARG = .FALSE. LREAL = Auto -90.000 -89.800 -89.600 -89.400 -89.200 -89.000 -88.800 -88.600 -88.400 0 0.2 0.4 0.6 0.8 1 Energy (eV) Distance (b) Generalized stacking fault curve ENCUT=250 ENCUT=280 ISMEAR = 1 ENCUT = 240.3 EDIFF = 1e-6 NSW=100 ISIF=2 IBRION=2 Once
One of the output files, CONTCAR, was then used as the POSCAR input file. The INCAR file was then adjusted from the relaxation run. The chart below shows the results for the ground state energy of the two systems.
Part c) Surface Formation Energy
When a crystal is split and a new surface is created the surface geometry will relax and possibly reconstruct to allow the surface atoms to reach equilibrium positions. Ideal FCC surfaces such as (111), (110), and (100) for aluminum were analyzed to determine their surface formation energy.
For (111), (110), and (100) surfaces the simulation box must be extended in the z, z, and x directions, respectively, to simulate a vacuum. The first step is to generate the simulation box using the following script found on the ICME course website:
#!/usr/bin/env python
# Purpose: Calculate FCC (111), (110), and (100) surface energies
# Author: Sungho Kim and Laalitha Liyanage
import os
import sys
import math
usage="""
Usage: ./gen_fcc_surface.py a surf vacuum nx ny nz adatom
Mandatory arguments
-------------------
a - equilibrium lattice constant
surf - Type of surface 100, 110 or 111
Optional arguments
-------------------
vacuum - length of vacuum; DEFAULT = 8.0 angstroms
nx,ny,nz - periodicity of supercell; DEFAULT (1,1,1)
adatom - 1/0 (True/False); DEFAULT = 0 (False)
"""
#Default setting
#--------------------------------------------------------------------
nx,ny,nz = 1,1,1
vacuum = 8.0
adatom = 0
#--------------------------Surface (100)-----------------------------
def gen_data_for_fcc(a,nx=2,ny=2,nz=4):
""" Generate datafile of FCC structure with lattice constant a """
xa=[]; ya=[]; za=[]
x0 = 0.0
bx,by,bz = a*nx,a*ny,a*nz+vacuum
x,y,z = bx,by,bz
for i in range(nx):
for j in range(ny):
for k in range(nz):
xa.append(0 + i*a); ya.append( 0 + j*a); za.append( 0 + k*a)
xa.append( 0 + i*a); ya.append(a/2 + j*a); za.append(a/2 + k*a)
xa.append(a/2 + i*a); ya.append( 0 + j*a); za.append(a/2 + k*a)
xa.append(a/2 + i*a); ya.append(a/2 + j*a); za.append( 0 + k*a)
if adatom != 0:
# xa.append(x0); ya.append(x0); za.append(x0+nz*a)
xa.append(bx/2.); ya.append(by/2.); za.append(x0+nz*a)
return xa,ya,za,bx,by,bz
#--------------------------Surface (110)-----------------------------
def gen_data_for_110_fcc(a,nx=4,ny=2,nz=1):
""" 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
z3 = math.sqrt(3)/3. * a
z5 = math.sqrt(3)*2./3. * a
bx,by,bz = ax*nx + vacuum, ay*ny, az*nz
for i in range(nx):
for j in range(ny):
for k in range(nz):
xa.append(x0+i*ax); ya.append(x0+j*ay); za.append(x0+k*az)
xa.append(x2+i*ax); ya.append(y2+j*ay); za.append(x0+k*az)
xa.append(x0+i*ax); ya.append(y3+j*ay); za.append(z3+k*az)
xa.append(x2+i*ax); ya.append(y4+j*ay); za.append(z3+k*az)
xa.append(x0+i*ax); ya.append(y5+j*ay); za.append(z5+k*az)
xa.append(x2+i*ax); ya.append(y6+j*ay); za.append(z5+k*az)
if adatom != 0:
xa.append(x0+nx*ax); ya.append(by/2.); za.append(bz/2.)
return xa,ya,za,bx,by,bz
#--------------------------Surface (111)-----------------------------
def gen_data_for_111_fcc(a,nx=2,ny=2,nz=4):
""" 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
bx,by,bz = ax*nx, ay*ny, az*nz+vacuum
for i in range(nx):
for j in range(ny):
layer = 0
for k in range(nz):
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 adatom != 0:
xa.append(bx/2.); ya.append(by/2.); za.append(x0+nz*az)
return xa,ya,za,bx,by,bz
#----------------------------POSCAR generation------------------------------------------------
def gen_poscar(xa,ya,za,bx,by,bz):
fout = open("POSCAR","w")
fout.write("Fe\n")
fout.write("1.0\n")
fout.write(" %22.16f %22.16f %22.16f\n"%(bx,0,0))
fout.write(" %22.16f %22.16f %22.16f\n"%(0,by,0))
fout.write(" %22.16f %22.16f %22.16f\n"%(0,0,bz))
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\n"%(xa[i],ya[i],za[i]))
# 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) > 2:
a_latt = float(sys.argv[1])
surf = sys.argv[2]
if len(sys.argv) == 4:
vacuum = float(sys.argv[3])
elif len(sys.argv) == 7:
vacuum = float(sys.argv[3])
nx = int(sys.argv[4])
ny = int(sys.argv[5])
nz = int(sys.argv[6])
elif len(sys.argv) == 8:
vacuum = float(sys.argv[3])
nx = int(sys.argv[4])
ny = int(sys.argv[5])
nz = int(sys.argv[6])
adatom = sys.argv[7]
if surf == '100' :
xa,ya,za,bx,by,bz = gen_data_for_fcc(a_latt,nx,ny,nz)
gen_poscar(xa,ya,za,bx,by,bz)
elif surf == '110':
xa,ya,za,bx,by,bz = gen_data_for_110_fcc(a_latt,nx,ny,nz)
gen_poscar(xa,ya,za,bx,by,bz)
elif surf == '111':
xa,ya,za,bx,by,bz = gen_data_for_111_fcc(a_latt,nx,ny,nz)
gen_poscar(xa,ya,za,bx,by,bz)
else:
print "Error: wrong number of arguments!!!"
print usage
This script required several inputs in order to run from the command line: the equilibrium lattice constant, FCC surface, vacuum length, simulation box dimensions, and additional atoms may be optionally specified. The command below creates the POSCAR file needed to run the simulation:
./script 4.0464027 111 10.0 2 2 4
This ensured there was a 10.0 angstrom vacuum length and there were at least 5 layers of atoms in the direction of the vacuum. Other VASP inputs included the POTCAR file, which is provided for aluminum, and the KPOINTS file and the INCAR file, which needed to be adusted. A relaxation run was also required. The relaxation run was done in the same way as the one described in part b. The pre-relaxation and post-relaxation simulation boxes are shown below.
The simulation results were tabulated as seen below.
Part d) Stacking Fault Energies
Overview of Stacking Fault Energy
Stacking faults are formed to connect partial dislocations within a single dislocation. The stacking fault region modifies the properties of the dislocation and is therefore important for the understanding of properties like dislocation mobility. This assignment focused on two types of stacking faults, the intrinsic and the extrinsic stacking faults. The intrinsic stacking fault is defined by stacking sequence, ... ABCACABC ... , where one ( 111) plane is missing in the FCC stacking sequence. This is the defect, that appears between the two partials of a split edge dislocation which bas a (111) plane as its slip plane. The 'A', 'B', and 'C' refer to the different positions of the layers in the FCC structure. The extrinsic stacking fault is defined by the sequence ... ABCBABC ... , where one extra (111) plane bas been inserted in the FCC stacking sequence. These two stacking faults are illustrated by the image below, where the intrinsic stacking fault is on the left and the extrinsic on the right.
Calculations
The procedure to calculate the energies is similar to procedures used previously in this assignment. First the total energy of the un-faulted crystal was calculated. This was then compared to the total energy of the intrinsic or extrinsic fault. The stacking fault energy was defined by the difference between the energy of the faulted and un-faulted systems divided by the area of the surface perpendicular to the fault. This was expressed in the equation below where 'N' is the number of atoms, 'e' is the energy per atom, A is the area described previously, and E is the energy of the system described by the subscript.
The energy-volume curve generated previously gave a value for e = -3.05 ev.
Results
The calculations resulted in the following plots for intrinsic and extrinsic stacking faults. Also provided is the data from the k-points convergence test which was also done for this simulation.
Part e) General Stacking Fault Curve
The assignment here was to generate a generalized stacking fault curve. To simulate the block shearing process, a supercell consisting of 26 atoms was used in these calculations. For the reciprocal space integration, we have used different special k-points in the irreducible Brillouin zone, using the scheme of Monkhorst and Pack. The VASP input decks were provided and a k-points convergence test resulted in stacking fault energy curves seen below. The convergence test showed a converged results for the k-points 13:1:13.
