VASP: Computing Bulk Modulus of Au
Modified by Yanming Wang (Sep, 2015)
Input files
Here we give an example of how to use VASP to compute the bulk modulus of LDA-Au. We performed this calculation on MC-CC in serial model in the ~/Codes/VASP/runs/Au/LDA/perfect.21x21x21 directory. This 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
POSCAR for FCC Au (created manually) 4.068 1.0 0 0 0 1.0 0 0 0 1.0 4 selective dynamics direct (real coordinates r) 0 0 0 T T T 0.5000 0.5000 0 T T T 0 0.5000 0.5000 T T T 0.5000 0 0.5000 T T T
POTCAR A POTCAR file that specifies the pseudopotential is required to be in the same directory. In this example, we will present the results from two calculations with US-LDA PP and PAW-LDA PP.
Run VASP
We can write one .pbs script to submit the VASP simulation request to the cluster.
Using MC2 cluster as an example, we provide the following vasp.pbs for the calculation of equilibrium lattice constant, cohesive energy and bulk modulus of Au. The goal of this script is to run vasp in iterations with different lattice constants.
#!/bin/bash #PBS -N vasp.Au #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 4.064 4.065 4.066 4.067 4.068 4.069 4.070 4.071 4.072 do cat > POSCAR << FIN POSCAR for FCC Au (created manually) $a 1.0 0 0 0 1.0 0 0 0 1.0 4 selective dynamics direct (real coordinates r) 0 0 0 T T T 0.5 0.5 0 T T T 0 0.5 0.5 T T T 0.5 0 0.5 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.B.dat p=`grep pressure OUTCAR | cut -b 25-34` echo $a $p >> platt.B.dat done
To submit the job, make sure POTCAR, INCAR and KPOINTS are placed in the same folder of your PBS script, and then type the following command.
qsub vasp.pbs
Analyze data
After the simulation is finished, it will create data files Elatt.B.dat and platt.B.dat.
For simulation with US-LDA PP, the result files should be similar to Elatt.B.dat
4.064 1 -.17546329E+02 -.17546306E+02 -.175463E+02 4.065 1 -.17546473E+02 -.17546449E+02 -.175465E+02 4.066 1 -.17546574E+02 -.17546550E+02 -.175466E+02 4.067 1 -.17546633E+02 -.17546607E+02 -.175466E+02 4.068 1 -.17546649E+02 -.17546622E+02 -.175466E+02 4.069 1 -.17546622E+02 -.17546594E+02 -.175466E+02 4.070 1 -.17546553E+02 -.17546524E+02 -.175466E+02 4.071 1 -.17546441E+02 -.17546412E+02 -.175464E+02 4.072 1 -.17546288E+02 -.17546257E+02 -.175463E+02
platt.B.dat
4.064 5.34 5.34 5.34 4.065 3.95 3.95 3.95 4.066 2.56 2.56 2.56 4.067 1.18 1.18 1.18 4.068 -0.20 -0.20 -0.20 4.069 -1.57 -1.57 -1.57 4.070 -2.93 -2.93 -2.93 4.071 -4.29 -4.29 -4.29 4.072 -5.63 -5.64 -5.64
Launch Matlab and run the following functions fit_a0EB.m and fit_a0B.m,
[a0, Ecoh, B] = fit_a0EB('Elatt.B.dat',4,4);
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. For simulation with US-LDA PP, the result is a0 = 4.068 
, Ecoh = -4.39 eV, B = 186 GPa.