VASP Computing Generalized Stacking Fault Energy of Au

VASP: Genearlized Stacking Fault Energy of Au

Input files

Here are the basic input files required for VASP calculation. Some of the files need to be changed since we need to perform a large number of calculations.

INCAR

System = fcc Au
LWAVE = .FALSE.
ENCUT  =  400
ISMEAR  = 1
SIGMA = 0.1
ISIF = 2

KPOINTS

test convergence
0        0 = automatic generation of k-points
Monkhorst
7 7 11
0 0 0

POSCAR

POSCAR for FCC Au (created manually)
4.0605
 1.22474487139159 0                0
 0                1.73205080756888 0
 0                0                0.70710678118655
6
Cartesian  (real coordinates r)
 0                0                0
 0.40824829046386 0.57735026918963 0
 0.81649658092773 1.15470053837925 0
 0.61237243569579 0                0.35355339059327
 1.02062072615966 0.57735026918963 0.35355339059327
 0.20412414523193 1.15470053837925 0.35355339059327

You also need to put the LDA pseudopotential file as POTCAR in this directory.

Results

Perfect crystal

The following table shows energy convergence with k-points. (Reference value from Au_bulk calculation: E = -4.39 eV/atom.)

KPOINTS	E (eV/atom)	optimal number of CPUs	computational time (second)
3x3x3	-4.4735468	--	--
5x5x5	-4.4121563	--	--
7x7x7	-4.4043722	--	--
7x7x11	-4.3981093	16	341
7x5x13	-4.3974777	16	291

Unrelaxed stacking fault energy

From now on, we fix KPOINTS at 7x5x13, for which the energy of perfect crystal is E0 = -4.3974777 (eV/atom).

To create a stacking fault, we shift the repeat vector b in the x-direction by sqrt(6)/6. The new POSCAR file is (notice the change in 3rd line).

POSCAR for FCC Au (created manually)
4.0605
 1.22474487139159 0                0
 0.40824829046386 1.73205080756888 0
 0                0                0.70710678118655
6
Cartesian  (real coordinates r)
 0                0                0
 0.40824829046386 0.57735026918963 0
 0.81649658092773 1.15470053837925 0
 0.61237243569579 0                0.35355339059327
 1.02062072615966 0.57735026918963 0.35355339059327
 0.20412414523193 1.15470053837925 0.35355339059327

Unrelaxed stacking fault energy
E1 (raw data)	-26.343752 (eV)
E0 (perfect crystal)	-4.3974777 (eV/atom)
Ex (excess) = (E1-E0*6)	0.0411142 (eV)
Area	14.27873 (angstrom^2)
Esf (unrelaxed) = Ex / Area	0.0028794 (eV/angstrom^2) = 46.128 (mJ/m^2)

The result can be compared with previous calculations Esf = 59 (mJ/m^2) (Rosengaard and Skriver, PRB, 47, 12865, 1993) and experimental value of Esf = 42 (mJ/m^2) (Devlin, J. Phys. F: Metal Phys. 4, 1865, 1974).

Relaxed stacking fault energy

We now displace the y-component of the repeat vector b and find the minimum energy. To do so we need to submit the following PBS script.

 relax.isf.pbs

...
ncpu=`cat $PBS_NODEFILE | wc -w`
echo "Number of processors = $ncpu "
...
./auto.relax.isf.par $ncpu
...

which calls the following shell script.

 auto.relax.isf.par

#!/bin/bash

ncpu=$1

rm WAVECAR 
for by in 1.730 1.735 1.740 1.745
do
cat > POSCAR << FIN
POSCAR for FCC Au (created by tcl)
4.0605
 1.22474487139159 0                0
 0.40824829046386 $by              0
 0                0                0.70710678118655
6
Cartesian  (real coordinates r)
 0                0                0
 0.40824829046386 0.57735026918963 0
 0.81649658092773 1.15470053837925 0
 0.61237243569579 0                0.35355339059327
 1.02062072615966 0.57735026918963 0.35355339059327
 0.20412414523193 1.15470053837925 0.35355339059327
FIN

echo "by=$by  ncpu=$ncpu"

mpiexec -np $ncpu ../../../../bin/vasp.mpotts.mva2

E=`tail -1 OSZICAR`
echo $by $E | sed -s 's/F=//; s/E0=//; s/d E =//;' >> E.by.dat

done

Here are the energy values at different $b_{y}$.

Energy minimization
$b_{y}$ (angstrom)	E1 (eV)
1.730	-26.341633
1.735	-26.345601
1.740	-26.345699
1.745	-26.342276

The following Matlab script finds the mininum value of E1 = -26.3461(eV) (at $b_{y}$ = 1.7378 angstrom).

data=[
1.730 1  -.26341633E+02  -.26341112E+02 -.156102E-02
1.735 1  -.26345601E+02  -.26345033E+02 -.170179E-02
1.740 1  -.26345699E+02  -.26345097E+02 -.180565E-02
1.745 1  -.26342276E+02  -.26341660E+02 -.184888E-02
];
x = [min(data(:,1)):0.0001:max(data(:,1))];
[p,s] = polyfit(data(:,1),data(:,3),2);
y = polyval(p,x);
[y_min ind] = min(y); 
x_min = x(ind);
figure(1); plot(data(1:4,1),data(1:4,3),'o', x,y,'-', x_min,y_min,'r*');
y_min, x_min

Relaxed stacking fault energy
E1 (min)	-26.3461 (eV)
E0 (perfect crystal)	-4.3974777 (eV/atom)
Ex (excess) = (E1-E0*6)	0.0387662 (eV)
Area	14.27873 (angstrom^2)
Esf (unrelaxed) = Ex / Area	0.002715 (eV/angstrom^2) = 43.49 (mJ/m^2)

Number of CPUs = 16. Total computation time = 20 minutes.

Unrelaxed generalized stacking fault energy

The unrelaxed generalized stacking fault energy is obtained by displacing the x-component of the repeat vector b (while keeping y-component of b fixed). This is done using the following shell script.

auto.unrelax.gsf.par 

#!/bin/bash

ncpu=$1

rm WAVECAR 
by=1.73205080756888
for bx in 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.408248 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.224745
do
cat > POSCAR << FIN
POSCAR for FCC Au (created by tcl)
4.0605
 1.22474487139159 0                0
 $bx              $by              0
 0                0                0.70710678118655
6
Cartesian  (real coordinates r)
 0                0                0
 0.40824829046386 0.57735026918963 0
 0.81649658092773 1.15470053837925 0
 0.61237243569579 0                0.35355339059327
 1.02062072615966 0.57735026918963 0.35355339059327
 0.20412414523193 1.15470053837925 0.35355339059327
FIN

echo "bx=$bx  by=$by  ncpu=$ncpu"

mpiexec -np $ncpu ../../../../bin/vasp.mpotts.mva2 >& vasp.log

E=`tail -1 OSZICAR`
echo $bx $by $E | sed -s 's/F=//; s/E0=//; s/d E =//;' >> E.bxby.dat

timeused=`grep "Total CPU time used" OUTCAR`
echo $bx $by $timeused  >> time.bxby.dat

done

The results can be plotted using the following Matlab script.

data2=[
0.00 1.73205080756888 1  -.26384866E+02  -.26384852E+02 -.419247E-04
0.05 1.73205080756888 1  -.26366785E+02  -.26366967E+02 0.547795E-03
0.10 1.73205080756888 1  -.26327184E+02  -.26327283E+02 0.294942E-03
0.15 1.73205080756888 1  -.26288320E+02  -.26288001E+02 -.957691E-03
0.20 1.73205080756888 1  -.26269732E+02  -.26269867E+02 0.405146E-03
0.25 1.73205080756888 1  -.26274866E+02  -.26274925E+02 0.176062E-03
0.30 1.73205080756888 1  -.26294139E+02  -.26293966E+02 -.517175E-03
0.35 1.73205080756888 1  -.26323923E+02  -.26323824E+02 -.295196E-03
0.408248 1.73205080756888 1  -.26343752E+02  -.26343211E+02 -.162154E-02
0.45 1.73205080756888 1  -.26331580E+02  -.26331671E+02 0.272932E-03
0.50 1.73205080756888 1  -.26252512E+02  -.26252146E+02 -.109923E-02
];
bx = data2(:,1); E1 = data2(:,4);
Egsf_unrlx = (E1 - (-4.3974777*6)) / 14.27873 * 1.602e-19*1e20 * 1e3;
pp = spline(bx, Egsf_unrlx);
x = [min(bx):0.001:max(bx)]; y = ppval(pp, x);
figure(2); plot(bx,Egsf_unrlx,'o', x,y,'-');
set(gca,'FontSize',16);
xlabel('b_x   (angstrom)'); ylabel('\gamma_{gsf}^{unrlx}   (mJ/m^2)');

Relaxed generalized stacking fault energy

The relaxed generalized stacking fault energy is obtained by displacing the x-component of the repeat vector b, while allowing the y-component of b to relax. This means for every value of $b_{x}$ we will need to a set of calculations with different values of $b_{y}$.

(Insert relaxed GSF curve here)

Relaxed GSF energy
bx (angstrom)	E1 (eV)	E_gsf (eV/angstrom^2)	E_gsf (mJ/m^2)
0.1	??	??	??
0.2	??	??	??
0.3	??	??	??
0.4	??	??	??
0.5	??	??	??
0.6	??	??	??