# ICME-QM1

Author(s): Laalitha Liyanage

## Contents

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.

Right [1]
Animation of simulation cell for generalized stacking fault curve for Al (111) surface. Black box indicates simulation cell size.

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.

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

1. 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)