はじめに
今回は、入力パラメータの与え方として簡単なtipsを書きます。あくまで経験則なので、あまり信用しすぎないように。今回はOpenMXのバージョン3.5を使います。おそらく3.6でも共通の話だと思います。タングステンのBccの計算
今回の計算対象はタングステンのBccのバルクとします。 この計算は非直交の辺をもつセルを使えばタングステン原子1つだけの系で計算可能ですが、後々空孔なども計算したいので直交系で行うことにします。すなわち、格子定数aの立方体を1セルとして、1セル内に2つのタングステンを配置したものとします。この1セルだけで構成(1x1x1セル構造)された2原子系がBccを計算するための最小の直交系です。これより大きな系として、x,y,z方向にセルを2列ずつ並べた(2x2x2セル構造)16原子の系や、3列ずつ並べた(3x3x3セル構造)54原子の系があります。 今回は今挙げた三つの系で、構造最適化の計算をしてみます。 まずは、入力ファイルを見てみましょう。例として1x1x1セル構造の入力ファイルを示します。
#
# File Name
#
System.CurrrentDirectory ./ # default=./
System.Name W2
DATA.PATH ../DFT_DATA06 # default=../DFT_DATA
level.of.stdout 1 # default=1 (1-3)
level.of.fileout 1 # default=1 (1-3)
#
# restart using a restart file, *.rst
#
scf.restart off # on|off,default=off
#
# Definition of Atomic Species
#
Species.Number 1
<Definition.of.Atomic.Species
W W8.0-s2p2d1f1 W_PBE
Definition.of.Atomic.Species>
#
# Atoms
#
Atoms.Number 2
Atoms.SpeciesAndCoordinates.Unit Ang # Ang|AU
<Atoms.SpeciesAndCoordinates # Unit=Ang.
1 W 0.0000 0.0000 0.0000 14.0 14.0
2 W 1.590 1.590 1.590 14.0 14.0
Atoms.SpeciesAndCoordinates>
Atoms.UnitVectors.Unit Ang # Ang|AU
<Atoms.UnitVectors # unit=Ang.
3.180 0.000000 0.000000
0.000000 3.180 0.000000
0.000000 0.000000 3.180
Atoms.UnitVectors>
#
# SCF or Electronic System
#
scf.XcType GGA-PBE # LDA|LSDA-CA|LSDA-PW|GGA-PBE
scf.SpinPolarization on # On|Off
scf.ElectronicTemperature 300.0 # default=300 (K)
scf.energycutoff 400.0 # default=150 (Ry)
scf.maxIter 40 # default=40
scf.EigenvalueSolver Band # Recursion|Cluster|Band
scf.lapack.dste dstevx # dstegr|dstedc|dstevx, default=dstegr
scf.Kgrid 6 6 6 # means 4x4x4
scf.Mixing.Type rmm-diisk # Simple|Rmm-Diis|Gr-Pulay
scf.Init.Mixing.Weight 0.01 # default=0.30
scf.Min.Mixing.Weight 0.001 # default=0.001
scf.Max.Mixing.Weight 0.200 # default=0.40
scf.Mixing.History 6 # default=5
scf.Mixing.StartPulay 6 # default=6
scf.criterion 1.0e-6 # default=1.0e-6 (Hartree)
#
# MD or Geometry Optimization
#
MD.Type Opt # Nomd|Opt|DIIS|NVE|NVT_VS|NVT_NH
MD.Opt.DIIS.History 7 # default=7
MD.Opt.StartDIIS 5 # default=5
MD.maxIter 100 # default=1
MD.TimeStep 1.0 # default=0.5 (fs)
MD.Opt.criterion 1.0e-4 # default=1.0e-4 (Hartree/bohr)
|
ここで、格子定数は3.18 angstromとしました。MD.TypeをOptにすることで、全エネルギー計算から得られ得る力の方向に原子を動かすことで、原子位置の最適化(構造最適化)を行います。原子を動かす回数はMD.maxIter で設定した値が上限となりまあすが、MD.Opt.criterionで設定した値よりも原子にかかる力が小さくなったときは、安定構造構造に達したと判断して計算が終了します。 さて、今計算対象としてセットした系はBcc構造のバルクですので、周期境界にする必要があります。scf.EigenvalueSolverをbandに設定することで周期境界条件が適用されます。
実際の計算
さて、上記の入力ファイルの計算条件で実際に計算を行ってみるとどうなるでしょうか。この場合は、Bcc構造が対象ですので、対称性から原子にかかる力が釣り合って常にキャンセルするはずです。格子定数aが最安定長と異なっていても、系のサイズが変わらない周期境界条件下ですので必ず力はキャンセルします。よって、構造安定化の計算を行っても、一度目のイタレーションでMD.Opt.criterion以下の力となり、計算が収束するはずです。この様な予測の元、2原子の1x1x1セル構造、16原子の2x2x2構造、54原子の3x3x3構造の三種の計算を行ったところ、実際に一度目のイタレーションで安定構造と判断され計算が収束します。このとき全ての原子にかかる力が0.000と表示されています。カットオフエネルギー(実空間グリッド)と対称性の破れ
上記では話の都合上、先にうまくいった例を示しました。しかし私が実際の作業で遭遇したのは、Bccの構造安定化計算にも関わらず54原子の3x3x3構造の場合に二度目以降のイタレーションに突入して収束しないという問題でした。2原子の1x1x1セル構造と16原子の2x2x2構造では正しく一度のイタレーションで計算が終了します。ログを見てみると、54原子の3x3x3構造の場合だけ、原子にかかる力がキャンセルアウトせぜに0.000以外の値を取っていました。 正直原因がわからずかなり悩みましたが、何らかの対称性の破れが生じて、力がキャンセルしないのだろうと予想されました。実は先のうまくいった例と、このうまくいかない系で異なるのは、カットオフエネルギーscf.energycutoffが前者では400.0、後者では300.0と与えていたことでした。結論から言うと、このうまくいかなかった場合は、カットオフエネルギーから算出される実空間グリッド数がセル数(今回は3)の倍数になっていないことが原因でした。 どういうことかというと、OpenMXの仕組みと関わってきます。OpenMXではカットオフエネルギーを与えても、その値を基に実空間グリッドを生成して、ポアソン方程式の解の算出など各種の処理を行っているようです。その時に、実空間グリッドのグリッド数と、原子配置させたBccのセル数がずれていると対称性が破れが生じて、数値上の誤差が出てきてしまうようです。 実際の計算ログを見てみましょう。 どの様な実空間グリッドが生成されるかは、計算の最初の方で行われています。 まずは、うまくいった場合(54原子の3x3x3構造カットオフエネルギー400.0)のログ出力の一部です。
Trial cutoff energies (a,b,c) = 7.848 (16), 7.848 (16), 7.848 (16)
Trial cutoff energies (a,b,c) = 31.392 (32), 31.392 (32), 31.392 (32)
Trial cutoff energies (a,b,c) = 70.631 (48), 70.631 (48), 70.631 (48)
Trial cutoff energies (a,b,c) = 125.566 (64), 125.566 (64), 125.566 (64)
Trial cutoff energies (a,b,c) = 196.197 (80), 196.197 (80), 196.197 (80)
Trial cutoff energies (a,b,c) = 282.524 (96), 282.524 (96), 282.524 (96)
Trial cutoff energies (a,b,c) = 384.546 (112), 384.546 (112), 384.546 (112)
Trial cutoff energies (a,b,c) = 502.264 (128), 502.264 (128), 502.264 (128)
UCell_Box: Cutoff=502.264028(128) 502.264028(128) 502.264028(128)
UCell_Box: (tuned) Cutoff=441.442993(120) 441.442993(120) 441.442993(120)
Required cutoff energy (Ryd) for 3D-grids = 400.0000
Used cutoff energy (Ryd) for 3D-grids = 441.4430, 441.4430, 441.4430
Num. of grids of a-, b-, and c-axes = 120, 120, 120
|
上から順番に実空間グリッド数を大きくし(16, 32, 48,・・・)、それに対してTrial cutoff energyも上昇していきます(7.848, 31.392, 70.631,・・・)。Trial cutoff energyが設定値であった400.0を超えたときに、内挿的に実空間グリッドを決定しています。今回は実空間グリッドは120となり、カットオフエネルギーにして441.442993となるものが採用されたようです。最終的な実空間グリッド数は"Num. of grids of a-, b-, and c-axes = 120, 120, 120"と表示されています。 今回採用された実空間グリッドの120という数値は、セルサイズである3で割り切れるので、グリッドとセルの間の対称性は破れていません。 次に、うまくいかなかった場合(54原子の3x3x3構造カットオフエネルギー300.0)のログを見てましょう。
Trial cutoff energies (a,b,c) = 7.848 (16), 7.848 (16), 7.848 (16)
Trial cutoff energies (a,b,c) = 31.392 (32), 31.392 (32), 31.392 (32)
Trial cutoff energies (a,b,c) = 70.631 (48), 70.631 (48), 70.631 (48)
Trial cutoff energies (a,b,c) = 125.566 (64), 125.566 (64), 125.566 (64)
Trial cutoff energies (a,b,c) = 196.197 (80), 196.197 (80), 196.197 (80)
Trial cutoff energies (a,b,c) = 282.524 (96), 282.524 (96), 282.524 (96)
Trial cutoff energies (a,b,c) = 384.546 (112), 384.546 (112), 384.546 (112)
UCell_Box: Cutoff=384.545896(112) 384.545896(112) 384.545896(112)
UCell_Box: (tuned) Cutoff=306.557634(100) 306.557634(100) 306.557634(100)
Required cutoff energy (Ryd) for 3D-grids = 300.0000
Used cutoff energy (Ryd) for 3D-grids = 306.5576, 306.5576, 306.5576
Num. of grids of a-, b-, and c-axes = 100, 100, 100
|
同様のトライアルがなされていますが、採用された実空間グリッド数は112でした。この数値はセルサイズである3の倍数にはなっていない為、グッリドとセルの間の対称性が破れます。これによって、各原子にかかる力に微妙な差が生まれ、構造安定化で収束しないという現象を生み出しているようです。 対称性の問題を回避する為には、実空間グリッドがセルサイズの倍数になるように、注意深くカットオフエネルギーを指定する必要があります。また、論文等に記述する際には、採用されたグリッド数に相当する採用されたカットオフエネルギーもこちらで指定した値とは変わっていることにも注意する必要があります。