晶格常数测试 (Equation of state method)
必要输入文件
run_a0.sh
POTCAR
INCAR 和 KPOINTS可在run_a0.sh中直接设置，也可以单独给出。
EOS.in三维立方晶格脚本示例：
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.5000000000
0.5000000000 0.5000000000 0.0000000000
0.7500000000 0.2500000000 0.7500000000
0.7500000000 0.7500000000 0.2500000000
0.2500000000 0.2500000000 0.2500000000
!
#优化计算
cp INCAR.relax INCAR
echo "a=$i"; time mpirun -n 16 vasp_std
cp CONTCAR POSCAR
rm INCAR
#静态计算
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提交脚本运算后得到ev.out。
利用vaspkit软件进行状态方程拟合。
vaspkit需要准EOS.in文件，文件格式及说明如下：
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 untis (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)参考例子
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
运行 vaspkit -task 205
得到主要的输出文件 PARAM.out 如下：
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通过平衡体积可以得到晶格常数 $a_0=5.4694803016$
建议对得到的平衡晶格常数再优化一次。
如果不想状态方程拟合，在较高精度下直接使用$ISIF = 3$ 的参数直接弛豫优化晶胞。
INCAR示例如下：
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)注：在新版本的vaspkit中，给出了更多状态方程及其整个新的拟合流程，具体可参考程序实例。
二维材料的六方晶格示例：
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
注：对于平面为矩形的二维材料，可以通过修改VASP源代码，直接使用ISIF=3对x、y方向优化。
详情见软件编译手册。
