10.3K Views
April 18, 21
スライド概要
水素原子に対する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方程式を,原点付近と原点から十分遠い点か ら数値的に解いた. それぞれの解を接合することにより,固有値及び 波動関数を得た. 計算によって得られた波動関数から,運動エネル ギー及びポテンシャルエネルギーを計算した. 計算で求めた固有値及び波動関数のいずれも, 解析的に求められる値とほとんど完全に一致して いた.