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

スライド概要

水素原子に対するSchrödinger方程式の数値解法について解説したスライドです。

profile-image

dc1394

@dc1394

作者について:

主にC++,C#,Juliaを使っています。数値解析が趣味。 博士後期課程中退。専門は第一原理計算。

スライド一覧
シェア
埋め込む

作成日

2021-04-18 09:54:51

各ページのテキスト

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

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

3. 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. で表す.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

35. 実行画面

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

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

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