第11章交換相関汎関数
前章でコーン・シャム(KS)方程式を導出したとき，多電子問題の困難はすべて交換相関エネルギー $E_{xc}[\rho]$ という一つの汎関数に押し込められた．形式は厳密であり，方程式の形も全エネルギーの表式も自己無撞着(SCF)の手続きも，$E_{xc}$ の中身を問わずに確定した．しかし $E_{xc}[\rho]$ の厳密な形は誰も知らない——知っていれば多体問題は解けている．したがって第一原理計算の精度は，この最後の未知項をどう近似するかだけで決まる．本章ではまず最も素朴な近似である局所密度近似(LDA)を構成し，続いて「なぜこんな乱暴な近似がこれほどうまくいくのか」という問いに，断熱接続と交換相関ホールという概念で答える．この節が本章の理論的核心であり，以後の汎関数開発のすべてがここに立脚している．さらに厳密汎関数が満たすべき数学的制約(スケーリング則・Lieb–Oxford下限・自己相互作用の消去)を導き，それらを「試験問題」として設計された一般化勾配近似(GGA)の代表であるPBE汎関数を，経験的パラメータを一つも使わずに組み立てる．最後に，ハイブリッド汎関数・meta-GGA・van der Waals汎関数・DFT+U へと続く「ヤコブの梯子」を概観し，KS法の変分的性質の応用としてハリス汎関数を扱う．第12章以降で固体の電子状態を計算するとき，実際に手を動かして選ぶのは本章で登場する汎関数の名前である．
- 局所密度近似 $E_{xc}^{\mathrm{LDA}}[\rho]=\int\rho\,\varepsilon_{xc}(\rho)\dd^3r$ の構成と，交換相関ポテンシャル $v_{xc}=\varepsilon_{xc}+\rho\,\dd\varepsilon_{xc}/\dd\rho$ の導出，スピン分極版(LSDA)
- 断熱接続(結合定数積分)の完全な導出:Hellmann–Feynman定理を $\lambda$ 微分に適用して $E_{xc}=\int_0^1\dd\lambda\braket{\Psi_\lambda|\hat V_{ee}|\Psi_\lambda}-E_{\mathrm{H}}$ を得る
- 交換相関ホール $\bar n_{xc}(\rr,\rr')$ による $E_{xc}$ の表現，総和則 $\int n_x\dd^3r'=-1$，$\int n_c\dd^3r'=0$，そして $E_{xc}$ がホールの球対称平均にしか依存しないことの導出——これがLDAの成功の理由である
- 厳密な制約:一様スケーリング $E_x[\rho_\gamma]=\gamma E_x[\rho]$ の証明，相関のスケーリング不等式，Lieb–Oxford下限，Levy–Perdewのビリアル型関係式，自己相互作用の消去条件
- 勾配展開(GEA)の失敗と実空間カットオフによるGGAの再生，無次元換算勾配 $s$，PBEの交換増強因子 $F_x(s)=1+\kappa-\kappa/(1+\mu s^2/\kappa)$ の各パラメータが厳密条件から決まること
- LDAとGGAの数値的比較(原子の交換相関エネルギー，分子の過結合，固体の格子定数，bcc鉄の基底状態)とGGAが失敗する例
- ヤコブの梯子:meta-GGA，ハイブリッド汎関数(断熱接続からの混合比の根拠)，HSE，vdW-DFとDFT-D，DFT+U
- ハリス汎関数 $E_{\mathrm{Harris}}[\rho_{\mathrm{in}}]$ の誤差が $\delta\rho$ の2次であることの証明と，非自己無撞着計算・局所力定理
本章でもハートリー原子単位系($\hbar=m_e=e=4\pi\varepsilon_0=1$，第1章参照)を用いる．エネルギーの単位はハートリー(Ha)で，$1\ \mathrm{Ha}=27.2114\ \mathrm{eV}$，長さの単位はボーア半径 $a_0=0.529177\ \text{Å}$ である．固体物性の文献では格子定数を Å，体積弾性率を GPa で表すのが慣習なので，数値比較の節ではそれらを併記する．
11.1 局所密度近似(LDA)
これから，第7章で厳密に(あるいは数値的に極めて高精度に)求めた一様電子ガスの交換相関エネルギーを，不均一系に「局所的に」持ち込むという近似を構成する．得られるのは $E_{xc}$ の具体形と，KS方程式に入れるべき交換相関ポテンシャル $v_{xc}(\rr)$ の閉じた表式である．
11.1.1 発想:密度を局所的に一様電子ガスで置き換える
第7章で扱ったジェリウム(一様電子ガス)は，密度 $\rho$ が空間的に一定という理想化された系であった．この系については，電子1個あたりの交換相関エネルギー
$$ \begin{equation} \varepsilon_{xc}(\rho)=\varepsilon_{x}(\rho)+\varepsilon_{c}(\rho) \label{eq:11-eps-xc-def} \end{equation} $$が，密度 $\rho$ というただ一つの数の関数として分かっている．交換部分 $\varepsilon_x$ は解析的に，相関部分 $\varepsilon_c$ は量子モンテカルロ法による数値計算とその解析的なフィットとして与えられる(11.1.3節)．
いま，実在の原子・分子・固体のように密度が場所によって変わる系を考える．もし密度の空間変化が十分ゆっくりであれば，点 $\rr$ のまわりの小さな体積要素 $\Delta V$ の中では密度はほぼ一定 $\rho(\rr)$ とみなせる．そこでこの体積要素の中の電子は「密度 $\rho(\rr)$ の一様電子ガスの中にいる」と近似し，その交換相関エネルギーを
$$ \begin{equation} \Delta E_{xc}\simeq \underbrace{\rho(\rr)\Delta V}_{\text{この体積中の電子数}}\times\underbrace{\varepsilon_{xc}\bigl(\rho(\rr)\bigr)}_{\text{電子1個あたり}} \label{eq:11-lda-cell} \end{equation} $$と見積もる．系全体ではこれを足し上げればよく，$\Delta V\to0$ の極限で積分になる．
定義:局所密度近似(LDA: Local Density Approximation)
$$ \begin{equation} E_{xc}^{\mathrm{LDA}}[\rho]=\int \dd^3r\;\rho(\rr)\,\varepsilon_{xc}\bigl(\rho(\rr)\bigr) \label{eq:11-lda-def} \end{equation} $$ここで $\varepsilon_{xc}(\bar\rho)$ は密度 $\bar\rho$ の一様電子ガスの，電子1個あたりの交換相関エネルギーである．$\varepsilon_{xc}$ は「関数の関数」ではなく1変数の普通の関数であることに注意せよ．汎関数としての非自明さは，その引数に $\rho(\rr)$ を代入して積分するところにしかない．
注意:「局所」の意味と，$\varepsilon_{xc}$ と「エネルギー密度」の区別
$E_{xc}^{\mathrm{LDA}}$ は局所汎関数である．すなわち点 $\rr$ における被積分関数は，その点の密度 $\rho(\rr)$ だけで決まり，隣の点の密度にも密度の勾配にも依存しない．これに対し，たとえばハートリーエネルギー $E_{\mathrm{H}}$ は $\rho(\rr)$ と $\rho(\rr')$ を掛け合わせて二重積分するので非局所汎関数である．真の $E_{xc}$ も本来は非局所である(11.2節で見るように，$\rr$ の電子は $\rr'$ に他の電子が来ないよう「穴」を掘っており，その穴は有限の広がりを持つ)．LDAはその非局所性をまるごと捨てる近似である．
また，$\varepsilon_{xc}(\rho)$ は電子1個あたりのエネルギー(単位:Ha)であり，単位体積あたりのエネルギー(エネルギー密度，単位:Ha/$a_0^3$)は $\rho\,\varepsilon_{xc}(\rho)$ である．文献によっては後者を $\varepsilon_{xc}$ と書くこともあるので，式の形($\rho$ が掛かっているかどうか)で見分ける習慣をつけるとよい．
11.1.2 交換部分:ディラックの交換エネルギー密度
まず交換部分を書き下す．第7章で，一様電子ガスの電子1個あたりの交換エネルギーが
$$ \begin{equation} \varepsilon_x(\rho)=-\frac{3k_F}{4\pi}, \qquad k_F=\bigl(3\pi^2\rho\bigr)^{1/3} \label{eq:11-eps-x-kf} \end{equation} $$と求まっていた($k_F$ はフェルミ波数，第3章)．これを密度だけの式に直す．$k_F$ を代入して
\begin{align} \varepsilon_x(\rho) &=-\frac{3}{4\pi}\bigl(3\pi^2\rho\bigr)^{1/3} \label{eq:11-epsx-1}\\ &=-\frac{3}{4}\left(\frac{3\pi^2}{\pi^3}\right)^{1/3}\rho^{1/3} \label{eq:11-epsx-2}\\ &=-\frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3}\rho^{1/3} \;\equiv\;-C_x\,\rho^{1/3}. \label{eq:11-epsx-3} \end{align}\eqref{eq:11-epsx-1} から \eqref{eq:11-epsx-2} へは，前の因子 $1/\pi$ を3乗根の中に入れた($1/\pi=(1/\pi^3)^{1/3}$)．\eqref{eq:11-epsx-3} では $3\pi^2/\pi^3=3/\pi$ と約分した．定数は
$$ \begin{equation} C_x=\frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3}=0.738559\ldots \label{eq:11-cx} \end{equation} $$である．したがってLDAの交換エネルギーは
$$ \begin{equation} E_x^{\mathrm{LDA}}[\rho]=\int \dd^3r\,\rho(\rr)\varepsilon_x(\rho(\rr)) =-C_x\int \dd^3r\;\rho^{4/3}(\rr). \label{eq:11-ex-lda} \end{equation} $$この形は歴史的にはDiracが1930年にThomas–Fermi模型の補正として導いたもので，ディラック交換あるいは(係数を可変にした場合)Slaterの $X\alpha$ 交換と呼ばれる．$\rho^{4/3}$ という指数は，$\varepsilon_x\propto\rho^{1/3}$($\propto k_F\propto$ 逆長さ)に密度 $\rho$ を掛けたことから来ている．
物理的意味:なぜ $\rho^{1/3}$ か
交換エネルギーは「同スピンの電子が互いを避けることによるクーロンエネルギーの利得」である．その大きさはおよそ(1電子あたり)$-1/(\text{電子間の典型距離})$ のオーダーであり，電子間距離は $\rho^{-1/3}$ に比例する．したがって $\varepsilon_x\sim-\rho^{1/3}$ となる．式 \eqref{eq:11-epsx-3} の係数 $-\tfrac34(3/\pi)^{1/3}$ は，この次元解析的な見積もりをフェルミ球の計算で厳密化した結果である．密度パラメータ $r_s$(半径 $r_s$ の球が電子1個分の体積を持つ，第3章)を使うと $\rho=3/(4\pi r_s^3)$ だから
$$ \varepsilon_x=-\frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3}\left(\frac{3}{4\pi}\right)^{1/3}\frac{1}{r_s} =-\frac{0.458165}{r_s}\ \mathrm{Ha} $$となり，「電子間距離の逆数」という描像がそのまま見える．
11.1.3 相関部分:量子モンテカルロとそのフィット
相関エネルギー $\varepsilon_c(\rho)$ には閉じた解析形がない．第7章で見たように，高密度極限($r_s\to0$)と低密度極限($r_s\to\infty$)では摂動論的・古典的に漸近形が分かる．以下はすべてハートリー単位である——第7章では同じ結果をリュードベリ単位で($\varepsilon_c\simeq0.0622\ln r_s-0.094$ Ry のように)書いたので，係数はちょうど半分になっていることに注意してほしい:
\begin{align} \varepsilon_c(r_s)&\to A\ln r_s+B+O(r_s\ln r_s), \qquad A=\frac{1-\ln 2}{\pi^2}=0.031091,\quad B\simeq-0.0469 &&(r_s\to0)\quad[\mathrm{Ha}] \label{eq:11-ec-high}\\ \varepsilon_c(r_s)&\to-\frac{d_0}{r_s}+\frac{d_1}{r_s^{3/2}}+\cdots, \qquad d_0\simeq0.4400,\quad d_1\simeq1.325 &&(r_s\to\infty)\quad[\mathrm{Ha}] \label{eq:11-ec-low} \end{align}高密度側 \eqref{eq:11-ec-high} はGell-Mann–Brueckner のRPA(乱雑位相近似)による結果，低密度側 \eqref{eq:11-ec-low} は電子が結晶化する(ウィグナー結晶)極限のマーデルング的評価である(第7章のウィグナー結晶の全エネルギー $E/N\simeq-1.792/r_s+2.66/r_s^{3/2}$ Ry から交換項 $-0.916/r_s$ Ry を差し引いてHaに直すと，$d_0\simeq0.438$，$d_1\simeq1.33$ が得られる)．実在物質の価電子は $r_s\simeq2$〜$6$ という中間領域にあり，どちらの漸近形も使えない．
この中間領域を埋めたのがCeperleyとAlderによる量子モンテカルロ(QMC)計算である．彼らはSlater–Jastrow型の試行波動関数
$$ \begin{equation} \Psi_{\mathrm{QMC}}=\exp\!\left[-\sum_{i\lt j}^{N}u\bigl(\abs{\rr_i-\rr_j}\bigr)\right]\; \det\bigl|\varphi_1\varphi_2\cdots\varphi_N\bigr| \label{eq:11-slater-jastrow} \end{equation} $$を用いた．行列式部分は第5章のスレーター行列式(反対称性とパウリ相関を担う)であり，指数因子(Jastrow因子)は2電子が近づくほど波動関数を小さくして，行列式では表せないクーロン相関を取り込む．この試行関数を出発点に確率的に基底状態エネルギーを求め，そこから既知の運動エネルギー・交換エネルギーを差し引いて $\varepsilon_c(r_s)$ を数点で決定した．
数点の数値データをKS方程式の中で使うには，$r_s$ の滑らかな関数として補間する必要がある．代表的な補間形が2つある．
定義:相関エネルギー密度の解析的フィット
(a) Perdew–Zunger (PZ81) 形——低密度側と高密度側で式を切り替える簡潔な形．スピン無分極の場合(Ha単位):
$$ \begin{equation} \varepsilon_c(r_s)= \begin{cases} \dfrac{\gamma}{1+\beta_1\sqrt{r_s}+\beta_2 r_s}, & r_s\ge1,\\[2ex] A\ln r_s+B+C r_s\ln r_s+D r_s, & r_s\lt1, \end{cases} \label{eq:11-pz} \end{equation} $$$\gamma=-0.1423$，$\beta_1=1.0529$，$\beta_2=0.3334$，$A=0.0311$，$B=-0.048$，$C=0.0020$，$D=-0.0116$．$r_s\ge1$ 側の形は低密度極限 \eqref{eq:11-ec-low} の $-d_0/r_s$ の振る舞いを，$r_s\lt1$ 側は高密度極限 \eqref{eq:11-ec-high} の対数を，それぞれ正しく再現するように選ばれている．
(b) Vosko–Wilk–Nusair (VWN) 形——RPAの解析構造をそのまま関数形に反映させたもの．$x=\sqrt{r_s}$，$X(x)=x^2+bx+c$，$Q=\sqrt{4c-b^2}$ として
$$ \begin{equation} \varepsilon_c=A\left\{\ln\frac{x^2}{X(x)}+\frac{2b}{Q}\arctan\frac{Q}{2x+b} -\frac{bx_0}{X(x_0)}\left[\ln\frac{(x-x_0)^2}{X(x)}+\frac{2(b+2x_0)}{Q}\arctan\frac{Q}{2x+b}\right]\right\} \label{eq:11-vwn} \end{equation} $$パラメータは $A=0.0310907$，$b=3.72744$，$c=12.9352$，$x_0=-0.10498$(スピン無分極)．第1項の $\ln(x^2/X)$ が $r_s\to0$ で $\ln r_s$ の対数発散を，$\arctan$ 項が中間領域の曲がりを担う．
両者は数値的にはほぼ同じ曲線を与える(差は $10^{-4}$ Ha 程度)．より新しい Perdew–Wang (PW92) 形も広く使われるが，構造は同種である．「LDA」と書かれた計算がPZを使っているかVWNを使っているかで結果が有意に変わることはまずない．
| $r_s$ ($a_0$) | $\varepsilon_x$ (Ha) | $\varepsilon_c$ (Ha) | $\abs{\varepsilon_c/\varepsilon_x}$ |
|---|---|---|---|
| 1 | $-0.4582$ | $-0.0596$ | 13.0% |
| 2 | $-0.2291$ | $-0.0451$ | 19.7% |
| 3 | $-0.1527$ | $-0.0372$ | 24.4% |
| 4 | $-0.1145$ | $-0.0321$ | 28.0% |
| 6 | $-0.0764$ | $-0.0255$ | 33.4% |
| 10 | $-0.0458$ | $-0.0186$ | 40.5% |
11.1.4 交換相関ポテンシャル $v_{xc}^{\mathrm{LDA}}$ の導出
KS方程式に入るのは $E_{xc}$ そのものではなく，その汎関数微分 $v_{xc}(\rr)=\delta E_{xc}/\delta\rho(\rr)$ である(第10章)．LDAの形 \eqref{eq:11-lda-def} に対してこれを実行する．
導出:局所汎関数の汎関数微分
$E_{xc}^{\mathrm{LDA}}$ は
$$ \begin{equation} E_{xc}^{\mathrm{LDA}}[\rho]=\int \dd^3r\;f\bigl(\rho(\rr)\bigr), \qquad f(\rho)\equiv\rho\,\varepsilon_{xc}(\rho) \label{eq:11-f-def} \end{equation} $$という「1変数関数 $f$ を密度に代入して積分しただけ」の形をしている．密度を $\rho\to\rho+\delta\rho$ と変化させ，$\delta\rho$ について1次までを拾う:
\begin{align} \delta E_{xc}^{\mathrm{LDA}} &=\int \dd^3r\;\Bigl[f\bigl(\rho(\rr)+\delta\rho(\rr)\bigr)-f\bigl(\rho(\rr)\bigr)\Bigr] \label{eq:11-fd-1}\\ &=\int \dd^3r\;\Bigl[f'\bigl(\rho(\rr)\bigr)\,\delta\rho(\rr) +\tfrac12 f''\bigl(\rho(\rr)\bigr)\,\delta\rho(\rr)^2+\cdots\Bigr] \label{eq:11-fd-2}\\ &=\int \dd^3r\;f'\bigl(\rho(\rr)\bigr)\,\delta\rho(\rr)+O(\delta\rho^2). \label{eq:11-fd-3} \end{align}\eqref{eq:11-fd-1} は定義通りの差．\eqref{eq:11-fd-2} では，各点 $\rr$ ごとに1変数関数 $f$ を $\rho(\rr)$ のまわりでテイラー展開した(積分の中で点ごとに独立に展開できるのは，$f$ の引数がその点の密度だけだからである——これが「局所」であることの利点である)．\eqref{eq:11-fd-3} で2次以上を捨てた．汎関数微分の定義 $\delta E=\int\dd^3r\,\frac{\delta E}{\delta\rho(\rr)}\delta\rho(\rr)$(第4章)と見比べると
$$ \begin{equation} v_{xc}^{\mathrm{LDA}}(\rr)=\frac{\delta E_{xc}^{\mathrm{LDA}}}{\delta\rho(\rr)} =f'(\rho)\Big|_{\rho=\rho(\rr)} =\frac{\dd}{\dd\rho}\bigl[\rho\,\varepsilon_{xc}(\rho)\bigr]_{\rho=\rho(\rr)}. \label{eq:11-vxc-fprime} \end{equation} $$積の微分(ライプニッツ則)を実行して
$$ \begin{equation} v_{xc}^{\mathrm{LDA}}(\rr)=\varepsilon_{xc}\bigl(\rho(\rr)\bigr) +\rho(\rr)\,\frac{\dd\varepsilon_{xc}}{\dd\rho}\bigg|_{\rho=\rho(\rr)}. \label{eq:11-vxc-lda} \end{equation} $$∎
注意すべきは，$v_{xc}\ne\varepsilon_{xc}$ であることである．$\varepsilon_{xc}$ は「電子1個あたりのエネルギー」，$v_{xc}$ は「電子を1個足したときのエネルギー変化率」であり，後者には $\rho\,\dd\varepsilon_{xc}/\dd\rho$ という項——電子を足すと密度が上がり，既にいる電子の $\varepsilon_{xc}$ も変わる効果——が加わる．
交換部分については具体的に微分を実行できる．$\varepsilon_x=-C_x\rho^{1/3}$ だから $f(\rho)=\rho\varepsilon_x=-C_x\rho^{4/3}$ であり，
\begin{align} v_x^{\mathrm{LDA}}(\rr) &=\frac{\dd}{\dd\rho}\left(-C_x\rho^{4/3}\right) \label{eq:11-vx-1}\\ &=-\frac{4}{3}C_x\rho^{1/3} \label{eq:11-vx-2}\\ &=-\frac{4}{3}\cdot\frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3}\rho^{1/3} \label{eq:11-vx-3}\\ &=-\left(\frac{3}{\pi}\right)^{1/3}\rho^{1/3}(\rr) =-0.984745\,\rho^{1/3}(\rr). \label{eq:11-vx-4} \end{align}\eqref{eq:11-vx-1}→\eqref{eq:11-vx-2} は冪の微分，\eqref{eq:11-vx-3} で $C_x$ の定義 \eqref{eq:11-cx} を代入，\eqref{eq:11-vx-4} で $\tfrac43\cdot\tfrac34=1$ と約分した．同時に
$$ \begin{equation} v_x^{\mathrm{LDA}}=\frac{4}{3}\,\varepsilon_x \label{eq:11-vx-43} \end{equation} $$という簡潔な関係も得られる($\varepsilon_x=-C_x\rho^{1/3}$ の $4/3$ 倍が \eqref{eq:11-vx-2} だから)．交換ポテンシャルはエネルギー密度より $33\%$ 深い．
例:$X\alpha$ 法との関係
Slaterは1951年，ハートリー・フォック法の非局所交換演算子を「フェルミ球平均」で置き換えて $v_x^{\mathrm{Slater}}=-\tfrac32(3\rho/\pi)^{1/3}$ を得た．これは \eqref{eq:11-vx-4} の $3/2$ 倍である．この食い違いは，Slaterがポテンシャルを平均したのに対し，DFTではエネルギーを先に決めてから微分することに由来する．両者を補間するために係数を可変にしたのが $X\alpha$ 法
$$ v_{x\alpha}=-\frac{3}{2}\alpha\left(\frac{3\rho}{\pi}\right)^{1/3} $$である．$\alpha=1$ がSlater，$\alpha=2/3$ がKS-LDAの交換に対応する．経験的には $\alpha\simeq0.7$ が原子の全エネルギーをよく再現し，DFT成立以前の標準的手法として広く使われた．
11.1.5 スピン分極版(LSDA)とスピン内挿
磁性体・開殻分子・奇数電子系では，上向きと下向きのスピン密度 $\rho^\uparrow,\rho^\downarrow$ を独立変数にとるスピン密度汎関数法(第10.6節)が必要である．このときLDAは局所スピン密度近似(LSDA)
$$ \begin{equation} E_{xc}^{\mathrm{LSDA}}[\rho^\uparrow,\rho^\downarrow] =\int \dd^3r\;\rho(\rr)\,\varepsilon_{xc}\bigl(\rho^\uparrow(\rr),\rho^\downarrow(\rr)\bigr), \qquad \rho=\rho^\uparrow+\rho^\downarrow \label{eq:11-lsda-def} \end{equation} $$へ拡張される．$\varepsilon_{xc}(\rho^\uparrow,\rho^\downarrow)$ は「スピン分極した一様電子ガス」の電子1個あたりの交換相関エネルギーである．慣習として，全密度 $\rho$ とスピン分極率
$$ \begin{equation} \zeta\equiv\frac{\rho^\uparrow-\rho^\downarrow}{\rho}\in[-1,1] \label{eq:11-zeta} \end{equation} $$を変数に選ぶ($\zeta=0$ が無分極，$\zeta=\pm1$ が完全分極)．対応する交換相関ポテンシャルは，スピンチャネルごとの偏微分
$$ \begin{equation} v_{xc}^{\sigma}(\rr)=\frac{\delta E_{xc}^{\mathrm{LSDA}}}{\delta\rho^\sigma(\rr)} =\frac{\partial}{\partial\rho^\sigma}\Bigl[\rho\,\varepsilon_{xc}(\rho^\uparrow,\rho^\downarrow)\Bigr] \label{eq:11-vxc-sigma} \end{equation} $$である(11.1.4節の導出をそのままスピンごとに繰り返せばよい)．$v_{xc}^\uparrow\ne v_{xc}^\downarrow$ となることが，KS方程式におけるスピン分裂——磁性の起源——を生む．
交換部分については，$\zeta$ 依存性が厳密に決まる．これを導く．
導出:交換のスピンスケーリング関係
ステップ1:交換はスピンチャネルごとに分離する．第5章・第10章で見たように，交換相互作用は同じスピンを持つ電子の間にしか働かない．KS行列式に対する交換エネルギーは
$$ \begin{equation} E_x=-\frac{1}{2}\sum_{\sigma=\uparrow,\downarrow}\iint \dd^3r\,\dd^3r'\, \frac{\abs{\rho_1^\sigma(\rr,\rr')}^2}{\abs{\rr-\rr'}}, \qquad \rho_1^\sigma(\rr,\rr')=\sum_{i}^{\mathrm{occ},\sigma}\phi_{i\sigma}(\rr)\phi_{i\sigma}^*(\rr') \label{eq:11-ex-spin-split} \end{equation} $$と書ける．$\sigma=\uparrow$ の項は上向き軌道だけ，$\sigma=\downarrow$ の項は下向き軌道だけからできており，両者は完全に独立である．したがって
$$ \begin{equation} E_x[\rho^\uparrow,\rho^\downarrow]=E_x[\rho^\uparrow,0]+E_x[0,\rho^\downarrow]. \label{eq:11-ex-additive} \end{equation} $$ステップ2:片スピンだけの汎関数を無分極汎関数で書く．無分極系($\rho^\uparrow=\rho^\downarrow=n$)に \eqref{eq:11-ex-additive} を適用すると
$$ \begin{equation} E_x[n,n]=E_x[n,0]+E_x[0,n]=2E_x[n,0] \label{eq:11-ex-nn} \end{equation} $$(最後の等号は，上向きだけの系と下向きだけの系が物理的に同じだから)．ここで「全密度 $\rho$ の無分極系の交換エネルギー」を $E_x^{\mathrm{unpol}}[\rho]\equiv E_x[\rho/2,\rho/2]$ と定義すると，$n=\rho/2$ すなわち $\rho=2n$ とおいて $E_x[n,n]=E_x^{\mathrm{unpol}}[2n]$．これを \eqref{eq:11-ex-nn} に入れると
$$ \begin{equation} E_x[n,0]=\tfrac12E_x^{\mathrm{unpol}}[2n]. \label{eq:11-ex-half} \end{equation} $$ステップ3:合成する．\eqref{eq:11-ex-half} を \eqref{eq:11-ex-additive} に代入して
$$ \begin{equation} E_x[\rho^\uparrow,\rho^\downarrow] =\tfrac12E_x^{\mathrm{unpol}}[2\rho^\uparrow]+\tfrac12E_x^{\mathrm{unpol}}[2\rho^\downarrow]. \label{eq:11-spin-scaling} \end{equation} $$これが交換のスピンスケーリング関係である．近似は一切していない——厳密汎関数でもLDAでもGGAでも成り立つ．
ステップ4:LDAに適用する．$E_x^{\mathrm{unpol}}[\rho]=-C_x\int\rho^{4/3}$ を \eqref{eq:11-spin-scaling} に入れる:
\begin{align} E_x^{\mathrm{LSDA}} &=-\frac{C_x}{2}\int \dd^3r\Bigl[(2\rho^\uparrow)^{4/3}+(2\rho^\downarrow)^{4/3}\Bigr] \label{eq:11-lsda-x-1}\\ &=-\frac{C_x}{2}\,2^{4/3}\int \dd^3r\Bigl[(\rho^\uparrow)^{4/3}+(\rho^\downarrow)^{4/3}\Bigr] \label{eq:11-lsda-x-2}\\ &=-\frac{C_x}{2}\,2^{4/3}\left(\frac{1}{2}\right)^{4/3}\int \dd^3r\;\rho^{4/3} \Bigl[(1+\zeta)^{4/3}+(1-\zeta)^{4/3}\Bigr] \label{eq:11-lsda-x-3}\\ &=-\frac{C_x}{2}\int \dd^3r\;\rho^{4/3}\Bigl[(1+\zeta)^{4/3}+(1-\zeta)^{4/3}\Bigr]. \label{eq:11-lsda-x-4} \end{align}\eqref{eq:11-lsda-x-2} で $2^{4/3}$ をくくり出し，\eqref{eq:11-lsda-x-3} では $\rho^\uparrow=\rho(1+\zeta)/2$，$\rho^\downarrow=\rho(1-\zeta)/2$(定義 \eqref{eq:11-zeta} と $\rho=\rho^\uparrow+\rho^\downarrow$ から従う)を代入して $(\rho/2)^{4/3}$ をくくり出した．\eqref{eq:11-lsda-x-4} で $2^{4/3}(1/2)^{4/3}=1$．
検算:$\zeta=0$ で括弧内は $1+1=2$ となり $E_x=-C_x\int\rho^{4/3}$(無分極のLDA)に戻る．$\zeta=1$ で括弧内は $2^{4/3}$ となり $E_x=-2^{1/3}C_x\int\rho^{4/3}$，すなわち完全分極系の交換は無分極系の $2^{1/3}=1.2599$ 倍だけ深い．∎
この結果を「無分極値と完全分極値の内挿」の形に整理しておく．
$$ \begin{equation} \varepsilon_x(\rho,\zeta)=\varepsilon_x(\rho,0) +\bigl[\varepsilon_x(\rho,1)-\varepsilon_x(\rho,0)\bigr]f(\zeta), \qquad f(\zeta)=\frac{(1+\zeta)^{4/3}+(1-\zeta)^{4/3}-2}{2^{4/3}-2} \label{eq:11-fzeta} \end{equation} $$と書けば，$f(0)=0$，$f(1)=1$ であり，$\varepsilon_x(\rho,1)=2^{1/3}\varepsilon_x(\rho,0)$ を使えば
\begin{align} \varepsilon_x(\rho,\zeta) &=\varepsilon_x(\rho,0)\Bigl[1+\bigl(2^{1/3}-1\bigr)f(\zeta)\Bigr] \label{eq:11-fz-check1}\\ &=\varepsilon_x(\rho,0)\left[1+\frac{(2^{1/3}-1)\bigl[(1+\zeta)^{4/3}+(1-\zeta)^{4/3}-2\bigr]}{2(2^{1/3}-1)}\right] \label{eq:11-fz-check2}\\ &=\varepsilon_x(\rho,0)\cdot\frac{(1+\zeta)^{4/3}+(1-\zeta)^{4/3}}{2} \label{eq:11-fz-check3} \end{align}となって \eqref{eq:11-lsda-x-4} と一致する．\eqref{eq:11-fz-check2} では分母の $2^{4/3}-2=2(2^{1/3}-1)$ を使い，\eqref{eq:11-fz-check3} で約分と整理を行った．
記号の注意:第7章では同じ記号 $f(\zeta)$ を，\eqref{eq:11-fz-check3} の右端に現れる交換の増大因子 $\bigl[(1+\zeta)^{4/3}+(1-\zeta)^{4/3}\bigr]/2$($f(0)=1$，$f(1)=2^{1/3}$)の意味で使っている．本章の $f(\zeta)$ は $f(0)=0$，$f(1)=1$ に規格化した内挿関数であり，両者は「第7章の $f$」$=1+(2^{1/3}-1)\times$「本章の $f$」で結ばれる．
注意:相関のスピン内挿は近似である
相関には \eqref{eq:11-spin-scaling} のような厳密なスピンスケーリング関係がない．異なるスピン同士も相関するからである．実際のLSDA汎関数では，QMCで別々に計算された無分極値 $\varepsilon_c(\rho,0)$ と完全分極値 $\varepsilon_c(\rho,1)$ を，交換と同じ内挿関数 $f(\zeta)$ を借用して
$$ \varepsilon_c(\rho,\zeta)\simeq\varepsilon_c(\rho,0)+\bigl[\varepsilon_c(\rho,1)-\varepsilon_c(\rho,0)\bigr]f(\zeta) $$と補間する．これはPZ形の処方であり，$f(\zeta)$ を使う理由は「交換で厳密なのだから相関でもそこそこだろう」という以上のものではない．VWN形はここをもう一歩進め，スピン剛性(スピン磁化率に対応する $\partial^2\varepsilon_c/\partial\zeta^2|_{\zeta=0}$)も正しく再現するように3つの関数を組み合わせる．磁性体の計算で「VWN5 を使え」といった注意が出るのは，この内挿の詳細が磁気モーメントの大きさに効くからである．
この節で得られたもの．一様電子ガスの結果を局所的に貼り付けるだけで，$E_{xc}[\rho]$ の具体形 \eqref{eq:11-lda-def} と，KS方程式に代入できる交換相関ポテンシャル \eqref{eq:11-lda-key} が手に入った．これでDFT計算は原理的に実行可能になる．しかし，これほど乱暴な近似がなぜ機能するのかはまったく明らかでない．次節がその答えである．
11.2 なぜLDAは(不当なほど)うまくいくのか
本節は本章の理論的核心である．まずLDAが本来要求する条件が実在系では全く成り立っていないことを確認し，それでも機能する理由を交換相関ホールという概念で説明する．そのために，(a) 断熱接続(結合定数積分)によって $E_{xc}$ を厳密に書き直し，(b) それがホールとのクーロン相互作用エネルギーであることを示し，(c) ホールが満たす総和則を証明し，(d) $E_{xc}$ がホールの球対称平均にしか依存しないことを導く．この4段構えが揃うと，LDAの成功が「誤差の系統的な相殺」として理解できる．
11.2.1 LDAの適用条件は破れている
LDAの導出(11.1.1節)で暗黙に使ったのは，「微小セルの中で密度が一定とみなせる」という仮定である．これを定量化しよう．セルの大きさとして意味があるのは，電子の量子力学的な「大きさ」，すなわち局所フェルミ波長
$$ \begin{equation} \lambda_F(\rr)=\frac{2\pi}{k_F(\rr)}, \qquad k_F(\rr)=\bigl(3\pi^2\rho(\rr)\bigr)^{1/3} \label{eq:11-lambdaF} \end{equation} $$である．この距離を動く間に密度の相対変化 $\lambda_F\abs{\nabla\rho}/\rho$ が小さいことが要請される:
$$ \begin{equation} \left|\frac{2\pi}{k_F(\rr)}\cdot\frac{\nabla\rho(\rr)}{\rho(\rr)}\right|\ll1. \label{eq:11-lda-condition} \end{equation} $$実在の系でこれが成り立つか確かめる．水素様の裾 $\rho(r)\propto e^{-2\alpha r}$ を持つ原子の外側では $\abs{\nabla\rho}/\rho=2\alpha$ で一定である一方，$k_F\propto\rho^{1/3}\propto e^{-2\alpha r/3}$ は指数的にゼロへ向かう．したがって \eqref{eq:11-lda-condition} の左辺は $r$ とともに指数関数的に発散する．原子核の近くでも密度は激しく変化する．つまりLDAの適用条件は，原子・分子・固体のどこでも——価電子領域ですら——満たされていない．にもかかわらずLDAは，結合長を数%，凝集エネルギーを $0.5$ eV 程度の誤差で再現する．この「理由なき成功」の説明が本節の目標である．
11.2.2 対密度と交換相関ホール
まず，多体波動関数が持つ「電子の避け合い」の情報を，密度に準ずる量として取り出す道具を用意する．
定義:対密度(pair density)
規格化された $N$ 電子状態 $\Psi$ に対して
$$ \begin{equation} n_2(\rr,\rr')\equiv\braket{\Psi\Bigl|\sum_{i\ne j}\delta(\rr-\rr_i)\,\delta(\rr'-\rr_j)\Bigr|\Psi} \label{eq:11-pair-def} \end{equation} $$を対密度と呼ぶ．$n_2(\rr,\rr')\dd^3r\,\dd^3r'$ は，「$\rr$ のまわりの $\dd^3r$ に1個，かつ $\rr'$ のまわりの $\dd^3r'$ に別の1個の電子を見出す確率」に $N(N-1)$ を掛けたものである(スピンについては和をとった)．$i\ne j$ の制限により，同じ電子を二度数えることはない．
対密度には2つの基本性質がある．第一に，電子間相互作用の期待値がこれで書ける:
\begin{align} \braket{\Psi|\hat V_{ee}|\Psi} &=\braket{\Psi\Bigl|\sum_{i\lt j}\frac{1}{\abs{\rr_i-\rr_j}}\Bigr|\Psi} =\frac{1}{2}\braket{\Psi\Bigl|\sum_{i\ne j}\frac{1}{\abs{\rr_i-\rr_j}}\Bigr|\Psi} \label{eq:11-vee-1}\\ &=\frac{1}{2}\iint \dd^3r\,\dd^3r'\; \frac{n_2(\rr,\rr')}{\abs{\rr-\rr'}}. \label{eq:11-vee-2} \end{align}\eqref{eq:11-vee-1} では $i\lt j$ の和を $i\ne j$ の和の半分に書き換えた(各対 $\{i,j\}$ が2回数えられるため)．\eqref{eq:11-vee-2} では，恒等式 $\frac{1}{\abs{\rr_i-\rr_j}}=\iint\dd^3r\dd^3r'\frac{\delta(\rr-\rr_i)\delta(\rr'-\rr_j)}{\abs{\rr-\rr'}}$ を挿入して，デルタ関数の組を \eqref{eq:11-pair-def} の定義にまとめた．
第二に，片方の変数について積分すると密度に戻る:
\begin{align} \int \dd^3r'\,n_2(\rr,\rr') &=\braket{\Psi\Bigl|\sum_{i\ne j}\delta(\rr-\rr_i)\Bigr|\Psi} \label{eq:11-pair-sum-1}\\ &=\braket{\Psi\Bigl|\sum_{i}(N-1)\,\delta(\rr-\rr_i)\Bigr|\Psi} \label{eq:11-pair-sum-2}\\ &=(N-1)\,\rho(\rr). \label{eq:11-pair-sum-3} \end{align}\eqref{eq:11-pair-sum-1} は $\int\dd^3r'\delta(\rr'-\rr_j)=1$ による．\eqref{eq:11-pair-sum-2} では，$i$ を固定したとき $j\ne i$ となる項が $N-1$ 個あることを使った．\eqref{eq:11-pair-sum-3} は密度の定義 $\rho(\rr)=\braket{\Psi|\sum_i\delta(\rr-\rr_i)|\Psi}$ である．
もし電子が互いに完全に独立なら，$\rr$ に1個いることと $\rr'$ に1個いることは無関係で $n_2=\rho(\rr)\rho(\rr')$ となる(厳密には $N/(N-1)$ の因子だけずれるが，それも含めて以下の定義に吸収される)．実際にはパウリ原理とクーロン反発により，$\rr$ に電子がいると $\rr'$ に他の電子が来にくくなる．そのずれを測るのが交換相関ホールである．
定義:交換相関ホール(exchange-correlation hole)
$$ \begin{equation} n_{xc}(\rr,\rr')\equiv\frac{n_2(\rr,\rr')}{\rho(\rr)}-\rho(\rr') \label{eq:11-hole-def} \end{equation} $$すなわち，$\rr$ に電子が1個いるという条件のもとで見た他の電子の密度 $n_2(\rr,\rr')/\rho(\rr)$ から，条件なしの密度 $\rho(\rr')$ を引いたものである．$\rr$ の電子が他の電子を追い払った分だけ $n_{xc}\lt0$ となるので，「電子がまわりに掘った穴(ホール)」と呼ぶ．$\rr$ を基準点，$\rr'$ を観測点と呼ぶことにする．
定理:交換相関ホールの総和則
$$ \begin{equation} \int \dd^3r'\;n_{xc}(\rr,\rr')=-1 \qquad\text{(すべての }\rr\text{ について)} \label{eq:11-sumrule-xc} \end{equation} $$証明．定義 \eqref{eq:11-hole-def} を積分し，\eqref{eq:11-pair-sum-3} と密度の規格化 $\int\rho\,\dd^3r'=N$ を使う:
\begin{align} \int \dd^3r'\,n_{xc}(\rr,\rr') &=\frac{1}{\rho(\rr)}\int \dd^3r'\,n_2(\rr,\rr')-\int \dd^3r'\,\rho(\rr') \label{eq:11-sr-1}\\ &=\frac{(N-1)\rho(\rr)}{\rho(\rr)}-N \label{eq:11-sr-2}\\ &=(N-1)-N=-1. \label{eq:11-sr-3} \end{align} ∎物理的意味:ホールはちょうど電子1個分
総和則 \eqref{eq:11-sumrule-xc} は電荷保存則そのものである．系全体には $N$ 個の電子がいるが，そのうち1個を「いま $\rr$ にいる基準電子」として取り出すと，残りは $N-1$ 個しかない．条件付き密度 $n_2/\rho$ を全空間で積分すると $N-1$，条件なしの $\rho$ を積分すると $N$，その差が $-1$ である．つまりホールは形は系や近似によって違っても，汲み出す電荷は必ずちょうど電子1個分でなければならない．この「全電荷が固定されている」という性質が，後で見るLDAの誤差相殺の鍵になる．
11.2.3 断熱接続(結合定数積分)
ここまでの議論は「ある波動関数 $\Psi$ に対して」の話であった．ところが $E_{xc}$ の定義(第10章)は
$$ \begin{equation} E_{xc}[\rho]=\underbrace{\bigl(T[\rho]-T_s[\rho]\bigr)}_{\text{運動エネルギーの差}} +\underbrace{\bigl(V_{ee}[\rho]-E_{\mathrm{H}}[\rho]\bigr)}_{\text{相互作用の非古典部分}} \label{eq:11-exc-def-recall} \end{equation} $$であり，運動エネルギーの差という「対密度では書けない項」を含んでいる．ここで登場するのが断熱接続である．電子間相互作用の強さを連続的に変えるパラメータ $\lambda$ を導入し，その積分で $E_{xc}$ 全体——運動エネルギー部分も含めて——を相互作用エネルギーだけで表す．これが本章で最も重要な技術的道具である．
定義:断熱接続の族
$0\le\lambda\le1$ をパラメータとして，ハミルトニアンの族
$$ \begin{equation} \hat H_\lambda=\hat T+\lambda\hat V_{ee}+\sum_{i=1}^{N}v_\lambda(\rr_i) \label{eq:11-Hlambda} \end{equation} $$を考える．ここで外部ポテンシャル $v_\lambda(\rr)$ は，$\hat H_\lambda$ の基底状態密度が $\lambda$ によらず与えられた $\rho(\rr)$ に等しくなるように，$\lambda$ ごとに選ぶ．この条件が族全体を貫く拘束である．両端は
- $\lambda=1$:$v_1=v$(実在系の外部ポテンシャル)，$\hat H_1=\hat H$，基底状態 $\Psi_1=\Psi$(真の多体波動関数)
- $\lambda=0$:相互作用のない系．$v_0=v_s=v_{\mathrm{eff}}$(KS有効ポテンシャル)，基底状態 $\Psi_0=\Phi$(KS行列式)
である．$\lambda=0$ の端でちょうどKS系になるのは，「同じ密度を与える非相互作用系」というKS系の定義そのものだからである(第10章)．$\lambda$ を $0$ から $1$ へゆっくり上げる過程で，系はKS系から実在系へ密度を保ったまま連続的に変形される．
数学ノート:Hellmann–Feynman定理
パラメータ $\lambda$ に依存するエルミート演算子 $\hat H_\lambda$ の，規格化された固有状態 $|\Psi_\lambda\rangle$ と固有値 $E_\lambda$ を考える($\hat H_\lambda|\Psi_\lambda\rangle=E_\lambda|\Psi_\lambda\rangle$，$\braket{\Psi_\lambda|\Psi_\lambda}=1$)．このとき
$$ \begin{equation} \frac{\dd E_\lambda}{\dd\lambda}=\braket{\Psi_\lambda\Bigl|\frac{\partial\hat H_\lambda}{\partial\lambda}\Bigr|\Psi_\lambda} \label{eq:11-hf-thm} \end{equation} $$が成り立つ．証明:$E_\lambda=\braket{\Psi_\lambda|\hat H_\lambda|\Psi_\lambda}$ を $\lambda$ で微分すると，積の微分則から3項出る:
\begin{align} \frac{\dd E_\lambda}{\dd\lambda} &=\braket{\partial_\lambda\Psi_\lambda|\hat H_\lambda|\Psi_\lambda} +\braket{\Psi_\lambda|\partial_\lambda\hat H_\lambda|\Psi_\lambda} +\braket{\Psi_\lambda|\hat H_\lambda|\partial_\lambda\Psi_\lambda} \label{eq:11-hf-1}\\ &=E_\lambda\braket{\partial_\lambda\Psi_\lambda|\Psi_\lambda} +E_\lambda\braket{\Psi_\lambda|\partial_\lambda\Psi_\lambda} +\braket{\Psi_\lambda|\partial_\lambda\hat H_\lambda|\Psi_\lambda} \label{eq:11-hf-2}\\ &=E_\lambda\frac{\dd}{\dd\lambda}\braket{\Psi_\lambda|\Psi_\lambda} +\braket{\Psi_\lambda|\partial_\lambda\hat H_\lambda|\Psi_\lambda} \label{eq:11-hf-3}\\ &=\braket{\Psi_\lambda|\partial_\lambda\hat H_\lambda|\Psi_\lambda}. \label{eq:11-hf-4} \end{align}\eqref{eq:11-hf-2} では固有値方程式を2通りに使った．第1項では右側のケットに $\hat H_\lambda|\Psi_\lambda\rangle=E_\lambda|\Psi_\lambda\rangle$ を適用して $\braket{\partial_\lambda\Psi_\lambda|\hat H_\lambda|\Psi_\lambda}=E_\lambda\braket{\partial_\lambda\Psi_\lambda|\Psi_\lambda}$ とした．第3項では左側のブラに作用させる必要があるが，$\hat H_\lambda$ がエルミートで固有値 $E_\lambda$ が実数であることから $\langle\Psi_\lambda|\hat H_\lambda=E_\lambda\langle\Psi_\lambda|$ が成り立ち，$\braket{\Psi_\lambda|\hat H_\lambda|\partial_\lambda\Psi_\lambda}=E_\lambda\braket{\Psi_\lambda|\partial_\lambda\Psi_\lambda}$ となる．\eqref{eq:11-hf-3} では，こうして出た2項が $\lambda$ 微分の積の規則により $E_\lambda\,\dd\braket{\Psi_\lambda|\Psi_\lambda}/\dd\lambda$ にまとまることを使った．\eqref{eq:11-hf-4} では規格化 $\braket{\Psi_\lambda|\Psi_\lambda}=1$ が $\lambda$ によらない定数なので，その微分がゼロであることを使った．
要点は，波動関数の変化による寄与が規格化条件だけで消えることである．だから「ハミルトニアンの $\lambda$ 依存部分の期待値」だけ計算すればよい．第16章のヘルマン・ファインマン力(原子核座標での微分)も同じ定理の応用である．
導出:断熱接続公式 $E_{xc}=\int_0^1\dd\lambda\,\braket{\Psi_\lambda|\hat V_{ee}|\Psi_\lambda}-E_{\mathrm{H}}$
ステップ1:$\lambda$ 微分にHellmann–Feynman定理を適用する．$\hat H_\lambda$ の $\lambda$ 依存部分は，明示的な $\lambda\hat V_{ee}$ と，密度固定の要請から $\lambda$ 依存する外部ポテンシャル $v_\lambda$ の2箇所である:
$$ \begin{equation} \frac{\partial\hat H_\lambda}{\partial\lambda}=\hat V_{ee}+\sum_{i=1}^{N}\frac{\partial v_\lambda(\rr_i)}{\partial\lambda}. \label{eq:11-dH-dlambda} \end{equation} $$これを \eqref{eq:11-hf-thm} に入れる:
\begin{align} \frac{\dd E_\lambda}{\dd\lambda} &=\braket{\Psi_\lambda|\hat V_{ee}|\Psi_\lambda} +\braket{\Psi_\lambda\Bigl|\sum_i\frac{\partial v_\lambda(\rr_i)}{\partial\lambda}\Bigr|\Psi_\lambda} \label{eq:11-ac-1}\\ &=\braket{\Psi_\lambda|\hat V_{ee}|\Psi_\lambda} +\int \dd^3r\,\frac{\partial v_\lambda(\rr)}{\partial\lambda}\,\rho_\lambda(\rr) \label{eq:11-ac-2}\\ &=\braket{\Psi_\lambda|\hat V_{ee}|\Psi_\lambda} +\int \dd^3r\,\frac{\partial v_\lambda(\rr)}{\partial\lambda}\,\rho(\rr). \label{eq:11-ac-3} \end{align}\eqref{eq:11-ac-2} では，1体演算子 $\sum_i u(\rr_i)$ の期待値が $\int u(\rr)\rho_\lambda(\rr)\dd^3r$ で与えられること(密度の定義)を使った．\eqref{eq:11-ac-3} が断熱接続の全体を支える一歩である:族の作り方から $\rho_\lambda(\rr)=\rho(\rr)$ がすべての $\lambda$ で成り立つので，この積分に現れる密度を $\lambda$ に依存しない $\rho(\rr)$ で置き換えられる．密度を固定するという不自然にも見える拘束を課したのは，まさにこの一歩を可能にするためである．
ステップ2:$\lambda$ について $0$ から $1$ まで積分する．
\begin{align} E_1-E_0 &=\int_0^1 \dd\lambda\,\frac{\dd E_\lambda}{\dd\lambda} \label{eq:11-ac-4}\\ &=\int_0^1 \dd\lambda\,\braket{\Psi_\lambda|\hat V_{ee}|\Psi_\lambda} +\int_0^1 \dd\lambda\int \dd^3r\,\frac{\partial v_\lambda(\rr)}{\partial\lambda}\rho(\rr) \label{eq:11-ac-5}\\ &=\int_0^1 \dd\lambda\,\braket{\Psi_\lambda|\hat V_{ee}|\Psi_\lambda} +\int \dd^3r\,\rho(\rr)\int_0^1 \dd\lambda\,\frac{\partial v_\lambda(\rr)}{\partial\lambda} \label{eq:11-ac-6}\\ &=\int_0^1 \dd\lambda\,\braket{\Psi_\lambda|\hat V_{ee}|\Psi_\lambda} +\int \dd^3r\,\rho(\rr)\bigl[v_1(\rr)-v_0(\rr)\bigr]. \label{eq:11-ac-7} \end{align}\eqref{eq:11-ac-4} は微積分学の基本定理．\eqref{eq:11-ac-6} では積分順序を交換した——これが許されるのは，$\rho(\rr)$ が $\lambda$ に依存しないので $\lambda$ 積分の外に出せるからである(ここでも密度固定の拘束が効いている)．\eqref{eq:11-ac-7} で内側の $\lambda$ 積分を実行した:$\int_0^1\partial_\lambda v_\lambda\,\dd\lambda=v_1-v_0$．
ステップ3:両端のエネルギーを書き下す．$\lambda=1$ の端は実在系そのものだから
$$ \begin{equation} E_1=\braket{\Psi|\hat T+\hat V_{ee}|\Psi}+\int \dd^3r\,v(\rr)\rho(\rr) =T+V_{ee}+\int v\rho. \label{eq:11-ac-E1} \end{equation} $$$\lambda=0$ の端は相互作用のないKS系だから，その基底状態はKS行列式 $\Phi$ であり
$$ \begin{equation} E_0=\braket{\Phi|\hat T|\Phi}+\int \dd^3r\,v_s(\rr)\rho(\rr)=T_s+\int v_s\rho. \label{eq:11-ac-E0} \end{equation} $$ステップ4:代入して整理する．$v_1=v$，$v_0=v_s$ に注意して \eqref{eq:11-ac-E1}，\eqref{eq:11-ac-E0} を \eqref{eq:11-ac-7} に入れると
\begin{align} \Bigl(T+V_{ee}+\int v\rho\Bigr)-\Bigl(T_s+\int v_s\rho\Bigr) &=\int_0^1 \dd\lambda\,\braket{\Psi_\lambda|\hat V_{ee}|\Psi_\lambda} +\int \rho\,(v-v_s) \label{eq:11-ac-8}\\ T-T_s+V_{ee} &=\int_0^1 \dd\lambda\,\braket{\Psi_\lambda|\hat V_{ee}|\Psi_\lambda}. \label{eq:11-ac-9} \end{align}\eqref{eq:11-ac-8} から \eqref{eq:11-ac-9} へは，両辺に共通して現れる $\int\rho(v-v_s)$ を消去しただけである．左辺は $E_{xc}$ の定義 \eqref{eq:11-exc-def-recall} とほぼ同じ形をしており，$E_{\mathrm{H}}$ を引けば完成する:
$$ \begin{equation} E_{xc}=(T-T_s)+(V_{ee}-E_{\mathrm{H}}) =\int_0^1 \dd\lambda\,\braket{\Psi_\lambda|\hat V_{ee}|\Psi_\lambda}-E_{\mathrm{H}}. \label{eq:11-ac-final} \end{equation} $$ ∎物理的意味:運動エネルギーはどこへ消えたか
式 \eqref{eq:11-ac-key} の右辺には運動エネルギーが一切現れない．しかし左辺の $E_{xc}$ には $T-T_s$(運動エネルギーの相関寄与)が含まれている．矛盾ではない．$T-T_s$ は$\lambda$ 積分の中に潜んでいるのである．もし $W_\lambda$ が $\lambda$ によらない定数なら $\int_0^1W_\lambda\dd\lambda=W_1=V_{ee}-E_{\mathrm{H}}$ となり，運動エネルギー寄与はゼロになる．実際には $\lambda$ が大きくなるほど電子は強く避け合い，$W_\lambda$ は単調に深くなる(図11.3)．その「$\lambda$ 依存性」が $T-T_s$ を表している．断熱接続の巧妙さは，計算しにくい運動エネルギー項を，計算しやすい(対密度で書ける)相互作用項の $\lambda$ 依存性に変換した点にある．
両端の値には明確な意味がある．$\lambda=0$ ではKS行列式に対する相互作用の期待値からハートリー項を引いたものだから
$$ W_0=\braket{\Phi|\hat V_{ee}|\Phi}-E_{\mathrm{H}}=E_x $$すなわち厳密交換エネルギー(KS軌道で組んだハートリー・フォック型の交換積分，第5章)である．これは11.6節でハイブリッド汎関数を正当化するときの出発点になる．$\lambda=1$ では $W_1=V_{ee}-E_{\mathrm{H}}$ で，真の波動関数における非古典的相互作用エネルギーである．
11.2.4 結合定数平均されたホールによる $E_{xc}$ の表現
断熱接続公式 \eqref{eq:11-ac-key} の右辺は相互作用の期待値だけでできているから，11.2.2節の対密度・ホールの言葉に翻訳できる．
各 $\lambda$ の状態 $\Psi_\lambda$ に対して対密度 $n_2^\lambda(\rr,\rr')$ とホール $n_{xc}^\lambda(\rr,\rr')$ を \eqref{eq:11-pair-def}，\eqref{eq:11-hole-def} で定義する．定義を移項すると $n_2^\lambda=\rho(\rr)\bigl[\rho(\rr')+n_{xc}^\lambda(\rr,\rr')\bigr]$ であり(密度は $\lambda$ によらないので $\rho$ に添字は不要)，これを \eqref{eq:11-vee-2} に入れて
\begin{align} \braket{\Psi_\lambda|\hat V_{ee}|\Psi_\lambda} &=\frac{1}{2}\iint \dd^3r\,\dd^3r'\,\frac{\rho(\rr)\bigl[\rho(\rr')+n_{xc}^\lambda(\rr,\rr')\bigr]}{\abs{\rr-\rr'}} \label{eq:11-holeE-1}\\ &=\underbrace{\frac{1}{2}\iint \dd^3r\,\dd^3r'\,\frac{\rho(\rr)\rho(\rr')}{\abs{\rr-\rr'}}}_{=\,E_{\mathrm{H}}[\rho]} +\frac{1}{2}\iint \dd^3r\,\dd^3r'\,\frac{\rho(\rr)\,n_{xc}^\lambda(\rr,\rr')}{\abs{\rr-\rr'}} \label{eq:11-holeE-2} \end{align}となる．\eqref{eq:11-holeE-2} では単に和を2つに分け，第1項がハートリーエネルギーの定義そのものであることを認めた．したがって
$$ \begin{equation} W_\lambda=\braket{\Psi_\lambda|\hat V_{ee}|\Psi_\lambda}-E_{\mathrm{H}} =\frac{1}{2}\iint \dd^3r\,\dd^3r'\,\frac{\rho(\rr)\,n_{xc}^\lambda(\rr,\rr')}{\abs{\rr-\rr'}}. \label{eq:11-Wlambda-hole} \end{equation} $$これを \eqref{eq:11-ac-key} に代入し，$\lambda$ 積分を内側に押し込む($\rho$ もクーロン核も $\lambda$ に依存しないので，$\lambda$ 積分は $n_{xc}^\lambda$ にだけ掛かる):
\begin{align} E_{xc} &=\int_0^1 \dd\lambda\;\frac{1}{2}\iint \dd^3r\,\dd^3r'\,\frac{\rho(\rr)\,n_{xc}^\lambda(\rr,\rr')}{\abs{\rr-\rr'}} \label{eq:11-exc-hole-1}\\ &=\frac{1}{2}\iint \dd^3r\,\dd^3r'\,\frac{\rho(\rr)}{\abs{\rr-\rr'}} \underbrace{\int_0^1 \dd\lambda\;n_{xc}^\lambda(\rr,\rr')}_{\displaystyle\equiv\,\bar n_{xc}(\rr,\rr')}. \label{eq:11-exc-hole-2} \end{align}定義:結合定数平均された交換相関ホール
$$ \begin{equation} \bar n_{xc}(\rr,\rr')\equiv\int_0^1 \dd\lambda\;n_{xc}^\lambda(\rr,\rr') \label{eq:11-holebar-def} \end{equation} $$これは驚くほど古典的な式である．交換相関エネルギーは，電子密度 $\rho(\rr)$ と，各電子が自分のまわりに引き連れている「電荷 $-1$ の穴」 $\bar n_{xc}(\rr,\rr')$ との，静電相互作用エネルギーにほかならない．量子力学的な多体効果のすべてが，この「穴の形」に押し込められている．
11.2.5 交換ホールと相関ホール，その総和則
ホールを交換部分と相関部分に分けておく．$\lambda=0$ のホール，すなわちKS行列式に対するホールを交換ホール $n_x$ と呼び，残りを相関ホール $n_c$ と呼ぶ:
$$ \begin{equation} n_x(\rr,\rr')\equiv n_{xc}^{\lambda=0}(\rr,\rr'), \qquad n_c(\rr,\rr')\equiv\bar n_{xc}(\rr,\rr')-n_x(\rr,\rr'). \label{eq:11-nx-nc-def} \end{equation} $$この分け方は $E_{xc}=E_x+E_c$ の分割と整合している．実際 \eqref{eq:11-exc-hole-key} に代入すれば
$$ \begin{equation} E_x=\frac{1}{2}\iint \dd^3r\,\dd^3r'\,\frac{\rho(\rr)n_x(\rr,\rr')}{\abs{\rr-\rr'}}, \qquad E_c=\frac{1}{2}\iint \dd^3r\,\dd^3r'\,\frac{\rho(\rr)n_c(\rr,\rr')}{\abs{\rr-\rr'}} \label{eq:11-ex-ec-hole} \end{equation} $$であり，$E_x=W_0$ という前節の結論と一致する．
導出:交換ホールの具体形と総和則 $\int n_x\dd^3r'=-1$
$\lambda=0$ の状態はKS行列式 $\Phi$ である．第5章・第6章で示したように，単一スレーター行列式の対密度は1体密度行列だけで書ける:
$$ \begin{equation} n_2^{\mathrm{KS}}(\rr,\rr')=\rho(\rr)\rho(\rr')-\sum_{\sigma=\uparrow,\downarrow}\abs{\rho_1^\sigma(\rr,\rr')}^2, \qquad \rho_1^\sigma(\rr,\rr')=\sum_{i}^{\mathrm{occ},\sigma}\phi_{i\sigma}(\rr)\phi_{i\sigma}^*(\rr'). \label{eq:11-n2-ks} \end{equation} $$第1項は「独立粒子」的な寄与，第2項がパウリ原理による交換の補正である(異なるスピン同士は行列式でも無相関なので，和はスピンチャネルごとに閉じている)．ホールの定義 \eqref{eq:11-hole-def} に入れると
$$ \begin{equation} n_x(\rr,\rr')=\frac{n_2^{\mathrm{KS}}(\rr,\rr')}{\rho(\rr)}-\rho(\rr') =-\frac{1}{\rho(\rr)}\sum_{\sigma}\abs{\rho_1^\sigma(\rr,\rr')}^2. \label{eq:11-nx-explicit} \end{equation} $$まず読み取れるのは$n_x(\rr,\rr')\le0$ がすべての $\rr,\rr'$ で成り立つことである(絶対値の2乗の和にマイナスが付いている)．交換ホールは決して正にならない．この符号条件から直ちに $E_x\le0$ も従う(\eqref{eq:11-ex-ec-hole} の被積分関数が至るところ $\le0$ だから)．
総和則を確かめる．$\rr'$ について積分すると
\begin{align} \int \dd^3r'\,\abs{\rho_1^\sigma(\rr,\rr')}^2 &=\int \dd^3r'\sum_{i}^{\mathrm{occ},\sigma}\sum_{j}^{\mathrm{occ},\sigma} \phi_{i\sigma}(\rr)\phi_{i\sigma}^*(\rr')\phi_{j\sigma}^*(\rr)\phi_{j\sigma}(\rr') \label{eq:11-nxsum-1}\\ &=\sum_{i,j}^{\mathrm{occ},\sigma}\phi_{i\sigma}(\rr)\phi_{j\sigma}^*(\rr) \underbrace{\int \dd^3r'\,\phi_{i\sigma}^*(\rr')\phi_{j\sigma}(\rr')}_{=\,\delta_{ij}} \label{eq:11-nxsum-2}\\ &=\sum_{i}^{\mathrm{occ},\sigma}\abs{\phi_{i\sigma}(\rr)}^2=\rho^\sigma(\rr). \label{eq:11-nxsum-3} \end{align}\eqref{eq:11-nxsum-1} は $\abs{\rho_1^\sigma}^2=\rho_1^\sigma(\rho_1^\sigma)^*$ を展開しただけ．\eqref{eq:11-nxsum-2} で $\rr'$ に依存する因子だけを積分の中に残し，KS軌道の規格直交性 $\int\phi_{i\sigma}^*\phi_{j\sigma}\dd^3r'=\delta_{ij}$ を使った．\eqref{eq:11-nxsum-3} でクロネッカーデルタにより和をつぶした．したがって
$$ \begin{equation} \int \dd^3r'\,n_x(\rr,\rr')=-\frac{1}{\rho(\rr)}\sum_\sigma\rho^\sigma(\rr) =-\frac{\rho(\rr)}{\rho(\rr)}=-1. \label{eq:11-sumrule-x} \end{equation} $$ ∎注意:第5章・第9章のスピン分解した交換ホールとの対応
第5章5.8.4項(および第9章9.1.2節)では，交換ホールをスピンチャネルごとに $n_{\mathrm{x}}^\sigma(\rr,\rr')=-\abs{\rho_1^\sigma(\rr,\rr')}^2/\rho^\sigma(\rr)$ と定義し， $\int n_{\mathrm{x}}^\sigma\dd^3r'=-1$ と $E_x=\frac12\sum_\sigma\iint\dd^3r\,\dd^3r'\,\rho^\sigma(\rr)\,n_{\mathrm{x}}^\sigma(\rr,\rr')/\abs{\rr-\rr'}$ を示した． 本章の $n_x$ \eqref{eq:11-nx-explicit} は，これをスピンで重み付き平均したもの
$$ n_x(\rr,\rr')=\sum_\sigma\frac{\rho^\sigma(\rr)}{\rho(\rr)}\,n_{\mathrm{x}}^\sigma(\rr,\rr') $$である．重み $\rho^\sigma/\rho$ の和が $1$ なので総和則は $\int n_x\dd^3r'=\sum_\sigma(\rho^\sigma/\rho)(-1)=-1$ となって両者で一致し， $\rho\,n_x=\sum_\sigma\rho^\sigma n_{\mathrm{x}}^\sigma$ より $E_x$ の表式 \eqref{eq:11-ex-ec-hole} も第5章の結果と同一である． 記号が違うだけで，定義・符号・総和則のいずれにも食い違いはない． なお第5章では密度を $n$，密度行列を $\gamma^\sigma$ と書いたが，本章の $\rho$，$\rho_1^\sigma$ と同じ量である ($\gamma^\sigma=(\rho_1^\sigma)^*$ で，使うのは $\abs{\cdot}^2$ だけなのでどちらでも結果は変わらない)．
相関ホールの総和則はこれらの引き算で出る．まず，総和則 \eqref{eq:11-sumrule-xc} の証明は「$\Psi_\lambda$ が規格化されていること」と「密度が $\rho$ であること」しか使っていないので，すべての $\lambda$ について $\int n_{xc}^\lambda\dd^3r'=-1$ が成り立つ．したがって
\begin{align} \int \dd^3r'\,\bar n_{xc}(\rr,\rr') &=\int_0^1 \dd\lambda\int \dd^3r'\,n_{xc}^\lambda(\rr,\rr') =\int_0^1 \dd\lambda\,(-1)=-1, \label{eq:11-sumrule-bar}\\ \int \dd^3r'\,n_c(\rr,\rr') &=\int \dd^3r'\,\bar n_{xc}-\int \dd^3r'\,n_x=(-1)-(-1)=0. \label{eq:11-sumrule-c} \end{align}物理的意味:2つの総和則の対比
交換ホールは電荷を汲み出す．$\int n_x=-1$ は，パウリ原理により同スピン電子が完全に排除されること(電子1個分)を表す．したがって $E_x$ は必ず負で，しかも大きい．
相関ホールは電荷を運ぶだけである．$\int n_c=0$ は，クーロン相関が電子を「近くから遠くへ押しやる」だけで，正味の電荷を汲み出さないことを意味する．近距離では $n_c\lt0$(さらに避ける)，遠距離では $n_c\gt0$(押しやられた分が積み上がる)となり，両者が打ち消し合う．$E_c$ が $E_x$ に比べて1桁小さい(表11.1)のは，この打ち消しのためである．
この2つの制約は，近似汎関数を作るときの最も基本的な指針になる．「モデルホールを作り，総和則を課し，\eqref{eq:11-exc-hole-key} で積分する」というのが11.4節で見るGGA構成法の骨格である．
11.2.6 $E_{xc}$ はホールの球対称平均にしか依存しない
ここが本節の山場である．式 \eqref{eq:11-exc-hole-key} の二重積分を，基準点 $\rr$ からの相対座標に書き換えて角度積分を実行すると，ホールの詳細な形が驚くほど落ちることが分かる．
導出:球対称平均への帰着
ステップ1:相対座標に移る．$\bm{s}\equiv\rr'-\rr$ と置く．$\rr$ を固定して $\rr'$ を積分するとき，$\rr'=\rr+\bm{s}$ で $\dd^3r'=\dd^3s$(平行移動のヤコビアンは1)．また $\abs{\rr-\rr'}=\abs{\bm{s}}\equiv s$．したがって
$$ \begin{equation} E_{xc}=\frac{1}{2}\int \dd^3r\;\rho(\rr)\int \dd^3s\;\frac{\bar n_{xc}(\rr,\rr+\bm{s})}{s}. \label{eq:11-sph-1} \end{equation} $$ステップ2:球座標に分解する．$\bm{s}$ の積分を球座標で書く．$\dd^3s=s^2\,\dd s\,\dd\Omega_s$($\dd\Omega_s=\sin\theta\,\dd\theta\,\dd\varphi$ は $\bm{s}$ の向きの立体角要素)．すると
\begin{align} \int \dd^3s\;\frac{\bar n_{xc}(\rr,\rr+\bm{s})}{s} &=\int_0^\infty \dd s\;s^2\int \dd\Omega_s\;\frac{\bar n_{xc}(\rr,\rr+\bm{s})}{s} \label{eq:11-sph-2}\\ &=\int_0^\infty \dd s\;s\int \dd\Omega_s\;\bar n_{xc}(\rr,\rr+\bm{s}). \label{eq:11-sph-3} \end{align}\eqref{eq:11-sph-2} から \eqref{eq:11-sph-3} へは $s^2/s=s$ と約分しただけである．ここが決定的な点である:クーロン核 $1/s$ は $\bm{s}$ の大きさにしか依存せず向きには依存しないので，立体角積分の外に出せてしまう．その結果，角度積分は $\bar n_{xc}$ だけに掛かる．
ステップ3:球対称平均を定義する．
$$ \begin{equation} \bar n_{xc}^{\mathrm{SA}}(\rr,s)\equiv\frac{1}{4\pi}\int \dd\Omega_s\;\bar n_{xc}(\rr,\rr+\bm{s}) \label{eq:11-SA-def} \end{equation} $$と定義する(SA = spherical average)．これは「基準点 $\rr$ から距離 $s$ だけ離れた球面上でホールを平均した値」である．$\int\dd\Omega_s\,\bar n_{xc}=4\pi\,\bar n_{xc}^{\mathrm{SA}}$ だから \eqref{eq:11-sph-3} は
$$ \begin{equation} \int \dd^3s\;\frac{\bar n_{xc}}{s}=\int_0^\infty \dd s\;4\pi s\;\bar n_{xc}^{\mathrm{SA}}(\rr,s). \label{eq:11-sph-4} \end{equation} $$ステップ4:結果．\eqref{eq:11-sph-1} に戻して
$$ \begin{equation} E_{xc}=\frac{1}{2}\int \dd^3r\;\rho(\rr)\int_0^\infty 4\pi s\,\dd s\;\bar n_{xc}^{\mathrm{SA}}(\rr,s). \label{eq:11-sph-final} \end{equation} $$ ∎右の等式は総和則 \eqref{eq:11-sumrule-bar} を同じ球座標分解で書き直したものである($\int\dd^3s\,\bar n_{xc}=\int_0^\infty s^2\dd s\int\dd\Omega_s\,\bar n_{xc}=\int_0^\infty4\pi s^2\dd s\,\bar n_{xc}^{\mathrm{SA}}$)．
物理的意味:近似汎関数に課される要求は驚くほど緩い
式 \eqref{eq:11-sph-key} は次の3点を主張している．
- ホールの角度依存性は $E_{xc}$ に一切効かない．クーロン相互作用が等方的($1/s$)であるため，角度積分がホールを平均してしまう．ホールが基準電子のまわりでどれほど歪んでいようと，球面平均が同じなら $E_{xc}$ は同じである．
- 重みは $s$ であって $s^2$ ではない．体積要素の $s^2$ の一つがクーロン核の $1/s$ で相殺された結果，\eqref{eq:11-sph-key} 左の積分では重み $4\pi s$ が掛かる．一方，総和則(右)の重みは $4\pi s^2$ である．つまり $E_{xc}$ は総和則よりも近距離を重く見る．
- 総和則が全体の目盛りを固定する．球平均ホールの「総量」が $-1$ に釘付けされているので，形が多少違っても $E_{xc}$ の値は大きくは外れない．
近似汎関数に要求されるのは「厳密なホールを再現すること」ではなく，「球平均したホールを，近距離を重視した重みのもとでそれらしく再現し，かつ総和則を満たすこと」だけである．これがLDAの成功を理解する鍵になる．
11.2.7 LDAのホールはなぜ十分によいのか
LDAの $E_{xc}$ は，定義からして「基準点 $\rr$ における密度 $\rho(\rr)$ の一様電子ガスのホール」
$$ \begin{equation} n_{xc}^{\mathrm{LDA}}(\rr,\rr+\bm{s})=\bar n_{xc}^{\mathrm{unif}}\bigl(\rho(\rr);\,s\bigr) \label{eq:11-lda-hole} \end{equation} $$を使って \eqref{eq:11-exc-hole-key} を評価することに等しい．実際，一様電子ガスに対する \eqref{eq:11-sph-final} は $\rho$ が定数だから $E_{xc}^{\mathrm{unif}}=\rho V\cdot\frac12\int_0^\infty4\pi s\,\dd s\,\bar n_{xc}^{\mathrm{unif}}(\rho;s)$，すなわち電子1個あたり $\varepsilon_{xc}(\rho)=\frac12\int_0^\infty4\pi s\,\dd s\,\bar n_{xc}^{\mathrm{unif}}(\rho;s)$ となり，これを $\rho\to\rho(\rr)$ として \eqref{eq:11-sph-final} に代入すれば $E_{xc}^{\mathrm{LDA}}=\int\rho\varepsilon_{xc}(\rho)$ に戻る．
このLDAホールには2つの重大な長所がある．
- 総和則を自動的に満たす．一様電子ガスの厳密なホールは総和則 \eqref{eq:11-sumrules-key} を満たしている．LDAはそれを密度 $\rho(\rr)$ で評価して使うだけなので，$\int n_{xc}^{\mathrm{LDA}}\dd^3s=-1$ が任意の $\rr$ で成り立つ．「電子1個分の穴」という最も基本的な物理が保たれている．
- 符号条件と近距離の振る舞いが正しい．一様ガスのホールは $n_x\le0$ を満たす．また $s\to0$ での深さが物理的に正しい:スピン無分極では交換ホールは同スピン電子(全体の半分)しか排除しないので $n_x(s\to0)=-\rho/2$ となり(演習11.2)，これに相関ホールの寄与が加わってさらに深くなる．$s$ が大きくなるとゼロに近づく．$E_{xc}$ を強く支配する近距離領域(重み $4\pi s$)の振る舞いが妥当であることが重要である．
他方，LDAホールは球対称であり，しかも基準電子を中心に据えている．実在系の厳密なホールはそうではない．たとえば原子の中で基準電子が核から少し離れた位置にあると，他の電子は核に強く引き付けられているため，ホール(=他の電子が来られない領域)は基準電子を中心とせず核の方向へ引き伸ばされる．形としてはまったく別物である．
ところが11.2.6節で見たとおり，$E_{xc}$ が見るのは球平均だけである．歪んだ厳密ホールを基準電子まわりで球面平均すると，球対称なLDAホールと非常に近い曲線になる．しかも両者の面積(総和則)は $-1$ で一致している．これが「形は違うがエネルギーは合う」という誤差相殺の内実である．
物理的意味:LDAの成功を支える2つの機構
(1) ホールの球平均と総和則(本節)．$E_{xc}$ が球平均にしか依存せず，その総量が $-1$ に固定されているため，LDAホールの「形の間違い」は $E_{xc}$ にほとんど伝わらない．適用条件 \eqref{eq:11-lda-condition} が破れていても機能するのはこのためである．逆に言えば，総和則を破る近似は，条件 \eqref{eq:11-lda-condition} に近い状況でさえ悪い結果を与えうる．11.4節で見る勾配展開(GEA)の失敗がまさにそれである．
(2) エネルギー差における系統誤差の相殺．実際に物理量として意味を持つのは全エネルギーそのものではなく，その差である:結合エネルギー(分子と原子)，凝集エネルギー(固体と孤立原子)，構造相転移のエネルギー差，反応障壁．LDAの誤差は密度の絶対値に依存する系統的なものなので，似た密度環境どうしを比べる差では大きく打ち消される．第10章のΔSCF法が実用的に働くのも同じ理屈である．
この2つはしばしば混同されるが，独立な機構である．(1) は $E_{xc}$ 自体の誤差が小さいことを，(2) は残った誤差が差で消えることを説明している．
注意:LDAが破綻する場面は総和則からは救えない
上の議論は「一様電子ガスのホールが，実在系の球平均ホールの良い近似である」という前提に立っている．この前提が成り立たない状況ではLDAは実際に破綻する．代表的なものは:
- 密度が重なっていない2つのフラグメントの間(van der Waals力)．基準電子から遠く離れた別の分子に相関ホールの尾が伸びるべきなのに，局所汎関数はそれを表現できない(11.6節)．
- 1電子系・少数電子系．ホールは厳密には $n_{xc}=-\rho(\rr')$(自分自身の密度そのもの)でなければならないが，一様ガスのホールはそうならない．これが自己相互作用誤差である(11.3.5節)．
- 強く局在した $d$・$f$ 電子．一様ガスは電子が遍歴する系なので，局在極限のホールを表現できない(11.6節のDFT+U)．
この節で得られたもの．断熱接続 \eqref{eq:11-ac-key} により，$E_{xc}$ は「密度と交換相関ホールの静電相互作用」\eqref{eq:11-exc-hole-key} という古典的な形に厳密に書き直された．さらに，クーロン核の等方性から $E_{xc}$ はホールの球平均 \eqref{eq:11-sph-key} にしか依存せず，その総量は総和則 \eqref{eq:11-sumrules-key} で $-1$ に固定されている．LDAはホールの形を大きく誤るが，球平均と総和則は満たすので，$E_{xc}$ の誤差は小さく抑えられる．次節では，厳密汎関数が満たすべき制約をさらに列挙し，汎関数設計の「試験問題」を用意する．
11.3 厳密な制約
厳密な $E_{xc}[\rho]$ は誰も書き下せないが，それが満たさねばならない性質はいくつも証明できる．これらは近似汎関数に課す「試験問題」として使える:制約を満たすように関数形とパラメータを決めれば，経験的なフィットに頼らずに汎関数が作れる．本節ではその主要なものを導出する．
11.3.1 一様スケーリングと交換のスケーリング等式
定義:密度の一様スケーリング
正の数 $\gamma$ に対して
$$ \begin{equation} \rho_\gamma(\rr)\equiv\gamma^3\rho(\gamma\rr) \label{eq:11-scal-def} \end{equation} $$を $\rho$ の一様スケール密度と呼ぶ．$\gamma\gt1$ は密度を原点に向かって圧縮(高密度化)，$\gamma\lt1$ は膨張(低密度化)する操作である．前の因子 $\gamma^3$ は電子数を保つためのもので，実際
$$ \begin{equation} \int \dd^3r\,\rho_\gamma(\rr)=\gamma^3\int \dd^3r\,\rho(\gamma\rr) =\gamma^3\cdot\frac{1}{\gamma^3}\int \dd^3u\,\rho(\bm{u})=N \label{eq:11-scal-norm} \end{equation} $$である(変数変換 $\bm{u}=\gamma\rr$，$\dd^3r=\dd^3u/\gamma^3$ を行った)．また $(\rho_\gamma)_{\gamma'}=\rho_{\gamma\gamma'}$ が成り立つ:$\bigl[(\rho_\gamma)_{\gamma'}\bigr](\rr)=\gamma'^3\rho_\gamma(\gamma'\rr)=\gamma'^3\gamma^3\rho(\gamma\gamma'\rr)=\rho_{\gamma\gamma'}(\rr)$．
定理:交換エネルギーのスケーリング等式
$$ \begin{equation} E_x[\rho_\gamma]=\gamma\,E_x[\rho] \label{eq:11-ex-scaling} \end{equation} $$導出:$E_x[\rho_\gamma]=\gamma E_x[\rho]$
厳密な交換エネルギーは11.2.3節で見たとおり $E_x[\rho]=\braket{\Phi[\rho]|\hat V_{ee}|\Phi[\rho]}-E_{\mathrm{H}}[\rho]$ である($\Phi[\rho]$ は密度 $\rho$ を与えるKS行列式)．3つの部品に分けて示す．
ステップ1:スケールした行列式は，スケールした密度のKS行列式である．$\Phi$ を作る軌道 $\{\phi_i\}$ から
$$ \begin{equation} \phi_i^\gamma(\rr)\equiv\gamma^{3/2}\phi_i(\gamma\rr) \label{eq:11-scal-orb} \end{equation} $$を作る．(i) 規格直交性:$\int\phi_i^{\gamma*}\phi_j^\gamma\dd^3r=\gamma^3\int\phi_i^*(\gamma\rr)\phi_j(\gamma\rr)\dd^3r=\int\phi_i^*(\bm u)\phi_j(\bm u)\dd^3u=\delta_{ij}$(同じ変数変換)．(ii) 密度:$\sum_i\abs{\phi_i^\gamma(\rr)}^2=\gamma^3\sum_i\abs{\phi_i(\gamma\rr)}^2=\gamma^3\rho(\gamma\rr)=\rho_\gamma(\rr)$．(iii) 運動エネルギー:連鎖律により $\nabla^2\phi_i^\gamma(\rr)=\gamma^{3/2}\gamma^2(\nabla^2\phi_i)(\gamma\rr)$ だから $T_s[\{\phi_i^\gamma\}]=\gamma^2T_s[\{\phi_i\}]$．因子 $\gamma^2$ は全軌道に共通の正の数なので，$T_s$ を最小にする軌道の組はスケーリングで移り合う．ゆえに $\Phi[\rho_\gamma]=\Phi[\rho]_\gamma$，すなわちスケールした行列式が $\rho_\gamma$ のKS行列式である．
ステップ2:相互作用の期待値は $\gamma$ 倍になる．$\Phi_\gamma(\rr_1,\dots,\rr_N)=\gamma^{3N/2}\Phi(\gamma\rr_1,\dots,\gamma\rr_N)$ に対し
\begin{align} \braket{\Phi_\gamma|\hat V_{ee}|\Phi_\gamma} &=\int\!\cdots\!\int \dd^3r_1\cdots \dd^3r_N\; \gamma^{3N}\abs{\Phi(\gamma\rr_1,\dots,\gamma\rr_N)}^2\sum_{i\lt j}\frac{1}{\abs{\rr_i-\rr_j}} \label{eq:11-vee-scal-1}\\ &=\int\!\cdots\!\int \frac{\dd^3u_1\cdots \dd^3u_N}{\gamma^{3N}}\; \gamma^{3N}\abs{\Phi(\bm u_1,\dots,\bm u_N)}^2\sum_{i\lt j}\frac{\gamma}{\abs{\bm u_i-\bm u_j}} \label{eq:11-vee-scal-2}\\ &=\gamma\braket{\Phi|\hat V_{ee}|\Phi}. \label{eq:11-vee-scal-3} \end{align}\eqref{eq:11-vee-scal-2} では全電子座標を $\bm u_i=\gamma\rr_i$ と変換した．体積要素から $\gamma^{-3N}$，クーロン核 $1/\abs{\rr_i-\rr_j}=\gamma/\abs{\bm u_i-\bm u_j}$ から $\gamma$ が出る．\eqref{eq:11-vee-scal-3} で $\gamma^{3N}\cdot\gamma^{-3N}=1$ と約分し，$\gamma$ だけが残った．
ステップ3:ハートリー項も $\gamma$ 倍になる．
\begin{align} E_{\mathrm{H}}[\rho_\gamma] &=\frac{1}{2}\iint \dd^3r\,\dd^3r'\,\frac{\gamma^3\rho(\gamma\rr)\,\gamma^3\rho(\gamma\rr')}{\abs{\rr-\rr'}} \label{eq:11-eh-scal-1}\\ &=\frac{1}{2}\iint \frac{\dd^3u}{\gamma^3}\frac{\dd^3u'}{\gamma^3}\, \frac{\gamma^6\rho(\bm u)\rho(\bm u')\,\gamma}{\abs{\bm u-\bm u'}} \label{eq:11-eh-scal-2}\\ &=\gamma E_{\mathrm{H}}[\rho]. \label{eq:11-eh-scal-3} \end{align}同じ変数変換である:$\gamma^6$(密度から)$\times\gamma^{-6}$(2つの体積要素から)$\times\gamma$(クーロン核から)$=\gamma$．
ステップ4:合成する．
$$ \begin{equation} E_x[\rho_\gamma]=\braket{\Phi_\gamma|\hat V_{ee}|\Phi_\gamma}-E_{\mathrm{H}}[\rho_\gamma] =\gamma\braket{\Phi|\hat V_{ee}|\Phi}-\gamma E_{\mathrm{H}}[\rho]=\gamma E_x[\rho]. \label{eq:11-ex-scal-final} \end{equation} $$ ∎例:LDA交換とGGA交換はスケーリング等式を自動的に満たす
LDA:式 \eqref{eq:11-ex-lda} に $\rho_\gamma$ を入れる．
\begin{align} E_x^{\mathrm{LDA}}[\rho_\gamma] &=-C_x\int \dd^3r\,\bigl[\gamma^3\rho(\gamma\rr)\bigr]^{4/3} =-C_x\gamma^4\int \dd^3r\,\rho(\gamma\rr)^{4/3} \label{eq:11-ldascal-1}\\ &=-C_x\gamma^4\cdot\frac{1}{\gamma^3}\int \dd^3u\,\rho(\bm u)^{4/3} =\gamma E_x^{\mathrm{LDA}}[\rho]. \label{eq:11-ldascal-2} \end{align}GGA:11.4節で導入する形 $E_x^{\mathrm{GGA}}=-C_x\int\rho^{4/3}F_x(s)\dd^3r$ を考える．無次元換算勾配 $s=\abs{\nabla\rho}/\bigl(2(3\pi^2)^{1/3}\rho^{4/3}\bigr)$ は $\abs{\nabla\rho}/\rho^{4/3}$ に比例する．$\rho_\gamma$ に対しては，連鎖律から $\nabla\rho_\gamma(\rr)=\gamma^4(\nabla\rho)(\gamma\rr)$，また $\rho_\gamma^{4/3}(\rr)=\gamma^4\rho(\gamma\rr)^{4/3}$ なので，$\gamma^4$ が分子分母で相殺して
$$ s_\gamma(\rr)=s(\gamma\rr) $$となる．すなわち $s$ はスケーリングで値を変えずに座標だけが移される．よって \eqref{eq:11-ldascal-2} と同じ変数変換で $E_x^{\mathrm{GGA}}[\rho_\gamma]=\gamma E_x^{\mathrm{GGA}}[\rho]$ が従う．「$\rho^{4/3}$ に $s$ の関数を掛ける」という形を選んだ時点で，交換のスケーリング等式は自動的に満たされる——これがGGAの標準形が広く採用されている理由の一つである．
11.3.2 相関のスケーリング不等式
相関エネルギーには等式ではなく不等式が課される．証明には第9章のLevyの制約付き探索を使う．
定理:相関エネルギーのスケーリング不等式
$$ \begin{equation} E_c[\rho_\gamma]\le\gamma E_c[\rho]\quad(\gamma\le1), \qquad E_c[\rho_\gamma]\ge\gamma E_c[\rho]\quad(\gamma\ge1) \label{eq:11-ec-scaling} \end{equation} $$導出:相関のスケーリング不等式
ステップ1:$E_c$ を2つの制約付き最小化の差として書く．Levyの流儀(第9章)で
$$ \begin{equation} F[\rho]=\min_{\Psi\to\rho}\braket{\Psi|\hat T+\hat V_{ee}|\Psi}=\braket{\Psi[\rho]|\hat T+\hat V_{ee}|\Psi[\rho]}, \qquad T_s[\rho]=\min_{\Phi\to\rho}\braket{\Phi|\hat T|\Phi} \label{eq:11-levy-recall} \end{equation} $$と定義する($\Psi$ は $\rho$ を与える任意の反対称波動関数，$\Phi$ は $\rho$ を与える単一行列式)．$E_{xc}=F-T_s-E_{\mathrm{H}}$，$E_x=\braket{\Phi[\rho]|\hat V_{ee}|\Phi[\rho]}-E_{\mathrm{H}}$ だから
\begin{align} E_c[\rho]&=E_{xc}[\rho]-E_x[\rho] \notag\\ &=\braket{\Psi[\rho]|\hat T+\hat V_{ee}|\Psi[\rho]}-\braket{\Phi[\rho]|\hat T|\Phi[\rho]}-E_{\mathrm{H}} -\braket{\Phi[\rho]|\hat V_{ee}|\Phi[\rho]}+E_{\mathrm{H}} \notag\\ &=\braket{\Psi[\rho]|\hat T+\hat V_{ee}|\Psi[\rho]}-\braket{\Phi[\rho]|\hat T+\hat V_{ee}|\Phi[\rho]}. \label{eq:11-ec-as-diff} \end{align}$\Phi[\rho]$ も密度 $\rho$ を与える1つの波動関数だから，$\Psi[\rho]$ の最小性により $E_c[\rho]\le0$ が直ちに従う．以下では
$$ \begin{equation} T_c[\rho]\equiv\braket{\Psi[\rho]|\hat T|\Psi[\rho]}-T_s[\rho]\;(\ge0), \qquad U_c[\rho]\equiv\braket{\Psi[\rho]|\hat V_{ee}|\Psi[\rho]}-\braket{\Phi[\rho]|\hat V_{ee}|\Phi[\rho]} \label{eq:11-Tc-Uc} \end{equation} $$と置く($E_c=T_c+U_c$)．$T_c\ge0$ は $T_s$ が $\rho$ を与える行列式の中で運動エネルギーを最小化する量であり，$\Psi[\rho]$ から作った密度行列に対応する行列式を考えれば $\braket{\Psi|\hat T|\Psi}\ge T_s$ が言えることによる．
ステップ2:試行波動関数としてスケールした $\Psi[\rho]$ を使う．$\Psi[\rho]$ をスケールした $\Psi[\rho]_\gamma$ は密度 $\rho_\gamma$ を与える(11.3.1節ステップ1と同じ計算)．したがって $F[\rho_\gamma]$ の制約付き最小化の候補である．よって最小性から
$$ \begin{equation} \braket{\Psi[\rho_\gamma]|\hat T+\hat V_{ee}|\Psi[\rho_\gamma]} \le\braket{\Psi[\rho]_\gamma|\hat T+\hat V_{ee}|\Psi[\rho]_\gamma} =\gamma^2\braket{\Psi[\rho]|\hat T|\Psi[\rho]}+\gamma\braket{\Psi[\rho]|\hat V_{ee}|\Psi[\rho]}. \label{eq:11-ec-scal-1} \end{equation} $$最後の等号では，運動エネルギーが $\gamma^2$ 倍，相互作用が $\gamma$ 倍になること(11.3.1節ステップ1(iii)とステップ2)を使った．
ステップ3:行列式側は等式である．11.3.1節ステップ1で示したとおり $\Phi[\rho_\gamma]=\Phi[\rho]_\gamma$ が厳密に成り立つので
$$ \begin{equation} \braket{\Phi[\rho_\gamma]|\hat T+\hat V_{ee}|\Phi[\rho_\gamma]} =\gamma^2 T_s[\rho]+\gamma\braket{\Phi[\rho]|\hat V_{ee}|\Phi[\rho]}. \label{eq:11-ec-scal-2} \end{equation} $$ステップ4:差をとる．\eqref{eq:11-ec-as-diff} に \eqref{eq:11-ec-scal-1} と \eqref{eq:11-ec-scal-2} を入れると
\begin{align} E_c[\rho_\gamma] &\le\gamma^2\bigl(\braket{\Psi|\hat T|\Psi}-T_s\bigr) +\gamma\bigl(\braket{\Psi|\hat V_{ee}|\Psi}-\braket{\Phi|\hat V_{ee}|\Phi}\bigr) \label{eq:11-ec-scal-3}\\ &=\gamma^2 T_c[\rho]+\gamma U_c[\rho] \label{eq:11-ec-scal-4}\\ &=\gamma\Bigl(T_c+U_c\Bigr)+\gamma(\gamma-1)T_c \label{eq:11-ec-scal-5}\\ &=\gamma E_c[\rho]+\gamma(\gamma-1)T_c[\rho]. \label{eq:11-ec-scal-6} \end{align}\eqref{eq:11-ec-scal-5} では $\gamma^2T_c=\gamma T_c+\gamma(\gamma-1)T_c$ と分解した．
ステップ5:$\gamma\le1$ の場合．$\gamma\gt0$ かつ $\gamma-1\le0$，$T_c\ge0$ なので $\gamma(\gamma-1)T_c\le0$．よって \eqref{eq:11-ec-scal-6} から
$$ \begin{equation} E_c[\rho_\gamma]\le\gamma E_c[\rho]\qquad(\gamma\le1). \label{eq:11-ec-scal-le} \end{equation} $$ステップ6:$\gamma\ge1$ の場合．いま得た \eqref{eq:11-ec-scal-le} を，密度 $\sigma\equiv\rho_\gamma$ とスケール因子 $1/\gamma\;(\le1)$ に適用する．$(\rho_\gamma)_{1/\gamma}=\rho_{\gamma\cdot(1/\gamma)}=\rho$ だから
$$ \begin{equation} E_c[\rho]=E_c\bigl[\sigma_{1/\gamma}\bigr]\le\frac{1}{\gamma}E_c[\sigma]=\frac{1}{\gamma}E_c[\rho_\gamma]. \label{eq:11-ec-scal-ge1} \end{equation} $$両辺に $\gamma\;(\gt0)$ を掛けて
$$ \begin{equation} \gamma E_c[\rho]\le E_c[\rho_\gamma]\qquad(\gamma\ge1). \label{eq:11-ec-scal-ge} \end{equation} $$ ∎物理的意味:相関は交換より「線形性が弱い」
交換は $\gamma$ に厳密に比例する \eqref{eq:11-ex-scaling} のに対し，相関は $\gamma$ に比例するよりゆるやかにしか大きくならない．密度を圧縮する($\gamma\to\infty$)と，電子はフェルミ縮退圧に支配されるようになり，クーロン相関の相対的重要性が下がるからである．実際，有限系では
$$ \lim_{\gamma\to\infty}E_c[\rho_\gamma]=E_c^{\mathrm{GL2}}[\rho]\quad(\text{有限の負の定数}) $$となることが知られている(Görling–Levyの2次摂動論)．LDAはこの制約を破る．一様電子ガスでは $\varepsilon_c\simeq A\ln r_s+B$ で $r_s\propto\rho^{-1/3}\propto\gamma^{-1}$ だから，$E_c^{\mathrm{LDA}}[\rho_\gamma]\simeq -AN\ln\gamma\to-\infty$ と対数発散してしまう．高密度極限で相関を過大評価するというLDAの欠陥は，表11.5で見る「原子の相関エネルギーを2倍以上に見積もる」という結果に直結している．11.4.5節で見るPBE相関の関数形は，まさにこの対数発散を打ち消すように設計されている．
11.3.3 Lieb–Oxford 下限
次は「$E_{xc}$ はどこまでも負にはなれない」という下からの押さえである．
定理:Lieb–Oxford 下限
任意の $N$ 電子密度 $\rho$ に対して
$$ \begin{equation} E_{xc}[\rho]\ \ge\ -C_{\mathrm{LO}}\int \dd^3r\;\rho^{4/3}(\rr), \qquad C_{\mathrm{LO}}\le1.68 \label{eq:11-LO} \end{equation} $$が成り立つ．(証明はLiebとOxfordによる評価であり，多体波動関数から古典静電エネルギーの下限を評価する精密な議論を要する．ここでは結果を認めて使う．厳密な最良定数は未確定で，$1.43\le C_{\mathrm{LO}}\le1.68$ の範囲にあることが知られている．)
この不等式の意味を「交換増強因子」の言葉に翻訳しておく．LDA交換 \eqref{eq:11-ex-lda} が $-C_x\int\rho^{4/3}$，$C_x=0.7386$ だったから，\eqref{eq:11-LO} は
$$ \begin{equation} \frac{E_{xc}[\rho]}{E_x^{\mathrm{LDA}}[\rho]}\le\frac{C_{\mathrm{LO}}}{C_x}=\frac{1.68}{0.7386}=2.273 \label{eq:11-LO-ratio} \end{equation} $$と書ける($E_x^{\mathrm{LDA}}\lt0$ で割ったので不等号の向きが反転していることに注意)．すなわち交換相関エネルギーは，LDA交換の $2.273$ 倍より深くなれない．これは11.4.4節でPBEのパラメータ $\kappa$ を決める根拠になる．
11.3.4 Levy–Perdew のビリアル型関係式
スケーリング等式 \eqref{eq:11-ex-scaling} を $\gamma$ について微分すると，$E_x$ と $v_x$ を結ぶ有用な関係式が出る．第2章のビリアル定理と同じ「座標スケーリングによる微分」の技法である．
導出:Levy–Perdew関係式 $E_x[\rho]=-\int\rho(\rr)\,\rr\cdot\nabla v_x(\rr)\,\dd^3r$
ステップ1:$\gamma=1$ で微分する．\eqref{eq:11-ex-scaling} の両辺を $\gamma$ で微分し $\gamma=1$ と置くと，右辺は $\dd(\gamma E_x[\rho])/\dd\gamma=E_x[\rho]$．左辺は汎関数の連鎖律により
$$ \begin{equation} \frac{\dd}{\dd\gamma}E_x[\rho_\gamma]\bigg|_{\gamma=1} =\int \dd^3r\;\frac{\delta E_x}{\delta\rho(\rr)}\bigg|_{\rho}\cdot \frac{\partial\rho_\gamma(\rr)}{\partial\gamma}\bigg|_{\gamma=1} =\int \dd^3r\;v_x(\rr)\,\frac{\partial\rho_\gamma(\rr)}{\partial\gamma}\bigg|_{\gamma=1}. \label{eq:11-lp-1} \end{equation} $$ステップ2:$\partial\rho_\gamma/\partial\gamma$ を計算する．$\rho_\gamma(\rr)=\gamma^3\rho(\gamma\rr)$ を $\gamma$ で微分する．第1因子からは $3\gamma^2\rho(\gamma\rr)$，第2因子からは連鎖律で $\gamma^3\,\rr\cdot(\nabla\rho)(\gamma\rr)$ が出る:
$$ \begin{equation} \frac{\partial\rho_\gamma(\rr)}{\partial\gamma}=3\gamma^2\rho(\gamma\rr)+\gamma^3\,\rr\cdot(\nabla\rho)(\gamma\rr) \;\xrightarrow{\ \gamma=1\ }\;3\rho(\rr)+\rr\cdot\nabla\rho(\rr). \label{eq:11-lp-2} \end{equation} $$したがって
$$ \begin{equation} E_x[\rho]=\int \dd^3r\;v_x(\rr)\bigl[3\rho(\rr)+\rr\cdot\nabla\rho(\rr)\bigr]. \label{eq:11-lp-3} \end{equation} $$ステップ3:第2項を部分積分する．ベクトル場 $\bm{A}(\rr)=\rr\,v_x(\rr)\rho(\rr)$ の発散をとる．積の発散の公式と $\nabla\cdot\rr=3$ から
\begin{align} \nabla\cdot\bigl[\rr\,v_x\rho\bigr] &=(\nabla\cdot\rr)\,v_x\rho+\rr\cdot\nabla(v_x\rho) \label{eq:11-lp-4}\\ &=3v_x\rho+\rho\,\rr\cdot\nabla v_x+v_x\,\rr\cdot\nabla\rho. \label{eq:11-lp-5} \end{align}これを全空間で積分する．左辺はガウスの定理により無限遠の球面積分になるが，束縛系では $\rho$ が指数的に減衰するので(そして $v_x$ は高々 $1/r$ でしか増大しないので)表面項はゼロである:
$$ \begin{equation} 0=\int \dd^3r\Bigl[3v_x\rho+\rho\,\rr\cdot\nabla v_x+v_x\,\rr\cdot\nabla\rho\Bigr]. \label{eq:11-lp-6} \end{equation} $$これを $\int v_x\,\rr\cdot\nabla\rho\,\dd^3r$ について解くと
$$ \begin{equation} \int \dd^3r\;v_x\,\rr\cdot\nabla\rho =-3\int \dd^3r\;v_x\rho-\int \dd^3r\;\rho\,\rr\cdot\nabla v_x. \label{eq:11-lp-7} \end{equation} $$ステップ4:代入する．\eqref{eq:11-lp-7} を \eqref{eq:11-lp-3} に入れると
\begin{align} E_x[\rho] &=3\int v_x\rho\;-3\int v_x\rho-\int \rho\,\rr\cdot\nabla v_x \label{eq:11-lp-8}\\ &=-\int \dd^3r\;\rho(\rr)\,\rr\cdot\nabla v_x(\rr). \label{eq:11-lp-9} \end{align} ∎この関係式は実用的な意味を持つ．交換ポテンシャル $v_x$ さえ分かればエネルギー $E_x$ が計算できるので，たとえば「厳密な $v_{xc}$ を逆問題として求める」タイプの計算(第10章で触れた逆KS法)や，軌道依存汎関数の検証に使われる．相関についても類似の関係式
$$ \begin{equation} E_c[\rho]+T_c[\rho]=-\int \dd^3r\;\rho(\rr)\,\rr\cdot\nabla v_c(\rr) \label{eq:11-lp-corr} \end{equation} $$が成り立つ(スケーリングが等式でなく不等式なので，$T_c$ という余分な項が現れる)．
11.3.5 自己相互作用の消去条件
最も直感的で，しかも近似汎関数がしばしば破る制約である．
定理:1電子系における自己相互作用の消去
電子が1個しかない系($N=1$，密度 $\rho$)では
$$ \begin{equation} E_{\mathrm{H}}[\rho]+E_{xc}[\rho]=0, \qquad\text{すなわち}\qquad E_x[\rho]=-E_{\mathrm{H}}[\rho],\quad E_c[\rho]=0 \label{eq:11-sic-condition} \end{equation} $$が厳密に成り立たねばならない．
証明．電子が1個なら電子間相互作用は存在しない($\hat V_{ee}=0$)．ところがKS分解では，古典クーロン項 $E_{\mathrm{H}}[\rho]=\frac12\iint\rho\rho'/\abs{\rr-\rr'}\gt0$ が必ず正の値で入ってしまう．これは「電子が自分自身の作る電荷分布と反発する」という物理的に存在しない相互作用(自己相互作用)である．定義上 $E_{xc}$ は「残り全部」だから，この偽の項をちょうど打ち消さねばならない．第5章5.5.2項で見たように，ハートリー・フォック法ではこの打ち消しが $J_{ll}=K_{ll}$ という厳密な等式として成立していた(交換項が同じ積分をそのまま差し引く)．KS法でも厳密な $E_{xc}$ を使う限り事情は同じだが，$E_{xc}$ を近似した瞬間に打ち消しは不完全になる——これが以下で見る自己相互作用誤差である．
ホールの言葉でも同じことが見える．$N=1$ では対密度 $n_2\equiv0$($i\ne j$ の組がない)なので，\eqref{eq:11-hole-def} から $n_{xc}(\rr,\rr')=0-\rho(\rr')=-\rho(\rr')$．これを \eqref{eq:11-exc-hole-key} に入れると
$$ E_{xc}=\frac{1}{2}\iint \dd^3r\,\dd^3r'\frac{\rho(\rr)\bigl[-\rho(\rr')\bigr]}{\abs{\rr-\rr'}}=-E_{\mathrm{H}}[\rho]. $$総和則 $\int n_{xc}\dd^3r'=-\int\rho\,\dd^3r'=-N=-1$ も自動的に満たされている．∎
例:水素原子におけるLSDAの自己相互作用誤差
水素原子の基底状態密度は $\rho(r)=e^{-2r}/\pi$(1電子，完全スピン分極 $\zeta=1$)である．各項を計算する．
(1) ハートリーエネルギー．1s密度に対する古典的自己エネルギーは初等的に
$$ E_{\mathrm{H}}=\frac{1}{2}\iint \dd^3r\,\dd^3r'\frac{\rho(\rr)\rho(\rr')}{\abs{\rr-\rr'}}=\frac{5}{16}=0.3125\ \mathrm{Ha}. $$厳密汎関数ならこれをちょうど打ち消す $E_{xc}=-0.3125$ Ha が必要である．
(2) LSDA交換．完全分極なので \eqref{eq:11-spin-scaling} より $E_x^{\mathrm{LSDA}}=\tfrac12E_x^{\mathrm{unpol}}[2\rho]=-2^{1/3}C_x\int\rho^{4/3}\dd^3r$．積分を実行する:
\begin{align} \int \dd^3r\,\rho^{4/3} &=\frac{1}{\pi^{4/3}}\int \dd^3r\;e^{-8r/3} =\frac{4\pi}{\pi^{4/3}}\int_0^\infty \dd r\;r^2e^{-8r/3} \label{eq:11-sic-int1}\\ &=\frac{4\pi}{\pi^{4/3}}\cdot\frac{2}{(8/3)^3} =\frac{4\pi}{4.6011}\cdot\frac{2\cdot27}{512}=0.28805. \label{eq:11-sic-int2} \end{align}\eqref{eq:11-sic-int1} では球対称性から $\dd^3r=4\pi r^2\dd r$，\eqref{eq:11-sic-int2} では公式 $\int_0^\infty r^2e^{-ar}\dd r=2/a^3$($a=8/3$)を使った．したがって
$$ E_x^{\mathrm{LSDA}}=-2^{1/3}\times0.738559\times0.28805=-0.2680\ \mathrm{Ha}. $$(3) LSDA相関．厳密には $E_c=0$ でなければならないが，LSDAは $E_c^{\mathrm{LSDA}}\simeq-0.022$ Ha を与えてしまう(密度 $\rho(r)$ に対応する $r_s(r)$ で $\varepsilon_c(r_s,\zeta=1)$ を積分した値)．
(4) 差し引き．
$$ E_{\mathrm{H}}+E_{xc}^{\mathrm{LSDA}}=0.3125-0.2680-0.022=+0.022\ \mathrm{Ha}=0.60\ \mathrm{eV}. $$厳密にはゼロであるべき量が $0.6$ eV も残る．しかもこれは正，すなわちLSDAは電子を余計に不安定化している．誤差の内訳を見ると，交換が自己ハートリーの $86\%$ しか打ち消せない一方，あるはずのない相関が余分に加わっている．この2つは部分的に相殺するので，結果として $0.02$ Ha 程度に収まっているが，これも一種の偶然の相殺である．
実際，自己無撞着に解いたLSDAの水素原子の全エネルギーは $-0.479$ Ha であり，厳密値 $-0.5$ Ha より $0.021$ Ha 高い．上の見積もりとよく合う．
物理的意味:自己相互作用誤差がもたらす3つの症状
自己相互作用誤差(SIE: self-interaction error)は多電子系にも「各電子が少しずつ自分と反発している」という形で残り，次の症状を生む．
- 占有準位が浅くなる．電子は本来 $N-1$ 個の他電子を感じるべきなのに，$N$ 個分の反発を感じる．その結果，占有KS準位が押し上げられ，$\varepsilon_{\mathrm{HOMO}}$ は $-I$ より浅くなり，バンドギャップが過小評価される(第10章)．
- 局在した状態が不当に不安定になる．電子が狭い領域に集まるほど自己ハートリーが大きくなるので，SIEはそれを嫌う．$d$・$f$ 電子の局在，電荷移動状態，遷移状態(反応障壁)などが不当に高いエネルギーを与えられる．これが反応障壁の系統的な過小評価の原因である．
- 電荷の非局在化誤差．分数電子数に対するエネルギーが下に凸になり，電荷を複数のフラグメントに分散させたほうが得になってしまう($\mathrm{H}_2^+$ の解離の破綻，第10章の演習)．
SIEを消す処方には，Perdew–Zungerの自己相互作用補正(軌道ごとに自己ハートリー+自己交換を差し引く)，厳密交換を混ぜるハイブリッド汎関数，DFT+Uなどがあり，いずれも11.6節で扱う．
11.3.6 制約カタログ
ここまでの制約を，LDAとPBE-GGAがそれぞれ満たすかどうかとともに整理しておく．
| 制約 | 内容 | LDA | PBE |
|---|---|---|---|
| 交換ホールの符号 | $n_x(\rr,\rr')\le0$ | ◯ | ◯ |
| 交換ホールの総和則 | $\int n_x\dd^3r'=-1$ | ◯ | ◯ |
| 相関ホールの総和則 | $\int n_c\dd^3r'=0$ | ◯ | ◯ |
| 交換のスピンスケーリング | $E_x[\rho^\uparrow,\rho^\downarrow]=\frac12(E_x[2\rho^\uparrow]+E_x[2\rho^\downarrow])$ | ◯ | ◯ |
| 交換の一様スケーリング | $E_x[\rho_\gamma]=\gamma E_x[\rho]$ | ◯ | ◯ |
| 相関のスケーリング不等式 | 式 \eqref{eq:11-ec-scaling} | △ | ◯ |
| 高密度極限の有限性 | $E_c[\rho_\gamma]$ が $\gamma\to\infty$ で有限 | ×(対数発散) | ◯ |
| 低密度極限 | $\gamma\to0$ で $E_c[\rho_\gamma]\to$ 線形 | ◯ | ◯ |
| Lieb–Oxford下限 | $E_{xc}\ge-1.68\int\rho^{4/3}$ | ◯ | ◯ |
| 一様ガス極限 | $\rho=$ 一定で厳密値 | ◯ | ◯ |
| 線形応答($s\to0$) | 一様ガスの応答を再現 | ◯ | ◯ |
| 自己相互作用の消去 | $N=1$ で $E_{xc}=-E_{\mathrm{H}}$ | × | × |
| 汎関数微分の不連続 | $\Delta_{xc}\ne0$ | × | × |
| $v_{xc}$ の漸近形 | $v_{xc}\to-1/r$($r\to\infty$) | ×(指数的に消える) | × |
下の3つの×が，LDAとGGAに共通して残る本質的な欠陥である．いずれも局所汎関数であることそのものに起因しており，勾配を追加する程度では直らない．これらを直すには軌道を陽に使う汎関数(ハイブリッド，meta-GGA)へ上がる必要があり，それが11.6節の主題になる．
この節で得られたもの．厳密汎関数が満たす制約——スケーリング等式・不等式，Lieb–Oxford下限，ビリアル型関係式，自己相互作用の消去——を導いた．次節では，これらの一部を「必ず満たす」ように関数形を選ぶことで，経験的パラメータを持たない汎関数PBEを構成する．
11.4 一般化勾配近似(GGA)とPBE
LDAの次に自然な一手は「密度の勾配 $\nabla\rho$ も使う」ことである．しかし素朴な勾配展開は失敗し，そこから正しい教訓を汲み取るのに20年近くを要した．本節ではその歴史をたどり，現代の標準汎関数PBEを厳密条件から組み立てる．
11.4.1 勾配展開(GEA)とその失敗
密度が空間的にゆっくり変化する系を考え，$E_{xc}$ を勾配の冪で展開する:
$$ \begin{equation} E_{xc}^{\mathrm{GEA}}[\rho]=\int \dd^3r\Bigl[ \rho\,\varepsilon_{xc}(\rho) +C_{xc}(\rho)\frac{\abs{\nabla\rho}^2}{\rho^{4/3}} +O(\nabla^4)\Bigr]. \label{eq:11-gea} \end{equation} $$1次の項 $\nabla\rho$ が現れないのは，$E_{xc}$ がスカラーであり，また空間反転($\rr\to-\rr$)に対して不変でなければならないからである($\nabla\rho$ はベクトルなので単独では不変量を作れない)．これを勾配展開近似(GEA: Gradient Expansion Approximation)と呼ぶ．係数は一様電子ガスの線形応答から解析的に決まり，交換については後で見る $10/81$ が出る．
理屈のうえでは，LDAは展開の0次，GEAはその1段上である．したがってGEAはLDAより良いはずである．ところが実際に原子や分子に適用すると，GEAはLDAより悪い結果を与える．原子の交換エネルギーの誤差は改善するどころか悪化し，分子の結合エネルギーは大幅に過大評価される．
原因は11.2節の言葉で明快に説明できる．GEAの2次項は，対応するホールの2次補正として書き直すことができるが，その補正されたホールは総和則を破るのである:
$$ \begin{equation} \int \dd^3r'\;n_x^{\mathrm{GEA}}(\rr,\rr')\ne-1. \end{equation} $$さらに，GEAホールは遠距離で符号の振動する長い尾を持ち，$n_x\le0$ という符号条件も破る．密度変化が本当にゆっくりならこれらの破れは微小で問題にならないが，実在系のように $\abs{\nabla\rho}/\rho$ が大きい領域では破れが大きくなり，11.2節で見た「総和則がエネルギーの目盛りを固定する」という保護機構が失われてしまう．LDAが良かったのは0次だからではなく，総和則を満たしていたからである．
物理的意味:展開の次数より制約の充足
これは汎関数開発の歴史における最も重要な教訓である．「系統的な展開の高次項を足せば必ず良くなる」という直観は，$E_{xc}$ については成り立たない．$E_{xc}$ の値を守っているのは展開の収束ではなく，ホールが満たす厳密な積分制約だからである．近似を改良するときは，まず制約を壊さないことを確認しなければならない．この方針が徹底されたのが11.4.4節以降のPBEである．
11.4.2 実空間カットオフによるGGAの再生
PerdewとWangは1986年，この診断から直接的な処方を提案した．GEAのホールを捨てるのではなく，手で修理するのである．
- GEAのホール $n_x^{\mathrm{GEA}}(\rr,\bm{s})$ を出発点とする．
- 符号条件を回復するため，$n_x^{\mathrm{GEA}}\gt0$ となる領域では強制的に $n_x=0$ とおく．
- 総和則を回復するため，ある半径 $s_c(\rr)$ より遠方を切り捨てる．カットオフ半径 $s_c$ は $\int_{s\lt s_c}n_x\dd^3s=-1$ がちょうど成り立つように点ごとに決める．
- こうして作った「修理済みホール」を \eqref{eq:11-sph-final} に代入して $E_x$ を数値的に評価する．
この手続きの結果は解析式では書けないが，数値的に求めた値は「$\rho$ と $\abs{\nabla\rho}$ の局所的な関数」として整理できる．それを解析形にフィットしたものが一般化勾配近似(GGA: Generalized Gradient Approximation)である．「一般化」とは「勾配の冪級数ではなく，勾配の一般の関数を許す」という意味である．
$$ \begin{equation} E_{xc}^{\mathrm{GGA}}[\rho]=\int \dd^3r\;f\bigl(\rho(\rr),\nabla\rho(\rr)\bigr) \label{eq:11-gga-general} \end{equation} $$PW86(1986)，B88(Becke，1988)，PW91(1991)，PBE(1996)といった汎関数がこの系譜に属する．B88とLYPを組み合わせたBLYPは量子化学で，PBEは物性物理で標準的な地位を占めている．
11.4.3 無次元換算勾配 $s$
GGAを構成する前に，勾配の「大きさ」を測る適切な物差しを定義する．
定義:無次元換算勾配(reduced density gradient)
$$ \begin{equation} s(\rr)\equiv\frac{\abs{\nabla\rho(\rr)}}{2k_F(\rr)\,\rho(\rr)}, \qquad k_F=\bigl(3\pi^2\rho\bigr)^{1/3} \label{eq:11-s-def} \end{equation} $$$s$ が無次元であることを確認しよう．原子単位で長さの次元を $L$ とすると，$[\rho]=L^{-3}$，$[\nabla\rho]=L^{-4}$，$[k_F]=L^{-1}$ だから
$$ [s]=\frac{L^{-4}}{L^{-1}\cdot L^{-3}}=L^0. $$すなわち無次元である．分母の因子2は歴史的な慣習だが，$2k_F$ が「フェルミ球の直径」であることに対応している．
$s$ の物理的意味は3通りに読める．
- フェルミ波長あたりの密度変化．$1/(2k_F)$ はおよそフェルミ波長の $1/4\pi$ 倍であり，$s=\lambda_F\abs{\nabla\rho}/(4\pi\rho)$ と書ける．LDAの適用条件 \eqref{eq:11-lda-condition} は $4\pi s\ll1$，すなわち $s\ll0.08$ と等価である．
- スケーリング不変量．11.3.1節の例で示したとおり，一様スケーリングで $s$ の値は保たれる．したがって $s$ の関数で汎関数を作れば，交換のスケーリング等式が自動的に満たされる．
- 密度の「不均一さ」の目安．$s=0$ は一様電子ガス，$s\lesssim1$ は緩やかに変化する価電子領域，$s\gg1$ は密度が急変する領域(原子の裾，真空との界面)である．
例:実在系での $s$ の値
- 原子の裾．$\rho\propto e^{-2\alpha r}$ とすると $\abs{\nabla\rho}/\rho=2\alpha$(一定)，$k_F\propto\rho^{1/3}\propto e^{-2\alpha r/3}$．よって $s\propto e^{+2\alpha r/3}\to\infty$ と指数的に発散する．したがって $F_x(s)$ は $s\to\infty$ でも有界でなければならない——これが11.4.4節の条件(ii)である．
- 単純金属の価電子．ほぼ一様なので $s\simeq0$〜$0.3$．ここではGGAはLDAとほとんど変わらない．
- 共有結合の結合中心．$s\simeq0.3$〜$1$．GGAが効き始める領域．
- 原子核の近傍．密度が非常に大きいので $k_F$ も大きく，勾配が急でも $s$ は $0.5$ 程度に抑えられる．
- 分子間の隙間・表面の外側．$s\gtrsim3$．GGAが $F_x$ の飽和値に近づく領域で，汎関数間の差が最も大きくなる．
11.4.4 PBE交換:3つの条件から関数形を決める
これから，PerdewとBurkeとErnzerhofが1996年に提案した交換汎関数を組み立てる．出発点は，11.3.1節で見た「スケーリング等式を自動的に満たす」形である:
定義:交換増強因子
$$ \begin{equation} E_x^{\mathrm{GGA}}[\rho]=\int \dd^3r\;\rho(\rr)\,\varepsilon_x^{\mathrm{unif}}\bigl(\rho(\rr)\bigr)\,F_x\bigl(s(\rr)\bigr) =-C_x\int \dd^3r\;\rho^{4/3}(\rr)\,F_x\bigl(s(\rr)\bigr) \label{eq:11-Fx-def} \end{equation} $$$F_x(s)$ を交換増強因子と呼ぶ．$F_x\equiv1$ がLDAである．$F_x\gt1$ は「LDAより交換を深くする」ことを意味する．
あとは1変数関数 $F_x(s)$ を決めればよい．PBEはこれを次の3条件から決める．
導出:PBE交換増強因子
条件(i):一様極限 $F_x(0)=1$．密度が一様なら $s=0$ であり，そこではLDA(=厳密解)に戻らねばならない．
条件(ii):$s\to\infty$ で有界．上の例で見たとおり，原子の裾では $s\to\infty$ になる．そこで $F_x$ が発散すると(GEAの $1+\mu s^2$ はまさに発散する)，密度が指数的に小さい領域でエネルギー密度が発散的に大きくなり，数値的にも物理的にも破綻する．さらに，Lieb–Oxford下限 \eqref{eq:11-LO-ratio} を局所的に課すことで上限値が決まる．
ここで注意が要る．\eqref{eq:11-LO-ratio} は $E_{xc}$ 全体に対する制約なので，$F_x$ 単独に $2.273$ を課すのは正しくない．PBEは最も厳しい条件——完全スピン分極($\zeta=1$)かつ高密度($r_s\to0$，相関が交換に比べて無視できる)——で評価する．スピンスケーリング関係 \eqref{eq:11-spin-scaling} より，完全分極系では実質的に密度 $2\rho$ の無分極系の交換を半分にしたものが効く．$2\rho$ に対する換算勾配は
$$ s[2\rho]=\frac{\abs{\nabla(2\rho)}}{2(3\pi^2\cdot2\rho)^{1/3}(2\rho)} =\frac{2\abs{\nabla\rho}}{2\cdot2^{1/3}k_F\cdot2\rho}=2^{-1/3}s $$であり，また $E_x^{\mathrm{unpol}}[2\rho]\propto-(2\rho)^{4/3}=-2^{4/3}\rho^{4/3}$ である．したがって完全分極系の増強比は，無分極LDA交換を基準にすると
$$ \frac{E_x^{\zeta=1}}{E_x^{\mathrm{LDA},\zeta=0}} =\frac{\tfrac12\cdot2^{4/3}}{1}F_x\bigl(2^{-1/3}s\bigr)=2^{1/3}F_x\bigl(2^{-1/3}s\bigr). $$これが $2.273$ 以下であればよいから，任意の引数 $\tilde s=2^{-1/3}s$ について
$$ \begin{equation} F_x(\tilde s)\le\frac{2.273}{2^{1/3}}=\frac{2.273}{1.2599}=1.804. \label{eq:11-kappa-derive} \end{equation} $$そこでPBEは $F_x$ の上限を $1+\kappa$ と書き，
$$ \begin{equation} \kappa=0.804 \label{eq:11-kappa} \end{equation} $$と定める．これがLieb–Oxford下限から $\kappa$ が決まる筋道である．
条件(iii):$s\to0$ の展開係数．小さい $s$ での振る舞いを $F_x(s)=1+\mu s^2+O(s^4)$ と書く．GEAの値は $\mu_{\mathrm{GE}}=10/81=0.1235$ だが，PBEはあえてこれを採らない．理由は次の11.4.5節で相関と合わせて説明するが，結論だけ先に書くと
$$ \begin{equation} \mu=\beta\frac{\pi^2}{3}=0.066725\times3.28987=0.21951 \label{eq:11-mu} \end{equation} $$である($\beta$ は相関の勾配展開係数)．
関数形の決定．以上3条件，すなわち
$$ F_x(0)=1,\qquad F_x(s)\simeq1+\mu s^2\ (s\to0),\qquad F_x(\infty)=1+\kappa $$を満たす最も簡単な有理関数を選ぶ:
$$ \begin{equation} F_x^{\mathrm{PBE}}(s)=1+\kappa-\frac{\kappa}{1+\mu s^2/\kappa}. \label{eq:11-Fx-pbe} \end{equation} $$検算．(a) $s=0$:$F_x=1+\kappa-\kappa/1=1$ ◯．(b) $s\to\infty$:第2項の分母が発散して第2項がゼロになり $F_x\to1+\kappa$ ◯．(c) 小さい $s$:$x\equiv\mu s^2/\kappa$ として $1/(1+x)=1-x+x^2-\cdots$ を使うと
\begin{align} F_x&=1+\kappa-\kappa\left(1-\frac{\mu s^2}{\kappa}+\frac{\mu^2s^4}{\kappa^2}-\cdots\right) \label{eq:11-Fx-exp-1}\\ &=1+\kappa-\kappa+\mu s^2-\frac{\mu^2}{\kappa}s^4+\cdots \label{eq:11-Fx-exp-2}\\ &=1+\mu s^2+O(s^4)\ \text{◯}. \label{eq:11-Fx-exp-3} \end{align} ∎
11.4.5 PBE相関と，$\mu$ の由来
相関部分も同じ思想で作られる．まずLDA相関に「勾配補正 $H$」を加える形を置く:
$$ \begin{equation} E_c^{\mathrm{PBE}}[\rho^\uparrow,\rho^\downarrow] =\int \dd^3r\;\rho\Bigl[\varepsilon_c^{\mathrm{unif}}(r_s,\zeta)+H(r_s,\zeta,t)\Bigr] \label{eq:11-pbe-c} \end{equation} $$ここで $t$ は相関に適した換算勾配で，遮蔽波数 $k_s=\sqrt{4k_F/\pi}$(第8章でトーマス・フェルミ遮蔽波数 $\kappa_{\mathrm{TF}}$ と書いたものと同じ量．ハートリー原子単位で $\kappa_{\mathrm{TF}}^2=4k_F/\pi$)を物差しに使う:
$$ \begin{equation} t\equiv\frac{\abs{\nabla\rho}}{2\phi(\zeta)k_s\rho}, \qquad k_s=\sqrt{\frac{4k_F}{\pi}}, \qquad \phi(\zeta)=\frac{(1+\zeta)^{2/3}+(1-\zeta)^{2/3}}{2}. \label{eq:11-t-def} \end{equation} $$交換では $k_F$(パウリ相関の物差し)，相関では $k_s$(遮蔽の物差し)を使う——この使い分けが物理的に自然である．$\phi(\zeta)$ はスピン分極による遮蔽長の変化を補正する因子で，$\phi(0)=1$，$\phi(1)=2^{-1/3}$ である．
$H$ の関数形は次の3条件から決まる．以下，$\gamma_{\mathrm{PBE}}=(1-\ln2)/\pi^2=0.031091$，$\beta=0.066725$ とする．
導出:PBE相関の関数形と3条件の検証
PBEが採用する形は
$$ \begin{equation} H=\gamma_{\mathrm{PBE}}\,\phi^3\ln\left[1+\frac{\beta}{\gamma_{\mathrm{PBE}}}t^2 \frac{1+At^2}{1+At^2+A^2t^4}\right], \qquad A=\frac{\beta}{\gamma_{\mathrm{PBE}}} \left[\exp\!\left(-\frac{\varepsilon_c^{\mathrm{unif}}}{\gamma_{\mathrm{PBE}}\phi^3}\right)-1\right]^{-1} \label{eq:11-H-pbe} \end{equation} $$である．一見複雑だが，次の3条件をすべて満たす最小限の構造になっている．
条件(a):緩やかに変化する極限 $t\to0$ で勾配展開に戻る．$t$ が小さいとき，対数の引数の第2項は小さい．$\ln(1+x)\simeq x$($x\to0$)と，括弧内の分数が $t\to0$ で $1$ になることを使うと
$$ \begin{equation} H\to\gamma_{\mathrm{PBE}}\phi^3\cdot\frac{\beta}{\gamma_{\mathrm{PBE}}}t^2=\beta\phi^3t^2. \label{eq:11-H-smallt} \end{equation} $$これが相関の2次勾配展開の正しい形である($\beta=0.066725$ はその係数)．
条件(b):急激に変化する極限 $t\to\infty$ で相関が消える．密度が急変する領域では相関エネルギーはゼロに向かうべきである(電子は互いを避ける余裕がない)．$t\to\infty$ で分数は
$$ \frac{1+At^2}{1+At^2+A^2t^4}\simeq\frac{At^2}{A^2t^4}=\frac{1}{At^2} $$となり，対数の引数は
$$ 1+\frac{\beta}{\gamma_{\mathrm{PBE}}}t^2\cdot\frac{1}{At^2}=1+\frac{\beta}{\gamma_{\mathrm{PBE}}A}. $$ここで $A$ の定義 \eqref{eq:11-H-pbe} を代入する．$\dfrac{\beta}{\gamma_{\mathrm{PBE}}A}=\exp\!\left(-\dfrac{\varepsilon_c^{\mathrm{unif}}}{\gamma_{\mathrm{PBE}}\phi^3}\right)-1$ だから，対数の引数はちょうど
$$ 1+\left[\exp\!\left(-\frac{\varepsilon_c^{\mathrm{unif}}}{\gamma_{\mathrm{PBE}}\phi^3}\right)-1\right] =\exp\!\left(-\frac{\varepsilon_c^{\mathrm{unif}}}{\gamma_{\mathrm{PBE}}\phi^3}\right) $$となる．よって
$$ \begin{equation} H\to\gamma_{\mathrm{PBE}}\phi^3\cdot\left(-\frac{\varepsilon_c^{\mathrm{unif}}}{\gamma_{\mathrm{PBE}}\phi^3}\right) =-\varepsilon_c^{\mathrm{unif}}, \label{eq:11-H-larget} \end{equation} $$すなわち \eqref{eq:11-pbe-c} の被積分関数の括弧内が $\varepsilon_c^{\mathrm{unif}}-\varepsilon_c^{\mathrm{unif}}=0$ になり，相関が完全に消える．$A$ の一見奇妙な定義は，この条件(b)を満たすためだけに逆算されたものである．
条件(c):高密度極限($\gamma\to\infty$ のスケーリング)で $E_c$ が有限にとどまる．11.3.2節で見たとおり，LDAはここで対数発散する．PBEがどう救うかを見よう．高密度($r_s\to0$)では $\varepsilon_c^{\mathrm{unif}}\simeq\gamma_{\mathrm{PBE}}\ln r_s\to-\infty$ である．すると $-\varepsilon_c^{\mathrm{unif}}/(\gamma_{\mathrm{PBE}}\phi^3)\to+\infty$ となり，$A$ の定義の指数関数が発散するので
$$ A\to0. $$$A\to0$ では分数 $\dfrac{1+At^2}{1+At^2+A^2t^4}\to1$ なので
$$ \begin{equation} H\to\gamma_{\mathrm{PBE}}\phi^3\ln\left[1+\frac{\beta}{\gamma_{\mathrm{PBE}}}t^2\right]. \label{eq:11-H-highdens} \end{equation} $$いま一様スケーリング $\rho\to\rho_\gamma$ を行うと，$t^2\propto\abs{\nabla\rho}^2/(k_s^2\rho^2)\propto\abs{\nabla\rho}^2/(k_F\rho^2)\propto\rho^{1/3}\times(\text{無次元量})$ なので $t^2\propto\gamma$ と増大する．したがって $\ln$ の中が大きくなり
$$ H\simeq\gamma_{\mathrm{PBE}}\ln t^2+\text{const}\simeq\gamma_{\mathrm{PBE}}\ln\gamma+\text{const}. $$一方 $\varepsilon_c^{\mathrm{unif}}\simeq\gamma_{\mathrm{PBE}}\ln r_s=-\gamma_{\mathrm{PBE}}\ln\gamma+\text{const}$ である($r_s\propto1/\gamma$)．両者の和
$$ \varepsilon_c^{\mathrm{unif}}+H\simeq-\gamma_{\mathrm{PBE}}\ln\gamma+\gamma_{\mathrm{PBE}}\ln\gamma+\text{const}=\text{const} $$では対数がちょうど打ち消し合う．これでPBEは高密度極限で有限にとどまり，LDAが破っていた制約(表11.2)を満たす．$H$ が $\ln$ の形をしているのは，まさにこの相殺を仕込むためである．∎
最後に，条件(iii) で保留していた $\mu$ の値を決める．
導出:$\mu=\beta\pi^2/3$——なぜGEAの $10/81$ を使わないのか
ステップ1:方針．一様電子ガスに緩やかな摂動を加えたときの線形応答は，実はLDAの方がGEAより正確であることが知られている．これは偶然ではなく，交換のGEA補正と相関のGEA補正が互いに打ち消し合う傾向があるためである．PBEはこの事実を制約として採用する:「$s\to0$($t\to0$)の極限で，交換と相関の勾配補正の和がゼロになるように $\mu$ を選ぶ」．
ステップ2:$t^2$ を $s^2$ で表す．スピン無分極($\phi=1$)として，定義 \eqref{eq:11-s-def}，\eqref{eq:11-t-def} の比をとる:
\begin{align} \frac{t^2}{s^2} &=\frac{\bigl[\abs{\nabla\rho}/(2k_s\rho)\bigr]^2}{\bigl[\abs{\nabla\rho}/(2k_F\rho)\bigr]^2} =\frac{k_F^2}{k_s^2} \label{eq:11-mu-1}\\ &=\frac{k_F^2}{4k_F/\pi}=\frac{\pi k_F}{4}. \label{eq:11-mu-2} \end{align}\eqref{eq:11-mu-2} で $k_s^2=4k_F/\pi$ を代入した．すなわち $t^2=\dfrac{\pi k_F}{4}s^2$．
ステップ3:相関の勾配補正を書く．条件(a) \eqref{eq:11-H-smallt} より，勾配補正の相関エネルギーへの寄与は
$$ \begin{equation} \Delta E_c=\int \dd^3r\;\rho\,\beta t^2 =\int \dd^3r\;\rho\,\beta\,\frac{\pi k_F}{4}s^2. \label{eq:11-mu-3} \end{equation} $$ステップ4:交換の勾配補正を書く．$F_x\simeq1+\mu s^2$ より，LDAからのずれは
$$ \begin{equation} \Delta E_x=\int \dd^3r\;\rho\,\varepsilon_x^{\mathrm{unif}}(\rho)\,\mu s^2 =\int \dd^3r\;\rho\left(-\frac{3k_F}{4\pi}\right)\mu s^2. \label{eq:11-mu-4} \end{equation} $$ここで $\varepsilon_x^{\mathrm{unif}}=-3k_F/(4\pi)$(式 \eqref{eq:11-eps-x-kf})を使った．
ステップ5:和をゼロにする．
\begin{align} \Delta E_x+\Delta E_c &=\int \dd^3r\;\rho\,s^2\left[\frac{\beta\pi k_F}{4}-\frac{3\mu k_F}{4\pi}\right] \label{eq:11-mu-5}\\ &=\int \dd^3r\;\rho\,s^2\,\frac{k_F}{4}\left[\beta\pi-\frac{3\mu}{\pi}\right]. \label{eq:11-mu-6} \end{align}これが任意の密度でゼロになるためには，角括弧の中がゼロでなければならない:
$$ \begin{equation} \beta\pi=\frac{3\mu}{\pi} \quad\Longleftrightarrow\quad \mu=\frac{\beta\pi^2}{3}=\frac{0.066725\times9.8696}{3}=0.21951. \label{eq:11-mu-final} \end{equation} $$ ∎物理的意味:$\mu$ の選択は「どちらの極限を大事にするか」の選択である
$\mu$ には2つの候補がある．
- $\mu=10/81=0.1235$:交換だけの勾配展開に忠実．緩やかに変化する密度(単純金属，固体のバルク)を正しく記述する．
- $\mu=0.2195$:交換相関の和について線形応答をLDAに合わせる．原子・分子のように密度が急変する系で良好．
PBEは後者を採った．その結果，PBEは分子の結合エネルギーや原子の全エネルギーを非常によく再現する一方，固体の格子定数をやや過大評価する傾向を持つ(11.5節)．この欠点を直すために2008年に提案されたのがPBEsolで，$\mu=10/81$，$\beta=0.046$ に戻したものである(図11.5)．PBEsolは固体の格子定数と表面エネルギーを改善するが，分子の結合エネルギーは悪化する．どの制約を優先するかによって最適な汎関数が変わる——これが「万能の汎関数は存在しない」ことの具体的な現れである．
11.4.6 非経験性の意義
PBEに現れるパラメータを改めて並べる．
| パラメータ | 値 | 決定根拠 |
|---|---|---|
| $\kappa$ | 0.804 | Lieb–Oxford下限 \eqref{eq:11-LO} を完全分極・高密度で局所的に課す(式 \eqref{eq:11-kappa-derive}) |
| $\beta$ | 0.066725 | 一様電子ガスの相関の2次勾配展開係数(高密度極限での解析的な値) |
| $\mu$ | $\beta\pi^2/3=0.21951$ | 交換相関の和が一様ガスの線形応答をLDA通りに再現する条件(式 \eqref{eq:11-mu-final}) |
| $\gamma_{\mathrm{PBE}}$ | $(1-\ln2)/\pi^2=0.031091$ | 一様電子ガスの相関エネルギーの高密度極限の対数の係数 \eqref{eq:11-ec-high} |
4つとも，一様電子ガスの厳密な性質か，厳密な不等式かのどちらかから決まっている．PBEは「経験的パラメータを持たない(non-empirical)汎関数」である．これは単なる美学ではない．
物理的意味:なぜ非経験性が重要か
- 予測性．実験データにフィットした汎関数は，フィットに使ったデータセットに似た系(たとえば第2周期元素からなる有機分子)では優秀だが，そこから外れた系(遷移金属，高圧下の物質，未知の化合物)での信頼性が保証されない．厳密条件だけから作られた汎関数は，どの系についても同じ論理で構成されている．
- 誤りの診断可能性．PBEが特定の系で失敗したとき，「どの制約が足りないのか」を問うことができる．実際，vdWの欠如は「非局所相関の制約が入っていない」，ギャップの過小評価は「汎関数微分の不連続が入っていない」と診断でき，次の汎関数の設計指針になる．フィット型の汎関数では，失敗の原因を切り分けにくい．
- 系統的改良の道筋．より多くの厳密条件を満たす関数形を探す，という明確な目標が立つ．11.6節で紹介するSCAN(meta-GGA)は，この方針を徹底して「meta-GGAが満たしうる17個の既知の厳密条件をすべて満たす」ように構成されている．
一方，量子化学で広く使われるB3LYPのように，3つのパラメータを分子のデータセットにフィットした汎関数もあり，対象を限れば極めて高精度である．どちらが「正しい」かではなく，どこまで外挿できるかが違うと理解するのが正しい．
この節で得られたもの．勾配展開の失敗が総和則の破れに起因することを見たうえで，実空間カットオフによってGGAが再生された歴史をたどり，PBE汎関数を厳密条件だけから組み立てた．次節では，LDAとGGAが実際にどれだけ違う数値を与えるかを見る．
11.5 LDAとGGAを数値で比べる
ここまでは理論だけを追ってきた．実際に使うときに知りたいのは「どちらをどれだけ信用してよいか」である．本節では原子・分子・固体それぞれについて代表的な数値を並べ，それぞれの汎関数がどの方向に誤るかという系統性を読み取る．系統性さえ分かっていれば，計算結果を批判的に読むことができる．
11.5.1 原子:交換エネルギーと相関エネルギー
原子は $E_{xc}$ の近似を試す最も清潔な舞台である．厳密な(あるいは厳密に近い)$E_x$，$E_c$ の値が別途知られているからである．
| 原子 | 厳密 | LSDA | 誤差 | PBE | 誤差 |
|---|---|---|---|---|---|
| H | 0.3125 | 0.2680 | $-14.2\%$ | 0.3059 | $-2.1\%$ |
| He | 1.0258 | 0.8840 | $-13.8\%$ | 1.0136 | $-1.2\%$ |
| Be | 2.6658 | 2.3124 | $-13.3\%$ | 2.6358 | $-1.1\%$ |
| N | 6.6011 | 5.9008 | $-10.6\%$ | 6.5521 | $-0.7\%$ |
| Ne | 12.1050 | 11.0335 | $-8.9\%$ | 12.0667 | $-0.3\%$ |
| 原子 | 厳密 | LSDA | PBE |
|---|---|---|---|
| H | 0.0000 | 0.0222 | 0.0060 |
| He | 0.0420 | 0.1125 | 0.0420 |
| Be | 0.0950 | 0.2240 | 0.0856 |
| N | 0.1883 | 0.4268 | 0.1799 |
| Ne | 0.3939 | 0.7428 | 0.3513 |
物理的意味:2つの誤差は逆向きに効く
表11.4と表11.5を並べると，LSDAの誤りが2つあり，しかも互いに逆向きであることが分かる．
- 交換を浅く見積もる(Neで $-1.07$ Ha)．原子内では密度が急変するため，一様ガスのホールでは交換を十分深く掘れない．
- 相関を深く見積もる(Neで $+0.35$ Ha)．11.3.2節で見たとおり，高密度極限で対数発散する欠陥の現れである．
Neでは $-1.07+0.35=-0.72$ Ha が正味の誤差として残る．それでも，個別の誤差より小さくなっている．LDAの「実用的な精度」の一部は，この $E_x$ と $E_c$ の誤差相殺に支えられている．したがって「交換だけをより良い汎関数にする」といった中途半端な改良はしばしば結果を悪化させる．$E_x$ と $E_c$ は同じ論理で構成された組(たとえばPBExとPBEc)で使わねばならない．
PBEは両方の誤差を1桁縮めており，相殺に頼る度合いが減っている．これが「制約を満たすように作る」ことの実利である．
11.5.2 分子:LDAの過結合
分子の解離エネルギー(原子化エネルギー)は，DFTの実用性が最初に問われた量である．傾向ははっきりしている．
| 分子 | 実験 | LSDA | 誤差 | PBE | 誤差 |
|---|---|---|---|---|---|
| H$_2$ | 4.75 | 4.90 | $+0.15$ | 4.55 | $-0.20$ |
| N$_2$ | 9.91 | 11.60 | $+1.69$ | 10.55 | $+0.64$ |
| O$_2$ | 5.23 | 7.59 | $+2.36$ | 6.24 | $+1.01$ |
| F$_2$ | 1.66 | 3.35 | $+1.69$ | 2.31 | $+0.65$ |
| H$_2$O | 10.09 | 11.58 | $+1.49$ | 10.19 | $+0.10$ |
| CO | 11.24 | 12.94 | $+1.70$ | 11.66 | $+0.42$ |
なぜLDAは過結合するのか．原子化エネルギーは「分子のエネルギー」と「バラバラの原子のエネルギーの和」の差である．孤立原子の密度は，分子の密度よりも急峻に変化する(核のまわりに孤立して局在し，外側で急激に減衰する)．したがって $s$ が大きく，LDAの誤差も大きい．11.5.1節で見たようにLDAは原子の $E_{xc}$ を浅く見積もるので，原子側のエネルギーが実際より高くなる．その結果，分子との差である結合エネルギーが過大になる．GGAは $F_x(s)$ を通じて $s$ の大きい原子側をより深く安定化するので，この不均衡が緩和される．
同じ理由で，GGAは以下の系統的な改善をもたらす．
- 結合長は LDA が $1$〜$2\%$ 短く，PBE は実験値に近いか $1\%$ 程度長い．
- 振動数は LDA が高すぎ(結合が強すぎる)，PBE が改善する．
- 水素結合(水の二量体など)は LDA が大幅に過大，PBE がほぼ正しい．
11.5.3 固体:格子定数と体積弾性率
固体では傾向が反転する場面が出てくる．
| 物質 | $a_0^{\mathrm{LDA}}$ | $a_0^{\mathrm{PBE}}$ | $a_0^{\mathrm{exp}}$ | $B_0^{\mathrm{LDA}}$ | $B_0^{\mathrm{PBE}}$ | $B_0^{\mathrm{exp}}$ |
|---|---|---|---|---|---|---|
| C(ダイヤモンド) | 3.535 | 3.574 | 3.555 | 466 | 431 | 443 |
| Si | 5.410 | 5.479 | 5.422 | 96 | 89 | 99 |
| Al | 3.985 | 4.041 | 4.020 | 84 | 77 | 79 |
| Cu | 3.523 | 3.635 | 3.596 | 190 | 141 | 142 |
| MgO | 4.169 | 4.256 | 4.188 | 172 | 149 | 165 |
| NaCl | 5.474 | 5.699 | 5.565 | 32 | 23 | 27 |
物理的意味:「LDAは縮み，GGAは伸びる」
LDAが過結合するのは分子でも固体でも同じで，その結果として結合が強すぎ，平衡体積が小さくなる．これがLDAの過結合(overbinding)である．
PBEはこれを補正するが，しばしば行き過ぎる．原因は11.4.5節で見た $\mu=0.2195$ という選択にある．この値は「原子・分子のように密度が急変する系」に合わせたもので，固体のバルク(密度が緩やかに変化する領域)では交換増強が強すぎ，密度の不均一性を過度に好む．結果として原子を引き離す方向に働き，格子定数が伸びる．$\mu=10/81$ に戻したPBEsolが固体で優れているのはこのためである．
体積弾性率 $B_0=V\,\dd^2E/\dd V^2$ は格子定数と逆相関する．格子定数が小さすぎれば $B_0$ は大きすぎ，逆も然り．したがって表11.7で $a_0$ と $B_0$ の誤差の符号が逆になっているのは当然である．この相関を知っていれば，「格子定数がずれているのに弾性定数だけ合っている」といった結果を疑うことができる．
凝集エネルギー(固体と孤立原子のエネルギー差)も同様で，LDAは $15$〜$25\%$ 過大，PBEは $5\%$ 前後の誤差に収まる．たとえばSiでは LDA $5.28$ eV，PBE $4.55$ eV に対し実験値 $4.63$ eV，Cuでは LDA $4.55$ eV，PBE $3.51$ eV に対し実験値 $3.49$ eV である．
11.5.4 決定的な例:bcc鉄の基底状態
LDAとGGAの違いが「精度の程度問題」ではなく「定性的に正しいか誤りか」を分ける有名な例が，鉄の基底状態である．
例:鉄の基底状態を予測する
実験事実:常温常圧の鉄は強磁性の体心立方(FM-bcc)構造をとる．原子あたりの磁気モーメントは $2.2\mu_B$ である．
計算では，格子定数(あるいはWigner–Seitz半径 $R_{\mathrm{WS}}$)を変えながら，いくつかの構造・磁気状態について全エネルギー曲線を描き，最も低い極小を探す．候補は非磁性fcc(NM-fcc)，非磁性hcp(NM-hcp)，強磁性bcc(FM-bcc)である．
- LDAの答え:非磁性の稠密構造(fcc/hcp)が最も安定．FM-bccはこれより高い．実験と定性的に食い違う．
- GGA(PBE)の答え:FM-bccが最も安定で，平衡格子定数も磁気モーメントも実験とよく合う．正しい．
なぜこの差が出るのか．鍵はスピン分極のコストと利得のバランスである．
- 強磁性状態では $3d$ 電子のスピンが揃うため，同スピン電子が増え，交換エネルギーの利得が大きい．一方，スピンを揃えるとフェルミ準位が上がり運動エネルギーが損をする(ストーナー機構)．両者の差は $0.01$ Ha($0.3$ eV)程度の微妙な競合である．
- bcc構造はfcc/hcpより充填率が低く，原子まわりに「隙間」が多い．その分，$d$ 電子は局在的になり，密度の不均一性($s$ の値)が大きい．
- GGAは $F_x(s)$ を通じて $s$ の大きい局在的な状態を安定化する．すなわちGGAは局在化とスピン分極の両方を，LDAより有利に評価する．この $0.01$ Ha オーダーの補正が，上の微妙な競合をひっくり返す．
この例は，汎関数の選択が単なる「精度」ではなく「予測される物理そのもの」を変えうることを示している．磁性金属や磁気秩序を扱うときは，GGA以上を使うのが最低条件である．
11.5.5 GGAが常に良いわけではない
「GGAはLDAの上位互換」ではない．GGAが劣る場面を知っておく必要がある．
- 格子定数の過大評価．表11.7のとおり，重い元素や貴金属(Ag，Au)ではPBEの過大評価が $1.5$〜$2\%$ に達し，LDAの過小評価より悪いことがある．Agでは LDA $4.01$ Å，PBE $4.15$ Å に対し実験値 $4.07$ Å である．
- 表面エネルギー．金属表面のエネルギーは，LDAが実験値に近く，PBEは $10$〜$30\%$ 過小評価する．ジェリウム表面という厳密解の分かっている系でも，LDAの方が正確である．これは表面の密度勾配領域で $F_x$ が働きすぎるためで，PBEsolが開発された直接の動機の一つである．
- 強く圧縮された系．高圧下では密度が上がって $s$ が小さくなり，LDAとGGAの差は縮む．ただしPBEの相関が高密度極限で有限にとどまるという性質(条件(c))が効き，高圧相図ではGGAの方が信頼できる．
- 両者に共通する破綻．van der Waals力(11.6節)，バンドギャップ($40$〜$50\%$ 過小)，強相関系の $d$/$f$ 電子，反応障壁($30\%$ 過小)は，GGAでも直らない．これらは表11.2の下3行の×に対応する本質的な欠陥である．
注意:計算値を読むときのチェックリスト
- 汎関数名が書かれているか．「DFT計算によると」だけでは情報がない．LDA/PBE/PBEsol/SCAN/HSEでは結果が有意に違う．
- 誤差の向きが分かっているか．PBEの格子定数は大きめに出ると分かっていれば，実験値との $1\%$ のずれは想定内と判断できる．逆にPBEが実験より小さい格子定数を出したら，計算のどこかを疑うべきである．
- 比較しているのは差か絶対値か．全エネルギーの絶対値には $E_{xc}$ の誤差がまるごと乗るが，同種の系どうしの差では大きく相殺する．異なる汎関数で計算した全エネルギーを直接引き算してはならない．
- 実験値の補正．格子定数の実験値にはゼロ点振動と熱膨張が含まれる．$0$ K・振動なしの計算値と比べるには補正が必要で，典型的には $0.5\%$ 程度である．
この節で得られたもの．LDAは交換を浅く相関を深く見積もり，その相殺のうえに実用精度が立っている．GGA(PBE)は両方の誤差を1桁縮め，分子の過結合と磁性金属の記述を劇的に改善するが，固体の格子定数と表面エネルギーではやり過ぎる．次節では，これらの汎関数に共通して残る本質的欠陥を直すための，さらに上の段を見る．
11.6 さらに先へ:ヤコブの梯子
Perdewは，汎関数の発展を「ハートリー世界(地上，$E_{xc}=0$)から化学精度の天国へ向かう梯子」になぞらえた．各段は新しい材料(ingredient)を1つ加えることで定義される．段を上がるほど精度は上がるが計算コストも増え，しかも「上の段が必ず良い」保証はない——梯子は上るものであって，飛び越えるものではない．
11.6.1 出発点:自己相互作用誤差とバンドギャップ問題
第3段までの汎関数(LDA，GGA，meta-GGA)に共通して残る欠陥を，第10章の結果と接続して整理しておく．基本ギャップ $E_g=I-A$ とKSギャップの関係は
$$ \begin{equation} E_g=E_g^{\mathrm{KS}}+\Delta_{xc}, \qquad \Delta_{xc}=\lim_{\eta\to0^+}\left[ \frac{\delta E_{xc}}{\delta\rho(\rr)}\bigg|_{N+\eta}-\frac{\delta E_{xc}}{\delta\rho(\rr)}\bigg|_{N-\eta}\right] \label{eq:11-gap-recall} \end{equation} $$であった．$E_{xc}$ が $\rho$，$\nabla\rho$，$\tau$ のなめらかな関数の積分で書けている限り，$v_{xc}$ も密度のなめらかな関数であり，電子数が整数を横切っても跳ばない．したがって $\Delta_{xc}\equiv0$ となり，計算されるのは $E_g^{\mathrm{KS}}$ だけである．これに自己相互作用誤差(11.3.5節)による占有準位の押し上げが重なり，ギャップは $40$〜$50\%$ 過小評価される．
| 物質 | PBE | HF | HSE06 | 実験 |
|---|---|---|---|---|
| Si | 0.62 | 6.4 | 1.15 | 1.17 |
| GaAs | 0.49 | 7.0 | 1.21 | 1.52 |
| SiC | 1.35 | 7.9 | 2.36 | 2.42 |
| C(ダイヤ) | 4.15 | 12.3 | 5.36 | 5.48 |
| MgO | 4.75 | 16.2 | 6.50 | 7.83 |
物理的意味:HFとGGAは反対方向に誤る
ハートリー・フォック法では，占有軌道の電子は他の $N-1$ 個の電子だけを感じる(自己相互作用が交換で厳密に打ち消される)．一方，非占有軌道の電子は $N$ 個すべてを感じる．この非対称性により占有・非占有の間に大きな隔たりが生まれ，ギャップが過大になる．しかもHFには相関による遮蔽が全く入っていないため，隔たりを縮める効果もない．
GGAでは逆に，自己相互作用が残るため占有軌道の電子が「自分を含む $N$ 個」を感じてしまい，占有準位が浅くなってギャップが縮む．
2つの誤差が逆向きなら，混ぜれば中間になる．これがハイブリッド汎関数の直観的な動機である．だがそれだけでは混合比を決められない．次節で断熱接続から根拠を与える．
11.6.2 ハイブリッド汎関数:混合比の根拠
11.2.3節で導いた断熱接続公式 \eqref{eq:11-ac-key} をもう一度書く:
$$ E_{xc}=\int_0^1 \dd\lambda\,W_\lambda, \qquad W_0=E_x^{\mathrm{exact}}, \qquad W_1=V_{ee}-E_{\mathrm{H}}. $$ここで $W_0$ が厳密交換(KS軌道で組んだ交換積分)であることが決定的である．積分の端点の一方は厳密に分かっているのだから，これを使わない手はない．
導出:混合比 $a=1/4$ の根拠(Perdew–Ernzerhof–Burke)
ステップ1:台形則(Beckeのhalf-and-half)．最も素朴には，$W_\lambda$ が $\lambda$ の1次関数だと仮定して台形則を使う:
$$ \begin{equation} E_{xc}\simeq\frac{W_0+W_1}{2}=\frac{1}{2}E_x^{\mathrm{exact}}+\frac{1}{2}W_1. \label{eq:11-halfhalf-1} \end{equation} $$$W_1$ は分からないので，近似汎関数の値 $E_{xc}^{\mathrm{DFA}}$(DFA = density functional approximation，LSDAやPBE)で置き換えると
$$ \begin{equation} E_{xc}^{\text{half\&half}}=\frac{1}{2}E_x^{\mathrm{exact}}+\frac{1}{2}E_{xc}^{\mathrm{LSDA}}. \label{eq:11-halfhalf-2} \end{equation} $$これがBeckeが1993年に提案した最初のハイブリッド汎関数である．分子の解離エネルギーが劇的に改善した．
ステップ2:$W_\lambda$ の曲率を考慮する．図11.3で見たとおり，$W_\lambda$ は直線ではなく $\lambda$ について下に凸である．$\lambda$ が小さいうちは急に深くなり，大きくなると変化が鈍る．したがって台形則は $W_0$ の重みを過大に見積もる．曲率を取り込むため，モデル関数
$$ \begin{equation} W_\lambda\simeq W_1+(W_0-W_1)(1-\lambda)^{n-1} \label{eq:11-Wlambda-model} \end{equation} $$を仮定する($n\ge2$ の整数)．$\lambda=0$ で $W_0$，$\lambda=1$ で $W_1$ を正しく与え，$n$ が大きいほど $W_0$ から $W_1$ への移行が速い(=曲率が大きい)．積分すると
\begin{align} \int_0^1 W_\lambda\,\dd\lambda &=W_1+(W_0-W_1)\int_0^1(1-\lambda)^{n-1}\dd\lambda \label{eq:11-peb-1}\\ &=W_1+(W_0-W_1)\left[-\frac{(1-\lambda)^n}{n}\right]_0^1 \label{eq:11-peb-2}\\ &=W_1+\frac{1}{n}(W_0-W_1). \label{eq:11-peb-3} \end{align}\eqref{eq:11-peb-2} では $\int_0^1(1-\lambda)^{n-1}\dd\lambda$ を $u=1-\lambda$ の置換($\dd u=-\dd\lambda$，積分区間が $1\to0$ に反転)で $\int_0^1u^{n-1}\dd u=1/n$ と計算した．
ステップ3:近似汎関数で $W_1$ を消す．$W_1\simeq E_{xc}^{\mathrm{DFA}}$，$W_0=E_x^{\mathrm{exact}}$ とし，さらに近似汎関数側でも同じ関係 $E_{xc}^{\mathrm{DFA}}\simeq E_x^{\mathrm{DFA}}+\frac{1}{n}(\cdots)$ が成り立つとすると，整理して
$$ \begin{equation} E_{xc}^{\mathrm{hyb}}=E_{xc}^{\mathrm{DFA}} +\frac{1}{n}\Bigl(E_x^{\mathrm{exact}}-E_x^{\mathrm{DFA}}\Bigr) \label{eq:11-hybrid-general} \end{equation} $$を得る．すなわち，近似汎関数の交換の一部を厳密交換で置き換える形になる．
ステップ4:$n$ を決める．PerdewとErnzerhofとBurkeは，典型的な分子の原子化エネルギーに対する $W_\lambda$ の曲がり方を4次の摂動論的解析から見積もり，$n=4$ が妥当であると結論した．したがって混合比は
$$ \begin{equation} a=\frac{1}{n}=\frac{1}{4}=0.25. \label{eq:11-a-quarter} \end{equation} $$ ∎定義:代表的なハイブリッド汎関数
PBE0(PBE1PBE)——非経験的ハイブリッド．上の $n=4$ をPBEに適用したもの:
$$ \begin{equation} E_{xc}^{\mathrm{PBE0}}=\frac{1}{4}E_x^{\mathrm{exact}}+\frac{3}{4}E_x^{\mathrm{PBE}}+E_c^{\mathrm{PBE}} \label{eq:11-pbe0} \end{equation} $$B3LYP——3パラメータ経験的ハイブリッド．量子化学で最も広く使われる:
$$ \begin{equation} E_{xc}^{\mathrm{B3LYP}}=E_{xc}^{\mathrm{LSDA}} +a_0\bigl(E_x^{\mathrm{exact}}-E_x^{\mathrm{LSDA}}\bigr) +a_x\,\Delta E_x^{\mathrm{B88}} +a_c\,\Delta E_c^{\mathrm{LYP}} \label{eq:11-b3lyp} \end{equation} $$$a_0=0.20$，$a_x=0.72$，$a_c=0.81$ は分子の実験データ(生成熱・イオン化エネルギー・プロトン親和力)へのフィットで決められた．$a_0=0.20$ が上の $1/4$ に近いことは示唆的である．
HSE(Heyd–Scuseria–Ernzerhof)——遮蔽ハイブリッド．クーロン核を短距離と長距離に分割する:
$$ \begin{equation} \frac{1}{r}=\underbrace{\frac{\mathrm{erfc}(\omega r)}{r}}_{\text{短距離(SR)}} +\underbrace{\frac{\mathrm{erf}(\omega r)}{r}}_{\text{長距離(LR)}} \label{eq:11-range-sep} \end{equation} $$($\mathrm{erf}+\mathrm{erfc}=1$ なので恒等式である．$\mathrm{erfc}(\omega r)/r$ は $r\gg1/\omega$ でガウス的に急減衰し，$\mathrm{erf}(\omega r)/r$ は $r\to0$ で有限値 $2\omega/\sqrt{\pi}$ に留まる．)厳密交換は短距離部分にだけ混ぜる:
$$ \begin{equation} E_{xc}^{\mathrm{HSE}}=\frac{1}{4}E_x^{\mathrm{exact,SR}}(\omega) +\frac{3}{4}E_x^{\mathrm{PBE,SR}}(\omega)+E_x^{\mathrm{PBE,LR}}(\omega)+E_c^{\mathrm{PBE}} \label{eq:11-hse} \end{equation} $$HSE06では $\omega=0.11\ a_0^{-1}$(遮蔽長 $\simeq9\ a_0$)．
物理的意味:なぜ長距離の厳密交換を切るのか
金属では，フェルミ面近傍の状態に対する厳密交換の長距離部分が，状態密度をフェルミ準位でゼロにしてしまう(HF法の悪名高い「対数発散」)．実際の金属では，遍歴電子による遮蔽がこの長距離交換を打ち消すため，そのような異常は起こらない．$\mathrm{erfc}$ で長距離を切り落とすことは，この遮蔽を粗く模擬することに相当する．
実務上の利点も大きい．平面波基底(第14章)で厳密交換を評価するとき，長距離部分は逆空間で $1/q^2$ の特異性を持ち，収束が非常に遅い．$\mathrm{erfc}$ で切ると被積分関数が $\exp(-q^2/4\omega^2)$ で減衰し，計算量が桁で減る．固体の第一原理計算でHSEが標準になっているのは，この2つの理由による．表11.8のとおり，ギャップの再現は $GW$ 近似(第18章)に匹敵する．
ただしハイブリッド汎関数のコストはGGAの $10$〜$100$ 倍であり，計算できる系のサイズは大きく制限される．また混合比 $a$ は本来系依存(遮蔽の強い金属では小さく，絶縁体では大きくすべき)であり，$a=1/4$ を万能に使うことの理論的正当性は限定的である．
11.6.3 meta-GGA:運動エネルギー密度を材料にする
第3段は，KS軌道から作られる運動エネルギー密度
$$ \begin{equation} \tau(\rr)=\frac{1}{2}\sum_{i}^{\mathrm{occ}}\abs{\nabla\phi_i(\rr)}^2 \label{eq:11-tau-def} \end{equation} $$を新しい材料に加える．$\tau$ は軌道から作られるが，$\rr$ の1点での値しか使わないので「準局所」であり，コストはGGAとほとんど変わらない．これが $\tau$ を使う最大の実利である．
$\tau$ が有用なのは，密度と勾配だけでは区別できない状況を見分けられるからである．指標として
$$ \begin{equation} \alpha(\rr)\equiv\frac{\tau-\tau_W}{\tau^{\mathrm{unif}}}, \qquad \tau_W=\frac{\abs{\nabla\rho}^2}{8\rho}, \qquad \tau^{\mathrm{unif}}=\frac{3}{10}(3\pi^2)^{2/3}\rho^{5/3} \label{eq:11-alpha-def} \end{equation} $$を使う．$\tau_W$ はvon Weizsäcker運動エネルギー密度である．
導出:$\alpha$ の3つの目印
(i) 1軌道領域では $\alpha=0$．その点の密度が1本の実軌道 $\phi=\sqrt{\rho}$ だけで作られているとする(1電子系，共有結合の結合中心，原子の最外殻の裾など)．$\rho=\phi^2$ すなわち $\phi=\sqrt\rho$ なので，連鎖律より
$$ \nabla\phi=\frac{\nabla\rho}{2\sqrt\rho}, \qquad \abs{\nabla\phi}^2=\frac{\abs{\nabla\rho}^2}{4\rho}. $$したがって
$$ \tau=\frac{1}{2}\abs{\nabla\phi}^2=\frac{\abs{\nabla\rho}^2}{8\rho}=\tau_W, $$すなわち $\alpha=(\tau_W-\tau_W)/\tau^{\mathrm{unif}}=0$ となる．
(ii) 一様電子ガスでは $\alpha=1$．$\nabla\rho=0$ なので $\tau_W=0$，また $\tau$ はフェルミ球の運動エネルギー密度そのもの $\tau^{\mathrm{unif}}$ である(第3章)．よって $\alpha=1$．
(iii) 密度の重なりが弱い領域では $\alpha\gg1$．2つの閉殻フラグメントの隙間では，それぞれのフラグメントから来た軌道が浅く重なる．$\rho$ も $\nabla\rho$ も小さいので $\tau_W$ は小さいが，$\tau$ は各軌道の勾配の2乗の和なので相対的に大きく残る．$\tau^{\mathrm{unif}}\propto\rho^{5/3}$ が急速に小さくなるため $\alpha$ は大きくなる．
また $\tau\ge\tau_W$(1軌道のとき等号)が一般に成り立つので，$\alpha\ge0$ である．∎
この3つの目印により，meta-GGAは「いま自分がいるのは共有結合か，金属的な領域か，弱い相互作用の隙間か」を判定し，それぞれに適した増強因子を使い分けられる．GGAには $s$ しか手がかりがなく，たとえば「共有結合の中心」と「弱い重なり領域」がどちらも小さい $s$ を持つため区別できなかった．
定義:SCAN汎関数
SCAN(Strongly Constrained and Appropriately Normed，Sun–Ruzsinszky–Perdew 2015)は，$\rho$，$\nabla\rho$，$\tau$ を材料とするmeta-GGAが満たしうる既知の厳密条件17個をすべて満たすように構成された非経験的汎関数である．さらに，厳密条件だけでは決まらない自由度を，「厳密に近い値が分かっている規範系(norm)」——希ガス原子，ジェリウム表面，圧縮された希ガス二量体など，いずれも結合が関与しない系——に合わせて固定している．
形式的には，$\alpha$ を補間パラメータとして「1軌道極限($\alpha=0$)に適した汎関数」と「緩やかに変化する極限($\alpha=1$)に適した汎関数」をなめらかにつなぐ:
$$ E_x^{\mathrm{SCAN}}=\int \dd^3r\,\rho\,\varepsilon_x^{\mathrm{unif}} \Bigl[h_x^{1}(s)+f(\alpha)\bigl\{h_x^{0}-h_x^{1}(s)\bigr\}\Bigr]g_x(s) $$($f(\alpha)$ は $f(0)=1$，$f(1)=0$ の補間関数)．SCANは分子の結合エネルギー，固体の格子定数，氷の多形の相対安定性などを同時に高精度で与え，meta-GGAの実用性を確立した．ただし数値的に敏感で，細かい実空間グリッドを要する．この点を改良したのが r$^2$SCAN である．
11.6.4 van der Waals力:なぜ準局所汎関数では出ないのか
希ガス二量体，層状物質(グラファイト，遷移金属ダイカルコゲナイド)，分子結晶，物理吸着——これらを支配するのはvan der Waals(vdW，ロンドン分散)力である．LDAもGGAもこれを記述できない．理由を厳密に示す．
定理:準局所汎関数は密度の重ならない2体間に相互作用を与えない
2つのフラグメントAとBの密度 $\rho_A$，$\rho_B$ の台(値がゼロでない領域)$\Omega_A$，$\Omega_B$ が交わらない($\Omega_A\cap\Omega_B=\emptyset$)とする．このとき任意の準局所汎関数 $E_{xc}[\rho]=\int f(\rho,\nabla\rho,\tau)\dd^3r$ について
$$ \begin{equation} E_{xc}[\rho_A+\rho_B]=E_{xc}[\rho_A]+E_{xc}[\rho_B] \label{eq:11-vdw-additive} \end{equation} $$が厳密に成り立つ．すなわち交換相関からの相互作用エネルギーは恒等的にゼロである．
証明．$\rho=\rho_A+\rho_B$ とする．空間を3つに分ける．
- $\rr\in\Omega_A$:$\rho_B(\rr)=0$ かつ($\Omega_B$ から離れているので)$\nabla\rho_B(\rr)=0$，$\tau_B(\rr)=0$．よって $f(\rho,\nabla\rho,\tau)=f(\rho_A,\nabla\rho_A,\tau_A)$．
- $\rr\in\Omega_B$:同様に $f=f(\rho_B,\nabla\rho_B,\tau_B)$．
- それ以外:$\rho=0$，$\nabla\rho=0$，$\tau=0$ なので $f(0,0,0)=0$(交換相関エネルギー密度は電子がなければゼロ)．
したがって
$$ \int_{\text{全空間}}f\,\dd^3r=\int_{\Omega_A}f(\rho_A,\dots)\dd^3r+\int_{\Omega_B}f(\rho_B,\dots)\dd^3r =E_{xc}[\rho_A]+E_{xc}[\rho_B]. $$ ∎他の項も同様である．$T_s$ は，互いに重ならない領域に局在した軌道が自動的に直交するので加算的である．$E_{\mathrm{H}}+\int v\rho+E_{nn}$ は古典静電相互作用を与えるが，中性で球対称なフラグメント(希ガス原子など)ではこれもゼロである．結局，準局所DFTは中性フラグメント間に $R\to\infty$ で厳密にゼロの相互作用しか与えない．真のvdW引力 $-C_6/R^6$ はどこにも現れない．
物理的意味:vdW力は「遠く離れた電子密度ゆらぎの相関」である
vdW引力の起源は，フラグメントAの中で瞬間的に生じた電気双極子ゆらぎが，Bに誘起双極子を作り，両者が引き合うことである．これは$\rr\in\Omega_A$ の電子と $\rr'\in\Omega_B$ の電子の間の相関であり，11.2節の言葉では「相関ホールの尾が別のフラグメントまで伸びている」状況である．$\rr$ と $\rr'$ を結ぶ非局所的な情報が本質なので，点 $\rr$ の局所量だけを見る汎関数には原理的に表現できない．
なお，実際の計算で希ガス二量体を扱うとLDAはしばしば「それらしい」結合を与える．これは密度の裾がわずかに重なる領域でLDAの交換が誤って引力を作っているためで，$R^{-6}$ ではなく指数関数的に減衰する偽物である．PBEはほとんど結合を与えない．どちらも正しくない．
定義:vdWを取り込む2つの処方
(a) 非局所相関汎関数(vdW-DF)——Dionらが2004年に提案．相関を局所部分と非局所部分に分ける:
$$ \begin{equation} E_c=E_c^{\mathrm{LDA}}[\rho]+E_c^{\mathrm{nl}}[\rho], \qquad E_c^{\mathrm{nl}}[\rho]=\frac{1}{2}\iint \dd^3r\,\dd^3r'\; \rho(\rr)\,\phi\bigl(\rr,\rr'\bigr)\,\rho(\rr') \label{eq:11-vdwdf} \end{equation} $$核 $\phi$ は $\abs{\rr-\rr'}$ と，両点における局所量 $q_0(\rr)=q_0[\rho(\rr),\abs{\nabla\rho(\rr)}]$ に依存する．$\phi$ の形は，断熱接続・ゆらぎ散逸定理から出発し，誘電応答をプラズモン極モデルで近似して導かれる．$\abs{\rr-\rr'}\to\infty$ で $\phi\propto-1/\abs{\rr-\rr'}^6$ となるように作られているので，二重積分は正しく $-C_6/R^6$ を与える．密度だけから $C_6$ が出るのが最大の長所である．改良版として vdW-DF2，rev-vdW-DF2，optB88-vdW などがある．
(b) 経験的分散補正(DFT-D)——Grimmeらの処方．全エネルギーに原子対の項を足す:
$$ \begin{equation} E_{\mathrm{disp}}=-\sum_{A\lt B}\left[ s_6\frac{C_6^{AB}}{R_{AB}^6}f_{\mathrm{damp}}^{(6)}(R_{AB}) +s_8\frac{C_8^{AB}}{R_{AB}^8}f_{\mathrm{damp}}^{(8)}(R_{AB})\right] \label{eq:11-dftd} \end{equation} $$$C_6^{AB}$ は原子種の組ごとに与えられた分散係数，$f_{\mathrm{damp}}$ は短距離で発散を抑える減衰関数(短距離ではGGAが既に相関を含んでいるので二重計上を避ける)．DFT-D3では $C_6$ を各原子の配位数の関数にすることで化学環境依存性を取り込んでいる．計算コストはほぼゼロで，既存の計算に後付けできるのが利点である．
11.6.5 DFT+U:局在電子へのハバード補正
遷移金属酸化物(NiO，CoO，MnO)は実験的には反強磁性絶縁体だが，LDA/GGAでは金属または極端に小さいギャップの半導体と予測される．原因は $3d$ 電子の強い局在性と，それに伴う自己相互作用誤差である．$d$ 軌道は空間的に狭く，自己ハートリーが大きいため，SIEも大きい．その結果，本来なら整数占有(占有された $d$ 軌道と空の $d$ 軌道が分離)すべきところで，分数占有が不当に安定化されてしまう．
処方は，局在軌道に対してだけハバード型の相互作用を陽に加えることである．
導出:DFT+Uの簡約形(Dudarev形)と二重数え補正
ステップ1:全体の構造．
$$ \begin{equation} E^{\mathrm{DFT}+U}=E^{\mathrm{DFT}}[\rho]+E^{\mathrm{Hub}}\bigl[\{n_{m\sigma}\}\bigr]-E^{\mathrm{dc}}\bigl[\{n_{m\sigma}\}\bigr] \label{eq:11-dftu-structure} \end{equation} $$$n_{m\sigma}$ は着目原子の局在軌道($m$ は $d$ 軌道なら $m=1,\dots,5$)の占有数で，KS波動関数を局在軌道へ射影して求める．$E^{\mathrm{Hub}}$ が新たに加えるハバード項，$E^{\mathrm{dc}}$ が二重数え補正である．二重数え補正が必要なのは，$E^{\mathrm{DFT}}$ が既に $d$ 電子間の平均的なクーロン反発を(平均場的に)含んでいるからで，これを引かないと同じ相互作用を2回数えることになる．
ステップ2:ハバード項(球平均近似)．簡単のため1つのスピンチャネルと $J=0$ の場合を書く．異なる軌道にいる電子どうしが $U$ で反発するとして
$$ \begin{equation} E^{\mathrm{Hub}}=\frac{U}{2}\sum_{m\ne m'}n_m n_{m'} =\frac{U}{2}\left[\Bigl(\sum_m n_m\Bigr)^2-\sum_m n_m^2\right] =\frac{U}{2}\Bigl[N^2-\sum_m n_m^2\Bigr] \label{eq:11-dftu-hub} \end{equation} $$($N=\sum_m n_m$ は全占有数)．2番目の等号では，$\sum_{m,m'}n_mn_{m'}=N^2$ から対角項 $\sum_m n_m^2$ を引いた．
ステップ3:二重数え(fully localized limit)．DFTが含んでいる平均場的な反発は，「$N$ 個の電子が互いに $U$ で反発する古典的な値」
$$ \begin{equation} E^{\mathrm{dc}}=\frac{U}{2}N(N-1) \label{eq:11-dftu-dc} \end{equation} $$と見積もる($N$ 個から2個選ぶ組合せの数 $N(N-1)/2$ に $U$ を掛けたもの)．
ステップ4:差をとる．
\begin{align} E_U&=E^{\mathrm{Hub}}-E^{\mathrm{dc}} =\frac{U}{2}\Bigl[N^2-\sum_m n_m^2\Bigr]-\frac{U}{2}\bigl[N^2-N\bigr] \label{eq:11-dftu-diff-1}\\ &=\frac{U}{2}\Bigl[N-\sum_m n_m^2\Bigr] \label{eq:11-dftu-diff-2}\\ &=\frac{U}{2}\sum_m\bigl(n_m-n_m^2\bigr) =\frac{U}{2}\sum_m n_m\bigl(1-n_m\bigr). \label{eq:11-dftu-diff-3} \end{align}\eqref{eq:11-dftu-diff-2} で $N^2$ が相殺した．\eqref{eq:11-dftu-diff-3} では $N=\sum_m n_m$ を和の中に入れた．交換 $J$ も含めた場合，$U\to U_{\mathrm{eff}}=U-J$ と置き換えればよい．スピンを戻して行列表記すると
$$ \begin{equation} E_U=\sum_{I}\frac{U_{\mathrm{eff}}^{I}}{2}\sum_{\sigma} \mathrm{Tr}\Bigl[\bm{n}^{I\sigma}\bigl(\bm{1}-\bm{n}^{I\sigma}\bigr)\Bigr] \label{eq:11-dftu-key} \end{equation} $$($I$ は補正を加える原子，$\bm{n}^{I\sigma}$ は占有数行列)．∎
物理的意味:$E_U$ は分数占有に罰金を課す
式 \eqref{eq:11-dftu-diff-3} の $n(1-n)$ は，$n=0$ と $n=1$ でゼロ，$n=1/2$ で最大 $U_{\mathrm{eff}}/8$ をとる下に凸でない(上に凸な)放物線である．したがって $E_U\ge0$ であり，整数占有では罰金がゼロ，分数占有では正の罰金となる．これはSIEが引き起こす「分数占有の不当な安定化」をちょうど打ち消す方向に働く．
準位のシフト．ポテンシャルは占有数微分で与えられる:
$$ v_U^{m}=\frac{\partial E_U}{\partial n_m}=\frac{U_{\mathrm{eff}}}{2}\bigl(1-2n_m\bigr) $$占有軌道($n_m=1$)は $-U_{\mathrm{eff}}/2$ だけ下がり，空軌道($n_m=0$)は $+U_{\mathrm{eff}}/2$ だけ上がる．両者の間に $U_{\mathrm{eff}}$ のギャップが開く．これがNiOやCoOのギャップを回復する機構である．
軌道分極の駆動．3重縮退した軌道に2電子がいる場合を考える．平均占有 $(2/3,2/3,2/3)$ では
$$ E_U=3\times\frac{U_{\mathrm{eff}}}{2}\cdot\frac{2}{3}\cdot\frac{1}{3}=\frac{U_{\mathrm{eff}}}{3}, $$分極占有 $(1,1,0)$ では $E_U=0$ である．差 $U_{\mathrm{eff}}/3$ だけ分極状態が有利になる．DFT+Uが軌道秩序やヤーン・テラー歪みを再現できるのはこのためである．
注意:$U$ をどう決めるか，そして何が犠牲になるか
$U$ の決め方．(i) 実験のギャップや磁気モーメントに合わせる(経験的，最も一般的)．(ii) 拘束付きDFT:局在軌道の占有数を強制的に変えたときのエネルギー応答から $U=\partial^2E/\partial n^2$ を求める(Cococcioni–de Gironcoliの線形応答法)．(iii) 制限付きRPA(cRPA):遮蔽された相互作用を第一原理的に計算する．3d遷移金属酸化物では $U_{\mathrm{eff}}=3$〜$8$ eV が典型的である．
犠牲になるもの．DFT+Uは補正を加える原子と軌道を人間が指定する必要があり，その意味で第一原理性を一部失う．また，局在軌道の定義(射影の仕方)に結果が依存する．金属的な遍歴性と局在性が競合する系(動的相関が本質的な系)では，$U$ を一つの静的な数で表す近似そのものが不十分で，動的平均場理論(DMFT)のような枠組みが必要になる．
11.6.6 第5段:非占有軌道を使う汎関数
梯子の最上段では，非占有KS軌道とその固有値も材料に加える．代表はRPA(乱雑位相近似)相関エネルギーで，断熱接続とゆらぎ散逸定理から
$$ \begin{equation} E_c^{\mathrm{RPA}}=\frac{1}{2\pi}\int_0^\infty \dd\omega\; \mathrm{Tr}\Bigl[\ln\bigl(1-\chi_0(i\omega)v\bigr)+\chi_0(i\omega)v\Bigr] \label{eq:11-rpa} \end{equation} $$と書ける($\chi_0$ は非相互作用応答関数で，占有・非占有軌道の両方を使って作られる，第8章)．RPAはvdW力を第一原理的に含み，金属も絶縁体も同じ枠組みで扱えるという美点を持つ．一方，コストはGGAの $10^3$ 倍以上に達し，また短距離相関の記述が不十分で結合エネルギーは系統的に過小になる．
量子化学では二重ハイブリッド汎関数(B2PLYPなど)が使われる．ハイブリッド汎関数に，非占有軌道を使う2次のMøller–Plesset摂動論的相関を一部混ぜたもので，化学精度($1$ kcal/mol)に手が届く．
この節で得られたもの．ヤコブの梯子の5段を，それぞれ何を材料に加え，何を直すのかという観点から概観した．第4段以上でのみ自己相互作用の消去と $\Delta_{xc}\ne0$ が可能になること，vdWには非局所相関か経験的補正が必要なこと，局在 $d$/$f$ 電子にはDFT+Uが効くことを見た．
11.7 発展:KS法の変分的性質とハリス汎関数
本節は交換相関汎関数そのものの話から少し離れ，KS法が持つ変分的な性質を利用した実用技術を扱う．中心となるのは「自己無撞着に解かなくても，1回の対角化で高精度な全エネルギーが得られる」というハリス汎関数である．第14章の実装，第16章の力の計算，そして磁気異方性のような微小エネルギー差の評価に直結する．
11.7.1 動機:SCFは高価である
第10章で見たとおり，KS方程式は不動点問題 $\rho_{\mathrm{out}}=\mathcal{K}[\rho_{\mathrm{in}}]$ であり，収束まで数十回の対角化を要する．もし1回の対角化で十分な精度の全エネルギーが得られるなら，計算量は1桁減る．実際，原子密度を重ね合わせた $\rho_{\mathrm{in}}$ から出発して1回だけ対角化するという処方は古くから使われてきた．問題は「その値はどれだけ信用できるのか」である．
KS全エネルギー汎関数を非自己無撞着な密度 $\rho_{\mathrm{in}}$ で評価すること自体は，変分原理から意味がある:
$$ \begin{equation} E_{\mathrm{KS}}[\rho_{\mathrm{in}}]=T_s[\rho_{\mathrm{in}}]+\int v\rho_{\mathrm{in}} +E_{\mathrm{H}}[\rho_{\mathrm{in}}]+E_{xc}[\rho_{\mathrm{in}}]\ \ge\ E_0 \label{eq:11-eks-variational} \end{equation} $$で，誤差は $\delta\rho=\rho_{\mathrm{in}}-\rho_0$ の2次である($\rho_0$ は自己無撞着密度)．しかしこの式には $T_s[\rho_{\mathrm{in}}]$ が現れ，$T_s$ を密度から直接計算する術がない(第10章)．実際に計算できるのは「$\rho_{\mathrm{in}}$ から作った $v_{\mathrm{eff}}$ で対角化して得た固有値と軌道」である．この使える情報だけで作った汎関数がハリス汎関数である．
11.7.2 ハリス汎関数の定義
定義:ハリス汎関数(Harris functional)
入力密度 $\rho_{\mathrm{in}}$ から有効ポテンシャル
$$ \begin{equation} v_{\mathrm{eff}}[\rho_{\mathrm{in}}](\rr)=v(\rr)+v_{\mathrm{H}}[\rho_{\mathrm{in}}](\rr)+v_{xc}[\rho_{\mathrm{in}}](\rr) \label{eq:11-veff-in} \end{equation} $$を作り，1電子方程式
$$ \begin{equation} \hat h[\rho_{\mathrm{in}}]\,\phi_i=\Bigl(-\tfrac12\nabla^2+v_{\mathrm{eff}}[\rho_{\mathrm{in}}]\Bigr)\phi_i=\varepsilon_i\phi_i \label{eq:11-h-in} \end{equation} $$を1回だけ解いて固有値 $\{\varepsilon_i\}$ と軌道 $\{\phi_i\}$ を得る．このとき
$$ \begin{equation} E_{\mathrm{Harris}}[\rho_{\mathrm{in}}] =\sum_{i}^{\mathrm{occ}}\varepsilon_i[\rho_{\mathrm{in}}] -E_{\mathrm{H}}[\rho_{\mathrm{in}}] -\int \dd^3r\;v_{xc}[\rho_{\mathrm{in}}](\rr)\,\rho_{\mathrm{in}}(\rr) +E_{xc}[\rho_{\mathrm{in}}] \label{eq:11-harris-def} \end{equation} $$をハリス汎関数と呼ぶ．第10章で導いた自己無撞着な全エネルギーの表式
$$ E=\sum_i\varepsilon_i-E_{\mathrm{H}}[\rho]-\int v_{xc}\rho+E_{xc}[\rho] $$とまったく同じ形であり，違いは「$\rho_{\mathrm{in}}$ が自己無撞着でなくてよい」という点だけである．したがって $\rho_{\mathrm{in}}=\rho_0$ のとき $E_{\mathrm{Harris}}[\rho_0]=E_0$ が成り立つ．
注意:出力密度は使わない
\eqref{eq:11-harris-def} の右辺には，対角化から得られる情報として固有値の和だけが現れ，残りの3項はすべて入力密度 $\rho_{\mathrm{in}}$ で評価される．「固有値は $\rho_{\mathrm{in}}$ が作るポテンシャルの固有値，二重数え補正も $\rho_{\mathrm{in}}$」という一貫性がハリス汎関数の要点である．もし補正項に出力密度 $\rho_{\mathrm{out}}=\sum_i\abs{\phi_i}^2$ を混ぜると，以下で示す2次の性質は失われる．
11.7.3 誤差が2次であることの証明
これから，$E_{\mathrm{Harris}}$ が自己無撞着密度 $\rho_0$ のまわりで停留する——すなわち1次の変分がゼロになる——ことを示す．そこから誤差が $\delta\rho$ の2次であることが従う．
数学ノート:1次摂動論による固有値和の変分
ハミルトニアン $\hat h$ に微小な摂動 $\delta v_{\mathrm{eff}}$ が加わったとき，非縮退な固有値の1次変化は
$$ \begin{equation} \delta\varepsilon_i=\braket{\phi_i|\delta v_{\mathrm{eff}}|\phi_i} =\int \dd^3r\;\delta v_{\mathrm{eff}}(\rr)\,\abs{\phi_i(\rr)}^2 \label{eq:11-pert-1st} \end{equation} $$である(1次摂動論)．これは11.2.3節のHellmann–Feynman定理 \eqref{eq:11-hf-thm} を $\hat h$ に適用したものと同じで，固有関数の変化による寄与が規格化条件で消える．占有軌道について和をとると
$$ \begin{equation} \delta\Bigl(\sum_i^{\mathrm{occ}}\varepsilon_i\Bigr) =\int \dd^3r\;\delta v_{\mathrm{eff}}(\rr)\sum_i^{\mathrm{occ}}\abs{\phi_i(\rr)}^2 =\int \dd^3r\;\delta v_{\mathrm{eff}}(\rr)\,\rho_{\mathrm{out}}(\rr). \label{eq:11-band-var} \end{equation} $$ここで $\rho_{\mathrm{out}}=\sum_i^{\mathrm{occ}}\abs{\phi_i}^2$ は出力密度である．「固有値和の変分は，ポテンシャルの変化に出力密度を掛けて積分したもの」——この単純な関係が以下のすべてを支える．
導出:$\delta E_{\mathrm{Harris}}/\delta\rho_{\mathrm{in}}$ の計算
入力密度を $\rho_{\mathrm{in}}\to\rho_{\mathrm{in}}+\delta\rho_{\mathrm{in}}$ と変化させ，\eqref{eq:11-harris-def} の4項それぞれの1次変化を求める．以下，$v_{\mathrm{H}}$，$v_{xc}$ は特に断らない限り $\rho_{\mathrm{in}}$ で評価した値とする．
(1) 固有値和．$v$ は固定なので $\delta v_{\mathrm{eff}}=\delta v_{\mathrm{H}}+\delta v_{xc}$．\eqref{eq:11-band-var} より
$$ \begin{equation} \delta\Bigl(\sum_i\varepsilon_i\Bigr) =\int \dd^3r\;\bigl[\delta v_{\mathrm{H}}(\rr)+\delta v_{xc}(\rr)\bigr]\rho_{\mathrm{out}}(\rr). \label{eq:11-hv-1} \end{equation} $$ここで
$$ \delta v_{\mathrm{H}}(\rr)=\int \dd^3r'\frac{\delta\rho_{\mathrm{in}}(\rr')}{\abs{\rr-\rr'}}=v_{\mathrm{H}}[\delta\rho_{\mathrm{in}}](\rr), \qquad \delta v_{xc}(\rr)=\int \dd^3r'\,K_{xc}(\rr,\rr')\,\delta\rho_{\mathrm{in}}(\rr') $$である($K_{xc}=\delta^2E_{xc}/\delta\rho\delta\rho'$ は交換相関カーネル．2階の汎関数微分なので $K_{xc}(\rr,\rr')=K_{xc}(\rr',\rr)$ と対称である)．
(2) $-E_{\mathrm{H}}[\rho_{\mathrm{in}}]$．第10章の結果 $\delta E_{\mathrm{H}}/\delta\rho=v_{\mathrm{H}}$ から
$$ \begin{equation} \delta\bigl(-E_{\mathrm{H}}\bigr)=-\int \dd^3r\;v_{\mathrm{H}}(\rr)\,\delta\rho_{\mathrm{in}}(\rr). \label{eq:11-hv-2} \end{equation} $$(3) $E_{xc}[\rho_{\mathrm{in}}]$．定義より $\delta E_{xc}/\delta\rho=v_{xc}$ だから
$$ \begin{equation} \delta E_{xc}=\int \dd^3r\;v_{xc}(\rr)\,\delta\rho_{\mathrm{in}}(\rr). \label{eq:11-hv-3} \end{equation} $$(4) $-\int v_{xc}\rho_{\mathrm{in}}$．$v_{xc}$ も $\rho_{\mathrm{in}}$ も変化するので，積の変分で2項出る:
$$ \begin{equation} \delta\left(-\int v_{xc}\rho_{\mathrm{in}}\right) =-\int \dd^3r\;\delta v_{xc}(\rr)\,\rho_{\mathrm{in}}(\rr) -\int \dd^3r\;v_{xc}(\rr)\,\delta\rho_{\mathrm{in}}(\rr). \label{eq:11-hv-4} \end{equation} $$合計する．\eqref{eq:11-hv-1}〜\eqref{eq:11-hv-4} を足す．まず (3) の第1項と (4) の第2項が完全に相殺する($+\int v_{xc}\delta\rho_{\mathrm{in}}-\int v_{xc}\delta\rho_{\mathrm{in}}=0$)．残りは
\begin{align} \delta E_{\mathrm{Harris}} &=\int \delta v_{\mathrm{H}}\,\rho_{\mathrm{out}} +\int \delta v_{xc}\,\rho_{\mathrm{out}} -\int v_{\mathrm{H}}\,\delta\rho_{\mathrm{in}} -\int \delta v_{xc}\,\rho_{\mathrm{in}} \label{eq:11-hv-5}\\ &=\underbrace{\left[\int \delta v_{\mathrm{H}}\,\rho_{\mathrm{out}}-\int v_{\mathrm{H}}\,\delta\rho_{\mathrm{in}}\right]}_{(\mathrm{A})} +\underbrace{\int \delta v_{xc}\,\bigl(\rho_{\mathrm{out}}-\rho_{\mathrm{in}}\bigr)}_{(\mathrm{B})}. \label{eq:11-hv-6} \end{align}(A) の整理．クーロン核の対称性を使う:
\begin{align} \int \dd^3r\,\delta v_{\mathrm{H}}(\rr)\rho_{\mathrm{out}}(\rr) &=\iint \dd^3r\,\dd^3r'\;\frac{\delta\rho_{\mathrm{in}}(\rr')\,\rho_{\mathrm{out}}(\rr)}{\abs{\rr-\rr'}} \label{eq:11-hvA-1}\\ &=\int \dd^3r'\;\delta\rho_{\mathrm{in}}(\rr') \underbrace{\int \dd^3r\frac{\rho_{\mathrm{out}}(\rr)}{\abs{\rr-\rr'}}}_{=\,v_{\mathrm{H}}[\rho_{\mathrm{out}}](\rr')} \label{eq:11-hvA-2}\\ &=\int \dd^3r'\;v_{\mathrm{H}}[\rho_{\mathrm{out}}](\rr')\,\delta\rho_{\mathrm{in}}(\rr'). \label{eq:11-hvA-3} \end{align}\eqref{eq:11-hvA-1} は $\delta v_{\mathrm{H}}$ の定義を代入したもの．\eqref{eq:11-hvA-2} では積分順序を入れ替え，$\rr$ 積分を先に実行した(クーロン核が $\rr\leftrightarrow\rr'$ 対称なので，これは $\rho_{\mathrm{out}}$ が作るハートリーポテンシャルにほかならない)．よって
$$ \begin{equation} (\mathrm{A})=\int \dd^3r\;\Bigl(v_{\mathrm{H}}[\rho_{\mathrm{out}}]-v_{\mathrm{H}}[\rho_{\mathrm{in}}]\Bigr)\delta\rho_{\mathrm{in}} =\int \dd^3r\;v_{\mathrm{H}}\bigl[\rho_{\mathrm{out}}-\rho_{\mathrm{in}}\bigr]\,\delta\rho_{\mathrm{in}} \label{eq:11-hvA-4} \end{equation} $$($v_{\mathrm{H}}$ は密度について線形なので差をまとめられる)．
(B) の整理．$K_{xc}$ の対称性を使う:
\begin{align} (\mathrm{B})&=\iint \dd^3r\,\dd^3r'\;K_{xc}(\rr,\rr')\,\delta\rho_{\mathrm{in}}(\rr')\, \bigl(\rho_{\mathrm{out}}-\rho_{\mathrm{in}}\bigr)(\rr) \label{eq:11-hvB-1}\\ &=\int \dd^3r'\;\delta\rho_{\mathrm{in}}(\rr') \int \dd^3r\;K_{xc}(\rr',\rr)\bigl(\rho_{\mathrm{out}}-\rho_{\mathrm{in}}\bigr)(\rr). \label{eq:11-hvB-2} \end{align}結論．両者をまとめ，演算子 $K\equiv K_{\mathrm{H}}+K_{xc}$($K_{\mathrm{H}}(\rr,\rr')=1/\abs{\rr-\rr'}$)を導入すると
$$ \begin{equation} \frac{\delta E_{\mathrm{Harris}}}{\delta\rho_{\mathrm{in}}(\rr)} =\Bigl[K\bigl(\rho_{\mathrm{out}}-\rho_{\mathrm{in}}\bigr)\Bigr](\rr). \label{eq:11-harris-gradient} \end{equation} $$両方の項が共通因子 $(\rho_{\mathrm{out}}-\rho_{\mathrm{in}})$ を持つ．自己無撞着解では定義により $\rho_{\mathrm{out}}=\rho_{\mathrm{in}}=\rho_0$ だから
$$ \begin{equation} \frac{\delta E_{\mathrm{Harris}}}{\delta\rho_{\mathrm{in}}}\bigg|_{\rho_{\mathrm{in}}=\rho_0}=0. \label{eq:11-harris-stationary} \end{equation} $$すなわち $E_{\mathrm{Harris}}$ は自己無撞着密度で停留する．$E_{\mathrm{Harris}}[\rho_0]=E_0$ と合わせて
$$ \begin{equation} E_{\mathrm{Harris}}[\rho_0+\delta\rho]=E_0+O\bigl(\delta\rho^2\bigr). \label{eq:11-harris-2nd} \end{equation} $$ ∎11.7.4 誤差の符号:ハリス汎関数は下から近づく
2次の誤差の符号まで分かると，SCFの真値を上下から挟むことができる．線形応答の範囲で符号を決定しよう．
導出:2次項の符号
$\delta n\equiv\rho_{\mathrm{in}}-\rho_0$ とする．有効ポテンシャルの変化は $\delta v_{\mathrm{eff}}=K\,\delta n$(上の $K$ の定義そのもの)．この摂動に対する非相互作用系の密度応答は，応答関数 $\chi_0$ を使って
$$ \begin{equation} \rho_{\mathrm{out}}-\rho_0=\chi_0\,\delta v_{\mathrm{eff}}=\chi_0K\,\delta n \label{eq:11-chi0-response} \end{equation} $$と書ける(第8章の線形応答)．ここで $\chi_0$ は負定値である:ポテンシャルを上げれば電子は逃げるので，任意の $f$ に対し $\braket{f,\chi_0 f}\lt0$．また $K_{\mathrm{H}}$(クーロン核)は正定値である:$\braket{f,K_{\mathrm{H}}f}=\iint f(\rr)f(\rr')/\abs{\rr-\rr'}$ は電荷分布 $f$ の静電自己エネルギーの2倍であり正である．$K_{xc}$ は一般に小さいので，$K=K_{\mathrm{H}}+K_{xc}$ も正定値としてよい．
これを \eqref{eq:11-harris-gradient} に代入すると，勾配は
\begin{align} \frac{\delta E_{\mathrm{Harris}}}{\delta\rho_{\mathrm{in}}} &=K\bigl(\rho_{\mathrm{out}}-\rho_{\mathrm{in}}\bigr) =K\bigl[(\rho_{\mathrm{out}}-\rho_0)-(\rho_{\mathrm{in}}-\rho_0)\bigr] \label{eq:11-2nd-1}\\ &=K\bigl[\chi_0K\,\delta n-\delta n\bigr] =\bigl(K\chi_0K-K\bigr)\delta n\;\equiv\;M\,\delta n \label{eq:11-2nd-2} \end{align}と $\delta n$ について線形になる．$K$ と $\chi_0$ が対称なので $M=K\chi_0K-K$ も対称であり，したがって $E_{\mathrm{Harris}}$ は $\delta n$ の2次形式として
$$ \begin{equation} E_{\mathrm{Harris}}[\rho_0+\delta n]=E_0+\frac{1}{2}\braket{\delta n,\,M\,\delta n}+O(\delta n^3) \label{eq:11-2nd-3} \end{equation} $$と書ける．$M$ の符号を調べる．任意の $f\ne0$ について
\begin{align} \braket{f,Mf}&=\braket{f,K\chi_0Kf}-\braket{f,Kf} \label{eq:11-2nd-4}\\ &=\braket{Kf,\;\chi_0\,(Kf)}-\braket{f,Kf}. \label{eq:11-2nd-5} \end{align}\eqref{eq:11-2nd-5} では $K$ の対称性を使って左側の $K$ をブラ側へ移した．第1項は $\chi_0$ が負定値なので負，第2項は $K$ が正定値なので $-\braket{f,Kf}\lt0$．よって $\braket{f,Mf}\lt0$，すなわち $M$ は負定値である．したがって
$$ \begin{equation} E_{\mathrm{Harris}}[\rho_{\mathrm{in}}]\le E_0. \label{eq:11-harris-below} \end{equation} $$ ∎一方，KS汎関数を非自己無撞着密度で評価した \eqref{eq:11-eks-variational} は変分原理により $E_{\mathrm{KS}}[\rho_{\mathrm{in}}]\ge E_0$ である．両者を合わせて
$$ \begin{equation} E_{\mathrm{Harris}}[\rho_{\mathrm{in}}]\ \le\ E_0\ \le\ E_{\mathrm{KS}}[\rho_{\mathrm{in}}] \label{eq:11-sandwich} \end{equation} $$となり，真値を上下から挟むことができる．しかも両者の誤差はともに2次で，大きさもほぼ同じである．実務的には，2つを平均すればさらに精度が上がる．
例:重ね合わせ原子密度からの1回対角化
結晶や分子について，各原子の孤立原子密度を単純に足し合わせた $\rho_{\mathrm{in}}=\sum_I\rho_I^{\mathrm{atom}}$ を作る．これは自己無撞着密度から数%ずれているが，$\delta\rho$ の2次の効果しかないので，ハリス汎関数で評価した全エネルギーは $10^{-3}$ Ha 程度の精度に収まることが多い．$\mathrm{Si}$ の結晶では，SCF計算の全エネルギーとハリス汎関数の値の差は原子あたり $1$ meV 以下になる．この処方は
- 大規模系の予備計算・構造探索のスクリーニング
- 分子動力学における力の高速評価(Harris–Foulkes MD)
- タイトバインディング型DFT(DFTB)の理論的基礎
- 希ガスのGordon–Kim型ポテンシャル，埋め込み原子法(EAM)，結合次数ポテンシャル(BOP)といったモデルポテンシャルの導出
の基礎になっている．
11.7.5 局所力定理と第二変分法
ハリス汎関数の議論をもう一歩推し進めると，実務上きわめて有用な定理が得られる．
定理:局所力定理(local force theorem)
自己無撞着解の近傍で，有効ポテンシャルに微小な変化 $\delta v_{\mathrm{eff}}$ を加えたとき，全エネルギーの変化はバンドエネルギー(占有固有値の和)の変化だけで与えられる:
$$ \begin{equation} \Delta E_{\mathrm{KS}}=\sum_i^{\mathrm{occ}}\varepsilon_i\bigl[v_{\mathrm{eff}}+\delta v_{\mathrm{eff}}\bigr] -\sum_i^{\mathrm{occ}}\varepsilon_i\bigl[v_{\mathrm{eff}}\bigr]+O(\delta^2). \label{eq:11-local-force} \end{equation} $$論拠．$E_{\mathrm{KS}}$ を $\rho$ と $v_{\mathrm{eff}}$ の形式的に独立な2変数の関数とみなす．自己無撞着解では $\delta E/\delta\rho=0$ が成り立つので，密度の変化 $\delta\rho$ に由来する寄与は1次で消える．残るのは $v_{\mathrm{eff}}$ の陽な依存性だけであり，$v_{\mathrm{eff}}$ は全エネルギー表式 \eqref{eq:11-harris-def} のうち固有値和にしか入っていない．したがって \eqref{eq:11-local-force} が従う．誤差は $O(\delta\rho^2,\delta v_{\mathrm{eff}}^2,\delta\rho\,\delta v_{\mathrm{eff}})$ である．
例:局所力定理の応用
(a) 第二変分法による磁気異方性エネルギー(MAE)．磁化の向きを変えたときのエネルギー差は $0.1$ meV/atom というきわめて小さい量で，スピン軌道相互作用(SOI)に起因する．愚直にやるなら磁化方向ごとにSOI込みのSCFを収束させねばならないが，非常に高価であり，また数値誤差がMAEと同程度になりかねない．局所力定理を使えば，SOIなしでSCFを1回収束させ，その $v_{\mathrm{eff}}$ にSOI項を加えて1回だけ対角化し，固有値和の差からMAEを読み取れる．誤差は2次なので，この処方は完全SCFの結果とよく一致する．
(b) 交換相互作用定数 $J_{ij}$．ハイゼンベルグ模型 $H=-\sum_{i\ne j}J_{ij}\,\bm{e}_i\cdot\bm{e}_j$ の $J_{ij}$ を，強磁性状態のまわりで局所的なスピン回転を微小摂動として与え，バンドエネルギーの変化から評価する(Liechtensteinの方法)．得られた $J_{ij}$ を平均場近似で処理してキュリー温度を見積もると，bcc鉄で約 $1300$ K(実験 $1043$ K)，fccニッケルで約 $450$ K(実験 $627$ K)といった値が得られる．平均場近似が横磁化ゆらぎを無視するため系統的なずれが残るが，材料間の傾向は正しく再現される．
(c) フリーデル模型的な安定性の議論．第13章で扱う遷移金属の結晶構造安定性の議論——$d$ バンドの状態密度の形だけから構造安定性を論じる——も，局所力定理により正当化される．バンドエネルギーの差が全エネルギーの差を代表するからこそ，状態密度の形の議論が意味を持つ．
注意:2次だから安心，ではない
「誤差が2次」は「誤差が小さい」を保証しない．2次係数 $M$ が大きければ，$\delta\rho$ が小さくても誤差は無視できない．特に
- 電荷移動が大きい系(イオン結晶，金属/半導体界面)では，重ね合わせ原子密度が自己無撞着密度から大きくずれる
- 金属で $\chi_0$ が大きい系(状態密度がフェルミ準位で大きい)では応答が強く，$M$ が大きい
といった場合には，非自己無撞着計算は危険である．逆に，共有結合性の半導体や希ガス固体のように電荷移動の少ない系では，驚くほどよく機能する．使う前に，対象系でSCF計算と比較して検証しておくのが鉄則である．
この節で得られたもの．ハリス汎関数の誤差が入力密度の誤差の2次であること，しかも真値に下から近づくことを証明した．これにより，1回の対角化で高精度なエネルギーが得られる根拠と，その適用限界が明らかになった．局所力定理はこの性質を「ポテンシャルの摂動」に読み替えたもので，磁気異方性や交換相互作用定数の実用的な評価法を与える．
11.8 まとめ・演習・参考文献
11.8.1 実践的な汎関数選択の指針
本章の内容を，実際に計算するときの判断材料として整理する．
| 目的・対象 | まずこれ | 検証に | 理由(本章の該当節) |
|---|---|---|---|
| 分子の構造・結合エネルギー | PBE，B3LYP | SCAN，PBE0 | LDAは過結合(11.5.2) |
| 固体の格子定数・弾性定数 | PBEsol，SCAN | PBE，LDA(上下の目安) | PBEは過大，LDAは過小(11.5.3) |
| 磁性金属・磁気秩序 | PBE | SCAN，PBE+U | LDAはbcc鉄を誤る(11.5.4) |
| バンドギャップ | HSE06 | $GW$(第18章)，DFT+U | $\Delta_{xc}=0$ とSIE(11.6.1) |
| 層状物質・分子結晶・物理吸着 | vdW-DF系，PBE-D3 | rev-vdW-DF2，SCAN+rVV10 | 準局所汎関数はvdWを与えない(11.6.4) |
| 遷移金属酸化物・強相関系 | DFT+U | HSE06，DMFT | 局在 $d$/$f$ 電子のSIE(11.6.5) |
| 反応障壁・遷移状態 | ハイブリッド(PBE0など) | meta-GGA | SIEが障壁を下げる(11.3.5) |
| 金属表面・表面エネルギー | PBEsol，LDA | SCAN | PBEは表面エネルギーを過小(11.5.5) |
| 大規模系のスクリーニング | PBE(非自己無撞着でも可) | SCFで再計算 | ハリス汎関数の2次性(11.7) |
物理的意味:汎関数は道具である
本章で見てきたとおり，あらゆる系であらゆる物理量を高精度に与える汎関数は存在しない．梯子を1段上がるたびに新しい物理が取り込まれるが，それと引き換えに計算コストが増え，しばしば別の物理量が悪化する．PBEが固体の格子定数を過大評価し，PBEsolが分子の結合エネルギーを悪化させるという「あちら立てればこちら立たず」の関係は，パラメータ調整の未熟さではなく，限られた材料($\rho$ と $\nabla\rho$)で複数の厳密条件を同時に満たすことの原理的な難しさに由来する．
したがって実践者に求められる態度は3つである．
- 系に応じて選ぶ．表11.9は出発点にすぎない．対象系に近い既知の物質で，汎関数ごとの誤差の大きさと向きを確かめてから本番に進む．
- 誤差の向きを知ったうえで読む．「PBEの格子定数だから $1\%$ 大きめ」「PBEのギャップだから半分」と補正して解釈できれば，多くの場合それで十分な判断ができる．
- 複数の汎関数で検証する．結論が汎関数に依存しないなら信頼できる．依存するなら，その依存性こそが物理的に重要な情報である(たとえば「LDAとGGAで基底状態が変わる」なら，その系は磁性と局在性の微妙な競合の上にある)．
汎関数は自然法則ではなく，近似という道具である．道具の癖を知り，目的に応じて選び，結果を検証する——これが第一原理計算を使いこなすということである．
11.8.2 まとめ
- LDA．$E_{xc}^{\mathrm{LDA}}=\int\rho\,\varepsilon_{xc}(\rho)\dd^3r$(式 \eqref{eq:11-lda-def})．交換部分は解析的に $\varepsilon_x=-\tfrac34(3/\pi)^{1/3}\rho^{1/3}$，相関部分はジェリウムのQMC結果をPZ/VWN形にフィットしたもの．交換相関ポテンシャルは $v_{xc}=\varepsilon_{xc}+\rho\,\dd\varepsilon_{xc}/\dd\rho$(式 \eqref{eq:11-lda-key})であり，交換だけなら $v_x=-(3/\pi)^{1/3}\rho^{1/3}=\tfrac43\varepsilon_x$．スピン分極版(LSDA)では交換のスピンスケーリング関係 \eqref{eq:11-spin-scaling} が厳密に成り立つが，相関の $\zeta$ 内挿は近似である．
- 断熱接続．相互作用の強さ $\lambda$ を $0\to1$ と変えながら密度を固定する外部ポテンシャル $v_\lambda$ を考え，Hellmann–Feynman定理を $\lambda$ 微分に適用すると $E_{xc}=\int_0^1\dd\lambda\braket{\Psi_\lambda|\hat V_{ee}|\Psi_\lambda}-E_{\mathrm{H}}$(式 \eqref{eq:11-ac-key})．運動エネルギーの相関寄与 $T-T_s$ は $W_\lambda$ の $\lambda$ 依存性として吸収されている．$\lambda=0$ 端の値 $W_0$ は厳密交換 $E_x$ である．
- 交換相関ホール．$E_{xc}=\frac12\iint\rho(\rr)\bar n_{xc}(\rr,\rr')/\abs{\rr-\rr'}$(式 \eqref{eq:11-exc-hole-key})——密度と「電荷 $-1$ の穴」との静電相互作用．総和則 $\int n_x\dd^3r'=-1$，$\int n_c\dd^3r'=0$(式 \eqref{eq:11-sumrules-key})．
- 球平均への帰着．クーロン核 $1/s$ が等方的なので角度積分がホールを平均し，$E_{xc}=\frac12\int\dd^3r\,\rho\int_0^\infty4\pi s\,\dd s\,\bar n_{xc}^{\mathrm{SA}}$(式 \eqref{eq:11-sph-key})．$E_{xc}$ はホールの形ではなく球平均にしか依存しない．LDAホールは形を大きく誤るが球平均が近く総和則も満たすので，$E_{xc}$ の誤差は小さい．これがLDAの成功の内実である．
- 厳密な制約．交換の一様スケーリング $E_x[\rho_\gamma]=\gamma E_x[\rho]$(式 \eqref{eq:11-ex-scaling})，相関のスケーリング不等式 \eqref{eq:11-ec-scaling}，Lieb–Oxford下限 $E_{xc}\ge-1.68\int\rho^{4/3}$，Levy–Perdew関係式 $E_x=-\int\rho\,\rr\cdot\nabla v_x$(式 \eqref{eq:11-lp-key})，1電子系での自己相互作用の消去 $E_{xc}=-E_{\mathrm{H}}$．LDA/GGAは最後の条件を破り，水素原子で $0.6$ eV の誤差を残す．
- GEAの失敗とGGA．2次勾配展開はホールの総和則と符号条件を破るため，LDAより悪化する．実空間カットオフで総和則を回復させたものがGGAである．「展開の次数より制約の充足」が教訓．
- PBE．$F_x(s)=1+\kappa-\kappa/(1+\mu s^2/\kappa)$(式 \eqref{eq:11-pbe-x-key})．$\kappa=0.804$ はLieb–Oxford下限から，$\mu=\beta\pi^2/3=0.21951$ は交換と相関の勾配補正が打ち消し合ってLDAの線形応答を回復する条件から決まる．相関の $H$ 項は $t\to0$ で勾配展開，$t\to\infty$ で相関消失，高密度極限で $\ln$ の相殺により有限——3条件をすべて満たす．経験的パラメータはゼロ．
- 数値の傾向．LSDAは交換を $10\%$ 浅く，相関を $2$〜$3$ 倍深く見積もる(誤差が逆向きに相殺)．分子では過結合(結合1本あたり $1$ eV 超)，固体では格子定数を $1$〜$2\%$ 過小評価．PBEは誤差を1桁縮めるが，固体の格子定数は逆に $1$〜$2\%$ 過大になる．bcc鉄の基底状態はLDAでは誤り，GGAで正しく得られる．
- ヤコブの梯子．LDA($\rho$)→ GGA($+\nabla\rho$)→ meta-GGA($+\tau$)→ ハイブリッド($+$占有軌道=厳密交換)→ RPA($+$非占有軌道)．ハイブリッドの混合比 $a=1/4$ は断熱接続曲線を $(1-\lambda)^{n-1}$ でモデル化して $n=4$ とすることから出る．HSEは $1/r=\mathrm{erfc}(\omega r)/r+\mathrm{erf}(\omega r)/r$ と分割し短距離だけに厳密交換を混ぜる．
- vdWとDFT+U．密度の重ならない2体間で準局所汎関数は厳密に加算的になり，相互作用エネルギーがゼロになる——だからvdWが出ない．処方は非局所相関汎関数(vdW-DF)か経験的補正(DFT-D)．局在 $d$/$f$ 電子には $E_U=\frac{U_{\mathrm{eff}}}{2}\sum n(1-n)$(式 \eqref{eq:11-dftu-key})が分数占有に罰金を課し，ギャップを $U_{\mathrm{eff}}$ だけ開き，軌道分極を促す．二重数え補正が必須．
- ハリス汎関数．$E_{\mathrm{Harris}}[\rho_{\mathrm{in}}]=\sum_i\varepsilon_i-E_{\mathrm{H}}-\int v_{xc}\rho_{\mathrm{in}}+E_{xc}$(式 \eqref{eq:11-harris-key})．変分の勾配が $K(\rho_{\mathrm{out}}-\rho_{\mathrm{in}})$ という形になるので自己無撞着密度で停留し，誤差は $\delta\rho$ の2次．2次形式は負定値なので下から近づき，KS汎関数(上から)と真値を挟む．局所力定理はその応用で，磁気異方性や交換相互作用定数の実用的な評価法を与える．
次章への橋渡し．第9〜11章で，密度汎関数理論の枠組み(HK定理)，実用形式(KS方程式)，そして最後の未知項の近似($E_{xc}$)がすべて揃った．原理的には，任意の原子配置に対して電子密度と全エネルギーを計算できる．しかし実際の固体では原子が $10^{23}$ 個も並んでおり，そのままでは有限の計算に落ちない．ここで効いてくるのが周期性である．第12章では，周期ポテンシャル中の電子波動関数がブロッホの定理に従うこと，その結果として無限系の問題が単位胞1個分の有限の問題に還元されること，そしてエネルギー準位が連続的な「バンド」を作ることを示す．第11章までに作り上げた $v_{\mathrm{eff}}$ が周期的であるという，ただそれだけの事実から，固体物理の全体像が立ち上がる．
11.8.3 演習問題
演習11.1 LDA交換ポテンシャルと $X\alpha$ 法
(1) 一様電子ガスの交換エネルギー密度 $\varepsilon_x=-\tfrac34(3/\pi)^{1/3}\rho^{1/3}$ から出発し，$v_x^{\mathrm{LDA}}=-(3/\pi)^{1/3}\rho^{1/3}$ を導け．またこれを $r_s$ で表し，$v_x=-0.6109/r_s$ Ha となることを確かめよ．
(2) $X\alpha$ 法のポテンシャル $v_{x\alpha}=-\frac32\alpha(3\rho/\pi)^{1/3}$ が，$\alpha=2/3$ のとき (1) の結果に一致することを示せ．
(3) $X\alpha$ 法の交換エネルギー汎関数 $E_{x\alpha}[\rho]=-C\int\rho^{4/3}\dd^3r$ が，その汎関数微分として $v_{x\alpha}$ を与えるように定数 $C$ を $\alpha$ の関数として決めよ．$\alpha=2/3$ で $C=C_x$ に戻ることを確認せよ．
ヒント:(1) 11.1.4節の導出をなぞる．$r_s$ への変換は $\rho=3/(4\pi r_s^3)$，すなわち $\rho^{1/3}=(3/4\pi)^{1/3}/r_s=0.62035/r_s$．(3) $E=-C\int\rho^{4/3}$ の汎関数微分は $-\frac43C\rho^{1/3}$ だから，$-\frac43C=-\frac32\alpha(3/\pi)^{1/3}$ を解けばよい．
演習11.2 一様電子ガスの交換ホールと総和則
スピン無分極の一様電子ガスにおいて，交換ホールは
$$ n_x(s)=-\frac{9\rho}{2}\left[\frac{\sin(k_Fs)-k_Fs\cos(k_Fs)}{(k_Fs)^3}\right]^2 $$で与えられる($s=\abs{\rr-\rr'}$)．
(1) $s\to0$ の極限で $n_x\to-\rho/2$ となることを示せ．この値の物理的意味を述べよ(なぜ $-\rho$ ではないのか)．
(2) 総和則 $\int_0^\infty 4\pi s^2\,n_x(s)\,\dd s=-1$ が成り立つことを，$x=k_Fs$ と置換したうえで公式 $\int_0^\infty\frac{(\sin x-x\cos x)^2}{x^4}\dd x=\frac{\pi}{4}$ を使って確かめよ．
(3) $n_x(s)$ を式 \eqref{eq:11-sph-key} の左側に代入し，公式 $\int_0^\infty\frac{(\sin x-x\cos x)^2}{x^5}\dd x=\frac14$ を使って $\varepsilon_x=-3k_F/(4\pi)$ を導け．
ヒント:(1) $\sin y-y\cos y=y^3/3-y^5/30+\cdots$ を使うと，括弧内は $y\to0$ で $1/3$ に近づく．$-\frac92\rho\cdot\frac19=-\frac\rho2$．$-\rho/2$ になるのは，交換ホールが同スピン電子(全体の半分)しか排除しないからである．(2)(3) では $\rho=k_F^3/(3\pi^2)$ を使う．(3) の重みが $s$(すなわち $x$ の1乗少ない)なので，被積分関数の $x$ の冪が (2) より1つ下がることに注意．
演習11.3 PBE交換増強因子の性質
$F_x(s)=1+\kappa-\dfrac{\kappa}{1+\mu s^2/\kappa}$，$\kappa=0.804$，$\mu=0.21951$ とする．
(1) $F_x(1)$，$F_x(2)$，$F_x(\infty)$ を数値で求めよ．
(2) $\dd F_x/\dd s$ を計算し，$F_x$ が $s$ の単調増加関数であることを示せ．また $F_x$ の変曲点(2階微分がゼロになる $s$)を求めよ．
(3) PBEsol は $\mu=10/81$，$\kappa=0.804$ である．$s=1$ における $F_x^{\mathrm{PBE}}$ と $F_x^{\mathrm{PBEsol}}$ の差を求め，それが固体の格子定数にどう効くかを $50$ 字程度で述べよ．
(4) $F_x$ を $s\to\infty$ で $1+\kappa$ に飽和させることが必須である理由を，原子の裾での $s$ の振る舞い(11.4.3節の例)を使って説明せよ．
ヒント:(2) $u=\mu s^2/\kappa$ と置くと $F_x=1+\kappa-\kappa/(1+u)$，$\dd F_x/\dd s=\kappa\frac{\dd u/\dd s}{(1+u)^2}$．$\dd u/\dd s=2\mu s/\kappa\ge0$ なので単調増加．(3) $F_x$ が大きいほど密度の不均一性が有利になり，原子まわりに密度が集まって原子間距離が伸びる．
演習11.4 ハリス汎関数の1次変分
簡単のため $E_{xc}=0$(ハートリー近似)の場合を考える．このとき
$$ E_{\mathrm{Harris}}[\rho_{\mathrm{in}}]=\sum_i^{\mathrm{occ}}\varepsilon_i[\rho_{\mathrm{in}}]-E_{\mathrm{H}}[\rho_{\mathrm{in}}] $$である．
(1) 11.7.3節の導出をなぞって，$\delta E_{\mathrm{Harris}}/\delta\rho_{\mathrm{in}}=v_{\mathrm{H}}[\rho_{\mathrm{out}}-\rho_{\mathrm{in}}]$ を示せ．
(2) 比較のため，誤って補正項に出力密度を使った汎関数
$$ \tilde E[\rho_{\mathrm{in}}]=\sum_i^{\mathrm{occ}}\varepsilon_i[\rho_{\mathrm{in}}]-E_{\mathrm{H}}[\rho_{\mathrm{out}}] $$を考える．この汎関数の1次変分が一般にゼロにならないこと(すなわち2次の精度を失うこと)を示せ．
(3) 自己無撞着解のまわりで $\rho_{\mathrm{out}}-\rho_0=\chi_0K_{\mathrm{H}}\delta n$ と書けるとして，$E_{\mathrm{Harris}}$ の2次形式が負定値であることを再確認せよ．
ヒント:(1) 11.7.3節の (A) の計算だけを行えばよい．(2) $E_{\mathrm{H}}[\rho_{\mathrm{out}}]$ の変分には $\delta\rho_{\mathrm{out}}=\chi_0\delta v_{\mathrm{eff}}$ という応答が現れ，$\rho_{\mathrm{out}}=\rho_{\mathrm{in}}$ を代入しても消えない項が残る．(3) $\braket{f,K\chi_0Kf}\lt0$ と $-\braket{f,Kf}\lt0$ を使う．
演習11.5 DFT+Uによる軌道分極と準位シフト
1つの原子の $d$ 軌道($m=1,\dots,5$)に，あるスピンチャネルの電子が $N_d$ 個入っているとする．$E_U=\frac{U_{\mathrm{eff}}}{2}\sum_{m}n_m(1-n_m)$，$\sum_m n_m=N_d$ とする．
(1) $N_d=2$ のとき，$E_U$ を最小にする占有パターンを求めよ．また最大にする占有パターンを求め，両者のエネルギー差を $U_{\mathrm{eff}}$ で表せ．
(2) $N_d$ が整数のとき，$E_U$ の最小値が常にゼロであることを示せ．
(3) $v_U^m=\partial E_U/\partial n_m$ を求め，占有軌道と空軌道のエネルギー差が $U_{\mathrm{eff}}$ だけ広がることを示せ．
(4) 二重数え補正 $E^{\mathrm{dc}}=\frac{U}{2}N_d(N_d-1)$ を引かなかった場合，$E_U$ はどのような形になるか．そのとき (2) の結論はどう変わるか．
ヒント:(1) $n(1-n)$ は $n=0,1$ でゼロ，$n=1/2$ で最大．$\sum n_m=2$ の拘束のもとでは $(1,1,0,0,0)$ が最小，$(1/2\times4,0)$ すなわち4軌道に $1/2$ ずつが最大．(4) 二重数えを引かないと $E^{\mathrm{Hub}}=\frac{U}{2}[N_d^2-\sum n_m^2]$ のままで，$N_d^2$ という占有パターンによらない項が残るうえ，$-\frac U2\sum n_m^2$ は分数占有を安定化する向きに働く．
11.8.4 参考文献
- W. Kohn and L. J. Sham, Self-Consistent Equations Including Exchange and Correlation Effects, Phys. Rev. 140, A1133 (1965). ——KS方程式の原論文であると同時に，LDAを最初に提案した論文でもある(11.1節)．
- D. M. Ceperley and B. J. Alder, Ground State of the Electron Gas by a Stochastic Method, Phys. Rev. Lett. 45, 566 (1980). ——一様電子ガスの相関エネルギーのQMC計算．LDAの数値的基礎(11.1.3節)．
- S. H. Vosko, L. Wilk, and M. Nusair, Accurate spin-dependent electron liquid correlation energies for local spin density calculations, Can. J. Phys. 58, 1200 (1980). ——VWNフィット(式 \eqref{eq:11-vwn})．
- J. P. Perdew and A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Phys. Rev. B 23, 5048 (1981). ——PZフィット(式 \eqref{eq:11-pz})と自己相互作用補正(11.3.5節)．
- J. P. Perdew and Y. Wang, Accurate and simple analytic representation of the electron-gas correlation energy, Phys. Rev. B 45, 13244 (1992). ——PW92．現在最も広く使われるLDA相関のパラメータ化．
- O. Gunnarsson and B. I. Lundqvist, Exchange and correlation in atoms, molecules, and solids by the spin-density-functional formalism, Phys. Rev. B 13, 4274 (1976). ——断熱接続と交換相関ホールによるLDAの正当化(11.2節)．
- D. C. Langreth and J. P. Perdew, Exchange-correlation energy of a metallic surface: Wave-vector analysis, Solid State Commun. 17, 1425 (1975); Phys. Rev. B 15, 2884 (1977). ——結合定数積分の定式化と，ホールの球平均がエネルギーを決めることの明示(11.2.6節)．
- O. Gunnarsson, M. Jonson, and B. I. Lundqvist, Descriptions of exchange and correlation effects in inhomogeneous electron systems, Phys. Rev. B 20, 3136 (1979). ——Ne原子における厳密ホールとLDAホールの比較(図11.4の元になった解析)．
- E. H. Lieb and S. Oxford, Improved lower bound on the indirect Coulomb energy, Int. J. Quantum Chem. 19, 427 (1981). ——Lieb–Oxford下限(式 \eqref{eq:11-LO})．
- M. Levy and J. P. Perdew, Hellmann-Feynman, virial, and scaling requisites for the exact universal density functionals, Phys. Rev. A 32, 2010 (1985). ——スケーリング則とビリアル型関係式(11.3節)．
- J. P. Perdew and Y. Wang, Accurate and simple density functional for the electronic exchange energy: Generalized gradient approximation, Phys. Rev. B 33, 8800 (1986). ——実空間カットオフによるGGAの構成(11.4.2節)．
- J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77, 3865 (1996). ——PBE汎関数．11.4節の主題．表11.4〜11.6の原子・分子の数値もこの論文に基づく．
- J. P. Perdew et al., Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces, Phys. Rev. Lett. 100, 136406 (2008). ——PBEsol．固体の格子定数・表面エネルギーの改善(11.4.5節，11.5.5節)．
- A. D. Becke, A new mixing of Hartree-Fock and local density-functional theories, J. Chem. Phys. 98, 1372 (1993); Density-functional thermochemistry. III. The role of exact exchange, J. Chem. Phys. 98, 5648 (1993). ——ハイブリッド汎関数の提唱とB3LYP(11.6.2節)．
- J. P. Perdew, M. Ernzerhof, and K. Burke, Rationale for mixing exact exchange with density functional approximations, J. Chem. Phys. 105, 9982 (1996). ——混合比 $a=1/4$ の断熱接続による根拠(式 \eqref{eq:11-a-quarter})．
- J. Heyd, G. E. Scuseria, and M. Ernzerhof, Hybrid functionals based on a screened Coulomb potential, J. Chem. Phys. 118, 8207 (2003); 124, 219906(E) (2006). ——HSE汎関数(式 \eqref{eq:11-hse})．
- J. Sun, A. Ruzsinszky, and J. P. Perdew, Strongly Constrained and Appropriately Normed Semilocal Density Functional, Phys. Rev. Lett. 115, 036402 (2015). ——SCAN汎関数(11.6.3節)．
- M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist, Van der Waals Density Functional for General Geometries, Phys. Rev. Lett. 92, 246401 (2004). ——vdW-DF(式 \eqref{eq:11-vdwdf})．
- S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, J. Chem. Phys. 132, 154104 (2010). ——DFT-D3(式 \eqref{eq:11-dftd})．
- V. I. Anisimov, J. Zaanen, and O. K. Andersen, Band theory and Mott insulators: Hubbard U instead of Stoner I, Phys. Rev. B 44, 943 (1991). ——LDA+U法の原論文(11.6.5節)．
- S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study, Phys. Rev. B 57, 1505 (1998). ——簡約形 \eqref{eq:11-dftu-key}．
- M. Cococcioni and S. de Gironcoli, Linear response approach to the calculation of the effective interaction parameters in the LDA+U method, Phys. Rev. B 71, 035105 (2005). ——$U$ の第一原理的決定．
- J. Harris, Simplified method for calculating the energy of weakly interacting fragments, Phys. Rev. B 31, 1770 (1985). ——ハリス汎関数(11.7節)．
- W. M. C. Foulkes and R. Haydock, Tight-binding models and density-functional theory, Phys. Rev. B 39, 12520 (1989). ——ハリス汎関数の2次性の詳細な解析とタイトバインディング法との関係．
- A. R. Mackintosh and O. K. Andersen, in Electrons at the Fermi Surface, ed. M. Springford (Cambridge University Press, 1980). ——局所力定理(式 \eqref{eq:11-local-force})．
- A. I. Liechtenstein, M. I. Katsnelson, V. P. Antropov, and V. A. Gubanov, Local spin density functional approach to the theory of exchange interactions in ferromagnetic metals and alloys, J. Magn. Magn. Mater. 67, 65 (1987). ——局所力定理による交換相互作用定数 $J_{ij}$ の評価．
- T. Asada and K. Terakura, Cohesive properties of iron obtained by use of the generalized gradient approximation, Phys. Rev. B 46, 13599 (1992). ——bcc鉄の基底状態におけるLDAとGGAの差(11.5.4節)．
- P. Haas, F. Tran, and P. Blaha, Calculation of the lattice constant of solids with semilocal functionals, Phys. Rev. B 79, 085104 (2009). ——固体の格子定数・体積弾性率の系統的比較(表11.7)．
- R. M. Martin, Electronic Structure: Basic Theory and Practical Methods (Cambridge University Press, 2004). ——固体物理寄りの標準教科書．第8章が本章の全範囲をカバーする．
- R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules (Oxford University Press, 1989). ——化学寄りの標準教科書．断熱接続と交換相関ホールの記述が詳しい．
- J. P. Perdew and K. Schmidt, Jacob's ladder of density functional approximations for the exchange-correlation energy, AIP Conf. Proc. 577, 1 (2001). ——ヤコブの梯子(図11.6)の原典．