第14章コーン・シャム方程式の数値解法

第I部から第IV部までで，密度汎関数理論の原理(Hohenberg–Kohn定理，Kohn–Sham方程式，交換相関汎関数)と，周期系の電子構造(第12章のBlochの定理，第13章の状態密度とグリーン関数)を学んだ．理論の道具立ては揃った．しかしKohn–Sham(KS)方程式が解析的に解ける系はごくわずかであり，実在の分子や固体に適用するには計算機で数値的に解くしかない．本章から始まる第V部では，KS方程式が計算機の中で「どのように」解かれているのかを，基底関数の選択，行列の組み立て，全エネルギーの数値評価，自己無撞着(SCF)反復の収束加速まで分解して追いかける．章の後半では，球対称な原子のKS方程式を自分の手で解くための動径方程式ソルバを設計図のレベルまで具体化する．原子の計算は，第15章で学ぶ擬ポテンシャル構成の出発点でもある．

単位系: 本章では断りのない限りHartree原子単位系($\hbar=m_e=e=4\pi\varepsilon_0=1$)を用いる．エネルギーの単位は Hartree(1 Ha $=27.2114$ eV)，長さの単位は Bohr(1 Bohr $=0.529177$ Å)である．

この章で学ぶこと

KS方程式+SCFループを「非線形固有値問題」として定式化し，計算全体の流れを把握する
変分原理から一般化固有値問題 $HC=SC\varepsilon$ を一行も飛ばさずに導出し，平面波・局在基底・実空間グリッドの長所短所を比較する
局在基底実装の急所である全エネルギーのクーロン項再編成(中性原子分解)を恒等変形で確認する
SCF反復を線形化して単純混合の収束条件を求め，charge sloshing の起源，Kerker 混合，Pulay(DIIS)混合を導出する
球対称原子のKS方程式を変数分離し，対数メッシュとシューティング法で解く「自作原子DFTソルバ」の全体像を得る
コード間再現性の定量指標であるΔゲージの考え方を知る

14.1 何を計算するのか — 非線形固有値問題としてのKS方程式

まず，これから数値的に解こうとしている問題の構造を正確に見定めておこう．KS方程式は

\begin{equation} \left[-\frac{1}{2}\nabla^{2}+v_{\mathrm{eff}}[n](\rr)\right]\varphi_{i}(\rr) =\varepsilon_{i}\,\varphi_{i}(\rr) \label{eq:14-ks} \end{equation}

であり，有効ポテンシャルは電子密度 $n(\rr)$ の汎関数として

\begin{equation} v_{\mathrm{eff}}[n](\rr) = v_{\mathrm{ext}}(\rr) + \underbrace{\int \dd^3 r'\,\frac{n(\rr')}{\abs{\rr-\rr'}}}_{\displaystyle v_{\mathrm{H}}[n](\rr)} +\, v_{xc}[n](\rr) \label{eq:14-veff} \end{equation}

で与えられる．ここで $v_{\mathrm{ext}}$ は原子核(または第15章で導入する擬ポテンシャル)による外場，$v_{\mathrm{H}}$ はHartreeポテンシャル，$v_{xc}=\delta E_{xc}/\delta n$ は交換相関ポテンシャル(汎関数微分は第4章)である．密度は占有数 $f_i$ を使って

\begin{equation} n(\rr)=\sum_{i} f_i\,\abs{\varphi_i(\rr)}^{2}, \qquad \sum_i f_i = N_{\mathrm{e}} \label{eq:14-dens} \end{equation}

と軌道から組み立てられる(スピン分極系ではスピン $\sigma=\uparrow,\downarrow$ ごとに $n_\sigma$ を持つが，本章では記法の簡潔さのためスピン和を省略して書く)．

式 \eqref{eq:14-ks}–\eqref{eq:14-dens} をまとめて眺めると，この問題の特殊さが見える．ハミルトニアンが固有関数自身に依存しているのである．$v_{\mathrm{eff}}$ を作るには $n$ が要り，$n$ を作るには $\varphi_i$ が要り，$\varphi_i$ を求めるには $v_{\mathrm{eff}}$ が要る．すなわちKS方程式は非線形固有値問題であり，線形代数の固有値ソルバを1回呼べば済む問題ではない．標準的な解法は，この循環を反復で断ち切ることである．入力密度 $n_{\mathrm{in}}$ を仮定して $v_{\mathrm{eff}}[n_{\mathrm{in}}]$ を作り，固有値問題を(その反復の間だけ)線形な問題として解き，得られた軌道から出力密度 $n_{\mathrm{out}}$ を作る．この一連の写像を

\begin{equation} n_{\mathrm{out}} = F[n_{\mathrm{in}}] \label{eq:14-fixmap-intro} \end{equation}

と書けば，求めたい自己無撞着解は不動点 $n^{*}=F[n^{*}]$ である．$n_{\mathrm{in}}$ と $n_{\mathrm{out}}$ が(許容誤差の範囲で)一致するまで反復する手続きが自己無撞着場(SCF)ループである．どう反復すれば速く安定に不動点へ到達できるかは，それ自体が一つの理論問題であり，14.5節で詳しく解析する．

図14.1 KS方程式の数値解法の全体フロー．1回のSCF反復の中では $v_{\mathrm{eff}}$ を固定した線形な固有値問題を解き，外側のループで密度を更新する．収束しない場合は電荷混合(14.5節)を経て有効ポテンシャルの構築に戻る．

14.1.1 設計の三つの軸

図14.1のフローはどのDFTコードでも共通だが，各ステップの実現方法には大きな設計の自由度がある．実用コードの個性は，おおよそ次の三つの軸で分類できる．

(i) 内殻電子の扱い．すべての電子を陽に扱う全電子法と，化学的に不活性な内殻電子を核とまとめて有効ポテンシャルに置き換え，価電子だけを解く擬ポテンシャル法がある．擬ポテンシャルの理論と構成法は第15章で詳述する．本章では，外場の局所部分 $v_{\mathrm{core},I}(\rr)$ が遠方で $-Z_I/r$($Z_I$ は全電子法なら核電荷，擬ポテンシャル法なら価電子数)と振る舞う，という事実だけを使う．

(ii) 基底関数(離散化)の選択．波動関数を有限個の数で表現する方法である．平面波，原子軌道様の局在関数，実空間グリッドが三大流派であり，14.2節で変分原理から統一的に定式化して比較する．

(iii) 求解アルゴリズム．固有値問題の解き方(直接対角化か反復対角化か)，Poisson方程式の解き方，SCFの収束加速(電荷混合)などである．本章では特にSCF収束(14.5節)を掘り下げる．

表14.1 KS方程式の数値解法の代表的な組み合わせ
方式	内殻電子	基底	実装例	特徴
全電子・混合基底 ((L)APW+lo など)	全て扱う	格子間は平面波，原子球内は数値動径関数	WIEN2k, FLEUR, ELK	最高精度の参照値を与える．内殻準位も計算可能．コストは大きい
擬ポテンシャル・平面波	価電子のみ	平面波	VASP, Quantum ESPRESSO	カットオフ1パラメータで系統的収束．周期系向き．扱いやすい
擬ポテンシャル・局在基底	価電子のみ	原子軌道様の局在関数	OpenMX, SIESTA, CP2K	少数基底・疎行列・$O(N)$法と整合．基底設計に経験を要する
実空間グリッド	主に価電子	差分格子(基底を使わない)	Octopus, GPAW(FD)	格子間隔で系統的収束．大規模並列に向く

物理的意味

「基底を選ぶ」ことは「ヒルベルト空間のどの部分空間で変分するか」を選ぶことである．14.2節で示すように，有限基底のKS計算はすべてRayleigh–Ritzの変分計算であり，基底を増やせば全エネルギーは(同じハミルトニアンに対して)単調に真の値へ近づく．平面波はこの収束を1個のパラメータ(カットオフエネルギー)で制御できる点が強みであり，局在基底は少数の関数で化学的に意味のある部分空間を張れる点が強みである．どちらが「正しい」かではなく，どの物理系にどちらの部分空間が効率的か，という工学的な選択の問題である．

14.1.2 計算コストの見取り図

基底の数を $M$，電子数(または原子数)を $N$ とすると，図14.1の各ステップのコストは大まかに次のようになる．行列要素の組み立ては平面波ではFFTを使って $O(M\log M)$，局在基底では行列が疎になるため $O(N)$ で済む．一方，固有値問題の直接対角化は $O(M^{3})$ であり，系が大きくなるとここが律速になる．平面波では原子あたり $10^{2}$–$10^{3}$ 個の基底が必要なのに対し，よく設計された局在基底は原子あたり $10$–$50$ 個で済む．この差が，大規模系や電子輸送計算で局在基底実装が好まれる理由である．さらに局在基底ではハミルトニアンの疎性を利用した $O(N)$ 法への道が開ける(本書では扱わないが，14.3節の基底設計はその土台である)．

SCFが収束したら，全エネルギー $E_{\mathrm{tot}}$(14.4節)と原子に働く力 $\bm{F}_I=-\partial E_{\mathrm{tot}}/\partial\bm{\tau}_I$ を計算する．力の計算では，基底関数が原子位置に張り付いている場合(局在基底)，Hellmann–Feynman項に加えて基底の移動に由来する補正(Pulay力)が必要になることを注意しておく．

14.2 基底展開と一般化固有値問題

この節では，KS方程式 \eqref{eq:14-ks} を「有限個の数を決める問題」に落とす．方針は明快である．波動関数を有限個の基底関数で展開し，変分原理(Rayleigh商の停留条件)を課すと，展開係数に対する行列方程式 — 一般化固有値問題 — が得られる．これがあらゆる基底型DFTコードの心臓部である．

14.2.1 有限基底展開

互いに一次独立な $M$ 個の基底関数 $\{\chi_\mu(\rr)\}_{\mu=1}^{M}$ を用意し，KS軌道を

\begin{equation} \varphi_i(\rr)=\sum_{\mu=1}^{M} c_{\mu i}\,\chi_\mu(\rr) \label{eq:14-expand} \end{equation}

と展開する．基底は一般に直交していなくてよい(局在原子軌道は互いに重なる)．ここで2つの $M\times M$ 行列，ハミルトニアン行列 $H$ と重なり行列 $S$ を定義する:

\begin{equation} H_{\mu\nu}=\int \dd^3 r\,\chi_\mu^{*}(\rr)\,\hat{h}_{\mathrm{KS}}\,\chi_\nu(\rr), \qquad S_{\mu\nu}=\int \dd^3 r\,\chi_\mu^{*}(\rr)\,\chi_\nu(\rr), \label{eq:14-hs-def} \end{equation}

ただし $\hat{h}_{\mathrm{KS}}=-\frac{1}{2}\nabla^2+v_{\mathrm{eff}}$ である．$\hat{h}_{\mathrm{KS}}$ がエルミートなので $H_{\mu\nu}=H_{\nu\mu}^{*}$，また定義から明らかに $S_{\mu\nu}=S_{\nu\mu}^{*}$ であり，$H$ も $S$ もエルミート行列である．

1個の軌道 $\varphi$(添字 $i$ をしばらく省略)のエネルギー期待値とノルムを，展開 \eqref{eq:14-expand} を代入して係数で表そう．積分と和の順序交換(有限和なので常に許される)により

\begin{align} \braket{\varphi|\hat{h}_{\mathrm{KS}}|\varphi} &=\int \dd^3 r\,\Big(\sum_\mu c_{\mu}\chi_\mu(\rr)\Big)^{*}\hat{h}_{\mathrm{KS}}\Big(\sum_\nu c_{\nu}\chi_\nu(\rr)\Big) \label{eq:14-exp-h-1}\\ &=\sum_{\mu\nu} c_{\mu}^{*}c_{\nu}\int \dd^3 r\,\chi_\mu^{*}(\rr)\,\hat{h}_{\mathrm{KS}}\,\chi_\nu(\rr) \label{eq:14-exp-h-2}\\ &=\sum_{\mu\nu} c_{\mu}^{*}\,H_{\mu\nu}\,c_{\nu} =\bm{c}^{\dagger}H\bm{c}. \label{eq:14-exp-h} \end{align}

1行目から2行目へは，複素共役を各項に分配し($\hat{h}_{\mathrm{KS}}$ は係数に作用しない線形演算子であることを使い)，二重和を積分の外に出した．3行目は定義 \eqref{eq:14-hs-def} を代入しただけである．全く同じ手順で

\begin{equation} \braket{\varphi|\varphi}=\sum_{\mu\nu}c_\mu^{*}S_{\mu\nu}c_\nu=\bm{c}^{\dagger}S\bm{c} \label{eq:14-exp-s} \end{equation}

を得る．ここで $\bm{c}=(c_1,\dots,c_M)^{\mathrm T}$ は係数を並べた列ベクトルである．

14.2.2 Rayleigh商の停留条件から $H\bm{c}=\varepsilon S\bm{c}$ へ

変分原理によれば，規格化条件 $\braket{\varphi|\varphi}=1$ の下で $\braket{\varphi|\hat{h}_{\mathrm{KS}}|\varphi}$ を停留にする $\varphi$ が固有状態である．拘束条件付きの停留問題なので，Lagrangeの未定乗数 $\varepsilon$ を導入して

\begin{equation} L(\bm{c},\bm{c}^{*},\varepsilon) =\bm{c}^{\dagger}H\bm{c} -\varepsilon\left(\bm{c}^{\dagger}S\bm{c}-1\right) =\sum_{\mu\nu}c_\mu^{*}H_{\mu\nu}c_\nu -\varepsilon\Big(\sum_{\mu\nu}c_\mu^{*}S_{\mu\nu}c_\nu-1\Big) \label{eq:14-lagrangian} \end{equation}

を係数について停留化する．未定乗数 $\varepsilon$ を使う理由は，拘束(規格化)を満たす方向だけの変分を，拘束なしの変分に置き換えるためである(第4章のLagrange未定乗数法と同じ論理)．係数は複素数なので，微分の取り扱いを数学ノートで確認しておく．

数学ノート:複素係数による微分($c$ と $c^{*}$ を独立に扱ってよい理由)

複素変数 $c=a+ib$($a,b$ は実数)の実数値関数 $f(a,b)$ の停留条件は $\partial f/\partial a=0$ かつ $\partial f/\partial b=0$ である．ここで形式的な微分演算子

$$ \frac{\partial}{\partial c}\equiv\frac{1}{2}\left(\frac{\partial}{\partial a}-i\frac{\partial}{\partial b}\right), \qquad \frac{\partial}{\partial c^{*}}\equiv\frac{1}{2}\left(\frac{\partial}{\partial a}+i\frac{\partial}{\partial b}\right) $$

を定義する(Wirtinger微分)．この定義に従うと，$c=a+ib$，$c^{*}=a-ib$ に対して

$$ \frac{\partial c}{\partial c^{*}}=\frac{1}{2}\left(\frac{\partial}{\partial a}+i\frac{\partial}{\partial b}\right)(a+ib) =\frac{1}{2}(1+i\cdot i)=\frac{1}{2}(1-1)=0, \qquad \frac{\partial c^{*}}{\partial c^{*}}=\frac{1}{2}(1-i\cdot i)=1 $$

となる．つまりこの微分規則の下では $c$ と $c^{*}$ は互いに独立な変数のように振る舞う．さらに，$\partial f/\partial a$ と $\partial f/\partial b$ の同時消滅は $\partial f/\partial c^{*}=0$(1本の複素方程式=実2本)と同値である．したがって「$c^{*}$ で微分してゼロと置く」ことは，実部・虚部それぞれで停留化することと完全に等価であり，計算は圧倒的に簡単になる．

では \eqref{eq:14-lagrangian} を $c_\kappa^{*}$ で微分する．第1項は

\begin{align} \frac{\partial}{\partial c_\kappa^{*}}\sum_{\mu\nu}c_\mu^{*}H_{\mu\nu}c_\nu &=\sum_{\mu\nu}\frac{\partial c_\mu^{*}}{\partial c_\kappa^{*}}H_{\mu\nu}c_\nu =\sum_{\mu\nu}\delta_{\mu\kappa}H_{\mu\nu}c_\nu =\sum_{\nu}H_{\kappa\nu}c_\nu. \label{eq:14-dterm1} \end{align}

ここで，数学ノートの規則により $c_\nu$ は $c_\kappa^{*}$ と独立なので微分にかからず，$\partial c_\mu^{*}/\partial c_\kappa^{*}=\delta_{\mu\kappa}$ だけが生き残った．第2項も同じ構造なので

\begin{equation} \frac{\partial}{\partial c_\kappa^{*}}\left[\varepsilon\Big(\sum_{\mu\nu}c_\mu^{*}S_{\mu\nu}c_\nu-1\Big)\right] =\varepsilon\sum_{\nu}S_{\kappa\nu}c_\nu. \label{eq:14-dterm2} \end{equation}

停留条件 $\partial L/\partial c_\kappa^{*}=0$($\kappa=1,\dots,M$)は，\eqref{eq:14-dterm1} から \eqref{eq:14-dterm2} を引いてゼロと置くことにより

\begin{equation} \sum_{\nu}H_{\kappa\nu}c_\nu=\varepsilon\sum_{\nu}S_{\kappa\nu}c_\nu \qquad(\kappa=1,\dots,M) \label{eq:14-stationary} \end{equation}

となる($c_\kappa$ での微分は，この式のエルミート共役を与えるだけで新しい条件を生まない．$H$，$S$ がエルミートだからである)．これをベクトル・行列記法でまとめ，軌道の添字 $i$ を復活させると

\begin{equation} H\bm{c}_i=\varepsilon_i\,S\,\bm{c}_i \qquad\Longleftrightarrow\qquad HC=SC\,\mathrm{diag}(\varepsilon_1,\dots,\varepsilon_M) \label{eq:14-gevp} \end{equation}

を得る．ここで $C$ は固有ベクトル $\bm{c}_i$ を列として並べた $M\times M$ 行列である．重なり行列 $S$ が単位行列でない点だけが通常の固有値問題と異なり，これを一般化固有値問題と呼ぶ．これが得られたもの:KS方程式の数値解法とは，SCFの各反復で一般化固有値問題 \eqref{eq:14-gevp} を解くことである．

注意:2種類の $S$ を混同しないこと

第15章のウルトラソフト擬ポテンシャルとPAW法でも，$\hat H\left|\phi\right\rangle=\varepsilon\,\hat S\left|\phi\right\rangle$ という形の上では同じ一般化固有値問題が現れる(15.7.2項)．しかし両者の $S$ は起源がまったく違う．本節の $S_{\mu\nu}=\braket{\chi_\mu|\chi_\nu}$ は基底が直交していないことだけから来るもので，平面波基底なら $S=I$ に戻る．一方ウルトラソフトの $\hat S=1+\sum_{ij}q_{ij}\left|\beta_i\right\rangle\left\langle\beta_j\right|$ は，ノルム保存を捨てた擬変換そのものが持ち込む計量の変形であり，平面波のような完全直交基底を使っても $\hat S\neq1$ のまま残る．両者は独立な効果なので，非直交局在基底とウルトラソフト擬ポテンシャルを併用すれば，2つが重なった $S$ を扱うことになる．本章で以下 $S$ と書けば，つねに前者(重なり行列)を指す．

未定乗数 $\varepsilon$ の意味も確認しておこう．\eqref{eq:14-stationary} の両辺に $c_\kappa^{*}$ を掛けて $\kappa$ について和をとると

\begin{align} \sum_{\kappa\nu}c_\kappa^{*}H_{\kappa\nu}c_\nu &=\varepsilon\sum_{\kappa\nu}c_\kappa^{*}S_{\kappa\nu}c_\nu &&\text{(両辺に $\textstyle\sum_\kappa c_\kappa^{*}$ を作用)}\nonumber\\ \Longrightarrow\quad \varepsilon&=\frac{\bm{c}^{\dagger}H\bm{c}}{\bm{c}^{\dagger}S\bm{c}} =\frac{\braket{\varphi|\hat{h}_{\mathrm{KS}}|\varphi}}{\braket{\varphi|\varphi}} \label{eq:14-rayleigh-val} \end{align}

となる．2行目へは，$\bm{c}^{\dagger}S\bm{c}\neq 0$(後述のように $S$ は正定値)で割り，\eqref{eq:14-exp-h}, \eqref{eq:14-exp-s} を使った．すなわち未定乗数は軌道のRayleigh商そのもの，つまりKS軌道エネルギーである．

異なる固有値に属する解の直交性も導いておく．\eqref{eq:14-gevp} を状態 $i$ について書き，左から $\bm{c}_j^{\dagger}$ を掛ける:

\begin{align} \bm{c}_j^{\dagger}H\bm{c}_i&=\varepsilon_i\,\bm{c}_j^{\dagger}S\bm{c}_i. \label{eq:14-orth1} \end{align}

次に状態 $j$ の式 $H\bm{c}_j=\varepsilon_j S\bm{c}_j$ のエルミート共役をとる．$H^{\dagger}=H$，$S^{\dagger}=S$，また $\varepsilon_j$ が実数であること(\eqref{eq:14-rayleigh-val} でエルミート行列の期待値の商だから)を使うと $\bm{c}_j^{\dagger}H=\varepsilon_j\bm{c}_j^{\dagger}S$ となり，これを右から $\bm{c}_i$ に作用させて

\begin{align} \bm{c}_j^{\dagger}H\bm{c}_i&=\varepsilon_j\,\bm{c}_j^{\dagger}S\bm{c}_i. \label{eq:14-orth2} \end{align}

\eqref{eq:14-orth1} から \eqref{eq:14-orth2} を辺々引くと

\begin{equation} 0=(\varepsilon_i-\varepsilon_j)\,\bm{c}_j^{\dagger}S\bm{c}_i \quad\Longrightarrow\quad \bm{c}_j^{\dagger}S\bm{c}_i=0\quad(\varepsilon_i\neq\varepsilon_j). \label{eq:14-sorth} \end{equation}

つまり固有ベクトルは $S$ を計量とする内積で直交する．これは $\braket{\varphi_j|\varphi_i}=\bm{c}_j^{\dagger}S\bm{c}_i=0$，すなわち軌道自身の直交性に他ならない．

定理:Rayleigh–Ritz の変分性と基底拡大の単調性

(1) 一般化固有値問題 \eqref{eq:14-gevp} の最低固有値 $\varepsilon_1^{(M)}$ は，真の最低固有値 $\varepsilon_1$ の上界である:$\varepsilon_1^{(M)}\ge\varepsilon_1$．これはRayleigh商の最小値が部分空間の制限により真の最小値以上になることから直ちに従う．
(2) さらに，基底を1個追加して $M\to M+1$ とすると，各固有値は単調に下がるか不変である:$\varepsilon_k^{(M+1)}\le\varepsilon_k^{(M)}$($k=1,\dots,M$)．この「入れ子の部分空間での固有値のはさみうち」はHylleraas–Undheim–MacDonaldの定理として知られる(証明は文献 [MacDonald 1933] に譲る)．実用上の含意は重要で，同一のハミルトニアンに対しては基底を増やすほど固有値・全エネルギーは系統的に改善する．「基底を増やして結果が変わらなくなるまで確かめる」という収束テストの正当性はこの定理が保証している．

14.2.3 重なり行列の性質と標準固有値問題への変換

$S$ の最重要性質は正定値性である．任意の係数ベクトル $\bm{c}\neq\bm{0}$ に対して

\begin{equation} \bm{c}^{\dagger}S\bm{c} =\sum_{\mu\nu}c_\mu^{*}S_{\mu\nu}c_\nu =\int \dd^3 r\,\Big|\sum_\mu c_\mu\chi_\mu(\rr)\Big|^{2} =\int \dd^3 r\,\abs{\varphi(\rr)}^{2}\ \ge 0, \label{eq:14-spos} \end{equation}

であり(2番目の等号は \eqref{eq:14-exp-s} の導出をそのまま使った)，等号は $\varphi(\rr)=0$(ほとんど至るところ)のときに限る．基底が一次独立ならば $\varphi=0$ となる非自明な $\bm{c}$ は存在しないので，$\bm{c}^{\dagger}S\bm{c}\gt 0$，すなわち $S$ は正定値で可逆である．

正定値性を使うと \eqref{eq:14-gevp} は標準形に変換できる．$S$ はエルミートなのでユニタリ行列 $U$ で対角化でき，$S=U\,\mathrm{diag}(s_1,\dots,s_M)\,U^{\dagger}$，$s_\mu\gt0$ と書ける．そこで

\begin{equation} S^{-1/2}\equiv U\,\mathrm{diag}(s_1^{-1/2},\dots,s_M^{-1/2})\,U^{\dagger} \label{eq:14-shalf-def} \end{equation}

を定義する(Löwdin直交化)．$S^{-1/2}SS^{-1/2}=U\,\mathrm{diag}(s^{-1/2}s\,s^{-1/2})\,U^{\dagger}=UU^{\dagger}=1$ を確かめておこう(対角行列同士の積は成分ごとの積，$U^{\dagger}U=1$ を2回使った)．新しい変数 $\bm{d}_i\equiv S^{1/2}\bm{c}_i$ を導入し，\eqref{eq:14-gevp} の左から $S^{-1/2}$ を掛けると

\begin{align} S^{-1/2}H\bm{c}_i&=\varepsilon_i\,S^{-1/2}S\bm{c}_i &&\text{(左から $S^{-1/2}$)}\nonumber\\ S^{-1/2}H\underbrace{S^{-1/2}S^{1/2}}_{=1}\bm{c}_i&=\varepsilon_i\,S^{-1/2}S\underbrace{S^{-1/2}S^{1/2}}_{=1}\bm{c}_i &&\text{(恒等式 $S^{-1/2}S^{1/2}=1$ を挿入)}\nonumber\\ \big(S^{-1/2}HS^{-1/2}\big)\,\bm{d}_i&=\varepsilon_i\,\big(S^{-1/2}SS^{-1/2}\big)\,\bm{d}_i=\varepsilon_i\,\bm{d}_i \label{eq:14-loewdin-transform} \end{align}

となる．$\tilde{H}\equiv S^{-1/2}HS^{-1/2}$ はエルミート($\tilde H^\dagger = S^{-1/2}H^\dagger S^{-1/2}=\tilde H$，$S^{-1/2}$ がエルミートであることを使った)なので，通常のエルミート固有値ソルバがそのまま使える．実装ではコレスキー分解 $S=LL^{\dagger}$ を使う同等の変換が主流である．

物理的意味:$S$ の条件数と「基底の張りすぎ」

\eqref{eq:14-loewdin-transform} には $s_\mu^{-1/2}$ が現れる．もし基底がほとんど一次従属(例えば裾の長い局在関数を密に並べた場合)なら，$S$ の最小固有値 $s_{\min}$ が $10^{-8}$ のような微小値になり，$s_{\min}^{-1/2}\sim10^{4}$ が数値誤差を増幅する．これが「基底は多いほどよい」とは限らない実装上の理由である．局在基底コードでは，$s_\mu$ が閾値以下の成分を捨てる(部分空間を切り詰める)処理が組み込まれていることが多い．基底のカットオフ半径を伸ばしすぎたときに計算が不安定になるのは，たいていこの過完全性が原因である．

14.2.4 周期系では:Bloch和と $\kk$ 点ごとの固有値問題

結晶では第12章のBlochの定理に従い，局在基底 $\chi_\mu(\rr-\bm{\tau}_\mu)$ から Bloch 和

\begin{equation} \chi_{\mu}^{(\kk)}(\rr)=\frac{1}{\sqrt{N_{\mathrm{cell}}}}\sum_{\bm{R}}e^{i\kk\cdot\bm{R}}\,\chi_\mu(\rr-\bm{\tau}_\mu-\bm{R}) \label{eq:14-bloch-sum} \end{equation}

を作って展開する($\bm{R}$ は格子ベクトル)．ハミルトニアンが格子並進と可換なので，異なる $\kk$ は混ざらず，固有値問題は $\kk$ 点ごとにブロック対角化される:

\begin{equation} H(\kk)\,C(\kk)=S(\kk)\,C(\kk)\,\mathrm{diag}\big(\varepsilon_1(\kk),\dots\big). \label{eq:14-gevp-k} \end{equation}

行列の次元は「単位胞内の基底の数」であり，結晶が無限に大きくても固有値問題は有限次元で済む．これがBlochの定理の計算上の恩恵である．以下，本章では記法の簡潔さのため $\kk$ 依存性を省いて書くが，すべての式は $\kk$ 点ごとに成り立つ．

14.2.5 基底関数の三大流派

ここまでの定式化は基底の選び方によらない．では実際に何を $\chi_\mu$ に選ぶか．図14.2に三大流派を示す．

図14.2 基底(離散化)の三大流派の模式図．(a) 平面波:カットオフ以下の波数の正弦波をすべて使う．(b) 局在基底:原子に張り付いた関数の線形結合．重なり行列 $S$ が非自明になる．(c) 実空間グリッド:格子点上の値で波動関数を表し，微分は差分で近似する．

(a) 平面波基底．体積 $V$ の周期セルで，規格化された平面波 $\chi_{\bm{G}}(\rr)=V^{-1/2}e^{i(\kk+\bm{G})\cdot\rr}$($\bm{G}$ は逆格子ベクトル)を基底にとる．平面波同士の直交性から $S_{\bm{G}\bm{G}'}=\delta_{\bm{G}\bm{G}'}$，つまり重なり行列は単位行列で，\eqref{eq:14-gevp} は最初から標準固有値問題である．基底の大きさはカットオフエネルギー $E_{\mathrm{cut}}$ ただ1つで制御する:

\begin{equation} \frac{1}{2}\abs{\kk+\bm{G}}^{2}\le E_{\mathrm{cut}}. \label{eq:14-pw-cutoff} \end{equation}

使用する平面波の数を見積もろう．逆格子点は逆空間で密度 $V/(2\pi)^{3}$ で分布する(第3章のフェルミ球の状態数の数え方と同じ)．半径 $q_c=\sqrt{2E_{\mathrm{cut}}}$ の球内の点の数は

\begin{equation} N_{\mathrm{PW}}\simeq\frac{4}{3}\pi q_c^{3}\times\frac{V}{(2\pi)^{3}} =\frac{V\,(2E_{\mathrm{cut}})^{3/2}}{6\pi^{2}} \label{eq:14-pw-count} \end{equation}

である．$E_{\mathrm{cut}}$ を上げれば厳密解に系統的に収束する(前掲の定理の (2))．これが平面波の最大の美点である．弱点は，波動関数が原子核近傍で鋭く振動する場合(遷移金属の $3d$，酸素の $2p$ など)に $E_{\mathrm{cut}}$ が莫大になることで，擬ポテンシャル(第15章)による「軟化」が事実上必須になる．また真空領域にも一様に基底を張るため，孤立分子や表面では損をする．

例:シリコン結晶の平面波数

Siダイヤモンド構造の慣用格子定数は $a=5.431$ Å $=10.26$ Bohr，fcc単位胞(2原子)の体積は $V=a^{3}/4=270$ Bohr$^{3}$ である．標準的なカットオフ $E_{\mathrm{cut}}=12.5$ Ha(25 Ry)に対して $q_c=\sqrt{2\times12.5}=5.0$ Bohr$^{-1}$ であり，\eqref{eq:14-pw-count} より $$N_{\mathrm{PW}}\simeq\frac{270\times125}{6\pi^{2}}\approx570$$ つまり原子あたり約285個の平面波を使う．一方，後述の局在基底なら原子あたり十数個で同等の価電子帯が記述できる．ただし平面波の $N_{\mathrm{PW}}$ は $E_{\mathrm{cut}}$ を上げるだけで系統的に厳密解へ向かうのに対し，局在基底の改善には関数形の設計が要る．この対比が両流派の性格をよく表している．

(b) 局在基底(LCAO/LCPAO)．原子 $I$ の位置 $\bm{\tau}_I$ に中心を持つ原子軌道様の関数 $\chi_{I\alpha}(\rr-\bm{\tau}_I)$(角運動量 $lm$ と動径自由度をまとめて $\alpha$ と書く)の線形結合で展開する．量子化学のGauss基底，Slater基底，そして数値表の動径関数を使う数値原子軌道(擬原子軌道，LCPAO)がこの流派である．長所は，(i) 原子あたり10–50個という圧倒的に少ない基底数，(ii) 有限のカットオフ半径を持たせれば $H$，$S$ が疎行列になり，大規模系・$O(N)$法・非平衡Green関数輸送計算と整合すること，(iii) 係数 $c_{\mu i}$ がそのまま化学的な解釈(どの原子のどの軌道が寄与するか)を与えることである．短所は，系統的収束の単一パラメータが存在しないこと，そして次の BSSE である．

数学ノート:基底関数重ね合わせ誤差(BSSE)

局在基底では，基底が原子に張り付いているため，系の形が変わると変分空間そのものが変わる．2分子 A，B の結合エネルギーを $\Delta E=E_{AB}-E_A-E_B$ で見積もる場合を考える．二量体 AB の計算では，Aの電子はBの基底関数「も」借りて自分の波動関数を改善できる(逆も同様)．単量体の計算ではこの借用ができない．したがって二量体だけが余分に変分的に得をし，$\Delta E$ は結合を過大評価する方向にずれる．これが基底関数重ね合わせ誤差(basis set superposition error, BSSE)である．標準的な補正は counterpoise 法 [Boys–Bernardi 1970]:単量体 A を「Bの位置に基底だけ置いた(原子核も電子もない)ゴースト原子」付きで計算し直し，同じ変分空間同士で差をとる．BSSEは基底を完全系に近づければ消える誤差であり，平面波基底(空間に一様に張られる)には存在しない．

(c) 実空間グリッド．波動関数を格子点上の値 $\varphi(\rr_p)$ で表し，ラプラシアンを有限差分，例えば2次精度なら

\begin{equation} \frac{\dd^2\varphi}{\dd x^2}\bigg|_{x_p}\approx\frac{\varphi(x_{p+1})-2\varphi(x_p)+\varphi(x_{p-1})}{h^{2}} \label{eq:14-fd} \end{equation}

で近似する(実際は6次〜8次精度の広いステンシルを使う)．格子間隔 $h$ を小さくすれば系統的に収束し，通信が局所的なので超並列計算機と相性がよい．原子が格子点に対してどこにあるかでエネルギーがわずかに波打つ「eggbox誤差」に注意が要る．なお，平面波法も内部ではFFTで実空間グリッドを併用しており，両者は親戚である．

表14.2 基底三流派の比較
	平面波	局在基底(LCAO/LCPAO)	実空間グリッド
収束制御	$E_{\mathrm{cut}}$ 1個で系統的	基底数・カットオフ半径(設計が必要)	格子間隔 $h$ で系統的
基底数/原子	$10^{2}$–$10^{3}$	$10$–$50$	グリッド点 $10^{3}$–$10^{4}$
$S$ 行列	単位行列	非対角(正定値)	(ほぼ)単位行列
行列の疎性	密(ただしFFTで作用は速い)	疎(カットオフ半径)	疎(差分ステンシル)
BSSE	なし	あり(counterpoise補正)	なし
得意な系	周期結晶・金属	大規模系・分子・輸送・$O(N)$	孤立系・超並列
力の計算	Hellmann–Feynmanのみ	Pulay補正が必要	eggbox誤差に注意

以上で「KS方程式をどんな行列方程式に落とすか」が確定した．次節では，三流派のうち局在基底，特に数値擬原子軌道(LCPAO)の設計を詳しく見る．少数の関数で高精度を出すための工夫が詰まっているからである．

14.3 基底関数の選択 — 擬原子軌道(LCPAO)の設計

前節で三流派の性格を表14.2に並べたが，コードを設計する立場で本当に効いてくるのは「行列を作るコスト・掛けるコストが系のサイズ $N$ にどう依存するか」という定量的な事情である．この節ではまずその点を具体的な数値で確認し(14.3.1)，続いて数値擬原子軌道(LCPAO)基底をどう作り，どう最適化するかを，閉じ込めポテンシャルの必要性から変分最適化の勾配式まで具体化する(14.3.2〜14.3.4)．少数の関数で高精度を出すための工夫がここに集約されている．

14.3.1 三流派のコストを数値で比べる

(a) 平面波:行列を「作らずに」掛ける．平面波基底では \eqref{eq:14-hs-def} の行列要素が解析的に書ける．運動エネルギー項は平面波の直交性から対角，局所ポテンシャル項はそのFourier成分になり

\begin{equation} H_{\bm{G}'\bm{G}}=\frac{1}{2}\abs{\kk+\bm{G}}^{2}\,\delta_{\bm{G}'\bm{G}} +\tilde v_{\mathrm{eff}}(\bm{G}'-\bm{G}), \qquad \tilde v_{\mathrm{eff}}(\bm{G})\equiv\frac{1}{\Omega}\int_{\Omega}\dd^3 r\,v_{\mathrm{eff}}(\rr)\,e^{-i\bm{G}\cdot\rr} \label{eq:14-pw-hmat} \end{equation}

となる($\Omega$ は単位胞の体積)．しかしこの行列は密である．$N_{\mathrm{PW}}=10^{5}$ なら要素数は $10^{10}$ 個，倍精度複素数で $160$ GB になり，そもそも記憶できない．したがって平面波法では行列を陽に作らない．必要なのは「与えられたベクトル $\bm{c}$ に $H$ を掛けた結果」だけであり，それはFFTで計算できる．理由を確かめておこう．波動関数を

\begin{equation} \psi_{\kk}(\rr)=\frac{1}{\sqrt{\Omega}}\sum_{\bm{G}}c_{\bm{G}}\,e^{i(\kk+\bm{G})\cdot\rr}, \qquad v_{\mathrm{eff}}(\rr)=\sum_{\bm{G}}\tilde v_{\mathrm{eff}}(\bm{G})\,e^{i\bm{G}\cdot\rr} \label{eq:14-pw-expand} \end{equation}

と展開し，実空間での積 $v_{\mathrm{eff}}\psi_{\kk}$ を計算すると

\begin{align} v_{\mathrm{eff}}(\rr)\,\psi_{\kk}(\rr) &=\frac{1}{\sqrt{\Omega}}\sum_{\bm{G}}\sum_{\bm{G}'}\tilde v_{\mathrm{eff}}(\bm{G}')\,c_{\bm{G}}\,e^{i(\kk+\bm{G}+\bm{G}')\cdot\rr} \label{eq:14-conv-1}\\ &=\frac{1}{\sqrt{\Omega}}\sum_{\bm{G}''}\left[\sum_{\bm{G}}\tilde v_{\mathrm{eff}}(\bm{G}''-\bm{G})\,c_{\bm{G}}\right]e^{i(\kk+\bm{G}'')\cdot\rr} \label{eq:14-conv-2} \end{align}

となる．1行目は \eqref{eq:14-pw-expand} の2つの和を単に掛け合わせたもの，2行目では和の変数を $\bm{G}''=\bm{G}+\bm{G}'$ に取り替えた(逆格子ベクトルの集合は加法について閉じているので，$\bm{G}'$ を走らせることと $\bm{G}''$ を走らせることは同じである)．角括弧の中身は \eqref{eq:14-pw-hmat} のポテンシャル項と $\bm{c}$ の積そのものである．つまり

物理的意味:なぜFFTが効くのか

運動エネルギーは逆空間で対角，局所ポテンシャルは実空間で対角である．$H\bm{c}$ を計算するには，$\bm{c}$ をFFTで実空間に運び($O(M\log M)$)，格子点ごとに $v_{\mathrm{eff}}(\rr_p)$ を掛け($O(M)$)，FFTで逆空間に戻して($O(M\log M)$)運動エネルギー項を足せばよい．$O(M^{2})$ の行列–ベクトル積が $O(M\log M)$ になる．あとは Davidson 法や共役勾配法のような反復対角化を使えば，$M$ 本すべてではなく占有軌道 $N_{\mathrm{occ}}$ 本だけを求めればよく，全体のコストは $O(N_{\mathrm{occ}}M\log M)$ + 直交化の $O(N_{\mathrm{occ}}^{2}M)$ に落ちる．平面波法が $10^{5}$ 個規模の基底を扱えるのは，この「行列を作らない」戦略のおかげである．

FFT格子の目の粗さも決まっている．格子間隔 $h$ の実空間格子で表現できる最短波長は $2h$(1周期に2点)であり，対応する最大波数は

\begin{equation} q_{\max}=\frac{2\pi}{2h}=\frac{\pi}{h} \label{eq:14-nyquist} \end{equation}

である(Nyquist波数)．波動関数は $\abs{\kk+\bm{G}}\le q_c=\sqrt{2E_{\mathrm{cut}}}$ までしか成分を持たないが，密度 $n=\abs{\psi}^{2}$ とポテンシャルは畳み込みにより $2q_c$ まで成分を持つ．したがって折り返し誤差(エイリアシング)を避けるには $q_{\max}\ge 2q_c$，すなわち

\begin{equation} h\le\frac{\pi}{2q_c}=\frac{\pi}{2\sqrt{2E_{\mathrm{cut}}}} \label{eq:14-grid-spacing} \end{equation}

が必要である．$E_{\mathrm{cut}}=12.5$ Ha なら $q_c=5.0$ Bohr$^{-1}$，$h\le 0.31$ Bohr となる．前掲のSi(2原子，$\Omega=270$ Bohr$^{3}$)なら格子点数は $270/0.31^{3}\approx 9\times10^{3}$，原子あたり約 $4\times10^{3}$ 点であり，表14.2の見積もりと合う．

(b) 平面波の弱点は「真空にも基底を張る」こと．\eqref{eq:14-pw-count} が示すように $N_{\mathrm{PW}}\propto \Omega$ であり，平面波はセルのどこに原子があるかを知らない．厚さ $10$ Bohr のスラブを真空層 $20$ Bohr とともに周期セルに入れると，$\Omega$ の $2/3$ は真空であり，平面波の $2/3$ は「何もない場所で波動関数がゼロであること」を記述するために費やされる．孤立分子ではこの浪費はさらに激しい．局在基底は原子に張り付いているので，この意味での無駄が原理的に存在しない．

(c) 局在基底:カットオフ半径が疎行列性を生む．基底 $\chi_{I\alpha}$ が原子 $I$ から半径 $r_c$ の外で恒等的にゼロなら，

\begin{equation} \abs{\bm{\tau}_I-\bm{\tau}_J}\gt r_c^{I}+r_c^{J} \quad\Longrightarrow\quad S_{I\alpha,J\beta}=0,\quad H_{I\alpha,J\beta}=0 \label{eq:14-sparsity} \end{equation}

である(2つの関数の台が重ならなければ積分は恒等的にゼロ．$H$ については局所ポテンシャル項も同じ理由で消え，非局所擬ポテンシャル項も射影子の台が有限なので同様である)．1行あたりの非ゼロ要素数は「半径 $r_c^{I}+r_c^{J}$ の球に入る原子数 × 原子あたりの軌道数」であり，これは系の大きさ $N$ に依らない定数である．ゆえに $H$，$S$ の格納量と組み立てコストは $O(N)$ になる．

例:Si結晶での非ゼロ要素数

Siの慣用格子(1辺 $a=10.26$ Bohr，8原子)から原子数密度は $8/a^{3}=7.4\times10^{-3}$ 原子/Bohr$^{3}$ である．カットオフ半径 $r_c=6.0$ Bohr の基底を使うと，行列要素が非ゼロになりうるのは $\abs{\bm{\tau}_I-\bm{\tau}_J}\le 12$ Bohr の相手だけであり，その球の体積は $\frac{4}{3}\pi(12)^{3}=7.24\times10^{3}$ Bohr$^{3}$，含まれる原子数は $$7.24\times10^{3}\times7.4\times10^{-3}\approx 54\ \text{個}$$ である．原子あたり13軌道(後述の s2p2d1)なら1行あたり $54\times13\approx700$ 個の非ゼロ要素で，これは $N$ に依らない．$10^{4}$ 原子($M=1.3\times10^{5}$)の系なら，密行列として持てば $1.7\times10^{10}$ 要素(約 $270$ GB)だが，疎行列なら $9\times10^{7}$ 要素(約 $1.4$ GB)で済む．約200倍の差である．ただし注意:$H$ と $S$ が疎でも，\eqref{eq:14-gevp} を直接対角化して得る固有ベクトルは一般に密であり，直接対角化のコストは $O(M^{3})$ のままである．$O(N)$ を最後まで貫くには，固有ベクトルを経由せず密度行列を直接求める別のアルゴリズムが必要になる．疎性はその前提条件である．

14.3.2 なぜ「閉じ込める」のか — 擬原子軌道の生成

\eqref{eq:14-sparsity} は「基底が有限半径で厳密にゼロになる」ことを前提にしていた．ところが自由原子の軌道はそうではない．束縛状態の動径部分は遠方で

\begin{equation} u_{nl}(r)\ \propto\ e^{-\sqrt{2\abs{\varepsilon_{nl}}}\,r} \qquad(r\to\infty) \label{eq:14-free-tail} \end{equation}

と指数関数的に減衰する(この漸近形は14.6節で動径方程式から導く)．価電子準位 $\varepsilon\simeq-0.5$ Ha なら減衰長は $1$ Bohr であり，$r=15$ Bohr でも $e^{-15}\approx3\times10^{-7}$ でゼロではない．裾は決して切れない．

では自由原子の軌道を $r_c$ でぶつ切りにすればよいかというと，それは二重の意味でまずい．第一に，切断点で関数が不連続に飛ぶと，その2階微分がデルタ関数の微分を含み，運動エネルギー行列要素 $\braket{\chi_\mu|-\tfrac{1}{2}\nabla^{2}|\chi_\nu}$ が発散的に大きな誤差を持つ．第二に，原子が動いて相手の切断面をまたぐたびに行列要素が飛ぶので，全エネルギーが原子座標の連続関数でなくなり，力が定義できなくなる．

正しい処方は「切る」のではなく「はじめから短距離な関数が出てくるように原子問題を変える」ことである．すなわち，原子のKS方程式を解くときに核ポテンシャルを

\begin{equation} V_{\mathrm{core}}(r)= \begin{cases} -\dfrac{Z}{r} & (r\le r_1)\\[4pt] \displaystyle\sum_{n=0}^{3}b_n r^{n} & (r_1\lt r\le r_c)\\[4pt] h & (r_c\lt r) \end{cases} \label{eq:14-confine} \end{equation}

という閉じ込めポテンシャルに置き換えて解く．3次多項式を挟む理由は接続条件の数え上げから分かる．$V_{\mathrm{core}}$ が $r_1$ と $r_c$ で値と1階微分の両方について連続であることを要求すると

\begin{equation} \sum_n b_n r_1^{n}=-\frac{Z}{r_1},\quad \sum_n n\,b_n r_1^{n-1}=\frac{Z}{r_1^{2}},\quad \sum_n b_n r_c^{n}=h,\quad \sum_n n\,b_n r_c^{n-1}=0 \label{eq:14-confine-cond} \end{equation}

の4本の線形方程式が立ち，未知数もちょうど4個($b_0,\dots,b_3$)なので一意に決まる．多項式の次数はこの数え上げが決めている．$C^{1}$ 級を要求するのは，ポテンシャルに折れ目があると波動関数の3階微分が飛び，高次の数値積分公式の精度が落ちるためである．

この井戸の中で原子のKS方程式を自己無撞着に解くと何が起きるか．$r\gt r_c$ では $v_{\mathrm{eff}}\approx h$(定数)なので，14.6節で導く動径方程式 \eqref{eq:14-radial-u} の遠心力項を無視すれば

\begin{align} \frac{\dd^{2}u}{\dd r^{2}}&=2\left[h-\varepsilon\right]u \qquad(r\gt r_c) \label{eq:14-confine-decay1}\\ \Longrightarrow\quad u(r)&\propto e^{-\sqrt{2(h-\varepsilon)}\,r} \label{eq:14-confine-decay2} \end{align}

となる．1行目は \eqref{eq:14-radial-u} で $v_{\mathrm{eff}}\to h$ と置いただけ，2行目は定数係数の2階線形常微分方程式の減衰解を取ったものである(発散解は境界条件 $u(\infty)=0$ で排除される)．ここで $h-\varepsilon$ を数 Hartree にとれば減衰長 $1/\sqrt{2(h-\varepsilon)}$ は $0.5$ Bohr 以下になり，$r_c$ を1〜2 Bohr 超えた時点で振幅は $10^{-2}$〜$10^{-4}$ に落ちる．実装では $r\gt r_c$ で基底を厳密にゼロと定義するが，その打ち切りが持ち込む誤差は上のように指数関数的に小さい．閉じ込めは，疎行列性という数値的要請を，変分空間の定義の側に押し込む工夫である．

(a) 閉じ込めポテンシャル V_core(r) のグラフ．r1 までは裸の核ポテンシャル −Z/r（藍の実線，灰色の破線は自由原子の −Z/r），r1 から rc までは3次多項式で滑らかに持ち上がり，rc の外では定数 h となる． — 図14.3 (a) 閉じ込めポテンシャル \eqref{eq:14-confine}．$r_1$ までは裸の(擬)核ポテンシャル，$r_1$ から $r_c$ までを3次多項式で滑らかに持ち上げ，$r_c$ の外では定数 $h$ とする．(b) この井戸の中で得られる同じ角運動量 $l$ の動径関数．節の数 $0,1,2,\dots$ の順に並び，いずれも $r_c$ で実質的にゼロになる．これらを multiple-$\zeta$ 基底として使う．

物理的意味:$r_c$ と基底数は「変分パラメータ」として扱える

閉じ込め半径 $r_c$ を大きくすれば，基底関数はより広い空間を張るので全エネルギーは下がる．同じ $r_c$ のもとで動径関数の本数を増やす場合は，追加前の関数がそのまま残るので変分空間は入れ子であり，14.2節の定理(2)が厳密に適用でき，エネルギーは単調に下がる．一方 $r_c$ を変えると閉じ込めポテンシャル自体が変わり，関数の形も変わるので部分空間は入れ子にならない．したがって $r_c$ に関する単調性は定理の保証を受けない．それでも実際には $r_c$ を伸ばすとエネルギーは単調に低下することが経験的に知られており，「$r_c$ と基底数の2つを変分パラメータとみなして収束を確認する」という運用が可能になる．ただし $r_c$ を伸ばしすぎると 14.2.3 節で述べた $S$ の悪条件化が始まるので，下げ幅と安定性のトレードオフになる．

14.3.3 multiple-$\zeta$ と分極軌道

原子あたり何本の関数を用意すべきか．最小限の選択は，価電子の各 $(n,l)$ に動径関数を1本ずつ割り当てるsingle-$\zeta$(SZ)基底である．しかしこれでは，結合を作ったときに原子の電荷状態や電子雲の広がりが自由原子と変わることを表現できない．必要なのは「動径方向の伸び縮み」を表す自由度である．

そのために何本必要かは，次の議論で分かる．Slater型の $1s$ 動径関数

\begin{equation} R(r;\zeta)=2\zeta^{3/2}e^{-\zeta r} \label{eq:14-slater} \end{equation}

を考え，指数 $\zeta$(軌道の「締まり具合」)が環境によって $\zeta_0$ から少しずれるとしよう．$\zeta_0$ のまわりでTaylor展開すると

\begin{align} R(r;\zeta)&=R(r;\zeta_0)+(\zeta-\zeta_0)\left.\frac{\partial R}{\partial \zeta}\right|_{\zeta_0}+O\!\left((\zeta-\zeta_0)^{2}\right) \label{eq:14-zeta-taylor} \end{align}

である．この微分を実行すると，積の微分則により

\begin{align} \frac{\partial R}{\partial \zeta} &=2\cdot\frac{3}{2}\zeta^{1/2}e^{-\zeta r}+2\zeta^{3/2}\cdot(-r)e^{-\zeta r} \label{eq:14-zeta-deriv1}\\ &=\left(\frac{3}{2\zeta}-r\right)R(r;\zeta) \label{eq:14-zeta-deriv2} \end{align}

となる．1行目は $\zeta^{3/2}$ と $e^{-\zeta r}$ をそれぞれ微分したもの，2行目は共通因子 $2\zeta^{3/2}e^{-\zeta r}=R$ をくくり出したものである．ここが要点である:$\partial R/\partial\zeta$ は $R$ と $rR$ の線形結合であり，$rR$ は $R$ と直交する方向(節を1つ持つ関数)を含む．つまり

物理的意味:double-$\zeta$ は「軌道の伸縮」を1次まで表現する基底である

1本の関数 $R(r;\zeta_0)$ しか持たなければ，係数 $c$ を変えても $cR(r;\zeta_0)$ という「振幅のスケール」しか作れず，軌道の広がりは自由原子のまま固定される．2本目として $\partial R/\partial\zeta$(実質的に節を1つ持つ関数)を加えると，線形結合 $c_1R+c_2\,\partial_\zeta R$ が $\zeta$ の変化を1次まで再現できる．これが double-$\zeta$(DZ)基底の変分的な意味である．実装では $\partial_\zeta R$ をそのまま使うのではなく，図14.3(b)のように同じ閉じ込め井戸の中の節数 $0,1,2,\dots$ の固有関数を順に採用する．これらはSturm–Liouville問題の固有関数なので互いに直交しており，$S$ の条件数の観点でも有利である．$N_\zeta$ 本使えば「$\zeta$ の変化を $N_\zeta-1$ 次まで」表現できることになる．

もう一つ必要なのが分極軌道である．自由原子で占有されていない角運動量($\mathrm{Si}$ の $d$，$\mathrm{H}$ の $p$ など)の関数を加えるのだが，なぜそれが必要かは1次摂動論で明快に分かる．原子が非球対称な環境(隣の原子が作る電場)に置かれたとする．摂動を一様電場 $\mathcal{E}$ による $\hat{V}'=\mathcal{E}z$ とすると，$s$ 軌道 $\ket{s}$ の1次補正は

\begin{equation} \ket{\delta s}=\sum_{k\neq s}\frac{\braket{k|\mathcal{E}z|s}}{\varepsilon_s-\varepsilon_k}\ket{k} \label{eq:14-pol-pt} \end{equation}

である(量子力学の標準的な1次摂動論)．ここで $z=r\cos\theta$ であり，球面調和関数を使えば $\cos\theta=\sqrt{4\pi/3}\,Y_{10}(\hat{\rr})$ と書ける．したがって行列要素は

\begin{equation} \braket{k|z|s} \propto \int \dd\Omega\ Y_{l_k m_k}^{*}(\hat\rr)\,Y_{10}(\hat\rr)\,Y_{00}(\hat\rr) \label{eq:14-pol-angular} \end{equation}

という角度積分を含む．$Y_{00}$ は定数なので，これは $Y_{l_km_k}^{*}$ と $Y_{10}$ の直交性の判定に帰着し，$l_k=1,m_k=0$ 以外ではゼロになる．結論:$s$ 軌道を電場方向に歪ませるには $p$ 関数が必要であり，$p$ 軌道を歪ませるには $s$ と $d$ が要る．一般に「占有軌道の最大 $l$ より1つ大きい $l$ の関数を1組加える」という慣用則は，この選択則から来ている．

表14.3 LCPAO基底の表記と関数の本数(1原子あたり)
表記	意味	本数の内訳	合計	使いどころ
s1p1	SZ(最小基底)	$1\times1+1\times3$	4	予備計算・定性的傾向のみ
s2p2	DZ	$2\times1+2\times3$	8	動径の伸縮は表現できるが分極が不足
s2p2d1	DZP(標準)	$2\times1+2\times3+1\times5$	13	第2周期・第3周期元素の実用標準
s3p3d2f1	TZDP	$3+9+10+7$	29	高精度が要るとき・遷移金属

$l$ の関数1本は $2l+1$ 個の磁気量子数を伴うので，$s$ は1個，$p$ は3個，$d$ は5個，$f$ は7個と数えることに注意する．第2周期元素なら s2p2d1 の13関数が実用上の標準であり，平面波法の原子あたり数百個と比べて一桁以上少ない．

14.3.4 基底関数そのものを変分最適化する

図14.3(b)の関数群をそのまま使うのではなく，それらを「原始基底」とみなして少数本に縮約し，縮約係数を変分的に最適化する，という考え方がある．すなわち原始動径関数 $\{R_{l\zeta}(r)\}_{\zeta=1}^{N_\zeta}$ から

\begin{equation} \phi_{l\mu}(r)=\sum_{\zeta=1}^{N_\zeta}a_{l\mu\zeta}\,R_{l\zeta}(r), \qquad \mu=1,\dots,N_\mu\ (\ll N_\zeta) \label{eq:14-contraction} \end{equation}

を作り，係数 $a_{l\mu\zeta}$ を「代表的な参照系(二量体や結晶)の全エネルギーが最小になるように」決める．これが実行可能なのは，$E_{\mathrm{tot}}$ の $a$ に関する勾配が解析的に書けるからである．それを導こう．

導出:基底パラメータに関する全エネルギーの勾配

基底関数に含まれる任意のパラメータを $\lambda$ とする($\lambda$ は縮約係数でも，原子位置 $\bm{\tau}_I$ でもよい)．$\lambda$ は $H$ と $S$ の両方に入る．まず1本の軌道エネルギー(Rayleigh商)

$$ \varepsilon_i=\frac{\bm{c}_i^{\dagger}H\bm{c}_i}{\bm{c}_i^{\dagger}S\bm{c}_i} \label{eq:14-rayleigh-again} $$

を $\lambda$ で微分する．$\varepsilon_i$ は $\lambda$ に陽に依存する($H$，$S$ を通じて)ほか，$\bm{c}_i$ を通じて陰にも依存するので，全微分は

$$ \frac{\dd\varepsilon_i}{\dd\lambda} =\underbrace{\frac{\partial\varepsilon_i}{\partial\lambda}\bigg|_{\bm{c}_i}}_{\text{陽}} +\underbrace{\sum_\mu\left(\frac{\partial\varepsilon_i}{\partial c_{\mu i}}\frac{\dd c_{\mu i}}{\dd\lambda} +\frac{\partial\varepsilon_i}{\partial c_{\mu i}^{*}}\frac{\dd c_{\mu i}^{*}}{\dd\lambda}\right)}_{\text{陰}} \label{eq:14-total-deriv} $$

となる．ここで陰の部分は消える．なぜなら $\bm{c}_i$ は \eqref{eq:14-gevp} の解，すなわちRayleigh商を停留にするベクトルであり，定義により $\partial\varepsilon_i/\partial c_{\mu i}^{*}=0$ だからである．これが一般化固有値問題版のHellmann–Feynman定理である．残る陽の部分を，$\bm{c}_i^{\dagger}S\bm{c}_i=1$ と規格化したうえで商の微分則 $\dd(A/B)=\dd A/B-(A/B^{2})\dd B$ に従って計算すると

$$ \frac{\dd\varepsilon_i}{\dd\lambda} =\bm{c}_i^{\dagger}\frac{\partial H}{\partial\lambda}\bm{c}_i -\varepsilon_i\,\bm{c}_i^{\dagger}\frac{\partial S}{\partial\lambda}\bm{c}_i \label{eq:14-gen-hf} $$

を得る($B=1$，$A/B=\varepsilon_i$ を使った)．占有数 $f_i$ を掛けて和をとり，密度行列とエネルギー密度行列

\begin{equation} \rho_{\nu\mu}=\sum_i f_i\,c_{\nu i}c_{\mu i}^{*}, \qquad e_{\nu\mu}=\sum_i f_i\,\varepsilon_i\,c_{\nu i}c_{\mu i}^{*} \label{eq:14-densmat} \end{equation}

を導入すれば

$$ \frac{\dd}{\dd\lambda}\sum_i f_i\varepsilon_i =\sum_{\mu\nu}\left[\rho_{\nu\mu}\frac{\partial H_{\mu\nu}}{\partial\lambda} -e_{\nu\mu}\frac{\partial S_{\mu\nu}}{\partial\lambda}\right] \label{eq:14-basis-grad} $$

となる．第1項が素朴に予想される寄与，第2項が「基底が非直交で，しかも $\lambda$ とともに動く」ことに由来する補正である．$\lambda=\bm{\tau}_I$ と取れば，この第2項が 14.1.2 節で予告したPulay力(重なり項の寄与)にほかならない．平面波基底では $S$ が単位行列で $\lambda$ に依らないため第2項が消え，Hellmann–Feynman項だけが残る．

あとは $\partial H/\partial a$，$\partial S/\partial a$ が要る．\eqref{eq:14-contraction} を \eqref{eq:14-hs-def} に代入すると，行列要素は縮約係数について双線形であることが分かる．重なり行列で書けば

\begin{align} S_{(I\mu)(J\nu)} &=\int\dd^3r\,\Big(\sum_\zeta a^{*}_{I\mu\zeta}R_{I\zeta}\Big)^{*}\Big(\sum_{\zeta'}a_{J\nu\zeta'}R_{J\zeta'}\Big)\ \ (\text{角度部分は略記}) \label{eq:14-contr-s1}\\ &=\sum_{\zeta\zeta'}a_{I\mu\zeta}\,a_{J\nu\zeta'}\;S^{(0)}_{(I\zeta)(J\zeta')} \label{eq:14-contr-s2} \end{align}

となる．ここで $S^{(0)}$ は原始基底同士の重なり(縮約係数を含まない)であり，$a$ を実数にとった．したがって

\begin{equation} \frac{\partial S_{(I\mu)(J\nu)}}{\partial a_{I\mu\zeta}}=\sum_{\zeta'}a_{J\nu\zeta'}S^{(0)}_{(I\zeta)(J\zeta')} \label{eq:14-contr-grad} \end{equation}

という単純な式になり，$H$ についても同様である．原始基底の行列 $S^{(0)},H^{(0)}$ を一度計算しておけば，勾配は行列の掛け算だけで得られる．あとは第16章で扱う準Newton法をそのまま流用して $a$ を最適化すればよい．

注意:最適化基底の「移植性」

この最適化は，ある参照系(たとえば二量体，あるいは代表的な結晶)の全エネルギーを下げるように行われる．したがって得られた基底は，その参照系の化学環境に最適化されている．実際には，複数の参照系のエネルギーの重み付き和を最小化することで移植性(transferability)を確保する．ここが局在基底の難しさであり，同時に面白さでもある．平面波の $E_{\mathrm{cut}}$ のように「1つ上げれば必ず良くなる」パラメータが存在しない代わりに，系の物理を知っていれば少数の関数で高精度に到達できる．基底ライブラリを配布する際には，後述するΔゲージ(14.7節)のような客観的な精度指標で検証することが不可欠になる．

以上で「どんな関数で展開するか」が決まった．次節では，その基底の上で全エネルギーをどう数値的に評価するかに進む．ここには，局在基底実装の成否を分ける核心的な工夫が待っている．

14.4 全エネルギーの実装形 — クーロン項の再編成

第10章で導いた全エネルギーの表式は，数学的には完結している．しかしそれをそのまま計算機に書き下すと，まともな精度が出ない．原因は長距離クーロン力である．この節では，まず何が起きるのかを逆空間で正確に見(14.4.1)，次に「中性原子ポテンシャル+差電荷」という再編成によって発散も桁落ちも同時に消える様子を，恒等変形として1行ずつ追う(14.4.2〜14.4.3)．最後に差電荷のPoisson方程式をFFTで解く手順を示す(14.4.4)．

14.4.1 素朴な表式と，その数値的な破綻

擬ポテンシャルを使う実装では，全エネルギーは次の5項に分かれる(第10章で導いたKS全エネルギー表式に，擬ポテンシャルの非局所項と核間反発を書き加えたものである)．

\begin{equation} E_{\mathrm{tot}}=E_{\mathrm{kin}}+E_{\mathrm{ec}}+E_{\mathrm{ee}}+E_{xc}+E_{\mathrm{cc}} \label{eq:14-etot-5} \end{equation}

各項の定義は

\begin{align} E_{\mathrm{kin}}&=\sum_{\mu\nu}\rho_{\nu\mu}\int\dd^3r\,\chi_\mu^{*}(\rr)\left(-\tfrac{1}{2}\nabla^{2}\right)\chi_\nu(\rr), \label{eq:14-ekin}\\ E_{\mathrm{ec}}&=\underbrace{\int\dd^3r\;n(\rr)\sum_I V_{\mathrm{core},I}(\rr-\bm{\tau}_I)}_{\displaystyle E_{\mathrm{ec}}^{(\mathrm{L})}} +\underbrace{\sum_{\mu\nu}\rho_{\nu\mu}\braket{\chi_\mu|\hat V_{\mathrm{NL}}|\chi_\nu}}_{\displaystyle E_{\mathrm{ec}}^{(\mathrm{NL})}}, \label{eq:14-eec}\\ E_{\mathrm{ee}}&=\frac{1}{2}\iint\dd^3r\,\dd^3r'\,\frac{n(\rr)\,n(\rr')}{\abs{\rr-\rr'}}=\frac{1}{2}\int\dd^3r\;n(\rr)\,v_{\mathrm{H}}[n](\rr), \label{eq:14-eee}\\ E_{xc}&=\int\dd^3r\;n(\rr)\,\epsilon_{xc}\big(n(\rr)\big), \qquad E_{\mathrm{cc}}=\frac{1}{2}\sum_{I\neq J}\frac{Z_IZ_J}{\abs{\bm{\tau}_I-\bm{\tau}_J}} \label{eq:14-exc-ecc} \end{align}

である．$V_{\mathrm{core},I}$ は擬ポテンシャルの局所部分($r\to\infty$ で $-Z_I/r$，$Z_I$ は価電子数)，$\hat V_{\mathrm{NL}}$ は非局所部分(第15章)，$\rho_{\nu\mu}$ は \eqref{eq:14-densmat} の密度行列である．

ここで問題になるのは $E_{\mathrm{ec}}^{(\mathrm{L})}$，$E_{\mathrm{ee}}$，$E_{\mathrm{cc}}$ の3つ，すなわち $1/r$ の長距離テールを持つ項である．周期系ではこの3つは個別に発散する．逆空間で見ると理由が一目で分かる．

導出:三つのクーロン項の $\bm{G}=0$ 発散とその相殺

周期セル $\Omega$ でのFourier変換の規約を

\begin{equation} \tilde f(\bm{G})=\int_{\Omega}\dd^3r\,f(\rr)\,e^{-i\bm{G}\cdot\rr}, \qquad f(\rr)=\frac{1}{\Omega}\sum_{\bm{G}}\tilde f(\bm{G})\,e^{i\bm{G}\cdot\rr} \label{eq:14-ft-conv} \end{equation}

と定める．$2\pi$ をすべて逆格子ベクトルの定義側($\bm{b}_i\cdot\bm{a}_j=2\pi\delta_{ij}$，第12章)に押し込み，変換の前に $(2\pi)^{-3/2}$ のような因子を付けない，という点は第12章と共通である．ただしセル体積 $\Omega$ の置き場所だけが違うことに注意してほしい:第12章の $V_{\bm{G}}=\Omega^{-1}\int_\Omega V e^{-i\bm{G}\cdot\rr}\dd^3r$(および \eqref{eq:14-pw-hmat} の $\tilde v_{\mathrm{eff}}$)は順変換の側で $\Omega$ で割る規約であり，本節の $\tilde f$ とは $\tilde f(\bm{G})=\Omega\,f_{\bm{G}}$ の関係にある．以下 \eqref{eq:14-ft-conv} を使う箇所では $\Omega$ の因子がこの規約で入っている．Poisson方程式 $\nabla^{2}v_{\mathrm{H}}=-4\pi n$ の両辺をFourier変換すると($\bm{G}$ について1次同次なので，この式自体はどちらの規約でも同形である)，$\nabla^{2}\to-\abs{\bm{G}}^{2}$ より

\begin{equation} \tilde v_{\mathrm{H}}(\bm{G})=\frac{4\pi\,\tilde n(\bm{G})}{\abs{\bm{G}}^{2}} \label{eq:14-poisson-g} \end{equation}

である．Parsevalの関係 $\int_\Omega f^{*}g\,\dd^3r=\Omega^{-1}\sum_{\bm{G}}\tilde f^{*}(\bm{G})\tilde g(\bm{G})$ を使うと，Hartreeエネルギーは

$$ E_{\mathrm{ee}}=\frac{1}{2}\int_\Omega n\,v_{\mathrm{H}} =\frac{1}{2\Omega}\sum_{\bm{G}}\tilde n^{*}(\bm{G})\tilde v_{\mathrm{H}}(\bm{G}) =\frac{2\pi}{\Omega}\sum_{\bm{G}}\frac{\abs{\tilde n(\bm{G})}^{2}}{\abs{\bm{G}}^{2}} \label{eq:14-eee-g} $$

となる．$\bm{G}=0$ 項は $\tilde n(0)=\int_\Omega n\,\dd^3r=N_{\mathrm{e}}$ なので $+2\pi N_{\mathrm{e}}^{2}/(\Omega\abs{\bm{G}}^{2})\to+\infty$ である．同様に，$V_{\mathrm{core},I}$ の長距離部分 $-Z_I/\abs{\rr-\bm{\tau}_I}$ のFourier変換は $-4\pi Z_Ie^{-i\bm{G}\cdot\bm{\tau}_I}/\abs{\bm{G}}^{2}$ だから

$$ E_{\mathrm{ec}}^{(\mathrm{L})}\Big|_{\bm{G}=0} =-\frac{4\pi\,N_{\mathrm{e}}\,Z_{\mathrm{tot}}}{\Omega\abs{\bm{G}}^{2}}\to-\infty, \qquad Z_{\mathrm{tot}}\equiv\sum_I Z_I \label{eq:14-eec-g0} $$

であり，核間反発をEwald和で書いたときの $\bm{G}=0$ 項は $+2\pi Z_{\mathrm{tot}}^{2}/(\Omega\abs{\bm{G}}^{2})\to+\infty$ である．3つを足すと

$$ \frac{2\pi}{\Omega\abs{\bm{G}}^{2}}\Big[N_{\mathrm{e}}^{2}-2N_{\mathrm{e}}Z_{\mathrm{tot}}+Z_{\mathrm{tot}}^{2}\Big] =\frac{2\pi}{\Omega\abs{\bm{G}}^{2}}\big(N_{\mathrm{e}}-Z_{\mathrm{tot}}\big)^{2}=0 \label{eq:14-g0-cancel} $$

となる．最後の等号は電気的中性条件 $N_{\mathrm{e}}=Z_{\mathrm{tot}}$ による．すなわち3つの発散は完全に相殺するが，それは「3つを足した後」の話である．

この相殺を手作業で仕組むのが古典的なEwald和である．しかしそれで話が終わるわけではない．発散を取り除いても，残った有限の3項はそれぞれ原子あたり数 Hartree〜数十 Hartree の大きさを持ち，符号が交互で，和は原子あたり数 Hartree に落ちる．ところが我々が知りたい凝集エネルギーや結合エネルギーの差は $10^{-2}$〜$10^{-3}$ Ha である．大きな数の差として小さな数を得る計算は，各項に相対誤差 $10^{-6}$ 以下を要求する．実空間格子上の数値積分で各項を独立に $10^{-6}$ の相対精度で評価するのは，格子を極端に細かくしない限り不可能である．これが「クーロン項の再編成が数値実装の核心である」と言われる理由である．

14.4.2 中性原子ポテンシャルと差電荷

発想の出発点は「原子はもともと中性である」という当たり前の事実である．孤立した(閉じ込め付きの)原子 $I$ の球対称な価電子密度を $n_I^{(a)}(\rr-\bm{\tau}_I)$ と書く．中性なので

\begin{equation} \int\dd^3r\;n_I^{(a)}(\rr)=Z_I \label{eq:14-atomic-neutral} \end{equation}

である．これらを重ね合わせたものを $n^{(a)}\equiv\sum_I n_I^{(a)}$，実際の密度との差を

\begin{equation} \delta n(\rr)\equiv n(\rr)-\sum_I n_I^{(a)}(\rr-\bm{\tau}_I) \label{eq:14-dn-def} \end{equation}

と定義する．これが差電荷密度である．結合が作られると電子は原子間に少し移動するが，その移動量は全電子数に比べれば小さい．したがって $\delta n$ は $n$ よりずっと小さく，しかも重要なことに

\begin{equation} \int_\Omega\dd^3r\;\delta n(\rr)=N_{\mathrm{e}}-\sum_I Z_I=0 \qquad\Longleftrightarrow\qquad \widetilde{\delta n}(\bm{G}=0)=0 \label{eq:14-dn-neutral} \end{equation}

である．差電荷はそれ自体が電気的に中性なのである．これが後で効いてくる．

次に，原子 $I$ の価電子密度が作るHartreeポテンシャルを $v_{\mathrm{H},I}^{(a)}$ と書き(すなわち $\nabla^{2}v_{\mathrm{H},I}^{(a)}=-4\pi n_I^{(a)}$)，擬ポテンシャルの局所部分と足し合わせた

\begin{equation} V_{\mathrm{na},I}(\rr)\equiv V_{\mathrm{core},I}(\rr)+v_{\mathrm{H},I}^{(a)}(\rr) \label{eq:14-vna-def} \end{equation}

を中性原子ポテンシャル(neutral atom potential)と呼ぶ．これは「原子核(イオン芯)とその周りの電子雲を一緒にした，電気的に中性な原子が作るポテンシャル」である．中性なのだから遠方では消えるはずだ，というのが直観だが，実は「遠方で消える」どころか有限距離で厳密にゼロになる．それを保証するのが次の定理である．

数学ノート:球殻定理(球対称電荷分布の外部ポテンシャル)

球対称な電荷分布 $\rho(r)$ が作る静電ポテンシャル $v(r)$ を考える．対称性から $v$ も $r$ だけの関数であり，Poisson方程式 $\nabla^{2}v=-4\pi\rho$ の球対称部分は

$$ \frac{1}{r^{2}}\frac{\dd}{\dd r}\left(r^{2}\frac{\dd v}{\dd r}\right)=-4\pi\rho(r) \label{eq:14-radial-poisson} $$

である(球座標のラプラシアンについては14.6節の数学ノートを参照)．両辺に $r^{2}$ を掛けて $0$ から $r$ まで積分すると，左辺は完全微分なので

$$ \left[r'^{2}\frac{\dd v}{\dd r'}\right]_{0}^{r} =r^{2}\frac{\dd v}{\dd r} =-4\pi\int_{0}^{r}\rho(r')\,r'^{2}\dd r' \equiv -Q(r) \label{eq:14-gauss} $$

となる($v$ が原点で正則なら $r'^{2}v'\to0$，また $Q(r)$ は半径 $r$ の球内の全電荷である)．すなわち

\begin{equation} \frac{\dd v}{\dd r}=-\frac{Q(r)}{r^{2}}. \label{eq:14-dvdr} \end{equation}

いま $\rho$ が半径 $R$ の外でゼロならば，$r\gt R$ では $Q(r)=Q_{\mathrm{tot}}$(定数)なので $\dd v/\dd r=-Q_{\mathrm{tot}}/r^{2}$ であり，$v(\infty)=0$ の条件で積分すれば

\begin{equation} v(r)=\frac{Q_{\mathrm{tot}}}{r}\qquad(r\gt R). \label{eq:14-shell} \end{equation}

球対称分布は，その外側では中心に置いた点電荷と全く同じポテンシャルを作る．これがNewtonの球殻定理である．

この定理を \eqref{eq:14-vna-def} に適用しよう．$n_I^{(a)}$ は球対称で，閉じ込めにより半径 $r_c^{I}$ の外でゼロだから，$r\gt r_c^{I}$ では $v_{\mathrm{H},I}^{(a)}(r)=+Z_I/r$ である(\eqref{eq:14-atomic-neutral} により全電荷は $Z_I$)．一方 $V_{\mathrm{core},I}(r)=-Z_I/r$ は擬ポテンシャルの芯半径の外で厳密に成り立つ．両者を足すと

\begin{equation} V_{\mathrm{na},I}(r)=-\frac{Z_I}{r}+\frac{Z_I}{r}=0 \qquad (r\gt r_c^{I}) \label{eq:14-vna-short} \end{equation}

となる．「遠方で速く減衰する」のではなく，有限半径の外で恒等的にゼロである．したがって $V_{\mathrm{na},I}$ を含む積分はすべて短距離の二中心積分になり，逆空間で高精度に評価できる．ここに全ての鍵がある．

図14.4 中性原子分解の考え方．(a) 擬ポテンシャル局所部分 $V_{\mathrm{core}}$ の $-Z/r$ テールと，原子価電子密度が作る $v_{\mathrm{H}}^{(a)}$ の $+Z/r$ テールが厳密に相殺し，中性原子ポテンシャル $V_{\mathrm{na}}$ は $r_c$ の外で恒等的にゼロになる．(b) 実際の密度 $n$ と原子密度の重ね合わせ $\sum_I n_I^{(a)}$ の差 $\delta n$ は小さく，かつ全体で中性である．

14.4.3 恒等変形:三つの発散項を三つの短距離項へ

準備が整った．これから \eqref{eq:14-eec}〜\eqref{eq:14-exc-ecc} の3つの長距離項を，近似を一切使わずに書き換える．使う道具は2つだけである:(i) Poisson方程式の線形性，(ii) クーロン核の対称性 $\int n_1v_{\mathrm{H}}[n_2]=\int n_2v_{\mathrm{H}}[n_1]$．後者は

\begin{equation} \int\dd^3r\,n_1(\rr)\,v_{\mathrm{H}}[n_2](\rr) =\iint\dd^3r\,\dd^3r'\,\frac{n_1(\rr)n_2(\rr')}{\abs{\rr-\rr'}} \label{eq:14-kernel-sym} \end{equation}

が $n_1\leftrightarrow n_2$ の入れ替えについて対称であること($\abs{\rr-\rr'}$ が $\rr\leftrightarrow\rr'$ で不変で，積分変数の名前を付け替えればよい)から従う．

導出:$E_{\mathrm{ec}}^{(\mathrm{L})}+E_{\mathrm{ee}}+E_{\mathrm{cc}}=E_{\mathrm{na}}+E_{\delta\mathrm{ee}}+E_{\mathrm{scc}}$

ステップ1(記号の準備)．次の3つを定義する．

$$ V_{\mathrm{na}}(\rr)\equiv\sum_I V_{\mathrm{na},I}(\rr-\bm{\tau}_I), \qquad v_{\mathrm{H}}^{(a)}\equiv\sum_I v_{\mathrm{H},I}^{(a)}, \qquad \delta v_{\mathrm{H}}\equiv v_{\mathrm{H}}[\delta n]. \label{eq:14-defs-step1} $$

Poisson方程式は線形なので，$n=n^{(a)}+\delta n$ に対して

\begin{equation} v_{\mathrm{H}}[n]=v_{\mathrm{H}}[n^{(a)}]+v_{\mathrm{H}}[\delta n]=v_{\mathrm{H}}^{(a)}+\delta v_{\mathrm{H}} \label{eq:14-vh-split} \end{equation}

が厳密に成り立つ(線形方程式の解の重ね合わせ)．

ステップ2($E_{\mathrm{ec}}^{(\mathrm{L})}$ の書き換え)．定義 \eqref{eq:14-vna-def} より $V_{\mathrm{core}}=V_{\mathrm{na}}-v_{\mathrm{H}}^{(a)}$ である．これを $E_{\mathrm{ec}}^{(\mathrm{L})}$ に代入すると

\begin{equation} E_{\mathrm{ec}}^{(\mathrm{L})} =\int n\,V_{\mathrm{core}} =\int n\,V_{\mathrm{na}}-\int n\,v_{\mathrm{H}}^{(a)} \equiv E_{\mathrm{na}}-\int n\,v_{\mathrm{H}}^{(a)} \label{eq:14-step2} \end{equation}

となる．ここで中性原子エネルギー $E_{\mathrm{na}}\equiv\int\dd^3r\,n(\rr)V_{\mathrm{na}}(\rr)$ を定義した．$V_{\mathrm{na}}$ は各原子まわりの短距離関数なので，$E_{\mathrm{na}}$ は二中心積分の和として高精度に計算できる．

ステップ3($E_{\mathrm{ee}}$ の書き換え)．\eqref{eq:14-vh-split} を \eqref{eq:14-eee} に代入する．

\begin{equation} E_{\mathrm{ee}}=\frac{1}{2}\int n\left(v_{\mathrm{H}}^{(a)}+\delta v_{\mathrm{H}}\right) =\frac{1}{2}\int n\,v_{\mathrm{H}}^{(a)}+\frac{1}{2}\int n\,\delta v_{\mathrm{H}} \label{eq:14-step3a} \end{equation}

第2項の $n$ をさらに $n^{(a)}+\delta n$ と分解する:

\begin{equation} \frac{1}{2}\int n\,\delta v_{\mathrm{H}} =\frac{1}{2}\int n^{(a)}\delta v_{\mathrm{H}}+\frac{1}{2}\int \delta n\,\delta v_{\mathrm{H}} \label{eq:14-step3b} \end{equation}

ここで核の対称性 \eqref{eq:14-kernel-sym} を使う．$\delta v_{\mathrm{H}}=v_{\mathrm{H}}[\delta n]$，$v_{\mathrm{H}}^{(a)}=v_{\mathrm{H}}[n^{(a)}]$ だから

$$ \int n^{(a)}\,\delta v_{\mathrm{H}}=\int \delta n\,v_{\mathrm{H}}^{(a)} \label{eq:14-step3c} $$

である．これを \eqref{eq:14-step3b} に入れ，\eqref{eq:14-step3a} に戻すと

$$ E_{\mathrm{ee}}=\frac{1}{2}\int n\,v_{\mathrm{H}}^{(a)}+\frac{1}{2}\int \delta n\,v_{\mathrm{H}}^{(a)}+\frac{1}{2}\int\delta n\,\delta v_{\mathrm{H}} \label{eq:14-step3d} $$

を得る．第3項を差電荷Hartreeエネルギー $E_{\delta\mathrm{ee}}\equiv\frac{1}{2}\int\delta n\,\delta v_{\mathrm{H}}$ と名づける．前2項は $\delta n=n-n^{(a)}$ を使って1つにまとめられる:

$$ \frac{1}{2}\int n\,v_{\mathrm{H}}^{(a)}+\frac{1}{2}\int \delta n\,v_{\mathrm{H}}^{(a)} =\frac{1}{2}\int\big(n+n-n^{(a)}\big)v_{\mathrm{H}}^{(a)} =\int n\,v_{\mathrm{H}}^{(a)}-\frac{1}{2}\int n^{(a)}v_{\mathrm{H}}^{(a)} \label{eq:14-step3e} $$

したがって

\begin{equation} E_{\mathrm{ee}}=\int n\,v_{\mathrm{H}}^{(a)}-\frac{1}{2}\int n^{(a)}v_{\mathrm{H}}^{(a)}+E_{\delta\mathrm{ee}}. \label{eq:14-step3f} \end{equation}

ステップ4(足し合わせ)．\eqref{eq:14-step2} と \eqref{eq:14-step3f} を足す．

$$ E_{\mathrm{ec}}^{(\mathrm{L})}+E_{\mathrm{ee}} =E_{\mathrm{na}}\underbrace{-\int n\,v_{\mathrm{H}}^{(a)}+\int n\,v_{\mathrm{H}}^{(a)}}_{=0} -\frac{1}{2}\int n^{(a)}v_{\mathrm{H}}^{(a)}+E_{\delta\mathrm{ee}} \label{eq:14-step4} $$

長距離量 $\int n\,v_{\mathrm{H}}^{(a)}$ が厳密に消えた．これがこの変形の要である．残ったのは

$$ E_{\mathrm{ec}}^{(\mathrm{L})}+E_{\mathrm{ee}} =E_{\mathrm{na}}+E_{\delta\mathrm{ee}}-\frac{1}{2}\int n^{(a)}v_{\mathrm{H}}^{(a)} \label{eq:14-step4b} $$

である．

ステップ5(核間反発との合体)．最後の項を原子ごとに分解する．$n^{(a)}=\sum_In_I^{(a)}$，$v_{\mathrm{H}}^{(a)}=\sum_Jv_{\mathrm{H},J}^{(a)}$ なので

$$ \frac{1}{2}\int n^{(a)}v_{\mathrm{H}}^{(a)} =\frac{1}{2}\sum_{I}\sum_{J}U_{IJ}, \qquad U_{IJ}\equiv\int\dd^3r\;n_I^{(a)}(\rr-\bm{\tau}_I)\,v_{\mathrm{H},J}^{(a)}(\rr-\bm{\tau}_J) \label{eq:14-uij} $$

と書ける．$I=J$ と $I\neq J$ に分けて $E_{\mathrm{cc}}$ と合わせると

\begin{equation} E_{\mathrm{cc}}-\frac{1}{2}\int n^{(a)}v_{\mathrm{H}}^{(a)} =\underbrace{\frac{1}{2}\sum_{I\neq J}\left[\frac{Z_IZ_J}{\abs{\bm{\tau}_I-\bm{\tau}_J}}-U_{IJ}\right]}_{\text{位置に依存する部分}} -\underbrace{\frac{1}{2}\sum_I U_{II}}_{\text{原子ごとの定数}} \equiv E_{\mathrm{scc}} \label{eq:14-escc-def} \end{equation}

を得る．これを遮蔽核間反発エネルギーと呼ぶ．$U_{II}$ は孤立原子の量であり原子配置に依らないので，力にもストレスにも寄与しないが，全エネルギーの絶対値を合わせるためには含めなければならない．

$E_{\mathrm{scc}}$ が短距離であることを確認しよう．原子 $I$ と $J$ の電子雲が重ならない，すなわち $\abs{\bm{\tau}_I-\bm{\tau}_J}\gt r_c^{I}+r_c^{J}$ のとき，球殻定理 \eqref{eq:14-shell} により $v_{\mathrm{H},J}^{(a)}$ は $n_I^{(a)}$ の存在領域では点電荷ポテンシャル $Z_J/\abs{\rr-\bm{\tau}_J}$ に等しい．したがって

\begin{align} U_{IJ}&=\int\dd^3r\;n_I^{(a)}(\rr-\bm{\tau}_I)\,\frac{Z_J}{\abs{\rr-\bm{\tau}_J}} \label{eq:14-uij-far1}\\ &=Z_J\cdot\frac{Z_I}{\abs{\bm{\tau}_I-\bm{\tau}_J}} =\frac{Z_IZ_J}{\abs{\bm{\tau}_I-\bm{\tau}_J}} \label{eq:14-uij-far2} \end{align}

となる．2行目では，今度は $n_I^{(a)}$ の側に球殻定理を適用した(球対称分布 $n_I^{(a)}$ が点 $\bm{\tau}_J$ に作るポテンシャルは $Z_I/\abs{\bm{\tau}_I-\bm{\tau}_J}$ に等しく，$\int n_I^{(a)}v=$ その値に $Z_J$ を掛けたもの)．ゆえに \eqref{eq:14-escc-def} の角括弧は恒等的にゼロになり，$E_{\mathrm{scc}}$ の和は重なる原子対だけに制限される．

物理的意味:Ewald和が要らなくなる仕組み

もとの $E_{\mathrm{cc}}=\frac{1}{2}\sum_{I\neq J}Z_IZ_J/\abs{\bm{\tau}_{IJ}}$ は $1/r$ 和であり，周期系では条件収束(足す順序で値が変わる)する．Ewald和はこれを逆空間と実空間の2つの絶対収束和に分解する技法だった．中性原子分解では，この $1/r$ 和が「点電荷 $Z_IZ_J/r$」から「その原子の電子雲による遮蔽 $U_{IJ}$」を差し引いた形で現れ，両者は電子雲が重ならない距離で厳密に相殺する．すなわち各原子は自分の価電子雲を纏った中性の物体として振る舞い，離れた原子同士は互いに何の力も及ぼさない．残った長距離的な相互作用は差電荷 $\delta n$ が担うが，$\delta n$ は \eqref{eq:14-dn-neutral} により中性なので，その $\bm{G}=0$ 成分はゼロであり，Poisson方程式を逆空間で解いても $1/\abs{\bm{G}}^{2}$ の発散に出会わない．Ewald和が不要になるのではなく，Ewald和がやっていた「中性化」を，物理的に意味のある単位(中性原子)で最初から済ませているのである．

14.4.4 差電荷のPoisson方程式をFFTで解く

残った長距離項は $E_{\delta\mathrm{ee}}=\frac{1}{2}\int\delta n\,\delta v_{\mathrm{H}}$ だけである．これは逆空間で解く．規約 \eqref{eq:14-ft-conv} のもと，$\delta n$ を一様FFT格子上でサンプリングして $\widetilde{\delta n}(\bm{G})$ を得れば，\eqref{eq:14-poisson-g} と同じ手順で

\begin{equation} \widetilde{\delta v}_{\mathrm{H}}(\bm{G})=\frac{4\pi\,\widetilde{\delta n}(\bm{G})}{\abs{\bm{G}}^{2}}\quad(\bm{G}\neq0), \qquad \widetilde{\delta v}_{\mathrm{H}}(0)=0 \label{eq:14-dpoisson} \end{equation}

となる．$\bm{G}=0$ を単に $0$ と置いてよいのは，\eqref{eq:14-dn-neutral} により分子 $\widetilde{\delta n}(0)$ もゼロだからである(通常のPoisson方程式では $\tilde n(0)=N_{\mathrm{e}}\neq0$ なので $0/0$ ではなく $N_{\mathrm{e}}/0$ の真の発散になり，一様な補償背景電荷を導入する必要があった)．エネルギーはParsevalの関係で

\begin{align} E_{\delta\mathrm{ee}} &=\frac{1}{2}\int_\Omega\dd^3r\;\delta n(\rr)\,\delta v_{\mathrm{H}}(\rr) =\frac{1}{2\Omega}\sum_{\bm{G}}\widetilde{\delta n}^{*}(\bm{G})\,\widetilde{\delta v}_{\mathrm{H}}(\bm{G}) \label{eq:14-edee-g1}\\ &=\frac{2\pi}{\Omega}\sum_{\bm{G}\neq0}\frac{\abs{\widetilde{\delta n}(\bm{G})}^{2}}{\abs{\bm{G}}^{2}} \label{eq:14-edee-g2} \end{align}

と評価される(1行目でParseval，2行目で \eqref{eq:14-dpoisson} を代入した．$\delta n$ が実関数なので $\widetilde{\delta n}(-\bm{G})=\widetilde{\delta n}^{*}(\bm{G})$ であり，和は実数になる)．実際の手順は次の3ステップである:格子上の $\delta n(\rr_p)$ を順方向FFT($O(N_g\log N_g)$)，$\bm{G}$ ごとに $4\pi/\abs{\bm{G}}^{2}$ を掛ける($O(N_g)$)，逆方向FFTで $\delta v_{\mathrm{H}}(\rr_p)$ を得る($O(N_g\log N_g)$)．

14.4.5 実装形のまとめ

以上をまとめると，全エネルギーは次の形になる．

\begin{equation} E_{\mathrm{tot}} =E_{\mathrm{kin}} +E_{\mathrm{na}} +E_{\mathrm{ec}}^{(\mathrm{NL})} +E_{\delta\mathrm{ee}} +E_{xc} +E_{\mathrm{scc}} \label{eq:14-etot-final} \end{equation}

もとの \eqref{eq:14-etot-5} と項数は同じだが，中身は決定的に違う．発散する項も，他項との大きな相殺を要する項も，一つも残っていない．それぞれの項の性格と評価法を表14.4にまとめる．

表14.4 全エネルギー各項の性格と数値評価法
項	表式	距離	評価法
$E_{\mathrm{kin}}$	$\sum\rho_{\nu\mu}\braket{\chi_\mu\|-\frac{1}{2}\nabla^{2}\|\chi_\nu}$	短(基底の $r_c$)	二中心積分，逆空間で解析的
$E_{\mathrm{na}}$	$\int n\sum_I V_{\mathrm{na},I}$	短(\eqref{eq:14-vna-short})	二中心・三中心積分，逆空間
$E_{\mathrm{ec}}^{(\mathrm{NL})}$	$\sum\rho_{\nu\mu}\braket{\chi_\mu\|\hat V_{\mathrm{NL}}\|\chi_\nu}$	短(射影子の台)	分離型射影子の二中心積分
$E_{\delta\mathrm{ee}}$	$\frac{1}{2}\int\delta n\,\delta v_{\mathrm{H}}$	長だが小さい	一様FFT格子，\eqref{eq:14-edee-g2}
$E_{xc}$	$\int n\,\epsilon_{xc}(n)$	局所(GGAは準局所)	実空間の一様格子で数値積分
$E_{\mathrm{scc}}$	$\frac{1}{2}\sum_{I\neq J}[Z_IZ_J/\tau_{IJ}-U_{IJ}]+\text{const}$	短(\eqref{eq:14-uij-far2})	重なる原子対のみ，動径細格子

例:初期密度としての $n^{(a)}$

この分解にはおまけがある．SCFの初期密度として何を使うべきかという問いに，自然な答えが用意されているのである．$n^{(0)}=n^{(a)}=\sum_In_I^{(a)}$ と取れば，そのとき $\delta n=0$，$E_{\delta\mathrm{ee}}=0$ から出発することになる．すなわち「原子を並べただけの状態」から出発し，SCFループが $\delta n$ を育てていくという描像になる．図14.1のフローで初期密度を原子密度の重ね合わせとしたのはこの理由による．結合の形成による電荷移動は全電子数の数パーセントに過ぎないので，この初期値は真の解にかなり近く，SCFの反復回数を大きく減らす．

14.5 SCFと電荷混合

14.1節で述べたように，KS方程式は非線形固有値問題であり，その解は写像 $F$ の不動点である．ここではその不動点反復を線形化して収束条件を厳密に求め，金属で悪名高いcharge sloshing(電荷の揺さぶり)がなぜ起きるのかを第8章の誘電遮蔽の言葉で説明する．そのうえで，Kerker混合とPulay(DIIS)混合という2つの標準的処方を，それぞれ「前処理」と「残差最小化」の原理から導出する．

14.5.1 不動点問題と単純混合の収束条件

SCFの1反復は写像

\begin{equation} n_{\mathrm{out}}=F[n_{\mathrm{in}}] \label{eq:14-scf-map} \end{equation}

である(具体的には，$n_{\mathrm{in}}$ から $v_{\mathrm{eff}}$ を作り，\eqref{eq:14-gevp} を解き，\eqref{eq:14-dens} で密度を組み直す一連の操作)．求めたいのは不動点 $n^{*}=F[n^{*}]$ である．反復の $n$ 回目の誤差を

\begin{equation} e^{(n)}(\rr)\equiv n_{\mathrm{in}}^{(n)}(\rr)-n^{*}(\rr) \label{eq:14-err-def} \end{equation}

と定義し，$e$ が小さいとして $F$ を不動点のまわりで線形化する．汎関数のTaylor展開(第4章)により

\begin{align} F[n^{*}+e]&=F[n^{*}]+\int\dd^3r'\left.\frac{\delta F(\rr)}{\delta n(\rr')}\right|_{n^{*}}e(\rr')+O(e^{2}) \label{eq:14-linearize1}\\ &=n^{*}+\big(J e\big)(\rr)+O(e^{2}) \label{eq:14-linearize2} \end{align}

となる．ここで $J$ は写像 $F$ のヤコビ演算子(汎関数微分を積分核とする線形演算子)である．1行目は汎関数のTaylor展開の定義そのもの，2行目では $F[n^{*}]=n^{*}$(不動点の定義)を使い，積分演算子を $J$ と略記した．

最も単純な更新則である単純混合

\begin{equation} n_{\mathrm{in}}^{(n+1)}=(1-\alpha)\,n_{\mathrm{in}}^{(n)}+\alpha\,n_{\mathrm{out}}^{(n)}, \qquad 0\lt\alpha\le1 \label{eq:14-simple-mix} \end{equation}

を採用したときの誤差の伝播を調べよう．両辺から $n^{*}$ を引き，$n_{\mathrm{out}}^{(n)}=F[n_{\mathrm{in}}^{(n)}]=n^{*}+Je^{(n)}$ を代入する:

\begin{align} e^{(n+1)} &=(1-\alpha)\,n_{\mathrm{in}}^{(n)}+\alpha\,n_{\mathrm{out}}^{(n)}-n^{*} \label{eq:14-err-rec1}\\ &=(1-\alpha)\big(n^{*}+e^{(n)}\big)+\alpha\big(n^{*}+Je^{(n)}\big)-n^{*} \label{eq:14-err-rec2}\\ &=\big[(1-\alpha)+\alpha\big]n^{*}-n^{*}+\big[(1-\alpha)\hat{1}+\alpha J\big]e^{(n)} \label{eq:14-err-rec3}\\ &=\big[(1-\alpha)\hat{1}+\alpha J\big]\,e^{(n)} \label{eq:14-err-rec4} \end{align}

1行目は定義，2行目は $n_{\mathrm{in}}^{(n)}=n^{*}+e^{(n)}$ と線形化 \eqref{eq:14-linearize2} の代入，3行目は $n^{*}$ を含む項と $e^{(n)}$ を含む項に整理したもの，4行目は $(1-\alpha)+\alpha=1$ より $n^{*}$ の項が消えることによる．

したがって誤差は反復ごとに行列 $A\equiv(1-\alpha)\hat 1+\alpha J$ を掛けられて伝わる．$n$ 回反復した後の誤差は $e^{(n)}=A^{n}e^{(0)}$ であり，これが $n\to\infty$ でゼロに収束する必要十分条件は，$A$ のスペクトル半径(固有値の絶対値の最大)が1未満であることである．$J$ の固有値を $\lambda_j$ とすると $A$ の固有値は $1-\alpha+\alpha\lambda_j$ だから

\begin{equation} \big|1-\alpha\left(1-\lambda_j\right)\big|\lt1 \qquad\text{(すべての固有値 $\lambda_j$ について)} \label{eq:14-conv-cond} \end{equation}

が単純混合の収束条件である．$\lambda_j$ が実数のときこの不等式を解こう．絶対値を外すと $-1\lt1-\alpha(1-\lambda_j)\lt1$ であり，右側は $\alpha(1-\lambda_j)\gt0$，$\alpha\gt0$ より $\lambda_j\lt1$ を要求する．左側は $\alpha(1-\lambda_j)\lt2$，すなわち

\begin{equation} \alpha\lt\frac{2}{1-\lambda_j} \qquad\Longrightarrow\qquad \alpha\lt\frac{2}{1-\lambda_{\min}} \label{eq:14-alpha-bound} \end{equation}

である($\lambda_{\min}$ は最も負の固有値)．$J$ に大きな負の固有値があると，混合パラメータ $\alpha$ を極端に小さくせざるを得ない．次項では，そのような固有値がどこから来るのかを明らかにする．

14.5.2 ヤコビ演算子の正体 — 誘電遮蔽とcharge sloshing

$J=\delta n_{\mathrm{out}}/\delta n_{\mathrm{in}}$ を物理的に読み解く．写像 $F$ は2段階からなる．第1段:入力密度の変化がポテンシャルを変える．第2段:ポテンシャルの変化が出力密度を変える．

第1段は静電相互作用が支配する(交換相関の寄与は小さいのでここでは無視する)．\eqref{eq:14-poisson-g} より，逆空間では

\begin{equation} \delta v_{\mathrm{eff}}(\bm{q})=\frac{4\pi}{\abs{\bm{q}}^{2}}\,\delta n_{\mathrm{in}}(\bm{q}) \label{eq:14-step1-jac} \end{equation}

である．第2段は，KS方程式を「与えられた外部摂動に対する独立粒子応答」として扱うことに相当し，応答関数 $\chi_0$ を使って

\begin{equation} \delta n_{\mathrm{out}}(\bm{q})=\chi_0(\bm{q})\,\delta v_{\mathrm{eff}}(\bm{q}) \label{eq:14-step2-jac} \end{equation}

と書ける．両者をつなぐと，一様電子ガス的な近似のもとでヤコビ演算子は逆空間で対角化され，その固有値は

\begin{equation} \lambda(\bm{q})=\chi_0(\bm{q})\,\frac{4\pi}{\abs{\bm{q}}^{2}} \label{eq:14-lambda-q} \end{equation}

となる．ここで第8章で導入した誘電関数の定義 $\varepsilon(\bm{q})=1-(4\pi/\abs{\bm{q}}^{2})\chi_0(\bm{q})$ を思い出すと，驚くほど簡単な関係が得られる:

\begin{equation} \lambda(\bm{q})=1-\varepsilon(\bm{q}) \qquad\Longrightarrow\qquad e^{(n+1)}(\bm{q})=\big[1-\alpha\,\varepsilon(\bm{q})\big]\,e^{(n)}(\bm{q}) \label{eq:14-amp-eps} \end{equation}

右側の式は \eqref{eq:14-err-rec4} に $\lambda=1-\varepsilon$ を代入したものである．単純混合の誤差増幅率は $1-\alpha\varepsilon(\bm{q})$ であり，系の誘電関数がそのまま収束の良し悪しを決める．収束条件 \eqref{eq:14-conv-cond} は $0\lt\alpha\varepsilon(\bm{q})\lt2$，すなわち

\begin{equation} \alpha\lt\frac{2}{\varepsilon_{\max}}=\frac{2}{\varepsilon(q_{\min})} \label{eq:14-alpha-eps} \end{equation}

である($\varepsilon$ は $q$ の減少関数なので最大値は最小波数で取る)．周期セルの1辺を $L$ とすれば，セル内で表現できる最小の非ゼロ波数は $q_{\min}=2\pi/L$ である．

金属では第8章のThomas–Fermi遮蔽により

\begin{equation} \varepsilon(q)\simeq1+\frac{\kappa_{\mathrm{TF}}^{2}}{q^{2}}, \qquad \kappa_{\mathrm{TF}}^{2}=4\pi N(E_{\mathrm{F}}) \label{eq:14-tf-eps} \end{equation}

である($N(E_{\mathrm{F}})$ はフェルミ準位での単位体積あたり状態密度)．これを \eqref{eq:14-alpha-eps} に入れると

\begin{equation} \alpha\lt\frac{2}{1+\kappa_{\mathrm{TF}}^{2}L^{2}/(4\pi^{2})} \ \xrightarrow[L\to\infty]{}\ \frac{8\pi^{2}}{\kappa_{\mathrm{TF}}^{2}L^{2}}\ \propto\ \frac{1}{L^{2}} \label{eq:14-alpha-L} \end{equation}

となる．セルを大きくするほど，許される混合パラメータは $L^{-2}$ で小さくなる．これが charge sloshing の正体である．

物理的意味:なぜ長波長で暴れるのか

金属は極めて有効に電荷を遮蔽する．逆に言えば，わずかな電荷の偏りが大きなポテンシャル変化を生む．長い胞の左半分にごく少量の電子が余分に溜まったとしよう．その静電ポテンシャルの変化は $4\pi/q^{2}$ に比例するので，波長が長い($q$ が小さい)ほど巨大になる．次の反復で電子はそのポテンシャルを避けて右半分へ大量に逃げる．さらに次の反復では左へ戻る．系が長いほど振幅は大きくなり，反復は発散する．これが「電荷が胞の中で左右にざぶざぶ揺れる」という比喩の由来である．
一方，絶縁体では $q\to0$ で $\chi_0(q)\propto q^{2}$ となり，$\varepsilon(0)=\varepsilon_\infty$ は有限にとどまる(Si なら約 $12$，多くの分子系では $2$ 程度)．この場合 $\alpha\lesssim2/\varepsilon_\infty$ は $L$ に依らない定数であり，単純混合でも $\alpha=0.2$〜$0.4$ 程度で十分収束する．「絶縁体のSCFは素直に回るのに金属は回らない」という経験則は，誘電関数の $q\to0$ の振る舞いの違いそのものである．

例:単純金属スラブでの収束速度の見積もり

典型的な単純金属($k_{\mathrm{F}}=1.0$ Bohr$^{-1}$)では，自由電子ガスの状態密度 $N(E_{\mathrm{F}})=k_{\mathrm{F}}/\pi^{2}$(両スピン込み，原子単位)から $$\kappa_{\mathrm{TF}}^{2}=4\pi\cdot\frac{k_{\mathrm{F}}}{\pi^{2}}=\frac{4k_{\mathrm{F}}}{\pi}=1.27\ \text{Bohr}^{-2}, \qquad \kappa_{\mathrm{TF}}=1.13\ \text{Bohr}^{-1}$$ である．セルの1辺 $L=20$ Bohr なら $q_{\min}=2\pi/20=0.314$ Bohr$^{-1}$， $$\varepsilon(q_{\min})=1+\frac{1.27}{0.0987}=13.9\quad\Longrightarrow\quad\alpha\lt0.14$$ となる．$L=60$ Bohr なら $q_{\min}=0.105$，$\varepsilon(q_{\min})=1+115=116$，$\alpha\lt0.017$ である．しかも $\alpha$ を小さくすると，今度は短波長成分($\varepsilon\approx1$)の増幅率が $1-\alpha\approx0.983$ となり，こちらが律速になる．誤差を $10^{-4}$ にするのに要する反復数は $$n\simeq\frac{\ln 10^{-4}}{\ln 0.983}\approx 540\ \text{回}$$ である．1反復に対角化が入ることを思えば，これは実用にならない．長波長成分と短波長成分で必要な $\alpha$ が桁違いに違う——この一点が問題の本質であり，次項の解決策を導く．

14.5.3 Kerker混合 — 長波長成分の前処理

上の例が教えるのは，$\alpha$ を波数に依らない定数にしていることが誤りだ，ということである．$q$ ごとに異なる混合パラメータ $\alpha(q)$ を使えばよい．増幅率 \eqref{eq:14-amp-eps} を $q$ に依らない一定値にしたいのだから，要求は

\begin{equation} \alpha(q)\,\varepsilon(q)=\alpha=\text{const} \qquad\Longleftrightarrow\qquad \alpha(q)=\frac{\alpha}{\varepsilon(q)} \label{eq:14-precond-ideal} \end{equation}

である．すなわち理想的な混合とは，ヤコビ演算子の逆(=誘電関数の逆)を前処理として掛けることである．実際の $\varepsilon(q)$ は分からないが，金属では \eqref{eq:14-tf-eps} が良いモデルを与える．これを代入すると

\begin{align} \alpha(q)&=\frac{\alpha}{1+\kappa_{\mathrm{TF}}^{2}/q^{2}} \label{eq:14-kerker-derive1}\\ &=\alpha\,\frac{q^{2}}{q^{2}+\kappa_{\mathrm{TF}}^{2}} \label{eq:14-kerker-derive2} \end{align}

となる(2行目は分母分子に $q^{2}$ を掛けただけ)．これがKerker混合であり，実装では $\kappa_{\mathrm{TF}}$ を調整可能パラメータ $q_0$ に置き換えて

\begin{equation} \widetilde{\delta n}_{\mathrm{in}}^{(n+1)}(\bm{q}) =\big[1-\alpha\,w(q)\big]\,\widetilde{\delta n}_{\mathrm{in}}^{(n)}(\bm{q}) +\alpha\,w(q)\,\widetilde{\delta n}_{\mathrm{out}}^{(n)}(\bm{q}), \qquad w(q)=\frac{q^{2}}{q^{2}+q_0^{2}} \label{eq:14-kerker-mix} \end{equation}

と書く．$w(q)$ をKerker因子と呼ぶ．実際に増幅率が一定になることを確かめよう:

\begin{align} 1-\alpha\,w(q)\,\varepsilon(q) &=1-\alpha\cdot\frac{q^{2}}{q^{2}+q_0^{2}}\cdot\frac{q^{2}+\kappa_{\mathrm{TF}}^{2}}{q^{2}} \label{eq:14-kerker-check1}\\ &=1-\alpha\,\frac{q^{2}+\kappa_{\mathrm{TF}}^{2}}{q^{2}+q_0^{2}} \ \xrightarrow{\ q_0=\kappa_{\mathrm{TF}}\ }\ 1-\alpha \label{eq:14-kerker-check2} \end{align}

1行目では \eqref{eq:14-tf-eps} を通分した形 $\varepsilon(q)=(q^{2}+\kappa_{\mathrm{TF}}^{2})/q^{2}$ を使い，2行目で $q^{2}$ を約分した．$q_0$ を $\kappa_{\mathrm{TF}}$ に合わせれば，すべての波数で増幅率が $1-\alpha$ になり，$L$ 依存性が完全に消える．実際には $\kappa_{\mathrm{TF}}$ は既知でないので，$q_0=\gamma\,q_{\min}$($\gamma\sim1$〜$3$)や $q_0\sim1$ Bohr$^{-1}$ という経験的な選び方をする．$q_0$ を大きくしすぎると長波長成分がほとんど更新されなくなり，逆に収束が遅くなるので，これは調整すべきパラメータである．

注意:$w(0)=0$ が意味すること

$q\to0$ で $w(q)\to0$ である．すなわち $\bm{q}=0$ 成分(=全電子数)は決して混合されない．これは正しい振る舞いである．$\bm{q}=0$ 成分は電子数保存によって常に固定されており，入力と出力で必ず一致している(誤差はゼロ)．混合する必要も，してよい理由もない．14.4節の差電荷 $\delta n$ を混合の対象にとる実装が多いのは，$\delta n$ が最初から $\widetilde{\delta n}(0)=0$ を満たしていて，この点が自動的に保証されるためでもある．

14.5.4 Pulay(DIIS)混合 — 過去の履歴から最適な入力を作る

Kerker混合はヤコビ演算子のモデルを使った前処理だった．もう一つの方向は，$J$ をモデル化するのではなく，これまでの反復で実際に観測された入出力の関係から $J$ の情報を推定する，という準Newton的な発想である．その代表が Pulay の DIIS(direct inversion in the iterative subspace)法であり，電子構造計算では RMM-DIIS とも呼ばれる．

まず残差を

\begin{equation} R_m\equiv n_{\mathrm{out}}^{(m)}-n_{\mathrm{in}}^{(m)}=F[n_{\mathrm{in}}^{(m)}]-n_{\mathrm{in}}^{(m)} \label{eq:14-resid-def} \end{equation}

と定義する．$R_m=0$ が自己無撞着の条件そのものである．過去 $p$ 回の入力密度の線形結合

\begin{equation} \bar n_{\mathrm{in}}=\sum_{m=n-p+1}^{n}c_m\,n_{\mathrm{in}}^{(m)}, \qquad \sum_{m}c_m=1 \label{eq:14-affine} \end{equation}

を考える．拘束 $\sum_m c_m=1$ を課す理由は次の補題で分かる．

導出:アフィン結合の残差は残差のアフィン結合になる

線形化 \eqref{eq:14-linearize2} のもとで，任意の密度 $n$ の残差は

\begin{equation} R[n]=F[n]-n=n^{*}+J(n-n^{*})-n=(\hat1-J)(n^{*}-n) \label{eq:14-resid-lin} \end{equation}

と書ける(2番目の等号で $n=n^{*}+(n-n^{*})$ と分けて整理した)．これを \eqref{eq:14-affine} の $\bar n_{\mathrm{in}}$ に適用すると

$$ R[\bar n_{\mathrm{in}}]=(\hat1-J)\Big(n^{*}-\sum_m c_m n_{\mathrm{in}}^{(m)}\Big) =(\hat1-J)\Big(\sum_m c_m n^{*}-\sum_m c_m n_{\mathrm{in}}^{(m)}\Big) \label{eq:14-diis-lemma1} $$

となる．ここで2番目の等号で拘束 $\sum_m c_m=1$ を使った($n^{*}=\sum_mc_mn^{*}$ と書き換えた)．まとめれば

$$ R[\bar n_{\mathrm{in}}]=\sum_m c_m\,(\hat1-J)\big(n^{*}-n_{\mathrm{in}}^{(m)}\big)=\sum_m c_m R_m \label{eq:14-diis-lemma2} $$

である(最後は再び \eqref{eq:14-resid-lin} を各 $m$ に適用した)．すなわち過去の入力をアフィン結合すると，その残差は過去の残差の同じ係数によるアフィン結合になる．$R_m$ は既に計算済みの既知量だから，これは「まだ計算していない密度の残差を，係数だけで予測できる」ことを意味する．拘束 $\sum c_m=1$ は，この予測が成り立つための必要条件である($\sum c_m\neq1$ だと $n^{*}$ の項が残ってしまう)．

そこで，予測残差のノルムを最小にする係数を選ぶ．最小化すべきは

\begin{equation} \Big\|\sum_m c_m R_m\Big\|^{2}=\sum_{m m'}c_m\,B_{mm'}\,c_{m'}, \qquad B_{mm'}\equiv\braket{R_m|R_{m'}} \label{eq:14-diis-obj} \end{equation}

であり，拘束は $\sum_m c_m=1$ である(以下 $c_m$ と $R_m$ は実数・実関数とする．$B$ は実対称でグラム行列なので半正定値)．拘束付き最小化なのでLagrangeの未定乗数 $\mu$ を導入し

\begin{equation} \mathcal{L}(\{c_m\},\mu)=\sum_{mm'}c_mB_{mm'}c_{m'}-2\mu\Big(\sum_m c_m-1\Big) \label{eq:14-diis-lagrangian} \end{equation}

を停留化する(未定乗数を $2\mu$ と書いたのは，後で $2$ が消えて式がきれいになるからで，本質ではない)．$c_k$ で微分すると，第1項からは $B$ の対称性により

\begin{equation} \frac{\partial}{\partial c_k}\sum_{mm'}c_mB_{mm'}c_{m'} =\sum_{m'}B_{km'}c_{m'}+\sum_{m}c_mB_{mk} =2\sum_{m'}B_{km'}c_{m'} \label{eq:14-diis-deriv} \end{equation}

が出る(第1の和は $m=k$ の項から，第2の和は $m'=k$ の項から来る．$B_{mk}=B_{km}$ を使って足し合わせた)．第2項の微分は $-2\mu$ だから，停留条件 $\partial\mathcal{L}/\partial c_k=0$ は

\begin{equation} \sum_{m'}B_{km'}c_{m'}=\mu \qquad(k=n-p+1,\dots,n) \label{eq:14-diis-stationary} \end{equation}

となる．これと拘束条件を合わせると，$(p+1)$ 元の連立一次方程式

\begin{equation} \begin{pmatrix} B_{11}&\cdots&B_{1p}&-1\\ \vdots&\ddots&\vdots&\vdots\\ B_{p1}&\cdots&B_{pp}&-1\\ 1&\cdots&1&0 \end{pmatrix} \begin{pmatrix}c_1\\ \vdots\\ c_p\\ \mu\end{pmatrix} = \begin{pmatrix}0\\ \vdots\\ 0\\ 1\end{pmatrix} \label{eq:14-diis-eq} \end{equation}

が得られる(添字を $1,\dots,p$ に振り直した)．この縁付き連立方程式を解けば係数 $\{c_m\}$ が決まる．$p$ はせいぜい $5$〜$20$ なので，この求解のコストは対角化に比べて無視できる．

得られた係数から次の入力密度を作る．単に $\bar n_{\mathrm{in}}=\sum_mc_mn_{\mathrm{in}}^{(m)}$ とするだけでは，新しい密度は過去の入力が張るアフィン空間から出られない．そこで最適残差の一部を足して新しい方向に踏み出す:

\begin{equation} n_{\mathrm{in}}^{(n+1)}=\sum_m c_m\,n_{\mathrm{in}}^{(m)}+\beta\sum_m c_m R_m \label{eq:14-pulay-update} \end{equation}

$\beta$ は単純混合の $\alpha$ と同じ役割の小さな正数である．実用的な実装では，この $\sum_mc_mR_m$ に Kerker因子を掛けて長波長成分を抑え，さらに内積を

\begin{equation} \braket{R_m|R_{m'}}=\sum_{\bm{q}}\frac{R_m^{*}(\bm{q})\,R_{m'}(\bm{q})}{w(q)} \label{eq:14-kerker-metric} \end{equation}

というKerker計量で測る．$1/w(q)$ は小さい $q$ で大きいので，この計量は長波長の残差を重く評価する．危険な成分を優先的に潰す，という意図である．

物理的意味:DIISは「ヤコビアンを測りながら」進む準Newton法である

不動点問題 $R[n]=0$ をNewton法で解くなら，更新は $n^{(n+1)}=n^{(n)}-(\partial R/\partial n)^{-1}R^{(n)}$，すなわち \eqref{eq:14-resid-lin} より $n^{(n+1)}=n^{(n)}+(\hat1-J)^{-1}R^{(n)}$ である．$J$ を陽に計算するのは高価だが，DIISは過去 $p$ 回の $(n_{\mathrm{in}}^{(m)},R_m)$ の組を「$J$ の作用のサンプル」とみなし，その張る部分空間の中で $(\hat1-J)^{-1}$ を近似している．単純混合は $(\hat1-J)^{-1}\approx\alpha\hat1$ という最も粗い近似，Kerker混合は $(\hat1-J)^{-1}\approx\alpha\,\varepsilon^{-1}$ というモデル近似，DIISは履歴からの経験的近似——3者はすべて同じ準Newton法の枠組みの中にある．だからこそ，Kerker前処理とDIISは競合せず，組み合わせるのが最も強力なのである．同じ考え方が第16章の構造最適化(座標に対するDIIS/BFGS)にも現れる．

図14.5 (a) charge sloshing の模式図．長い胞では長波長の電荷の偏りが反復ごとに符号を変えて増幅する．(b) 残差ノルムの収束の模式図(縦軸は対数)．単純混合は長波長成分に律速されて停滞し，Kerker前処理でその律速が外れ，さらにPulay(DIIS)で履歴を使うと収束が加速する．

14.5.5 実践的なコツ

(i) 金属ではsmearingが必須である．絶対零度の占有数はステップ関数 $f_i=\theta(\mu-\varepsilon_i)$ であり，フェルミ準位をまたぐ状態が1つでもあると，$v_{\mathrm{eff}}$ の無限小の変化で占有数が $0$ と $1$ の間を飛ぶ．すると $F$ は不連続になり，\eqref{eq:14-linearize2} の線形化そのものが成立しない．どんな混合法を使っても収束しないのは当然である．処方は占有数を滑らかにすること，すなわち

\begin{equation} f_i=\frac{1}{1+\exp\!\big[(\varepsilon_i-\mu)/\sigma\big]} \label{eq:14-fermi-smear} \end{equation}

という有限幅 $\sigma$ のFermi分布(または Methfessel–Paxton 型の展開)を使うことである．これにより $\partial f/\partial\varepsilon$ が有限($\varepsilon=\mu$ で $-1/(4\sigma)$)になり，$F$ が微分可能になる．副次的な効果として，Brillouin域積分の $\kk$ 点収束も大幅に改善される(第13章の状態密度の議論と同じ理由:フェルミ面の鋭さが緩和される)．ただし $\sigma\neq0$ では変分量が全エネルギーではなく自由エネルギー $\Omega_{\mathrm{F}}=E-\sigma S_{\mathrm{el}}$ になるので，力やエネルギーを議論するときは $\sigma\to0$ 外挿(Fermi分布なら $E_0\simeq(E+\Omega_{\mathrm{F}})/2$ の関係が知られている [Gillan 1989])が必要である．記号の注意:ここだけ $S_{\mathrm{el}}=-k_{\mathrm{B}}\sum_i[f_i\ln f_i+(1-f_i)\ln(1-f_i)]$ は占有数のエントロピーであり，14.2節の重なり行列 $S$ とは別物である．同様に $\Omega_{\mathrm{F}}$ は自由エネルギーであって，14.4節以降で単位胞の体積を表す $\Omega$ とは無関係である(混同を避けるため本書では添字を付けて区別する)．

(ii) 混合パラメータの選び方．\eqref{eq:14-alpha-eps} が上限を与える．絶縁体・分子なら $\alpha=0.2$〜$0.5$，金属バルクなら $0.1$ 前後，長いスラブや磁性金属なら $0.02$〜$0.05$ から始めるのが安全である．Kerker前処理を入れれば $\alpha$ をセルサイズに依らず選べるようになるので，まずKerkerを有効にしてから $\alpha$ を上げていくのが実務的な順序である．DIISの履歴長 $p$ は $5$〜$20$ 程度．$B$ 行列が悪条件になったら(過去の残差がほぼ一次従属になったら)履歴を捨ててやり直す．

(iii) 初期密度．14.4節の例で述べたように，原子密度の重ね合わせ $n^{(a)}$ が標準である．加えて，磁性系では初期スピン密度に有限の分極を与えて対称性を破っておかないと，無偏極解($\zeta=0$)は\eqref{eq:14-scf-map} の不動点でもあるため，そこに留まってしまう．反強磁性秩序を狙うなら，サイトごとに反対向きの初期モーメントを与える必要がある．SCFは最も近い不動点を見つけるだけで，最低エネルギーの解を保証しない——これは初学者が最も陥りやすい落とし穴である．

表14.5 混合法の位置づけ
方法	$(\hat 1-J)^{-1}$ の近似	必要な情報	長所と限界
単純混合	$\alpha\hat 1$	なし	実装が容易．$\alpha\lt2/\varepsilon(q_{\min})$ に縛られる
Kerker混合	$\alpha\,q^{2}/(q^{2}+q_0^{2})$	モデル誘電関数	セルサイズ依存性を除去．$q_0$ の調整が要る
Pulay(DIIS)	履歴が張る部分空間内で最適	過去 $p$ 回の $(n_{\mathrm{in}},R)$	収束が速い．履歴が一次従属になると不安定
Kerker + DIIS	前処理付き履歴近似	両方	実用コードの標準的な組み合わせ

14.6 原子のDFT計算 — 動径方程式を解く

DFTの数値解法を本当に理解する最短経路は，原子のKS方程式を自分で解くコードを書いてみることである．球対称という強い仮定のおかげで問題は1次元常微分方程式に落ち，それでいて「固有値問題を解く」「Poisson方程式を解く」「交換相関ポテンシャルを作る」「SCFを回す」という要素がすべて登場する．しかも原子の計算は，14.3節の擬原子軌道の生成にも，第15章の擬ポテンシャル構成にも，そのまま必要になる．この節では，その設計図を最後まで書き下す．

14.6.1 球対称ポテンシャルでの変数分離

孤立原子(スピン非分極，球対称近似)のKS方程式は

\begin{equation} \left[-\frac{1}{2}\nabla^{2}+v_{\mathrm{eff}}(r)\right]\varphi_{nlm}(\rr)=\varepsilon_{nl}\,\varphi_{nlm}(\rr), \qquad v_{\mathrm{eff}}(r)=-\frac{Z}{r}+v_{\mathrm{H}}(r)+v_{xc}\big(n(r)\big) \label{eq:14-atomic-ks} \end{equation}

である．ポテンシャルが $r$ だけの関数であることが，以下すべての出発点である．

数学ノート:球座標のラプラシアンと球面調和関数

球座標 $(r,\theta,\phi)$ でのラプラシアンは

\begin{equation} \nabla^{2} =\frac{1}{r^{2}}\frac{\partial}{\partial r}\left(r^{2}\frac{\partial}{\partial r}\right) +\frac{1}{r^{2}\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial}{\partial\theta}\right) +\frac{1}{r^{2}\sin^{2}\theta}\frac{\partial^{2}}{\partial\phi^{2}} \label{eq:14-lap-sph} \end{equation}

である．第2項と第3項は角度変数だけを含み，$r^{-2}$ という共通因子を持つ．そこで角運動量演算子の2乗を(原子単位 $\hbar=1$ で)

$$ \hat L^{2}\equiv-\left[\frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial}{\partial\theta}\right) +\frac{1}{\sin^{2}\theta}\frac{\partial^{2}}{\partial\phi^{2}}\right] \label{eq:14-lsq-def} $$

と定義すると，\eqref{eq:14-lap-sph} は簡潔に

\begin{equation} \nabla^{2}=\frac{1}{r^{2}}\frac{\partial}{\partial r}\left(r^{2}\frac{\partial}{\partial r}\right)-\frac{\hat L^{2}}{r^{2}} \label{eq:14-lap-compact} \end{equation}

と書ける．$\hat L^{2}$ の固有関数が球面調和関数 $Y_{lm}(\theta,\phi)$ であり

\begin{equation} \hat L^{2}Y_{lm}=l(l+1)\,Y_{lm}, \qquad l=0,1,2,\dots,\quad m=-l,\dots,l \label{eq:14-ylm-eig} \end{equation}

が成り立つ．$Y_{lm}$ は球面上で完全正規直交系をなす:

\begin{equation} \int \dd\Omega\;Y_{lm}^{*}(\hat\rr)\,Y_{l'm'}(\hat\rr)=\delta_{ll'}\delta_{mm'} \label{eq:14-ylm-orth} \end{equation}

本節で後に使うもう一つの性質は，同じ $l$ の $\abs{Y_{lm}}^{2}$ を $m$ について足すと角度依存性が消えることである(Unsöldの定理):

\begin{equation} \sum_{m=-l}^{l}\abs{Y_{lm}(\hat\rr)}^{2}=\frac{2l+1}{4\pi} \label{eq:14-unsold} \end{equation}

これは，閉殻(ある $l$ の $2l+1$ 個の軌道がすべて等しく占有された状態)の電子密度が厳密に球対称になることを意味する．

波動関数を動径部分と角度部分に分離して

\begin{equation} \varphi_{nlm}(\rr)=R_{nl}(r)\,Y_{lm}(\theta,\phi) \label{eq:14-sep-ansatz} \end{equation}

と置き，\eqref{eq:14-lap-compact} と \eqref{eq:14-ylm-eig} を使って \eqref{eq:14-atomic-ks} に代入する:

\begin{align} -\frac{1}{2}\left[\frac{Y_{lm}}{r^{2}}\frac{\dd}{\dd r}\left(r^{2}\frac{\dd R_{nl}}{\dd r}\right) -\frac{R_{nl}}{r^{2}}\,l(l+1)\,Y_{lm}\right]+v_{\mathrm{eff}}R_{nl}Y_{lm} &=\varepsilon_{nl}R_{nl}Y_{lm} \label{eq:14-sep1} \end{align}

$\hat L^{2}$ は $R_{nl}(r)$ に作用しない(角度変数だけの微分だから)ので $Y_{lm}$ にだけ作用し，\eqref{eq:14-ylm-eig} により $l(l+1)$ に置き換わった．全体に共通因子 $Y_{lm}$ があるのでこれで割ると，角度変数が完全に消えて

\begin{equation} -\frac{1}{2r^{2}}\frac{\dd}{\dd r}\left(r^{2}\frac{\dd R_{nl}}{\dd r}\right) +\left[\frac{l(l+1)}{2r^{2}}+v_{\mathrm{eff}}(r)\right]R_{nl}=\varepsilon_{nl}R_{nl} \label{eq:14-radial-R} \end{equation}

という $r$ だけの常微分方程式になる．角度部分は解析的に処理され，$l(l+1)/(2r^{2})$ という遠心力ポテンシャルだけを残していった．

14.6.2 $u=rR$ の置換で標準形へ

\eqref{eq:14-radial-R} には1階微分の項が隠れており($\frac{1}{r^2}(r^2R')'=R''+\frac{2}{r}R'$)，このままでは数値積分に不便である．標準的な処方は

\begin{equation} u_{nl}(r)\equiv r\,R_{nl}(r) \qquad\Longleftrightarrow\qquad R_{nl}=\frac{u_{nl}}{r} \label{eq:14-u-def} \end{equation}

と置くことである．実際に計算してみよう．まず微分を実行する:

\begin{align} \frac{\dd R}{\dd r}&=\frac{\dd}{\dd r}\left(\frac{u}{r}\right)=\frac{u'}{r}-\frac{u}{r^{2}} \label{eq:14-usub1}\\ \frac{\dd^{2}R}{\dd r^{2}}&=\frac{u''}{r}-\frac{u'}{r^{2}}-\left(\frac{u'}{r^{2}}-\frac{2u}{r^{3}}\right) =\frac{u''}{r}-\frac{2u'}{r^{2}}+\frac{2u}{r^{3}} \label{eq:14-usub2} \end{align}

(1行目は商の微分則，2行目は1行目をもう一度 $r$ で微分し，$\dd(u'/r)/\dd r=u''/r-u'/r^{2}$ と $\dd(u/r^{2})/\dd r=u'/r^{2}-2u/r^{3}$ を使った)．次に \eqref{eq:14-radial-R} の運動エネルギー部分を展開形で書き直す:

\begin{align} \frac{1}{r^{2}}\frac{\dd}{\dd r}\left(r^{2}\frac{\dd R}{\dd r}\right) &=\frac{1}{r^{2}}\left(2r\frac{\dd R}{\dd r}+r^{2}\frac{\dd^{2}R}{\dd r^{2}}\right) =\frac{\dd^{2}R}{\dd r^{2}}+\frac{2}{r}\frac{\dd R}{\dd r} \label{eq:14-usub3} \end{align}

(積の微分則を使い，$r^{2}$ で割った)．ここに \eqref{eq:14-usub1}，\eqref{eq:14-usub2} を代入すると

\begin{align} \frac{\dd^{2}R}{\dd r^{2}}+\frac{2}{r}\frac{\dd R}{\dd r} &=\left(\frac{u''}{r}-\frac{2u'}{r^{2}}+\frac{2u}{r^{3}}\right) +\frac{2}{r}\left(\frac{u'}{r}-\frac{u}{r^{2}}\right) \label{eq:14-usub4}\\ &=\frac{u''}{r} \underbrace{-\frac{2u'}{r^{2}}+\frac{2u'}{r^{2}}}_{=0} \underbrace{+\frac{2u}{r^{3}}-\frac{2u}{r^{3}}}_{=0} =\frac{u''}{r} \label{eq:14-usub5} \end{align}

となり，余分な項がすべて相殺する．これが $u=rR$ という置換の御利益である．これを \eqref{eq:14-radial-R} に戻し，両辺に $r$ を掛けると($R=u/r$ も代入して)

\begin{equation} \left[-\frac{1}{2}\frac{\dd^{2}}{\dd r^{2}}+\frac{l(l+1)}{2r^{2}}+v_{\mathrm{eff}}(r)\right]u_{nl}(r) =\varepsilon_{nl}\,u_{nl}(r) \label{eq:14-radial-u} \end{equation}

を得る．これは1次元のSchrödinger方程式と全く同じ形であり，有効ポテンシャルは $W_l(r)=v_{\mathrm{eff}}(r)+l(l+1)/(2r^{2})$，定義域は $r\in[0,\infty)$ である．$l\ge1$ では遠心力障壁が原点で $+\infty$ に発散し，電子を原点から遠ざける．

規格化条件も簡単になる．\eqref{eq:14-ylm-orth} を使って角度積分を実行すると

\begin{equation} \int\dd^3r\,\abs{\varphi_{nlm}}^{2} =\int_0^\infty\!\!\dd r\,r^{2}\abs{R_{nl}}^{2}\int\dd\Omega\,\abs{Y_{lm}}^{2} =\int_0^\infty\!\!\dd r\,\abs{u_{nl}(r)}^{2}=1 \label{eq:14-u-norm} \end{equation}

となる($r^{2}R^{2}=u^{2}$ による)．$u$ は「1次元の規格化された波動関数」そのものである．

14.6.3 境界条件と漸近形

常微分方程式を数値積分するには両端での振る舞いが要る．

原点近傍．$r\to0$ では遠心力項 $l(l+1)/(2r^{2})$ が $v_{\mathrm{eff}}\sim-Z/r$ よりも速く発散するので，$l\ge1$ では最も特異な項どうしが釣り合う．\eqref{eq:14-radial-u} の主要部だけを残すと

\begin{equation} \frac{\dd^{2}u}{\dd r^{2}}\simeq\frac{l(l+1)}{r^{2}}u \label{eq:14-origin-eq} \end{equation}

である．冪級数の先頭を $u\propto r^{s}$ と置くと $\dd^{2}u/\dd r^{2}=s(s-1)r^{s-2}$，右辺は $l(l+1)r^{s-2}$ だから，指数方程式

\begin{equation} s(s-1)=l(l+1) \quad\Longrightarrow\quad s=l+1\ \ \text{または}\ \ s=-l \label{eq:14-indicial} \end{equation}

を得る(2次方程式 $s^{2}-s-l(l+1)=0$ の2解)．$s=-l$ は $R=u/r\propto r^{-l-1}$ を与え，$l\ge1$ では規格化積分 \eqref{eq:14-u-norm} が発散し，$l=0$ でも $R\propto1/r$ となって $\nabla^{2}(1/r)=-4\pi\delta(\rr)$ が余分なデルタ関数源を生むので，いずれも物理解として棄却される．したがって

\begin{equation} u_{nl}(r)\ \xrightarrow[r\to0]{}\ C\,r^{\,l+1}\left(1+a_1r+a_2r^{2}+\cdots\right) \label{eq:14-origin-power} \end{equation}

である．係数 $a_1,a_2,\dots$ は \eqref{eq:14-radial-u} に級数を代入して同じ冪の係数を比較すれば順に決まる(Frobenius法)．$l=0$ の場合の $a_1=-Z$ は，よく知られた核カスプ条件に対応する．

無限遠．$r\to\infty$ では $v_{\mathrm{eff}}\to0$(中性原子)，遠心力項も $r^{-2}$ で消えるので

\begin{equation} \frac{\dd^{2}u}{\dd r^{2}}\simeq-2\varepsilon_{nl}\,u=2\abs{\varepsilon_{nl}}u \qquad(\varepsilon_{nl}\lt0) \label{eq:14-asym-eq} \end{equation}

となる．定数係数の2階線形方程式なので一般解は $u=Ae^{-\kappa r}+Be^{+\kappa r}$，$\kappa=\sqrt{2\abs{\varepsilon_{nl}}}$ である．境界条件 $u(\infty)=0$ が $B=0$ を要求し

\begin{equation} u_{nl}(r)\ \xrightarrow[r\to\infty]{}\ A\,e^{-\sqrt{2\abs{\varepsilon_{nl}}}\,r} \label{eq:14-asym} \end{equation}

を得る．これが 14.3.2 節で先取りして使った \eqref{eq:14-free-tail} である．

14.6.4 対数メッシュ

数値積分の格子をどう取るか．ここには実質的な困難がある．重い原子では $1s$ 軌道の半径が $\sim a_0/Z$($Z=30$ なら $0.033$ Bohr)であるのに対し，価電子の裾は $30$〜$50$ Bohr まで伸びる．両方を等間隔格子で解像しようとすると，内側に合わせて $\Delta r\sim10^{-3}$ Bohr，外側まで $50$ Bohr なので $5\times10^{4}$ 点が必要になる．しかも大半の点は，波動関数がほとんど変化しない外側領域に費やされる．

解決策は対数メッシュである．新しい変数

\begin{equation} x=\ln r \qquad\Longleftrightarrow\qquad r=e^{x}, \qquad \frac{\dd r}{\dd x}=r \label{eq:14-logmesh} \end{equation}

を導入し，$x$ について等間隔の格子 $x_i=x_{\min}+i\,\Delta x$ を取る．このとき $r$ 空間での間隔は $\Delta r=r\,\Delta x$ であり，相対分解能 $\Delta r/r=\Delta x$ が全域で一定になる．原子の波動関数は「$r$ が2倍になるあいだに1回程度変化する」という自己相似的な構造を持つので，これがまさに欲しい格子である．$x\in[-11.5,\,3.9]$($r\in[10^{-5},50]$ Bohr)を $N=3000$ 点で刻めば $\Delta x=0.0051$，すなわちどこでも $0.5$ %の相対分解能が得られる．等間隔格子の $5\times10^{4}$ 点に対して一桁以上の節約である．

方程式を $x$ で書き換えよう．合成関数の微分則から

\begin{align} \frac{\dd u}{\dd r}&=\frac{\dd x}{\dd r}\frac{\dd u}{\dd x}=\frac{1}{r}\frac{\dd u}{\dd x} \label{eq:14-dudr}\\ \frac{\dd^{2}u}{\dd r^{2}} &=\frac{1}{r}\frac{\dd}{\dd x}\left(\frac{1}{r}\frac{\dd u}{\dd x}\right) =\frac{1}{r}\left[-\frac{1}{r}\frac{\dd u}{\dd x}+\frac{1}{r}\frac{\dd^{2}u}{\dd x^{2}}\right] =\frac{1}{r^{2}}\left[\frac{\dd^{2}u}{\dd x^{2}}-\frac{\dd u}{\dd x}\right] \label{eq:14-d2udr2} \end{align}

である．\eqref{eq:14-dudr} は $\dd x/\dd r=1/r$($x=\ln r$ の微分)による．\eqref{eq:14-d2udr2} では，まず $\dd/\dd r=(1/r)\dd/\dd x$ を \eqref{eq:14-dudr} の右辺に作用させ，積 $(1/r)(\dd u/\dd x)$ の微分に積の微分則を使った．ここで $\dd(1/r)/\dd x=-(1/r^{2})(\dd r/\dd x)=-1/r$ を用いている(\eqref{eq:14-logmesh} の $\dd r/\dd x=r$)．

これを \eqref{eq:14-radial-u} に代入し，両辺に $-2r^{2}$ を掛けて整理すると

\begin{equation} \frac{\dd^{2}u}{\dd x^{2}}-\frac{\dd u}{\dd x} =\Big[l(l+1)+2r^{2}\big(v_{\mathrm{eff}}(r)-\varepsilon\big)\Big]u \label{eq:14-radial-x} \end{equation}

となる．都合の悪いことに，1階微分の項 $-\dd u/\dd x$ が現れた．これがあると，2階微分専用の高精度公式(Numerov法など)がそのままでは使えない．対処法は2通りある．

方法A:1階微分項を消す置換．$u=e^{x/2}y=\sqrt{r}\,y$ と置く．微分を実行すると

\begin{align} \frac{\dd u}{\dd x}&=e^{x/2}\left(\frac{1}{2}y+\frac{\dd y}{\dd x}\right) \label{eq:14-yA1}\\ \frac{\dd^{2}u}{\dd x^{2}}&=e^{x/2}\left(\frac{1}{4}y+\frac{\dd y}{\dd x}+\frac{\dd^{2}y}{\dd x^{2}}\right) \label{eq:14-yA2} \end{align}

である(積の微分則を2回．$\dd e^{x/2}/\dd x=\frac{1}{2}e^{x/2}$ を使った)．これらを \eqref{eq:14-radial-x} に代入し，共通因子 $e^{x/2}$ で割ると

\begin{align} \frac{1}{4}y+\frac{\dd y}{\dd x}+\frac{\dd^{2}y}{\dd x^{2}}-\frac{1}{2}y-\frac{\dd y}{\dd x} &=\Big[l(l+1)+2r^{2}(v_{\mathrm{eff}}-\varepsilon)\Big]y \label{eq:14-yA3}\\ \frac{\dd^{2}y}{\dd x^{2}} &=\Big[l(l+1)+\frac{1}{4}+2r^{2}(v_{\mathrm{eff}}-\varepsilon)\Big]y \label{eq:14-yA4} \end{align}

を得る．1行目で $\dd y/\dd x$ の項が相殺し，2行目では $\frac{1}{4}y-\frac{1}{2}y=-\frac{1}{4}y$ を右辺に移した．最後に $l(l+1)+\frac{1}{4}=\left(l+\frac{1}{2}\right)^{2}$ と平方完成すれば

\begin{equation} \frac{\dd^{2}y}{\dd x^{2}} =\left[\left(l+\tfrac{1}{2}\right)^{2}+2r^{2}\big(v_{\mathrm{eff}}(r)-\varepsilon\big)\right]y, \qquad u=\sqrt{r}\,y \label{eq:14-numerov-form} \end{equation}

という美しい形になる．1階微分項がないので，Numerov法(局所誤差 $O(\Delta x^{6})$)がそのまま適用できる．

方法B:1階の連立系にする．もう一つの実装は，$u=r^{l+1}L(r)$ と置いて $L$ について解く方法である．同じ手順(\eqref{eq:14-radial-x} に $u=e^{(l+1)x}L$ を代入し，$e^{(l+1)x}$ で割る)を実行すると，$L$ の係数は $(l+1)^{2}-(l+1)-l(l+1)=(l+1)[(l+1)-1-l]=0$ で消え，$\dd L/\dd x$ の係数は $2(l+1)-1=2l+1$ となって

\begin{equation} \frac{\dd L}{\dd x}=M, \qquad \frac{\dd M}{\dd x}=-(2l+1)M+2r^{2}\big(v_{\mathrm{eff}}(r)-\varepsilon\big)L \label{eq:14-first-order-sys} \end{equation}

という1階連立系が得られる．この形の利点は，\eqref{eq:14-origin-power} により $L(r)\to C$(定数)，$M\to0$ が原点での自然な初期値になり，$r^{l+1}$ という強い冪をコードが扱わなくてよいことである．予測子–修正子法(Adams型)を使うのが標準である．以下では方法A・Bのどちらでもよいので，単に「動径方程式を数値積分する」と呼ぶ．

図14.6 (a) 対数メッシュ．$x=\ln r$ について等間隔にとった格子点を $r$ 軸に写すと原点付近に密に集まり，相対分解能が全域で一定になる．(b) 有効動径ポテンシャル $W_l(r)=v_{\mathrm{eff}}(r)+l(l+1)/(2r^{2})$．原点側は遠心力障壁，外側は $0$ に漸近する．固有値 $\varepsilon_{nl}$ と交わる点(古典的転回点)がマッチング点 $r_{\mathrm{mp}}$ の目安になる．

14.6.5 シューティング法 — 両側から積分してつなぐ

\eqref{eq:14-radial-u} は $\varepsilon$ を含む固有値問題である．$\varepsilon$ を勝手に決めて数値積分すると，一般には境界条件のどちらか一方しか満たせない．この事実を逆手に取るのがシューティング法である．

手順は次の通り．(i) 原点側から $r_{\min}$ で \eqref{eq:14-origin-power} の初期値を与え，外向きに $r_{\mathrm{mp}}$ まで積分する．これを $u^{\mathrm{out}}$ と呼ぶ．(ii) 外側から $r_{\max}$ で \eqref{eq:14-asym} の初期値を与え，内向きに $r_{\mathrm{mp}}$ まで積分する．これを $u^{\mathrm{in}}$ と呼ぶ．(iii) 両者を $r_{\mathrm{mp}}$ で接続する．

接続条件を正しく書くには一つ注意が要る．\eqref{eq:14-radial-u} は斉次線形方程式なので，解を定数倍しても解である．したがって $u^{\mathrm{out}}$ と $u^{\mathrm{in}}$ の絶対的な大きさには意味がない．そこで一方をスケールして $u^{\mathrm{out}}(r_{\mathrm{mp}})=u^{\mathrm{in}}(r_{\mathrm{mp}})$ を常に満たすようにしておき，残る条件として対数微分の連続性

\begin{equation} D(\varepsilon)\equiv \left.\frac{1}{u^{\mathrm{out}}}\frac{\dd u^{\mathrm{out}}}{\dd r}\right|_{r_{\mathrm{mp}}} -\left.\frac{1}{u^{\mathrm{in}}}\frac{\dd u^{\mathrm{in}}}{\dd r}\right|_{r_{\mathrm{mp}}} =0 \label{eq:14-logderiv} \end{equation}

を課す．対数微分 $u'/u$ を使う理由は明快である:$u\to cu$ と定数倍しても $u'/u$ は変わらない(分子分母に同じ $c$ が掛かる)．すなわち対数微分はスケール不変な，解の「形」だけを表す量である．$D(\varepsilon)=0$ を満たす $\varepsilon$ だけが，両側の境界条件を同時に満たす真の固有値である．

どの固有値を捕まえたのかを判定するにはノード数を数える．Sturmの振動定理により，動径方程式の固有関数は固有値の低い順に節が $0,1,2,\dots$ と増える．主量子数 $n$ と角運動量 $l$ の状態については

\begin{equation} (\text{$u_{nl}$ の $0\lt r\lt\infty$ での節の数})=n-l-1 \label{eq:14-nodes} \end{equation}

である($1s$:0個，$2s$:1個，$2p$:0個，$3d$:0個)．積分中に符号反転の回数を数えれば，試した $\varepsilon$ が目標の状態より高いか低いかが分かるので，二分法の区間を確実に括ることができる．

導出:固有値の1次補正(Newton法のステップ)

二分法は確実だが遅い．マッチング点での微分の食い違いから，$\varepsilon$ の補正量を直接見積もることができる．\eqref{eq:14-radial-u} を

$$ \frac{\dd^{2}u}{\dd r^{2}}=2\big[W_l(r)-\varepsilon\big]u, \qquad W_l(r)\equiv v_{\mathrm{eff}}(r)+\frac{l(l+1)}{2r^{2}} \label{eq:14-schr-1d} $$

と書く．試行値 $\varepsilon$ での解を $u$(両側から積分し，$r_{\mathrm{mp}}$ で値だけ一致させたもの)，真の固有値 $\hat\varepsilon=\varepsilon+\Delta\varepsilon$ での厳密解を $\hat u$ とする．それぞれ

$$ u''=2(W_l-\varepsilon)u, \qquad \hat u''=2(W_l-\varepsilon-\Delta\varepsilon)\hat u \label{eq:14-two-eqs} $$

を満たす．第1式に $\hat u$ を，第2式に $u$ を掛けて差をとると，$W_l$ を含む項が相殺して

\begin{equation} \hat u\,u''-u\,\hat u'' =-2\varepsilon u\hat u+2(\varepsilon+\Delta\varepsilon)u\hat u =2\,\Delta\varepsilon\;u\,\hat u \label{eq:14-wronskian1} \end{equation}

が残る．左辺は完全微分である:

$$ \frac{\dd}{\dd r}\big(\hat u\,u'-u\,\hat u'\big) =\hat u'u'+\hat u u''-u'\hat u'-u\hat u'' =\hat u u''-u\hat u'' \label{eq:14-wronskian2} $$

($\hat u'u'$ と $u'\hat u'$ が相殺する)．したがって \eqref{eq:14-wronskian1} を区間ごとに積分できる．まず $[0,r_{\mathrm{mp}}]$ で外向き解について:

$$ \Big[\hat u\,u'-u\,\hat u'\Big]_{0}^{r_{\mathrm{mp}}} =\hat u(r_{\mathrm{mp}})\,u^{\mathrm{out}\prime}(r_{\mathrm{mp}})-u(r_{\mathrm{mp}})\,\hat u'(r_{\mathrm{mp}}) =2\,\Delta\varepsilon\int_{0}^{r_{\mathrm{mp}}}u\hat u\,\dd r \label{eq:14-wronskian3} $$

($r=0$ では $u,\hat u$ ともにゼロなので下端の寄与は消える)．次に $[r_{\mathrm{mp}},\infty)$ で内向き解について:

$$ \Big[\hat u\,u'-u\,\hat u'\Big]_{r_{\mathrm{mp}}}^{\infty} =-\hat u(r_{\mathrm{mp}})\,u^{\mathrm{in}\prime}(r_{\mathrm{mp}})+u(r_{\mathrm{mp}})\,\hat u'(r_{\mathrm{mp}}) =2\,\Delta\varepsilon\int_{r_{\mathrm{mp}}}^{\infty}u\hat u\,\dd r \label{eq:14-wronskian4} $$

($r\to\infty$ で両者は指数的にゼロなので上端の寄与は消える)．2式を足すと $\hat u'$ を含む項が相殺し，$u_m\equiv u(r_{\mathrm{mp}})$，$\hat u(r_{\mathrm{mp}})\simeq u_m$(1次の精度では区別しない)として

$$ u_m\Big[u^{\mathrm{out}\prime}(r_{\mathrm{mp}})-u^{\mathrm{in}\prime}(r_{\mathrm{mp}})\Big] =2\,\Delta\varepsilon\int_{0}^{\infty}u^{2}\,\dd r \label{eq:14-wronskian5} $$

すなわち

\begin{equation} \Delta\varepsilon =\frac{u_m\left[u^{\mathrm{out}\prime}(r_{\mathrm{mp}})-u^{\mathrm{in}\prime}(r_{\mathrm{mp}})\right]} {2\displaystyle\int_{0}^{\infty}u^{2}\,\dd r} \label{eq:14-eps-correction} \end{equation}

を得る．マッチング点での微分の食い違いを，規格化積分で割るだけで，固有値の補正量が求まる．これはNewton法の1ステップに相当し，二分法で十分に括ったあとに数回適用すれば $10^{-12}$ Ha 程度まで一気に収束する．

図14.7 シューティング法．(a) 原点側から外向きに，無限遠側から内向きに積分し，$r_{\mathrm{mp}}$ で値を合わせる．試行固有値が正しくないと微分が食い違う．(b) 判定量 $D(\varepsilon)$ は固有値を横切るたびに符号を変える．ノード数で状態を識別し，二分法で括ってから \eqref{eq:14-eps-correction} で仕上げる．

14.6.6 動径Hartreeポテンシャル

固有関数が得られたら密度を作る．占有数 $f_{nl}$(その $nl$ 殻に入っている電子数，スピン込み)を使うと，Unsöldの定理 \eqref{eq:14-unsold} により

\begin{align} n(r)&=\sum_{nl}\frac{f_{nl}}{2l+1}\sum_{m=-l}^{l}\abs{R_{nl}(r)Y_{lm}}^{2} =\sum_{nl}\frac{f_{nl}}{2l+1}\,\abs{R_{nl}}^{2}\cdot\frac{2l+1}{4\pi} \label{eq:14-dens-radial1}\\ &=\sum_{nl}f_{nl}\,\frac{u_{nl}(r)^{2}}{4\pi r^{2}} \label{eq:14-dens-radial} \end{align}

となる(1行目では $f_{nl}$ を $2l+1$ 個の $m$ に等分し，2行目で $R=u/r$ を代入した)．この密度は厳密に球対称である．検算しておこう:$\int n\,\dd^3r=\int_0^\infty 4\pi r^{2}n(r)\dd r=\sum_{nl}f_{nl}\int_0^\infty u_{nl}^{2}\dd r=\sum_{nl}f_{nl}$ となり，\eqref{eq:14-u-norm} の規格化により確かに全電子数に一致する．

次にHartreeポテンシャルである．球対称なので，数学ノートで導いた \eqref{eq:14-dvdr} が使える:

\begin{equation} \frac{\dd v_{\mathrm{H}}}{\dd r}=-\frac{Q(r)}{r^{2}}, \qquad Q(r)=4\pi\int_0^r n(r')r'^{2}\dd r' \label{eq:14-vh-step1} \end{equation}

境界条件 $v_{\mathrm{H}}(\infty)=0$ のもとで $r$ から $\infty$ まで積分すると

\begin{align} v_{\mathrm{H}}(\infty)-v_{\mathrm{H}}(r)&=-\int_r^\infty\frac{Q(r')}{r'^{2}}\dd r' \label{eq:14-vh-step2}\\ \Longrightarrow\quad v_{\mathrm{H}}(r)&=\int_r^\infty\frac{Q(r')}{r'^{2}}\dd r' \label{eq:14-vh-step3} \end{align}

を得る．この形のままでも計算できるが，二重積分になっていて非効率である．部分積分で1重積分2本に直そう．$\dd(-1/r')/\dd r'=1/r'^{2}$ を使って

\begin{align} \int_r^\infty\frac{Q(r')}{r'^{2}}\dd r' &=\left[-\frac{Q(r')}{r'}\right]_r^\infty+\int_r^\infty\frac{1}{r'}\frac{\dd Q}{\dd r'}\dd r' \label{eq:14-vh-parts1}\\ &=\frac{Q(r)}{r}+\int_r^\infty\frac{1}{r'}\cdot4\pi n(r')r'^{2}\,\dd r' \label{eq:14-vh-parts2} \end{align}

となる．1行目は部分積分の公式そのもの．2行目では，上端で $Q(\infty)/\infty\to0$(全電荷は有限)，下端から $+Q(r)/r$ が出ること，および \eqref{eq:14-vh-step1} の定義から $\dd Q/\dd r'=4\pi n(r')r'^{2}$ であることを使った．整理すると

\begin{equation} v_{\mathrm{H}}(r)=\frac{1}{r}\int_0^{r}4\pi n(r')\,r'^{2}\,\dd r' +\int_r^{\infty}4\pi n(r')\,r'\,\dd r' \label{eq:14-vh-radial} \end{equation}

物理的意味:内側は点電荷，外側は殻ごとの寄与

第1項は「半径 $r$ の内側にある全電荷 $Q(r)$ を原点に集めた点電荷が作るポテンシャル $Q(r)/r$」である．これは球殻定理 \eqref{eq:14-shell} の言い換えにほかならない．第2項は「半径 $r'\gt r$ の球殻が，その内部のどこでも作る一定のポテンシャル $\dd Q/r'=4\pi n(r')r'\dd r'$ を，$r'$ について足し上げたもの」である．一様に帯電した球殻は内部に電場を作らないが，ポテンシャルは一定値 $\dd Q/r'$ を持つ，という初等静電気学の事実がそのまま現れている．$r$ が全電荷分布の外にあれば第2項はゼロ，第1項は $Q_{\mathrm{tot}}/r$ となり，球殻定理に帰着する．

対数メッシュ上では $\dd r=r\,\dd x$ を使って

\begin{equation} v_{\mathrm{H}}(x)=\frac{4\pi}{r}\int_{x_{\min}}^{x}n(x')\,r'^{3}\,\dd x' +4\pi\int_{x}^{x_{\max}}n(x')\,r'^{2}\,\dd x' \label{eq:14-vh-logmesh} \end{equation}

と書ける．等間隔格子上の台形則やSimpson則を累積和として1回走査すれば，全格子点での $v_{\mathrm{H}}$ が $O(N)$ で得られる．

14.6.7 SCFループの全体像

部品が揃った．原子DFTコードの骨格は次の通りである．

擬似コード:原子のDFT計算

入力: 原子番号 Z, 電子配置 {(n,l,f_nl)}, 交換相関汎関数
格子: x_i = x_min + i*dx  (i = 0..N-1),  r_i = exp(x_i)

# 1. 初期密度(Thomas-Fermi近似か水素様軌道の重ね合わせ)
n(r) <- n_initial(r)

repeat:                                   # ---- SCF ループ ----
    # 2. 有効ポテンシャルの構築
    v_H(r)   <- 動径Hartree公式(14.6.6節)を累積積分で評価
    v_xc(r)  <- LDA/GGA の公式を格子点ごとに評価
    v_eff(r) <- -Z/r + v_H(r) + v_xc(r)

    # 3. 各 (n,l) について固有値問題を解く
    for each (n, l) in 電子配置:
        # 3a. 固有値を粗くスキャンし，節の数 = n-l-1 となる区間を括る
        [e_lo, e_hi] <- bracket_by_node_count(v_eff, l, n-l-1)
        # 3b. 二分法で D(eps) の符号変化を追い込む
        eps <- bisection(D, e_lo, e_hi, tol = 1e-6)
        # 3c. Newton 仕上げ(1次補正公式)
        repeat:
            u_out <- integrate_outward (v_eff, l, eps, r_min .. r_mp)
            u_in  <- integrate_inward  (v_eff, l, eps, r_max .. r_mp)
            scale u_in so that u_in(r_mp) = u_out(r_mp)
            d_eps <- u_m (u_out' - u_in')|_r_mp / (2 ∫ u² dr)
            eps   <- eps + d_eps
        until |d_eps| < 1e-12
        u_nl <- normalize(u_out ∪ u_in)     # ∫u² dr = 1

    # 4. 新しい密度
    n_out(r) <- Σ_nl f_nl * u_nl(r)² / (4 pi r²)

    # 5. 収束判定と混合
    if  ∫|n_out - n| 4πr² dr < δ : break
    n(r) <- (1-α) n(r) + α n_out(r)        # 単純混合で十分(原子は小さい)

# 6. 全エネルギー(二重数え補正つき，第10章)
E_tot = Σ_nl f_nl eps_nl
        - (1/2)∫ n v_H 4πr² dr
        + ∫ n eps_xc(n) 4πr² dr
        - ∫ n v_xc(n) 4πr² dr

手順6の全エネルギーは，第10章で導いた「固有値の和から二重数えを差し引く」形である．念のため書き下すと

\begin{equation} E_{\mathrm{tot}} =\sum_{nl}f_{nl}\,\varepsilon_{nl} -\frac{1}{2}\int\dd^3r\,n\,v_{\mathrm{H}} +E_{xc}[n] -\int\dd^3r\,n\,v_{xc} \label{eq:14-etot-atom} \end{equation}

である．第1項は $\sum f\varepsilon=T_s+\int n v_{\mathrm{eff}}$ を含んでおり，$\int nv_{\mathrm{H}}$ を丸ごと数えてしまっている(Hartreeエネルギーは $\frac{1}{2}\int nv_{\mathrm{H}}$ なので1回分多い)ため第2項で $\frac{1}{2}\int nv_{\mathrm{H}}$ を引き，同様に $\int nv_{xc}$ を引いて正しい $E_{xc}$ を足し直している．

例:水素原子でコードを検証する

$Z=1$，電子1個，Hartree項と交換相関項を強制的にゼロにすれば，\eqref{eq:14-radial-u} は水素原子の厳密解を持つ:$\varepsilon_{nl}=-1/(2n^{2})$ Ha，基底状態は $u_{1s}(r)=2re^{-r}$ である．実際に代入して検算しよう．$u=2re^{-r}$ より $u'=2e^{-r}-2re^{-r}$，$u''=-2e^{-r}-2e^{-r}+2re^{-r}=-4e^{-r}+2re^{-r}$ である．$l=0$，$v_{\mathrm{eff}}=-1/r$ での \eqref{eq:14-radial-u} の左辺は $$-\frac{1}{2}\left(-4e^{-r}+2re^{-r}\right)-\frac{1}{r}\cdot2re^{-r} =2e^{-r}-re^{-r}-2e^{-r}=-re^{-r}=-\frac{1}{2}\,u$$ となり，確かに $\varepsilon=-1/2$ Ha の固有関数である．節の数は $0$ で \eqref{eq:14-nodes} の $n-l-1=0$ に一致する．自作コードはまずこれを $10^{-10}$ Ha の精度で再現できなければならない．次にHartree項を有効にしてHe原子を計算し，LDAで $E_{\mathrm{tot}}\approx-2.83$ Ha(実験値 $-2.90$ Ha)が出れば，Poissonソルバと交換相関ルーチンも正しく動いている．

14.7 コラム:計算の再現性とΔゲージ

ここまで，基底の選択，全エネルギーの再編成，SCFの収束加速，原子ソルバと，実装の自由度を見てきた．ここで自然に湧く疑問がある．同じ汎関数を使っていれば，どのコードで計算しても同じ答えが出るのだろうか．

答えは「原理的にはイエス，しかし現実には長らくノーだった」である．$E_{xc}$ の形(たとえばPBE)を固定しても，実際に得られる格子定数や体積弾性率は，基底(平面波か局在基底か，そのカットオフ)，擬ポテンシャル(あるいは全電子か)，実空間積分の格子，Brillouin域積分の $\kk$ 点数といった実装上の選択に依存する．1990年代から2010年代初頭まで，文献に報告されるSiのPBE格子定数は $0.02$ Å 程度ばらついていた．これは「PBEという理論の予言」と呼べるものが存在しなかったということである．

注意:precision と accuracy を分けて考える

この議論では2種類の誤差を厳密に区別しなければならない．

accuracy(正確さ):計算値が実験値からどれだけずれているか．これは主として汎関数の近似($E_{xc}$ の形)に由来する．第11章で見たように，PBEは格子定数を1%程度過大評価する傾向がある．
precision(精密さ，再現性):同じ汎関数を使った別々の計算どうしがどれだけ一致するか．これは実装(基底・擬ポテンシャル・数値積分)に由来する．

汎関数の改良を議論するには，まずprecisionが確保されていなければならない．実装誤差が汎関数間の差より大きければ，「PBEとSCANのどちらが良いか」という問い自体が意味をなさないからである．

14.7.1 状態方程式を合わせるという発想

コード間の一致を定量化するには，比較する「量」を決めなければならない．格子定数だけを比べるのでは不十分である．同じ平衡体積を与えても，その周りでのエネルギー曲線の曲がり具合(体積弾性率)が違えば，圧力下の振る舞いや弾性的性質は違ってしまう．そこで，単一の数値ではなくエネルギー–体積曲線 $E(V)$ 全体を比較対象にする．

$E(V)$ は少数の物理パラメータで表せる．標準的なのが3次のBirch–Murnaghan状態方程式である:

\begin{equation} E(V)=E_0+\frac{9V_0B_0}{16}\left\{ \left[\left(\frac{V_0}{V}\right)^{2/3}\!\!-1\right]^{3}\!\!B_0' +\left[\left(\frac{V_0}{V}\right)^{2/3}\!\!-1\right]^{2}\left[6-4\left(\frac{V_0}{V}\right)^{2/3}\right]\right\} \label{eq:14-birch} \end{equation}

ここで $E_0$ は平衡エネルギー，$V_0$ は平衡体積，$B_0=-V(\partial P/\partial V)_{V_0}$ は体積弾性率，$B_0'=(\partial B/\partial P)_{V_0}$ はその圧力微分である．いくつかの体積で全エネルギーを計算し，この4パラメータをフィットすれば，$E(V)$ 曲線が滑らかな関数として得られる．$E_0$ の絶対値はコードによって基準が違う(擬ポテンシャルか全電子かで内殻の寄与が異なる)ので，比較の前に各曲線をその最小値がゼロになるようシフトしておく．

14.7.2 Δ値の定義

2つのコード $a$，$b$ の曲線 $E_a(V)$，$E_b(V)$ の食い違いを1つの数にまとめる．自然な選択は，着目する体積範囲での二乗平均平方根である:

\begin{equation} \Delta_{ab}=\sqrt{\frac{1}{V_{\max}-V_{\min}}\int_{V_{\min}}^{V_{\max}} \Big[E_a(V)-E_b(V)\Big]^{2}\,\dd V} \label{eq:14-delta-gauge} \end{equation}

積分範囲は通常 $V_{\min}=0.94\,V_0$，$V_{\max}=1.06\,V_0$($\pm6$%)にとる．この幅は，常温常圧付近の熱膨張や中程度の圧力で実際に到達する体積範囲におおむね対応している．$\Delta_{ab}$ は原子あたりのエネルギー(meV/atom)の次元を持ち，値が小さいほど2つの計算はよく一致していることになる．

Δゲージの考え方．(a) 2つのコードの E(V) 曲線（藍が E_a，赤が E_b）を最小値がゼロになるようにシフトして重ね，0.94V0 から 1.06V0 の範囲に帯をかけた図．(b) 差の二乗 [E_a − E_b]² を同じ範囲で積分する様子で，その面積を幅で割って平方根をとったものが Δ である． — 図14.8 Δゲージの考え方．(a) 2つのコードの $E(V)$ 曲線をそれぞれ最小値がゼロになるようにシフトして重ねる．(b) 差の二乗を平衡体積の $\pm6$% の範囲で積分し，区間幅で割って平方根をとった値が $\Delta_{ab}$ である．平衡体積のずれだけでなく，曲率(体積弾性率)のずれも同時に評価できる．

14.7.3 何が分かったか

この指標を使った大規模な検証が2016年に報告された [Lejaeghere et al. 2016]．15の主要コード，71元素の単体結晶，GGA-PBE汎関数，スカラー相対論という条件を揃えて $E(V)$ 曲線を計算し，全ペアについてΔ値を評価したものである．結果の要点は3つある．

全電子コード同士は $0.3$〜$1.0$ meV/atom で一致する．(L)APW+lo など独立に開発された全電子法が互いにこの精度で一致したことにより，「PBEの答え」を定義する基準が確立した．
丁寧に構成された擬ポテンシャル/PAWのセットは，全電子基準に対して $1$〜$2$ meV/atom で一致する．つまり内殻を有効ポテンシャルに畳み込む近似は，正しく作れば実用上の誤差を生まない．一方，古い世代の擬ポテンシャルライブラリでは $10$ meV/atom を超えるものもあった．差は理論ではなく，擬ポテンシャルとその基底の作り込みにある．
PBEと実験の差(accuracy)は平均 $23$ meV/atom 程度である．これは実装のばらつき($1$〜$2$ meV/atom)の10倍以上大きい．

物理的意味:Δゲージが与えたもの

最後の点が本質的である．実装誤差が汎関数誤差より一桁以上小さいことが確認されて初めて，「PBEはこの物質の格子定数を1.2%過大評価する」という言明が，特定のコードの癖ではなく汎関数そのものの性質として意味を持つ．Δゲージは，DFTを「コードごとの職人芸」から「再現可能な予測科学」へ移すための計測器だったのである．
実務的な含意も大きい．新しい擬ポテンシャルライブラリや基底セットを公開するときには，全電子基準とのΔ値を示すことが事実上の標準になった．14.3.4節で述べた「局在基底には $E_{\mathrm{cut}}$ のような単一の収束パラメータがない」という弱点も，Δゲージという客観的な合格基準があれば，基底ライブラリの品質保証という形で補うことができる．読者が自分の計算を報告するときにも，収束テスト(カットオフ，$\kk$ 点，格子)の記載と併せて，使った擬ポテンシャル・基底の出自を明記する習慣を持ってほしい．

14.8 まとめと演習

14.8.1 本章のまとめ

KS方程式は非線形固有値問題であり，内側の線形固有値問題(基底展開+対角化)と外側の密度の不動点反復(SCF)という二重構造で解かれる(14.1節，図14.1)．
波動関数を有限基底 $\{\chi_\mu\}$ で展開しRayleigh商を停留化すると，一般化固有値問題 $H\bm{c}_i=\varepsilon_i S\bm{c}_i$ \eqref{eq:14-gevp} が得られる．$S$ は正定値なので $S^{-1/2}$ による標準形への変換が可能で，その最小固有値が小さいことが「基底の張りすぎ」による数値不安定の原因である(14.2節)．
平面波はカットオフ $E_{\mathrm{cut}}$ ただ1つで系統的に収束し，行列を作らずFFTで $H\bm{c}$ を計算できる．局在基底は原子あたり $10$–$50$ 個と少なく，有限カットオフ半径により $H$，$S$ が疎になる(14.3.1節)．
局在基底は閉じ込めポテンシャル \eqref{eq:14-confine} のもとで原子問題を解いて生成する．裾を切るのではなく，はじめから短距離な関数を作るのが要点である．multiple-$\zeta$ は軌道の伸縮を，分極軌道は角度方向の変形を担う．後者が $l+1$ でなければならないことは1次摂動論の選択則から従う(14.3.2〜14.3.3節)．
全エネルギーの3つの長距離クーロン項 $E_{\mathrm{ec}}^{(\mathrm{L})}$，$E_{\mathrm{ee}}$，$E_{\mathrm{cc}}$ は個別に発散し，和だけが有限である．中性原子ポテンシャル $V_{\mathrm{na},I}=V_{\mathrm{core},I}+v_{\mathrm{H},I}^{(a)}$ と差電荷 $\delta n$ による再編成 $E_{\mathrm{ec}}^{(\mathrm{L})}+E_{\mathrm{ee}}+E_{\mathrm{cc}}=E_{\mathrm{na}}+E_{\delta\mathrm{ee}}+E_{\mathrm{scc}}$ は厳密な恒等式であり，これにより発散も桁落ちも消え，Ewald和も不要になる(14.4節)．
差電荷は $\int\delta n\,\dd^3r=0$ を満たすので，そのPoisson方程式は逆空間で $\bm{G}=0$ の困難なく解ける(14.4.4節)．
SCFの単純混合の誤差増幅率は $1-\alpha\varepsilon(\bm{q})$ であり，収束条件は $\alpha\lt2/\varepsilon(q_{\min})$．金属では $\varepsilon(q)\simeq1+\kappa_{\mathrm{TF}}^{2}/q^{2}$ のため $\alpha\propto L^{-2}$ となる．これが charge sloshing である(14.5.1〜14.5.2節)．
Kerker混合の因子 $w(q)=q^{2}/(q^{2}+q_0^{2})$ はモデル誘電関数の逆数であり，増幅率を波数に依らず $1-\alpha$ に均す．Pulay(DIIS)混合は $\sum_mc_m=1$ の拘束下で予測残差 $\sum_mc_mR_m$ のノルムを最小化する係数を求める準Newton法であり，Lagrange未定乗数法により縁付き連立一次方程式 \eqref{eq:14-diis-eq} に帰着する(14.5.3〜14.5.4節)．
球対称原子ではKS方程式が変数分離でき，$u=rR$ の置換により1次元Schrödinger方程式 \eqref{eq:14-radial-u} になる．対数メッシュ $x=\ln r$ は相対分解能を全域で一定に保つ．シューティング法では対数微分の連続性 $D(\varepsilon)=0$ を条件とし，ノード数で状態を識別し，二分法と1次補正 \eqref{eq:14-eps-correction} で固有値を求める(14.6節)．
球対称密度のHartreeポテンシャルは $v_{\mathrm{H}}(r)=\frac{1}{r}\int_0^r4\pi n r'^{2}\dd r'+\int_r^\infty4\pi n r'\dd r'$ という1重積分2本で書ける(14.6.6節)．
コード間の再現性はΔゲージ \eqref{eq:14-delta-gauge} で定量化される．実装誤差($1$–$2$ meV/atom)が汎関数誤差($\sim23$ meV/atom)より一桁小さいことが確認され，DFTは再現可能な予測科学になった(14.7節)．

次章への橋渡し．本章では擬ポテンシャルを「外場の局所部分が遠方で $-Z_I/r$ になり，加えて非局所項がある」というブラックボックスとして扱ってきた．しかし14.3節の擬原子軌道も，14.4節の中性原子ポテンシャルも，その中身に決定的に依存している．第15章では，この箱を開ける．内殻電子を消し去りながら価電子の散乱特性を保存するという要請が，ノルム保存条件という具体的な数学的条件に翻訳される過程を追う．そこで使う道具は，本章14.6節で作った原子ソルバそのものである．

14.8.2 演習問題

演習14.1 基底の大きさを見積もる

1辺 $a=8.0$ Bohr の立方体セルに，厚さ $4.0$ Bohr のスラブを置き，残り $4.0$ Bohr を真空とする(セルは $8\times8\times8$ Bohr$^{3}$)．カットオフ $E_{\mathrm{cut}}=20$ Ha を用いる．

\eqref{eq:14-pw-count} を使って平面波の数 $N_{\mathrm{PW}}$ を求めよ．
\eqref{eq:14-grid-spacing} からFFT格子の間隔 $h$ の上限と，必要な格子点数(各方向 $2^{k}$ に切り上げなくてよい)を求めよ．
真空層を $12$ Bohr に増やしたとき，$N_{\mathrm{PW}}$ と格子点数はそれぞれ何倍になるか．局在基底ならこの増加はどうなるか，理由とともに述べよ．

ヒント:(1) $q_c=\sqrt{2E_{\mathrm{cut}}}$．(3) $N_{\mathrm{PW}}\propto\Omega$ であることと，局在基底の基底数が原子数だけで決まることを比較する．

演習14.2 混合パラメータの上限とKerker前処理

$\kappa_{\mathrm{TF}}=1.0$ Bohr$^{-1}$ の金属を，1辺 $L$ の立方セルで計算する．

$L=15,\,30,\,60$ Bohr のそれぞれについて，単純混合が収束するための $\alpha$ の上限 \eqref{eq:14-alpha-eps} を求めよ．
$L=60$ Bohr，$\alpha$ を上限の $0.8$ 倍にとったとき，短波長成分($\varepsilon\approx1$)の誤差を $10^{-5}$ にするのに必要な反復回数を見積もれ．
Kerker因子を $q_0=\kappa_{\mathrm{TF}}$ で導入すると，(2) の反復回数はどうなるか．$\alpha=0.3$ が使えるとして計算せよ．
次に，履歴を2個($R_1,R_2$)だけ使うDIISを考える．拘束 $c_1+c_2=1$ のもとで $\left\|c_1R_1+c_2R_2\right\|^{2}$ を最小化する係数が $$c_1=\frac{\braket{R_2|R_2-R_1}}{\left\|R_2-R_1\right\|^{2}},\qquad c_2=1-c_1$$ となることを，$c_2=1-c_1$ を代入して $c_1$ の1変数2次関数を最小化する方法で示せ．同じ結果が \eqref{eq:14-diis-eq} の $2\times2$ 版を解いても得られることも確かめよ．
密度が1変数 $n$ の場合を考え，(4) の係数から作られる $\bar n_{\mathrm{in}}=c_1n^{(1)}+c_2n^{(2)}$ が，方程式 $R(n)=0$ に対する割線法(secant method)の1ステップと一致することを示せ．

ヒント:誤差は反復ごとに $\abs{1-\alpha\varepsilon(q)}$ 倍される．必要回数は $n=\ln(10^{-5})/\ln\abs{1-\alpha\varepsilon}$．(3) では \eqref{eq:14-kerker-check2} により増幅率が $1-\alpha$ で波数に依らないことを使う．(4) は $\left\|R_2+c_1(R_1-R_2)\right\|^{2}$ を展開して $c_1$ で微分する．(5) 1変数では $R_1,R_2$ は数であり，割線法の更新は $n^{(3)}=n^{(2)}-R_2(n^{(2)}-n^{(1)})/(R_2-R_1)$ である．

演習14.3 動径方程式と動径Hartree公式の検算

水素様原子(核電荷 $Z$)の $2p$ 状態 $u_{2p}(r)=C\,r^{2}e^{-Zr/2}$ を \eqref{eq:14-radial-u} に代入し，$\varepsilon_{2p}=-Z^{2}/8$ となることを示せ．また節の数が \eqref{eq:14-nodes} と整合することを確かめよ．
半径 $R$ の球内に電荷 $Q$ が一様に分布した密度 $n(r)=3Q/(4\pi R^{3})\ (r\le R)$，$0\ (r\gt R)$ について，\eqref{eq:14-vh-radial} の両方の項を実行し $$v_{\mathrm{H}}(r)=\frac{Q}{2R}\left(3-\frac{r^{2}}{R^{2}}\right)\ (r\le R),\qquad v_{\mathrm{H}}(r)=\frac{Q}{r}\ (r\gt R)$$ を導け．$r=R$ で値と1階微分が連続であることも確認せよ．
この結果が球殻定理 \eqref{eq:14-shell} と矛盾しないことを説明せよ．

ヒント:(1) $u''$ を丁寧に計算し，$l=1$ の遠心力項 $1/r^{2}$ と $-Z/r$ を合わせる．$r^{-2}$，$r^{-1}$，定数の各冪の係数が別々に釣り合う．(2) 第1項は $r\le R$ で $Q(r)=Q(r/R)^{3}$，第2項は $\int_r^R$ の積分になる．

演習14.4 中性原子分解の恒等式を確かめる

原子が2個(位置 $\bm{\tau}_1,\bm{\tau}_2$，価電荷 $Z_1,Z_2$)しかない系で，14.4.3節の恒等式を最初から書き下してみよ．$\delta n=0$(密度が原子密度の重ね合わせに完全に等しい)という仮想的な場合に，$E_{\mathrm{ec}}^{(\mathrm{L})}+E_{\mathrm{ee}}+E_{\mathrm{cc}}$ が $E_{\mathrm{na}}+E_{\mathrm{scc}}$ に等しくなることを，\eqref{eq:14-step2}〜\eqref{eq:14-escc-def} の各ステップをたどって確認せよ．さらに2原子が十分離れているとき，この値が「孤立原子2個分のエネルギーの和」になること(相互作用がゼロになること)を示せ．

ヒント:$\delta n=0$ のとき $E_{\delta\mathrm{ee}}=0$．十分離れているときは \eqref{eq:14-uij-far2} により $E_{\mathrm{scc}}$ の $I\neq J$ 項が消え，$E_{\mathrm{na}}$ も $V_{\mathrm{na},I}$ の短距離性 \eqref{eq:14-vna-short} により各原子の自己項だけになる．

参考文献

教科書

R. M. Martin, Electronic Structure: Basic Theory and Practical Methods (Cambridge University Press, 2004) — 第12〜16章が本章の内容全般に対応する標準的な参考書．
J. Kohanoff, Electronic Structure Calculations for Solids and Molecules (Cambridge University Press, 2006) — 基底，SCF収束，擬ポテンシャルの実装的側面．
M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias and J. D. Joannopoulos, Rev. Mod. Phys. 64, 1045 (1992) — 平面波法と反復対角化の古典的レビュー．

変分性と基底

J. K. L. MacDonald, Phys. Rev. 43, 830 (1933) — 基底拡大に対する固有値の単調性(Hylleraas–Undheim–MacDonaldの定理)．
S. F. Boys and F. Bernardi, Mol. Phys. 19, 553 (1970) — counterpoise法によるBSSE補正．
O. F. Sankey and D. J. Niklewski, Phys. Rev. B 40, 3979 (1989) — 閉じ込めによる有限台の数値原子軌道．
J. Junquera, Ó. Paz, D. Sánchez-Portal and E. Artacho, Phys. Rev. B 64, 235111 (2001) — 局在基底の閉じ込め方と系統的な基底生成．
J. M. Soler et al., J. Phys.: Condens. Matter 14, 2745 (2002) — 数値原子軌道による $O(N)$ 実装(SIESTA)．
T. Ozaki, Phys. Rev. B 67, 155108 (2003) — 基底関数の変分最適化．
T. Ozaki and H. Kino, Phys. Rev. B 69, 195113 (2004) — H から Kr までの数値原子基底の系統的構成と検証．

全エネルギーとPoisson求解

T. Ozaki and H. Kino, Phys. Rev. B 72, 045121 (2005) — 中性原子ポテンシャルと差電荷によるクーロン項の再編成．
P. P. Ewald, Ann. Phys. 369, 253 (1921) — 周期系のクーロン和．

SCF収束と電荷混合

G. P. Kerker, Phys. Rev. B 23, 3082 (1981) — 長波長成分の前処理(Kerker混合)．
P. Pulay, Chem. Phys. Lett. 73, 393 (1980); J. Comput. Chem. 3, 556 (1982) — DIIS法．
D. G. Anderson, J. ACM 12, 547 (1965) — Anderson混合(DIISと等価な定式化)．
D. D. Johnson, Phys. Rev. B 38, 12807 (1988) — 修正Broyden法による密度混合．
G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996) — Kerker計量つきRMM-DIISと実装上の詳細．
N. D. Mermin, Phys. Rev. 137, A1441 (1965) — 有限温度DFTと自由エネルギー汎関数(smearingの理論的根拠)．
M. Methfessel and A. T. Paxton, Phys. Rev. B 40, 3616 (1989) — 高次のsmearing関数．
M. J. Gillan, J. Phys.: Condens. Matter 1, 689 (1989) — $\sigma\to0$ 外挿と $E_0\simeq(E+\Omega)/2$ の関係．
E. R. Davidson, J. Comput. Phys. 17, 87 (1975) — 反復対角化法．

原子のDFT計算

D. R. Hartree, The Calculation of Atomic Structures (Wiley, 1957) — 動径方程式の数値解法の古典．
S. E. Koonin and D. C. Meredith, Computational Physics (Addison-Wesley, 1990) — シューティング法とNumerov法の実装．
W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical Recipes, 3rd ed. (Cambridge University Press, 2007) — 常微分方程式の積分と根の探索．

再現性

F. Birch, Phys. Rev. 71, 809 (1947) — Birch–Murnaghan状態方程式．
K. Lejaeghere et al., Science 351, aad3000 (2016) — 固体のDFT計算の再現性検証とΔゲージ．
K. Lejaeghere, V. Van Speybroeck, G. Van Oost and S. Cottenier, Crit. Rev. Solid State Mater. Sci. 39, 1 (2014) — Δ値の定義とその統計的性質．