バンド計算を行うための手法は多数生み出されているが、その中でも最初期に登場したセルラー法（Wigner-Seitz法）に基づいてナトリウム金属の価電子帯の最低エネルギー準位を計算するコードをPythonで書いてみた。

#解くべき式
詳細については参考文献[1][2]を見てもらうことにして、解くべき式を原子単位系で書くと

\frac{d^2R(r)}{dr^2}=\left(-\frac{2Q(r)}{r^2}-E\right)R(r)\tag{1}

となる。ここで$R(r)=r\chi(r)$であり、求めたい波動関数が$\chi(r)$である。
$Q(r)$は文献[2]にあるように

\begin{eqnarray}
Q(r)=\left\{ \begin{array}{ll}
11r & (0\leq r<0.01) \\
-26.4 r^2+11.53r-0.00264 & (0.01\leq r <0.15) \\
-2.84r^2+4.46r+0.5274 & (0.15\leq r<1.00) \\
\vdots\\
\end{array} \right.
\end{eqnarray}

と非常に複雑な関数になっている(Naの3s軌道の電子が感じるポテンシャルになるように実験値から決定されたようだ)。

式(1)の微分方程式を4次のルンゲクッタ法で解けば$R(r)$が求まる。すなわち次の二つの微分方程式に分割して、これらを連立させて解く。

\frac{dR(r)}{dr} = p\\
\frac{dp}{dr} = \left(-\frac{2Q(r)}{r^2}-E\right)R(r) = F_1(r,R,p)

ただし、セルラー法の仮定から導かれる境界条件を課す必要がある。$r_0$をウィグナー・サイツ基本格子と同じ体積を持つ球体の半径であるとすると境界条件は以下のようになる。

\left.\frac{dR(r)}{dr}\right|_{r_0} = \frac{R(r_0)}{r_0} \tag{2}

この境界条件を満たす解を探すために論文[2]中では次のような戦略がとられた。

$E$に適当な値を代入して式(1)を解く。
得られた$R(r)$を$r$に対してプロットし、式の境界条件が満たされる$r_0$を作図によって探す。すなわち点($r_0$,$R(r_0)$)での$R(r)$の接線が原点を通るかどうかを判断すればよい。
1. 2.を繰り返して$r_0$とそのときのエネルギー$E$をプロットして曲線$E(r_0)$を得る。これは境界条件の導入によって、ある$r_0$に対してあるエネルギー固有値（そして波動関数）が決まることを示している。
金属Naを計算したい場合は$r_0=3.95$でのエネルギーを読み取り、再度式(1)に代入することで波動関数$\chi(r)$を得る。

ただし、2. のプロセスについては$r_i$での微分を差分に置き換えて、それとその点での傾きとの差の絶対値${\rm abs}$が小さい点を抽出するようにした。$r_i$は離散化した$i$番目の$r$の値である。

{\rm abs}=\left|\frac{R(r_{i+1})-R(r_{i-1})}{2h} - \frac{R(r_i)}{r_i}\right|\tag{3}

#プログラム

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from scipy.optimize import curve_fit

h = 0.002
hinv = 1 / h
xmin = h
xmax = 10.0
h2 = h * h
dh = h * 0.5
num = int((xmax - xmin)/h)

pot = []
def Make_pot(h):
    mw = 0
    w = int(0.01*hinv)
    for i in range(mw+1,w+1):
        pot.append(11.0*i*h)

    mw = w
    w = mw + int(0.14*hinv)
    for i in range(mw+1,w+1):
        pot.append(-26.4*i*i*h2+11.53*i*h-0.00264)

    mw = w
    w = mw + int(0.85*hinv)
    for i in range(mw+1,w+1):
        pot.append(-2.84*i*i*h2+4.46*i*h+0.5275)

    mw = w
    w = mw + int(0.55*hinv)
    for i in range(mw+1,w+1):
        pot.append(1.508*i*i*h2-4.236*i*h+4.876)

    mw = w
    w = mw + int(1.75*hinv)
    for i in range(mw+1,w+1):
        pot.append(0.1196*i*i*h2+0.2072*i*h+1.319)

    mw = w
    w = mw + int(3.44*hinv)
    for i in range(mw+1,w+1):
        pot.append(0.0005*i*i*h2+0.9933*i*h+0.0222)

    mw = w
    w = mw + int(3.26*hinv)
    for i in range(mw+1,w+1):
        pot.append(i*h)

    print(len(pot))


Make_pot(h)
pa = np.array(pot)
xa = np.arange(xmin,xmax+h,h)
v = -pa/xa/xa
print(v)
def Solv(e):
    def F1(x,y,p):
        nd = int(x*hinv)-1
        return (2.0*v[nd]-e)*y

    x = xmin
    y = x
    p = 1.0

    for i in range(num):
        dp1 = h * F1(x,y,p)
        dy1 = h * p

        dp2 = h * F1(x+dh,y+dy1*0.5,p+dp1*0.5)
        dy2 = h * (p+dp1*0.5)

        dp3 = h * F1(x+dh,y+dy2*0.5,p+dp2*0.5)
        dy3 = h * (p+dp2*0.5)

        dp4 = h * F1(x+h,y+dy3,p+dp3)
        dy4 = h * (p+dp3)

        xl.append(x)
        yl.append(y)

        x = x + h
        y = y + (dy1 + 2.0*dy2 + 2.0*dy3 + dy4)/6.0
        p = p + (dp1 + 2.0*dp2 + 2.0*dp3 + dp4)/6.0


results = []
for j in range(30,76):
    xl = []
    yl = []

    ex =  -j * 0.01
    Solv(ex)
    mdif1 = 50.0
    mdif2 = 51.0
    mem1 = 0.0
    mem2 = 0.0
    for k in range(200,num-1):
        dif = (yl[k+1]-yl[k-1])/h*0.5
        tes = abs(dif-yl[k]/xl[k])
        if tes < mdif1:
            mdif2 = mdif1
            mem2 = mem1
            mdif1 = tes
            mem1 = xl[k]

        elif tes < mdif2 and xl[k]-mem1 > 0.02:
            mdif2 = tes
            mem2 = xl[k]

    if abs(mem1-mem2) > 0.02:
        results.append([mem1,ex,mdif1])
        results.append([mem2,ex,mdif2])
    else:
        results.append([mem1,ex,mdif1])

$h$は微分方程式を解く際の離散化の幅である。

プログラムのはじめの部分ではポテンシャル関数を生成している（Make_pot）。
つづいてあるエネルギーの値に対して4次のルンゲクッタ法によって式(1)を解く関数Solvが定義されている。

$r=r_{\rm min}(=h)$から解き始めるが、初期値については$R(r)$が原点近傍で

R(r) \sim r

となることから求めた（参考）。

その後の$j$に関するループではエネルギーを-0.3から-0.75の間で0.01ずつ変えながら式(1)を解いている。続く処理では式(3)で評価される絶対値が小さい点を抽出している。ただし、$r$の値が近い2点が選ばれないようにif文による条件分岐がある。

結果はresultsというリストに

[r_0,\ E,\ {\rm abs}]

の形で格納される。

#解析
まず$E$と$r_0$のグラフを見てみる。きれいな曲線関係が見て取れる。

rl = np.array(results)
sc = plt.scatter(rl[:,0],rl[:,1],c=rl[:,2],norm=matplotlib.colors.LogNorm())
plt.colorbar(sc)
plt.show()

この図では$\log({\rm abs})$の大きさによって点の色分けがなされている。$E\leq-0.69$の点については絶対値がほかの点に比べて大きいため境界条件が満たされていないと判断し、以下の解析では除くことにする。

ちなみに、いくつかの$E$について$R(r)$を$r$に対してプロットした図が下図である。上の図と比較してみると境界条件(2)が満たされているところでの$r$が$r_0$となっていることがわかる。

続いて得られた点に対して曲線あてはめを行う。いろいろな関数式を試した結果

f(r\ ;\ e,r_0,n_1,n_2,n_3,a,b) = e\left\{ \left( \frac{r_0}{r} \right)^{n_1} - a \left( \frac{r_0}{r} \right)^{n_2} - b\left(\frac{r_0}{r} \right)^{n_3} \right\}\tag{4}

で比較的よくフィットできることがわかった。ただし変数を$r_0$から$r$に変えた。以下はカーブフィッティングのコード。

rl2 = rl[rl[:,1]>-0.69]
fitx = rl2[:,0]
fity = rl2[:,1]

def Fit_func(x,ep,x0,n1,n2,n3,a,b):
    return ep*((x0/x)**n1 - a*(x0/x)**n2 - b*(x0/x)**n3)

popt, pcov = curve_fit(Fit_func,fitx,fity,maxfev=10000)

plt.scatter(fitx,fity,s=15,zorder=2)

rr = plt.axis()
ra = np.arange(rr[0],rr[1],0.01)
plt.ylim(rr[2],rr[3])

plt.xlabel(r"$r$")
plt.ylabel(r"$E$")

plt.plot(ra,Fit_func(ra,*popt),c="orange",zorder=1,label="f(r)")
plt.legend()
plt.show()

フィッティングの結果を示す。点が非常によく曲線上にのっていることがわかる。

最後に$r=3.95$について、求められたパラメータのもとで式(4)を計算すると$E=-0.6175$となった。

これを再び式(1)に代入して波動関数$\chi(r)$を計算する。

r0 = 3.95
xl = []
yl = []
ex = Fit_func(r0,*popt)
print(ex)
Solv(-0.6175)
xa = np.array(xl)
ya = np.array(yl)

plt.plot(xa,ya/xa)
plt.xlim(0,r0)
plt.ylim(-0.2,0.2)
plt.hlines(0,0,r0,"black",linewidth=0.5)
plt.xlabel(r"$r$")
plt.ylabel(r"$\chi(r)$")

plt.show()

規格化をしていないので絶対値は合わないが、[1]や[3]のグラフと非常によく似たものが得られた。

#最後に
文献[2]に

The calculation of a wave function took about two afternoons, and five wave functions were calculated on the whole, giving the ten points of the figure.

という一節がある。コンピュータなしでこの問題を解くのはさぞ骨の折れる作業だっただろう。それにも関わらず正しく結果が求められていてすごいなと思った。

現代ではこのような手計算をする必要はないにせよ、自分の計算力が圧倒的に不足しているような気がしてならない、、、

#参考文献
[1] アシュクロフト・マーミン，固体物理の基礎上・IIー固体のバンド構造ー
[2] E. Wigner and F. Seitz, Phys. Rev. 43:804, 1933.
[3] E. Wigner and F. Seitz, Phys. Rev. 46:509, 1934.

また、ルンゲクッタ法の実装については
https://nemuneko-gensokyou.blog.so-net.ne.jp/2015-04-14
を参考にした。