JP2015141533A - 電子状態算出装置、電子状態算出方法、及び電子状態算出プログラム - Google Patents

Abstract

【課題】トライ・アンド・エラーを繰り返さずに、より広範な物質に対して、実験により検出された電子の状態を再現することができる電子状態算出装置を提供すること。
【解決手段】物質に含まれる電子のスピン状態を示す値を、仮想磁場を示す値を利用して繰り返し算出する算出部と、前記スピン状態を示す値が所定値に収束するか否かを判定する判定部と、前記判定部による判定結果に基づいて、前記仮想磁場に関する値を変更する変更部と、を備える電子状態算出装置である。
【選択図】図１

Description

この発明は、電子状態算出装置、電子状態算出方法、及び電子状態算出プログラムに関する。

物性の解析や予測を行うため、物質（結晶や分子だけではなく、原子も含む）中の電子状態（例えば、電子の波動関数や電子密度、電子の持つエネルギー等）を知るための電子状態計算の研究が行われている。

これに関連し、近年のコンピュータの性能の向上により、第一原理計算と呼ばれる精度の高い電子状態計算が、広く物質の研究に用いられている。例えば、第一原理計算は、分子の安定構造・反応機構の解析や、固体材料の力学的特性・光学的特性・電気的特性の予測に用いられる（非特許文献１参照）。第一原理計算では、物質の電子状態(すなわち電子の波動関数またはコーン・シャム軌道、電子密度、電子の持つエネルギー、スピン状態)が算出される。

より具体的には、第一原理計算は、初期値として与えた電子密度に基づいて有効ポテンシャルを算出し、算出された有効ポテンシャルに基づいて電子状態を算出し、算出した電子状態から再度電子密度を算出するという演算を繰り返し行う。この繰り返しの演算は、自己無撞着場計算と呼ばれ、算出される電子状態が収束するまで行う。このような演算によって算出された電子状態（例えば、電子の波動関数等）に基づいて、分子の安定構造・反応機構の解析や、固体材料の力学的特性・光学的特性・電気的特性の予測等を行うことができる。

ところで、電子は上向きスピンと下向きスピンとのうちいずれか一方のスピンを有しており、物質ごとに上向きスピンを有する電子と下向きスピンを有する電子の数が異なる場合がある。物質中の上向きスピンを有する電子と下向きスピンを有する電子の数との個数差を、電子のスピン状態と称する。電子のスピン状態は、電子の状態の一例である。物質の物性は、物質に含まれる電子のスピン状態の変化に応じて変化する。物質中の電子のスピン状態は、例えば、ＥＳＲ（Electron Spin Resonance）等の実験により検出することができる。

ここで、第一原理計算でスピン分極した物質を扱う場合、上向き及び下向きスピンの電子の状態をそれぞれ独立に計算し、上述の実験により検出された電子のスピン状態を再現する。しかし、物質によっては、第一原理計算によって実験により検出されたスピン状態を再現できない場合がある。

この問題を解決するため、第一原理計算において、物質に仮想的な磁場が印加されていると仮定し、電子のスピン状態を計算する手法が知られている（非特許文献２参照）。この手法では、自己無撞着場計算の初回に、物質に対して仮想的な磁場を印加し、自己無撞着場計算による繰り返し演算の数が増える毎に、仮想的な磁場を減衰させる。そして、電子のスピン状態が収束する時に、この仮想的な磁場を無視できるほど小さくする。

"A Chemist’s Guide to Density Functional Theory", W.Koch, M.C.Holthausen, WILEY-VCH "Elk code manual", http://elk.sourceforge.net/elk.pdf

しかしながら、非特許文献２に記載の手法では、所望のスピン分極が起こらず、実験により検出される電子のスピン状態を再現することができない場合があった。また、最初に印加すべき仮想的な磁場の大きさや、適切な磁場の減衰速度が物質によって変化することが知られている。このため、算出する電子のスピン状態を所望の状態に収束させるには、トライ・アンド・エラーを繰り返さなければならず、時間のロスが生じるという問題もあった。さらに、物質によっては、所望の状態を再現する最適値が見つからない場合もあった。

そこで本発明は、上記従来手法の問題に鑑みてなされたものであり、トライ・アンド・エラーを繰り返さずに、より広範な物質に対して、実験により検出される電子の状態を再現することができる電子状態算出装置、電子状態算出方法、及び電子状態算出プログラムを提供する。

本発明の一態様は、物質に含まれる電子のスピン状態を示す値を、仮想磁場を示す値を利用して繰り返し算出する算出部と、前記スピン状態を示す値が所定値に収束するか否かを判定する判定部と、前記判定部による判定結果に基づいて、前記仮想磁場に関する値を変更する変更部と、を備える電子状態算出装置である。

そこで本発明は、トライ・アンド・エラーを繰り返さずに、より広範な物質に対して、実験により検出された電子の状態を再現することができる電子状態算出装置、電子状態算出方法、及び電子状態算出プログラムを提供することができる。

第１の実施形態に係る電子状態算出装置１の機能構成の一例を示す図である。 電子状態算出装置１のハードウェア構成の一例を示す図である。 算出制御部１０により実行される、算出制御部１０が計算スピン状態を算出し、目標スピン状態に収束させる処理の流れの一例を示すフローチャートである。 第１の実施形態に係る電子状態算出装置１との比較対象となる電子状態算出装置Ｘによる計算スピン状態を算出する処理の流れの一例を説明するフローチャートである。 仮想磁場を印加した場合に、計算モデル中の電子に生じる変化の様子を示す図である。 第２の実施形態に係る電子状態算出装置１の機能構成の一例を示す図である。 第２の実施形態に係る電子状態算出装置１によってＥｕＯと呼ばれるＮａＣｌ型結晶に含まれる電子の計算スピン状態を算出されるまでの様子を示す図である。 第３の実施形態に係る電子状態算出装置１によってＳｒ３１Ｓ３２：Ｅｕ２＋のと呼ばれるＮａＣｌ型結晶のＳｒＳのうち、Ｓｒの一部をＥｕで置換した物質に含まれる電子の計算スピン状態を算出されるまでの様子を示す図である。

＜概要＞
まず、以下に示す実施形態に係る電子状態算出装置１の概要を説明し、その後により詳細な実施形態を説明する。
電子状態算出装置１は、物質に含まれる電子のスピン状態を示す値を、仮想磁場を示す値を利用して繰り返し算出し、スピン状態を示す値が所定の範囲内に存在するか否かを判定し、その判定結果に基づいて、仮想磁場に関する値を変更する。このため、電子状態算出装置１は、トライ・アンド・エラーを繰り返さずに、より広範な物質に対して、実験により検出される電子の状態を再現することができる。

本実施形態において、物質に含まれる電子とは、分子や固体材料等に含まれた電子を示す。また、再現される電子状態とは、電子のスピン状態を含むが、電子の波動関数や電子密度、電子の持つエネルギー等も含む。また、電子のスピン状態とは、物質中の上向きスピンを有する電子（以下、上向きスピン電子と称する）の個数と、下向きスピンを有する電子（以下、下向きスピン電子と称する）の個数の差分として定義される量である。なお、同じ構造が周期的に並ぶ結晶などの周期的境界条件を課すことができる物質では単位格子中の全電子の上向きスピン電子と下向きスピン電子の差をスピン状態としても良い。また、本実施形態において、仮想磁場に関する値とは、仮想磁場の値であるが、これに代えて、仮想磁場に結合する結合定数等であってもよい。

より具体的には、電子状態算出装置１は、第一原理計算による自己無撞着場計算によって算出される電子のスピン状態（以下、計算スピン状態と称する）と実験により検出される電子のスピン状態（以下、目標スピン状態と称する）とが所定の範囲に存在するような電子状態に収束させる。電子状態算出装置１は、繰り返される自己無撞着場計算の１回の演算ごとに、計算スピン状態が目標スピン状態に近い値となるまで、仮想磁場の値を増大させる。目標スピン状態に近い値とは、後述する所定の第１範囲に収まる値である。

そして、電子状態算出装置１は、計算スピン状態が目標スピン状態に近い値になった後、目標スピン状態に収束するまでの間、仮想磁場の値を減少させる。この仮想磁場を減少させる処理に関して、電子状態算出装置１は、最終的に計算スピン状態と目標スピン状態とが所定の範囲に存在する時、仮想磁場の値が無視できるほど小さくなるように減少させる。
このように、電子状態算出装置１は、計算スピン状態を常にチェックしながら、磁場を逐次変動させることで、計算スピン状態と目標スピン状態とが所定の範囲に存在するような電子状態に収束させることができる。

＜第１の実施形態＞
以下、本発明の第１の実施形態について、図面を参照して説明する。図１は、第１の実施形態に係る電子状態算出装置１の機能構成の一例を示す図である。電子状態算出装置１は、例えば、算出制御部１０と、記憶部３０と、入力受付部３２と、表示部３４とを備える。第１の実施形態において、電子状態算出装置１は、孤立系を仮定したスピン分極している分子に含まれる電子に関する計算スピン状態を算出する。孤立系とは、真空中に１原子や１分子だけが孤立して存在することを仮定したモデルであり、例えば、有機低分子等の１分子単独の電子状態を計算する場合に用いられるモデルである。

なお、電子状態算出装置１は、孤立系を仮定したスピン分極している分子に含まれる電子に代えて、孤立系を仮定したスピン分極している原子等に含まれる電子でもよい。また、孤立系を仮定した分子や原子ではない固体材料中に含まれる電子の計算スピン状態の算出については、第２の実施形態及び第３の実施形態において説明する。

電子状態算出装置１は、密度汎関数理論（Density Functional Theory：DFT）を用いた第一原理計算によって電子の計算スピン状態を算出する。密度汎関数理論は、ホーエンベルク・コーンの第１の定理及び第２の定理に基づいて、ピエール・ホーエンベルクとヴァルター・コーンによって提案された理論である。第１の実施形態において、電子状態算出装置１は、孤立系を仮定したスピン分極している原子に含まれる電子の計算スピン状態を算出する。そのため、電子状態算出装置１は、スピン分極を考慮した密度汎関数理論によって、電子状態を算出する。

スピン分極を考慮した密度汎関数理論において、上向きスピン電子及び下向きスピン電子の状態は、電子の波動関数に代えて、コーン・シャム軌道と呼ばれる関数によってそれぞれ表される。また、前述した通り、第１の実施形態において、電子は、孤立系を仮定したスピン分極している原子に含まれている。そのため、上向きスピン電子及び下向きスピン電子のコーン・シャム軌道は、原子軌道による線形結合法（Linear Combination of Atomic Orbital：LCAO）や、原子に強く束縛された近似（Tight Binding Approximation）に基づいた原子軌道の線形和によってそれぞれ表すことができる。電子状態算出装置１は、この原子軌道の線形和によって表された上向きスピン電子及び下向きスピン電子のコーン・シャム軌道に基づいて、スピン分極を考慮した電子状態を算出する。

より具体的には、ＬＣＡＯやＴＢＡを仮定した場合、上向きスピン電子及び下向きスピン電子のコーン・シャム軌道は、以下に示した式（１）及び式（２）によってそれぞれ表される。なお、以下の文中では、便宜上、Ψ_{KS}^{up,n}(r)によって式（１）の左辺に示される上向きスピン電子のコーン・シャム軌道を示し、Ψ_{KS}^{dn,n}(r)によって式（２）の左辺に示される下向きスピン電子のコーン・シャム軌道を示す。ここで、rは、計算モデル中の座標系における位置を表す変数である。これらの記載方法と同様に、以下では、「A_{…}」と記載した場合、中括弧内の文字が先頭のアルファベット（この例では「A」）又はギリシャ文字に対する下付き添え字を表すものとし、「^{…}」と記載した場合、中括弧内の文字が先頭のアルファベット又はギリシャ文字に対する上付き添え字を表すものとする。また、式（１）及び式（２）の右辺のギリシャ文字ファイを、以下の文中では、φと記載する。なお、以降全てハートリーの原子単位系を用いる。

式（１）及び式（２）のコーン・シャム軌道の添え字nは、上向きスピン電子及び下向きスピン電子のそれぞれについて、エネルギーの小さいものから順に付与した電子の番号である。ここで、φ_{i_{atom},i_{orbital}}は、i_{atom}番目の原子のi_{orbital}番目の原子軌道である。また、C_{i_{atom},i_{orbital}}^{up,n}やC_{i_{atom},i_{orbital}}^{dn,n}は、n番目の電子に対応する各基底の展開係数（重み）である。また、N_{orbital}^{i_{atom}}は、各原子の原子軌道のうち、基底として用いるものの数である。また、N_{atom}は、計算モデル中の原子数である。式（１）及び式（２）に示したように、ＬＣＡＯやＴＢＡを仮定した場合、上向きスピン電子及び下向きスピン電子のコーン・シャム軌道は、計算モデル（孤立系）の各原子の原子軌道φ_{i_{atom},i_{orbital}}の線形和で表される。電子状態算出装置１は、これらのコーン・シャム軌道に基づいて、計算スピン状態を算出する。

ここより、電子状態算出装置１の各機能部について説明する。算出制御部１０は、初期値算出部１２と、有効ポテンシャル算出部１４と、密度算出部１８と、スピン状態算出部２０と、スピン状態判定部２２と、収束判定部２４と、仮想磁場変更部２６とを備える。

初期値算出部１２は、入力受付部３２により受け付けられたユーザーからの初期条件の入力に基づいて、原子中の上向きスピン電子の密度（以下、上向きスピン電子密度と称する）及び下向きスピン電子の密度（以下、下向きスピン電子密度と称する）をそれぞれ第一原理計算（自己無撞着場計算）に用いる初期値として算出する。初期条件とは、例えば、式（１）及び式（２）に含まれるN_{orbital}^{i_{atom}}やN_{atom}、計算モデルにおける原子の位置を示す情報（例えば、i_{orbital}と原子とを対応付ける情報等）、スピン状態をコントロールしたい原子（ターゲット原子）を示す情報（例えば、i_{atom}とターゲット原子とを対応づける情報等）、ターゲット原子の目標スピン状態、密度汎関数の種類を示す情報や収束条件等を含む。

有効ポテンシャル算出部１４は、初期値算出部１２により算出された上向きスピン電子密度及び下向きスピン電子密度に基づいて、それぞれに対応する有効ポテンシャルを算出する。より具体的には、有効ポテンシャル算出部１４は、記憶部３０により予め記憶された複数種類の密度汎関数を表す情報のうち、初期条件によって指定された密度汎関数を表す情報を読み込む。そして、有効ポテンシャル算出部１４は、読み込んだ密度汎関数と、初期値算出部１２により算出された上向きスピン電子密度及び下向きスピン電子密度に基づいて、それぞれに対する有効ポテンシャルを算出する。また、有効ポテンシャル算出部１４は、仮想磁場が電子に印加されているとして、有効ポテンシャルを算出する。有効ポテンシャルを算出する処理についての詳細は後述する。

密度算出部１８は、有効ポテンシャル算出部１４により算出された上向きスピン電子及び下向きスピン電子に対応する有効ポテンシャルに基づいて、上向きスピン電子及び下向きスピン電子のコーン・シャム軌道に対するコーン・シャム方程式を生成する。そして、密度算出部１８は、これらのコーン・シャム方程式を解き、コーン・シャム軌道とその固有エネルギーを算出する。密度算出部１８は、算出されたこれらのコーン・シャム軌道とその固有エネルギーに基づいて、上向きスピン電子密度と下向きスピン電子密度とをそれぞれ算出する。コーン・シャム方程式に基づく上向きスピン電子密度と下向きスピン電子密度とを算出するための処理についての詳細は後述する。

スピン状態算出部２０は、密度算出部１８により算出された上向きスピン電子密度及び下向きスピン電子密度に基づいて、原子中の上向きスピン電子の個数及び下向きスピン電子の個数をそれぞれ算出する。なお、密度汎関数理論において、上向きスピン電子の個数及び下向きスピン電子の個数は、後述するようにコーン・シャム軌道の絶対値の２乗に基づいて定義された量であり、必ずしも整数ではない。スピン状態算出部２０は、算出された個数に基づいて、計算スピン状態を算出する。計算スピン状態を算出するための処理についての詳細は後述する。
スピン状態判定部２２は、スピン状態算出部２０により算出された計算スピン状態が、所定の第１範囲内に含まれるか否かを判定するための第１判定処理を行う。所定の第１範囲とは、例えば、目標スピン状態から±０．５程度の範囲とするが、これに代えて、より広い範囲でも狭い範囲でもよい。スピン状態判定部２２は、判定部の一例である。

収束判定部２４は、スピン状態判定部２２により計算スピン状態が第１範囲内であると判定された場合、当該演算ルーチンにおいて計算スピン状態を算出するために使用した上向きスピン電子密度と、前回の演算ルーチンにおいて計算スピン状態を算出するために使用した上向きスピン電子密度との差分が所定の閾値Ｙ未満か否かを判定する。

同様に、収束判定部２４は、当該演算ルーチンにおいて計算スピン状態を算出するために使用した下向きスピン電子密度と、前回の演算ルーチンにおいて計算スピン状態を算出するために使用した下向きスピン電子密度との差分が所定の閾値Ｙ未満か否かを判定する。以下、この判定を、収束判定処理と称する。なお、収束判定部２４は、収束判定処理を、上向きスピン電子密度と下向きスピン電子密度とのうちいずれか一方又は両方に対して行ってもよい。以下、収束判定部２４は、収束判定処理を上向きスピン電子密度と下向きスピン電子密度との両方に対して行うものとして説明する。

仮想磁場変更部２６は、算出制御部１０に含まれる他の各機能部の連携によって電子のスピン状態が算出される演算ルーチンごとに、スピン状態判定部２２による判定結果に基づいて、有効ポテンシャル算出部１４が有効ポテンシャルを算出するために使用する仮想磁場の値を増大させる。また、仮想磁場変更部２６は、スピン状態判定部２２により計算スピン状態が第１範囲内であると判定された場合、算出制御部１０に含まれる他の各機能部の連携によって電子のスピン状態が算出される演算ルーチンごとに、有効ポテンシャル算出部１４が有効ポテンシャルを算出するために使用する仮想磁場の値を減少させる。仮想磁場の値を変更する処理についての詳細は後述する。

記憶部３０は、例えば、ＨＤＤ（Hard Disk Drive）やＳＳＤ（Solid State Drive）、ＥＥＰＲＯＭ（Electrically Erasable Programmable Read-Only Memory）、ＲＯＭ（Read-Only Memory）、ＲＡＭ（Random Access Memory）などを含み、電子状態算出装置１が処理する情報やプログラムを格納する。なお、記憶部３０は、電子状態算出装置１に内蔵されるものに代えて、外付け型の記憶装置でもよい。

入力受付部３２は、キーボードやマウス、タッチパッド、その他の入力装置である。なお、入力受付部３２は、表示部として機能してもよく、さらに、タッチパネルとして構成されてもよい。
表示部３４は、例えば、液晶ディスプレイパネル、あるいは、有機ＥＬ（ElectroLuminescence）ディスプレイパネルである。

次に、図２を参照することで、電子状態算出装置１のハードウェア構成について説明する。図２は、電子状態算出装置１のハードウェア構成の一例を示す図である。電子状態算出装置１は、例えば、記憶部３０と、入力受付部３２と、表示部３４と、ＣＰＵ（Central Processing Unit）３６とを備える。ＣＰＵ３６は、記憶部３０に格納された各種プログラムを実行する。なお、電子状態算出装置１は、通信部を備えることで、他の装置と無線や有線によって通信可能に接続されていてもよい。

以下、図３を参照し、算出制御部１０が計算スピン状態を算出し、目標スピン状態に収束させる処理について説明する。図３は、算出制御部１０により実行される、算出制御部１０が計算スピン状態を算出し、目標スピン状態に収束させる処理の流れの一例を示すフローチャートである。

まず、初期値算出部１２は、入力受付部３２から、ユーザーが入力した前述の初期条件を取得する（ステップＳ１００）。次に、初期値算出部１２は、取得した初期条件に基づいて、原子中の上向きスピン電子密度の初期値ρ_{initial}^{up}及び下向きスピン電子密度の初期値ρ_{initial}^{dn}を算出する（ステップＳ１１０）。初期値算出部１２は、これらの初期値を、オクテット則等に基づいた酸化数から算出してもよいし、Gaussian等の密度汎関数理論計算用ソフトウェアと同じ方法を用いて算出してもよいし、他の既知の方法で算出してもよい。

次に、有効ポテンシャル算出部１４は、初期値算出部１２により算出された初期値に基づいて、有効ポテンシャルを算出する（ステップＳ１２０）ここで、有効ポテンシャルを算出するための処理について説明する。有効ポテンシャル算出部１４は、有効ポテンシャルを算出するためにまず、初期条件に基づいて記憶部３０からユーザーに指定された密度汎関数E_{XC}[ρ^{up}(r),ρ^{dn}(r)]を表す情報を読み込む。

密度汎関数E_{XC}[ρ^{up}(r),ρ^{dn}(r)]は、パウリの排他率による電子間の反発、量子揺らぎ、配置間相互作用等に起因する電子の交換相関エネルギーを表す。以下、有効ポテンシャル算出部１４は、密度汎関数E_{XC}[ρ^{up}(r),ρ^{dn}(r)]として、局所スピン密度近似（Local Spin Density Approximation）や一般化された密度勾配近似（Generalized Gradient Approximation）によって表される密度汎関数を表す情報を読み込むものとして説明する。有効ポテンシャル算出部１４は、これに代えて、例えば、ハイブリッド汎関数であるＢ３ＬＹＰ等を読み込んでもよい。その場合、電子状態算出装置１は、ハートリー・フォック理論と密度汎関数理論を組み合わせることによって、計算スピン状態を算出する。

有効ポテンシャル算出部１４は、密度汎関数E_{XC}[ρ^{up}(r),ρ^{dn}(r)]を表す情報を読み込むと、密度汎関数E_{XC}[ρ^{up}(r),ρ^{dn}(r)]の上向きスピン電子密度ρ^{up}による変分ν_{xc}^{up}[ρ^{up}(r),ρ^{dn}(r)]を算出する。また、有効ポテンシャル算出部１４は、密度汎関数E_{XC}[ρ^{up}(r),ρ^{dn}(r)]の下向きスピン電子密度ρ^{dn}による変分ν_{xc}^{dn}[ρ^{up}(r),ρ^{dn}(r)]を算出する。なお、有効ポテンシャル算出部１４は、計算機のメモリ・計算処理軽減のため密度汎関数E_{XC}[ρ^{up}(r),ρ^{dn}(r)]の変分ν_{xc}^{up}[ρ^{up}(r),ρ^{dn}(r)]とν_{xc}^{dn}[ρ^{up}(r),ρ^{dn}(r)]を記憶部３０から読み込み、上向きスピン電子密度ρ^{up}(r)と下向きスピン電子密度ρ^{dn}(r)を変分ν_{xc}^{up}とν_{xc}^{dn}に代入しても良い。以降、変分ν_{xc}^{up}とν_{xc}^{dn}とをまとめて交換相関ポテンシャルと呼ぶ。

上向きスピン電子及び下向きスピン電子に対する有効ポテンシャルは交換相関ポテンシャルと、上向きスピン電子及び下向きスピン電子によるクーロンポテンシャル（式（３）及び式（４）における右辺第１項）と、原子核の静電ポテンシャルν_{nucl}(r)と、内殻電子を静電ポテンシャルとして近似した擬ポテンシャルν_{pseudo}(r)と、計算モデル中の原子に印加した仮想磁場B_{initial}^{ext}(r)とによって、以下の式（３）及び（４）に示したように表される。なお、仮想磁場の印加は、計算モデル中のすべての原子に印加してもよいし、その中の一部であるターゲット原子にのみ印加してもよい。以下では、有効ポテンシャル算出部１４は、仮想磁場をターゲット原子にのみ印加するものとして説明する。

有効ポテンシャル算出部１４は、初期値算出部１２により算出された原子中の上向きスピン電子密度の初期値ρ_{initial}^{up}及び下向きスピン電子密度の初期値ρ_{initial}^{dn}を式（３）及び式（４）に代入することで、上向きスピン電子及び下向きスピン電子に対する有効ポテンシャルを算出する。

次に、密度算出部１８は、有効ポテンシャル算出部１４により算出された有効ポテンシャルに基づいて、新たな上向きスピン電子密度及び下向きスピン電子密度を算出する（ステップＳ１３０）。ここで、上向きスピン電子密度及び下向きスピン電子密度を算出するための処理について説明する。上向きスピン電子密度は、スピン分極を考慮した密度汎関数理論において、コーン・シャム軌道Ψ_{KS}^{up,n}に基づいて算出される。同様に、下向きスピン電子密度は、コーン・シャム軌道Ψ_{KS}^{dn,n}に基づいて算出される。これらのコーン・シャム軌道は、有効ポテンシャルに基づいて、以下の式（５）及び式（６）に示したコーン・シャム方程式における固有値問題を解くことによって算出される。

ここで、ε_{n}^{up}は、コーン・シャム軌道Ψ_{KS}^{up,n}を固有関数とする固有値方程式における固有値であり、コーン・シャム軌道Ψ_{KS}^{up,n}の固有エネルギーを表す。また、ε_{n}^{dn}は、コーン・シャム軌道Ψ_{KS}^{dn,n}を固有関数とする固有値方程式における固有値であり、コーン・シャム軌道Ψ_{KS}^{dn,n}の固有エネルギーを表す。これらの固有エネルギーも、コーン・シャム軌道と同様に、式（５）及び式（６）における固有値問題を解くことによって算出される。

このように算出されたコーン・シャム軌道に基づいて、上向きスピン電子密度及び下向きスピン電子密度は、以下に示した式（７）及び（８）から算出される。

ここで、a^{up,n}は、各コーン・シャム軌道における上向きスピン電子の占有数である。また、a^{dn,n}は、各コーン・シャム軌道における下向きスピン電子の占有数である。Neは、分子内の電子の総数である。ただし、擬ポテンシャルに含まれるとして露わには計算しない電子を除いた数とする。a^{up,n}及びa^{dn,n}は、対応する固有エネルギーε_{n}^{up}またはε_{n}^{dn}が小さいコーン・シャム軌道に対して下から順に上向きスピン電子を詰めていき、Ｎ番目までのコーン・シャム軌道を占有軌道とすることで算出される。各コーン・シャム軌道が、電子に占有されている場合は、a^{up,n}またはa^{dn,n}の値は１であり、非占有の場合は０となる。スピンの向きに関係なくエネルギーが小さい軌道から順に電子を詰めていくことで、上下スピンの占有数が同じになることもあるが、違いが出ることもある。すなわち、上向きスピン電子密度ρ^{up}（ｒ）と下向きスピン電子密度ρ^{dn}（ｒ）は、式（７）及び式（８）に基づいて算出する。

密度算出部１８は、式（５）及び式（６）によって上向きスピン電子及び下向きスピン電子のコーン・シャム軌道を算出する。そして、密度算出部１８は、算出された上向きスピン電子及び下向きスピン電子のコーン・シャム軌道に基づいて、式（７）及び式（８）から上向きスピン電子密度及び下向きスピン電子密度を算出する。

次に、スピン状態算出部２０は、電子の計算スピン状態を算出する（ステップＳ１５０）。ここで、電子の計算スピン状態を算出するための処理について説明する。電子の計算スピン状態は、ターゲット原子に含まれる上向きスピン電子の数n_{i_{target}}^{up}と下向きスピン電子の数n_{i_{target}}^{dn}の差S_C= n_{i_{target}}^{up}- n_{i_{target}}^{dn}として定義される。ターゲット原子に含まれる上向きスピン電子の数n_{i_{target}}^{up}と下向きスピン電子の数n_{i_{target}}^{dn}は、以下の式（９）及び式（１０）から算出される。

スピン状態算出部２０は、上記の式（９）及び式（１０）に基づいて、ターゲット電子に含まれる計算スピン状態S_C=n_{i_{target}}^{up}-n_{i_{target}}^{dn}を算出する。

次に、スピン状態判定部２２は、スピン状態算出部２０により算出された計算スピン状態が、前述した所定の第１範囲内に含まれるか否かを判定するスピン状態判定処理を行う（ステップＳ１５０）。仮想磁場変更部２６は、スピン状態判定部２２によって計算スピン状態が所定の第１範囲内に含まれていないと判定された場合（ステップＳ１５０−Ｎｏ）、仮想磁場に第１所定値を乗算することで、仮想磁場の値を増大させる（ステップＳ１７０）。その後、算出制御部１０は、ステップＳ１２０に戻り有効ポテンシャル算出部１４によって当該演算ルーチンで算出された上向きスピン電子密度及び下向きスピン電子密度に基づいて有効ポテンシャルを算出し直す。第１所定値とは、例えば、１．２程度である。

収束判定部２４は、スピン状態判定部２２によって計算スピン状態が所定の第１範囲内に含まれていると判定された場合（ステップＳ１５０−Ｙｅｓ）、収束判定処理を行う（ステップＳ１６０）。ここで、収束判定処理について説明する。収束判定部２４は、当該演算ルーチンにおいて計算スピン状態を算出するために使用した上向きスピン電子密度と、前回の演算ルーチンにおいて計算スピン状態を算出するために使用した上向きスピン電子密度との差分が所定の閾値Ｙ未満か否かを判定する。また、収束判定部２４は、当該演算ルーチンにおいて計算スピン状態を算出するために使用した下向きスピン電子密度と、前回の演算ルーチンにおいて計算スピン状態を算出するために使用した下向きスピン電子密度との差分が所定の閾値Ｙ未満か否かを判定する。

収束判定部２４は、これらの判定結果が両方とも所定の閾値Ｙ未満であると判定した場合（ステップＳ１６０−Ｙｅｓ）、計算スピン状態が収束したと判断して処理を終了する。仮想磁場変更部２６は、収束判定部２４によるこれらの判定結果のうち少なくともいずれか一方が所定の閾値Ｙ以上であると判定された場合（ステップＳ１６０−Ｎｏ）、仮想磁場に第２所定値を乗算することで、仮想磁場の値を減少させる（ステップＳ１８０）。その後、算出制御部１０は、ステップＳ１２０に戻り有効ポテンシャル算出部１４によって当該演算ルーチンで算出された上向きスピン電子密度及び下向きスピン電子密度に基づいて有効ポテンシャルを算出し直す。第２所定値は、例えば、０．８程度である。

なお、収束判定部２４は、電子密度に代えて、電子のトータルエネルギーの変化量で収束判定処理を行ってもよいし、電子密度と電子のトータルエネルギーとの両方で収束判定処理を行ってもよい。また、仮想磁場変更部２６は、ステップＳ１７０において、第１所定値を乗算するのに代えて、第１所定値を加算してもよい。また、仮想磁場変更部２６は、ステップＳ１７０において、仮想磁場を増大させる割合や加算する値を演算ルーチンごとに変化させてもよい。同様に、仮想磁場変更部２６は、ステップＳ１８０において、第２所定値を乗算するのに代えて、第２所定値を減算してもよい。また、仮想磁場変更部２６は、ステップＳ１８０において、仮想磁場を減少させる割合や減算する値を演算ルーチンごとに変化させてもよい。

以上のように、電子状態算出装置１は、ステップＳ１５０によって計算スピン状態を常にチェックしながら、ステップＳ１７０において仮想磁場を逐次変動させることで、計算スピン状態を目標スピン状態に収束させることができる。ここで、図４を参照することで、第１の実施形態に係る電子状態算出装置１との比較対象となる電子状態算出装置Ｘによる計算スピン状態を算出する処理の流れの一例を説明する。図４は、第１の実施形態に係る電子状態算出装置１との比較対象となる電子状態算出装置Ｘによる計算スピン状態を算出する処理の流れの一例を説明するフローチャートである。

図４に示したフローチャートにおいて、ステップＳ２００〜ステップＳ２４０までは、図３に示したフローチャートにおけるステップＳ１００〜ステップＳ１４０と同様の処理であるため、説明を省略する。電子状態算出装置Ｘは、ステップＳ２４０において算出された計算スピン状態に対して、収束判定処理を行う（ステップＳ２５０）。電子状態算出装置Ｘは、上向きスピン電子密度と下向きスピン電子密度に対する収束判定の結果が両方とも所定の閾値Ｙ未満であると判定した場合（ステップＳ２５０−Ｙｅｓ）、計算スピン状態が収束したと判断して処理を終了する。電子状態算出装置Ｘは、上向きスピン電子密度と下向きスピン電子密度に対する収束判定の判定結果のうち少なくともいずれか一方が所定の閾値Ｙ以上であると判定された場合（ステップＳ２５０−Ｎｏ）、仮想磁場に第２所定値を乗算することで、仮想磁場の値を減少させる（ステップＳ２６０）。

このように、電子状態算出装置Ｘは、第１の実施形態に係る電子状態算出装置１が行うスピン状態判定処理（ステップＳ１５０の処理）とそれに伴う仮想磁場の値を増大させる処理（ステップＳ１７０の処理）を行わないため、計算スピン状態が目標スピン状態に収束しない場合がある。一方で、第１の実施形態に係る電子状態算出装置１は、電子状態算出装置Ｘとは異なり、図３に示したステップＳ１５０において行われるスピン状態判定部２２によるスピン状態判定処理と、それに伴うステップＳ１７０の仮想磁場を増大させる処理のために、より確実に計算スピン状態を目標スピン状態に収束させることができる。

ここで、図５を参照することで、電子状態算出装置１が行う仮想磁場を逐次変動させる処理を、計算モデル中の電子に生じる変化の様子を示すことで説明する。図５は、仮想磁場を印加した場合に、計算モデル中の電子に生じる変化の様子を示す図である。図５（Ａ）には、仮想磁場が印加されていない場合のターゲット原子中の電子の様子を示す。仮想磁場が印加されていない場合、スピンＳ１を有する電子のエネルギーＥ１及びスピンＳ２を有する電子のエネルギーＥ２は、図５（Ａ）に示したように縮退する。この縮退は、仮想磁場を印加することによるゼーマン効果によって解かれる。なお、図５（Ａ）において、スピンＳ１は上向き、スピンＳ２は下向きであったとする。

図５（Ｂ）には、仮想磁場を印加した場合の計算モデル中の電子の様子を示す。図５（Ｃ）には、図５（Ｂ）より強い仮想磁場を印加した場合の計算モデル中の電子の様子を示す。図５（Ｂ）及び（Ｃ）に示したように、計算モデル中のスピンＳ１を有する電子のエネルギーＥ１とスピンＳ２を有する電子のエネルギーＥ２は、仮想磁場に比例するゼーマンエネルギーΔＥの分だけ縮退が解かれる。そして、スピンＳ１及びスピンＳ２は、エネルギーの低い電子のスピンに揃いやすくなるため、図５（Ｂ）及び図５（Ｃ）に示したように、両方とも上向きとなる。ここで、図５（Ｂ）及び図５（Ｃ）では、エネルギーＥ１＜エネルギーＥ２であるものとする。この結果として、仮想磁場が印加された計算モデル中の電子は、高スピン状態になりやすい。

従って、図４に示した電子状態算出装置Ｘでは、計算スピン状態を算出する際、最初に印加した仮想磁場から減少させていくのみであったため、計算スピン状態を目標スピン状態に収束させることができない場合があった。第１の実施形態に係る電子状態算出装置１では、計算スピン状態が目標スピン状態に近づくまで最初に印加した仮想磁場を増大させ、その後、計算スピン状態が目標スピン状態に収束するように仮想磁場を減少させていくため、より確実に計算スピン状態を目標スピン状態に収束させることができる。

以上説明したように、第１の実施形態に係る電子状態算出装置１は、原子に含まれる電子の計算スピン状態を、仮想磁場の値を利用して繰り返し算出し、計算スピン状態が目標スピン状態に収束するか否かを判定し、その判定結果に基づいて、仮想磁場の値を変更する。このため、電子状態算出装置１は、トライ・アンド・エラーを繰り返さずに、より広範な物質に対して、実験により検出された電子の状態を再現することができる。

また、電子状態算出装置１は、計算スピン状態を示す値が目標スピン状態に収束するか否かを、計算スピン状態が第１範囲内に含まれるか否かを指標として判定し、判定結果に基づいて、仮想磁場の値を増大させるため、この判定結果に基づく仮想磁場の増大を行う処理が無い場合と比べて、より確実に計算スピン状態を目標スピン状態に収束させることができる。

なお、仮想磁場は分子を構成する原子全てに一様にかかっているとしたが、一部の原子のみにかかっているとしてもよい。特にターゲット原子にのみ仮想磁場がかかっているとして有効ポテンシャルを計算することで、よりターゲット原子のスピン状態を制御しやすくなる。これは、式（３）と式（４）における仮想磁場が有効ポテンシャルに寄与する以下の式（１１）の項を式（９）と式（１０）を使って算出した計算スピン状態S_C=n_{i_{target}}^{up}-n_{i_{target}}^{dn}を使って代用することが望ましい。

これにより、ターゲット原子のみに磁場をかけることができるため非常に効率よく計算スピン状態を目標スピン状態に収束させられる。当然、分子中の特定原子にのみ磁場がかかった状態は現実では考えられないが、本計算では計算途中にのみ特定原子に仮想磁場を印加し、最終的な収束状態では仮想磁場を無視できるほど小さくするため、算出した電子状態は、実験より得たデータと比較・議論すべき妥当性を持つ。

＜第２の実施形態＞
以下、本発明の第２の実施形態について、図面を参照して説明する。第２の実施形態に係る電子状態算出装置１は、孤立系における原子に含まれる電子の計算スピン状態を算出するのに代えて、固体材料等（とくに結晶）の周期境界条件を課すことができる物質に含まれる電子の計算スピン状態を算出する。全体の構成については、図１を援用し、同じ機能部に対して同一の符号を付して説明する。また、同じ機能を有する機能部の処理については、説明を省略する。

まず、周期境界条件を課すことができる物質に含まれる電子に対応するコーン・シャム軌道と、それを利用した計算スピン状態の算出方法について説明する。周期境界条件を課すことができる物質に含まれる電子のコーン・シャム軌道は、以下に示した式（１２）及び式（１３）のように、平面波基底の線形和によって表わされる。

以下では、式中に登場する太字のアルファベットを、アルファベットに「（太字）」と付記することで表す。k（太字）は、ｋ番目の固体結晶格子の点（以下、ｋ点と称する）を表すベクトルであり、Ｇ（太字）は、平面波基底の波数ベクトルを表す（詳細な理論的背景は、「”Inorganic Chemistry”, Shriver&Atkins, Oxford University Press」に記載された内容に従うものとする）。第１の実施形態と異なり、第２の実施形態では、コーン・シャム軌道は、ｋ点毎に定義される。

そして、式（１２）及び式（１３）に示したコーン・シャム軌道が従うコーン・シャム方程式は、以下に示した式（１４）及び式（１５）のように表される。

コーン・シャム軌道とその固有エネルギーは、上記の式（１４）及び式（１５）を解くことで、ｋ点毎に算出される。そして、ｋ点毎の上向きスピン電子密度及び下向きスピン電子密度は、算出されたｋ点毎のコーン・シャム軌道に基づいて、以下に示した式（１６）及び式（１７）から算出される。

ここで、w_{k(太字)}は、各ｋ点における重みである。また、a^{up,n,k(太字)}及びa^{dn,n,k(太字)}は、上向きスピン電子及び下向きスピン電子によるコーン・シャム軌道の占有数である。これらの電子密度に基づいて、計算スピン状態は、以下に示した式（１８）から算出される。

ここで、R(太字)_{T}は、ターゲット原子の原子核の位置を示す。また、b_Tは、ターゲット原子の原子半径を示すが、例えば、ターゲット原子と最近接原子との間の距離の１／２でもよいし、上述の「”Inorganic Chemistry”, Shriver&Atkins, Oxford University Press」に記載されたinonic radiiやatomic radii等でもよい。

以下、第２の実施形態に係る電子状態算出装置１の機能構成について説明する。図６は、第２の実施形態に係る電子状態算出装置１の機能構成の一例を示す図である。電子状態算出装置１は、例えば、算出制御部１０ａと、記憶部３０と、入力受付部３２と、表示部３４とを備える。算出制御部１０ａは、例えば、初期値算出部１２と、有効ポテンシャル算出部１４と、密度算出部１８ａと、スピン状態算出部２０ａと、スピン状態判定部２２と、収束判定部２４と、仮想磁場変更部２６とを備える。

密度算出部１８ａは、上記の式（１４）及び式（１５）に基づいて、ｋ点毎のコーン・シャム軌道とその固有値を算出する。そして、算出されたコーン・シャム軌道に基づいて、上記の式（１６）及び式（１７）から、上向きスピン電子密度及び下向きスピン電子密度を算出する。
スピン状態算出部２０ａは、密度算出部１８ａによって算出された上向きスピン電子密度及び下向きスピン電子密度に基づいて、上記の式（１８）から計算スピン状態を算出する。

なお、第２の実施形態に係る電子状態算出装置１が計算スピン状態を算出し、目標スピン状態に収束させる処理の流れは、図３に示したフローチャートと同様であるため、説明を省略する。

ここで、図７を参照することで、第２の実施形態に係る電子状態算出装置１によって計算スピン状態を算出した具体例を示す。図７は、第２の実施形態に係る電子状態算出装置１によってＮａＣｌ型結晶である酸化ユウロピウムＥｕＯ内の電子状態が算出されるまで、計算スピン状態の様子を示す図である。なお、ユウロピウムＥｕは、４ｆ軌道に７つの電子を持ち、そのすべてが平行なスピン（例えば、すべて上向きスピン）を有することが知られている。すなわち、実験により検出されるスピン状態は７となる。従って、酸化ユウロピウムＥｕＯに含まれる電子の計算スピン状態を算出する場合、電子状態算出装置１による計算スピン状態を収束させたい目標スピン状態は７となる。

図７には、自己無撞着場計算の演算ルーチン毎に算出された計算スピン状態をプロットしたグラフを示した。このグラフでは、縦軸が計算スピン状態を表し、横軸が自己無撞着場計算の演算ルーチンの回数を表す。実線で示される曲線は、第２の実施形態に係る電子状態算出装置１によって算出された計算スピン状態のプロットを曲線で結んだものである。また、点線で示される曲線は、前述した電子状態算出装置Ｘによって算出された計算スピン状態のプロットを曲線で結んだものである。

図７に示したように、電子状態算出装置Ｘによって算出された計算スピン状態は、ほぼ０に収束している。一方、電子状態算出装置１によって算出された計算スピン状態は、ほぼ７に収束している。なお、図７に示した演算において、電子状態算出装置１は、図３に示したフローチャートにおけるステップＳ１７０において、仮想磁場を一定の割合（ここでは、１．２倍）で増大させ、ステップＳ１８０において、仮想磁場を一定の割合（ここで、０．８倍）で減少させた。より具体的には、電子状態算出装置１は、８回目の演算ルーチンまでは仮想磁場を一定の割合で増大させ、９回目以降の演算ルーチンにおいて仮想磁場を一定の割合で減少させた。図７に示した結果から、電子状態算出装置１は、電子状態算出装置Ｘと比べて、より確実に計算スピン状態を目標スピン状態に収束させることができることが分かる。

以上説明したように、第２の実施形態に係る電子状態算出装置１は、ｋ点毎にコーン・シャム軌道を算出し、算出されたコーン・シャム軌道に基づいて計算スピン状態を算出するため、固体結晶中の電子に対しても、より確実に計算スピン状態を目標スピン状態に収束させることができ、第１の実施形態と同様の効果を得ることができる。

＜第３の実施形態＞
以下、本発明の第３の実施形態について、図面を参照して説明する。第３の実施形態に係る電子状態算出装置１は、第２の実施形態と同様に固体結晶等の周期境界条件を課すことができる物質に含まれる電子の計算スピン状態を算出するが、第２の実施形態とは異なる基底関数を用いて計算スピン状態を算出する。全体の構成については、図６を援用し、同じ機能部に対して同一の符号を付して説明する。

第３の実施形態において、電子状態算出装置１は、ＦＬＡＰＷ（Full-potential Linearized Augumented Plane Wave）という計算方法を用いて、計算スピン状態を算出する。ＦＬＡＰＷでは、コーン・シャム軌道の基底を、各原子の原子核を中心とした球状の領域（Muffin-Tin:MT）と、固体結晶中の領域ＭＴに属さないすべての領域（Interstitial Region:I）に分けて定義される。以下では、固体結晶中に含まれる原子αの領域ＭＴを、領域ＭＴαと称して説明する。また、以下では、固体結晶中の領域ＭＴに属さないすべての領域を領域Ｉと称して説明する。

ＦＬＡＰＷにおいて、コーン・シャム軌道の基底は、以下に示した式（１９）のように各領域ＭＴαと、領域Ｉとで別々に定義される。以下では、式中に登場するハット付のアルファベットを、アルファベットに「（ハット）」と付記することで表す。また、式中に登場するドット付のアルファベットを、アルファベットに「（ドット）」と付記することで表す。

ここで、Y_{L}(r(ハット)^{α})は、球面調和関数を表す。r^{α}は、原子核を中心とする極座標の動径方向成分を表す。r(ハット)^{α}は、原子核を中心とする極座標の角度成分を表す。L(=(l,m))は、球面調和関数における方位量子数lと磁気量子数mを表す。ν_{1}^{α}(r^{α})は、動径方向の原子軌道成分を表す。ν(ドット)_{1}^{α}(r^{α})は、ν_{1}^{α}(r^{α})をr^{α}で微分した導関数を表す。a_{LG}^{α}とb_{LG}^{α}は、領域Ｉと領域ＭＴαの境界において、波動関数が微分可能になるように決められた複素数係数を表す。

このように、ＦＬＡＰＷにおいて、コーン・シャム軌道の基底は、領域ＭＴでは原子軌道で表され、その他の領域Ｉでは平面波で表される。そして、コーン・シャム軌道は、それら２つ（原子軌道と平面波）を連続的に接続されるような（微分可能となるような）基底によって表される。上記の式（１９）によって表されたコーン・シャム軌道は、領域Ｉでは緩やかに変化し、原子核近傍で激しい変化を起こすような軌道を高精度で表現できる。

さらに、ＦＬＡＰＷでは、特定の原子に局在化したコーン・シャム軌道が、以下に示した局所軌道（local orbital:lo）と呼ばれる基底によって表される。

なお、局所軌道とは、上記の式（２０）に示したように、ある１つの原子αの領域ＭＴα内でのみ値を持ち、他の領域では０となる基底である。局所軌道は、係数a_{Llo}^{α}とb_{Llo}^{α}を調整することで、各領域の境界が微分可能となるように接続される。また、第１の実施形態や第２の実施形態において擬ポテンシャルで置き換えていた内殻電子についても、局所軌道によって各電子の軌道を１つの基底で表すことができる。従って、ＦＬＡＰＷでは、基底の数をほぼ増やすことなく、擬ポテンシャル法を用いずに、すべての電子について計算するフルポテンシャル計算を行うことができる。

上向きスピン電子及び下向きスピン電子に対するコーン・シャム軌道は、上記の式（１９）及び式（２０）に示した基底を用いて、以下の式（２１）及び式（２２）のように表される。

C_{k(太字),G(太字)}^{up,n}、C_{k(太字),G(太字)}^{dn,n}、D_{k(太字),α,lo}^{up,n}、及びD_{k(太字),α,lo}^{dn,n}は、各基底の係数である。第３の実施形態に係る電子状態算出装置１は、上記の式（２１）及び式（２２）に示したコーン・シャム軌道を用いて、第２の実施形態に係る電子状態算出装置１と同様の機能部により、同様の流れで計算スピン状態を算出する。なお、第３の実施形態に係る電子状態算出装置１では、図３に示したステップＳ１４０において計算スピン状態を算出する際、原子半径b_Tを各原子の領域ＭＴの半径とするのが好ましい。

ここで、図８を参照することで、第３の実施形態に係る電子状態算出装置１によって計算スピン状態を算出した具体例を示す。図８は、第３の実施形態に係る電子状態算出装置１によってＳｒ_３１Ｓ_３２：Ｅｕ^２＋の組成で表されるＮａＣｌ型結晶の硫化ストロンチウムＳｒＳのうち、ストロンチウムＳｒの一部をユウロピウムＥｕで置換した物質に含まれる電子の計算スピン状態を算出されるまでの様子を示す図である。前述したようにユウロピウムＥｕは、４ｆ軌道に７つの電子を持ち、そのすべてが平行なスピン（例えば、すべて上向きスピン）を有することが知られている。そのため実験により検出されるスピン状態は７となる。従って、Ｓｒ_３１Ｓ_３２：Ｅｕ^２＋に含まれる電子の計算スピン状態を算出する場合、電子状態算出装置１による計算スピン状態を収束させたい目標スピン状態は７である。

図８には、自己無撞着場計算の演算ルーチン毎に算出された計算スピン状態をプロットしたグラフを示した。このグラフでは、縦軸が計算スピン状態を表し、横軸が自己無撞着場計算の演算ルーチンの回数を表す。実線で示される曲線は、第３の実施形態に係る電子状態算出装置１によって算出された計算スピン状態のプロットを曲線で結んだものである。また、点線で示される曲線は、前述した電子状態算出装置Ｘによって算出された計算スピン状態のプロットを曲線で結んだものである。

図８に示したように、電子状態算出装置Ｘによって算出された計算スピン状態は、ほぼ０に収束している。一方、電子状態算出装置１によって算出された計算スピン状態は、ほぼ７に収束している。なお、図８に示した演算において、電子状態算出装置１は、図３に示したフローチャートにおけるステップＳ１７０において、仮想磁場を一定の割合（ここでは、１．２倍）で増大させ、ステップＳ１８０において、仮想磁場を一定の割合（ここで、０．８倍）で減少させた。より具体的には、電子状態算出装置１は、３５回目の演算ルーチンまでは仮想磁場を一定の割合で増大させ、３６回目以降の演算ルーチンにおいて仮想磁場を一定の割合で減少させた。図８に示した結果から、電子状態算出装置１は、電子状態算出装置Ｘと比べて、より確実に計算スピン状態を目標スピン状態に収束させることができることが分かる。

以上説明したように、第３の実施形態に係る電子状態算出装置１は、ＦＬＡＰＷにおけるコーン・シャム軌道に基づいて計算スピン状態を算出するため、第１の実施形態及び第２の実施形態と同様の効果を得ることができる。

なお、図１、６における電子状態算出装置１を構成する各部の機能を実現するためのプログラムを、コンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することにより電子状態算出装置１の実施を行ってもよい。なお、ここでいう「コンピュータシステム」とは、ＯＳ（Operation System）や周辺機器等のハードウェアを含むものとする。

また、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ−ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間の間、動的にプログラムを保持するもの、その場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリのように、一定時間プログラムを保持しているものも含むものとする。また上記プログラムは、前述した機能の一部を実現するためのものであっても良く、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであっても良い。

以上、この発明の実施形態を、図面を参照して詳述してきたが、具体的な構成はこの実施形態に限られるものではなく、この発明の要旨を逸脱しない限り、変更、置換、削除等されてもよい。

（１）本発明の一態様は、物質に含まれる電子のスピン状態を示す値を、仮想磁場を示す値を利用して繰り返し算出する算出部と、前記スピン状態を示す値が所定値に収束するか否かを判定する判定部と、前記判定部による判定結果に基づいて、前記仮想磁場に関する値を変更する変更部と、を備える電子状態算出装置である。
（２）また、本発明の他の態様は、前記算出部は、前記スピン状態を示す値が前記所定値に収束するか否かを判定するための指標値をさらに算出する、（１）に記載の電子状態算出装置である。

（３）また、本発明の他の態様は、前記判定部は、前記指標値が、所定範囲内の値であるか否かを判定し、前記変更部は、当該判定結果に基づいて前記仮想磁場に関する値を変更する、（２）に記載の電子状態算出装置である。
（４）また、本発明の他の態様は、前記仮想磁場に関する値は、前記電子に印加された仮想磁場の大きさであり、前記変更部は、前記指標値が、前記所定範囲内の値ではないと判定した場合、前記仮想磁場の大きさを増大させる、（３）に記載の電子状態算出装置である。

（５）また、本発明の他の態様は、前記指標値は、前記電子のスピン状態を示す値である、（２）から（４）のうちいずれか一項に記載の電子状態算出装置である。
（６）また、本発明の他の態様は、前記電子のスピン状態を示す値は、上向きのスピンを有する前記電子の個数と、下向きのスピンを有する前記電子の個数との差分値である、（１）から（５）のうちいずれか一項に記載の電子状態算出装置である。
（７）また、本発明の他の態様は、前記物質とは、所定の原子であり、前記電子のスピン状態を示す値は、前記所定の原子に含まれる上向きのスピンを有する電子の個数と、前記所定の原子に含まれる下向きのスピンを有する電子の個数との差分値である、（１）から（６）のうちいずれか一項に記載の電子状態算出装置である。

（８）また、本発明の他の態様は、前記算出部は、コーン・シャム方程式に基づいて、前記スピン状態を示す値を算出する、（１）から（７）のうちいずれか一項に記載の電子状態算出装置である。
（９）また、本発明の他の態様は、物質に含まれる電子のスピン状態を示す値を、仮想磁場を示す値を利用して繰り返し算出し、前記スピン状態を示す値が所定値に収束するか否かを判定し、前記判定した結果に基づいて、前記仮想磁場に関する値を変更する、電子状態算出方法である。

（１０）また、本発明の他の態様は、コンピュータに物質に含まれる電子のスピン状態を示す値を、仮想磁場を示す値を利用して繰り返し算出させ、前記スピン状態を示す値が所定値に収束するか否かを判定させ、前記判定させた結果に基づいて、前記仮想磁場に関する値を変更させる、電子状態算出プログラムである。

１ 電子状態算出装置、１０ 算出制御部、１２ 初期値算出部、１４ 有効ポテンシャル算出部、１６ 方程式算出部、１８ 密度算出部、２０ スピン状態算出部、２２ スピン状態判定部、２４ 収束判定部、２６ 仮想磁場変更部、３０ 記憶部、３２ 入力受付部、３４ 表示部、３６ ＣＰＵ

Claims (10)

1. 物質に含まれる電子のスピン状態を示す値を、仮想磁場を示す値を利用して繰り返し算出する算出部と、
   前記スピン状態を示す値が所定値に収束するか否かを判定する判定部と、
   前記判定部による判定結果に基づいて、前記仮想磁場に関する値を変更する変更部と、
   を備える電子状態算出装置。
2. 前記算出部は、前記スピン状態を示す値が前記所定値に収束するか否かを判定するための指標値をさらに算出する、
   請求項１に記載の電子状態算出装置。
3. 前記判定部は、前記指標値が、所定範囲内の値であるか否かを判定し、
   前記変更部は、当該判定結果に基づいて前記仮想磁場に関する値を変更する、
   請求項２に記載の電子状態算出装置。
4. 前記仮想磁場に関する値は、前記電子に印加された仮想磁場の大きさであり、
   前記変更部は、前記指標値が、前記所定範囲内の値ではないと判定した場合、前記仮想磁場の大きさを増大させる、
   請求項３に記載の電子状態算出装置。
5. 前記指標値は、前記電子のスピン状態を示す値である、
   請求項２から４のうちいずれか一項に記載の電子状態算出装置。
6. 前記電子のスピン状態を示す値は、上向きのスピンを有する前記電子の個数と、下向きのスピンを有する前記電子の個数との差分値である、
   請求項１から５のうちいずれか一項に記載の電子状態算出装置。
7. 前記物質とは、所定の原子であり、
   前記電子のスピン状態を示す値は、前記所定の原子に含まれる上向きのスピンを有する電子の個数と、前記所定の原子に含まれる下向きのスピンを有する電子の個数との差分値である、
   請求項１から６のうちいずれか一項に記載の電子状態算出装置。
8. 前記算出部は、コーン・シャム方程式に基づいて、前記スピン状態を示す値を算出する、
   請求項１から７のうちいずれか一項に記載の電子状態算出装置。
9. 物質に含まれる電子のスピン状態を示す値を、仮想磁場を示す値を利用して繰り返し算出し、
   前記スピン状態を示す値が所定値に収束するか否かを判定し、
   前記判定した結果に基づいて、前記仮想磁場に関する値を変更する、
   電子状態算出方法。
10. コンピュータに
    物質に含まれる電子のスピン状態を示す値を、仮想磁場を示す値を利用して繰り返し算出させ、
    前記スピン状態を示す値が所定値に収束するか否かを判定させ、
    前記判定させた結果に基づいて、前記仮想磁場に関する値を変更させる、
    電子状態算出プログラム。

Images