非相対論的な水素様原子の電子の波動関数を求めるfortranコードを書きます。

解く問題

電荷\(+Ze\)の原子核の周りに存在する、電荷\(-e\)を持つ電子が満たす時間依存しないシュレーディンガー方程式は、
\(
\displaystyle \left[-\frac{1}{2}\Delta-\frac{Z}{r}\right]\psi({\bf r})=E\psi({\bf r})
\)
です。
\(
\left. \psi({\bf r})\right|_{|{\bf r}|\to \infty}\to 0
\)
の境界条件を満たす正規直交化された解は、
\(
\begin{align}
\psi({\bf r})&=R_{n,l}(r)\cdot Y_{l,m}(\theta,\varphi)\\
&=\sqrt{\left(\frac{2}{n}\right)^3\frac{(n-l-1)!}{2n[(n+l)!]}} e^{-Zr/n} \cdot \left(2\frac{Zr}{n}\right)^l\cdot L_{n-l-1}^{(2l+1)}\left(2\frac{Zr}{n}\right)\cdot Y_{l,m}(\theta,\varphi)
\end{align}
\)
と表されます。ここで、

• \(n,l,m\)  主量子数,方位量子数,磁気量子数
• \(r,\theta,\varphi\)  動径方向,z軸方向の角度,xy平面の角度
• \(L_n^{(\alpha)}(x)\)  ラゲール陪多項式[1]
• \(Y_{l,m}(\theta,\varphi)\)  球面調和関数[2]

を意味します。

動径関数\(R_{n,l}(r)\)は微分方程式
\(
\displaystyle \left[-\frac{1}{2}\left(\frac{d^2}{dr^2}+\frac{2}{r}\frac{d}{dr}\right)+\frac{l(l+1)}{2r^2}-\frac{Z}{r}\right]R_{n,l}(r)=E_n R_{n,l}(r)
\)
を満たす正規直交化された関数で、固有値は
\(
\displaystyle E_n=-\frac{1}{2n^2}
\)
です。

波動関数を出力するFortranコードはこれ。

計算結果

3次元の波動関数を表示するには3つの座標と波動関数の値の計4つが必要です。
しかし、4次元空間を分かりやすく描画する方法はありません。

そこで、ここでは\(\mathrm{Re}{\psi(x,0,z,t)}\)のように、\(y=0\)にしてグラフを描画します。
\(z\)軸周りの構造がうまく表現できませんが、良しとします。

\(n=1\)

\(n=2\)

\(n=3\)

”color”と書いてあるのは、波動関数の形がはっきり映るようにカラーバーの範囲を調整したものです。
なので、色が飽和しているところはそのカラーバーの範囲で表現することはできない領域です。

\(n=20\)まで同じような画像を作ったのですが、画像の容量が大きくなりすぎました。
なので、見たい方は以下のリンクから元ファイルをダウンロードしてください。

\(n=1\)  Hatomc_n1.png（48KB）

\(n=2\)  Hatomc_n2.png（223KB）

\(n=3\)  Hatomc_n3.png（612KB）

\(n=4\)  Hatomc_n4.png（1MB）

\(n=5\)  Hatomc_n5.png（2MB）

\(n=6\)  Hatomc_n6.png（3MB）

\(n=7\)  Hatomc_n7.png（4MB）

\(n=8\)  Hatomc_n8.png（3MB）

\(n=9\)  Hatomc_n9.png（6MB）

\(n=10\)  Hatomc_n10.png（6MB）

\(n=11\)  Hatomc_n11.png（9MB）

\(n=12\)  Hatomc_n12.png（9MB）

\(n=13\)  Hatomc_n13.png（10MB）

\(n=14\)  Hatomc_n14.png（14MB）

\(n=15\)  Hatomc_n15.png（15MB）

\(n=16\)  Hatomc_n16.png（16MB）

\(n=17\)  Hatomc_n17.png（18MB）

\(n=18\)  Hatomc_n18.png（19MB）

\(n=19\)  Hatomc_n19.png（21MB）

\(n=20\)  Hatomc_n20.png（21MB）

n=1～n=20までのファイル hatom.zip（164MB）

例えば、n=13, n=20はこんな感じです（荒くしてます）。

特殊関数に関するプログラミングは特殊関数へ。

How to plot in gnuplot?
Just for note for myself.

You need 2 files which named : multi.plt and atommulti.plt

Before ploting, you must make directory “dat” and move every data files to “dat”.

and in gnuplot, type

then, you can get figure.

multi.plt

atommulti.plt

参考

[1]ラゲール陪多項式\(L_m^{(\alpha)}(x)\)

• 定義域

  \(\;\;\;\;
  0 \lt x \lt \infty
  \)

• 微分方程式

  \(\;\;\;\;
  \displaystyle
  \left[x\frac{d^2}{dx^2}+(\alpha+1-x)\frac{d}{dx}+n\right]L_m^{(\alpha)}(x)=0
  \)

• 具体例

  \(\;\;\;\;
  \begin{align}
  & L^{(0)}_0(x)=1 \\
  & L^{(0)}_1(x)=-x+1 \\
  & L^{(0)}_2(x)=\frac{1}{2}(x^2-4x+2) \\
  & L^{(0)}_3(x)=\frac{1}{6}(-x^3+9x^2-18x+6) \\
  & L^{(1)}_0(x)=1 \\
  & L^{(1)}_1(x)=-x+2 \\
  & L^{(1)}_2(x)=\frac{1}{2}\left(x^2-6x+6\right) \\
  & L^{(1)}_3(x)=\frac{1}{6}\left(-x^3+12x^2-36x+24\right)
  \end{align}
  \)

• 漸化式

  \(\;\;\;\;
  \begin{align}
  L_n^{(\alpha)}(x)&=\binom{n+\alpha}{n}a_0(x) \\
  &\hspace{1em}a_{m-1}(x)=1-\frac{n-m+1}{m(\alpha+m)}\cdot x\cdot a_m(x) \\
  &\hspace{5em}(m=n,n-1,n-2,\cdots,1,\;\;a_n(x)=1)
  \end{align}
  \)
  [3]
  1階微分
  \(\;\;\;\;
  \begin{align}
  \frac{d}{dx}L_n^{(\alpha)}(x)=-L_{n-1}^{(\alpha+1)}(x)
  \end{align}
  \)
  2階微分
  \(\;\;\;\;
  \begin{align}
  \frac{d^2}{dx^2}L_n^{(\alpha)}(x)=L_{n-2}^{(\alpha+2)}(x)
  \end{align}
  \)

• 直交性

  \(\;\;\;\;
  \displaystyle
  \int_{0}^{\infty} e^{-x}x^{\alpha}L^{(\alpha)}_n(x)L^{(\alpha)}_{n’}(x) dx =\frac{\Gamma(\alpha+n+1)}{n!}\delta_{n,n’},\;\;(\alpha>-1)
  \)

[2]球面調和関数\(Y_{l,m}(\theta,\varphi)\)
ルジャンドル多項式(\(P_l^m(x)\),[3])を用いて
\(
Y_{l,m}(\theta,\varphi)=(-1)^{(m+|m|)/2}\sqrt{\frac{2l+1}{4\pi}\frac{(l-|m|)!}{(l+|m|)!}}P_l^{|m|}(\cos\theta)e^{im\phi}
\)

[3]ルジャンドル陪多項式\(P_l^m(x)\)

• 定義域

  \(\;\;\;\;
  1 \lt x \lt -1
  \)

• 微分方程式

  \(\;\;\;\;
  \displaystyle (1-x^2)\frac{d^2}{dx^2}P_l^m(x)-2x\frac{d}{dx}P_l^m(x)+\left[l(l+1)-\frac{m^2}{1-x^2}\right]P_l^m(x)=0
  \)

• 具体例

  \(\;\;\;\;
  \begin{align}
  & P_0^0(x)=1 \\
  & P_{1}^{-1}(x)=\frac{1}{2}(1-x^2)^{1/2} \\
  & P_{1}^{0}(x)=x \\
  & P_{1}^{1}(x)=-(1-x^2)^{1/2} \\
  & P_{2}^{-2}(x)=\frac{1}{8}(1-x^2) \\
  & P_{2}^{-1}(x)=\frac{1}{2}x(1-x^2)^{1/2} \\
  & P_{2}^{0}(x)=\frac{1}{2}(3x^2-1) \\
  & P_{2}^{1}(x)=-3(1-x^2)^{1/2} \\
  & P_{2}^{2}(x)=3(1-x^2)
  \end{align}
  \)

• 漸化式

  用いている式
  \(\;\;\;\;
  \begin{align}
  & P_{|m|}^{|m|}(x)=(-1)^{|m|}(2|m|-1)!!(1-x^2)^{\frac{|m|}{2}} \\
  & P_{l}^{-m}(x)=(-1)^m\frac{(l-m)!}{(l+m)!}P_l^m(x) \\
  & P_{|m|}^{|m|+1}(x)=(2|m|+1)xP_{|m|}^{|m|}(x) \\
  & P_{|m|}^{|m|+q}(x)=\left(\frac{2|m|-1}{q}+2\right)xP_{|m|}^{|m|+q-1}(x)-\left(\frac{2|m|-1}{q}+1\right)P_{|m|}^{|m|+q-2}(x)
  \end{align}
  \)[1]
  1階微分
  \(\;\;\;\;
  \begin{align}
  m=\pm 1&\; \\
  &:\hspace{1em} \frac{d}{dx}P_l^m(x)=\frac{lx}{(x^2-1)}P_l^m(x)-\frac{l+m}{(x^2-1)} P_{l-1}^m(x) \\
  m\ne\pm 1\; \\
  &:\hspace{1em} \frac{d}{dx}P_l^m(x)=c_{2}P^{m+2}_{l-1}(x)+c_0 P^{m}_{l-1}(x)+c_{-2} P^{m-2}_{l-1}(x) \\
  &\hspace{1em} c_{2}=\frac{1}{4(m+1)}\\
  &\hspace{1em} c_{0}=\frac{l+m}{2}\left(1+\frac{l}{1-m^2}\right)\\
  &\hspace{1em} c_{-2}=-\frac{(l+m)(l+m-1)(l+m-2)(l-m+1)}{4(m-1)}
  \end{align}
  \)
       ※\(m=\pm 1\)の時、\(x=\pm 1\)で発散します。
  2階微分
  \(\;\;\;\;
  \begin{align}
  m=\pm 1,\pm 3&\; \\
  &:\hspace{1em} \frac{d^2}{dx^2}P_l^m(x)=\frac{1}{1-x^2}\left[2x\frac{dP_l^m(x)}{dx}+\frac{(l+1)(l+m)}{2l+1}\frac{dP_{l-1}^m(x)}{dx}-\frac{l(l-m+1)}{2l+1}\frac{dP_{l+1}^m(x)}{dx}\right] \\
  m\ne\pm 1,\pm 3&\; \\
  &:\hspace{1em} \frac{d^2}{dx^2}P_l^m(x)=c_{2}\frac{dP^{m+2}_{l-1}(x)}{dx}+c_0 \frac{dP^{m}_{l-1}(x)}{dx}+c_{-2} \frac{dP^{m-2}_{l-1}(x)}{dx}
  \end{align}
  \)
       ※\(m=\pm 1,\pm 3\)の時、\(x=\pm 1\)で発散します。

• 直交性

  \(l\)に対する直交性
  \(\;\;\;\;
  \displaystyle
  \int_{-1}^{1} P_{m}^{l}(x)P_{m}^{l’}(x) dx =\frac{2}{2l+1}\cdot \frac{(l+m)!}{(l-m)!}\delta_{l,l’}
  \)

[3] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 22.18, Page 789)
[4] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 22.2, Page 774,775)