格子定数テスト (状態方程式法)#
必要な入力ファイル
run_a0.sh
POTCAR
INCAR と KPOINTS は run_a0.sh で直接設定できますが、別々に指定することもできます。
EOS.in
三次元立方格子スクリプトの例：
Si
#!/bin/sh
cat > INCAR.relax <<!
グローバルパラメータ
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 (最大ノード数; ハイブリッドには設定しない)
電子緩和
ISMEAR = 0 (ガウススミアリング; 金属:1)
SIGMA = 0.05 (スミアリング値, eV単位; 金属:0.2)
NELM = 60 (最大電子SCFステップ)
NELMIN = 4 (最小電子SCFステップ)
EDIFF = 1E-08 (SCFエネルギー収束; eV単位)
GGA = PE (PBEsol交換相関)
イオン緩和
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 <<!
グローバルパラメータ
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
静的計算
ISMEAR = -5 (DOSのための四面体法)
#LORBIT = 11 (投影DOSのためのPAW半径)
#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 : 結晶の名前 (最大256文字)
natoms : 単位セル内の原子数
etype : 状態方程式のタイプ (以下参照)
vplt1, vplt2, nvplt : エネルギー、圧力などをプロットするための体積間隔
およびプロット内のポイント数
nevpt : 入力されるエネルギー-体積ポイントの数
vpt(i) ept(i) : エネルギー-体積ポイント (VASP単位)
入力単位はVASPのデフォルト単位 (すなわち、A^3およびeV) であることに注意してください。
現在実装されている状態方程式は次の通りです：
1. ユニバーサルEOS (Vinet P et al., J. Phys.: Condens. Matter 1, p1941 (1989))
2. マーナハンEOS (Murnaghan F D, Am. J. Math. 49, p235 (1937))
3. バーチ-マーナハン3次EOS (Birch F, Phys. Rev. 71, p809 (1947))
4. バーチ-マーナハン4次EOS
5. 自然ひずみ3次EOS (Poirier J-P and Tarantola A, Phys. Earth
Planet Int. 109, p1 (1998))
6. 自然ひずみ4次EOS
7. (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
バーチ-マーナハン3次EOS
Birch F, Phys. Rev. 71, p809 (1947)
(デフォルトのVASP単位: eV, アングストロームなど)
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 の例は以下の通り：
グローバルパラメータ
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 (最大ノード数; ハイブリッドには設定しない)
電子緩和
ISMEAR = 0 (ガウススミアリング; 金属:1)
SIGMA = 0.05 (スミアリング値, eV単位; 金属:0.2)
NELM = 60 (最大電子SCFステップ)
NELMIN = 4 (最小電子SCFステップ)
EDIFF = 1E-08 (SCFエネルギー収束; eV単位)
GGA = PE (PBEsol交換相関)
イオン緩和
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 <<!
グローバルパラメータ
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 (最大ノード数; ハイブリッドには設定しない)
電子緩和
ISMEAR = 0 (ガウススミアリング; 金属:1)
SIGMA = 0.05 (スミアリング値, eV単位; 金属:0.2)
NELM = 60 (最大電子SCFステップ)
NELMIN = 4 (最小電子SCFステップ)
EDIFF = 1E-08 (SCFエネルギー収束; eV単位)
GGA = PE (PBEsol交換相関)
イオン緩和
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 <<!
グローバルパラメータ
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
静的計算
ISMEAR = -5 (DOSのための四面体法)
#SIGMA = 0.05
#LORBIT = 11 (投影DOSのためのPAW半径)
#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 方向を最適化できます。