VASP Computing Elastic Constants of Au

VASP: Computing Elastic Constants of Au

Input files

(Under construction!!!)

Here we give an example of how to use VASP to compute the elastic constants C11 and C44 of LDA-Au. We performed this calculation on MC-CC in PARALLEL model in the ~/Codes/VASP/runs/Au/LDA/perfect.21x21x21 directory. This directory contains the following files.

INCAR

ENCUT  =  400 or 500??
ISIF = 2 ! It relaxes atoms but keep cell volume fixed!
ISMEAR  = 1
SIGMA = 0.1

KPOINTS

21x21x21
0        0 = automatic generation of k-points
Monkhorst
21 21 21
0 0 0

POSCAR for C11

POSCAR for FCC Au (created manually)
 POSCAR for FCC Au
4.062241         from acurate VASP lattice calculation
 1.000  0.0 0.0     first  Bravais lattice vector
 0.0 1.0 0.0     second Bravais lattice vector
 0.0 0.0 1.0     third  Bravais lattice vector
 4               number of atoms per species
selective dynamics
direct           direct or cart (only first letter is significant)
  0.00000  0.00000  0.00000 T T T
  0.50000  0.50000  0.00000 T T T
  0.00000  0.50000  0.50000 T T T
  0.50000  0.00000  0.50000 T T T

To do this calculation, you also need to put the LDA pseudopotential file as POTCAR in this directory.

POSCAR for C44

POSCAR for FCC Au
4.062241         universal scaling factor
  1.000  1.000  0.0     first  lattice vector
 -1.0 1.0 0.0     second lattice vector
  0.0 0.0 1.0     third  lattice vector
 8               number of atoms per species
selective dynamics
direct           direct or cart (only first letter is significant)
         0    0.5000         0  T T T
    0.2500    0.7500    0.5000  T T T
         0         0         0  T T T
    0.5000         0         0  T T T
    0.2500    0.2500    0.5000  T T T
    0.5000    0.5000         0  T T T
    0.7500    0.7500    0.5000  T T T
    0.7500    0.2500    0.5000  T T T

Run VASP

Now we are ready to run vasp

To compute the elastic constants we use the following script auto.B.serial to run vasp repeated with different values for the first C11 component of the first Bravais lattice vector. Make sure you have previously calculated the equilibrium lattice constant for this supercell with 4 atoms, otherwise errors in the calculations could arise.

Script for C11

#!/bin/bash
rm WAVECAR 
for a in 0.997 0.998 0.999 1.000 1.001 1.002 1.003
do
cat > POSCAR << FIN
POSCAR for FCC Au
4.062241         from acurate VASP lattice calculation
$a  0.0 0.0     first  Bravais lattice vector
0.0 1.0 0.0     second Bravais lattice vector
0.0 0.0 1.0     third  Bravais lattice vector
4               number of atoms per species
selective dynamics
direct           direct or cart (only first letter is significant)
0.00000  0.00000  0.00000 T T T
0.50000  0.50000  0.00000 T T T
0.00000  0.50000  0.50000 T T T
0.50000  0.00000  0.50000 T T T
FIN

echo "a=$a"
./vasp

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

p=`grep pressure OUTCAR | cut -b 25-34`
echo $a $p >> platt.C11.dat

done

Script for C44

#!/bin/bash
rm WAVECAR 
for a in 0.997 0.998 0.999 1.000 1.001 1.002 1.003
do
cat > POSCAR << FIN
POSCAR for FCC Au
4.062241         universal scaling factor
$a  $a  0.0     first  lattice vector
-1.0 1.0 0.0     second lattice vector
0.0 0.0 1.0     third  lattice vector
8               number of atoms per species
selective dynamics
direct           direct or cart (only first letter is significant)
0    0.5000         0  T T T
0.2500    0.7500    0.5000  T T T
0         0         0  T T T
0.5000    0         0  T T T
0.2500    0.2500    0.5000  T T T
0.5000    0.5000         0  T T T
0.7500    0.7500    0.5000  T T T
0.7500    0.2500    0.5000  T T T
FIN
echo "a=$a"
./vasp

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

p=`grep pressure OUTCAR | cut -b 25-34`
echo $a $p >> platt.C44.dat

done

Analyze data

After running the two scripts above for C11 and C44 s ./auto.B.serial, it will create data files: Elatt.C11.dat and platt.C11.dat Elatt.C44.dat and platt.C44.dat.

Launch octave and run the following functions fit_a0EB.m and fit_a0B.m, cal_VASP.m

 [a0, Ecoh, B] = fit_a0EB('Elatt.B.dat',4,4);
 [C11,C12] = fit_C11EB('Elatt.C11.dat',a0,B);
 C44 = fit_C44EB('Elatt.C44.dat',a0,C11,C12);

The first line fits the energy data to a quadratic curve and computes the equilibrium lattice constant, cohesive energy and bulk modulus. (Refer to previous calculation?!) The second and third lines compute the C11 and C12 elastic constants and the third line calculates C44. You shuold obtain C11= 213.0 GPa, C12= 179.2 GPa and C44= 40.4 GPa.

Parallel computation

To run vasp in parallel, you need to submit vasp.pbs as

qsub vasp.pbs

You will need the following two files to do this calculation in parallel on SU-AHPCRC: auto.B.par and B.pbs.