VASP Computing Bulk Modulus of Si

VASP: Computing Bulk Modulus of Si

Adriano Sanchez，Yanming Wang and Wei Cai

Input files

Here we give an example of how to use VASP to compute the equilibrium lattice constant, cohesive energy and bulk modulus of Si. We performed this calculation on MC2 in parallel mode in the ~/Codes/VASP/runs/Si/LDA/perfect.21x21x21 directory. This directory should contain 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 example for Si

POSCAR for DC Si (created manually)
5.431           universal scaling factor 
0.0    0.5     0.5
0.5    0.0     0.5
0.5    0.5     0.0
2
Direct
-0.125 -0.125 -0.125
 0.125  0.125  0.125

POTCAR

It requires to put one pseudopotential (PP) file POTCAR in this directory. There are different choices of PP and these potential files are provided by VASP. In this example, we choose the ultrasoft-LDA PP. The following is an article explains the differences of different PPs.

https://www.vasp.at/vasp-workshop/slides/pseudoppdatabase.pdf

Run VASP

Now we are ready to run

vasp

To compute the equilibrium lattice constant, cohesive energy and bulk modulus, we use the following script auto.B.serial to run vasp repeated with different lattice constants.

#!/bin/bash
#PBS -N vasp.Si
#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 5.382 5.383 5.384 5.385 5.386 5.387 5.388 5.389 5.390 5.391 5.392 5.393 5.394 5.395 5.396 5.397
do
cat > POSCAR << FIN
cubic diamond
$a             universal scaling factor
0.0    0.5     0.5
 0.5    0.0     0.5
 0.5    0.5     0.0
  2
Direct
 -0.125 -0.125 -0.125
  0.125  0.125  0.125
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.B.dat

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

done

Analyze data

After running it as ./auto.B.serial, it will create data files Elatt.B.dat and platt.B.dat.

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

fit_a0EB('Elatt.B.dat');
fit_a0B ('platt.B.dat');

The first line fits the energy data to a quadratic curve and computes the equilibrium lattice constant, cohesive energy and bulk modulus. The second line fits the pressure data to a linear curve and computes the equilibrium lattice constant and bulk modulus. In this example, the result is a0 = 5.389 angstrom, Ecoh = -5.977 eV, B = 95 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.