VASP Computing Elastic Constants of Au

VASP: Computing Elastic Constants of Au

Modified by Yanming Wang (Sep, 2015)

Input files

Here we give an example of how to use VASP to compute the elastic constants C11, C12 and C44 of Au. We performed this calculation on MC2 in PARALLEL mode in ~/Codes/VASP/runs/Au_bulk_C11 and ~/Codes/VASP/runs/Au_bulk_C44 two directories. Each directory contains the following files.

INCAR

PREC = High

ISTART = 0 
ICHARG = 2  
ISMEAR = 1
SIGMA = 0.1

EDIFF  = 1E-09  
NELM = 40      
ENMAX = 500     
ENCUT = 500    

ISIF = 2 
NSW = 100  
IBRION = 2  

KPOINTS

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

POSCAR

Here's an example for Au C11 and C12 calculation, which will be changed later in PBS script every iteration.

POSCAR for FCC Au (created manually)
POSCAR for FCC Au
4.067877         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.0000  0.0000  0.0000 T T T
0.5000  0.5000  0.0000 T T T
0.0000  0.5000  0.5000 T T T
0.5000  0.0000  0.5000 T T T

Here's an example for Si C44 calculation, which will be changed later in PBS script every iteration.

POSCAR for FCC Au
4.067877         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.0000    0.5000    0.0000  T T T
0.2500    0.7500    0.5000  T T T
0.0000    0.0000    0.0000  T T T
0.5000    0.0000    0.0000  T T T
0.2500    0.2500    0.5000  T T T
0.5000    0.5000    0.0000  T T T
0.7500    0.7500    0.5000  T T T
0.7500    0.2500    0.5000  T T T

POTCAR

In this example, we currently choose the US-LDA PP of Au provided by VASP. We will provide results of simulations with Au PAW-LDA PP in the later section.

Run VASP

To compute the elastic constants, we prepared the following PBS scripts to submit the jobs to MC2 cluster. vasp.pbs.C11 is used to to run vasp repeated with different values for the first component of the first Bravais lattice vector, $a. Here we put 4.067877 as the length of the lattice constant in the script, which is obtained from our previous VASP calculation. In generally, you are expected to put a reliable value from your lattice constant calculation with the same PP. The script automatically creates the POSCAR file for every a specified in the script in the range from 0.997 to 1.003. For C44 calculation, the PBS script vasp.pbs.C44 changes the value of $a for the first and second component of the first Bravais lattice vector.

PBS Script for C11

#!/bin/bash
#PBS -N vasp.Au.C11
#PBS -j oe
#PBS -l nodes=1:ppn=8,walltime=48:00:00
#PBS -V

### ---------------------------------------
### BEGINNING OF EXECUTION
### ---------------------------------------

echo The master node of this job is `hostname`
echo The working directory is `echo $PBS_O_WORKDIR`
echo This job runs on the following nodes:
echo `cat $PBS_NODEFILE`

ncpu=`cat $PBS_NODEFILE | wc -w`
echo "Number of processors = $ncpu "

### end of information preamble

cd $PBS_O_WORKDIR

echo $PWD

cd $PBS_O_WORKDIR

module load vasp

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.067877        accurate equilibrium lattice constant
$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.0000  0.0000  0.0000 T T T
0.5000  0.5000  0.0000 T T T
0.0000  0.5000  0.5000 T T T
0.5000  0.0000  0.5000 T T T
FIN

echo "a=$a  ncpu=$ncpu"

cmd="mpiexec -np $ncpu vasp"
$cmd >& vasp.log

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

done

PBS script for C44

#!/bin/bash
#PBS -N vasp.Au.C44
#PBS -j oe
#PBS -l nodes=1:ppn=8,walltime=48:00:00
#PBS -V

### ---------------------------------------
### BEGINNING OF EXECUTION
### ---------------------------------------

echo The master node of this job is `hostname`
echo The working directory is `echo $PBS_O_WORKDIR`
echo This job runs on the following nodes:
echo `cat $PBS_NODEFILE`

ncpu=`cat $PBS_NODEFILE | wc -w`
echo "Number of processors = $ncpu "

### end of information preamble

cd $PBS_O_WORKDIR

echo $PWD

cd $PBS_O_WORKDIR

module load vasp

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.067877        accurate equilibrium lattice constant 
$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.0000    0.5000    0.0000  T T T
0.2500    0.7500    0.5000  T T T
0.0000    0.0000    0.0000  T T T
0.5000    0.0000    0.0000  T T T
0.2500    0.2500    0.5000  T T T
0.5000    0.5000    0.0000  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  ncpu=$ncpu"

cmd="mpiexec -np $ncpu vasp"
$cmd >& vasp.log

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

done

Analyze data

After running the two scripts above for C11 and C44, the data files: Elatt.C11.dat and Elatt.C44.dat will be created.

If you choose your POTCAR for US-LDA PP, the date should be similar to

Elatt.C11.dat

0.997 1  -.17546258E+02  -.17546235E+02 -.175463E+02
0.998 1  -.17546476E+02  -.17546452E+02 -.175465E+02
0.999 1  -.17546606E+02  -.17546580E+02 -.175466E+02
1.000 1  -.17546649E+02  -.17546622E+02 -.175466E+02
1.001 1  -.17546605E+02  -.17546577E+02 -.175466E+02
1.002 1  -.17546476E+02  -.17546446E+02 -.175465E+02
1.003 1  -.17546261E+02  -.17546229E+02 -.175463E+02

Elatt.C44.dat

0.997 1  -.35092701E+02  -.35092739E+02 -.350927E+02
0.998 1  -.35093172E+02  -.35093210E+02 -.350932E+02
0.999 1  -.35093440E+02  -.35093474E+02 -.350934E+02
1.000 1  -.35093511E+02  -.35093539E+02 -.350935E+02
1.001 1  -.35093392E+02  -.35093411E+02 -.350934E+02
1.002 1  -.35093087E+02  -.35093098E+02 -.350931E+02
1.003 1  -.35092601E+02  -.35092607E+02 -.350926E+02

Download the following Matlab functions. fit_a0EB.m, fit_C11EB.m, fit_C44EB.m,

Place the data files of Elatt.C11.dat, Elatt.C44.dat, Elatt.B.dat from Bulk modulus calculation of Au

Launch Matlab, and type the following commands.

[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([a0, Ecoh, B]) respectively. The first 4 specifies the number of atoms in the computation cell and the second 4 is the number of atoms in a conventional unit cell. The second lines computes the C11 and C12 elastic constants and the third line calculates C44. Suggested results are C11 = 206 GPa, C12 = 176 GPa and C44 = 36 GPa for calculations with US-LDA PP.

C44 of 36 GPa is a bit lower than the paper's result of 47 GPa (The experimental value is 45 GPa).