Lattice Constant Testing (Equation of state method)#
Necessary input files
run_a0.sh
POTCAR
INCAR and KPOINTS can be set directly in run_a0.sh or provided separately.
EOS.in
Three-dimensional cubic lattice script example:
Si
#!/bin/sh
cat > INCAR.relax <<!
Global Parameters
ISTART = 0 (Read existing wavefunction; if there)
ICHARG = 2 (Non-self-consistent: GGA/LDA band structures)
LREAL = F (Projection operators: automatic)
ENCUT = 500 (Cut-off energy for plane wave basis set, in eV)
PREC = Accurate (Precision level)
LWAVE = .FALSE. (Write WAVECAR or not)
LCHARG = .FALSE. (Write CHGCAR or not)
ADDGRID= .TRUE. (Increase grid; helps GGA convergence)
# LVTOT = .TRUE. (Write total electrostatic potential into LOCPOT or not)
# LVHAR = .TRUE. (Write ionic + Hartree electrostatic potential into LOCPOT or not)
# NELECT = (No. of electrons: charged cells; be careful)
# LPLANE = .TRUE. (Real space distribution; supercells)
# NPAR = 4 (Max is no. nodes; don't set for hybrids)
Electronic Relaxation
ISMEAR = 0 (Gaussian smearing; metals:1)
SIGMA = 0.05 (Smearing value in eV; metals:0.2)
NELM = 60 (Max electronic SCF steps)
NELMIN = 4 (Min electronic SCF steps)
EDIFF = 1E-08 (SCF energy convergence; in eV)
GGA = PE (PBEsol exchange-correlation)
Ionic Relaxation
NELMIN = 6 (Min electronic SCF steps)
NSW = 60 (Max electronic SCF steps)
IBRION = 2 (Algorithm: 0-MD; 1-Quasi-New; 2-CG)
ISIF = 4 (Stress/relaxation: 2-Ions, 3-Shape/Ions/V, 4-Shape/Ions)
EDIFFG = -0.001 (Ionic convergence; eV/AA)
# ISYM = 2 (Symmetry: 0=none; 2=GGA; 3=hybrids)
!
cat > INCAR.static <<!
Global Parameters
ISTART = 0 (Read existing wavefunction; if there)
ICHARG = 2 (Non-self-consistent: GGA/LDA band structures)
LREAL = F (Projection operators: automatic)
ENCUT = 500 (Cut-off energy for plane wave basis set, in eV)
PREC = Accurate (Precision level)
LWAVE = .FALSE. (Write WAVECAR or not)
LCHARG = .FALSE. (Write CHGCAR or not)
ADDGRID= .TRUE. (Increase grid; helps GGA convergence)
# LVTOT = .TRUE. (Write total electrostatic potential into LOCPOT or not)
# LVHAR = .TRUE. (Write ionic + Hartree electrostatic potential into LOCPOT or not)
# NELECT = (No. of electrons: charged cells; be careful)
# LPLANE = .TRUE. (Real space distribution; supercells)
# NPAR = 4 (Max is no. nodes; don't set for hybrids)
GGA = PE
Static Calculation
ISMEAR = -5 (tetrahedron method for DOS)
#LORBIT = 11 (PAW radii for projected DOS)
#NEDOS = 2001 (DOSCAR points)
NELM = 60 (Max electronic SCF steps)
EDIFF = 1E-08 (SCF energy convergence; in eV)
EDIFFG = -0.001
!
cat > KPOINTS <<!
A
0
M
9 9 9
0 0 0
!
echo 'a0' 'volume' 'free_energy(eV)' >ev.out
for i in $(seq 5.00 0.05 5.90)
do
cat > POSCAR <<!
Si8
1.0000000000
$i 0.0000000000 0.0000000000
0.0000000000 $i 0.0000000000
0.0000000000 0.0000000000 $i
Si
8
Direct
0.0000000000 0.0000000000 0.0000000000
0.2500000000 0.7500000000 0.7500000000
0.5000000000 0.0000000000 0.5000000000
0.0000000000 0.5000000000 0.5000000
0.5000000000 0.5000000000 0.0000000000
0.7500000000 0.2500000000 0.7500000000
0.7500000000 0.7500000000 0.2500000000
0.2500000000 0.2500000000 0.2500000000
!
#Optimization calculation
cp INCAR.relax INCAR
echo "a=$i"; time mpirun -n 16 vasp_std
cp CONTCAR POSCAR
rm INCAR
#Static calculation
cp INCAR.static INCAR
echo "a=$i"; time mpirun -n 16 vasp_std
V=$(grep "volume" OUTCAR | tail -1 | awk '{printf "%12.9f \n", $5 }')
E=$(grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.9f \n", $5 }')
echo $i $V $E >> ev.out
rm INCAR
done
Submit the script calculation to obtain ev.out.
Use vaspkit software for equation of state fitting.
vaspkit requires the quasi http://EOS.in file, the file format and description are as follows:
cname : name of crystal up to 256 characters
natoms : number of atoms in unit cell
etype : equation of state type (see below)
vplt1, vplt2, nvplt : volume interval over which to plot energy, pressure etc.
as well as the number of points in the plot
nevpt : number of energy-volume points to be inputted
vpt(i) ept(i) : energy-volume points (VASP units)
Note that the input units are VASP default units (i.e., A^3 and eV).
The equations of state currently implemented are:
1. Universal EOS (Vinet P et al., J. Phys.: Condens. Matter 1, p1941 (1989))
2. Murnaghan EOS (Murnaghan F D, Am. J. Math. 49, p235 (1937))
3. Birch-Murnaghan 3rd-order EOS (Birch F, Phys. Rev. 71, p809 (1947))
4. Birch-Murnaghan 4th-order EOS
5. Natural strain 3rd-order EOS (Poirier J-P and Tarantola A, Phys. Earth
Planet Int. 109, p1 (1998))
6. Natural strain 4th-order EOS
7. Cubic polynomial in (V-V0)
Reference example
Si
Si
8
3
124.00 206 500
19
125.000000000 -39.404864520
128.790000000 -40.307996560
132.650000000 -41.074573950
136.590000000 -41.714702140
140.610000000 -42.237893500
144.700000000 -42.653098650
148.880000000 -42.968706720
153.130000000 -43.192650580
157.460000000 -43.332358630
161.880000000 -43.394826990
166.380000000 -43.386635660
170.950000000 -43.313964240
175.620000000 -43.182635890
180.360000000 -42.998104870
185.190000000 -42.765493980
190.110000000 -42.489611560
195.110000000 -42.174957610
200.200000000 -41.825750520
205.380000000 -41.445910010
Run vaspkit -task 205
Obtain the main output file PARAM.out as follows:
Si
Birch-Murnaghan 3rd-order EOS
Birch F, Phys. Rev. 71, p809 (1947)
(Default VASP units: eV, Angstrom etc.)
V0 (A^3) = 163.6206779
E0 (eV) = -43.39626828
B0 = 0.2982734905E-02
B0' = 4.260966253
B0 (GPa) = 87.75507591
Through the equilibrium volume, the lattice constant can be obtained as $a_0=5.4694803016$
If you do not want to fit the equation of state, you can directly use the parameters of ISIF = 3 for cell relaxation optimization at a higher precision.
INCAR example is as follows:
Global Parameters
ISTART = 0 (Read existing wavefunction; if there)
ICHARG = 2 (Non-self-consistent: GGA/LDA band structures)
LREAL = F (Projection operators: automatic)
ENCUT = 500 (Cut-off energy for plane wave basis set, in eV)
PREC = Accurate (Precision level)
LWAVE = .FALSE. (Write WAVECAR or not)
LCHARG = .FALSE. (Write CHGCAR or not)
ADDGRID= .TRUE. (Increase grid; helps GGA convergence)
# LVTOT = .TRUE. (Write total electrostatic potential into LOCPOT or not)
# LVHAR = .TRUE. (Write ionic + Hartree electrostatic potential into LOCPOT or not)
# NELECT = (No. of electrons: charged cells; be careful)
# LPLANE = .TRUE. (Real space distribution; supercells)
# NPAR = 4 (Max is no. nodes; don't set for hybrids)
Electronic Relaxation
ISMEAR = 0 (Gaussian smearing; metals:1)
SIGMA = 0.05 (Smearing value in eV; metals:0.2)
NELM = 60 (Max electronic SCF steps)
NELMIN = 4 (Min electronic SCF steps)
EDIFF = 1E-08 (SCF energy convergence; in eV)
GGA = PE (PBEsol exchange-correlation)
Ionic Relaxation
NELMIN = 6 (Min electronic SCF steps)
NSW = 60 (Max electronic SCF steps)
IBRION = 2 (Algorithm: 0-MD; 1-Quasi-New; 2-CG)
ISIF = 3 (Stress/relaxation: 2-Ions, 3-Shape/Ions/V, 4-Shape/Ions)
EDIFFG = -0.001 (Ionic convergence; eV/AA)
# ISYM = 2 (Symmetry: 0=none; 2=GGA; 3=hybrids)
Note: In the new version of vaspkit, more equations of state and a whole new fitting process are provided, which can be referred to in the program examples.
Example of hexagonal lattice for two-dimensional materials:
Tl2O
#!/bin/sh
cat > INCAR.relax <<!
Global Parameters
ISTART = 0 (Read existing wavefunction; if there)
ICHARG = 2 (Non-self-consistent: GGA/LDA band structures)
LREAL = F (Projection operators: automatic)
ENCUT = 600 (Cut-off energy for plane wave basis set, in eV)
PREC = Accurate (Precision level)
LWAVE = .FALSE. (Write WAVECAR or not)
LCHARG = .FALSE. (Write CHGCAR or not)
ADDGRID= .TRUE. (Increase grid; helps GGA convergence)
# LVTOT = .TRUE. (Write total electrostatic potential into LOCPOT or not)
# LVHAR = .TRUE. (Write ionic + Hartree electrostatic potential into LOCPOT or not)
# NELECT = (No. of electrons: charged cells; be careful)
# LPLANE = .TRUE. (Real space distribution; supercells)
# NPAR = 4 (Max is no. nodes; don't set for hybrids)
Electronic Relaxation
ISMEAR = 0 (Gaussian smearing; metals:1)
SIGMA = 0.05 (Smearing value in eV; metals:0.2)
NELM = 60 (Max electronic SCF steps)
NELMIN = 4 (Min electronic SCF steps)
EDIFF = 1E-08 (SCF energy convergence; in eV)
GGA = PE (PBEsol exchange-correlation)
Ionic Relaxation
NELMIN = 6 (Min electronic SCF steps)
NSW = 100 (Max electronic SCF steps)
IBRION = 2 (Algorithm: 0-MD; 1-Quasi-New; 2-CG)
ISIF = 2 (Stress/relaxation: 2-Ions, 3-Shape/Ions/V, 4-Shape/Ions)
EDIFFG = -0.001 (Ionic convergence; eV/AA)
# ISYM = 2 (Symmetry: 0=none; 2=GGA; 3=hybrids)
!
cat > INCAR.static <<!
Global Parameters
ISTART = 0 (Read existing wavefunction; if there)
ICHARG = 2 (Non-self-consistent: GGA/LDA band structures)
LREAL = F (Projection operators: automatic)
ENCUT = 600 (Cut-off energy for plane wave basis set, in eV)
PREC = Accurate (Precision level)
LWAVE = .FALSE. (Write WAVECAR or not)
LCHARG = .FALSE. (Write CHGCAR or not)
ADDGRID= .TRUE. (Increase grid; helps GGA convergence)
# LVTOT = .TRUE. (Write total electrostatic potential into LOCPOT or not)
# LVHAR = .TRUE. (Write ionic + Hartree electrostatic potential into LOCPOT or not)
# NELECT = (No. of electrons: charged cells; be careful)
# LPLANE = .TRUE. (Real space distribution; supercells)
# NPAR = 4 (Max is no. nodes; don't set for hybrids)
GGA = PE
Static Calculation
ISMEAR = -5 (tetrahedron method for DOS)
#SIGMA = 0.05
#LORBIT = 11 (PAW radii for projected DOS)
#NEDOS = 2001 (DOSCAR points)
NELM = 60 (Max electronic SCF steps)
EDIFF = 1E-08 (SCF energy convergence; in eV)
EDIFFG = -0.001
!
cat > KPOINTS <<!
A
0
G
12 12 1
0 0 0
!
for i in $(seq 3.10 0.05 4.10)
do
j=`echo "$i * 0.5000000000"|bc`
k=`echo "$i * 0.8660254038"|bc`
cat > POSCAR <<!
Tl2O
1.0
$i 0.00000000000 0.00000000000
-$j $k 0.00000000000
0.00000000000 0.00000000000 25.00000000000
Tl O
2 1
Direct
0.3333330099999969 0.6666670059999973 0.5606771909028483
0.6666669459999994 0.3333329709999973 0.4393227710971557
0.0000000000000000 0.0000000000000000 0.5000000000000000
!
cp INCAR.relax INCAR
echo "a=$j"; time mpirun -np 16 vasp_std
cp CONTCAR POSCAR
rm INCAR
cp INCAR.static INCAR
echo "a=$j"; time mpirun -np 16 vasp_std
V=$(grep "volume" OUTCAR | tail -1 | awk '{printf "%12.9f \n", $5 }')
E=$(grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.9f \n", $5 }')
echo $i $V $E >> ev.dat
rm POSCAR
done
Note: For two-dimensional materials with rectangular planes, the VASP source code can be modified to directly use ISIF=3 for optimization in the x and y directions.