HOME > ナノ材料科学研究グループ > RMM-DIIS

>>	グループTOP
>>	Staff
>>	研究テーマ
>>	研究成果
>>	代表連絡先
>>	HOMEへ戻る

RMM-DIIS

１　はじめに

現在の第一原理分子動力学法^１）は、原子の動きの度毎に密度汎関数理論^２）に基づいて平面波基底の第一原理擬ポテンシャル法^３）で電子状態を求め、全エネルギーや原子に働く力を計算するものである。第一原理擬ポテンシャル法の電子状態計算では、周期境界条件のセルをバンド計算における単位胞として用い（スーパーセル法）、入力の電子密度と出力の電子密度が自己無撞着に一致するように巨大な固有値方程式（Kohn-Sham方程式）を解く。Car-Parrinello^４）に始まる高速化アルゴリズムは、この自己無撞着な固有状態計算を大きく高速化するものであり、Car-Parrinello法以降、最急降下法や共役勾配法など、類似手法が各種提案されている。
これらのアルゴリズムでは、適当なインプットの波動関数のセットから始めて、ハミルトニアンＨ×波動関数φの演算を繰り返しながら波動関数を固有状態の方向に少しずつ変化させ、最終的に基底状態の電子状態（自己無撞着な固有関数）のセットに収束させる。これらのアルゴリズムは、（１）波動関数の更新とハミルトニアンのポテンシャル（charge）の更新とを同時に行う手法（Car-Parrinello法、Teterらの共役勾配法^５））と、（２）ハミルトニアンの対角化とポテンシャル（charge）の更新とを分けて行い、前者に新しいアルゴリズムを用いる手法（Bylanderらの共役勾配法^６）、Block-Davidson法^７）、RMM-DIIS法^８，９）など）とがある。最近は、金属も半導体も同じ立場で扱える手法として、（２）を用いることが多い。
我々のグループでは（２）のうちで、金属系の大きなセルを扱うには、Bylanderらの共役勾配法の方がBlock-Davidson法よりも効率的であることを示した^１０）。しかし、RMM-DIIS法（以下、簡単のためDIIS法）についてはあまり検討されていない。一方、今後の計算機技術はメモリー分散の並列機が主流と考えられ、そうした場合の並列化に対応できるアルゴリズムが有利である。従来の共役勾配法等では、「ｋ点並列」は容易だが、各波動関数を別々のCPUで扱う「バンド間並列」については、波動関数の互いの直交化のためのCPU間の波動関数データの転送のため、効率が落ちることが知られている。特に大規模系では、k=0（Γ点）のみを用いることも多く、そうした場合には「バンド間並列」が不可避となる。DIIS法は、共役勾配法等に比べ、原理的に状態間の直交化処理の回数が少なくてすむため、「バンド間並列」が効率的に行える可能性がある。従って、今後の大規模系のための効率的な第一原理計算手法として、DIIS法の詳細な検討が必要である。

２　RMM-DIIS法の概要

RMM-DIIS法（residual minimization / direct inversion in the iterative subspace）は、Pulay^８）により、固有値問題のみならず、charge-mixing法や格子緩和法^１１）など、非線形最適化問題の一般的解法として
提案され、Pulay法とも呼ばれる。平面波基底擬ポテンシャル法の固有状態計算では、Wood-Zunger^１２）やHutterら^１３）による検討を経て、Kresse-Furthmuller^９）が詳細に検討している。今回はKresse-Furthmullerのやり方に沿って検討する。
上述のように?の方法論では、ポテンシャル（charge）を固定したハミルトニアンの固有状態計算をDIIS法で行い、その出力によるchargeの更新はKerker法^１４）など適当な方法を用い、この両者を交互に繰り返す。DIIS法による固有状態計算は、各k点の各バンドの状態|φ_ｍ>（ｍがk点とバンドを指定）毎に以下のように行う。
まず、現在の波動関数|φ_ｍ^０>に対して、残差ベクトル|Ｒ_ｍ^０>は、

|Ｒ_ｍ^０>＝（Ｈ－ε_ｍ^０）|φ_ｍ^０>          (1)
ε_ｍ^０＝<φ_ｍ^０|Ｈ|φ_ｍ^０>          (2)

で与えられる。残差ベクトルがゼロになれば、|φ_ｍ^０>はＨの固有状態に収束していることになる。
次にpreconditioningとして、Ｋ|Ｒ_ｍ^０>を計算する（Ｋは対角行列、後述）。ここまでの過程は、共役勾配法や最急降下法と基本的に同様であり、|Ｒ_ｍ^０>にマイナスをつければ、共役勾配法の勾配ベクトルになる。
新しい波動関数を次のように組み立てる（trial stepと呼ぶ）。

|φ_ｍ^１>＝|φ_ｍ^０>＋λＫ|Ｒ_ｍ^０>          (3)

λは期待値<φ_ｍ^１|Ｈ|φ_ｍ^１>／<φ_ｍ^１|φ_ｍ^１>が最小になるように決める（詳細は後述）。この過程は、波動関数の配位空間のエネルギー曲面の勾配方向に変化させることであり、最急降下法と同様である。共役勾配法では、共役方向ベクトルを組み立ててその方向に変化させる。
この|φ_ｍ^１>に対し(1)式のように残差ベクトル|Ｒ_ｍ^１>を計算する。次に、|φ_ｍ^０>、|φ_ｍ^１>と|Ｒ_ｍ^０>、|Ｒ_ｍ^１>に対して２×２のDIIS過程を行う。すなわち、|φ_ｍ^０>と
|φ_ｍ^１>の線形結合で新しい波動関数を次のように作る。

|φ_ｍ^Ｍ>＝Σα_ｉ|φ_ｍ^ｉ>　(i=0～M, M=1)          (4)

係数｛α_ｉ｝は、|φ_ｍ^Ｍ>の規格化条件のもとで線形性を想定した以下の残差

|Ｒ_ｍ^Ｍ>＝Σα_ｉ|Ｒ_ｍ^ｉ>　(i=0～M, M=1)          (5)

のノルムの二乗<Ｒ_ｍ^Ｍ|Ｒ_ｍ^Ｍ>を最小化する条件で決める。ラグランジュの未定係数法から、{α_ｉ}は、

[Ｒ]_ｉｊ{α_ｊ}＝ε[Ｓ]_ｉｊ{α_ｊ}          (6)
[Ｒ]_ｉｊ＝<Ｒ_ｍ^ｉ|Ｒ_ｍ^ｊ>, [Ｓ]_ｉｊ＝<φ_ｍ^ｉ|φ_ｍ^ｊ>          (7)

の固有値方程式の最低解になる。これを解けば新しい波動関数|φ_ｍ^Ｍ>が残差を最小化する方向で決まる。
さらに、その残差ベクトルを求め、trial step

|φ_ｍ^Ｍ＋１>＝|φ_ｍ^Ｍ>＋λＫ|Ｒ_ｍ^Ｍ>          (8)

から、新たな波動関数|φ_ｍ^Ｍ＋１>が生まれ、残差ベクトル|Ｒ_ｍ^Ｍ＋１>が計算される。今度は、|φ_ｍ^Ｍ＋１>=|φ_ｍ^２>を基底に加え、|φ_ｍ^０>、|φ_ｍ^１>、|φ_ｍ^２>について３×３のDIIS過程(4)～(7)を行うことになる。

このようにバンド毎にtrial stepとDIIS過程を繰り返しながら、波動関数の更新を行う。DIIS過程の次元は次第に増えるので、適当な回数で打ち切る。全体の状態を一巡した後、chargeの更新を行い、ハミルトニアン内のポテンシャルを更新して、また初めからDIIS法を行う。
重要なことは、DIIS法では、k点毎の各状態間での直交化処理を原理的には行わなくてもよいことである（実際には必要だが）。最急降下法や共役勾配法では、エネルギー期待値を最低化するように、波動関数の配位空間のエネルギー曲面の勾配に従って波動関数を変化させる。その際、各状態間の直交化の拘束条件を厳密に課していないと、全ての状態が一番下の同じレベルに収束してしまう。一方、DIIS法では、各状態が各々ある程度固有状態（残差をゼロにする状態）に近づいていれば、過去の波動関数の線形結合で各々近傍の残差ゼロの状態に近づくので、直交化の拘束条件は直接には不要であり、固有状態に収束すれば自然に直交化することになる。

３　具体的なアルゴリズム

３．１　DIISループの詳細
今回テストした具体的な計算過程を図１に示す。インプットのchargeからのポテンシャルを固定したハミルトニアンの対角化計算を、各k点の各状態毎に、一番内側のループでDIIS法で行う。ループの詳細を説明する。

１巡目は、インプットの|φ_ｍ^０>に対し、残差|Ｒ_ｍ^０>を組み立て、trial step|φ_ｍ^１>＝|φ_ｍ^０>＋λＫ|Ｒ_ｍ^０>で、期待値を最小化する方向でλを決定。λは最小値λ_ｍｉｎ、最大値λ_ｍａｘの間に制限する。|φ_ｍ^１>に対し、|Ｒ_ｍ^１>を計算する（|φ_ｍ^１>は規格化しておく）。
ここで、ループの打ち切り条件の判定を行う。期待値<φ_ｍ^１|Ｈ|φ_ｍ^１>の前回の値（<φ_ｍ^０|Ｈ|φ_ｍ^０>）からの低下値がある値以下、または、残差のノルムの二乗<Ｒ_ｍ^１|Ｒ_ｍ^１>の初期値<Ｒ_ｍ^０|Ｒ_ｍ^０>との比がある値以下、のどちらかが満たされれば打ち切る。

２巡目は、|φ_ｍ^０>、|φ_ｍ^１>と|Ｒ_ｍ^０>、|Ｒ_ｍ^１>に対して２×２のDIIS過程を行う。新しい波動関数|φ_ｍ^１>に対し残差ベクトルを求め、trial step|φ_ｍ^２>＝|φ_ｍ^１>＋λＫ|Ｒ_ｍ^１>を行う。λは１巡目と同じものを用いる。|φ_ｍ^２>に対し|Ｒ_ｍ^２>を計算（|φ_ｍ^２>は規格化）。
ここで、期待値<φ_ｍ^２|Ｈ|φ_ｍ^２>の前回からの低下値の打ち切り条件、または、残差のノルム二乗<Ｒ_ｍ^２|Ｒ_ｍ^２>の初期値との比の打ち切り条件の判定。

３巡目は、|φ_ｍ^０>、|φ_ｍ^１>、|φ_ｍ^２>と|Ｒ_ｍ^０>、|Ｒ_ｍ^１>、|Ｒ_ｍ^２>に対して３×３のDIIS過程を行う。新しい波動関数|φ_ｍ^２>に対し残差ベクトルを求め、trial step|φ_ｍ^３>＝|φ_ｍ^２>＋λＫ|Ｒ_ｍ^２>。λは１巡目と同じ。|φ_ｍ^３>に対し|Ｒ_ｍ^３>を計算（|φ_ｍ^３>は規格化）。
同様に期待値<φ_ｍ^３|Ｈ|φ_ｍ^３>、または、残差のノルム二乗<Ｒ_ｍ^３|Ｒ_ｍ^３>での打ち切り条件の判定。

４巡目は、|φ_ｍ^０>、|φ_ｍ^１>、|φ_ｍ^２>、|φ_ｍ^３>と|Ｒ_ｍ^０>、|Ｒ_ｍ^１>、|Ｒ_ｍ^２>、|Ｒ_ｍ^３>に対して４×４のDIIS過程。新しい波動関数|φ_ｍ^３>に対し残差ベクトルを求め、trial step|φ_ｍ^４>＝|φ_ｍ^３>＋λＫ|Ｒ_ｍ^３>。λは１巡目と同じ。|φ_ｍ^４>は規格化。

各状態毎のDIIS法のループは、最大４巡目で終了とする。また、占有数が空の状態は２巡目で打ち切る。この間、状態間の直交化処理は行わないので、各状態のDIIS法を各CPUで独立に行うことができ、効率的なバンド間並列が可能である。全体の状態で一巡した後、k点毎の各状態間の規格直交化処理（Gram-Schmidtの直交化）をまとめて行う。その後、k点毎の状態間でsubspace対角化を行う。直交化とsubspace対角化は、k点毎の並列化は可能だが、バンド間並列は不可能である。
これらの後、全ての状態の占有数の更新とchargeの更新、ハミルトニアンのポテンシャルの更新があり、はじめに戻って全状態のDIIS法を再び行う。この全体のループの一巡を１ステップとする（後述のテスト参照）。
以上は、基本的にRef.9に従ったものだが、Ref.9では、subspace対角化をループのはじめに行うなど少し異なっている。
Bylanderらの共役勾配法^６）との違いは、内側のループの中味と直交化をまとめて行う点だけである。残りの部分は基本的に同じプログラムが使える。各巡目でのハミルトニアンＨ×波動関数φの演算の回数（計算時間の目安になる）は各２回（１巡目のみ３回）で、共役勾配法と変わらない。しかし、直交化の総回数は、同じk点の他状態との直交化を内側のループの各巡目で複数回行わねばならない共役勾配法に比べ、DIIS法では数分の一以下になる。大きなサイズ（平面波基底数が大きい）ほど直交化処理の計算時間に占める割合は増えるので、DIIS法の方が有利になる。

３．２　Preconditioning関数と混合比λの計算
Preconditioningは、以下のような考え方から出てくる。現在の波動関数|φ_ｍ^０>と真の波動関数|φ_ｍ>とが

|φ_ｍ>＝|φ_ｍ^０>＋δφ_ｍ          (9)

の関係にあるとき、|φ_ｍ^０>＝|φ_ｍ>－δφ_ｍで、残差は

|Ｒ_ｍ^０>＝(Ｈ－ε_ｍ^０)|φ_ｍ^０>＝－(Ｈ－ε_ｍ^０)δφ_ｍ          (10)

となる（近似的に(Ｈ－ε_ｍ^０)|φ_ｍ>＝０）。とすると、現在の|φ_ｍ^０>に与えるべき変位δφ_ｍは、

δφ_ｍ＝－(Ｈ－ε_ｍ^０)^－１|Ｒ_ｍ^０>          (11)

となる。
Teterらは、この－(Ｈ－ε_ｍ^０)^－１（通常の最適化法のヘッセ行列^１５）の逆行列に対応）を対角行列Ｋで以下のように代用することを提案した^５）（平面波基底での表現）。

Ｋ_ＧＧ’＝－δ_ＧＧ’×<
(27+18x+12x^２+8x^３)／(27+18x+12x^２+8x^３+16x^４)          (12)
x=(ħ^２/2m)|k+G|^２／Ｔ          (13)
Ｔ＝<φ_ｍ|-(ħ^２/2m)∇^２|φ_ｍ>          (14)

Ｔは、その時点の波動関数φ_ｍの運動エネルギー期待値である。Bylanderらの共役勾配法でも同じものが用いられている。なお、Teterらの論文では、勾配ベクトル（残差ベクトルにマイナスをかけた形）に作用させる形なのでマイナスがついていないことに注意。
一方、Kresse-Furthmuller^９）は、以下の改良形を提案している。

Ｋ_ＧＧ’＝－δ_ＧＧ’×（4/(3Ｔ)）×
(27+18x+12x^２+8x^３)／(27+18x+12x^２+8x^３+16x^４)          (15)
x=(ħ^２/2m)|k+G|^２／(3Ｔ/2)          (16)

Ｔは(14)と同じである（Ref.9では、Ｔが残差ベクトルの運動エネルギー期待値と書かれているが、誤りと思われる）。
(3)式のtrial stepでのλの決定法は、以下と推定される。まず、期待値をλの２次までで展開する。

Ｅ[λ]＝<φ_ｍ^１|Ｈ|φ_ｍ^１>／<φ_ｍ^１|φ_ｍ^１>          (17)
＝Ｅ[０]＋Ａλ＋Ｂλ^２          (18)

Ａ，Ｂが求まれば、Ｅ[λ]を最小にするλ＝－Ａ／2Ｂが求まる。係数Ａ，Ｂを決めるための条件式を組み立てよう。(2)式からＥ[０]＝ε_ｍ^０である。適当なλ_１で(3)式から実際に波動関数を構築し、(17)式からＥ[λ_１]の値が計算できる。これは、(18)式からＥ[０]＋Ａλ_１＋Ｂλ_１^２に等しい。次に、Ｅ[λ]のλでの微分を考える。λ＝０での微分値は、波動関数の微分を用いた計算から、

∂Ｅ[λ]/∂λ|_λ＝０＝｛<ＫＲ_ｍ^０|Ｈ|φ_ｍ^０>
－<φ_ｍ^０|Ｈ|φ_ｍ^０><ＫＲ_ｍ^０|φ_ｍ^０>｝＋｛複素共役｝          (19)

となり、値が計算できる。これは、(18)式からＡに等しくなる。以上の条件式からＡ，Ｂが求まり、λが決まる。
λ_１は0.5を用いた。実際に求まるλは0.1～10程度の範囲で、λ_１にはそれほど依存しなかった。

３．３　Ｓを含む固有値方程式の解法
実際のプログラミングは、従来のBylanderらの共役勾配法^６）のプログラムの修正として比較的容易に行える。しかし、DIIS過程の基底間の重なりＳを含む固有値方程式(6)の解法は注意が必要である。小規模なので、通常の数学ライブラリのエルミート行列対角化のサブルーチンを用いて解くが、量子化学の教科書^１６）にあるように、まず、Ｓ行列の固有値・固有ベクトルを求め、それらをs_ｉ、{Ｕ_１ｉ, Ｕ_２ｉ, Ｕ_３ｉ, , , }として、変換行列Ｘは、s_ｉ^１／２{Ｕ_１ｉ, Ｕ_２ｉ, Ｕ_３ｉ, , , }を列ベクトルとする行列である。次にＲ行列をＸで変換したＲ'＝Ｘ^†ＲＸの固有値・固有ベクトルを求める。その固有ベクトルに行列Ｘを作用させれば最終的な解{α_ｉ}が求まる。
このとき、Ｓ行列の固有値に負やゼロがあると、変換行列Ｘの構築時に取り除いて行列の次数を下げて進める。問題は、DIIS法では、収束に近くなると、どうしてもＳの固有値にゼロに近いものが生じ、うまく取り除かないとゼロ割が発生し、また、DIIS過程そのものが意味をなさなくなることである。
例えば、収束に近くなると、DIIS過程の基底{|φ_ｍ^ｉ>}はどれも同じものに近くなり、Ｓ行列はどの要素も１に近くなる。そうすると、Ｓの固有値は、２×２の場合、０、２、３×３の場合、０、０、３、４×４の場合、０、０、０、４に近づく。ゼロを取り除いてＸ行列を組み立てて解いた解{α_ｉ}からの波動関数(4)式は、結局、基底の均等な混合となる。そうなると、DIIS過程を行うことでかえって波動関数が収束から遠ざかってしまう（最新の基底|φ_ｍ^Ｍ>が最も良いはずなので）。
今回は、Ｓの固有値が10^－１４以下でゼロと見なし、上記のように均等な混合になる場合、DIIS過程の出力の(4)式を最新の|φ_ｍ^Ｍ>で置き換え、次のtrial stepに進むようにした。

４　計算実行時の諸条件

実際の計算の実行についは、Ref.9に習って、(i)初めの数回のステップは、共役勾配法^６）で行い、その後にDIIS法を行った。(ii)また、半導体の場合も金属の場合も占有される状態に加えて上の空の状態を幾つか用意した。これらは、DIIS法の性質上、必ずしも全ての固有状態を自然に下から全て覆うわけではないこと、変分の自由度を増す目的などからである。
検討すべき計算条件やパラメータは、（１）trial stepでのλの制限範囲、特に最大値λ_ｍａｘ（最小値λ_ｍｉｎはRef.9と同様0.1とした。Ref.9ではλ_ｍａｘは1.0を用いている）、（２）DIIS法ループの打ち切り条件：期待値の前回値からの低下の値、（３）DIIS法ループの打ち切り条件：残差ベクトルのノルムの二乗値の初期値に対する比（Ref.9では0.3としている）、（４）Preconditioning関数としてTeterらのものを用いるか、Kresseらのものを用いるか、（５）DIIS法のループで１巡目に期待値での打ち切り判定を行うかどうか、等である。
最後の（５）の条件は、１巡目に期待値の判定を入れてしまうと、収束に近づくとほとんど１巡目でループが打ち切られることが起き（後述）、実質的に共役勾配法等と同様の方法になってしまうためである（２巡目に行かない限りDIIS過程は行わないので）。
今回、比較対照とするBylanderらの共役勾配法では、期待値の低下が10^－８Ry以下、あるいは最初の低下値の30％以下になることを各状態のループの打ち切り条件とした。全体のchargeの更新を含めたステップの収束条件は、chargeのself-consistencyにつき、メッシュ点のinputとoutputの差の最大値が10^－７未満、且つステップでのエネルギー低下が10^－８Ry以下である。
DIIS法での収束テストについて、ステップの収束条件は共役勾配法と同じとした。DIISループの打ち切り条件については、（２）の期待値の低下値は共役勾配法と同様に10^－８Ryとした。残りの（１）、（３）、（４）、（５）については、各種のテストを行い、最適条件を探った。
今回の収束テストでは、初期のcharge分布を中性原子のchargeの重ね合わせから作り、初期の波動関数は乱数で作った。chargeの更新は、simple-mixing法かKerker-mixing法^１４）を用いた。後者の場合、電子密度分布のフーリエ成分ρ(G)に対して、

ρ_ｉｎ^ｎｅｗ(G)＝ρ_ｉｎ^ｏｌｄ(G)
＋(ＡG^２/(G^２+G_０^２))(ρ_ｏｕｔ^ｏｌｄ(G)-ρ_ｉｎ^ｏｌｄ(G))          (20)

により、新しいinput chargeが構築される。ρ_ｉｎ^ｏｌｄ(G)、ρ_ｏｕｔ^ｏｌｄ(G)は、各々前回の入力と出力の電子密度分布のフーリエ成分で、ρ_ｉｎ^ｎｅｗ(G)は新しい入力の電子密度分布のフーリエ成分である。G_０^２は、Gの大きさ毎に混合比(ＡG^２/(G^２+G_０^２))を変えるためのパラメータで、Gが大きいと混合比はＡに近づき、小さいほど混合比は小さくなる。これは、振動の原因になる小さなGの混合比を小さくするための方法であり、G_０^２＝０なら混合比は全てのGでＡであり、simple-mixingになる。

５　完全結晶でのテスト

完全結晶のテストでは、（５）の条件については、１巡目では期待値の判定は行わないモードを用いた。

図２は、Si完全結晶セルのテスト結果である。chargeの更新はsimple-mixing（A=0.8）を用いた。DIIS法は共役勾配法（図のCG(BKL)）と同等の効率である。TeterらのPreconditioning関数を用いた場合、残差のノルム二乗値の比の打ち切り条件が0.2で、λ_ｍａｘ＝3.0（図のDIIS-2）、あるいはλ_ｍａｘ＝2.0のもの（図のDIIS-1）が他のパラメータ値の結果に比して良好であった。KresseらのPreconditioning関数もテストしてみた（DIIS-3、λ_ｍａｘ＝3.0、ノルム二乗値の比0.2）。

Fig. 2　Test of the RMM-DIIS scheme for the Si crystal (two atom cell).　
Ｅ_ｃｕｔ=20Ry and 10 k points per irreducible zone are used.　
Six bands for each k point are treated.　
The conjugate- gradient (CG) scheme is used for the first two iterations.

Al完全結晶セルのテスト結果もSiの場合と同様にDIIS法と共役勾配法でそれほど差はなかった。しかし、図３のfcc-Ti結晶では、DIIS法の方が収束が速かった。これはcharge-mixing法も関係すると考えられるが、共役勾配法で探索した範囲では、Kerker法でA=0.8、G_０^２=1.5がベストで、DIIS法も同じmixing法を用いている。TeterらのPreconditioning関数を用いた場合、λ_ｍａｘ＝3.0、残差のノルム二乗値の比の条件0.1が速かった（DIIS-1）。Kresseらのものを用いた場合、λ_ｍａｘ＝1.0、ノルム二乗値の比の条件0.1（DIIS-2）、あるいはλ_ｍａｘ＝2.0、ノルム二乗値の比の条件0.1（DIIS-3）が速かった。全体で は、KresseらのPreconditioning関数の方が良い。DIIS法の方が収束が速い理由は、後述のTi表面の場合と同じと考えられる。

Fig. 3　Test of the RMM-DIIS scheme for the fcc Ticrystal (one atom cell).　
Ｅ_ｃｕｔ=40Ry and 20 k points per irreducible zone are used.　
Five bands for each k point are treated with the Gaussian broadening scheme^１７）(σ=0.35eV).　
The CG scheme is used for the first five iterations.

６　表面でのテスト

比較的大きな系として表面スーパーセルのテストを行った。図４はSi{110}表面の結果である（Kerker法パラメータA=1.0、G_０^２=0.6）。KresseらのPreconditioning関数について、Ti結晶のテストで良好であった条件（残差のノルム二乗値の比の条件0.1で、λ_ｍａｘ＝1.0、あるいはλ_ｍａｘ＝2.0）をテストした。いづれも共役勾配法とあまり変わらなかった。ただし、上記（５）の条件、すなわちDIIS法のループで１巡目に期待値での打ち切り判定を行うかどうかで、（結晶のテストと違って）１巡目に判定を行うようにした方が収束が若干速く、計算時間も若干短い。期待値の低下条件を１巡目の打ち切り判定に入れた場合、収束に近づくと、DIIS法の各状態のループはほとんど１巡目で打ち切られ、共役勾配法とあまり変わらなくなる（DIIS過程を行わない）。しかし、Si表面は共役勾配法で良好に収束するので、問題はない。むしろ、無理やり２巡目のDIIS過程をやらせる方が計算時間が無駄なようである（ノルム二乗値の比の判定値を0.1と小さくしているので、ループの打ち切りはほとんど期待値の低下条件で起きる）。図４のDIIS-1は、１巡目の打ち切り判定に期待値を入れた場合のλ_ｍａｘ＝1.0の結果である。

Fig. 4　Test of the RMM-DIIS scheme for the Si {110} surface slab
(five atomic layers and vacant region corresponding to three atomic layers).　
Ｅ_ｃｕｔ=18Ry and 4 k points per irreducible zone are used.　
26 bands for each k point are treated with the Gaussian broadening scheme (σ=0.2eV).　
The CGscheme is used for the first three iterations.

ところが、Ti表面のテストでは、（５）の条件で、１巡目に期待値での打ち切り判定を行わないようにして、収束に近づいてもなるべく２巡目のDIIS過程をやらせるようにした方が収束が速いことが判明した。
図５に結果を示す。もともと共役勾配法でかなり収束の遅い系である（Kerker法パラメータA=0.8、G_０^２=2.0）。収束に近づくと、chargeの小さい振動のため収束スピードが急に落ちる（図のCG(BKL)のプロット）。こういう現象は、遷移金属／セラミックス界面などの大きなセルで頻繁に見られ、従来、対策に苦労してきた^１８）。図のDIIS-3は、図４のSi表面のDIIS-1と同じ条件での結果で、１巡目に期待値での打ち切り判定を行うようにしている。そうすると、収束に近づくと１巡目で打ち切られることが多く、ほとんど共役勾配法と同じ過程になり、共役勾配法と同様、収束スピードの鈍化が見られる。一方、図のDIIS-1、DIIS-2は、１巡目に期待値での打ち切り判定を行わないようにし、毎回２巡目のDIIS過程をなるべくやらせたもので、収束スピードの鈍化が押さえられている。DIIS-1、DIIS-2は、KresseらのPreconditioning関数で、残差のノルム二乗値の比の条件0.1、各々λ_ｍａｘ＝1.0とλ_ｍａｘ＝2.0である。

Fig. 5　Test of the RMM-DIIS scheme for the Ti{110} surface slab
(five atomic layers and vacant region corresponding to three atomic layers).
Ｅ_ｃｕｔ=30Ry and 8 k points per irreducible zone are used.　
18 bands for each k point are treated with the Gaussian broadening scheme (σ=0.3eV).　
The CG scheme is used for the first five iterations.

charge-mixing法の改良も重要と考えられるが、少なくとも、共役勾配法で問題になってきた遷移金属を含む大きなセルでの収束スピードの鈍化の問題^１８）は、DIIS法により克服の可能性がある。この点は、これまでのDIIS法の各種テスト^９）でも見出されていない点である。共役勾配法で収束が遅くなる原因として、共役勾配法では波動関数の配位空間でのエネルギー曲面の極小値（固有関数）に向けて、基本的に「勾配」に従って低い方へ波動関数を変化させるため、極小点付近の曲面の曲率が小さい場合、波動関数がポテンシャルと合わせて極小点付近を振動することが考えられる。DIIS法では、「勾配」や曲率に関わらず、数値的に強制的に極小点（勾配ゼロ＝残差ゼロの地点）に持っていくので、そうした振動が起きにくいと考えられる。

７　まとめ

DIIS法のプログラムを構築し、各種の計算条件、パラメータの検討を行った。従来の共役勾配法と同等の収束性が得られ、特に遷移金属を含む大きなセルでは、従来法での収束の鈍化の問題を克服できる可能性が示された。なお、並列化プログラムのテスト、charge-mixingや格子緩和へのDIIS法の適用も検討中である。

参考文献
1) 藤原毅夫, "固体電子構造", (1999) 朝倉書店.
2) P.Hohenberg and W.Kohn, Phys.Rev., 136, B864(1964); W.Kohn and L.J.Sham, ibid., 140, A1133(1965).
3) J.R.Chelikowsky and M.L.Cohen, "Handbook on Semiconductors, Vol.1", edited by P.T.Landsberg, p.60 (1992) Elsevier.
4) R.Car and M.Parrinello, Phys.Rev.Lett., 55, 2471 (1985).
5) M.P.Teter et al., Phys.Rev.B, 40, 12255(1989); M.C.Payne et al., Rev.Mod.Phys., 64, 1045 (1992).
6) D.M.Bylander et al., Phys.Rev.B, 42, 1394 (1990).
7) D.Singh, Phys.Rev.B, 40, 5428 (1989).
8) P.Pulay, Chem.Phys.Lett., 73, 393 (1980).
9) G.Kresse and J.Furthmuller, Phys.Rev.B, 54, 11169 (1996).
10) M.Kohyama, Model.Simul.Mater.Sci.Eng., 4, 397 (1996); 「共役勾配法」を参照
11) P.Pulay, J.Comp.Chem., 3, 556 (1982); P.Csaszar and P.Pulay, J.Mol.Struc., 114, 31 (1984).
12) D.M.Wood and A.Zunger, J.Phys.A, 18, 1343(1985).
13) J.Hutter et al., Comp.Mater.Sci., 2, 244(1994).
14) G.P.Kerker, Phys.Rev.B, 23, 3082 (1981).
15) 大野豊・磯田和男監修, "新版数値計算ハンドブック", p.792 (1990) オーム社.
16) A.Szabo and N.S.Ostlund, 大野公男他訳, "新しい量子化学", 上巻 p.156 (1987) 東大出版会.
17) C.L.Fu and K.M.Ho, Phys.Rev.B, 28, 5480(1983).
18) S.Tanaka and M.Kohyama, Phys.Rev.B, 64, 235308 (2001).