(57)【要約】  電子状態計算装置は、結晶構造、原子番号の入力を受付け（Ｓ１１）、通常のコーン・シャム理論の電子状態計算から揺らぎを発生し得る原子集団を特定することにより、参照系を算出する（Ｓ１２）。次に、拡張コーン・シャム理論に基づき、有効多体系について自己無撞着計算を行う（Ｓ１３）。参照系について得られた密度揺らぎと、自己無撞着計算により得られた密度揺らぎとが一致するか否かを判断し（Ｓ１４）、一致する場合の交換相関エネルギ及び局所相互作用のパラメータを取得する。取得したパラメータにより有効ハミルトニアンを決定し、既知結晶構造に対する電子状態の最適化を行う（Ｓ１６）。
2007141942.tif 000002
【特許請求の範囲】
【請求項１】
第一原理計算により物質の電子状態を計算する方法において、
前記物質の物性を特徴付ける物理量の揺らぎを再現するモデルを決定し、決定したモデルを初期条件とした第一原理計算により前記物質の電子状態を求めることを特徴とする電子状態計算方法。
【請求項２】
次式で表される残差交換相関エネルギ汎関数を含んだ表式により定式化される拡張コーン・シャム理論に基づく最適化手法を用いて前記電子状態を求めることを特徴とする請求項１に記載の電子状態計算方法。
【数１】
2007141942.tif 000024
ここで、
ｎΨは一電子密度、
Ｅxc［ｎΨ］は交換相関エネルギ汎関数、
Φは相関が発生していない状態の多粒子波動関数、
Ψ’は相関が発生した状態の多粒子波動関数、
Ｔは運動エネルギ演算子、
Ｖred は物理量の組｛Ｘi ｝により与えられる次の相関関数により定義される演算子である。
【数２】
2007141942.tif 000025
ただし、
ｉ＝１，２，３，…，Ｎ（Ｎは計算すべき物理量の個数）、
Ｗi は重みを表すパラメータ、
小括弧内は物理量Ｘi と物理量Ｘi の平均値との差を表す。
【請求項３】
電子間の有効相互作用を可変として前記揺らぎの大きさを変化させることにより、同じ電子密度を再現する揺らぎがある状態及び揺らぎがない状態の２つの状態における運動エネルギ及び全エネルギの評価を行うことを特徴とする請求項２に記載の電子状態計算方法。
【請求項４】
前記モデルを決定する際、多体電子系の計算手法を用いて前記物理量に対する揺らぎを計算し、電子間の交換相関ポテンシャル及びクーロン相互作用を含む拡張コーン・シャム理論に基づいて前記物理量に対する揺らぎを計算し、得られた２つの揺らぎが略等しくなるように前記モデルを記述するパラメータを決定することを特徴とする請求項１に記載の電子状態計算方法。
【請求項５】
前記拡張コーン・シャム理論は、次式で表される波動関数汎関数の最適化問題として定式化されることを特徴とする請求項４に記載の電子状態計算方法。
【数３】
2007141942.tif 000026
ここで、
Ψは多粒子波動関数、
Ｔは運動エネルギ演算子、
ｎΨは一電子密度、
Ｆ［ｎ］はユニバーサルエネルギ汎関数、
Ψ’は一電子密度ｎΨを再現する条件付きで＜Ｔ＋Ｖred ＞を最小化する多粒子波動関数、
ｖext は外部相互作用ポテンシャル、
Ｖred は数式２により定義される演算子である。
【請求項６】
前記パラメータは、電子間の交換相関エネルギ及び局所相互作用に係るパラメータを含むことを特徴とする請求項４に記載の電子状態計算方法。
【請求項７】
前記計算手法は、量子モンテカルロ法、トランスコリレイティッド法、配置間相互作用法、摂動計算・グリーン関数法、又は有効ポテンシャル法であることを特徴とする請求項４に記載の電子状態計算方法。
【請求項８】
前記拡張コーン・シャム理論に基づく相関関数の計算は、平面波基底展開法、補強平面波基底展開法、実空間、数値基底展開法、混合基底展開法、局在基底展開法、又は線型化ＭＴ軌道法を用いることを特徴とする請求項２又は請求項５に記載の電子状態計算方法。
【請求項９】
前記揺らぎは、正定値かつ有界な値域を持つ２次の相関関数により記述されることを特徴とする請求項１に記載の電子状態計算方法。
【請求項１０】
前記物理量は、局所軌道上の電子に関する数演算子又は生成・消滅演算子により記述される量であることを特徴とする請求項１に記載の電子状態計算方法。
【請求項１１】
前記物理量は、前記物質の構造因子、光学応答係数、又は誘電率、透磁率、導電率、若しくは弾性率を決定する物理量であることを特徴とする請求項１に記載の電子状態計算方法。
【請求項１２】
第一原理計算により物質の電子状態を計算する装置において、
前記物質の物性を特徴付ける物理量の揺らぎを再現するモデルを決定する手段と、決定したモデルを初期条件とした第一原理計算により前記物質の電子状態を求める手段とを備えることを特徴とする電子状態計算装置。
【請求項１３】
次式で表される残差交換相関エネルギ汎関数を含んだ表式により定式化される拡張コーン・シャム理論に基づく最適化手法を用いて前記電子状態を求めるようにしてあることを特徴とする請求項１２に記載の電子状態計算装置。
【数４】
2007141942.tif 000027
ここで、
ｎΨは一電子密度、
Ｅxc［ｎΨ］は交換相関エネルギ汎関数、
Φは相関が発生していない状態の多粒子波動関数、
Ψ’は相関が発生した状態の多粒子波動関数、
Ｔは運動エネルギ演算子、
Ｖred は物理量の組｛Ｘi ｝により与えられる次の相関関数により定義される演算子である。
【数５】
2007141942.tif 000028
ただし、
ｉ＝１，２，３，…，Ｎ（Ｎは計算すべき物理量の個数）、
Ｗi は重みを表すパラメータ、
小括弧内は物理量Ｘi と物理量Ｘi の平均値との差を表す。
【請求項１４】
電子間の有効相互作用を可変として前記揺らぎの大きさを変化させることにより、同じ電子密度を再現する揺らぎがある状態及び揺らぎがない状態の２つの状態における運動エネルギ及び全エネルギの評価する手段を備えることを特徴とする請求項１３に記載の電子状態計算装置。
【請求項１５】
多体電子系の計算手法を用いて前記物理量に対する揺らぎを計算する手段と、電子間の交換相関ポテンシャル及びクーロン相互作用を含む拡張コーン・シャム理論に基づいて前記物理量に対する揺らぎを計算する手段と、得られた２つの揺らぎが略等しくなるように前記モデルを記述するパラメータを決定する手段とを備えることを特徴とする請求項１２に記載の電子状態計算装置。
【請求項１６】
次式で表される波動関数汎関数の最適化問題として前記拡張コーン・シャム理論を定式化してあることを特徴とする請求項１５に記載の電子状態計算装置。
【数６】
2007141942.tif 000029
ここで、
Ψは多粒子波動関数、
Ｔは運動エネルギ演算子、
ｎΨは一電子密度、
Ｆ［ｎ］はユニバーサルエネルギ汎関数、
Ψ’は一電子密度ｎΨを再現する条件付きで＜Ｔ＋Ｖred ＞を最小化する多粒子波動関数、
ｖext は外部相互作用ポテンシャル、
Ｖred は数式５により定義される演算子である。
【請求項１７】
コンピュータに、物質の電子状態を第一原理計算により計算させるコンピュータプログラムにおいて、
コンピュータに、前記物質の物性を特徴付ける物理量の揺らぎを再現するモデルを決定させる第１ステップと、コンピュータに、決定させたモデルを初期条件とした第一原理計算により前記物質の電子状態を計算させる第２ステップとを有することを特徴とするコンピュータプログラム。
【請求項１８】
前記第一原理計算として、次式で表される残差交換相関エネルギ汎関数を含んだ表式により定式化される拡張コーン・シャム理論に基づく最適化手法を用いることを特徴とする請求項１７に記載のコンピュータプログラム。
【数７】
2007141942.tif 000030
ここで、
ｎΨは一電子密度、
Ｅxc［ｎΨ］は交換相関エネルギ汎関数、
Φは相関が発生していない状態の多粒子波動関数、
Ψ’は相関が発生した状態の多粒子波動関数、
Ｔは運動エネルギ演算子、
Ｖred は物理量の組｛Ｘi ｝により与えられる次の相関関数により定義される演算子である。
【数８】
2007141942.tif 000031
ただし、
ｉ＝１，２，３，…，Ｎ（Ｎは計算すべき物理量の個数）、
Ｗi は重みを表すパラメータ、
小括弧内は物理量Ｘi と物理量Ｘi の平均値との差を表す。
【請求項１９】
コンピュータに、電子間の有効相互作用を可変として前記揺らぎの大きさを変化させるステップと、コンピュータに、同じ電子密度を再現する揺らぎがある状態及び揺らぎがない状態の２つの状態における運動エネルギ及び全エネルギの評価を行わせるステップとを有することを特徴とする請求項１８に記載のコンピュータプログラム。
【請求項２０】
前記第１ステップは、多体電子系の計算手法を用いて前記物理量に対する揺らぎをコンピュータに計算させるステップと、電子間の交換相関ポテンシャル及びクーロン相互作用を含む拡張コーン・シャム理論に基づいて前記物理量に対する揺らぎをコンピュータに計算させるステップと、得られた２つの揺らぎが略等しくなるように前記モデルを記述するパラメータを決定させるステップとを含むことを特徴とする請求項１７に記載のコンピュータプログラム。
【請求項２１】
前記第一原理計算として、次式で表される波動関数汎関数を含んだ表式により定式化される拡張コーン・シャム理論に基づく最適化手法を用いることを特徴とする請求項２０に記載のコンピュータプログラム。
【数９】
2007141942.tif 000032
ここで、
Ψは多粒子波動関数、
Ｔは運動エネルギ演算子、
ｎΨは一電子密度、
Ｆ［ｎ］はユニバーサルエネルギ汎関数、
Ψ’は一電子密度ｎΨを再現する条件付きで＜Ｔ＋Ｖred ＞を最小化する多粒子波動関数、
ｖext は外部相互作用ポテンシャル、
Ｖred は数式８により定義される演算子である。
【請求項２２】
コンピュータに、物質の電子状態を第一原理計算により計算させるコンピュータプログラムが記録されているコンピュータでの読取りが可能な記録媒体において、
コンピュータに、前記物質の物性を特徴付ける物理量に関して、該物理量の揺らぎを再現するモデルを決定させるステップと、コンピュータに、決定させたモデルを初期条件とした第一原理計算により前記物質の電子状態を計算させるステップとを有するコンピュータプログラムが記録されていることを特徴とするコンピュータでの読取りが可能な記録媒体。
【請求項２３】
前記第一原理計算として、次式で表される残差交換相関エネルギ汎関数を含んだ表式により定式化される拡張コーン・シャム理論に基づく最適化手法を用いることを特徴とする請求項２２に記載の記録媒体。
【数１０】
2007141942.tif 000033
ここで、
ｎΨは一電子密度、
Ｅxc［ｎΨ］は交換相関エネルギ汎関数、
Φは相関が発生していない状態の多粒子波動関数、
Ψ’は相関が発生した状態の多粒子波動関数、
Ｔは運動エネルギ演算子、
Ｖred は物理量の組｛Ｘi ｝により与えられる次の相関関数により定義される演算子である。
【数１１】
2007141942.tif 000034
ただし、
ｉ＝１，２，３，…，Ｎ（Ｎは計算すべき物理量の個数）、
Ｗi は重みを表すパラメータ、
小括弧内は物理量Ｘi と物理量Ｘi の平均値との差を表す。
【請求項２４】
前記第一原理計算として、次式で表される波動関数汎関数を含んだ表式により定式化される拡張コーン・シャム理論に基づく最適化手法を用いることを特徴とする請求項２２に記載の記録媒体。
【数１２】
2007141942.tif 000035
ここで、
Ψは多粒子波動関数、
Ｔは運動エネルギ演算子、
ｎΨは一電子密度、
Ｆ［ｎ］はユニバーサルエネルギ汎関数、
Ψ’は一電子密度ｎΨを再現する条件付きで＜Ｔ＋Ｖred ＞を最小化する多粒子波動関数、
ｖext は外部相互作用ポテンシャル、
Ｖred は数式１１により定義される演算子である。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、第一原理計算により物質の電子状態を計算する電子状態計算方法、電子状態計算装置、コンピュータプログラム、及び該コンピュータプログラムが記録されている記録媒体に関する。
【背景技術】
【０００２】
従来、物質の物理的又は化学的性質（以下、物性と呼ぶ）を量子力学の基礎法則に則って予測する第一原理計算理論と呼ばれる方法のうち、全エネルギ及び一電子密度の再現を厳密に可能にすることで、弾性を含む力学特性、超伝導特性を含む伝導特性、誘電性、磁性などの物性を比較的高精度に再現し、しかも計算規模が実施可能な範囲である、密度汎関数理論に基づく第一原理電子状態計算理論が知られている。この計算理論は、既に物質設計に応用され、その予測精度が実験で実証された例が多く存在するものである（例えば、非特許文献１参照）。
【非特許文献１】「計算機マテリアルデザイン入門」、笠井秀明、赤井久純、吉田博編、２００５年、大阪大学出版会
【発明の開示】
【発明が解決しようとする課題】
【０００３】
従来では、通常コーン・シャム方程式と呼ばれる独立粒子系の補助方程式を導入して計算を実施するが、全エネルギ及び一電子密度以外の物理量に対して再現性は保証されておらず、これまで以下の問題点が指摘されてきた。【０００４】
まず、コーン・シャム方程式自身の与える励起エネルギはもとの物質系の電子・正孔励起を必ずしも与えず、本来電荷ギャップ形成がある物質においてギャップが存在しないと結論する場合があるという問題点を有している。【０００５】
また、電子密度とスピン密度（又は、４元電流密度）とを秩序変数とする拡張を用いた場合であっても、独立フェルミオン系の解を参照系として得ることができるのみであり、磁性及び強相関効果による超伝導、軌道秩序形成を議論する際に表現自体の不完全性から、物性の再現性に限界があるとされてきた。なお、コーン・シャム方程式の実用的応用には、さらに交換相関エネルギ汎関数の近似形が必要とされ、通常は局所密度近似（ＬＤＡ : Local Density Approximation）、密度勾配補正近似（ＧＧＡ: Generarized Gradient Approximation）が採用されてきた。【０００６】
本発明は斯かる事情に鑑みてなされたものであり、物質の物性を特徴付ける物理量の揺らぎを再現するモデルを決定し、決定したモデルを初期条件とした第一原理計算により電子状態を求める構成とすることにより、全エネルギ及び一電子密度以外の物理量、例えば、誘電応答、透磁率、系の反応過程、ダイナミクスなどを求めることができる電子状態計算方法、電子状態計算装置、コンピュータプログラム、及び記録媒体を提供することを目的とする。
【課題を解決するための手段】
【０００７】
第１発明に係る電子状態計算方法は、第一原理計算により物質の電子状態を計算する方法において、前記物質の物性を特徴付ける物理量の揺らぎを再現するモデルを決定し、決定したモデルを初期条件とした第一原理計算により前記物質の電子状態を求めることを特徴とする。【０００８】
第１発明にあっては、物質の物性を特徴付ける物理量の揺らぎを再現するモデルを決定し、決定したモデルを初期条件とした第一原理計算により電子状態を求めるようにしている。したがって、縮退のない基底状態を示す電子系のシュレーディンガー方程式と等価なエネルギ及び一電子密度を与えることができるモデルであって、物理量の揺らぎを同時再現するモデルを、ユニバーサルエネルギ汎関数を用いて波動関数に関する汎関数として与えることができた場合、モデルに内在するパラメータが一意に定まる。【０００９】
第２発明に係る電子状態計算方法は、次式で表される残差交換相関エネルギ汎関数を含んだ表式により定式化される拡張コーン・シャム理論に基づく最適化手法を用いて前記電子状態を求めることを特徴とする。【００１０】
【数１】
2007141942.tif 000003
【００１１】
ここで、ｎΨは一電子密度、Ｅxc［ｎΨ］は交換相関エネルギ汎関数、Φは相関が発生していない状態の多粒子波動関数、Ψ’は相関が発生した状態の多粒子波動関数、Ｔは運動エネルギ演算子、Ｖred は物理量の組｛Ｘi ｝により与えられる次の相関関数により定義される演算子である。【００１２】
【数２】
2007141942.tif 000004
【００１３】
ただし、ｉ＝１，２，３，…，Ｎ（Ｎは計算すべき物理量の個数）、Ｗi は重みを表すパラメータ、小括弧内は物理量Ｘi と物理量Ｘi の平均値との差を表す。【００１４】
第２発明にあっては、拡張コーン・シャム方程式において、残差交換相関エネルギを厳密な表式からの恒等変形により、通常の交換相関エネルギとそこからの差分として表記する。この差分の部分を、同じ電子密度を再現する揺らぎのある状態、無い状態の２つの状態における運動エネルギの差及びＶred の期待値（揺らぎ）として与え、この運動エネルギの変化を拡張コーン・シャム理論の枠内で与えることにより、量子モンテカルロ法、配置間相互作用法などの高精度計算手法を用いることなく電子状態の計算が可能となる。【００１５】
第３発明に係る電子状態計算方法は、電子間の有効相互作用を可変として前記揺らぎの大きさを変化させることにより、同じ電子密度を再現する揺らぎがある状態及び揺らぎがない状態の２つの状態における運動エネルギ及び全エネルギの評価を行うことを特徴とする。【００１６】
第３発明にあっては、同じ電子密度を再現する揺らぎのある状態、無い状態の２つの状態における運動エネルギ及び全エネルギを評価することによって、拡張コーン・シャム理論の枠内で電子状態の計算が可能となる。【００１７】
第４発明に係る電子状態計算方法は、前記モデルを決定する際、多体電子系の計算手法を用いて前記物理量に対する揺らぎを計算し、電子間の交換相関ポテンシャル及びクーロン相互作用を含む拡張コーン・シャム理論に基づいて前記物理量に対する揺らぎを計算し、得られた２つの揺らぎが略等しくなるように前記モデルを記述するパラメータを決定することを特徴とする。【００１８】
第４発明にあっては、参照系として計算した物理量の揺らぎを再現するようにモデルを決定するため、このモデルを決定する過程において定まるパラメータを用いて多体系の波動関数が求められる。【００１９】
第５発明に係る電子状態計算方法は、前記拡張コーン・シャム理論は、次式で表される波動関数汎関数の最適化問題として定式化されることを特徴とする。【００２０】
【数３】
2007141942.tif 000005
【００２１】
ここで、Ψは多粒子波動関数、Ｔは運動エネルギ演算子、ｎΨは一電子密度、Ｆ［ｎ］はユニバーサルエネルギ汎関数、Ψ’は一電子密度ｎΨを再現する条件付で＜Ｔ＋Ｖred ＞を最小化する多粒子波動関数、ｖext は外部相互作用ポテンシャル、Ｖred は数式２により定義される演算子である。【００２２】
第５発明にあっては、数式３により表される波動関数汎関数の最適化問題によって、クーロン系の基底状態エネルギ及びその一粒子密度が与えられる。【００２３】
第６発明に係る電子状態計算方法は、前記パラメータが、電子間の交換相関エネルギ及び局所相互作用に係るパラメータを含むことを特徴とする。【００２４】
第６発明にあっては、交換エネルギ及び相関エネルギを含めたエネルギ汎関数のルジャンドル変換を密度変分して与えられる有効ポテンシャルの定義が可能になり、ハミルトニアンの運動エネルギ部分とこの有効ポテンシャルとからなる有効一体ハミルトニアン中で、相互作用する仮想的電子系が運動するという表現が得られる。【００２５】
第７発明に係る電子状態計算方法は、前記計算手法が、量子モンテカルロ法、トランスコリレイティッド法、配置間相互作用法、摂動計算・グリーン関数法、又は有効ポテンシャル法であることを特徴とする。【００２６】
第７発明にあっては、参照系としての多体電子系の密度・密度相関を求める際に量子モンテカルロ法、トランスコリレイティッド法、配置間相互作用法、摂動計算・グリーン関数法、又は有効ポテンシャル法が用いられる。【００２７】
第８発明に係る電子状態計算方法は、前記拡張コーン・シャム理論に基づく相関関数の計算が、平面波基底展開法、補強平面波基底展開法、実空間、数値基底展開法、混合基底展開法、局在基底展開法、又は線型化ＭＴ軌道法を用いることを特徴とする。【００２８】
第８発明にあっては、拡張コーン・シャム理論に基づいて相関関数を求める際に平面波基底展開法、補強平面波基底展開法、実空間、数値基底展開法、混合基底展開法、局在基底展開法、又は線型化ＭＴ軌道法が用いられる。【００２９】
第９発明に係る電子状態計算方法は、前記揺らぎが、正定値かつ有界な値域を持つ２次の相関関数により記述されることを特徴とする。【００３０】
第９発明にあっては、物理量の揺らぎを正定値かつ有界な値域を持つ２次の相関関数により記述するため、縮体がない基底状態を示す電子系のシュレーディンガー方程式と等価なエネルギ及び一電子密度を与えることができるモデルであって、しかも正定値有界な相関関数を同時再現するモデルを、ユニバーサルエネルギ汎関数を用いて与えることができた場合、内在するパラメータが一意に定まる。【００３１】
第１０発明に係る電子状態計算方法は、前記物理量が、局所軌道上の電子に関する数演算子又は生成・消滅演算子により記述される量であることを特徴とする。【００３２】
第１０発明にあっては、物理量を局所軌道上の電子に関する数演算子又は生成・消滅演算子として記述しているため、対応する拡張コーン・シャム模型は、第一原理的不純物アンダーソンモデル及び第一原理的ハバードモデルを含む強相関電子モデルとして得られる。【００３３】
第１１発明に係る電子状態計算方法は、前記物理量が、前記物質の構造因子、光学応答係数、又は誘電率、透磁率、導電率、若しくは弾性率を決定する物理量であることを特徴とする。【００３４】
第１１発明にあっては、物理量の揺らぎを再現するモデルを決定することにより、物質の構造因子（すなわち、原子レベルでの原子配置）が定まるため、他の物理量が全て無矛盾かつ前提条件なしで定まる。【００３５】
第１２発明に係る電子状態計算装置は、第一原理計算により物質の電子状態を計算する装置において、前記物質の物性を特徴付ける物理量の揺らぎを再現するモデルを決定する手段と、決定したモデルを初期条件とした第一原理計算により前記物質の電子状態を求める手段とを備えることを特徴とする。【００３６】
第１２発明にあっては、物質の物性を特徴付ける物理量の揺らぎを再現するモデルを決定し、決定したモデルを初期条件とした第一原理計算により電子状態を求めるようにしている。したがって、縮退のない基底状態を示す電子系のシュレーディンガー方程式と等価なエネルギ及び一電子密度を与えることができるモデルであって、物理量の揺らぎを同時再現するモデルを、ユニバーサルエネルギ汎関数を用いて波動関数に関する汎関数として与えることができた場合、モデルに内在するパラメータが一意に定まる。【００３７】
第１３発明に係る電子状態計算装置は、次式で表される残差交換相関エネルギ汎関数を含んだ表式により定式化される拡張コーン・シャム理論に基づく最適化手法を用いて前記電子状態を求めるようにしてあることを特徴とする。【００３８】
【数４】
2007141942.tif 000006
【００３９】
ここで、ｎΨは一電子密度、Ｅxc［ｎΨ］は交換相関エネルギ汎関数、Φは相関が発生していない状態の多粒子波動関数、Ψ’は相関が発生した状態の多粒子波動関数、Ｔは運動エネルギ演算子、Ｖred は物理量の組｛Ｘi ｝により与えられる次の相関関数により定義される演算子である。【００４０】
【数５】
2007141942.tif 000007
ただし、ｉ＝１，２，３，…，Ｎ（Ｎは計算すべき物理量の個数）、Ｗi は重みを表すパラメータ、小括弧内は物理量Ｘi と物理量Ｘi の平均値との差を表す。【００４１】
第１３発明にあっては、拡張コーン・シャム方程式において、残差交換相関エネルギを厳密な表式からの恒等変形により、通常の交換相関エネルギとそこからの差分として表記する。この差分の部分を、同じ電子密度を再現する揺らぎのある状態、無い状態の２つの状態における運動エネルギの差及びＶred の期待値として与え、この運動エネルギの変化を拡張コーン・シャム理論の枠内で与えることにより、量子モンテカルロ法、配置間相互作用法などの高精度計算手法を用いることなく電子状態の計算が可能となる。【００４２】
第１４発明に係る電子状態計算装置は、電子間の有効相互作用を可変として前記揺らぎの大きさを変化させることにより、同じ電子密度を再現する揺らぎがある状態及び揺らぎがない状態の２つの状態における運動エネルギ及び全エネルギを評価する手段を備えることを特徴とする。【００４３】
第１４発明にあっては、同じ電子密度を再現する揺らぎのある状態、無い状態の２つの状態における運動エネルギ及び全エネルギを評価することによって、拡張コーン・シャム理論の枠内で電子状態の計算が可能となる。【００４４】
第１５発明に係る電子状態計算装置は、多体電子系の計算手法を用いて前記物理量に対する揺らぎを計算する手段と、電子間の交換相関ポテンシャル及びクーロン相互作用を含む拡張コーン・シャム理論に基づいて前記物理量に対する揺らぎを計算する手段と、得られた２つの揺らぎが略等しくなるように前記モデルを記述するパラメータを決定する手段とを備えることを特徴とする。【００４５】
第１５発明にあっては、参照系として計算した物理量の揺らぎを再現するようにモデルを決定するため、このモデルを決定する過程において定まるパラメータを用いて多体系の波動関数が求められる。【００４６】
第１６発明に係る電子状態計算装置は、次式で表される波動関数汎関数の最適化問題として前記拡張コーン・シャム理論を定式化してあることを特徴とする。【００４７】
【数６】
2007141942.tif 000008
【００４８】
ここで、Ψは多粒子波動関数、Ｔは運動エネルギ演算子、ｎΨは一電子密度、Ｆ［ｎ］はユニバーサルエネルギ汎関数、Ψ’は一電子密度ｎΨを再現する条件付で＜Ｔ＋Ｖred ＞を最小化する多粒子波動関数、ｖext は外部相互作用ポテンシャル、Ｖred は数式５により定義される演算子である。【００４９】
第１６発明にあっては、数式６により表される波動関数汎関数の最適化問題によって、クーロン系の基底状態エネルギ及びその一粒子密度が与えられる。【００５０】
第１７発明に係るコンピュータプログラムは、コンピュータに、物質の電子状態を第一原理計算により計算させるコンピュータプログラムにおいて、コンピュータに、前記物質の物性を特徴付ける物理量の揺らぎを再現するモデルを決定させる第１ステップと、コンピュータに、決定させたモデルを初期条件とした第一原理計算により前記物質の電子状態を計算させる第２ステップとを有することを特徴とする。【００５１】
第１７発明にあっては、物質の物性を特徴付ける物理量の揺らぎを再現するモデルを決定し、決定したモデルを初期条件とした第一原理計算により電子状態を求めるようにしている。したがって、縮退のない基底状態を示す電子系のシュレーディンガー方程式と等価なエネルギ及び一電子密度を与えることができるモデルであって、物理量の揺らぎを同時再現するモデルを、ユニバーサルエネルギ汎関数を用いて波動関数に関する汎関数として与えることができた場合、モデルに内在するパラメータが一意に定まる。【００５２】
第１８発明に係るコンピュータプログラムは、前記第一原理計算として、次式で表される残差交換相関エネルギ汎関数を含んだ表式により定式化される拡張コーン・シャム理論に基づく最適化手法を用いることを特徴とする。【００５３】
【数７】
2007141942.tif 000009
【００５４】
ここで、ｎΨは一電子密度、Ｅxc［ｎΨ］は交換相関エネルギ汎関数、Φは相関が発生していない状態の多粒子波動関数、Ψ’は相関が発生した状態の多粒子波動関数、Ｔは運動エネルギ演算子、Ｖred は物理量の組｛Ｘi ｝により与えられる次の相関関数により定義される演算子である。【００５５】
【数８】
2007141942.tif 000010
ただし、ｉ＝１，２，３，…，Ｎ（Ｎは計算すべき物理量の個数）、Ｗi は重みを表すパラメータ、小括弧内は物理量Ｘi と物理量Ｘi の平均値との差を表す。【００５６】
第１８発明にあっては、拡張コーン・シャム方程式において、残差交換相関エネルギを厳密な表式からの恒等変形により、通常の交換相関エネルギとそこからの差分として表記する。この差分の部分を、同じ電子密度を再現する揺らぎのある状態、無い状態の２つの状態における運動エネルギの差及びＶred の期待値として与え、この運動エネルギの変化を拡張コーン・シャム理論の枠内で与えることにより、量子モンテカルロ法、配置間相互作用法などの高精度計算手法を用いることなく電子状態の計算が可能となる。【００５７】
第１９発明に係るコンピュータプログラムは、コンピュータに、電子間の有効相互作用を可変として前記揺らぎの大きさを変化させるステップと、コンピュータに、同じ電子密度を再現する揺らぎがある状態及び揺らぎがない状態の２つの状態における運動エネルギ及び全エネルギの評価を行わせるステップとを有することを特徴とする。【００５８】
第１９発明にあっては、同じ電子密度を再現する揺らぎのある状態、無い状態の２つの状態における運動エネルギ及び全エネルギを評価することによって、拡張コーン・シャム理論の枠内で電子状態の計算が可能となる。【００５９】
第２０発明に係るコンピュータプログラムは、前記第１ステップが、多体電子系の計算手法を用いて前記物理量に対する揺らぎをコンピュータに計算させるステップと、電子間の交換相関ポテンシャル及びクーロン相互作用を含む拡張コーン・シャム理論に基づいて前記物理量に対する揺らぎをコンピュータに計算させるステップと、得られた２つの揺らぎが略等しくなるように前記モデルを記述するパラメータを決定させるステップとを含むことを特徴とする。【００６０】
第２０発明にあっては、参照系として計算した物理量の揺らぎを再現するようにモデルを決定するため、このモデルを決定する過程において定まるパラメータを用いて多体系の波動関数が求められる。【００６１】
第２１発明に係るコンピュータプログラムは、前記第一原理計算として、次式で表される波動関数汎関数を含んだ表式により定式化される拡張コーン・シャム理論に基づく最適化手法を用いることを特徴とする。【００６２】
【数９】
2007141942.tif 000011
【００６３】
ここで、Ψは多粒子波動関数、Ｔは運動エネルギ演算子、ｎΨは一電子密度、Ｆ［ｎ］はユニバーサルエネルギ汎関数、Ψ’は一電子密度ｎΨを再現する条件付で＜Ｔ＋Ｖred ＞を最小化する多粒子波動関数、ｖext は外部相互作用ポテンシャル、Ｖred は数式８により定義される演算子である。【００６４】
第２１発明にあっては、数式９により表される波動関数汎関数の最適化問題によって、クーロン系の基底状態エネルギ及びその一粒子密度が与えられる。【００６５】
第２２発明に係る記録媒体は、コンピュータに、物質の電子状態を第一原理計算により計算させるコンピュータプログラムが記録されているコンピュータでの読取りが可能な記録媒体において、コンピュータに、前記物質の物性を特徴付ける物理量に関して、該物理量の揺らぎを再現するモデルを決定させるステップと、コンピュータに、決定させたモデルを初期条件とした第一原理計算により前記物質の電子状態を計算させるステップとを有するコンピュータプログラムが記録されていることを特徴とする。【００６６】
第２２発明にあっては、物質の物性を特徴付ける物理量の揺らぎを再現するモデルを決定し、決定したモデルを初期条件とした第一原理計算により電子状態を求めるようにしている。したがって、縮退のない基底状態を示す電子系のシュレーディンガー方程式と等価なエネルギ及び一電子密度を与えることができるモデルであって、物理量の揺らぎを同時再現するモデルを、ユニバーサルエネルギ汎関数を用いて波動関数に関する汎関数として与えることができた場合、モデルに内在するパラメータが一意に定まる。【００６７】
第２３発明に係る記録媒体は、前記第一原理計算として、次式で表される残差交換相関エネルギ汎関数を含んだ表式により定式化される拡張コーン・シャム理論に基づく最適化手法を用いることを特徴とする。【００６８】
【数１０】
2007141942.tif 000012
【００６９】
ここで、ｎΨは一電子密度、Ｅxc［ｎΨ］は交換相関エネルギ汎関数、Φは相関が発生していない状態の多粒子波動関数、Ψ’は相関が発生した状態の多粒子波動関数、Ｔは運動エネルギ演算子、Ｖred は物理量の組｛Ｘi ｝により与えられる次の相関関数により定義される演算子である。【００７０】
【数１１】
2007141942.tif 000013
ただし、ｉ＝１，２，３，…，Ｎ（Ｎは計算すべき物理量の個数）、Ｗi は重みを表すパラメータ、小括弧内は物理量Ｘi と物理量Ｘi の平均値との差を表す。【００７１】
第２３発明にあっては、拡張コーン・シャム方程式において、残差交換相関エネルギを厳密な表式からの恒等変形により、通常の交換相関エネルギとそこからの差分として表記する。この差分の部分を、同じ電子密度を再現する揺らぎのある状態、無い状態の２つの状態における運動エネルギの差及びＶred の期待値として与え、この運動エネルギの変化を拡張コーン・シャム理論の枠内で与えることにより、量子モンテカルロ法、配置間相互作用法などの高精度計算手法を用いることなく電子状態の計算が可能となる。【００７２】
第２４発明に係る記録媒体は、前記第一原理計算として、次式で表される波動関数汎関数を含んだ表式により定式化される拡張コーン・シャム理論に基づく最適化手法を用いることを特徴とする。【００７３】
【数１２】
2007141942.tif 000014
【００７４】
ここで、Ψは多粒子波動関数、Ｔは運動エネルギ演算子、ｎΨは一電子密度、Ｆ［ｎ］はユニバーサルエネルギ汎関数、Ψ’は一電子密度ｎΨを再現する条件付で＜Ｔ＋Ｖred ＞を最小化する多粒子波動関数、ｖext は外部相互作用ポテンシャル、Ｖred は数式１１により定義される演算子である。【００７５】
第２４発明にあっては、上述した式により表される波動関数汎関数の最適化問題によって、クーロン系の基底状態エネルギ及びその一粒子密度が与えられる。
【発明の効果】
【００７６】
本発明による場合は、物質の物性を特徴付ける物理量の揺らぎを再現するモデルを決定し、決定したモデルを初期条件とした第一原理計算により電子状態を求めるようにしている。したがって、縮退のない基底状態を示す電子系のシュレーディンガー方程式と等価なエネルギ及び一電子密度を与えることができるモデルを記述するパラメータを一意に定めることができ、物質の電子状態を計算することができる。【００７７】
本発明による場合は、拡張コーン・シャム方程式において、残差交換相関エネルギを厳密な表式からの恒等変形により、通常の交換相関エネルギとそこからの差分として表記する。この差分の部分を、同じ電子密度を再現する揺らぎのある状態、無い状態の２つの状態における運動エネルギの差及びＶred の期待値として与え、この運動エネルギの変化を拡張コーン・シャム理論の枠内で与えることにより、量子モンテカルロ法、配置間相互作用法などの高精度計算手法を用いることなく電子状態の計算が可能となる。その結果、量子モンテカルロ法、配置間相互作用法と同程度の精度を確保しながら、局所密度近似法を越える高速計算手法を得ることができる。【００７８】
本発明による場合は、同じ電子密度を再現する揺らぎのある状態、無い状態の２つの状態における運動エネルギ及び全エネルギを評価することによって、拡張コーン・シャム理論の枠内で電子状態を計算することができる。【００７９】
本発明による場合は、参照系として計算した物理量の揺らぎを再現するようにモデルを決定する。したがって、このモデルを決定する過程において定まるパラメータを用いて多体系の波動関数を求めることができる。【００８０】
本発明による場合は、上述した式により表される波動関数汎関数の最適化問題によって、クーロン系の基底状態エネルギ及びその一粒子密度を与えることができる。【００８１】
本発明による場合は、交換エネルギ及び相関エネルギを含めたエネルギ汎関数のルジャンドル変換を密度変分して与えられる有効ポテンシャルの定義が可能になり、ハミルトニアンの運動エネルギ部分とこの有効ポテンシャルとからなる有効一体ハミルトニアン中で、相互作用する仮想的電子系が運動するという表現が得られる。したがって、有効ハミルトニアンの解として、例えば結晶を考えるときにブロッホ波を得た後に適当なバンドインデックスの組から一般化された逆フーリエ変換により直ちに導かれるワニエ関数から、ある局在軌道の表現を得ることができる。【００８２】
本発明による場合は、量子モンテカルロ法、トランスコリレイティッド法、配置間相互作用法、摂動計算・グリーン関数法、有効ポテンシャル法を用いて、参照系としての多体電子系における物理量の揺らぎを求めることができる。【００８３】
本発明による場合は、平面波基底展開法、補強平面波基底展開法、実空間、数値基底展開法、混合基底展開法、局在基底展開法、又は線型化ＭＴ軌道法を用いて、拡張コーン・シャム理論に基づいて物理量の揺らぎを求めることができる。【００８４】
本発明による場合は、物理量の揺らぎを正定値かつ有界な値域を持つ２次の相関関数により記述するため、縮体がない基底状態を示す電子系のシュレーディンガー方程式と等価なエネルギ及び一電子密度を与えることができるモデルであって、正定値有界な相関関数を同時再現するモデルを、ユニバーサルエネルギ汎関数を用いて与えることができた場合、内在するパラメータを一意に定めることができる。【００８５】
本発明による場合は、物理量を局所軌道上の電子に関する数演算子又は生成・消滅演算子として記述しているため、対応する拡張コーン・シャム模型は、第一原理的不純物アンダーソンモデル及び第一原理的ハバードモデルを含む強相関電子モデルとして得られる。【００８６】
本発明による場合は、物理量の揺らぎを再現するモデルを決定することにより、物質の構造因子（すなわち、原子レベルでの原子配置）が定まるため、他の物理量を全て無矛盾かつ前提条件なしで定めることができる。
【図面の簡単な説明】
【００８７】
【図１】実施の形態１に係る電子状態計算装置の内部構成を示すブロック図である。
【図２】電子状態計算装置が実行する処理の手順を示すフローチャートである。
【図３】参照系の算出処理を示すフローチャートである。
【図４】有効多体系の自己無撞着計算において電子状態計算装置が実行する処理の手順を説明するフローチャートである。
【図５】電子状態計算装置が実行する処理の手順を示すフローチャートである。
【図６】バルク系に対するＭＤ計算において電子状態計算装置が実行する処理の手順を説明するフローチャートである。
【図７】実施の形態３に係る電子状態計算装置の内部構成を示すブロック図である。
【図８】電子状態計算装置が実行する処理の手順を示すフローチャートである。
【図９】水素分子における１ｓ軌道上の揺らぎと運動エネルギとの関係を示すグラフである。
【図１０】水素分子における１ｓ軌道上の揺らぎとハートレーエネルギの期待値との関係を示すグラフである。
【符号の説明】
【００８８】
１０，２０  電子状態計算装置
１１，２１  ＣＰＵ
１２，２２  バス  
１３，２３  ＲＯＭ
１４，２４  ＲＡＭ
１５，２５  入出力ＩＦ
１５ａ，２５ａ  入力デバイス
１５ｂ，２５ｂ  出力デバイス
１６，２６  補助記憶装置
１７，２７  記憶装置
１７ａ，２７ａ  計算プログラム記憶領域
Ｍ１，Ｍ２  記録媒体  
【発明を実施するための最良の形態】
【００８９】
以下、実施の形態を図面を用いて具体的に説明する。
実施の形態１．
図１は実施の形態１に係る電子状態計算装置の内部構成を示すブロック図である。電子状態計算装置１０は、ＣＰＵ１１、ＲＯＭ１３、ＲＡＭ１４、入出力ＩＦ１５、補助記憶装置１６、及び記憶装置１７を備えており、これらの各ハードウェアがバス１２を介して互いに接続されている。ＲＯＭ１３には、ハードウェア各部を制御するための制御プログラムが格納されている。ＣＰＵ１１は、ＲＯＭ１３に格納されている制御プログラムをＲＡＭ１４上にロードして実行することにより、前述したハードウェア各部の制御を行う。入出力ＩＦ１５には、マウス、キーボードなどの入力デバイス１５ａ、及びＣＲＴ、液晶ディスプレイなどの出力デバイス１５ｂが接続される。入出力ＩＦ１５は、入力デバイス１５ａを通じて入力される情報、及び出力デバイス１５ｂへ出力すべき情報の入出力に関する制御を行う。【００９０】
補助記憶装置１６は、本実施の形態で説明する電子状態計算方法をコンピュータに実行させるためのコンピュータプログラムを記録したＦＤ、ＣＤ−ＲＯＭなどの記録媒体Ｍ１から前記コンピュータプログラムを読取るためのＦＤドライブ、ＣＤ−ＲＯＭドライブなどの読取装置を備える。補助記憶装置１６によって読取られたコンピュータプログラムは記憶装置１７に記憶される。ＣＰＵ１１は、記憶装置１７に記憶されたコンピュータプログラムをＲＡＭ１４上にロードして実行することにより、以下で説明する電子状態計算方法を実現する。【００９１】
なお、前記コンピュータプログラムを記録する記録媒体Ｍ１としては、前述したＦＤ及びＣＤ−ＲＯＭの他に、ＭＯ、ＭＤ、ＤＶＤ−ＲＯＭ等の光ディスク、ハードディスク等の磁気記録媒体、ＩＣカード、メモリカード、光カード等のカード型記録媒体、マスクＲＯＭ、ＥＰＲＯＭ（Erasable Programmable Read Only Memory）、ＥＥＰＲＯＭ（Electrically Erasable Programmable Read Only Memory）、フラッシュＲＯＭ等による半導体メモリの利用も可能である。また、インターネットを含む通信ネットワークを接続可能なシステムを構成することにより、通信ネットワークから前述のコンピュータプログラムをダウンロードするようにしてもよい。更に、ＲＯＭ１３の内部に前述のコンピュータプログラムを予め格納しておく構成であってもよい。【００９２】
記憶装置１７は、例えばハードディスクドライブであり、電子状態の計算に必要な初期データ、電子状態の計算により得られたデータなどが記憶されると共に、その記憶容量の一部は、各種計算手法に基づく計算プログラムが記憶される計算プログラム記憶領域１７ａとして利用される。この計算プログラム記憶領域１７ａに記憶される計算プログラムとしては、量子モンテカルロ法、トランスコリレイティッド法、配置間相互作用法、摂動計算・グリーン関数法、有効ポテンシャル法、平面波基底展開法、補強平面波基底展開法、線型化ＭＴ軌道法などに基づく計算プログラムが挙げられる。なお、電子状態の計算に用いられる計算手法としては、上記の計算手法に限定する必要はなく、将来的に電子状態計算分野で得られる各種精密計算手法を用いることも可能である。【００９３】
以下、本実施の形態に係る電子状態計算方法について説明する。なお、本実施の形態では、既知結晶構造に対する電子状態の最適化計算について説明する。図２は電子状態計算装置１０が実行する処理の手順を示すフローチャートである。電子状態計算装置１０は、まず、計算対象となる結晶の結晶構造、及びこの結晶構造を構成する原子の原子番号の入力を入力デバイス１５ａを通じて受付ける（ステップＳ１１）。次いで、電子状態計算装置１０は、参照系の算出処理を行う（ステップＳ１２）。この参照系の算出処理では、通常のコーン・シャム理論に基づく電子状態計算により、揺らぎを発生し得る原子とその周辺の原子集団とを特定する。【００９４】
図３は参照系の算出処理を示すフローチャートである。電子状態計算装置１０は、第一原理に基づいて電子状態の計算を行い（ステップＳ２１）、最適クラスター構造を決定する（ステップＳ２２）。すなわち、通常のコーン・シャム法による電子状態計算から揺らぎを発生し得る原子と、その周辺の原子集団を特定する。【００９５】
次いで、電子状態計算装置１０は、計算プログラム記憶領域１７ａに記憶されている計算プログラムを選択することにより、計算手法の選択を行い（ステップＳ２３）、選択した計算手法を用いて一電子密度及び軌道上の揺らぎを算出する（ステップＳ２４）。コーン・シャム方程式が与える量子状態は、相互作用のない独立フェルミオン系の解である単一スレータ積により表現される。しかし、コーン・シャム方程式が与える一粒子スペクトル構造には本来は意味がないとされており、むしろ仮想的な独立フェルミオン系の波動関数に関するエネルギ汎関数の最適化問題として定式化することが、厳密な表記法の一つを与えるということが知られている。ここでは、以下の定式化を採用する。【００９６】
【数１３】
2007141942.tif 000015
【００９７】
ここで、いわゆる交換相関エネルギ汎関数は、上記の等式から次のように定義される。【００９８】
【数１４】
2007141942.tif 000016
【００９９】
Ｇ0 ［Ψ］を用いた電子状態計算は既に広く方法論が知られており、各種物理量を直ちに計算することができる。参照系としては、現実物質を用いた実験により揺らぎが測定できる場合にはその結果を用いることもできる。実験が不可能な場合には、シュレーディンガー方程式を直接解くことができる高精度計算を用いて揺らぎの評価を行う。高精度計算はステップＳ２３で選択された手法が用いられる。その際、量子モンテカルロ法、トランスコリレイティッド法、配置間相互作用法、摂動計算・グリーン関数法などが使用される。また、実際の計算では、ステップＳ２２で決定された最適クラスター構造のデータ、並びに後述する有効多体系の自己無撞着計算により得られるモデルに内在するパラメータ、及び局在軌道のデータが使用される。【０１００】
コーン・シャム理論は密度汎関数理論に立脚して導入されるものであるため、そこで定まるものは一粒子密度であって、それ以外の物理量には個別に定式化しておく必要がある。なお、物理量Ｘi として局所軌道φi （ｒ）上の数演算子を用いた場合には、対応する拡張コーン・シャム模型は、第一原理的不純物アンダーソン模型、第一原理的ハバード模型を含む強相関電子模型として得られる。【０１０１】
説明を図２のフローチャートに戻す。電子状態計算装置１０は、参照系の算出処理の後、有効多体系について自己無撞着計算を行う（ステップＳ１３）。この自己無撞着計算では、本願発明者によって提案された拡張コーン・シャム理論に基づいて揺らぎ又は相関関数を求め、有効多体波動関数を決定する。なお、計算手法の詳細については後に詳述することとする。【０１０２】
次いで、電子状態計算装置１０は、参照系について算出した揺らぎと、自己無撞着計算を用いて有効多体系について計算した揺らぎとを比較し、両者が一致するか否かを判断する（ステップＳ１４）。両者が一致しないと判断した場合（Ｓ１４：ＮＯ）、有効モデルを記述するパラメータの探索範囲の再設定を行い（ステップＳ１５）、処理をステップＳ１３へ戻す。【０１０３】
両者が一致すると判断した場合（Ｓ１４：ＹＥＳ）、すなわち、揺らぎを再現するように決定される有効モデルを、第一原理計算における初期条件として利用することによって有効多体波動関数が求まった場合、既知結晶構造に対する電子状態の最適化を行う（ステップＳ１６）。【０１０４】
図４は有効多体系の自己無撞着計算において電子状態計算装置１０が実行する処理の手順を説明するフローチャートである。電子状態計算装置１０は、通常のコーン・シャム理論に基づいて初期電子密度を推定し（ステップＳ３１）、初期有効ポテンシャル及び初期局在軌道を算出する（ステップＳ３２，Ｓ３３）。【０１０５】
密度汎関数法の処方に基づく場合、交換エネルギ、相関エネルギまでも含めたエネルギ汎関数の密度変分として与えられる有効ポテンシャルが定義可能になる。すなわち、ハミルトニアンの運動エネルギ部分と有効ポテンシャルとからなる有効一体ハミルトニアン中で、仮想的電子系が運動するという表現が得られる。この場合、有効一体ハミルトニアンの解として、例えば、結晶を考えるときにはブロッホ波を得た後に適当なバンドインデックスの組から一般化された逆フーリエ変換により直ちに得られるワニエ関数から、ある局在軌道の表現が得られる。一般に、こうした局在軌道として、有効一体ハミルトニアンに現れる量子化されたエネルギ構造の中で占有数が揺らぐことが許される軌道の一群が発生する。本実施の形態では、揺らぎを評価する数演算子ｎi は、フェルミ順位近傍の準位を表現する規格直交化された局在軌道の一群φi に付随して定義する。したがって、φi （ｒ）上の密度・密度相関関数として次の量を定義する。【０１０６】
【数１５】
2007141942.tif 000017
【０１０７】
ここで、規格化された軌道φi （ｒ）に付随して生成・消滅演算子が与えられる。数演算子ｎi,σ  はこの軌道上の電子数を測定する。量子状態｜Ψ＞が与えられることによって、ｎi の期待値が＜Ψ｜ｎi,σ｜Ψ＞のように与えられる。【０１０８】
次いで、電子状態計算装置１０は、拡張コーン・シャム理論に基づいて電子密度、有効ポテンシャル、局所軌道を算出する（ステップＳ３４〜Ｓ３６）。第一原理計算のうち、密度汎関数法に基づくコーン・シャム理論に基づいて電子状態計算を与えるＧ0 ［Ψ］なる汎関数の精度を向上させるために、本実施の形態では、拡張コーン・シャム理論に基づくＧWi ［Ψ］なる汎関数と、残差交換相関エネルギ汎関数Ｅrxc ［ｎ］とを用いる。このとき、適宜の物理量Ｘi の揺らぎ、【０１０９】
【数１６】
2007141942.tif 000018
【０１１０】
について、適当に定めた参照系の値を再現するようにモデルのもつパラメータＷi を定める方法を用いる。通常のコーン・シャム理論による場合、Ｇ0 ［Ψ］によって与えられるモデルでは揺らぎを再現することが不可能であるが、拡張コーン・シャム理論では可能になるため、該当するモデルの精度は各種物理量の再現という点において著しく良くなる。ここで、揺らぎとして導入するものは正定値かつ有界な値域をもつ２次の相関関数を用いることができる。再現を要請する揺らぎとして導入する物理量の数と種類とを増やすことにより、モデルの精度を向上させることができる。拡張コーン・シャム理論に基づく汎関数ＧWi ［Ψ］は以下の式により与えられる。【０１１１】
【数１７】
2007141942.tif 000019
【０１１２】
ここで、Ｖred は次式で定義される物理量の組｛Ｘi ｝の与える相関関数を用いて、次のように定義される。【０１１３】
【数１８】
2007141942.tif 000020
【０１１４】
対応する残差交換相関関数は以下のように与えられる。【０１１５】
【数１９】
2007141942.tif 000021
【０１１６】
計算を実行する際には、残差交換相関エネルギ汎関数Ｅrxc ［ｎ］に近似を導入する必要がある。最終結果は、その近似精度に依存するため、残差交換相関エネルギ汎関数Ｅrxc ［ｎ］の評価の精度を向上させることが、計算精度を向上させることに繋がる。【０１１７】
以下、残差交換相関エネルギ汎関数Ｅrxc ［ｎ］の評価方法について説明する。残差交換相関エネルギ汎関数Ｅrxc ［ｎ］は、数式１９により厳密な値が与えられる。数式１９において、ｍｉｎ（Ψ’→ｎΨ）なる最適化過程は、波動関数の位相空間において定義される制限付き探索であるが、これは空間の各点において未定係数を設定したラグランジュ未定係数法により実行することができる。あるいは、ルジャンドル変換を行ってポテンシャルの最適化問題として再定義してから解いてもよい。【０１１８】
次いで、電子状態計算装置１０は収束判定を行い、計算が収束したか否かを判断する（ステップＳ３７）。多配置参照密度汎関数法に相当するＧWi［Ψ］の解は、スレータ積の線型結合で表される配置間相互作用の波動関数を解としてもつ。このとき、有効モデルの一体ハミルトニアン部分がΨから得られる一粒子密度を通して再決定されるため、一体部分が計算の途上で変化するが、与えられる位相空間の構造が不変となった段階で、有効相互作用が生じる部分空間を与えるアルゴリズムを全位相空間から部分空間を定めるユニタリ変換として与えることで位相空間を設定すると、有効相互作用の結果配置間の重みを決定しているものはＷi の組のみとなる。したがって、Ｗi の組とΨの構造とには一対一対応を見出すことが計算の収束段階において可能になる。【０１１９】
ステップＳ３７で計算が収束していないと判断した場合（Ｓ３７：ＮＯ）、電子状態計算装置１０は、処理をステップＳ３４へ戻す。また、ステップＳ３７で計算が収束したと判断した場合（Ｓ３７：ＹＥＳ）、決定される配置を最低エネルギ状態としてもつ有効模型がＷi の組を得ることで定められる。この解を先に得た上で、その自己無撞着解を固有解としてもつハミルトニアンを逆に決定してゆく方法を本実施の形態では採用している。今回の多配置参照密度汎関数法でこの方法が実施できる理由は、任意の揺らぎとそれに付随して発生する有効相互作用をもつ模型が少なくとも密度汎関数理論として成立するという本願発明者によって提案された理論に基づく。【０１２０】
このようにして最終的にＷi の組に対して対応する揺らぎが決定される。この物理量の特定された拡張コーン・シャム系を適切な拡張コーン・シャムモデルとして複数採用しておくことができる。Ｗi の組を変化させた場合の揺らぎの変化が、このようにして定まるため、最終的な参照系にフィットしたモデルが得られる。【０１２１】
また、自己無撞着解を決定してゆくプロセスにおいて、有効相互作用を設定する軌道の組とその上での相互作用の形態を決定するものは、すべて考慮する揺らぎの種類と数によって決定される。考慮する揺らぎの種類と数を変化させた場合に与えられるモデルが示す物理的性質（エネルギ、一電子密度、各種秩序変数、揺らぎ等）が変化しないところまで達していることが、有効模型として充分なものを得ているか否かの検証を行う指針となる。【０１２２】
次いで、電子状態計算装置１０は、計算プログラム記憶領域１７ａに記憶されている計算プログラムを選択することにより、計算手法の選択を行い（ステップＳ３８）、選択した計算手法を用いて有効多体波動関数の係数行列を決定することにより、バルク秩序を決定する（ステップＳ３９）。ここで選択される計算手法としては、厳密対角化法、量子モンテカルロ法、数値繰り込み群法、動的平均場近似法などが挙げられる。【０１２３】
実施の形態２．
実施の形態１では、既知結晶構造に対する電子状態の最適化計算について説明したが、分子動力学の手法を用いて系のダイナミクスを計算するようにしてもよい。なお、電子状態計算装置１０の内部構成については実施の形態１と全く同様であるため、その説明を省略することとする。【０１２４】
図５は電子状態計算装置１０が実行する処理の手順を示すフローチャートである。電子状態計算装置１０は、まず、計算対象となる結晶の結晶構造、及びこの結晶構造を構成する原子の原子番号の入力を入力デバイス１５ａを通じて受付ける（ステップＳ５１）。次いで、電子状態計算装置１０は、参照系の算出処理を行う（ステップＳ５２）。この参照系の算出処理は、図３のフローチャートにより説明される処理内容と全く同様であり、通常のコーン・シャム理論に基づく電子状態計算により、揺らぎの発生し得る原子とその周辺の原子集団とを特定する。【０１２５】
次いで、電子状態計算装置１０は、有効多体系について自己無撞着計算を行う（ステップＳ５３）。この自己無撞着計算は、図４のフローチャートにより説明される処理内容と全く同様であり、本願発明者によって提案された拡張コーン・シャム理論に基づいて密度・密度相関関数を求め、有効多体波動関数を決定する。【０１２６】
次いで、電子状態計算装置１０は、参照系について算出した揺らぎと、自己無撞着計算を用いて有効多体系について計算した揺らぎとを比較し、両者が一致するか否かを判断する（ステップＳ５４）。両者が一致しないと判断した場合（Ｓ５４：ＮＯ）、有効モデルを記述するパラメータの探索範囲の再設定を行い（ステップＳ５５）、処理をステップＳ５３へ戻す。【０１２７】
両者が一致すると判断した場合（Ｓ５４：ＹＥＳ）、すなわち、揺らぎを再現するように決定される有効モデルを第一原理計算における初期条件として利用することにより、ある時間発展ステップにおける有効多体波動関数が求まった場合、バルク系に対するＭＤ計算（MD : Molecular Dynamics）を行う（ステップＳ５６）。なお、計算手法については後に詳述することとする。【０１２８】
次いで、電子状態計算装置１０は、時系列に対する計算が完了したか否かを、ＭＤ計算の結果に基づいて判断する（ステップＳ５８）。計算が完了していないと判断した場合（Ｓ５８：ＮＯ）、電子状態計算装置１０は、参照系決定用の構造データを更新した上で（ステップＳ５９）、処理をステップＳ５２へ戻す。時系列に対する計算が完了したと判断した場合（Ｓ５８：ＹＥＳ）、本フローチャートによる処理を終了する。【０１２９】
図６はバルク系に対するＭＤ計算において電子状態計算装置１０が実行する処理の手順を説明するフローチャートである。電子状態計算装置１０は、拡張コーン・シャム理論に基づいて得られた電子密度、有効ポテンシャル、及び局在軌道を用いて原子間力及び応力を算出する（ステップＳ６１，Ｓ６２）。【０１３０】
次いで、電子状態計算装置１０は、算出した原子間力及び応力に基づいて物質構造変化率を決定する（ステップＳ６３）。そして、電子状態計算装置１０は、微小時間後の物質構造を決定し、ＭＤの時間発展ステップを更新する（ステップＳ６４）。【０１３１】
実施の形態３．
実施の形態１及び実施の形態２では、量子モンテカルロ法、配置間相互作用法などの高精度計算手法を利用して参照系について計算した揺らぎと、拡張コーン・シャム理論に基づく自己無撞着計算を利用して計算した揺らぎとを比較し、両者が一致するように有効モデルを記述するパラメータを決定する構成とした。しかしながら、このような手法では、量子モンテカルロ法、配置間相互作用法などの計算と比較することによって初めて電子状態の計算精度が高まるため、拡張コーン・シャム理論に基づく自己無撞着計算を単独で用いた場合には計算精度を上げる指標が存在しない。そこで、本実施の形態では、拡張コーン・シャム理論に基づく自己無撞着計算を単独で用いて電子状態を計算する手法について説明する。【０１３２】
図７は実施の形態３に係る電子状態計算装置の内部構成を示すブロック図である。電子状態計算装置２０は、実施の形態１で説明した電子状態計算装置１０と同様の構成であり、ＣＰＵ２１、ＲＯＭ２３、ＲＡＭ２４、入出力ＩＦ２５、補助記憶装置２６、及び記憶装置２７を備え、これらの各ハードウェアがバス２２を介して互いに接続されている。また、入出力ＩＦ２５には、マウス、キーボードなどの入力デバイス２５ａ、及びＣＲＴ、液晶ディスプレイなどの出力デバイス２５ｂが接続される。入出力ＩＦ２５は、入力デバイス２５ａを通じて入力される情報、及び出力デバイス２５ｂへ出力すべき情報の入出力に関する制御を行う。【０１３３】
補助記憶装置２６は、本実施の形態で説明する電子状態計算方法をコンピュータに実行させるためのコンピュータプログラムを記録したＦＤ、ＣＤ−ＲＯＭなどの記録媒体Ｍ２から前記コンピュータプログラムを読取るためのＦＤドライブ、ＣＤ−ＲＯＭドライブなどの読取装置を備える。補助記憶装置２６によって読取られたコンピュータプログラムは記憶装置２７に記憶される。ＣＰＵ２１は、記憶装置２７に記憶されたコンピュータプログラムをＲＡＭ２４上にロードして実行することにより、以下で説明する電子状態計算方法を実現する。【０１３４】
なお、前記コンピュータプログラムを記録する記録媒体Ｍ２としては、前述したＦＤ及びＣＤ−ＲＯＭの他に、ＭＯ、ＭＤ、ＤＶＤ−ＲＯＭ等の光ディスク、ハードディスク等の磁気記録媒体、ＩＣカード、メモリカード、光カード等のカード型記録媒体、マスクＲＯＭ、ＥＰＲＯＭ（Erasable Programmable Read Only Memory）、ＥＥＰＲＯＭ（Electrically Erasable Programmable Read Only Memory）、フラッシュＲＯＭ等による半導体メモリの利用も可能である。また、インターネットを含む通信ネットワークを接続可能なシステムを構成することにより、通信ネットワークから前述のコンピュータプログラムをダウンロードするようにしてもよい。更に、ＲＯＭ２３の内部に前述のコンピュータプログラムを予め格納しておく構成であってもよい。【０１３５】
記憶装置２７は、例えばハードディスクドライブであり、電子状態の計算に必要な初期データ、電子状態の計算により得られたデータなどが記憶されると共に、その記憶容量の一部は、拡張コーン・シャム理論による自己無撞着計算法、物理量評価法、従来型のコーン・シャム理論による物理量評価法などに基づく計算プログラムが記憶される計算プログラム記憶領域２７ａとして利用される。【０１３６】
以下、本実施の形態に係る電子状態計算方法について説明する。図８は電子状態計算装置２０が実行する処理の手順を示すフローチャートである。電子状態計算装置２０は、まず、揺らぎの大きさを可変とした拡張コーン・シャム理論に基づき自己無撞着計算を行う（ステップＳ７１）。【０１３７】
拡張コーン・シャム方程式において、基底状態のエネルギは以下のように記述することができる。【０１３８】
【数２０】
2007141942.tif 000022
【０１３９】
ここで、ｎGSは基底状態における一電子密度、Ｆはユニバーサルエネルギ汎関数、ｖext は外部相互作用ポテンシャルである。また、Ｔは運動エネルギ演算子、Ｖred は数式１８で定義される揺らぎ、ｎi,σは電子に関する数演算子、Ｅrxc ［ｎΨ ］は残差交換相関エネルギ汎関数である。【０１４０】
対応する残差交換相関エネルギ汎関数は、数式１９に示した厳密な表式からの恒等変形により、通常の交換相関エネルギとそこからの差分として次式により表記することができる。【０１４１】
【数２１】
2007141942.tif 000023
【０１４２】
ここで、Ｅxc［ｎΨ］は交換相関エネルギ汎関数であり、従来の近似理論であるＬＤＡ、ＧＧＡ等を用いて評価可能な部分である。また、Φは相関が発生していない状態の波動関数、Ψ’は相関が発生した状態の波動関数である。すなわち、数式２１の右辺第２項及び第３項は、同じ電子密度を再現する揺らぎのある状態及び揺らぎのない状態の２つの状態における運動エネルギの差及びＶred の期待値を表す。【０１４３】
このように、Ψを有効多電子方程式により決定したうえで、交換相関エネルギ汎関数の部分をＬＤＡ，ＧＧＡを用いて評価し、数式２１の運動エネルギのずれを評価することにより基底状態のエネルギＥ0 を評価することができる。すなわち、数式２１の評価において、Ψ’の与える運動エネルギがΦの与える運動エネルギより大きい場合、相関の発生していない状態よりもエネルギの低い状態が相関の発生により生じることを示している。拡張コーン・シャム方程式により、相関発生のある状態により系の表現が適切に行われる。【０１４４】
このような状況が起こり得ることは、次の水素分子に関する拡張コーン・シャム方程式の計算結果から示される。図９は水素分子における１ｓ軌道上の揺らぎと運動エネルギとの関係を示すグラフである。横軸は１ｓ軌道上の揺らぎの大きさを表し、縦軸は運動エネルギの期待値を表す。このグラフは、原子間距離を０．６Åから２．０Åまで０．２Åずつ変えた場合のエネルギ変化をシンボルを異ならせて表したものである。【０１４５】
また、図１０は水素分子における１ｓ軌道上の揺らぎとハートレーエネルギの期待値との関係を示すグラフである。図９のグラフと同様に、横軸は１ｓ軌道上の揺らぎの大きさ、縦軸はハートレーエネルギの期待値を表し、原子間距離を０．６Åから２．０Åまで０．２Åずつ変えた場合のエネルギ変化をシンボルを異ならせて表したものである。【０１４６】
図１０のグラフからは、原子間距離が２．０Åでは密度がほぼ揺らぎによらず、ハートレーエネルギがほぼ一定であることが分かる。また、図９のグラフからは、運動エネルギの相対値は原子間距離で傾向が逆転することが分かる。すなわち、原子間距離が短い場合、揺らぎが小さくになるにつれて運動エネルギの期待値が小さくなるのに対し、原子間距離が１．０Å付近でその傾向が逆転し、揺らぎが小さくなるにつれて運動エネルギが上昇する傾向が見られ、揺らぎの抑制が発生することが分かる。【０１４７】
したがって、上述の数式２１において揺らぎの大きさを変化させて運動エネルギの期待値を算出することにより、揺らぎの抑制が発生したか否かを判断することができる。電子状態計算装置２０は、数式２１の有効相互作用を可変として自己無撞着計算を行うことにより、揺らぎの抑制が発生したか否かを判断し（ステップＳ７２）、揺らぎの抑制が発生したと判断した場合（Ｓ７２：ＹＥＳ）、拡張コーン・シャム理論による物理量の評価を行い（ステップＳ７３）、揺らぎの抑制が発生していないと判断した場合（Ｓ７２：ＮＯ）、従来型のコーン・シャム理論による物理量の評価を行う（ステップＳ７４）。【０１４８】
このように、本実施の形態では、量子モンテカルロ法、配置間相互作用法などの高精度計算手法による参照系の計算手続きをスキップすることができるため、計算速度の著しい高速化を図ることができる。また、参照系を必要とせずに最終ターゲットに直接適用できるため、大規模系の計算に直接応用することが可能となる。その結果、著しい高精度化を図ることができ、従来では困難であったフェルミ面近傍の電子状態の記述が可能となった。更に、本実施の形態は、残差交換相関エネルギ汎関数法の評価法を与えるため、密度汎関数法における主要な近似法であったＬＤＡを越え、かつ計算ステップとして事実上同規模の計算手順に抑えた実用的計算方法として、高精度かつ高速計算方法を実現した。
【図１】
2007141942.tif 000036
【図２】
2007141942.tif 000037
【図３】
2007141942.tif 000038
【図４】
2007141942.tif 000039
【図５】
2007141942.tif 000040
【図６】
2007141942.tif 000041
【図７】
2007141942.tif 000042
【図８】
2007141942.tif 000043
【図９】
2007141942.tif 000044
【図１０】
2007141942.tif 000045
【国際調査報告】
2007141942.tif 000046
2007141942.tif 000047
2007141942.tif 000048
2007141942.tif 000049