VASP Computing Generalized Stacking Fault Energy of Au
Jump to navigation
Jump to search
VASP: Generalized 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 5 13 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 .
| Energy minimization | |
|---|---|
| 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 = 1.7378).
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
Number of CPUs = 16. Total computation time = 2 hours.
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
0.55 1.73205080756888 1 -.26085728E+02 -.26085280E+02 -.134254E-02
0.60 1.73205080756888 1 -.25829191E+02 -.25828980E+02 -.633646E-03
0.65 1.73205080756888 1 -.25506381E+02 -.25506199E+02 -.546431E-03
0.70 1.73205080756888 1 -.25167205E+02 -.25167286E+02 0.242977E-03
0.75 1.73205080756888 1 -.24887268E+02 -.24887237E+02 -.929289E-04
0.80 1.73205080756888 1 -.24745853E+02 -.24745735E+02 -.354033E-03
0.85 1.73205080756888 1 -.24780442E+02 -.24780426E+02 -.490752E-04
0.90 1.73205080756888 1 -.24974975E+02 -.24974773E+02 -.604431E-03
0.95 1.73205080756888 1 -.25276553E+02 -.25276421E+02 -.395463E-03
1.00 1.73205080756888 1 -.25615702E+02 -.25615797E+02 0.282608E-03
1.05 1.73205080756888 1 -.25927655E+02 -.25927281E+02 -.112246E-02
1.10 1.73205080756888 1 -.26168023E+02 -.26167870E+02 -.459552E-03
1.15 1.73205080756888 1 -.26319686E+02 -.26319608E+02 -.233191E-03
1.20 1.73205080756888 1 -.26387386E+02 -.26387561E+02 0.525098E-03
1.224745 1.73205080756888 1 -.26395070E+02 -.26395345E+02 0.826715E-03
];
bx_unrlx = data2(:,1); E1 = data2(:,4);
Egsf_unrlx = (E1 - (-4.3974777*6)) / 14.27873 * 1.602e-19*1e20 * 1e3;
pp = spline(bx_unrlx, Egsf_unrlx);
x_unrlx = [min(bx_unrlx):0.001:max(bx_unrlx)]; y_unrlx = ppval(pp, x_unrlx);
figure(2); plot(bx_unrlx,max(Egsf_unrlx,0),'.', x_unrlx,y_unrlx,'-');
set(gca,'FontSize',16);
xlabel('b_x / a_0'); ylabel('\gamma_{gsf}^{unrlx} (mJ/m^2)');
xlim([0 max(bx_unrlx)]); ylim([0 2000]);
% print out results
nx = length(bx_unrlx);
for i=1:nx, disp(sprintf('%12.6f %24.12f', bx_unrlx(i), Egsf_unrlx(i))); end
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 we will need to a set of calculations with different values of . This is done using the following shell script. We need to be careful that for every the miniminum energy is sampled within the range of . For example, when , we need to run the script again with = 1.775:0.005:1.815; when , = 1.800:0.005:1.840.
auto.relax.gsf.par
#!/bin/bash
ncpu=$1
rm WAVECAR
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
for by in 1.730 1.735 1.740 1.745 1.750 1.755 1.760 1.765 1.770
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
done
Number of CPUs = 16. Total computation time = 19 hours.
The following Matlab script extracts the minimized energy and plots the relaxed GSF curve.
data3 = load('E.bxby.dat');
N = length(data3(:,1)); ny = 9; nx = floor(N/ny);
bx_rlx = data3(1:ny:nx*ny,1); by = data3(1:ny,2);
E1 = reshape(data3(1:nx*ny,4),ny,nx);
figure(3); plot(E1,'.-');
E1_min = zeros(1,nx);
for i=1:nx,
E1_min(i) = min(spline(by,E1(:,i),[min(by):0.0001:max(by)]));
end
Egsf_rlx = (E1_min - (-4.3974777*6)) / 14.27873 * 1.602e-19*1e20 * 1e3;
pp = spline(bx_rlx, Egsf_rlx);
x_rlx = [min(bx_rlx):0.001:max(bx_rlx)]; y_rlx = ppval(pp, x_rlx);
figure(4); plot(bx_rlx,max(Egsf_rlx,0),'o', x_rlx,y_rlx,'-');
set(gca,'FontSize',16);
xlabel('b_x / a_0'); ylabel('\gamma_{gsf}^{rlx} (mJ/m^2)');
% print out results
for i=1:nx, disp(sprintf('%12.6f %24.12f', bx_rlx(i), Egsf_rlx(i))); end
Plot final results
The following Matlab script summarizes the results and plots the relaxed and unrelaxed GSF curves together.
data_unrlx = [
0.000000 0.000224389703
0.050000 20.286175591244
0.100000 64.716458956782
0.150000 108.319866262614
0.200000 129.174645364115
0.250000 123.414561659193
0.300000 101.791247821056
0.350000 68.375133082561
0.408248 46.128015866956
0.450000 59.784373260085
0.500000 148.494598889395
0.550000 335.617660954438
0.600000 623.438968591740
0.650000 985.615170536872
0.700000 1366.153181970664
0.750000 1680.228084990751
0.800000 1838.888435035886
0.850000 1800.081357655754
0.900000 1581.825346091701
0.950000 1243.470355136624
1.000000 862.962636312892
1.050000 512.967429456257
1.100000 243.286907449047
1.150000 73.128828964477
1.200000 -2.827085882290
1.224745 -11.448138314825
];
data_rlx = [
0.000000 -0.020670060764
0.050000 19.360018794024
0.100000 56.053066178278
0.150000 86.173859662894
0.200000 99.250871762213
0.250000 96.073197595982
0.300000 85.496344853148
0.350000 61.842558692410
0.408248 43.478576333529
0.450000 54.685416520609
0.500000 117.020469910803
0.550000 207.709070871746
0.600000 294.919590647244
0.650000 369.101339702239
0.700000 421.696425371603
0.750000 445.021797830292
0.800000 450.785540359893
0.850000 442.667257101993
0.900000 423.309809152175
0.950000 395.882657649786
% 1.000000 404.582714569149
1.050000 276.895997333094
1.100000 167.991003752519
1.15 62.813035785537
1.20 -2.921550847183
1.224745 -11.448138314825
];
bx_unrlx = data_unrlx(:,1); Egsf_unrlx = data_unrlx(:,2);
% apply offset
Egsf_unrlx = Egsf_unrlx - Egsf_unrlx(1);
Egsf_unrlx = Egsf_unrlx-Egsf_unrlx(end)/bx_unrlx(end)*bx_unrlx;
pp = spline(bx_unrlx, Egsf_unrlx);
x_unrlx = [min(bx_unrlx):0.001:max(bx_unrlx)]; y_unrlx = ppval(pp, x_unrlx);
bx_rlx = data_rlx(:,1); Egsf_rlx = data_rlx(:,2);
% apply offset
Egsf_unrlx = Egsf_unrlx - Egsf_unrlx(1);
Egsf_rlx = Egsf_rlx-Egsf_rlx(end)/bx_rlx(end)*bx_rlx;
pp = spline(bx_rlx, Egsf_rlx);
x_rlx = [min(bx_rlx):0.001:max(bx_rlx)]; y_rlx = ppval(pp, x_rlx);
% partial Burgers vector
bp = sqrt(6)/6;
figure(4);
subplot(1,2,1);
p41 = plot(x_unrlx/bp,y_unrlx,'r-', x_rlx/bp, y_rlx,'b-', ...
bx_unrlx/bp,max(Egsf_unrlx,0),'m.', bx_rlx/bp, max(Egsf_rlx, 0),'k.');
set(gca,'FontSize',12);
xlim([0 max(bx_unrlx/bp)]); ylim([0 2000]);
xlabel(' b_x / b_p ');
ylabel('\gamma_{GSF} (mJ/m^2)');
lg1=legend('unrelaxed', 'relaxed', 2); set(lg1,'FontSize',10);
subplot(1,2,2);
p42 = plot(x_unrlx/bp,y_unrlx,'r-', x_rlx/bp, y_rlx,'b-', ...
bx_unrlx/bp,max(Egsf_unrlx,0),'m.', bx_rlx/bp, max(Egsf_rlx, 0),'k.');
set(gca,'FontSize',12);
xlim([0 1.5]); ylim([0 200]);
xlabel(' b_x / b_p ');
ylabel('\gamma_{GSF} (mJ/m^2)');
lg2=legend('unrelaxed', 'relaxed', 2); set(lg2,'FontSize',10);
