軌道の最適化

OpenMXで用いられる基底の動径関数は、軌道最適化法 [30]を用いて変分的に最適化することができます。 メタン分子(入力ファイル「Methane_OO.dat」)を例として、軌道最適化の手順を説明します。 軌道最適化法では、最適化動径関数は、プリミティブ動径関数の線形結合で表されます。線形結合の係数は縮約係数 と呼ばれ、この縮約係数が変分原理に基づき最適化されます。 軌道最適化法でのプリミティブ動径関数と最適化動径関数の数は、次の様にして指定します。

    <Definition.of.Atomic.Species
      H   H5.0-s4>1            H_CA13
      C   C5.0-s4>1p4>1        C_CA13
    Definition.of.Atomic.Species>
水素原子(H)のs軌道に対しては、4つのプリミティブ動径関数の線形結合から1つの最適化動径関数が得られます。 同様に、炭素原子(C)の場合、4つのプリミティブ動径関数の線形結合から1つのs(p)軌道の最適化動径関数が得られます。 さらに、以下のキーワードが軌道最適化法に関連します。
    orbitalOpt.Method          species     # Off|Species|Atoms
    orbitalOpt.Opt.Method        EF        # DIIS|EF
    orbitalOpt.SD.step        0.001        # default=0.001
    orbitalOpt.HistoryPulay     30         # default=15
    orbitalOpt.StartPulay       10         # default=1
    orbitalOpt.scf.maxIter      60         # default=40
    orbitalOpt.Opt.maxIter     140         # default=100
    orbitalOpt.per.MDIter       20         # default=1000000
    orbitalOpt.criterion      1.0e-4       # default=1.0e-4 

    CntOrb.fileout               on        # on|off, default=off
    Num.CntOrb.Atoms             2         # default=1
    <Atoms.Cont.Orbitals
     1
     2
    Atoms.Cont.Orbitals>
入力ファイル「Methane_OO.dat」を用いてOpenMXを通常実行します。
    % ./openmx Methane_OO.dat
  
この計算が正常に終了すると、ファイル「met_oo.out」に軌道最適化の履歴が記録されています。
***********************************************************
***********************************************************
         History of orbital optimization   MD= 1
*********     Gradient Norm ((Hartree/borh)^2)     ********
              Required criterion=  0.000100000000
***********************************************************

   iter=   1  Gradient Norm=  0.057098961101  Uele= -3.217161102876
   iter=   2  Gradient Norm=  0.044668461503  Uele= -3.220120116009
   iter=   3  Gradient Norm=  0.034308306321  Uele= -3.223123238394
   iter=   4  Gradient Norm=  0.025847573248  Uele= -3.226177980300
   iter=   5  Gradient Norm=  0.019106400842  Uele= -3.229294858054
   iter=   6  Gradient Norm=  0.013893824906  Uele= -3.232489198284
   iter=   7  Gradient Norm=  0.010499500005  Uele= -3.235304178159
   iter=   8  Gradient Norm=  0.008362635043  Uele= -3.237652870812
   iter=   9  Gradient Norm=  0.006959703539  Uele= -3.239618540761
   iter=  10  Gradient Norm=  0.005994816379  Uele= -3.241268535418
   iter=  11  Gradient Norm=  0.005298095979  Uele= -3.242657118263
   iter=  12  Gradient Norm=  0.003059655878  Uele= -3.250892948269
   iter=  13  Gradient Norm=  0.001390201488  Uele= -3.255123241210
   iter=  14  Gradient Norm=  0.000780925380  Uele= -3.255179362845
   iter=  15  Gradient Norm=  0.000726631072  Uele= -3.255263012792
   iter=  16  Gradient Norm=  0.000390930576  Uele= -3.250873416989
   iter=  17  Gradient Norm=  0.000280785975  Uele= -3.250333677139
   iter=  18  Gradient Norm=  0.000200668585  Uele= -3.252345643243
   iter=  19  Gradient Norm=  0.000240367596  Uele= -3.254238199726
   iter=  20  Gradient Norm=  0.000081974594  Uele= -3.258146794679
多くの場合、20〜50回の反復計算で収束に達します。 プリミティブ基底関数および最適化基底関数で計算されたメタン分子の全エネルギーの比較結果を以下に示します。
    Primitive basis orbitals
       Utot  =      -7.992569945749 (Hartree) 

    Optimized orbitals by the orbital optimization 
       Utot  =      -8.133746986502 (Hartree)
図 16: プリミティブ基底関数および最適化基底関数で計算された炭素2量体(C$_2$)、メタン分子(CH$_4$)、 ダイヤモンド構造の炭素とケイ素、エタン分子(C$_2$H$_6$)、ヘキサフルオロエタン(C$_2$F$_6$) の全エネルギー。 全エネルギーと全基底数は、C$_2$、ダイヤモンド構造の炭素およびケイ素の場合は1原子についての値、 CH$_4$、C$_2$H$_6$およびC$_2$F$_6$の場合は1分子についての値。
\begin{figure}\begin{center}
\epsfig{file=OrbOpt.eps,width=16cm}
\end{center}
\end{figure}

軌道最適化によって、少ない基底数で高精度な基底関数が得られることが分かります。 図 16 に、プリミティブ基底および最適化基底を使って得られた分子およびバルクの全エネルギーの収束の様子を示します。 ここで扱った全ての系に対して、最適化基底の収束特性が優れていることが分かります。 上記の例でのメタン分子の場合、最適化された基底関数は「C_1.pao」と「H_2.pao」の2つのファイルに出力されます。 これらのファイル「C_1.pao」と「H_2.pao」は擬原子基底関数の入力データとして、OpenMXの計算に対してそのまま使用できます。 最適化された基底関数はファイルに出力されますので、計算しようとする系の基底関数を予め最適化しておくと便利です。 この際、軌道最適化法を適用する系としては、化学的に類似した小さな系を選択することを推奨します。

キーワード「orbitalOpt.Method」には、次の2つのオプションが用意されています。 (1) それぞれの原子上の基底関数が完全に最適化される「atoms」、 (2) それぞれの原子種の基底関数が最適化された「species」。

「入力ファイル」の章でも同様な情報が記載されていますが、ユーザーの利便性のため、 関連するキーワードの詳細を以下に列挙します。
orbitalOpt.scf.maxIter
軌道最適化におけるSCF反復の最大回数を「orbitalOpt.scf.maxIter」キーワードで指定します。

orbitalOpt.Opt.maxIter
軌道最適化の反復の最大回数を「orbitalOpt.Opt.maxIter」キーワードで指定します。軌道最適化の反復は、収束条件が達成しなかった場合でも、同キーワードで設定した回数で終了します。

orbitalOpt.Opt.Method
軌道最適化の収束方法として、2つの手法がサポートされています。 「EF」は固有ベクトル追跡法、「DIIS」は反復部分空間における直接反転法です。それぞれのアルゴリズムは構造最適化のそれと同じです。「orbitalOpt.Opt.Method」キーワードでは「EF」あるいは「DIIS」を指定して下さい。

orbitalOpt.StartPulay
「orbitalOpt.StartPulay」キーワードで指定した最適化ステップから、準ニュートン法である「EF」または「DIIS」法を開始します。

orbitalOpt.HistoryPulay
準ニュートン法である「EF」および「DIIS」法において、次ステップでの縮約係数を推定するために参照する過去の ステップ数を「orbitalOpt.HistoryPulay」キーワードで指定します。

orbitalOpt.SD.step
準ニュートン法である「EF」および「DIIS」法を開始するまでの最適化ステップは最急降下法が適用されます。 最急降下法で使用する前因子は「orbitalOpt.SD.step」キーワードで指定します。 多くのケースにおいて、「orbitalOpt.SD.step」の適切な値は0.001程度となります。

orbitalOpt.criterion
軌道最適化の収束条件((Hartree/bohr)$^2$)を「orbitalOpt.criterion」キーワードで指定します。 「微分のノルム$<$orbitalOpt.criterion」という条件が満たされた時に反復ループが終了します。

CntOrb.fileout
最適化動径関数をファイルに出力したい場合は、「CntOrb.fileout」キーワードを「ON」にする必要があります。

Num.CntOrb.Atoms
最適化動径関数をファイルに出力する際の原子数を「Num.CntOrb.Atoms」キーワードで指定します。

Atoms.Cont.Orbitals
「Atoms.SpeciesAndCoordinates」キーワードの第1列で定義した原子の通し番号を用いて、 最適化基底を出力する原子を「Atoms.Cont.Orbitals」キーワードで次のように指定します。

    <Atoms.Cont.Orbitals
     1
     2
    Atoms.Cont.Orbitals>
記述は「$<$Atoms.Cont.Orbitals」で始め、「Atoms.Cont.Orbitals$>$」で終わります。 行の数は、
「Atoms.Cont.Orbitals」で記述する数字と整合性がなければなりません。 例えば、最適化擬原子軌道が「C_1.pao」と「H_2.pao」などとして保存された場合、その原子種名はキーワード
「Definition.of.Atomic.Species」の 設定における第一列の記号に対応し、記号の後の数字はキーワード「Atoms.SpeciesAndCoordinates」の設定における第一列の数値を意味します。 これらの出力ファイル「C_1.pao」と「H_2.pao」は基底関数の入力データとして使用可能です。



2017-03-07