水素原子に対するSchrödinger方程式の数値解法

水素原子に対するSchrödinger方程式 の数値解法 @dc1394

概要       Schrödinger方程式と用いる単位系について 3次元のSchrödinger方程式（偏微分方程式）を変 数分離し，1次元の常微分方程式に還元 境界値の推定 常微分方程式の数値解法 Brentアルゴリズムを用いた固有値の探索 波動関数の構成と規格化（正規化）

Hartree原子単位系     このスライドでは，Schrödinger方程式の表式を簡 潔にするために，Hartree原子単位系を使用する． この単位系では，長さの単位はBohr半径a0 (1 [a0] = 5.29×10-11 [m]), 質量の単位は電子の質量me, 電荷は電気素量e, エネルギーはHartree (1 [Hartree] = 4.36×10-18 [J] = 27.2 [eV])を用いる． この単位系では，Dirac定数ℏと，Coulombポテン シャルの比例定数1 / (4πε0)が1となる． 単位を表す記号として，atomic unit の省略形であ る a.u. で表す．

水素原子のSchrödinger方程式  最も簡単な水素原子について，定常状態におけ るSchrödinger方程式（以下Sch方程式と書く）を以 下に示す． 電子の運動エネルギーポテンシャル    ここで， Coulombポテンシャル である． この方程式は（少なくとも見かけ上は）単純であり， また解析的に解くことができる（ただし，非常に面 倒である）． 今回は，この方程式を数値的に解くことを考える．

水素原子のSchrödinger方程式の数 値的解法の種類  純粋な数値的解法  修正狙い撃ち法（本スライドで解説） ◼  特殊な常微分方程式の二点境界値問題に帰着 行列を用いた解法（数値線形代数による解法） 基底で展開する方法  有限差分法（cf. https://qiita.com/cometscome_phys/items/f2dc91364f8a8 3020495 ）  有限要素法（cf. https://qiita.com/dc1394/items/c0d3d02fa738b040128c ）  ◼ いずれも（一般化）固有値問題に帰着

素朴な疑問  なぜ水素原子のSchrödinger方程式は解析的に解 けるのに，数値的に解こうとするのか？  多電子原子のSchrödinger方程式（と言うよりHartree- Fock方程式やKohn-Sham方程式）を解くための練習

Sch方程式の変数分離  水素原子に対するSch方程式を，以下のように書 く．  ここで，極座標におけるラプラシアンΔは，  であり，V(r)は球対称であるので，  と変数分離が可能である．ここで，n, l, mはそれぞ れ主量子数，方位量子数，磁気量子数である．

素朴な疑問2  変数分離せずに直接，3次元直交座標または3次 元球座標で解くことはできないか？  一応可能．ただし，水素原子のみならず，たいていの 原子で，ポテンシャルが球対称の条件を使えるので， 系の有する対称性を用いた方が良い．

Sch方程式の変数分離    この変数分離により，以下の二つの微分方程式 が得られる． なお，この計算は， https://github.com/dc1394/fukui_sono2/blob/mas ter/separation_of_variables.pdf を参照のこと． 第二式の解は，球面調和関数として解析的に得 られる．従って，数値的に解くべき方程式は第一 式である．

Sch方程式の数値解法上の困難     ここで，この常微分方程式の境界条件は， Lnl(0)=有限, Lnl(∞)=0である． しかし，これらの境界条件は，この常微分方程式 を数値的に解く上で，何も言っていないのと同じ である． さらに，この常微分方程式は固有方程式であり，E も未知数であるが，Eが固有値以外の値では解が 発散する．

Sch方程式を二点境界値問題にする     原点は確定特異点であるので，原点から微分方 程式を数値的に解くことができない． また，コンピュータでは無限大を扱えないので， 無限遠点から微分方程式を数値的に解くことも不 可能である． 従って，原点に十分近い点と，原点から十分離れ た点で， Lnl(r)とその微分dLnl(r)/drの値を調べ，そ の二つの点から微分方程式を数値的に解く． 最終的に，二点境界値問題に帰着する．

Sch方程式の固有値Eの求め方    適当なあるEを用いて，原点に十分近い点および 原点から十分離れた点から，常微分方程式を数 値的に解くと，Eが固有値でない限り，二つの常微 分方程式の解は一致しない． 逆に，もしEが固有値であれば，原点に十分近い 点と原点から十分離れた点の間の点（マッチング ポイント）で，二つの常微分方程式の解は，一致 はしないもののある関係式を満たす． 従って，これにより固有値Eを見つけることができ る．

対数メッシュの導入   ポテンシャルV(r)は原点付近で大きく変化する． このような空間の異方性を考慮するために，対数 メッシュr = exp(x)を導入（変数変換）して微分方 程式を変形する．

変数変換した微分方程式   この変数変換により， Lnl(r)の一階微分と二階微分 は以下の式となる． これらを目的の微分方程式に代入することによっ て，次式が得られる．

連立一階常微分方程式に書き直す     二階の常微分方程式を直接，数値的に解くことは 難しい． 従って，二階の常微分方程式を，二元連立一階 常微分方程式に書き直す．すると，次の二つの式 が得られる． これらの式が，数値的に解くべき方程式である． 次に，Lnl(x)とM(x)の境界値について考える．

V(r)とLnl(r)を級数展開する    まずは原点に十分近い点で，関数V(r)と関数Lnl(r) がどう振る舞うか調べてみる（Frobeniusの方法を 用いる）． 関数V(r)と関数Lnl(r)は，以下のように級数展開で きるはずである． 従って， Lnl(r)の一階微分と二階微分は以下となる．

級数展開した式を代入する  前ページの結果を上式に代入すると， 左辺=  右辺=  が得られる． 

Lnl(r)の級数展開の係数を求める   左辺と右辺のrのべき乗の係数を比較することによっ て， Lnl(r)の級数展開の係数b0～b4は以下のように得 られる． ここで，V(r)の級数展開の係数a0～a2（とE）は，未だ未 知数であるので，別に求める必要がある．

原点付近のLnl(x)とM(x)の近似値    V(r)の級数展開の係数a0～a2は，以下の連立一 次方程式から求められる． この連立一次方程式の解は，連立一次方程式の ソルバーを用いれば，簡単に得られる． その解を用いて，以下のように，原点付近(x = x0) でのLnl(x0)とM(x0)の近似値が求められる．

原点から十分離れた点での波動関数 の振る舞い    波動関数を以下のように書くと， 次式が成り立つ（この計算は， https://github.com/dc1394/fukui_sono2/blob/master/ separation_of_variables_2.pdf を参照のこと） 左辺中括弧の中の第二項と第三項の寄与は，原点か ら十分離れたところでは無視できる．従って，  となる．この微分方程式の解析解は容易に分かり，  である．

原点から十分離れた点でのLnl(x)と M(x)の近似値  ここで，  である．  従って，上式に前ページの結果を代入すると，原 点から十分離れた点(x = xmax)でのLnl(xmax)と M(xmax)の近似値は次式となる．

常微分方程式の解法    これらの近似値を初期値として用いて，あるエネ ルギーEを仮定し，原点に十分近い点と，原点か ら十分離れた点から，常微分方程式を数値的に 解いていく． 常微分方程式を数値的に解くには，数値計算ライ ブラリのソルバーを利用すると良い（もちろんソル バーを自作しても良い）． ソルバーのアルゴリズムは，Adams-BashforthMoulton法（予測子修正子法），Burilrsh-Stoer法（補 外法），Controlled Runge-Kutta法（誤差とステップ 幅を制御したRunge-Kutta法）のいずれかを使うと 良い．

Lnl(x)とM(x)のグラフ  l = 0, E = -0.5 (Hartree)に対して，x = -8.0（緑線） 及びx = 6.0（青線）から解くと，以下のようになる．

マッチングポイント マッチングポイント  二つのLnl(x)あるいはM(x)を比較する点を，マッチ ングポイントと呼ぶ．

関数ΔD(E)を定義する    もし選んだEが固有値であるならば，次式が成り立 つ． ここで，以下の関数ΔD(E)を定義する． この関数ΔD(E)がゼロになるE（つまり非線形方程 式ΔD(E) = 0の根）が，固有値である．

固有値Eを探索する   ΔD(E)は固有値の前後で符号が変化する． 従って，固有値を探索するアルゴリズムは以下の ようになる．  (1) Eをスキャンし，ΔD(E)の符号が変化する領域を探 す．  (2) 符号が変化する領域が見つかったら，その領域内 で，Brent法を用いて，ΔD(E) = 0の根を求める．なお， Brent法は，非線形方程式の根を非常に効率的に求 めるアルゴリズムである．  (3) 求まった根が目的の固有値Eである．

波動関数の構成と正規化  見つかった固有値Eに対して，以下の順序で波動 関数を求める．  (1) Lnl,O(x)とLnl,I(x)を接合してLnl(x)を構成する．  (2) 以下の式により，Lnl(x)を規格化（正規化）する．

波動関数のノード数   それぞれの波動関数のノード数はn – l – 1となる． 重複して解を求めてしまわないように，ノード数が 適切なものになっているか調べる必要がある． ノード（節）

最終的に得られる波動関数    最終的に得られる波動関数は（変数をxからrに戻 した）， 及び，それにrを乗じた形 である．これらをrのメッシュの値と共にファイルに 出力する．

厳密な固有値と計算で求めた固有値 の比較    水素原子の場合，厳密な固有値は解析的に求め られ， となる．上式より，1s軌道の場合，厳密な固有値 は-0.5 (Hartree)である． 著者のコードでの計算値は（インプットファイルで 与えるパラメータにもよるが），-0.49999998 Hartree)であり，非常に良く一致している．

厳密なRnl(r)と，計算で求めたRnl(r)を 比較したプロット（H原子の1s軌道）

厳密なRnl(r)と，計算で求めたRnl(r)を 比較したプロット（y軸対数目盛）

厳密なPnl(r)と，計算で求めたPnl(r)を 比較したプロット（H原子の1s軌道）

厳密なPnl(r)と，計算で求めたPnl(r)を 比較したプロット（y軸対数目盛）

実行画面

ソースコードへのリンク    このプログラムのソースコードは，GitHub上で公 開しています https://github.com/dc1394/schrac ライセンスは修正BSDライセンスとします．

schracで計算した水素原子の5dxy軌 道の3次元プロット  プロットには， SchracVisualize_ Direct3D_11とい う自作ソフトを 使用 （ https://github. com/dc1394/Sc hracVisualize_Di rect3D_11/rele ases/tag/v0.2 からダウンロー ド可能）

まとめ     Sch方程式を，原点付近と原点から十分遠い点か ら数値的に解いた． それぞれの解を接合することにより，固有値及び 波動関数を得た． 計算によって得られた波動関数から，運動エネル ギー及びポテンシャルエネルギーを計算した． 計算で求めた固有値及び波動関数のいずれも， 解析的に求められる値とほとんど完全に一致して いた．