この章で選ぶ $\hat H$ と当てる道具：Kohn–Sham の1電子のハミルトニアンと自己無撞着ループ

ハミルトニアンの選び方で固体を考える　第2部 第10章 §2.10.1　★★★ 全電子の $\hat H$ を、密度に依存する有効ポテンシャルの中の1電子の $\hat H_{\mathrm{KS}}$ に書き換え、密度とポテンシャルを互いに合うまで回す

第2部のここまでの章では、周期ポテンシャルの中の1電子の問題を、tight-binding の飛び移り $t$ やほぼ自由な電子の Fourier 成分 $V_{\bm{G}}$ という、外から与えたパラメータで解いてきた。第2部第9章の終わりには Si と GaAs のバンド図を読んだが、その図は多くの軌道を残した計算から借りたもので、どう計算されたかは述べなかった。この章ではその計算の手順を扱う。実在の物質のバンド図の大半は、密度汎関数理論（DFT）の Kohn–Sham 方程式を平面波の基底で解いて作られている。

DFT の定理の側、すなわち基底状態の密度が外部ポテンシャルを決めるという Hohenberg–Kohn の定理、密度についての変分原理を作り直す Levy–Lieb の制約付き探索、相互作用のない電子の問題へ書き換える Kohn–Sham の構成は、近似の数理 第9章の密度汎関数理論の節で証明とともに扱った。この章ではそれらを結果として引き、導き直さない。そのうえで、計算を実際に回すときに決めなければならない4つのこと、すなわち自己無撞着ループの回し方、基底（平面波）と擬ポテンシャルとカットオフと $\bm{k}$ 点、交換相関汎関数の選び方、そして得られた固有値をバンドとして読むときの注意（バンドギャップ問題）を順に述べる。

2.10.1節で Kohn–Sham の1電子のハミルトニアンを選び、自己無撞着ループと混合の安定性を扱う。1次元の小さな模型で、ループを実際に数値で回す。2.10.2節で平面波の基底と擬ポテンシャルを導入し、カットオフと $\bm{k}$ 点の数による収束を確かめる。2.10.3節で一様電子気体から局所密度近似（LDA）を作り、勾配近似（GGA）と混成汎関数へ進む。2.10.4節で、Kohn–Sham の固有値が準粒子のエネルギーでないことから、ギャップが実験より小さく出る理由を述べる。2.10.5節で何を削ったかを書き出して近似の裏付けの札を置き、2.10.6節で実験の量と授業で習った形に対応させ、2.10.7節でまとめる。

前の章から使う結果は次の4つである。第2部第4章で、周期ポテンシャルの中の1電子の固有状態が Bloch の定理で波数 $\bm{k}$ をラベルに持つことを示した。第2部第5章で、ほぼ自由な電子の Brillouin 域の境界に幅 $2\lvert V_{\bm{G}}\rvert$ のギャップが開くことを見た。第2部第7章で、ギャップ $E_{\mathrm g}$ の半導体の真性キャリア密度が $n_{\mathrm i}\propto\mathrm{e}^{-E_{\mathrm g}/(2 k_{\mathrm B} T)}$ となることを出した。第5部第1章で、一様電子気体の交換エネルギーと Thomas–Fermi の遮蔽を求めた。この章は第2部の中に置くが、第5部第1章の結果を先回りして借りる。

出発点の全体のハミルトニアン

第1部第1章で書いた固体の全体のハミルトニアンから出発する。原子核の位置 $\bm{R}_I$ は止め（Born–Oppenheimer 近似）、電子 $N$ 個だけを量子力学で扱う。

$$\hat H=\sum_{i=1}^{N}\Bigl[-\frac{\hbar^2}{2m}\nabla_i^2+v_{\mathrm{ext}}(\bm{r}_i)\Bigr]+\frac12\sum_{i\ne j}\frac{e^2}{4\pi\epsilon_0\lvert \bm{r}_i-\bm{r}_j\rvert },\qquad v_{\mathrm{ext}}(\bm{r})=-\sum_I\frac{Z_Ie^2}{4\pi\epsilon_0\lvert \bm{r}-\bm{R}_I\rvert }\tag{2.10.1}$$

ここで $v_{\mathrm{ext}}$ は電子1個が原子核から受ける位置エネルギーで、この章では「ポテンシャル」をすべてエネルギーの単位で書く。第2部のこれまでの章は、この式の電子間の反発を、何らかの平均の周期ポテンシャルに置き換えたものを出発点にしていた。この章で問うのは、その平均のポテンシャルを、模型のパラメータとしてではなく式 (2.10.1) から決める方法である。

Kohn–Sham の1電子のハミルトニアン

近似の数理 第9章の結果を、この巻の記号（SI）に書き直して引く。基底状態のエネルギーは電子密度 $n(\bm{r})$ の汎関数として

$$E[n]=T_{\mathrm s}[n]+\int v_{\mathrm{ext}}(\bm{r})n(\bm{r})\,\mathrm{d}^3r+E_{\mathrm H}[n]+E_{\mathrm{xc}}[n],\qquad E_{\mathrm H}[n]=\frac12\iint\frac{e^2\,n(\bm{r})n(\bm{r}')}{4\pi\epsilon_0\lvert \bm{r}-\bm{r}'\rvert }\,\mathrm{d}^3r\,\mathrm{d}^3r'\tag{2.10.2}$$

と書け、正しい密度でこれが最小になる。$T_{\mathrm s}$ は同じ密度を持つ相互作用のない電子の運動エネルギー、$E_{\mathrm H}$ は密度を古典的な電荷の雲とみなした静電エネルギー（Hartree エネルギー）で、残りをすべてまとめたものが交換相関エネルギー $E_{\mathrm{xc}}$ である。$E_{\mathrm{xc}}$ はこの式で定義されるので、ここまでは近似を含まない。$n=\sum_i\lvert \psi_i\rvert ^2$ と軌道で書いて変分をとると、次の1電子の固有値問題が得られる。

$$\hat H_{\mathrm{KS}}\psi_{n\bm{k}}=\varepsilon _{n\bm{k}}\psi_{n\bm{k}},\qquad \hat H_{\mathrm{KS}}=-\frac{\hbar^2}{2m}\nabla^2+v_{\mathrm{eff}}(\bm{r}),\qquad v_{\mathrm{eff}}=v_{\mathrm{ext}}+v_{\mathrm H}+v_{\mathrm{xc}}\tag{2.10.3}$$

$$v_{\mathrm H}(\bm{r})=\int\frac{e^2\,n(\bm{r}')}{4\pi\epsilon_0\lvert \bm{r}-\bm{r}'\rvert }\,\mathrm{d}^3r',\qquad v_{\mathrm{xc}}(\bm{r})=\frac{\delta E_{\mathrm{xc}}}{\delta n(\bm{r})},\qquad n(\bm{r})=\sum_{n\bm{k}\sigma}f_{n\bm{k}}\lvert \psi_{n\bm{k}}(\bm{r})\rvert ^2\tag{2.10.4}$$

$f_{n\bm{k}}$ は占有数（絶縁体では満ちたバンドで 1、空のバンドで 0）である。結晶では $v_{\mathrm{eff}}$ が格子の周期を持つので、第2部第4章の Bloch の定理がそのまま当たり、固有状態は $n$ と $\bm{k}$ でラベルされる。これが、この章で選ぶ $\hat H$ である。見た目は第2部の1電子の問題と同じだが、ポテンシャル $v_{\mathrm{eff}}$ が解くべき軌道から作った密度 $n$ に依存する点だけが違う。

この書き換えで何が得られ、何が残ったかを確かめておく。得られたものは、$3N$ 変数の多体の波動関数の代わりに、3変数の1電子の方程式を $N$ 本解けばよくなったことである。運動エネルギーの大部分 $T_{\mathrm s}$ を軌道で正確に扱うので、原子の殻構造や化学結合が正しく出る。残ったものは、$E_{\mathrm{xc}}$ の形が分からないことだけである。全エネルギーの中で $E_{\mathrm{xc}}$ は運動エネルギーや静電エネルギーより1桁ほど小さく、その小さな部分だけを近似すればよい、というのが Kohn と Sham の書き換えの要点である（W. Kohn and L. J. Sham (1965), Phys. Rev. 140, A1133）。

当てる道具：自己無撞着ループ

式 (2.10.3) は、$v_{\mathrm{eff}}$ が分からないと解けず、解けないと $n$ が分からず、$n$ が分からないと $v_{\mathrm{eff}}$ が作れない。そこで、入力のポテンシャルを仮に決めて解き、出てきた密度からポテンシャルを作り直し、入力と出力が一致するまで繰り返す（図 2.10.1）。これを自己無撞着場（SCF）の計算と呼ぶ。

出力をそのまま次の入力にする（$\alpha=1$）と、ループはたいてい収束しない。理由は線形化して見ると分かる。正しい解 $v^*$ からのずれ $\delta v_{\mathrm{in}}$ が小さいとき、出力のずれはそれに比例し、$\delta v_{\mathrm{out}}=\hat{\mathcal J}\,\delta v_{\mathrm{in}}$ と書ける。$\hat{\mathcal J}$ は「ポテンシャルを少し変えたとき、電子が動いて作り直されるポテンシャルがどれだけ変わるか」を表す線形の写像である。電子どうしは反発するので、ポテンシャルの谷を深くすると電子が集まり、その電子の作る Hartree ポテンシャルは谷を浅くする向きに加わる。このため $\hat{\mathcal J}$ の固有値 $\lambda$ は負になる。混合を入れると、次のずれは

$$\delta v_{\mathrm{in}}^{\mathrm{new}}=\bigl[1-\alpha(1-\hat{\mathcal J})\bigr]\delta v_{\mathrm{in}}\tag{2.10.5}$$

となる。固有値 $\lambda$ の向きのずれは毎回 $1-\alpha(1-\lambda)$ 倍される。$\lambda<0$ なので $1-\lambda>1$ で、$\alpha=1$ のときこの倍率は $\lambda$ そのものになる。$\lvert \lambda\rvert >1$ の向きがあれば $\alpha=1$ では発散し、収束させるにはすべての固有値について $\lvert 1-\alpha(1-\lambda)\rvert <1$、すなわち $0<\alpha<2/(1-\lambda_{\min})$ が要る。$1-\hat{\mathcal J}$ は、外から入れたポテンシャルが電子の応答でどれだけ弱められるかを表すので、固体の誘電関数（第5部第1章の $\varepsilon _{\mathrm r}$）と同じものである。混合の係数の上限は、その系の誘電関数のいちばん大きな値で決まる。

全エネルギー

ループが収束したら全エネルギーを求める。固有値の和 $\sum f\varepsilon$ は、運動エネルギーと $\int v_{\mathrm{eff}}n$ の和なので、式 (2.10.2) と比べると Hartree と交換相関を数え違えている。正しくは

$$E=\sum_{n\bm{k}\sigma}f_{n\bm{k}}\varepsilon _{n\bm{k}}-\int\bigl(v_{\mathrm H}+v_{\mathrm{xc}}\bigr)n\,\mathrm{d}^3r+E_{\mathrm H}[n]+E_{\mathrm{xc}}[n]+E_{\mathrm{II}}\tag{2.10.6}$$

である（$E_{\mathrm{II}}$ は原子核どうしの反発）。固有値の和は電子どうしの反発を「各電子が他のすべての電子から受ける」形で2回数えているので、その分を引き戻している。この式から、Kohn–Sham の固有値の和は全エネルギーでなく、固有値そのものにも全エネルギーの差という意味は直ちには無いことが分かる。固有値の意味は 2.10.4節で改めて扱う。