7. 計算精度など補足事項
7.1. 計算精度、収束性
7.1.1. カットオフエネルギーと計算精度
平面波基底を採用している利点の1つとして,カットオフエネルギーは大きくすればするほど必ず全エネルギーは小さくなり,密度汎関数理論の厳密解に近づく,という点が挙げられます。具体例として,面心立方格子のアルミニウム結晶を利用したテスト例を紹介します。 図 7.1 にカットオフエネルギーと全エネルギーの関係を示します。
図 7.1 アルミニウム結晶の場合の,カットオフエネルギーと全エネルギーの関係
図から明らかなように,カットオフエネルギーを大きくすると全エネルギーが小さくなり,一定の値に収束しています。 この振る舞いは利用している擬ポテンシャルに依存します。この例では,36 Rydbergほどで原子あたり1 meV程度まで収束しています。どの程度の収束を目指すべきかは対象とする問題によって異なってきますが,通常10meV程度の収束が得られていれば十分であると考えられます。また,全エネルギーを絶対エネルギーで評価するのではなく,(ふたつの構造の全エネルギーの差など)相対エネルギーで評価する場合はより小さなカットオフエネルギーで収束することが期待できます。
7.1.2. k 点サンプリングと計算精度
PHASEは平面波基底を採用しているので,扱える問題は周期系に限られます。したがって,すべての物理量は最終的には第一ブリユアンゾーン内で積分する必要があります。この第一ブリユアンゾーン内の積分の細かさを指定するのが\(k\)点サンプリングです。 \(k\)点サンプリング数と全エネルギーの関係を, 図 7.2 に示します。
図 7.2 アルミニウム結晶の場合の,k点数と全エネルギーの関係
\(k\)点数に関しては変分原理が成立するわけではないので,\(k\)点数に応じて全エネルギーが単純減少するわけではない点には 注意が必要です。図 7.2 の例でも,途中全エネルギーが大きくなってから収束へ至っていることが分かります。
なお,カットオフエネルギーの場合と同様ここでみた全エネルギーの絶対エネルギーではなく相対エネルギーの場合はより少ない\(k\)点サンプリング数で収束することが期待できます。
7.1.3. 収束判定と計算精度
SCF計算の収束判定を厳しくすると,原子に働く力をより精度よく計算することが可能となります。通常の構造最適化の場合\(10^{- 8}\) hartree程度の収束判定を採用すれば多くの場合問題なく収束します。他方,分子動力学シミュレーションにおいて保存量を保存させるには,さらに厳しい収束判定を採用する必要があります。
図 7.3 に,SiO2に対して収束判定を変化させながら力の計算を行った結果を例としてしめします。この図から, 力を収束させるためには\(10^{- 10}\) hartree以上の,比較的厳しい収束判定が要求されることが分かります。
図 7.3 SiO2 の,収束判定と力の最大値の関係
7.1.4. ベンチマーク計算例(波動関数ソルバーの比較)
7.1.4.1. FCC-Cu
FCC-Cu について,各波動関数ソルバーがどの程度のパフォーマンスを示すか確かめてみましょう。 なお,ここで紹介する計算の入力データは,以下のディレクトリのサブディレクトリにあります。
samples/sol_cmix_test/Cu
入力データ
まず,波動関数ソルバー以外の入力データです。
Control{
condition = initial
cpumax = 1 day
}
accuracy{
cutoff_wf = 25.0 rydberg
cutoff_cd = 225.0 rydberg
num_bands = 10
ksampling{
mesh{
nx = 10
ny = 10
nz = 10
}
}
scf_convergence{
delta_total_energy = 1.e-10 hartree
succession = 3
}
initial_wavefunctions = atomic_orbitals
initial_charge_density = atomic_charge_density
}
structure{
unit_cell_type = primitive
unit_cell{
!#units bohr
a_vector = 0.0000000 3.4704637 3.4704637
b_vector = 3.4704637 0.0000000 3.4704637
c_vector = 3.4704637 3.4704637 0.0000000
}
symmetry{
method = automatic
tspace{
lattice_system = fcc
}
sw_inversion = on
}
atom_list{
atoms{
#tag rx ry rz weight element mobile
0.000 0.000 0.000 1 Cu 0
}
}
element_list{
#tag element atomicnumber
Cu 29
}
}
wavefunction_solver{
次節参照
}
charge_mixing{
mixing_methods{
#tag method rmxs rmxe itr var prec istr nbmix update
broyden2 0.60 0.60 * * on 3 15 RENEW
}
}
printlevel{
base = 1
}
以下に波動関数ソルバーの設定例を示します。
行列対角化
wavefunction_solver{
solvers{
#tag id sol till_n dts dte itr var prec cmix submat
1 matrixdiagon -1 * * * * on 1 off
}
}
lm+msd, 部分空間対角化は波動関数更新後
wavefunction_solver{
solvers{
#tag id sol till_n dts dte itr var prec cmix submat
1 lm+msd 1 * * * * on 1 on
}
submat{
before_renewal=off
}
}
lm+msd, 部分空間対角化は波動関数更新前
wavefunction_solver{
solvers{
#tag id sol till_n dts dte itr var prec cmix submat
1 lm+msd 1 * * * * on 1 on
}
submat{
before_renewal=on
}
}
lm+msd → rmm3, 部分空間対角化は波動関数更新後
wavefunction_solver{
solvers{
#tag id sol till_n dts dte itr var prec cmix submat
1 lm+msd 1 * * * * on 1 on
2 rmm3 -1 * * * * on 1 on
}
rmm{
edelta_change_to_rmm = 5.0e-3
}
submat{ before_renewal=off }
}
lm+msd → rmm3, 部分空間対角化は波動関数更新前
wavefunction_solver{
solvers{
#tag id sol till_n dts dte itr var prec cmix submat
1 lm+msd 1 * * * * on 1 on
2 rmm3 -1 * * * * on 1 on
}
rmm{
edelta_change_to_rmm = 5.0e-3
}
submat{
before_renewal=on
}
}
Davidson → rmm3, 部分空間対角化は波動関数更新後
wavefunction_solver{
solvers{
#tag id sol till_n dts dte itr var prec cmix submat
1 davidson 1 * * * * off 1 off
2 rmm3 -1 * * * * on 1 on
}
rmm{
edelta_change_to_rmm = 5.0e-3
}
submat{
before_renewal=off
}
}
Davidson → rmm3, 部分空間対角化は波動関数更新前
wavefunction_solver{
solvers{
#tag id sol till_n dts dte itr var prec cmix submat
1 davidson 1 * * * * off 1 off
2 rmm3 -1 * * * * on 1 on
}
rmm{
edelta_change_to_rmm = 5.0e-3
}
submat{
before_renewal=on
}
}
Davidson
wavefunction_solver{
solvers{
#tag id sol till_n dts dte itr var prec cmix submat
1 davidson -1 * * * * off 1 off
}
}
結果
ベンチマークテストの計算結果を 図 7.4 に,計算時間を 表 7.1 に示します。 Intel Fortran Compiler 11.1 for Linux でコンパイルして, 2.4GHzのOpteron280 プロセッサを搭載したクラスターマシンにおいてk点4並列で計算を行いました。
図 7.4 波動関数ソルバーによる収束の速さの比較
表 7.1 は波動関数ソルバーによる計算時間の比較です。繰り返し回数は,電荷の混合の回数です。 また,計算時間はOpteron 280 2.4GHz のクラスターマシンで,k 点4 並列で行った場合の結果です。 あくまで参照値とお考えください。
手法 |
繰り返し回数 |
計算時間 |
---|---|---|
行列対角化法 |
13回 |
19.2秒 |
lm+msd, 部分空間対角化は波動関数更新後 |
67回 |
22.2秒 |
lm+msd, 部分空間対角化は波動関数更新前 |
105回 |
32.4秒 |
lm+msd \(\rightarrow\) rmm3, 部分空間対角化は波動関数更新後 |
34回 |
12.4秒 |
lm+msd \(\rightarrow\) rmm3, 部分空間対角化は波動関数更新前 |
16回 |
8.4秒 |
Davidson \(\rightarrow\) rmm3, 部分空間対角化は波動関数更新後 |
23回 |
11.2秒 |
Davidson \(\rightarrow\) rmm3, 部分空間対角化は波動関数更新前 |
15回 |
9.5秒 |
Davidson |
17回 |
11.8秒 |
横軸は繰り返し回数で,縦軸は収束値からの相対的なエネルギーです。 SCF 計算では,変分原理が成り立つので,エネルギーが低いほど,正確な値となります。 行列対角化法は,繰り返し回数で見た収束は速いのですが,1回あたりの計算量は一般に大きく, とくに,系のサイズが大きくなると,この手法は実際上適用できなくなります。 全体的には,途中でrmm3法へ移行するケースが速く収束しています。また,rmm3法は部分空間対角化を波動関数更新前に適用する方がより速く収束しています。
計算する系に依存して,収束の速さ,安定性は変化します。 そこで,その都度,最適な計算手法を選択する事をお勧めいたします。 比較的安定で収束の速い手法として, LM+MSD→RMM3法, Davidson法, Davidson→RMM3法が多くの場合推奨されます。 また,RMM3法を利用する場合部分空間対角化は波動関数更新前に適用した方がよい場合が多いです。 Davidson法を利用する場合,前処理(precon)はoffとしておいた方がよいでしょう。
7.1.4.2. Fe(100) 表面
次に,スピンを考慮した計算の例としてFe(100)表面の計算例を紹介します。 この問題では,波動関数ソルバーは固定とし,電荷密度混合法を変更したテストを行います。 ここで紹介する計算の入力データは,以下のディレクトリーのサブディレクトリーにあります。
samples/sol_cmix_test/Fe100
入力データ
まず,電荷密度混合法以外の入力データです。
control{
condition = initial
max_iteration = 200
}
accuracy{
num_bands = 52
ksampling{
method=monk
mesh{
nx = 6
ny = 6
nz = 1
}
}
cutoff_wf = 30 rydberg
cutoff_cd = 300 rydberg
initial_wavefunctions = atomic_orbitals
initial_charge_density = atomic_charge_density
scf_convergence{
delta_total_energy = 1e-9
succession = 3
}
force_convergence{
max_force = 0.0005 hartree/bohr
}
}
structure{
atom_list{
atoms{
#tag element rx ry rz mobile weight
Fe 0.5 0.5 0 off 1
Fe 0 0 0.0948333333333 off 2
Fe 0 0 0.2845 off 2
Fe 0.5 0.5 0.189666666667 off 2
}
}
ferromagnetic_state{
sw_fix_total_spin=off
total_spin=14
spin_fix_period=5
}
unit_cell{
a_vector = 5.3762704477 0.0 0.0
b_vector = 0.0 5.3762704477 0.0
c_vector = 0.0 0.0 28.3458898822
}
element_list{
#tag element atomicnumber mass zeta deviation
Fe 26 101802.230406 0.375 1.83
}
symmetry{
method = automatic
sw_inversion = on
}
magnetic_state=ferro
}
structure_evolution{
method = gdiis
gdiis{
initial_method = cg
c_forc2gdiis = 0.005 hartree/bohr
}
}
wavefunction_solver{
solvers{
#tag sol till_n prec cmix submat
davidson 1 off 1 off
rmm3 -1 on 1 on
}
rmm{
edelta_change_to_rmm = 5e-3 hartree
}
submat{
before_renewal = on
}
}
charge_mixing{
次節参照
}
printoutlevel{
base = 1
}
以下に電荷密度混合法の設定例を示します。
電荷密度の和と差で異なる混合比を採用,混合アルゴリズムはBroyden2法(case0)
charge_mixing{
spin_density_mixfactor=4
mixing_methods{
#tag no method rmxs rmxe itr var prec istr nbmix update
1 broyden2 0.1 0.1 40 linear on 3 5 renew
}
}
電荷密度の和と差で異なる混合比を採用,混合アルゴリズムはPulay 法(case1)
charge_mixing{
spin_density_mixfactor=4
mixing_methods{
#tag no method rmxs rmxe itr var prec istr nbmix update
1 pulay 0.1 0.1 40 linear on 3 15 renew
}
}
電荷密度の和と差で同じ混合比を採用,混合アルゴリズムはBroyden2 法(case2)
charge_mixing{
spin_density_mixfactor=1
mixing_methods{
#tag no method rmxs rmxe itr var prec istr nbmix update
1 broyden2 0.1 0.1 40 linear on 3 15 renew
}
}
電荷密度の和と差で同じ混合比を採用,混合アルゴリズムはPulay 法(case3)
charge_mixing{
spin_density_mixfactor=1
mixing_methods{
#tag no method rmxs rmxe itr var prec istr nbmix update
1 pulay 0.1 0.1 40 linear on 3 15 renew
}
}
結果
本ベンチマークの結果を,表 7.2 にまとめました。同じ電子状態へ収束していることを 確認するため,得られた全エネルギーもあわせて表示しています。 この例では,「電荷密度の和と差で同じ混合比を採用,混合法はPulay法 」(case3)で最少の回数で収束解を得ることができました。多くの場合case3の 設定によって少ない計算回数で収束解が得られますが,問題によってはBroyden2法を採用した方が良い場合もありますし,電荷密度の差の混合比を大きくした方がよい場合もあります。 さらに,スピンを考慮した計算の場合,初期スピン分極やスピンを固定するか否かなどによっても収束性は変化します。収束しづらい(しない)場合,上記の設定例などを参考に,最適な電荷密度混合法の設定を行ってください。
SCF回数 |
全エネルギー (ha.) |
|
---|---|---|
case0 |
36 |
-153.877775988322 |
case1 |
32 |
-153.877775991437 |
case2 |
34 |
-153.877775825592 |
case3 |
29 |
-153.877775990755 |
7.1.5. 波動関数ソルバーおよび電荷密度混合法の自動設定
PHASEに搭載されている波動関数ソルバーには、MSD法、lm+MSD法、Davidson法、CG法、RMM法、直接対角化法などの基本ソルバーと補助ソルバーとしてのsubspace rotationがあります。さらに、電荷密度ミキサーとして単純混合法、Pulay法、Broydenによる2番目の方法などを搭載しています。これらを、問題に応じて適切に組み合わせることによって高速な収束が期待できます。しかし、このように問題に応じて適切に組み合わせるのは非常に手間がかかる作業です。そこで、PHASEには、適切な波動関数ソルバーや電荷密度混合法をプログラムが自動的に選択する機能が備わっています。ここでは、この機能の利用方法を説明します。この機能によって様々な系を収束させることができていますが、もしなかなか収束させられない場合は、手動で波動関数ソルバーや電荷密度ミキサーの設定を行うようにしてください。
この機能を利用することによって、様々な系で安定に収束解が得られていることを確認しています。例として、鉄(100)表面の収束性を示したものを次の図に示します。
入力パラメーターファイルに波動関数ソルバーおよび電荷密度混合法の指定がない場合の収束性をPHASE v1100とPHASE/0 2014で比較しました。計算は磁性を考慮して行いました。結果は下図の通りです。
図 7.5 波動関数ソルバーおよび電荷密度ミキサーの自動設定
赤線が本機能の結果、緑線がPHASE v1100を利用した結果です。
図から明らかなように、PHASE v1100を利用すると波動関数ソルバー・電荷密度混合法の設定をしない場合は本例題を収束させることができませんが、本機能を利用すると収束させることができています
7.2. 構造最適化
7.2.1. 構造最適化手法
7.2.1.1. 計算例
構造最適化の各アルゴリズムがどのように振る舞うかを調べるため,以下に挙げる系の構造最適化を搭載されているアルゴリズムを利用して実施しました。
ケース1 : cis 型のジクロロシクロヘキサン
ケース2 : ルチル型TiO2
ケース3 : SiO2
ケース4 : Si(001) 表面
これらの入力は,以下のディレクトリー下のサブディレクトリーにあります。
samples/strevl_test
力の収束判定は,すべてのケースで\(10^{- 4}\) hartree/bohrとしました(この条件は,比較的厳しい収束判定条件です)。 また,電子状態計算の収束判定はすべてのケースで\(10^{- 10}\) hartreee 1回としました。 原子配置の更新は, 最大200回行い,それでも収束しないケースは未収束と見なしました。
各々の構造最適化のアルゴリズムは,以下のように設定しました。
quenched MD法
structure_evolution{
method = quench
}
cg 法
structure_evolution{
method = cg
}
gdiis 法
structure_evolution{
method = gdiis
gdiis{
initial_method = cg
c_forc2gdiis = 0.01 hartree/bohr
}
}
このように設定すると,まずはcg法のアルゴリズムに従って最適化が進行します。原子に働く力の最大値がc_forc2gdiisで指定する値よりも小さくなった時点でgdiis法に切り替わります。ただし,最初の3回は原子に働く力の最大値に関わらずにcg法が採用されます。
bfgs 法
structure_evolution{
method = bfgs
gdiis{
initial_method = cg
c_forc2gdiis = 0.01 hartree/bohr
}
}
このように設定すると,まずはcg法のアルゴリズムに従って最適化が進行します。原子に働く力の最大値がc_forc2gdiisで指定する値よりも小さくなった時点でbfgs法に切り替わります。ただし,最初の3回は原子に働く力の最大値に関わらずにcg法が採用されます。bfgs法の詳細設定はgdiisブロックで行い,また各変数の意味もgdiis法の場合と全く同じです。
7.2.1.2. 計算結果
ベンチマークの結果を表 7.3にまとめました。この結果,quenched MD法は収束が遅いことが分かります。 今回のケースでは時間刻みとしてデフォルト値(100 au)を採用しましたが,この値を調整することによって収束が改善される可能性はあります。 GDIIS法はSiO2のケースでは最も速く収束していますが,それ以外のケースではあまり有効に作用していません。 gdiis法の振る舞いは,c_forc2gdiisパラメーターをより小さなものにしたり,電子状態計算の収束判定をより厳しくすることによって改善される可能性はあります。 cg法は比較的安定に最適化ができています。bfgs法は,すべてのケースで収束し,また平均的に少ない回数で最適化を行うことができました。
ケース 1 |
ケース 2 |
ケース 3 |
ケース 4 |
|
---|---|---|---|---|
quenched MD |
未収束 |
115 |
166 |
未収束 |
cg |
195 |
62 |
28 |
124 |
GDIIS |
未収束 |
71 |
13 |
176 |
BFGS |
87 |
38 |
16 |
67 |
7.2.2. 構造更新時の電荷密度、波動関数の予測更新(収束性の向上)
PHASEには、構造最適化や分子動力学シミュレーションを行っている際に、波動関数や電荷密度を原子配置の変化に合わせて“補外”することによって収束性を向上させる機能が備わっています。補外は [Arias94] で紹介されている方法によって行っています。いくつかの系において本機能を適用した結果を示します。いずれのケースにおいても、おおむね高速化されます。
鉄粒界
no prediction |
charge0 |
charge1 |
|
---|---|---|---|
原子配置更新回数 |
11 |
11 |
11 |
SCF回数 |
245 |
185 |
163 |
ジクロロシクロヘキサン分子
no prediction |
charge0 |
charge1 |
|
---|---|---|---|
原子配置更新回数 |
25 |
25 |
25 |
SCF回数 |
481 |
330 |
310 |
SiO2結晶
no prediction |
charge0 |
charge1 |
|
---|---|---|---|
原子配置更新回数 |
16 |
16 |
16 |
SCF回数 |
193 |
144 |
136 |
Si (100) 表面
no prediction |
charge0 |
charge1 |
charge1+wf |
|
---|---|---|---|---|
iteration_ionic |
21 |
21 |
21 |
21 |
SCF回数 |
324 |
200 |
187 |
157 |
上記表中、charge0はsw_charge_predictor = on, charge1はsw_charge_predictor = on とsw_extrapolate_charge = on, charge1+wfはsw_charge_predictor = on , sw_extrapolate_charge = on, sw_wf_predictor = onに相当
Si 64原子300K MD 50回
no prediction |
charge0 |
charge1 |
charge1+wf |
|
---|---|---|---|---|
SCF回数 |
613 |
413 |
387 |
231 |
上記表中、charge0はsw_charge_predictor = on, charge1はsw_charge_predictor = on とsw_extrapolate_charge = on, charge1+wfはsw_charge_predictor = on , sw_extrapolate_charge = on, sw_wf_predictor = onに相当
7.3. PHASE/0 の単位系
PHASE/0において利用される単位は,原則としてハートリー原子単位系です。ここでは,ハートリー原子単位系からそのほかの単位に変換する際の変換係数を記述します。結果の解析の際にご活用ください。
エネルギー |
1 hartree = 2 rydberg = 27.21139615 eV = 4.359745836 \(\times 10^{- 18}\) J |
長さ |
1 bohr = 0.5291772480 Å= 0.5291772480 \(\times 10^{- 10}\)m |
質量 |
1 au mass = 電子の質量 = 9.1094\(\times 10^{- 31}\) kg |
体積 |
1 au volume = 0.1481847426 \(\mathring{\mathrm{A}}^{3}\) = 1.48184726 \(\times 10^{- 29}\) \(m^{3}\) |
速度 |
1 au velocity = 2.187691417 \(\times 10^{- 2}\) Å/s = 2.187691417 \(\times 10^{8}\) m/s |
力 |
1 hartree/bohr = 51.42208259 eV/Å= 8.238725025 \(\times 10^{- 8}\) N |
時間 |
1 au time = 2.418884327 \(\times 10^{- 2}\) fs = 2.418884327 \(\times 10^{- 17}\) s |
ストレス |
1 au stress = 2.903628623 \(\times 10^{9}\) atm = 2.942101703 \(\times 10^{13}\) Pa |
密度 |
1 au density = 1.23013834 \(\times 10^{4}\) amu/\(\mathring{\mathrm{A}}^{3}\) = 9.1093897\(\times 10^{- 4}\) g/\(\text{cm}^{3}\) = 9.1093897 \(\times 10^{- 1}\) kg/\(m^{3}\) |