7. 計算精度など補足事項

7.1. 計算精度、収束性

7.1.1. カットオフエネルギーと計算精度

平面波基底を採用している利点の1つとして,カットオフエネルギーは大きくすればするほど必ず全エネルギーは小さくなり,密度汎関数理論の厳密解に近づく,という点が挙げられます。具体例として,面心立方格子のアルミニウム結晶を利用したテスト例を紹介します。 図 7.1 にカットオフエネルギーと全エネルギーの関係を示します。

../_images/suppl_image1.svg

図 7.1 アルミニウム結晶の場合の,カットオフエネルギーと全エネルギーの関係

図から明らかなように,カットオフエネルギーを大きくすると全エネルギーが小さくなり,一定の値に収束しています。 この振る舞いは利用している擬ポテンシャルに依存します。この例では,36 Rydbergほどで原子あたり1 meV程度まで収束しています。どの程度の収束を目指すべきかは対象とする問題によって異なってきますが,通常10meV程度の収束が得られていれば十分であると考えられます。また,全エネルギーを絶対エネルギーで評価するのではなく,(ふたつの構造の全エネルギーの差など)相対エネルギーで評価する場合はより小さなカットオフエネルギーで収束することが期待できます。

7.1.2. k 点サンプリングと計算精度

PHASEは平面波基底を採用しているので,扱える問題は周期系に限られます。したがって,すべての物理量は最終的には第一ブリユアンゾーン内で積分する必要があります。この第一ブリユアンゾーン内の積分の細かさを指定するのが\(k\)点サンプリングです。 \(k\)点サンプリング数と全エネルギーの関係を, 図 7.2 に示します。

../_images/suppl_image2.svg

図 7.2 アルミニウム結晶の場合の,k点数と全エネルギーの関係

\(k\)点数に関しては変分原理が成立するわけではないので,\(k\)点数に応じて全エネルギーが単純減少するわけではない点には 注意が必要です。図 7.2 の例でも,途中全エネルギーが大きくなってから収束へ至っていることが分かります。

なお,カットオフエネルギーの場合と同様ここでみた全エネルギーの絶対エネルギーではなく相対エネルギーの場合はより少ない\(k\)点サンプリング数で収束することが期待できます。

7.1.3. 収束判定と計算精度

SCF計算の収束判定を厳しくすると,原子に働く力をより精度よく計算することが可能となります。通常の構造最適化の場合\(10^{- 8}\) hartree程度の収束判定を採用すれば多くの場合問題なく収束します。他方,分子動力学シミュレーションにおいて保存量を保存させるには,さらに厳しい収束判定を採用する必要があります。

図 7.3 に,SiO2に対して収束判定を変化させながら力の計算を行った結果を例としてしめします。この図から, 力を収束させるためには\(10^{- 10}\) hartree以上の,比較的厳しい収束判定が要求されることが分かります。

../_images/suppl_image3.svg

図 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並列で計算を行いました。

../_images/suppl_image4.svg

図 7.4 波動関数ソルバーによる収束の速さの比較

表 7.1 は波動関数ソルバーによる計算時間の比較です。繰り返し回数は,電荷の混合の回数です。 また,計算時間はOpteron 280 2.4GHz のクラスターマシンで,k 点4 並列で行った場合の結果です。 あくまで参照値とお考えください。

表 7.1 波動関数ソルバーによる計算時間の比較。

手法

繰り返し回数

計算時間

行列対角化法

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法を採用した方が良い場合もありますし,電荷密度の差の混合比を大きくした方がよい場合もあります。 さらに,スピンを考慮した計算の場合,初期スピン分極やスピンを固定するか否かなどによっても収束性は変化します。収束しづらい(しない)場合,上記の設定例などを参考に,最適な電荷密度混合法の設定を行ってください。

表 7.2 Fe (100) 表面の収束に至るまでのSCF 計算回数と到達した全エネルギー

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で比較しました。計算は磁性を考慮して行いました。結果は下図の通りです。

../_images/suppl_image5.svg

図 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法は,すべてのケースで収束し,また平均的に少ない回数で最適化を行うことができました。

表 7.3 構造緩和法の違いによる収束にいたる回数の比較。200 回の更新を経ても力の最大値が10−4 以下とならなかったケースは未収束とした。ケース1 がジクロロシクロヘキサン,ケース2 がルチル型TiO2, ケース3 がSiO2, ケース4 がSi(001) 表面に相当する。

ケース 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}\)