JP5447674B2 - 電子状態計算方法、電子状態計算装置及びコンピュータプログラム - Google Patents

Description

本発明は、物質の電子状態を計算により求める電子状態計算方法、電子状態計算装置及びコンピュータプログラムに関する。

従来、物質の物理的又は化学的性質（以下、物性と呼ぶ）を量子力学の基礎法則に則って予測する第一原理計算理論と呼ばれる計算理論が存在する。この計算理論のうち、弾性を含む力学特性、超伝導特性を含む伝導特性、誘電性、磁性などの物性を近似的に再現し、計算規模が実施可能な範囲である、密度汎関数理論に基づく計算理論が存在する。この計算理論は、既に物質設計に応用され、その予測精度及び精度限界が実験で実証された例が多く存在する（例えば、非特許文献１を参照）。

密度汎関数法に基づく計算理論には、局所密度近似（LDA : Local Density Approximation）による自己無撞着計算理論が含まれるが、計算結果は必ずしも実験と同一の物性とはならないことが知られている。

密度汎関数法に基づく計算手法には、ＬＤＡにおける問題点を克服する手法として位置づけられる一般化勾配近似（GGA : Generalized Gradient Approximation）、ＧＷ近似、ＧＷ＋Γ近似、ＬＤＡ＋Ｕ近似など複数の近似計算手法が知られているが、これらの近似計算手法についても、計算結果は必ずしも実験と同一の物性とはならないことが知られている。

「計算機マテリアルデザイン入門」、笠井秀明、赤井久純、吉田博編、２００５年、大阪大学出版会

国際公開第２００７／１４１９４２号パンフレット 特願２００８−２２３８１３号明細書

前述した近似計算手法は、いずれも真の解（厳密解）が示す物性に有限回数の計算で常に到達する方法を与えていないという問題点を有していた。すなわち、ＬＤＡが与える高精度の一電子基底を活用しながら、常に真の解の物性を再現することが保障されている自己無撞着計算理論は、存在していなかった。

一方、本願発明者の一名は、多配置参照密度汎関数法に基づく有効多電子計算として揺らぎ参照決定法を提案している（例えば、特許文献１を参照）。この方法により、一電子密度に加えて、局所密度揺らぎなどの正定値性を備えた揺らぎ変数（相関関数）を再現する拡張コーン・シャム方程式を与えて、物理量の再現精度を任意の精度で高めることが可能である。この計算手法では、ハミルトニアンの最低エネルギ状態の決定を通して、物質の全エネルギ、一電子密度、特定した揺らぎ変数を同時に再現した。更に、拡張コーン・シャム方程式が定めるモデルとその安定解を初期値として用いることで、全エネルギ、一電子密度、正準相関関数を再現する時間発展方程式を与えて物質の外場応答（力学変形、電磁場印加に対する応答など）を再現するための、高精度な第一原理電子状態計算の方法を与えた（特許文献１を参照）。

しかしながら、前記の多配置参照密度汎関数法に基づく有効多電子計算は、量子モンテカルロ法、トランスコリレイティッド法、配置間相互作用法、摂動計算・グリーン関数法、有効ポテンシャル法などの参照計算が不可避な計算手法であった。そのため、物理量の再現精度は参照計算の計算精度に依存することとなり、参照計算の計算精度を越えて真の解に到達する計算手法を与えるものではなかった。

そこでさらに、本願発明者の一名は、密度汎関数を用いる量子力学的変分計算理論（密度汎関数変分法）及び厳密解を含むバナッハ空間となるモデル空間中で厳密解に到達する近似モデルの列を与える方法を与え、計算機資源の許す範囲内で安定した数値解を求めることで、厳密解の評価ができる電子状態計算方法、電子状態計算装置、コンピュータプログラム、及び記録媒体を提供した（特許文献２を参照）。

しかしながら、前記モデル空間中で厳密解に到達する近似モデルの列により厳密解に到達する方法は、物質の物理的又は化学的性質のみを再現することに目的を限定すれば、より少ない有限回数の計算で厳密解の物性に到達することを用いてはいなかった。

本発明は斯かる事情に鑑みてなされたものであり、局所密度近似法、一般化勾配近似法、ハートレーフォック法、配置間相互作用法、又は摂動ＭＰ（Moeller-Plesset）法などの既存の近似モデルを出発点として、モデル空間中で与えられる密度汎関数変分法の変分原理に従って、複数の演算経路を探索しながら厳密解が示す物性に有限回の計算手続きで到達することができる電子状態計算方法、電子状態計算装置及びコンピュータプログラム提供することを目的とする。

本発明に係る電子状態計算方法は、計算装置を用いて物質の電子状態を計算する方法において、前記計算装置は、物質の電子状態に対する近似解を与える複数の演算モデルの夫々を構成要素として含む集合を設定し、複数の電子軌道上に存在する少なくとも２つ以上の電子を含む電子系に働く有効相互作用を含む有効ハミルトニアンを、各演算モデルを特定するために定め、前記有効ハミルトニアンの自己無撞着解を前記集合内の各演算モデルを用いて計算する過程で、計算した自己無撞着解同士が連続的に変化する方向を特定しながら、前記集合がなす空間内で距離が近接する複数の演算モデルのうちで最適な演算モデルを量子力学的変分法に基づいて定め、前記最適な演算モデルを逐次更新する際に、前記有効ハミルトニアンの自己無撞着解による電子系の変分エネルギを評価し、評価した変分エネルギが計算すべき厳密解のエネルギに近接してゆき、しかも変分エネルギが単調減少凸関数をなすように演算モデルを更新して、一又は複数の変分エネルギの系列から前記電子状態の厳密解を計算することを特徴とする。

本発明に係る電子状態計算方法は、局所密度近似法、一般化勾配近似法、又はハートレーフォック法を基に、有効相互作用として揺らぎを示す非局所型演算子を定める手続きを通じて前記有効ハミルトニアンを決定することを特徴とする。

本発明に係る電子状態計算方法は、前記揺らぎが、クーロン相互作用とハートレー平均場項のずれである密度揺らぎを含むことを特徴とする。

本発明に係る電子状態計算方法は、前記厳密解に近接してゆく変分エネルギの系列を与える一又は複数の演算モデルの系列のうち、電子状態を示す相が転移せずに厳密解と同一視される相を与え、かつ最小の計算ステップで前記有効ハミルトニアンの自己無撞着解に至る演算モデルを含むものを定めることを特徴とする。

本発明に係る電子状態計算方法は、各演算モデルによる自己無撞着解をＬＡＰＷ法（Linearized Augmented Plane Wave 法）、ＰＡＷ法（Projector Augmented Wave 法）、又は数値基底展開法を用いた並列計算により求めることを特徴とする。

本発明に係る電子状態計算装置は、物質の電子状態を計算する装置において、物質の電子状態に対する近似解を与える複数の演算モデルの夫々を構成要素として含む集合を設定する手段と、複数の電子軌道上に存在する少なくとも２つ以上の電子を含む電子系に働く有効相互作用を含む有効ハミルトニアンを、各演算モデルを特定するために定める手段と、前記有効ハミルトニアンの自己無撞着解を前記集合内の各演算モデルを用いて計算する過程で、計算した自己無撞着解同士が連続的に変化する方向を特定しながら、前記集合がなす空間内で距離が近接する複数の演算モデルのうちで最適な演算モデルを量子力学的変分法に基づいて定める手段と、前記最適な演算モデルを逐次更新する際に、前記有効ハミルトニアンの自己無撞着解による電子系の変分エネルギを評価する手段と、評価した変分エネルギが計算すべき厳密解のエネルギに近接してゆき、しかも変分エネルギが単調減少凸関数をなすように演算モデルを更新する手段と、一又は複数の変分エネルギの系列から前記電子状態の厳密解を計算する手段とを備えることを特徴とする。

本発明に係る電子状態計算装置は、局所密度近似法、一般化勾配近似法、又はハートレーフォック法を基に、有効相互作用として揺らぎを示す非局所型演算子を定める手続きを通じて前記有効ハミルトニアンを決定するようにしてあることを特徴とする。

本発明に係る電子状態計算装置は、前記揺らぎが、クーロン相互作用とハートレー平均場項のずれである密度揺らぎを含むことを特徴とする。

本発明に係る電子状態計算装置は、前記厳密解に近接してゆく変分エネルギの系列を与える一又は複数の演算モデルの系列のうち、電子状態を示す相が転移せずに厳密解と同一視される相を与え、かつ最小の計算ステップで前記有効ハミルトニアンの自己無撞着解に至る演算モデルを含むものを定めることを特徴とする。

本発明に係る電子状態計算装置は、各演算モデルによる自己無撞着解をＬＡＰＷ法（Linearized Augmented Plane Wave 法）、ＰＡＷ法（Projector Augmented Wave 法）、又は数値基底展開法を用いた並列計算により求めるようにしてあることを特徴とする。

本発明に係るコンピュータプログラムは、コンピュータに、物質の電子状態を計算させるコンピュータプログラムにおいて、物質の電子状態に対する近似解を与える複数の演算モデルの夫々を構成要素として含む集合を設定させるステップと、複数の電子軌道上に存在する少なくとも２つ以上の電子を含む電子系に働く有効相互作用を含む有効ハミルトニアンを、各演算モデルを特定するために定めさせるステップと、前記有効ハミルトニアンの自己無撞着解を前記集合内の各演算モデルを用いて計算させる過程で、計算させた自己無撞着解同士が連続的に変化する方向を特定しながら、前記集合がなす空間内で距離が近接する複数の演算モデルのうちで最適な演算モデルを量子力学的変分法に基づいて定めさせるステップと、前記最適な演算モデルを逐次更新させる際に、前記有効ハミルトニアンの自己無撞着解による電子系の変分エネルギを評価させるステップと、評価させた変分エネルギが計算すべき厳密解のエネルギに近接してゆき、しかも変分エネルギが単調減少凸関数をなすように演算モデルを更新させ、一又は複数の変分エネルギの系列から前記電子状態の厳密解を計算させるステップとをコンピュータに実行させることを特徴とする。

本発明にあっては、計算装置を用いて物質の電子状態を計算する方法、装置、コンピュータプログラムにおいて、物質の電子状態が示す物性を再現する有限回数の計算手続きを多配置参照密度汎関数理論に基づいて与える。

本発明にあっては、多配置参照密度汎関数理論による演算モデルが与える決定方程式の自己無撞着解を決定することを通して、物質の厳密解が示す電子状態がもつ物性を、計算装置が自動的に決定する。

本発明にあっては、厳密解が示す物性の再現に必要となる局所密度揺らぎなどの正定値性を備えた揺らぎ変数（相関関数）又は多重秩序変数が有限個であることが、密度汎関数理論における原理により保証される。

本発明による場合は、物質の電子状態に対する近似解を与える演算モデルを出発点とし、電子間に働く有効相互作用を示す有効ハミルトニアンの自己無撞着解を求める過程で、モデル間の距離が接近していくように、計算過程で使用する演算モデルを決定する。これにより、新しい演算子空間の最適化手続きが与えられ、繰り込み思想に基づくダウンフォールディングの方法では不可能であった上位モデルからなる系列を形成することで、厳密解が与える物性の評価方法が有限回の計算手続きにより与えられる。

本発明による場合、多種多様の演算モデルを含むモデル空間内で、厳密解に到達しうる複数の経路を与えることができるため、電子相関を取込む経路を適宜設定することができる。例えば、強相関効果をフェルミ準位から順次取り込み、エネルギスケール又は空間スケールに依存した全ての相関を考慮することでクーロン相互作用に漸近する演算子を導入することができる。

本発明による場合は、自己無撞着解を求める演算を並列して実行することが可能であるため、並列計算を実現することにより計算速度が向上する。また、計算の各ステップで密度汎関数法を含む既存の第一原理計算法で実現されている並列化計算技法を用いることで、計算速度が向上する。

本発明による場合は、局所密度近似法、一般化密度勾配近似法、ハートレーフォック法、配置間相互作用法、及び摂動ＭＰ法を含む既存の演算モデルを出発点として厳密解が示す物性に有限回の計算ステップで到達する手法を与えることができる。

本発明による場合は、有効ハミルトニアンの表式が、フェルミ準位近傍の電子間のクーロン相互作用に漸近する有効相互作用を演算子として含む手法が与えられる。

本発明による場合は、各演算モデルによる自己無撞着解をＬＡＰＷ法（Linearized Augmented Plane Wave 法）、ＰＡＷ法（Projector Augmented Wave 法）、又は数値基底展開の方法などを用いて、基底展開を行った高精度な一電子波動関数の表現方法を用いて求める手法が与えられる。

本発明による場合は、モデル系とその安定解を初期値として用いることで、全エネルギ、一電子密度、正準相関関数を再現する時間発展方程式をさらに与え、物質の外場応答（力学変形、電磁場印加に対する応答など）を再現する高精度な第一原理電子状態計算の方法を与えることができる。

本発明による場合は、与えられた計算機資源を利用する範囲内で効率良く厳密解が示す物性に到達し得る計算手続きを決定することができる。
その結果、物質設計を計算機シミュレーションによる量子シミュレータを用いて実施する新しい量子デザイン手法において、スピン・エレクロニクス、分子・エレクトロニクスなどの量子素子形成要素の設計法の提供、エレクトロニクス応用によるエネルギ問題回避法の提供、低環境負荷又は戦略的元素選択による環境問題・資源エネルギ問題回避型物質解の提供、新型センサー設計、生体親和材料設計、薬剤設計などの計算装置による物質設計が可能となる。

計算原理を説明する説明図である。 本実施の形態に係る電子状態計算装置の内部構成を示すブロック図である。 電子状態計算装置が実行する処理の手順を示すフローチャートである。 変分エネルギの評価結果の一例を示す図である。 最適モデル決定のプロセスを与える原理を説明する説明図である。 電子状態計算装置が最適モデル決定のために実行する処理の手順を示すフローチャートである。

以下、本発明をその実施の形態を示す図面に基づいて具体的に説明する。
実施の形態１．
本実施の形態に係る電子状態計算方法における計算原理は以下の通りである。本実施の形態では、電子状態の厳密解が示す物性を再現する最適化された演算モデルを決定する方法として、クーロン多体系が基底状態で発生する電子密度なる秩序変数を数値的に求めることを通して、物性の再現を行うことを目標とする。電子密度は実験で観測できる物理量であり、存在することが知られている。クーロン多体系における基底状態のエネルギＥ₀ は以下の表式を満たすため、数値的に評価を行うことができる。

（数式１）

エネルギ汎関数Ｇバー_{Xi, εi, gi}は、次の式により定義される。

（数式２）

また、ΔＥバー_{Xi, εi, gi}は、次の式により定義される。

（数式３）

ここで、Ｇバー_{Xi, εi, gi}は演算モデルが表現する電子の波動関数Ψに関するエネルギ汎関数を表し、ΔＥバー_{Xi, εi, gi}はクーロン系の変分エネルギを評価する際に発生するエネルギのずれを補正するエネルギ汎関数を表す。

ここで、Ψは多粒子波動関数、Ｔは運動エネルギに係る演算子、Ｖ_red ^Xi は以下の数式４によって表される作用素である。なお、数式中では、演算子Ｔ及び作用素Ｖ_red ^Xi はハット付の文字で表記しているが、明細書中においてはそれぞれハットを付さずに表記する。また、数式２において、ｅは電子の電荷、ｎ_Ψ （ｒ）は位置ベクトルｒでのΨが与える電子密度、v_ext （ｒ）は外部スカラーポテンシャルを表す。

（数式４）

ここで、Ｘ_i はモデルを定めるパラメータと演算子との組からなる集合を表す。Ｐ⁽ⁱ⁾ は射影演算子である。また、：：はノーマルオーダと呼ばれる演算子順序を指定するものであり、生成演算子を消滅演算子の前にフェルミオン演算子の交換規則に沿って順序交換してから用いることを指定するものである。また、Ｙ_n ⁽ⁱ⁾、Ｚ_n ⁽ⁱ⁾は、一体波動関数の完全系が与える軌道ｌ上での電子の消滅演算子ｃ_lσ ⁽ⁱ⁾ から与えられる演算子である。ｆ_{n, +, l, σ}及びｆ_{n, -, l, σ}は複素定数である。α_n ⁽ⁱ⁾は自己相互作用補正を制御するパラメータである。Ξ_n ⁽ⁱ⁾は、演算子Ｙ_n ⁽ⁱ⁾、Ｚ_n ⁽ⁱ⁾に関する関数であり、期待値が有限な下界をもつことが分かっていることが望ましいが、特に多項式関数であってもよい。

また、数式２におけるＥ_εi ^local ［Ψ］及びＥ_gi ^non-local ［Ψ］は局所性をもつモデル・エネルギ汎関数及び非局所性をもつモデル・エネルギ汎関数であり、例えば、以下のように定義されるものを用いてよい。

（数式５）

ここで、ε_i （ｎ）は有界単調減少連続関数、ｇ_i （ｌ₁ ，ｌ₂ ，σ）は実数値係数であり、入力の際に与えられる。これらは計算プロセスの途中で定まるコーン・シャム方程式の解を用いて最適化しなおすこともできるが、充分な計算規模を確保できたことにより数式４を通して取り込まれる揺らぎが充分であれば結果には影響しない。

数式２，３，４及び５の有効相互作用によるモデルを導入していることが本願発明の特徴の１つとなる。演算子に関する関数Ξ_n ⁽ⁱ⁾及び有界単調減少連続関数ε_i （ｎ）及びｇ_i （ｌ₁ ，ｌ₂ ，σ）を最適化することで、計算資源の範囲内で、高効率で高速かつ高精度な計算を可能とする最適化モデル系列を探索することができるようになる。

数式４は、クーロン系の密度揺らぎをもとにチャンネル分解を行った有効多粒子相互作用を導入することによって、以下の式を具体的な表式として用いることが可能である。

（数式６）

数式１は、モデル空間において成立する変分原理を与えている。この原理は、Ｅ_gi ^non-local ［Ψ］が特別な場合に限れば証明が公開されている（例えば、K. Kusakabe, J. Phys. Soc. Jpn 78, 114716(2009)を参照）。そのため、不等式として数式１に類似の評価式が与えられる。等号を成立させる条件の存在は示されているが、有限の計算機資源による計算では厳密な等号は達成されない。しかしながら、変分エネルギに関する評価値の数列を得ることを通して、十分な計算結果の集積により基底状態のエネルギＥ₀ を数値的に評価することが可能となる（特許文献２を参照）。本出願では、この変分原理が広範なエネルギ汎関数に対して成立することを用いる。さらに、この計算過程を複数用意することによって、厳密解が示す物性を導出する有限回数の計算過程を複数与え、計算資源の範囲内で、高効率で高速かつ高精度な計算方法を見つけ出すことができる。

数式１を与える原理は次のように厳密に示される。残差交換相関エネルギ汎関数ΔＥ_{Xi, εi, gi}［Ψ］を、自己相互作用補正項を分離せずにクーロン揺らぎを与える表式を用いて、次式により導入する。

（数式７）

ここで、λに関する積分は、クーロン相互作用をλ倍に変化させて得られる制限付き探索法によるエネルギ汎関数（数式８を参照）のλ微分に関するルベーグ積分、及びF_λ [Ψ]のλディニ微分が有限個の飛びを生じている点での補正をもとに定義される量であり、存在が示されている。密度を変動させない範囲で、λディニ微分が飛びを生じさせる点において、残る相関による相転移が発生する。

（数式８）

このとき、数式１のΔＥバー_{Xi, εi, gi}［Ψ］を、数式７で定義されたΔＥ_{Xi, εi, gi}［Ψ］に置き換えることによって、以下の表式を得ることが可能となる。

（数式９）

この数式９は、Ψ_{Xi, εi, gi}が変分原理を満たすことを通して得られるが、さらに次の不等式の成立を保証する。

（数式１０）

図１は計算原理を説明する説明図である。数式１の表式に従い、まず、演算モデルを固定した上でエネルギ汎関数Ｇバー_{Xi, εi, gi}を最小化し、次に、最小化されたエネルギ汎関数Ｇバー_{Xi, εi, gi}と、この最小値を与える波動関数Ψ_{Xi, εi, gi}によって定められるΔＥバー_{Xi, εi, gi}との和である変分エネルギの評価値を得る。続いて、演算モデルを変化させて同じ計算ステップを実行することで、変分エネルギの評価値を複数個得る。複数の演算モデルの夫々を構成要素として含む集合に、電子密度間の絶対値ノルムを用いて距離を導入することで得られる空間内で、演算モデルがもつ相互作用パラメータを連続的に変化させることを用いて、近接したモデルの系列を得る。その結果、相互作用パラメータに対して変分エネルギの変化が作る漸近曲線を得ることができる。この漸近曲線の収束値を数値的に評価することで、基底状態のエネルギＥ₀ の評価値を求める。図１は、この方法が複数個存在することを表している。演算モデルとは、Ｇバー_{Xi, εi, gi}と、この最小値を与える波動関数Ψ_{Xi, εi, gi}により構成される。波動関数Ψ_{Xi, εi, gi}は、自己無撞着方程式を形成するＧバー_{Xi, εi, gi}の最小化問題を定める決定方程式に対して、自己無撞着性を備えるように決定することができる。
後述するように、あるモデル（例えば、図１のＸ_i に関するモデル）、及び他のモデル（例えば、図１のＸ_j に関するモデル）の収束性について判定する。このとき、モデル系列｛Ｘ_i ｝内でモデル列Ｘ_i が電子密度及び変分エネルギに関して収束しているか否かを判定する。また、モデル系列｛Ｘ_j ｝内のモデル列Ｘ_j についても、電子密度及び変分エネルギに関して収束しているか否かを判定する。本実施の形態では、２種類の演算モデルを含むモデル列が共通の収束点として与える有効モデルを用いて物理量の算出を行う。

演算モデルを構成するときには、物質の電子状態に対する近似解を与える既存のモデルである、局所密度近似法（ＬＤＡ）、ＬＤＡ＋Ｕ、一般化勾配近似法（ＧＧＡ）、スピン依存ＧＧＡ、メタＧＧＡ、ハートレーフォック法、配置間相互作用法、及び摂動ＭＰ（Moeller-Plesset）法などを、初期条件として用いることができる。

必要充分な揺らぎの寄与を数式４が表現する非局所演算子を用いてＧバー_{Xi, εi, gi}に取り込むことにより、計算が与えている近似解が、真の解（厳密解）の充分近傍にまで接近すると、ユニバーサルエネルギ汎関数F_λ=1 [Ψ]からＧバー_{Xi, εi, gi}を定めるエネルギ汎関数（数式１１を参照）まで、制限付き探索法によるエネルギ汎関数のパラメータ変化に関する微分が連続になり、演算モデルが与えるn_Ψ (r)が、クーロン系が与える密度と同一相に含まれるようになる。このとき、定めている数値的精度限界の範囲内で波動関数Ψ_{Xi, εi, gi}は厳密解が示す物性と同一の物性を与える。

（数式１１）

次に、演算モデルが与えている近似解が、真の解（厳密解）の充分近傍にまで接近していることの確認を行う方法を開示する。複数の初期条件と、複数の揺らぎの導入順序を組み合わせて与えることにより、厳密解に漸近する演算モデル系列を複数与えることができる。このとき、系列ごとにモデル空間内でコーシー列を発生させる演算モデルの系列を常に構成することができる。それぞれの演算モデルにおいて明確に定義される一電子軌道が量子化を発生することを通して、上位モデル系列の発生手順が定義される。そこで、全ての系列が充分な揺らぎの取り込みがなされていることが系列ごとに定められる。同一系に対する複数系列が同じ物性を発生することは、電子密度が近接していくことを数値的に確認することでいつでも評価できる。よって、数値的に発生させられる全ての演算モデルが与えている近似解が、密度に関して同じ極限を与えるコーシー列を発生しているとき、厳密解に充分近傍にまで漸近していることが確認される。

すなわち、本実施の形態では、多種多様の演算モデルを含むモデル空間内で、厳密解に到達しうる複数の経路を与えることができるため、電子相関を取込む経路を適宜設定することができる。例えば、強相関効果をフェルミ準位から順次取り込み、全スケールに依存した相関を考慮することでクーロン相互作用に近づけることができる。
そこで、Ｘ_i ，Ｘ_i+1 ，Ｘ_i+2 ，…のようにパラメータを変化させたときに表れる変分エネルギの点列に対して数値的収束を調べることにより、厳密解に漸近する方向を見出すことが可能となる。
厳密解としての基底状態における物性が再現されていることは、電子密度がＸ_i を変化させる全ての方向について変動しなくなることで定められる。このとき、変分エネルギの最小化が完了した時点の電子密度が得られている。

量子力学系の定常状態に対する決定方程式を、連続座標空間を設定して、その上の微分方程式として得る際に、共変普遍性を満たす場の理論の表現から自由度の分離を行うことによる簡略化がなされた有効方程式を得る手続きが定義可能となるためには、時間的に不変な物理量期待値が与えるＣ数化された変量を与える自由度が見出されている必要がある。その例は、電子系のシュレディンガー方程式として採用される表現では、古典質点系として扱われている原子核配置であり、同時にマクスウェル方程式の解である古典場として扱われている静電場である。この課題設定によって電子系に対して付随して電子密度という秩序変数発生と、有効ハミルトニアンの決定が起こっている。このことは、現実の問題において自発的対称性の破れを伴う秩序変数が発生している状況では常に現れる。このようなＣ数化された物理変数発生を既知とし得ない純粋な多粒子系の量子状態決定問題に対しては、自由度の分離は不可能である。現実の物性を議論する多電子系では、低温での原子核の量子流体化が発生する場合を除くと、凝縮相において常に自由度の分離が可能であることが、電子密度の測定実験結果から確定されている。

この不可能性を除いた量子力学系の定常状態に対する決定方程式を得る際にさらに別の時間的に不変な物理量期待値が与えるＣ数化されたある変量を与える自由度が見出されれば、密度だけでなくより一般化された秩序変数などＣ数化された物理変数を用いることで有効方程式を改善することが直ちに可能となる。そこで、本出願における密度汎関数変分法では、一般的に多体量子系に対して成立する方法論を与えることができる。

密度汎関数の発生は外部スカラーポテンシャル中で多電子系の運動決定が行われることによる。エネルギ密度汎関数の評価をする際、従来の局所密度近似、局所スピン密度近似、一般化勾配近似、スピン依存一般化勾配近似、を含む汎用性が高いものの物質系によらずに関数形が汎関数として与えられている従来の方法を用いるのであれば、本願で開示する手法の適用が直ちに可能となる。

本願発明者は、国際公開第２００７／１４１９４２号パンフレットにて、多配置参照密度汎関数法に基づく有効多体電子計算として揺らぎ参照決定法を与えた。この決定法により、一電子密度に加えて、局所密度揺らぎなどの正定値性を備えた変数（相関関数）を再現するエネルギ汎関数を定義し（数式１２を参照）、物理量の再現を行うことができる。

（数式１２）

揺らぎ参照決定法を用いると、数式１２の残差交換・相関エネルギ汎関数に局所密度近似したモデルを用いることで、ハートレーフォック法、配置間相互作用法、及び摂動ＭＰ法などの近似的電子状態計算が与える局所揺らぎを再現する密度汎関数理論に基づいた演算モデルを直ちに構成することができる。

したがって、ハートレーフォック法、配置間相互作用法、及び摂動ＭＰ法などの電子状態の近似解を与える演算モデルを初期条件として、演算モデルの系列を取り入れる揺らぎを定める演算子空間を拡張してゆくことにより与え、当該系列からコーシー列を見出すことで真の解の物性に到達する手続きを与えることができる。さらに、この決定プロセスにおいて、数値的精度限界を設定すると、有限の計算ステップで収束する計算プロセスを、点列｛Ｘ_i ｝の発生を通して定めることができる。

図２は本実施の形態に係る電子状態計算装置の内部構成を示すブロック図である。電子状態計算装置１０は、ＣＰＵ１１、ＲＯＭ１３、ＲＡＭ１４、入力ＩＦ１５、出力ＩＦ１６、補助記憶装置１７、及び記憶装置１８を備えており、これらの各ハードウェアがバス１２を介して互いに接続されている。ＲＯＭ１３には、ハードウェア各部を制御するための制御プログラムが格納されている。ＣＰＵ１１は、ＲＯＭ１３に格納されている制御プログラムをＲＡＭ１４上にロードして実行することにより、前述したハードウェア各部の制御を行う。入力ＩＦ１５には、マウス、キーボードなどの入力デバイス２１が接続される。一方、出力ＩＦ１６には、ＣＲＴ、液晶ディスプレイなどの出力デバイス２２が接続される。

補助記憶装置１７は、本実施の形態で説明する電子状態計算方法をコンピュータで実現するためのコンピュータプログラムを記録したＦＤ、ＣＤ−ＲＯＭなどの記録媒体Ｍから前記コンピュータプログラムを読取るためのＦＤドライブ、ＣＤ−ＲＯＭドライブなどの読取装置を備える。補助記憶装置１７によって読取られたコンピュータプログラムは記憶装置１８に記憶される。演算手段としてのＣＰＵ１１は、記憶装置１８に記憶されたコンピュータプログラムを記憶手段としてのＲＡＭ１４上にロードして実行することにより、装置全体を本発明に係る電子状態計算装置として機能させる。ＲＡＭ１４には、前記コンピュータプログムの他、入力デバイス２１を通じて入力された各種情報、ＣＰＵ１１による演算の途中結果及び最終結果等が記憶される。

なお、前記コンピュータプログラムを記録する記録媒体Ｍとしては、前述したＦＤ及びＣＤ−ＲＯＭの他に、ＭＯ、ＭＤ、ＤＶＤ−ＲＯＭ等の光ディスク、ハードディスク等の磁気記録媒体、ＩＣカード、メモリカード、光カード等のカード型記録媒体、マスクＲＯＭ、ＥＰＲＯＭ（Erasable Programmable Read Only Memory）、ＥＥＰＲＯＭ（Electrically Erasable Programmable Read Only Memory）、フラッシュＲＯＭ等による半導体メモリの利用も可能である。また、インターネットを含む通信ネットワークを接続可能なシステムを構成することにより、通信ネットワークから前述のコンピュータプログラムをダウンロードするようにしてもよい。更に、ＲＯＭ１３の内部に前述のコンピュータプログラムを予め格納しておく構成であってもよい。

記憶装置１８は、例えばハードディスクドライブであり、補助記憶装置１７によって読み取られたコンピュータプログラム、電子状態の計算に必要な初期データ、電子状態の計算により得られた途中結果、最終結果などが記憶される。また、記憶装置１８が持つ記憶領域の一部は、物質の電子状態に対する近似解を与える既存のモデルを記憶するモデル記憶領域として利用される。この既存のモデルには、前述したように、ＬＤＡ、ＬＤＡ＋Ｕ、ＧＧＡ、スピン依存ＧＧＡ、メタＧＧＡ、ハートレーフォック法、配置間相互作用法、摂動ＭＰ法などが含まれる。
記憶装置１８に記憶されているコンピュータプログラム、初期データ、演算モデル等は、電子状態を計算する際に読み出され、ＲＡＭ１３に一時的に格納される。電子状態計算装置１０のＣＰＵ１１は、ＲＡＭ１３に格納されたコンピュータプログラムを実行することにより、電子状態の計算を行う。

なお、本実施の形態では、記録媒体Ｍに記録された本発明に係るコンピュータプログラムを記憶装置１８に記憶させて利用する構成としたが、予めＲＯＭ１３に組み込まれていてもよい。また、電子状態計算装置１０が通信手段を備え、前記コンピュータプログラムを通信により取得する構成としてもよいことは勿論のことである。

以下、具体的な計算手順について説明する。図３は電子状態計算装置１０が実行する処理の手順を示すフローチャートである。電子状態計算装置１０は、まず、初期値の設定を行う（ステップＳ１１）。初期値の設定は、計算対象となる物質の原子座標Ｒ_I を設定し、有効相互作用の形式Ｘ_i と組み合わせるε_i ，ｇ_i を定めることによって行う。初期値の設定を行うことにより、モデルが１つ定まる。原子座標Ｒ_I の設定は、入力デバイス２１を通じて原子座標Ｒ_I のデータを受付け、ＲＡＭ１４に格納することによって行う。または、予め記憶装置１８に原子座標Ｒ_I のデータを記憶しておき、設定時に記憶装置１８からそのデータを読み込み、ＲＡＭ１４に格納することによって行う。
また、ε_i ，ｇ_i の関数形などはサブルーチンとして与えられているものとする。Ｘ_i は揺らぎ項の組を示す。揺らぎ項は、相互作用パラメータＸ_n ⁽ⁱ⁾とα_n ⁽ⁱ⁾とを含むように、例えば、数式６によって与えられる。

次いで、電子状態計算装置１０は、Ｘ_i ＝０（ｉ＝１，２，３，…）にセットする（ステップＳ１２）。ここで、Ｘ_i ＝０にセットするとは、数式６における相互作用パラメータＸ_n ⁽ⁱ⁾をゼロにすることを意味する。

次いで、電子状態計算装置１０のＣＰＵ１１は、設定された初期値を用いて、外部スカラーポテンシャルｖ_ext （ｒ）を算出する（ステップＳ１３）。 外部スカラーポテンシャルｖ_ext （ｒ）の算出方法は公知であり、原子座標ＲI から一通りに決定される。

次いで、電子状態計算装置１０のＣＰＵ１１は、有効系自己無撞着解を決定し、局在軌道の算出を行う（ステップＳ１４）。Ｘ_i ＝０における有効系自己無撞着解は、通常のＬＤＡなどの計算ルーチンにより決定される。一体の展開基底がこの段階で初期値として与えられる。ユニタリ変換が原子座標Ｒ_I に応じて決定され、局在軌道が算出される。

次いで、電子状態計算装置１０のＣＰＵ１１は、ループカウンタｉを１にセットした上で（ステップＳ１５）、作用素Ｖ_red ^Xi を算出する（ステップＳ１６）。このステップで、射影演算子Ｐ⁽ⁱ⁾ とユニタリ変換の形式及び操作とを一連のものとして与える。すなわち、ループによって射影演算子とユニタリ変換の形式及び操作とが更新されてもよいこと意味する。更新によって、上位モデルに至る。

次いで、電子状態計算装置１０のＣＰＵ１１は、Ｘ_i におけるエネルギ汎関数の最適化を行う（ステップＳ１７）。自己無撞着解を持つモデルの列であるモデル列において、モデル同士が数値的に誤差範囲で同じ密度を与え、かつクーロン系の変分エネルギが最小化した値を持つことを確認したとき、収束したモデルＸ_i が得られたと判定する。一体の展開基底がこの段階で再決定される。また、ユニタリ変換が原子座標Ｒ_I に応じて決定され、局在軌道が算出される。

次いで、電子状態計算装置１０のＣＰＵ１１は、電子密度ｎ_ΨXi を算出する（ステップＳ１８）。電子密度ｎ_ΨXi は、収束モデルに対してその波動関数Ψ_Xiから算出される。

次いで、電子状態計算装置１０のＣＰＵ１１は、Ｘ_i+1 におけるエネルギ汎関数の最適化を行う（ステップＳ１９）。本願では、射影演算子を変更してＰ⁽ⁱ⁺¹⁾ とし、ユニタリ変換も同時に変更することにより、アップコンバージョンを行って、上位モデルＸ_i+1 を作成する。この上位モデルの列において、収束したモデルＸ_i+1 を得る。

次いで、電子状態計算装置１０のＣＰＵ１１は、電子密度ｎ_ΨXi+1 を算出する（ステップＳ２０）。電子密度ｎ_ΨXi+1 は、収束モデルに対してその波動関数Ψ_Xi+1から算出される。

次いで、電子状態計算装置１０のＣＰＵ１１は、モデルが収束したか否かを判定する（ステップＳ２１）。すでに下位モデルＸ_i から上位モデルＸ_i+1 までのモデル系列｛Ｘ_i ｝を与えている。このモデル系列｛Ｘ_i ｝内でモデル列Ｘ_i が電子密度及び変分エネルギに関して収束しているときに、収束したと判定する。
モデル系列が収束していないと判定した場合（Ｓ２１：ＮＯ）、電子状態計算装置１０のＣＰＵ１１は、ループカウンタｉを１だけインクリメントした上で（ステップＳ２２）、処理をステップＳ１６へ戻す。

モデルが収束したと判定した場合（Ｓ２１ＹＥＳ）、電子状態計算装置１０のＣＰＵ１１は、未計算のモデルが存在するか否かを判断する（ステップＳ２３）。未計算のモデル（例えば、図１に示すＸ_j に関するモデル）が存在する場合には（Ｓ２３：ＹＥＳ）、処理をステップＳ１１へ戻し、未計算のモデルについて初期値を設定する。そして、ステップＳ１２〜ステップＳ２１の処理を新たに設定したモデルについて実行することにより、例えば、図１に示すＸ_j に関するモデルについての収束性を判定する。

未計算のモデルが存在しない場合（Ｓ２３：ＮＯ）、すでに複数のモデル系列｛Ｘ_i ｝、｛Ｘ_j ｝が得られている。電子状態計算装置１０のＣＰＵ１１は、これらのモデル系列を含むモデル列が収束したか否かを判定する（ステップＳ２４）。この異なるモデル系列が収束先の上位モデルに関して、それらの与える電子密度が収束した値を得ているときに、モデル列が収束したと判定する。
モデル列が収束していない場合（Ｓ２４：ＮＯ）、まだ充分に揺らぎの取り込みに至っていないことから、モデル系列内の上位モデルを定める射影演算子とユニタリ変換とについて更新を行って、再びモデル系列をアップコンバージョンすべく、処理をステップＳ１６へ戻す。

モデル列が収束したと判定した場合（Ｓ２４：ＹＥＳ）、得られた２種類の演算モデルを含むモデル列が共通の収束点として与える有効モデルを用いて物理量の算出を行う（ステップＳ２５）。

以上のように、本実施の形態では、複数の一体有効モデル（ＬＤＡ，ＧＧＡなど）から出発し、クーロン解に漸近するモデル系列を複数構成することができる。
密度汎関数法が与えるエネルギ汎関数法が示す状態のレベル交差点が有限数である。また、一体有効ハミルトニアンが必ず量子化を示すことから、上位モデルを形成したときに有限回数でモデル系列が収束する。したがって、上記モデル系列は、有限回数の計算でクーロン解の示す物性に至る。

また、本実施の形態では、モデル系列間の比較を電子密度に関する距離を用いて行うことができる点を利用しており、複数のモデル系列の中で最も計算効率が高い方法を有限の計算回数で決定することができる。

なお、本実施の形態では、２種類の演算モデルを用いてモデル列の収束性を判定したが、３種類以上の演算モデルを用いてモデル列の収束性を判定するようにしてもよい。

本実施の形態では、特定の演算モデル（例えば、局所密度近似）を出発点として、演算モデルを逐次更新しながら、演算を実行する構成としたが、各演算モデルによる自己無撞着解をＬＡＰＷ法（LinearizedAugmented Plane Wave Method法）を用いた並列計算、及び数値規定展開法により求める構成としてもよい。

実施の形態２．
実施の形態１で説明した電子状態計算方法に、アップコンバージョンのプロセスを導入することも可能である。
実施の形態２では、実施の形態１の電子状態計算方法においてアップコンバージョンのプロセスを導入した形態について説明を行う。

以下の数式１３は、揺らぎ効果が重要になる一電子バンドや一電子軌道に対して射影演算子を導入し、その揺らぎ項が与える有効相互作用の強さを調整するパラメータとして遮蔽相互作用型の相互作用強度を導入したモデルによる、アップコンバージョン法の実施形態２に用いられるモデルの定義式である。

（数式１３）

アップコンバージョンのプロセスは、初期モデルをＬＤＡなどに取った場合には、２段階のプロセスとして纏めることができる。まず、揺らぎを導入する重要な一電子バンドや一電子軌道が、フェルミ準位近傍から順に取り入れられていく。
ここで、遮蔽相互作用を定めるパラメータκを正の有限値に導入すると、中距離相関効果がまず取り込まれる。

実施の形態１において説明したように、モデル列の収束点を示すモデルが得られたとき、このモデルを特徴付けるブロッホ軌道またはワニエ軌道を得ているため、これらの軌道により定められる射影演算子Ｐ_A ハットを導入する。ここで、λは、上記揺らぎを取り込んだモデルから最終的なクーロン系に至る単一連続パラメータによるアップコンバージョンの第二プロセスを示しているパラメータで、０から１までの数値をとる。

具体的な実施の形態として、Ｓｒ₂ ＣｕＯ₃ の電子状態計算を例にとる。
ＧＧＡによるコーン・シャム方程式の解を確認すると、フェルミ準位近傍にＣｕ一原子当たり１つの１次元性バンドが生じていることが直ちに確認できる。このバンド上のワニエ軌道に対して射影演算子を導入する。遮蔽相互作用を取り込むことは、ワニエ軌道上の対角型相関項のみからアップコンバージョンの第一プロセスを始めることに相当する。そこで、相関パラメータとしては、遮蔽相互作用の積分が与えるＵを単一パラメータとして採用することになる。相関パラメータを導入した多配置参照密度汎関数法計算は、自己無撞着に解を与えることが出来る。このパラメータＵを変化させていくと、多配置参照密度汎関数法による多電子状態が、量子相関発生した状態ベクトルを自己無撞着解として、密度を伴って与えられる。密度汎関数変分法の変分原理に従って、ΔＥバーを評価すると、変分エネルギがＵの関数として得られる。

図４は変分エネルギの評価結果の一例を示す図である。図４において、横軸は、相関パラメータＵ（単位はエレクトロンボルト）、縦軸は変分エネルギ（単位はリュードベリ）を表す。ここで、モデルが与える変分エネルギは、相関パラメータＵの連続関数となることが分かっている。そこで、この図から明らかなように、密度汎関数変分法の原理によって、変分エネルギを最小化する最適な相関パラメータを決定することができる。

決定された相関パラメータのモデルをさらに初期条件として、アップコンバージョンの第二プロセスを実行することもできる。その場合には、λを変化させて残差相関を取り込んでいくことによって、モデル系列の変動を起こさせ、その変化を数値的に追跡することで、モデル空間での相転移発生の有無を検証することにより行うことができる。

ここで、λを変化させたときのモデル系列が与えるエネルギの評価は、摂動展開グリーン関数法、モンテカルロ法などによって行うことができる。摂動展開グリーン関数法においては、相関パラメータＵが有限のときに定まった局所相関効果を含むグリーン関数Ｇ₀ を求め、それをもとにして、長距離相互作用効果をダイソン方程式の解析を通して行うことができる。

この場合の近似的評価法としては、ＧＷ近似における入力グリーン関数にＧ₀ を用いて、遮蔽効果を取り込むことによる、ワンショット計算が、密度を変化させないことを通して、アップコンバージョン・プロセスの収束が確認される。密度変動が見えた場合には、相関項をより加えることを射影演算子Ｐ_A ハットの再定義を通して行い、再び多配置参照密度汎関数法計算を行うことで、アップコンバージョンの第一プロセスを再度行う。

モンテカルロ法においては、相関パラメータＵが有限のときに定まった局所相関効果を含む状態ベクトルを求め、それを基にした多体波動関数を拡散モンテカルロ法の節固定のための試行関数に採用することで、極短距離相関効果を取り込んだ計算を実施することができる。この場合の近似的評価法としては、節固定近似した拡散モンテカルロ法、変分モンテカルロ法等を採用することができる。

実施の形態３．
実施の形態１及び２で説明した電子状態計算方法を用いることにより、厳密解の示す物性を再現しながら、計算回数の最も少ない、最適化された演算モデルを定めることができる。
実施の形態３では、最適化された演算モデルの選定手法について説明する。

以下に示す数式１４は、クーロン系と単一連続パラメータで凸性をもって繋がる連続モデル系列を与える、モデル汎関数を示したものである。Ｖハット・チルダ（λ）は、λ微分から正定値演算子を発生する相関項として与えられたものである。λは０から１までの値をとり、λ＝１において、もとのクーロン相互作用する系が再現されることになる。

（数式１４）

λ＝０においては、相関項に対して、凸性を保証するために設定された射影演算子を用いて定義されるモデル相互作用項が含まれることになる。この項は、数式１３のλ＝０での解を元にして構築することが出来る。

数式１４のλ＝０でのモデルの最適解は、そのエネルギの評価を摂動展開グリーン関数法によって行うことができるものである。これが、射影演算子Ｐ_A の再定義を用いるアップコンバージョンの第一プロセスである。このとき、密度を変動させずにエネルギ評価が可能な場合を探索することができる。

次に、λを０から１まで変化させたときの効果を、変分計算によって評価する。この数式１４が与えるアップコンバージョンの第二プロセスにおいては、ｎハット・チルダ（ｒ）のチルダ付きの演算子が与えるモデル汎関数のλ微分が正定値であることを通して、相転移回数の有限性定理が成立する条件が与えられる。この条件は、アップコンバージョン・プロセスにおける密度の不変性である。

したがって、計算実施が可能な第二プロセスにおける密度の不変性確認によって、それ以上のアップコンバージョン・プロセスが必要でないことが確認され、全プロセスの収束判定が、有限の数値計算プロセス内で定めた誤差範囲において可能となる。

図５は最適モデル決定のプロセスを与える原理を説明する説明図である。既に、実施の形態１を通して、複数の収束モデル列が与えられているとする。このとき、λ変形されたモデルを数式１４に従って構成することができる。このλ変形プロセスにおいて、既に充分に密度に関する収束が得られている。密度を変動させない範囲で残る相関による相転移点の発生は、高々有限回数である。なお、図５では明示されていないが、相転移点はパラメータλの変動に対して現れるグラフの屈曲点として生じ、この点においてλディニ微分が飛びを示す。このとき、厳密解の示す物性を再現するモデルにおいては、内部秩序変数の変化によるモデル空間における相転移発生が生じないことを用いている。

収束した密度を再現し、充分な相関を取り込んだモデルの収束点から始めて、この相転移点が発生せずにλ＝０まで至ることができるモデルを選択することによって、厳密解の示す物性を再現しながら、計算回数の最も少ない、最適化されたモデルを選び出すことができる。

図５の原理からは、相転移発生を通した確認が出来ることになるが、これを数値計算で実現するための方法が、図６において示される、電子状態計算装置が最適モデル決定のために実行する処理の手順を示すフローチャートである。実施の形態１の図３により示したフローチャートのステップＳ１１からステップＳ２４により、収束モデル列が複数得られているとする。次に、電子状態計算装置は、収束密度を与える複数のモデル間でλ＝０での相の差異を確認することにより、密度が収束していると判定されていながら、内部秩序変数の変動により相転移が発生しているか否かを判断する（ステップＳ３１）。λ変形において相転移発生があると確認した場合、例えば、図５の２つのモデル系列が得られている場合、（Ｓ３１：ＹＥＳ）、それらのモデル列を棄却し、収束したモデル系列が同一の相を与えることを確認するまで、ステップＳ１６に戻って揺らぎを更に取り込んだモデル系列の探索を行う。ステップＳ３１で相転移発生がないと確認した場合（Ｓ３１：ＮＯ）、同一相を与える複数のモデル系列を含むモデル列が得られる。このとき、λ変形に対する相転移発生は、アルゴリズム上、摂動展開グリーン関数法によるアップコンバージョンの第二プロセスが、密度変動を生じないことに加えて、各種応答の変動が無いことによって確認されることになる。

このとき、複数のモデル系列が、λ変形の上位において同一の物性を示すことが系の真の解を得たことの確認を与えることになる。この条件が満たされたとき、λ変形プロセスの下位モデルまでにおいて、相転移点の発生が無いものを選択することで、最適モデル判定を行う（ステップＳ３２）。

最後に、最適モデルを用いた物理量の算出を行い（ステップＳ３３）、計算規模に対して最も適切に選択されたモデル系による物性評価が可能となる。

以上の実施の形態では、既知結晶構造に対する電子状態の最適化計算について説明したが、分子動力学の手法を用いて系のダイナミックスを計算することも可能である。系のダイナミックスを計算するためには、ある時間発展ステップにおける有効多体波動関数を求め、バルク秩序についてＭＤ計算（MD : Molecular Dynamics）を行う。このＭＤ計算では、まず、拡張コーン・シャム方程式により得られた電子密度、外部スカラーポテンシャル、及び局在軌道を用いて、原子間力及び応力を求め、得られた原子間力及び応力から物質構造変化率を算出し、微小時間後の物質の構造を順次決定する手法が採用される。

また、本実施の形態、または分子動力学法の形態により得られる高精度に決定された基底状態を、引き続いて時間依存問題に対する初期条件として用いることができる。この方法によって、電子密度と電流密度を基本となる秩序変数として用いた時間依存電流密度汎関数法を実現することができる。

１０ 電子状態計算装置
１１ ＣＰＵ
１２ バス
１３ ＲＯＭ
１４ ＲＡＭ
１５ 入力ＩＦ
１６ 出力ＩＦ
１７ 補助記憶装置
１８ 記憶装置
２１ 入力デバイス
２２ 出力デバイス

Claims (11)

1. 計算装置を用いて物質の電子状態を計算する方法において、
   前記計算装置は、
   物質の電子状態に対する近似解を与える複数の演算モデルの夫々を構成要素として含む集合を設定し、
   複数の電子軌道上に存在する少なくとも２つ以上の電子を含む電子系に働く有効相互作用を含む有効ハミルトニアンを、各演算モデルを特定するために定め、
   前記有効ハミルトニアンの自己無撞着解を前記集合内の各演算モデルを用いて計算する過程で、計算した自己無撞着解同士が連続的に変化する方向を特定しながら、前記集合がなす空間内で距離が近接する複数の演算モデルのうちで最適な演算モデルを量子力学的変分法に基づいて定め、
   前記最適な演算モデルを逐次更新する際に、前記有効ハミルトニアンの自己無撞着解による電子系の変分エネルギを評価し、
   評価した変分エネルギが計算すべき厳密解のエネルギに近接してゆき、しかも変分エネルギが単調減少凸関数をなすように演算モデルを更新して、一又は複数の変分エネルギの系列から前記電子状態の厳密解を計算する
   ことを特徴とする電子状態計算方法。
2. 局所密度近似法、一般化勾配近似法、又はハートレーフォック法を基に、有効相互作用として揺らぎを示す非局所型演算子を定める手続きを通じて前記有効ハミルトニアンを決定することを特徴とする請求項１に記載の電子状態計算方法。
3. 前記揺らぎは、クーロン相互作用とハートレー平均場項のずれである密度揺らぎを含むことを特徴とする請求項２に記載の電子状態計算方法。
4. 前記厳密解に近接してゆく変分エネルギの系列を与える一又は複数の演算モデルの系列のうち、電子状態を示す相が転移せずに厳密解と同一視される相を与え、かつ最小の計算ステップで前記有効ハミルトニアンの自己無撞着解に至る演算モデルを含むものを定めることを特徴とする請求項１に記載の電子状態計算方法。
5. 各演算モデルによる自己無撞着解をＬＡＰＷ法（Linearized Augmented Plane Wave 法）、ＰＡＷ法（Projector Augmented Wave 法）、又は数値基底展開法を用いた並列計算により求めることを特徴とする請求項１から請求項４の何れか１つに記載の電子状態計算方法。
6. 物質の電子状態を計算する装置において、
   物質の電子状態に対する近似解を与える複数の演算モデルの夫々を構成要素として含む集合を設定する手段と、
   複数の電子軌道上に存在する少なくとも２つ以上の電子を含む電子系に働く有効相互作用を含む有効ハミルトニアンを、各演算モデルを特定するために定める手段と、
   前記有効ハミルトニアンの自己無撞着解を前記集合内の各演算モデルを用いて計算する過程で、計算した自己無撞着解同士が連続的に変化する方向を特定しながら、前記集合がなす空間内で距離が近接する複数の演算モデルのうちで最適な演算モデルを量子力学的変分法に基づいて定める手段と、
   前記最適な演算モデルを逐次更新する際に、前記有効ハミルトニアンの自己無撞着解による電子系の変分エネルギを評価する手段と、
   評価した変分エネルギが計算すべき厳密解のエネルギに近接してゆき、しかも変分エネルギが単調減少凸関数をなすように演算モデルを更新する手段と、
   一又は複数の変分エネルギの系列から前記電子状態の厳密解を計算する手段と
   を備えることを特徴とする電子状態計算装置。
7. 局所密度近似法、一般化勾配近似法、又はハートレーフォック法を基に、有効相互作用として揺らぎを示す非局所型演算子を定める手続きを通じて前記有効ハミルトニアンを決定するようにしてあることを特徴とする請求項６に記載の電子状態計算装置。
8. 前記揺らぎは、クーロン相互作用とハートレー平均場項のずれである密度揺らぎを含むことを特徴とする請求項７に記載の電子状態計算装置。
9. 前記厳密解に近接してゆく変分エネルギの系列を与える一又は複数の演算モデルの系列のうち、電子状態を示す相が転移せずに厳密解と同一視される相を与え、かつ最小の計算ステップで前記有効ハミルトニアンの自己無撞着解に至る演算モデルを含むものを定めることを特徴とする請求項６に記載の電子状態計算装置。
10. 各演算モデルによる自己無撞着解をＬＡＰＷ法（Linearized Augmented Plane Wave 法）、ＰＡＷ法（Projector Augmented Wave 法）、又は数値基底展開法を用いた並列計算により求めるようにしてあることを特徴とする請求項６から請求項９の何れか１つに記載の電子状態計算装置。
11. コンピュータに、物質の電子状態を計算させるコンピュータプログラムにおいて、
    物質の電子状態に対する近似解を与える複数の演算モデルの夫々を構成要素として含む集合を設定させるステップと、
    複数の電子軌道上に存在する少なくとも２つ以上の電子を含む電子系に働く有効相互作用を含む有効ハミルトニアンを、各演算モデルを特定するために定めさせるステップと、
    前記有効ハミルトニアンの自己無撞着解を前記集合内の各演算モデルを用いて計算させる過程で、計算させた自己無撞着解同士が連続的に変化する方向を特定しながら、前記集合がなす空間内で距離が近接する複数の演算モデルのうちで最適な演算モデルを量子力学的変分法に基づいて定めさせるステップと、
    前記最適な演算モデルを逐次更新させる際に、前記有効ハミルトニアンの自己無撞着解による電子系の変分エネルギを評価させるステップと、
    評価させた変分エネルギが計算すべき厳密解のエネルギに近接してゆき、しかも変分エネルギが単調減少凸関数をなすように演算モデルを更新させ、一又は複数の変分エネルギの系列から前記電子状態の厳密解を計算させるステップと
    をコンピュータに実行させることを特徴とするコンピュータプログラム。

Images