晶格常數測試 (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 (讀取現有波函數；如果有的話)
ICHARG = 2 (非自洽：GGA/LDA能帶結構)
LREAL = F (投影算子：自動)
ENCUT = 500 (平面波基組的截斷能量，以eV為單位)
PREC = Accurate (精度級別)
LWAVE = .FALSE. (是否寫入WAVECAR)
LCHARG = .FALSE. (是否寫入CHGCAR)
ADDGRID= .TRUE. (增加網格；有助於GGA收斂)
# LVTOT = .TRUE. (是否將總靜電勢寫入LOCPOT)
# LVHAR = .TRUE. (是否將離子 + 哈特里靜電勢寫入LOCPOT)
# NELECT = (電子數：帶電單元；小心)
# LPLANE = .TRUE. (實空間分佈；超晶胞)
# NPAR = 4 (最大節點數；不要為混合物設置)
Electronic Relaxation
ISMEAR = 0 (高斯展寬；金屬：1)
SIGMA = 0.05 (展寬值，以eV為單位；金屬：0.2)
NELM = 60 (最大電子SCF步驟)
NELMIN = 4 (最小電子SCF步驟)
EDIFF = 1E-08 (SCF能量收斂；以eV為單位)
GGA = PE (PBEsol交換-相關)
Ionic Relaxation
NELMIN = 6 (最小電子SCF步驟)
NSW = 60 (最大電子SCF步驟)
IBRION = 2 (算法：0-MD；1-準牛頓；2-CG)
ISIF = 4 (應力/放鬆：2-離子，3-形狀/離子/V，4-形狀/離子)
EDIFFG = -0.001 (離子收斂；eV/AA)
# ISYM = 2 (對稱性：0=無；2=GGA；3=混合物)
!
cat > INCAR.static <<!
Global Parameters
ISTART = 0 (讀取現有波函數；如果有的話)
ICHARG = 2 (非自洽：GGA/LDA能帶結構)
LREAL = F (投影算子：自動)
ENCUT = 500 (平面波基組的截斷能量，以eV為單位)
PREC = Accurate (精度級別)
LWAVE = .FALSE. (是否寫入WAVECAR)
LCHARG = .FALSE. (是否寫入CHGCAR)
ADDGRID= .TRUE. (增加網格；有助於GGA收斂)
# LVTOT = .TRUE. (是否將總靜電勢寫入LOCPOT)
# LVHAR = .TRUE. (是否將離子 + 哈特里靜電勢寫入LOCPOT)
# NELECT = (電子數：帶電單元；小心)
# LPLANE = .TRUE. (實空間分佈；超晶胞)
# NPAR = 4 (最大節點數；不要為混合物設置)
GGA = PE
Static Calculation
ISMEAR = -5 (四面體方法計算DOS)
#LORBIT = 11 (PAW半徑用於投影DOS)
#NEDOS = 2001 (DOSCAR點)
NELM = 60 (最大電子SCF步驟)
EDIFF = 1E-08 (SCF能量收斂；以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 需要準http://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 (讀取現有波函數；如果有的話)
ICHARG = 2 (非自洽：GGA/LDA能帶結構)
LREAL = F (投影算子：自動)
ENCUT = 500 (平面波基組的截斷能量，以eV為單位)
PREC = Accurate (精度級別)
LWAVE = .FALSE. (是否寫入WAVECAR)
LCHARG = .FALSE. (是否寫入CHGCAR)
ADDGRID= .TRUE. (增加網格；有助於GGA收斂)
# LVTOT = .TRUE. (是否將總靜電勢寫入LOCPOT)
# LVHAR = .TRUE. (是否將離子 + 哈特里靜電勢寫入LOCPOT)
# NELECT = (電子數：帶電單元；小心)
# LPLANE = .TRUE. (實空間分佈；超晶胞)
# NPAR = 4 (最大節點數；不要為混合物設置)
Electronic Relaxation
ISMEAR = 0 (高斯展寬；金屬：1)
SIGMA = 0.05 (展寬值，以eV為單位；金屬：0.2)
NELM = 60 (最大電子SCF步驟)
NELMIN = 4 (最小電子SCF步驟)
EDIFF = 1E-08 (SCF能量收斂；以eV為單位)
GGA = PE (PBEsol交換-相關)
Ionic Relaxation
NELMIN = 6 (最小電子SCF步驟)
NSW = 60 (最大電子SCF步驟)
IBRION = 2 (算法：0-MD；1-準牛頓；2-CG)
ISIF = 3 (應力/放鬆：2-離子，3-形狀/離子/V，4-形狀/離子)
EDIFFG = -0.001 (離子收斂；eV/AA)
# ISYM = 2 (對稱性：0=無；2=GGA；3=混合物)
注：在新版本的 vaspkit 中，給出了更多狀態方程及其整個新的擬合流程，具體可參考程序實例。
二維材料的六方晶格示例：
Tl2O
#!/bin/sh
cat > INCAR.relax <<!
Global Parameters
ISTART = 0 (讀取現有波函數；如果有的話)
ICHARG = 2 (非自洽：GGA/LDA能帶結構)
LREAL = F (投影算子：自動)
ENCUT = 600 (平面波基組的截斷能量，以eV為單位)
PREC = Accurate (精度級別)
LWAVE = .FALSE. (是否寫入WAVECAR)
LCHARG = .FALSE. (是否寫入CHGCAR)
ADDGRID= .TRUE. (增加網格；有助於GGA收斂)
# LVTOT = .TRUE. (是否將總靜電勢寫入LOCPOT)
# LVHAR = .TRUE. (是否將離子 + 哈特里靜電勢寫入LOCPOT)
# NELECT = (電子數：帶電單元；小心)
# LPLANE = .TRUE. (實空間分佈；超晶胞)
# NPAR = 4 (最大節點數；不要為混合物設置)
Electronic Relaxation
ISMEAR = 0 (高斯展寬；金屬：1)
SIGMA = 0.05 (展寬值，以eV為單位；金屬：0.2)
NELM = 60 (最大電子SCF步驟)
NELMIN = 4 (最小電子SCF步驟)
EDIFF = 1E-08 (SCF能量收斂；以eV為單位)
GGA = PE (PBEsol交換-相關)
Ionic Relaxation
NELMIN = 6 (最小電子SCF步驟)
NSW = 100 (最大電子SCF步驟)
IBRION = 2 (算法：0-MD；1-準牛頓；2-CG)
ISIF = 2 (應力/放鬆：2-離子，3-形狀/離子/V，4-形狀/離子)
EDIFFG = -0.001 (離子收斂；eV/AA)
# ISYM = 2 (對稱性：0=無；2=GGA；3=混合物)
!
cat > INCAR.static <<!
Global Parameters
ISTART = 0 (讀取現有波函數；如果有的話)
ICHARG = 2 (非自洽：GGA/LDA能帶結構)
LREAL = F (投影算子：自動)
ENCUT = 600 (平面波基組的截斷能量，以eV為單位)
PREC = Accurate (精度級別)
LWAVE = .FALSE. (是否寫入WAVECAR)
LCHARG = .FALSE. (是否寫入CHGCAR)
ADDGRID= .TRUE. (增加網格；有助於GGA收斂)
# LVTOT = .TRUE. (是否將總靜電勢寫入LOCPOT)
# LVHAR = .TRUE. (是否將離子 + 哈特里靜電勢寫入LOCPOT)
# NELECT = (電子數：帶電單元；小心)
# LPLANE = .TRUE. (實空間分佈；超晶胞)
# NPAR = 4 (最大節點數；不要為混合物設置)
GGA = PE
Static Calculation
ISMEAR = -5 (四面體方法計算DOS)
#SIGMA = 0.05
#LORBIT = 11 (PAW半徑用於投影DOS)
#NEDOS = 2001 (DOSCAR點)
NELM = 60 (最大電子SCF步驟)
EDIFF = 1E-08 (SCF能量收斂；以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 方向優化。