test

17 minute read

密度汎関数理論 (DFT) における交換相関汎関数の理論的基盤と数学的定式化

1. Kohn-Sham理論と交換相関エネルギーの定義

1.1 多電子ハミルトニアンから非相互作用参照系へ

多電子系におけるシュレーディンガー方程式の厳密な解法は、電子数の増加に伴い次元の呪い（exponential scaling）に直面する。密度汎関数理論（Density Functional Theory: DFT）は、多電子系の波動関数 $\Psi(\mathbf{r}_1, \dots, \mathbf{r}_N)$ の代わりに、3次元空間の物理量である電子密度 $\rho(\mathbf{r})$ を基本変数とすることで、この困難を回避するアプローチである。

Hohenberg-Kohnの第1定理により、基底状態の全エネルギーは電子密度 $\rho(\mathbf{r})$ の一意な汎関数 $E[\rho]$ として記述できることが保証されている。しかし、電子密度のみから相互作用する電子系の運動エネルギーを直接高精度に近似することは困難であった（Thomas-Fermiモデルの限界）。これを解決したのが、1965年にKohnとShamによって提案されたKohn-Sham（KS）法である。

Kohn-Sham法では、実際の相互作用する電子系と同じ基底状態密度 $\rho(\mathbf{r})$ を与える、相互作用のない仮想的な独立粒子系（Kohn-Sham参照系）を導入する。この非相互作用系のハミルトニアンは以下のように与えられる。 \(\hat{H}_s = \sum_{i=1}^N \left( -\frac{1}{2}\nabla_i^2 + v_{\text{eff}}(\mathbf{r}_i) \right)\) ここで、ハートリー原子単位系（$\hbar = m_e = e = 4\pi\varepsilon_0 = 1$）を用いている。このハミルトニアンに対する独立粒子方程式がKohn-Sham方程式である。 \(\left( -\frac{1}{2}\nabla^2 + v_{\text{eff}}(\mathbf{r}) \right) \phi_i(\mathbf{r}) = \varepsilon_i \phi_i(\mathbf{r})\) ここで、$\phi_i(\mathbf{r})$ はKohn-Sham軌道であり、全電子密度 $\rho(\mathbf{r})$ は占有軌道の二乗和として再現される。 \(\rho(\mathbf{r}) = \sum_{i=1}^N |\phi_i(\mathbf{r})|^2\)

1.2 交換相関エネルギー $E_{xc}[\rho]$ の定義

Kohn-Sham法において、相互作用する系の全エネルギー汎関数 $E[\rho]$ は、以下のように各成分に分割して定義される。 \(E[\rho] = T_s[\rho] + \int \rho(\mathbf{r}) v_{\text{ext}}(\mathbf{r}) d\mathbf{r} + E_H[\rho] + E_{xc}[\rho]\) ここで、

  • $T_s[\rho]$ は非相互作用系における運動エネルギーであり、Kohn-Sham軌道を用いて以下のように厳密に計算される。 \(T_s[\rho] = -\frac{1}{2} \sum_{i=1}^N \int \phi_i^*(\mathbf{r}) \nabla^2 \phi_i(\mathbf{r}) d\mathbf{r}\)
  • $v_{\text{ext}}(\mathbf{r})$ は原子核による静電ポテンシャルなどの外部ポテンシャルである。
  • $E_H[\rho]$ は電子密度の古典的な静電反発エネルギー（Hartreeエネルギー）である。 \(E_H[\rho] = \frac{1}{2} \iint \frac{\rho(\mathbf{r})\rho(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r} d\mathbf{r}'\)
  • $E_{xc}[\rho]$ が交換相関エネルギー（Exchange-Correlation Energy）である。

系の真の運動エネルギーを $T[\rho]$、電子間の全相互作用エネルギーを $V_{ee}[\rho]$ とすると、真の全エネルギーは $E[\rho] = T[\rho] + \int \rho(\mathbf{r})v_{\text{ext}}(\mathbf{r})d\mathbf{r} + V_{ee}[\rho]$ と書けるため、交換相関エネルギー $E_{xc}[\rho]$ の厳密な定義式は以下のように導かれる。 \(E_{xc}[\rho] = (T[\rho] - T_s[\rho]) + (V_{ee}[\rho] - E_H[\rho])\) この式は、交換相関エネルギーが単なる電子間の非古典的な静電相互作用エネルギー（$V_{ee}[\rho] - E_H[\rho]$）だけでなく、実系と非相互作用系の運動エネルギーの差（$T[\rho] - T_s[\rho]$、すなわち運動エネルギー相関効果）をも内包していることを意味する。

交換相関ポテンシャル $v_{xc}(\mathbf{r})$ は、このエネルギー汎関数の電子密度に対する関数微分として定義され、Kohn-Sham方程式の実効ポテンシャル $v_{\text{eff}}(\mathbf{r}) = v_{\text{ext}}(\mathbf{r}) + v_H(\mathbf{r}) + v_{xc}(\mathbf{r})$ に寄与する。 \(v_{xc}(\mathbf{r}) = \frac{\delta E_{xc}[\rho]}{\delta \rho(\mathbf{r})}\)

交換効果と相関効果について

交換相関ポテンシャル/エネルギーは，交換部分と相関部分に分けて考えられる．物理的には「パウリの排他原理（スピンの都合）」によるものを交換，「クーロン反発（電気的な都合）」によるものを相関と呼ぶ．前者は同じスピンを持つ電子同士が反発し，後者は全ての電子同士が反発するから，いずれにしても交換相関効果を考慮すると電子の周りに他の電子が近づきにくくなる．このことを次の節で説明する「交換相関孔」と呼ぶ．

参考までに，量子化学における波動関数理論（HF法など）における交換相関との定義との違いを確認しておこう．HF法では電子の多体波動関数を1つのスレーター行列式（パウリの排他原理を満たす形）で表すため，交換効果が厳密に含まれている．これをExact Exchangeと呼ぶ．一方で電子相関についてはクーロン力を平均場で扱っているため瞬間的な相関効果を無視している．そこで量子化学では正確なエネルギーと、HF法で得られたエネルギーの差をを相関エネルギー$E_{\text{corr}}$と定義する．

一方でDFTではあくまで交換相関ポテンシャルは人間の手による近似であり，交換エネルギーも相関エネルギーも正確な値ではない．特に交換エネルギーは自己相互作用誤差（自分自身とクーロン反発してしまうエラー）を含んでしまうという欠点が知られている．また，DFTの相関項にはHF法が無視していた「クーロン反発によるリアルタイムの回避運動」に加えて「仮想的な非相互作用系から、現実の相互作用系へ変換するときに生じる、運動エネルギーのズレ」も不純物として一緒に含まれている．

項目 交換 相関
物理的原因 パウリの排他律 クーロン反発
HF法 数式的に厳密に考慮 完全に無視
DFT 汎函数による近似 汎函数による近似

DFTにおける数式的な定義は以下のようになる．

交換エネルギー（$E_{\text{x}}$）の定義交換エネルギーは、「もし電子が『スピンの都合（パウリの排他原理）』だけで動いていて、クーロン反発（マイナス同士の斥力）が一切なかったとしたら？」という仮定のもとで定義されます。クーロン反発がないということは、電子はリアルタイムでお互いを避け合わないため、運動エネルギーのズレ（$\Delta T$）は発生しません。波動関数は単一のスレーター行列式のままで変化しないからです。そのため、DFTにおける交換エネルギー $E_{\text{x}}$ は、KS軌道（単一スレーター行列式）から計算される純粋な静的（ポテンシャル）エネルギーの差分のみと定義されます。\(E_{\text{x}}[\rho] = \langle \Phi_{\text{KS}} \vert{} \hat{V}_{\text{ee}} \vert{} \Phi_{\text{KS}} \rangle - J[\rho]\)ここには運動エネルギーの項（$T$）は一切登場しません。2. 相関エネルギー（$E_{\text{c}}$）の定義全体の箱（$E_{\text{xc}}$）から、上記の「純粋な交換（$E_{\text{x}}$）」を引き算した残りのすべてが「相関（$E_{\text{c}}$）」と定義されます。\(E_{\text{c}}[\rho] = E_{\text{xc}}[\rho] - E_{\text{x}}[\rho]\)数式を整理すると、以下のようになります。\(E_{\text{c}}[\rho] = (T - T_{\text{s}}) + \left( V_{\text{ee}} - \langle \Phi_{\text{KS}} \vert{} \hat{V}_{\text{ee}} \vert{} \Phi_{\text{KS}} \rangle \right)\)ご覧の通り、運動エネルギーのズレ（$T - T_{\text{s}}$）がそっくりそのまま相関エネルギーの第一項として居座ることになります。補足：「なぜ交換に含まれる気がしてしまうのか？」の正体質問者様が「一部は交換に含まれるのでは？」と感じられた感覚は、物理的には非常に鋭いです。なぜなら、パウリの排他原理（交換）によっても、電子の運動（波の形）は制限を受けるからです。しかし、KS-DFTの巧妙なところは、「パウリの排他原理による運動エネルギーの変化分は、最初から $T_{\text{s}}$（仮想系の運動エネルギー）の中に100%完璧に織り込み済みである」という点にあります。KS法では、相互作用がなくても「電子がフェルミ粒子である（スレーター行列式に従う）」という前提で $T_{\text{s}}$ を計算します。つまり、スピンの都合による運動の制約はすでに $T_{\text{s}}$ が引き受けているため、差額である $\Delta T = T - T_{\text{s}}$ には「スピンの都合」は残っておらず、「クーロン反発のせいで余計に歪んだ分の運動エネルギー」だけが残るのです。したがって、運動エネルギーのズレは一部たりとも交換には混ざらず、すべて相関の責任（$E_{\text{c}}$）になります。

1.3 交換相関孔（Exchange-Correlation Hole）と総和規則

交換相関エネルギーの物理的内容を理解するための極めて強力な概念が交換相関孔（Exchange-Correlation Hole） $n_{xc}(\mathbf{r}, \mathbf{r}’)$ である。 電子間の量子力学的な多体相関により、ある位置 $\mathbf{r}$ に電子が存在するとき、その周囲の他の位置 $\mathbf{r}’$ における電子の存在確率は、単純な独立確率密度 $\rho(\mathbf{r})\rho(\mathbf{r}’)$ から減少（または変形）する。この効果を考慮した2電子確率密度は、ペア相関関数 $g(\mathbf{r}, \mathbf{r}’)$ を用いて以下のように表される。 \(\rho_2(\mathbf{r}, \mathbf{r}') = \rho(\mathbf{r})\rho(\mathbf{r}') g(\mathbf{r}, \mathbf{r}')\) $g(\mathbf{r}, \mathbf{r}’)$ の具体的な形は当然系によって異なるから，ここではなんらかの2変数の関数と考えておけば良い．一番極端な例としては，電子間に一切の相関（パウリの排他原理もクーロン斥力も）がなければ $g(\mathbf{r}, \mathbf{r}’) = 1$ となり、2体確率は単純な独立確率の積になる．$g(\mathbf{r}, \mathbf{r}’) \neq 1$ となる偏りが電子間の「交換（排他原理）」と「相関（静電排斥）」を表す．

ここで、交換相関孔 $n_{xc}(\mathbf{r}, \mathbf{r}’)$ は以下のように定義される。 \(n_{xc}(\mathbf{r}, \mathbf{r}') = \rho(\mathbf{r}') [g(\mathbf{r}, \mathbf{r}') - 1]\) このとき、電子間の全相互作用エネルギー $V_{ee}[\rho]$ は静電的な表示を用いて次のように書ける。 \(V_{ee}[\rho] = \frac{1}{2} \iint \frac{\rho_2(\mathbf{r}, \mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r} d\mathbf{r}' = E_H[\rho] + \frac{1}{2} \iint \frac{\rho(\mathbf{r}) n_{xc}(\mathbf{r}, \mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r} d\mathbf{r}'\) したがって、交換相関エネルギーのうち相互作用項（非古典静電相互作用）は、位置 $\mathbf{r}$ にある電子が、自身によって排除された電子密度の分布（すなわち交換相関孔）のつくるポテンシャルと静電相互作用しているものと解釈できる。

また、ある参照電子は自分自身を排除するため、排除される全電荷量は正確に $-1$ （電子1個分）でなければならない。これが交換相関孔の総和規則（Sum Rule）である。 \(\int n_{xc}(\mathbf{r}, \mathbf{r}') d\mathbf{r}' = -1\) この総和規則は、任意の電子密度 $\rho(\mathbf{r})$ や相互作用強度に対して厳密に成立する強固な制約条件であり、局所密度近似（LDA）をはじめとする近似汎関数が実用上頑健なエネルギーを与える背景要因となっている。

2. 局所密度近似 (LDA) とディラック交換汎関数の数学的導出

2.1 局所密度近似 (LDA) の基本概念

局所密度近似（Local Density Approximation: LDA）は、交換相関汎関数の最も基本となる近似であり、ジェイコブの梯子の第1段に相当する。LDAでは、不均一な系（分子や固体など）の局所的な空間点 $\mathbf{r}$ における交換相関エネルギー密度が、その点と同じ一定の電子密度 $\rho = \rho(\mathbf{r})$ を持つ均質電子ガス（Uniform Electron Gas: UEG）の交換相関エネルギー密度に等しいと仮定する。

全体の交換相関エネルギーは以下のように定式化される。 \(E_{xc}^{LDA}[\rho] = \int \rho(\mathbf{r}) \varepsilon_{xc}^{UEG}(\rho(\mathbf{r})) d\mathbf{r}\) ここで、$\varepsilon_{xc}^{UEG}(\rho)$ は密度 $\rho$ の均質電子ガスにおける1粒子あたりの交換相関エネルギーである。通常、これは交換項（Exchange）と相関項（Correlation）に分離して扱われる。 \(\varepsilon_{xc}^{UEG}(\rho) = \varepsilon_x^{UEG}(\rho) + \varepsilon_c^{UEG}(\rho)\)

2.2 均質電子ガスの第1次簡約密度行列

交換エネルギー密度 $\varepsilon_x^{UEG}(\rho)$ を導出するために、相互作用のない均質電子ガス（UEG）モデルの第1次簡約密度行列（first-order reduced density matrix） $\rho(\mathbf{r}_1, \mathbf{r}_2)$ を考える。

非相互作用UEGの基底状態は、全空間にわたり一定の正電荷背景の中で運動する電子の Slater 行列式で表され、一粒子軌道は以下の平面波 $\psi_{\mathbf{k}}(\mathbf{r})$ である。 \(\psi_{\mathbf{k}}(\mathbf{r}) = \frac{1}{\sqrt{V}} e^{i \mathbf{k} \cdot \mathbf{r}}\) ここで $V$ は系の体積である。各軌道はスピンの縮退度 2 を持つ。第1次簡約密度行列 $\rho(\mathbf{r}_1, \mathbf{r}_2)$ は、占有軌道（フェルミ波数 $k_F$ 以下の波動ベクトル $\mathbf{k}$）の和として定義される。 \(\rho(\mathbf{r}_1, \mathbf{r}_2) = 2 \sum_{|\mathbf{k}| \le k_F} \psi_{\mathbf{k}}^*(\mathbf{r}_1) \psi_{\mathbf{k}}(\mathbf{r}_2) = \frac{2}{V} \sum_{|\mathbf{k}| \le k_F} e^{-i \mathbf{k} \cdot (\mathbf{r}_1 - \mathbf{r}_2)}\)

熱力学的極限（$V \to \infty$）において、$\mathbf{k}$ に関する和は以下のように積分へ移行する。 \(\sum_{\mathbf{k}} \to \frac{V}{(2\pi)^3} \int d\mathbf{k}\) したがって、密度行列は次のように書き直される。 \(\rho(\mathbf{r}_1, \mathbf{r}_2) = \frac{2}{(2\pi)^3} \int_{|\mathbf{k}| \le k_F} e^{-i \mathbf{k} \cdot (\mathbf{r}_1 - \mathbf{r}_2)} d\mathbf{k}\) ここで、相対座標ベクトルを $\mathbf{s} = \mathbf{r}_1 - \mathbf{r}_2$、その大きさを $s = |\mathbf{s}|$ とする。$\mathbf{k}$ の積分を、$\mathbf{s}$ の方向を極軸にとった球座標系（$d\mathbf{k} = k^2 \sin\theta dk d\theta d\phi$、$\mathbf{k} \cdot \mathbf{s} = k s \cos\theta$）で実行する。 \(\begin{aligned} \rho(\mathbf{r}_1, \mathbf{r}_2) &= \frac{2}{(2\pi)^3} \int_0^{k_F} k^2 dk \int_0^{2\pi} d\phi \int_0^\pi e^{-i k s \cos\theta} \sin\theta d\theta \\ &= \frac{4\pi}{(2\pi)^3} \int_0^{k_F} k^2 \left[ \frac{e^{-i k s \cos\theta}}{i k s} \right]_0^\pi dk \\ &= \frac{1}{2\pi^2} \int_0^{k_F} k^2 \left( \frac{e^{i k s} - e^{-i k s}}{i k s} \right) dk \\ &= \frac{1}{\pi^2 s} \int_0^{k_F} k \sin(k s) dk \end{aligned}\)

この定積分を部分積分を用いて実行する。 \(\begin{aligned} \int_0^{k_F} k \sin(k s) dk &= \left[ -k \frac{\cos(k s)}{s} \right]_0^{k_F} + \int_0^{k_F} \frac{\cos(k s)}{s} dk \\ &= -\frac{k_F \cos(k_F s)}{s} + \frac{\sin(k_F s)}{s^2} \\ &= \frac{\sin(k_F s) - k_F s \cos(k_F s)}{s^2} \end{aligned}\) これより、第1次簡約密度行列は以下のように求まる。 \(\rho(\mathbf{r}_1, \mathbf{r}_2) = \frac{\sin(k_F s) - k_F s \cos(k_F s)}{\pi^2 s^3}\)

ここで、均質電子ガスの空間密度 $\rho$ は $\rho(\mathbf{r}, \mathbf{r})$ に等しく、上記の極限 $s \to 0$ をとるか、または単純に運動量空間の占有球の体積から以下のように与えられる。 \(\rho = \frac{2}{V} \sum_{|\mathbf{k}| \le k_F} 1 = \frac{2}{(2\pi)^3} \left( \frac{4}{3}\pi k_F^3 \right) = \frac{k_F^3}{3\pi^2}\) したがって、$\pi^2 = \frac{k_F^3}{3\rho}$ である。これを密度行列の分母に代入することで、目的の表式が得られる。 \(\rho(\mathbf{r}_1, \mathbf{r}_2) = 3\rho \left( \frac{\sin(k_F s) - k_F s \cos(k_F s)}{(k_F s)^3} \right)\) この密度行列は、非相互作用のUEGにおける電子密度の空間的コヒーレンス（重なり）を表しており、ディラック交換のすべての数学的情報を含んでいる。

2.3 ディラック交換エネルギーの導出と定積分評価

ハートリー・フォック理論より、スピン制限閉殻系における交換エネルギー $E_x$ は、第1次簡約密度行列を用いて以下のように定義される。 \(E_x = -\frac{1}{4} \iint \frac{|\rho(\mathbf{r}_1, \mathbf{r}_2)|^2}{|\mathbf{r}_1 - \mathbf{r}_2|} d\mathbf{r}_1 d\mathbf{r}_2\) ここで、係数 $-1/4$ は、同スピン間（$\uparrow\uparrow$ および $\downarrow\downarrow$）のみで交換相互作用が働くこと、および全電子密度行列 $\rho = \rho_\uparrow + \rho_\downarrow$ において $\rho_\uparrow = \rho_\downarrow = \rho/2$ であること（$\sum_{\sigma} |\rho_\sigma|^2 = 2 \times |\rho/2|^2 = |\rho|^2/2$）に由来し、ハミルトニアンの交換項に含まれる $1/2$ と組み合わさって生じる。

この重積分を実行するため、以下の座標変換を行う。

  • 相対座標： $\mathbf{s} = \mathbf{r}_1 - \mathbf{r}_2$
  • 重心座標： $\mathbf{R} = \frac{\mathbf{r}_1 + \mathbf{r}_2}{2}$

この変換のヤコビアン行列式は 1 であり、微小体積要素は $d\mathbf{r}_1 d\mathbf{r}_2 = d\mathbf{R} d\mathbf{s}$ となる。系の均一性（均質性）より、被積分関数は重心座標 $\mathbf{R}$ に依存しないため、$\mathbf{R}$ に関する積分は単に系の全体積 $V$ を与える。 \(\int d\mathbf{R} = V\) したがって、交換エネルギーは相対座標 $\mathbf{s}$ のみの積分に帰着される。 \(E_x = -\frac{V}{4} \int \frac{|\rho(\mathbf{s})|^2}{s} d\mathbf{s}\) 被積分関数は球対称であるため、3次元球座標表示（$d\mathbf{s} = s^2 \sin\theta_s ds d\theta_s d\phi_s$）を用いて角度積分を実行すると、立体角 $4\pi$ が得られる。 \(E_x = -\frac{V}{4} \cdot 4\pi \int_0^\infty \frac{|\rho(s)|^2}{s} s^2 ds = -V\pi \int_0^\infty s |\rho(s)|^2 ds\)

ここで、先に導出した $\rho(s)$ を代入する。 \(E_x = -V\pi \int_0^\infty s \left[ 9\rho^2 \frac{(\sin(k_F s) - k_F s \cos(k_F s))^2}{(k_F s)^6} \right] ds\) 無次元の積分変数 $t = k_F s$ （$s = \frac{t}{k_F}$, $ds = \frac{dt}{k_F}$）を導入して整理する。 \(\begin{aligned} E_x &= -9\pi V \rho^2 \int_0^\infty \left(\frac{t}{k_F}\right) \frac{(\sin t - t\cos t)^2}{t^6} \left(\frac{dt}{k_F}\right) \\ &= -\frac{9\pi V \rho^2}{k_F^2} \int_0^\infty \frac{(\sin t - t\cos t)^2}{t^5} dt \end{aligned}\) 1粒子あたりの交換エネルギー $\varepsilon_x = \frac{E_x}{N} = \frac{E_x}{\rho V}$ は以下のように表される。 \(\varepsilon_x = -\frac{9\pi \rho}{k_F^2} I\) ここで、$I$ は以下の定積分である。 \(I = \int_0^\infty \frac{(\sin t - t\cos t)^2}{t^5} dt\)

この定積分 $I$ を、部分積分を繰り返し用いて厳密に評価する。

ステップ 1：第1回部分積分 $u(t) = (\sin t - t\cos t)^2$, $dv(t) = t^{-5} dt$ と置く。このとき $v(t) = -\frac{1}{4}t^{-4}$ である。 $u’(t)$ を計算すると、 \(u'(t) = 2(\sin t - t\cos t)(\cos t - \cos t + t\sin t) = 2t\sin t(\sin t - t\cos t) = 2t\sin^2 t - 2t^2\sin t\cos t\) となる。部分積分公式 $\int u dv = [uv] - \int v du$ を適用する。 境界項 $[u(t)v(t)]_0^\infty$ について評価する。

  • $t \to 0$ において、$\sin t - t\cos t \approx \frac{t^3}{3}$ より $u(t) \approx \frac{t^6}{9}$ となり、 \(\lim_{t \to 0} u(t)v(t) = \lim_{t \to 0} \left( \frac{t^6}{9} \right) \left( -\frac{1}{4t^4} \right) = \lim_{t \to 0} \left( -\frac{t^2}{36} \right) = 0\)
  • $t \to \infty$ において、$(\sin t - t\cos t)^2 \le (1+t)^2 \approx t^2$ より、 \(\lim_{t \to \infty} |u(t)v(t)| \le \lim_{t \to \infty} \frac{t^2}{4t^4} = 0\) よって、境界項は 0 となる。したがって、 \(I = -\int_0^\infty \left( -\frac{1}{4t^4} \right) (2t\sin^2 t - 2t^2\sin t\cos t) dt = \frac{1}{2} \int_0^\infty \frac{\sin^2 t - t\sin t\cos t}{t^3} dt\)

ステップ 2：第2回部分積分 $J = \int_0^\infty \frac{\sin^2 t - t\sin t\cos t}{t^3} dt$ と置き、さらに部分積分を行う。 $U(t) = \sin^2 t - t\sin t\cos t$, $dV(t) = t^{-3} dt$ と置く。このとき $V(t) = -\frac{1}{2}t^{-2}$ である。 $U’(t)$ を計算すると、 \(U'(t) = 2\sin t\cos t - \left( \sin t\cos t + t(\cos^2 t - \sin^2 t) \right) = \sin t\cos t - t\cos(2t) = \frac{1}{2}\sin(2t) - t\cos(2t)\) となる。境界項 $[U(t)V(t)]_0^\infty$ を評価する。

  • $t \to 0$ において、$\sin^2 t - t\sin t\cos t \approx t^2 - t\left(t - \frac{2}{3}t^3\right) = \frac{2}{3}t^4$ より、 \(\lim_{t \to 0} U(t)V(t) = \lim_{t \to 0} \left( \frac{2}{3}t^4 \right) \left( -\frac{1}{2t^2} \right) = 0\)
  • $t \to \infty$ において、$|U(t)V(t)| \le \frac{1+t}{2t^2} \to 0$ となる。 したがって、境界項は 0 となり、次式を得る。 \(J = -\int_0^\infty \left( -\frac{1}{2t^2} \right) \left( \frac{1}{2}\sin(2t) - t\cos(2t) \right) dt = \frac{1}{4} \int_0^\infty \frac{\sin(2t) - 2t\cos(2t)}{t^2} dt\)

ステップ 3：変数変換と厳密解の導出 積分変数 $y = 2t$ （$dt = \frac{dy}{2}$）を導入する。 \(J = \frac{1}{4} \int_0^\infty \frac{\sin y - y\cos y}{(y/2)^2} \frac{dy}{2} = \frac{1}{2} \int_0^\infty \frac{\sin y - y\cos y}{y^2} dy\) ここで、被積分関数が以下の微分の形に完全に一致することに注目する。 \(\frac{d}{dy} \left( -\frac{\sin y}{y} \right) = -\frac{\cos y \cdot y - \sin y}{y^2} = \frac{\sin y - y\cos y}{y^2}\) したがって、この積分は極めて単純に実行できる。 \(\int_0^\infty \frac{\sin y - y\cos y}{y^2} dy = \left[ -\frac{\sin y}{y} \right]_0^\infty = 0 - (-1) = 1\) ここで、$\lim_{y \to 0} \frac{\sin y}{y} = 1$ を用いた。 よって、 \(J = \frac{1}{2} \times 1 = \frac{1}{2}\) 最終的に、求める定積分 $I$ の値は以下のように決定される。 \(I = \frac{1}{2} J = \frac{1}{4}\)

この結果を $\varepsilon_x$ の式に代入する。 \(\varepsilon_x = -\frac{9\pi \rho}{k_F^2} \left(\frac{1}{4}\right) = -\frac{9\pi \rho}{4 k_F^2}\) ここに $k_F = (3\pi^2\rho)^{1/3}$ を代入する。 \(\varepsilon_x^{LDA}(\rho) = -\frac{9\pi \rho}{4 (3\pi^2\rho)^{2/3}} = -\frac{9\pi}{4 (3\pi^2)^{2/3}} \rho^{1/3}\) プレファクターを整理する。 \(\frac{9\pi}{4 (3\pi^2)^{2/3}} = \frac{9\pi}{4 \cdot 3^{2/3} \pi^{4/3}} = \frac{3^{2 - 2/3}}{4 \pi^{1/3}} = \frac{3}{4} \left( \frac{3}{\pi} \right)^{1/3}\) したがって、以下のディラック交換エネルギーの表式が得られる。 \(\varepsilon_x^{LDA}(\rho) = -C_x \rho^{1/3}\) ここで、ディラック交換定数 $C_x$ は以下のように定義される。 \(C_x = \frac{3}{4} \left( \frac{3}{\pi} \right)^{1/3} \approx 0.73855877\)

スピン極性への拡張：局所スピン密度近似 (LSDA)

実際の電子系では上向きスピン（$\uparrow$）と下向きスピン（$\downarrow$）の数が異なるスピン極性状態が存在する。スピン制限のない局所スピン密度近似（Local Spin Density Approximation: LSDA）では、交換エネルギーは各スピン密度の独立な寄与の和となる。 \(E_x^{LSDA}[\rho_\uparrow, \rho_\downarrow] = E_x[\rho_\uparrow] + E_x[\rho_\downarrow]\) ここで、スピン分極率 $\zeta = \frac{\rho_\uparrow - \rho_\downarrow}{\rho}$ （ただし $\rho = \rho_\uparrow + \rho_\downarrow$）を導入すると、各スピン密度は $\rho_\uparrow = \frac{1}{2}\rho(1+\zeta)$、$\rho_\downarrow = \frac{1}{2}\rho(1-\zeta)$ と書ける。 これらを代入して整理すると、LSDAの1粒子あたり交換エネルギーは、スピン分極のないLDAの式にスピンスケーリング関数 $f(\zeta)$ を乗じた形に定式化される。 \(\varepsilon_x^{LSDA}(\rho, \zeta) = \varepsilon_x^{LDA}(\rho) f(\zeta)\) \(f(\zeta) = \frac{(1+\zeta)^{4/3} + (1-\zeta)^{4/3}}{2}\) このスケーリング関数は、完全にスピン分極した強磁性極限（$\zeta = 1$）において、交換エネルギーが $2^{1/3} \approx 1.2599$ 倍に増大することを示している。

2.4 VWN (Vosko-Wilk-Nusair) 相関汎関数の定式化

均質電子ガスの相関エネルギー $\varepsilon_c^{UEG}(\rho)$ は、交換エネルギーのように単純な解析解を得ることができない。量子モンテカルロ（QMC）法などの高精度な数値計算結果を再現するために、様々なパラメータフィッティングが行われてきた。その代表であり、実質的な業界標準となっているのが Vosko-Wilk-Nusair (VWN) 相関汎関数（具体的には1980年に提案された彼らの第5式、VWN5）である。

VWN汎関数では、まず電子密度 $\rho$ を無次元のウィグナー・ザイツ半径 $r_s$ に変換する。 \(r_s = \left( \frac{3}{4\pi\rho} \right)^{1/3}\) スピン極性のない場合（$\zeta = 0$、パラ磁性状態）および完全に極化した場合（$\zeta = 1$、強磁性状態）の各々に対して、変数 $x = \sqrt{r_s}$ を用いた以下の関数形式で相関エネルギーをフィッティングする。 \(\varepsilon_c(r_s) = \frac{A}{2} \left\{ \ln\left(\frac{x^2}{X(x)}\right) + \frac{2b}{\sqrt{4c - b^2}} \arctan\left(\frac{\sqrt{4c-b^2}}{2x + b}\right) - \frac{bx_0}{X(x_0)} \left[ \ln\left(\frac{(x-x_0)^2}{X(x)}\right) + \frac{2(b + 2x_0)}{\sqrt{4c - b^2}} \arctan\left(\frac{\sqrt{4c-b^2}}{2x + b}\right) \right] \right\}\) ここで、$X(x) = x^2 + b x + c$ である。フィッティングパラメータ $A$, $x_0$, $b$, $c$ の値は以下の通りである。

状態 $A$ $x_0$ $b$ $c$
パラ磁性 ($\zeta=0$) $0.0621814$ $-0.10498$ $3.72744$ $12.9352$
強磁性 ($\zeta=1$) $0.0310907$ $-0.32500$ $7.06042$ $18.0578$

任意のスピン分極率 $\zeta$ に対する相関エネルギー $\varepsilon_c(r_s, \zeta)$ は、以下の補間式を用いて算出される。 \(\varepsilon_c(r_s, \zeta) = \varepsilon_c^P(r_s) + f_c(\zeta) \left[ \varepsilon_c^F(r_s) - \varepsilon_c^P(r_s) \right] (1 - \zeta^4) + \dots\) あるいは、スピン剛性（spin stiffness）$\alpha_c(r_s)$ を含めた高精度な補間式（VWN5の真の形式）が用いられる。 \(\varepsilon_c(r_s, \zeta) = \varepsilon_c^P(r_s) + \alpha_c(r_s) \frac{f(\zeta)}{f''(0)} \left(1 - \zeta^4\right) + \left[\varepsilon_c^F(r_s) - \varepsilon_c^P(r_s)\right] f(\zeta) \zeta^4\) ここで、スピン補間関数 $f(\zeta)$ は次のように定義される。 \(f(\zeta) = \frac{(1+\zeta)^{4/3} + (1-\zeta)^{4/3} - 2}{2^{4/3} - 2}\) また、$\alpha_c(r_s)$ はスピン剛性であり、$\varepsilon_c$ と同様の関数形式でパラメータフィッティングされている。VWN5はこのように理論的漸近挙動とQMCデータを極めて精密に繋ぐことで、実用上非常に信頼性の高い相関エネルギーのベースラインを提供している。

3. 一般化勾配近似 (GGA) と PBE 汎関数の理論とパラメータ設計

3.1 一般化勾配近似 (GGA) の基本概念と無次元還元密度勾配 $s$

局所密度近似（LDA）は電子密度が緩やかに変化する均質電子ガスに近い系では良好な結果を与えるが、実際の分子や表面のように電子密度が空間的に激しく変化する不均一な系では、結合エネルギーの過大評価（過大結合）などの顕著な問題が生じる。この密度の「不均一性」を取り入れるために、各空間点における電子密度 $\rho(\mathbf{r})$ に加えて、密度の一次勾配 $\nabla \rho(\mathbf{r})$ を局所的な変数として取り入れたのが一般化勾配近似（Generalized Gradient Approximation: GGA）である。

GGAの交換相関エネルギーは一般に以下の形式で定式化される。 \(E_{xc}^{GGA}[\rho] = \int \rho(\mathbf{r}) \varepsilon_x^{LDA}(\rho(\mathbf{r})) F_{xc}(\rho(\mathbf{r}), \nabla \rho(\mathbf{r})) d\mathbf{r}\) ここで、$\varepsilon_x^{LDA}(\rho) = -C_x \rho^{1/3}$ はLDAの交換エネルギー密度であり、$F_{xc}$ は増強因子（Enhancement Factor）と呼ばれる無次元の関数である。増強因子は、不均一性によるLDAからのズレ（補正）を表現する。

密度の不均一性を定量化する無次元の指標として、無次元還元密度勾配（dimensionless reduced density gradient） $s(\mathbf{r})$ が導入される。 \(s(\mathbf{r}) = \frac{|\nabla \rho(\mathbf{r})|}{2 k_F(\mathbf{r}) \rho(\mathbf{r})}\) ここで、$k_F(\mathbf{r}) = (3\pi^2 \rho(\mathbf{r}))^{1/3}$ は局所的なフェルミ波数である。この $k_F$ を用いて分母を書き下すと、次のように整理される。 \(s(\mathbf{r}) = \frac{|\nabla \rho(\mathbf{r})|}{2 (3\pi^2)^{1/3} \rho(\mathbf{r})^{4/3}} \approx \frac{|\nabla \rho(\mathbf{r})|}{6.187 \rho(\mathbf{r})^{4/3}}\) 還元勾配 $s$ は、密度の空間変化の特性長さである $|\nabla \rho|/\rho$ と、局所的なフェルミ波長 $\lambda_F = 2\pi/k_F$ （または平均電子間距離 $r_s \sim 1/k_F$）の比と見なすことができる。

  • $s = 0$ は、空間変化が一切ない均質電子ガス（UEG）の限界に対応し、増強因子は $F_{xc} \to 1$ （または相関部分を含めた対応値）となってLDAに帰着する。
  • $s$ が大きくなるにつれ、密度の空間的変化が激しく、より不均一な系であることを示す。

3.2 PBE 交換汎関数の定式化と増強因子 $F_x^{PBE}(s)$

1996年にPerdew、Burke、Ernzerhofによって提案されたPBE汎関数は、実験データを用いたパラメータフィッティングを排除し、厳密な物理条件のみを満たすように構築された非経験的（ab initio）なGGAである。

PBEの交換増強因子 $F_x^{PBE}(s)$ は以下の数学的表式で与えられる。 \(F_x^{PBE}(s) = 1 + \kappa - \frac{\kappa}{1 + \frac{\mu}{\kappa} s^2}\) この表式は、還元密度勾配 $s$ の値に応じて以下のように極限挙動を制御する。

  1. 緩やかに変化する限界（$s \to 0$）： テイラー展開を行うことで、次のように $s^2$ の2次の項で近似される。 \(F_x^{PBE}(s) \approx 1 + \mu s^2 + \mathcal{O}(s^4)\) ここで、$\mu$ は $s^2$ に比例する補正項の初期の立ち上がりを決定するパラメータである。
  2. 急激に変化する限界（$s \to \infty$）： 分母の $s^2$ の項が無限大に発散するため、増強因子は一定の上限値に漸近する。 \(\lim_{s \to \infty} F_x^{PBE}(s) = 1 + \kappa\) ここで、$\kappa$ は増強因子の最大値を支配するパラメータである。

3.3 PBE パラメータ $\kappa$ と Lieb-Oxford 限界の物理的導出

パラメータ $\kappa$ の値は、量子力学的な全交換相関エネルギー（および交換エネルギー）の上限を定めるLieb-Oxford限界（Lieb-Oxford bound）を厳密に満たすように決定される。

1981年にLiebとOxfordによって証明された定理によると、任意の電子密度 $\rho(\mathbf{r})$ を持つ多電子系において、負の交換相関エネルギー $E_{xc}[\rho]$ の大きさは、以下の不等式によって制限される。 \(E_{xc}[\rho] \ge E_x[\rho] \ge -C_{LO} \int \rho(\mathbf{r})^{4/3} d\mathbf{r}\) ここで、厳密な上限定数として $C_{LO} \le 1.68$ （Hartree原子単位系）が与えられている。

この不等式を、スピン分極のないLDAの交換エネルギー \(E_x^{LDA}[\rho] = -C_x \int \rho^{4/3} d\mathbf{r} \quad \left( C_x = \frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3} \approx 0.7386 \right)\) を用いて書き直すと、次のようになる。 \(E_{xc}[\rho] \ge 2.275 E_x^{LDA}[\rho] \quad \left( \because \frac{1.68}{0.7386} \approx 2.275 \right)\)

ここで、PBEは相関エネルギー $E_c[\rho]$ を含む全交換相関エネルギーにおいてLieb-Oxford限界を満足する必要がある。一般に交換エネルギーが支配的である極限を考えると、交換エネルギー単独においてもこの制約を満たさなければならない。したがって、任意の局所点において交換エネルギー密度に対する増強因子 $F_x(s)$ は、以下の条件を満たす必要がある。 \(F_x(s) \le 2.275\)

しかし、PBEではさらに厳しいスピン極性状態（完全にスピン分極した1電子系など）におけるLieb-Oxford限界も同時に満たすように設計されている。 スピンスケーリング関係式より、完全にスピン分極した系（$\zeta = 1$）における交換エネルギーは次のように表される。 \(E_x[\rho_\uparrow, 0] = \frac{1}{2} E_x[2\rho_\uparrow] = -2^{1/3} C_x \int \rho_\uparrow^{4/3} F_x(2^{-1/3} s) d\mathbf{r}\) ここで、$2^{1/3} \approx 1.2599$ である。このスピン分極した状態に対しても Lieb-Oxford 不等式 $E_x \ge -1.68 \int \rho^{4/3} d\mathbf{r}$ が成立するためには、交換増強因子 $F_x(s)$ は以下の限界を満たさなければならない。 \(2^{1/3} F_x(s) \le 2.275 \implies F_x(s) \le \frac{2.275}{2^{1/3}} \approx 1.804\)

PBEでは、この最も厳しい限界値を漸近限界として採用している。 \(1 + \kappa = 1.804 \implies \kappa = 0.804\) これにより、PBE交換汎関数は任意の電子密度分布および任意のスピン極性において、Lieb-Oxford限界を決して破らないことが数学的に保証される。

3.4 PBE パラメータ $\mu$ と均質電子ガスの線形応答

もう一つのパラメータ $\mu = 0.21951$ は、均質電子ガス（UEG）における線形応答理論（linear response theory）の漸近挙動と、交換・相関の gradient expansion 間の物理的整合性から決定される。

不均一性が極めて小さい緩やかな変化の限界（$s \to 0$）において、交換相関エネルギー密度の勾配補正項は、以下のような2次の勾配拡張（gradient expansion approximation: GEA）で表される。 \(E_{xc}^{GEA}[\rho] = E_{xc}^{LDA}[\rho] + \int \left[ C_x \sigma_x \frac{|\nabla \rho|^2}{\rho^{4/3}} + C_c \sigma_c(\rho) \frac{|\nabla \rho|^2}{\rho^{4/3}} \right] d\mathbf{r}\) ここで、交換項の2次勾配膨張係数 $\sigma_x$ は、均質電子ガスの線形応答理論（より具体的にはスピン分極のない一様電子ガスの誘電関数 $\varepsilon(q)$ の小波数極限）から厳密に決定され、以下のように与えられる。 \(\sigma_x = \frac{10}{81} \approx 0.12346\)

しかし、PBEでは交換パラメータ $\mu$ を単独の交換GEA係数である $10/81$ に設定するのではなく、交換と相関の勾配補正項の相殺効果を考慮して決定している。 PBE相関汎関数における勾配依存項の立ち上がりを決めるパラメータ $\beta$ は、高密度極限での勾配補正を正しく再現するよう $\beta \approx 0.066725$ に定められている。 高密度極限（$r_s \to 0$）において、相関の勾配修正項と交換の勾配修正項が、線形応答理論に基づいて整合性を保つためには、交換増強因子の立ち上がりパラメータ $\mu$ は以下の関係式を満たさなければならない。 \(\mu = \frac{\pi^2 \beta}{3}\) この関係式に $\beta \approx 0.066725$ を代入することで、パラメータ $\mu$ の値が導出される。 \(\mu = \frac{\pi^2 \times 0.066725}{3} \approx 0.21951\) この設計により、PBEは緩やかに変化する限界において、相関汎関数と交換汎関数の gradient correction が線形応答の性質を保ちつつ調和して機能し、高密度限界での正しい全エネルギー応答を与える。 （注：固体物理向けに設計された PBEsol 汎関数では、固体の格子定数や表面エネルギーの再現性を向上させるために、相関パラメータを $\beta = 0.046$ とし、交換パラメータを厳密な交換単独のGEA限界値である $\mu = 10/81 \approx 0.1235$ に戻している）。

3.5 経験的 GGA の代表例：Becke の B88 汎関数との比較

PBE汎関数の物理的な整合性と対比されるのが、量子化学分野で絶大な成功を収めた経験的（empirical）なGGAである Becke’s B88交換汎関数 である。

1988年にA. D. Beckeによって提案されたB88交換エネルギー汎関数は、以下の形式を持つ。 \(E_x^{B88} = E_x^{LDA} - \beta \int \rho^{4/3} \frac{x^2}{1 + 6 \beta x \sinh^{-1} x} d\mathbf{r}\) ここで、$x$ はBeckeが定義した無次元の密度勾配変数であり、PBEの還元密度勾配 $s$ との間には次の関係がある。 \(x = \frac{|\nabla \rho|}{\rho^{4/3}} = 2 (3\pi^2)^{1/3} s \approx 6.187 s\) また、プレファクターの $\beta$ は、希ガス原子（He, Ne, Ar, Kr, Xe）のハートリー・フォック交換エネルギーの厳密値に対してフィッティングされた経験的パラメータであり、その値は $\beta = 0.0042 \text{ a.u.}$ である。

B88汎関数とPBE汎関数は、設計思想と極限挙動において以下の明確なコントラストを示す。

  1. パラメータ決定の思想：
    • PBE：実験値や厳密な原子構造のエネルギーへのフィッティングを一切行わず、厳密な物理制約（Lieb-Oxford限界、線形応答）のみからパラメータ（$\kappa, \mu$）を決定している。
    • B88：少数の物理制約を満足させつつ、実際の原子の参照データにフィッティングして $\beta$ を決定しており、「化学的精度（結合エネルギーの精密記述）」の向上に重きを置いている。
  2. 大勾配限界（$s, x \to \infty$）での挙動：
    • PBE：$s \to \infty$ において $F_x(s) \to 1.804$ となり、Lieb-Oxford限界を常に満たす。
    • B88：$x \to \infty$ において、$\sinh^{-1} x \approx \ln(2x)$ となるため、補正項は以下のように振る舞う。 \(\Delta E_x^{B88} \approx - \beta \int \rho^{4/3} \frac{x^2}{6 \beta x \ln(2x)} d\mathbf{r} = - \frac{1}{6} \int \rho^{4/3} \frac{x}{\ln(2x)} d\mathbf{r}\) これは $x \to \infty$ で発散し、非常に大きい勾配を持つ系（あるいは電子密度が極限まで希薄な領域）において、Lieb-Oxford限界を破る可能性がある。
  3. 小勾配限界（$s, x \to 0$）での挙動：
    • PBE：一様電子ガスの誘電応答と整合する $\mu = 0.2195$ を採用している。
    • B88：$x \to 0$ で $E_x^{B88} \approx E_x^{LDA} - \beta \int \rho^{4/3} x^2 d\mathbf{r}$ となり、勾配補正の係数が $\beta = 0.0042$ となるが、これは誘電応答理論から得られる厳密な値とは一致しない。

B88は分子の結合エネルギーの予測に優れる一方、PBEはその物理的頑健性から、特に固体の記述や多様な不均一電子系において優れた汎用性と信頼性を示す。

4. メタ一般化勾配近似 (meta-GGA) と SCAN 汎関数における物理制約の充足

4.1 メタ一般化勾配近似 (meta-GGA) の基本設計と運動エネルギー密度 $\tau$

一般化勾配近似（GGA）は電子密度の空間的勾配 $\nabla \rho$ を導入することで、LDAから大きな進歩を遂げた。しかし、密度とその勾配情報だけでは、多様な化学結合（共有結合、イオン結合、金属結合、および弱い非共有結合相互作用）の極所的な違いを完全に識別し、それぞれに適した相関効果を自律的に記述することには限界があった。

この限界を克服するために、ジェイコブの梯子の第3段であるメタ一般化勾配近似（meta-GGA）では、電子密度 $\rho(\mathbf{r})$ とその勾配 $\nabla \rho(\mathbf{r})$ に加え、Kohn-Sham占有軌道の勾配から定義される正値の局所軌道運動エネルギー密度（positive-definite local kinetic energy density） $\tau(\mathbf{r})$ を追加の局所変数として導入する。 \(\tau(\mathbf{r}) = \frac{1}{2} \sum_{i=1}^{occ} |\nabla \phi_i(\mathbf{r})|^2\) ここで、$\phi_i(\mathbf{r})$ はKohn-Sham占有軌道の波動関数である。運動エネルギー密度 $\tau(\mathbf{r})$ は軌道の情報を陽に含むため、局所的な電子密度の重なりや、パウリの排他原理による電子排除の度合い（軌道の局在化状態）を非常に敏感に感知することができる。

4.2 Weizsäcker 運動エネルギー密度 $\tau_W$ と von Weizsäcker 限界

運動エネルギー密度 $\tau(\mathbf{r})$ の解釈および局所的な物理情報の抽出において極めて重要な役割を果たすのが、von Weizsäcker 運動エネルギー密度（von Weizsäcker kinetic energy density） $\tau_W(\mathbf{r})$ である。 \(\tau_W(\mathbf{r}) = \frac{|\nabla \rho(\mathbf{r})|^2}{8 \rho(\mathbf{r})}\) この $\tau_W(\mathbf{r})$ は、すべての電子が単一の軌道（ボース・アインシュタイン凝縮状態、あるいはスピン二重占有された単一の空間軌道）に収容されている仮想的な状態における運動エネルギー密度に正確に一致する。

量子力学における一般定理（Weizsäcker の不等式）により、任意の多電子系において、全運動エネルギー密度 $\tau(\mathbf{r})$ は、von Weizsäcker 運動エネルギー密度 $\tau_W(\mathbf{r})$ を下限とすることが数学的に証明されている。 \(\tau(\mathbf{r}) \ge \tau_W(\mathbf{r})\) この不等式は、軌道運動エネルギー密度 $\tau(\mathbf{r})$ が von Weizsäcker 限界 $\tau_W$ を下回ることができないことを意味し、両者の差 $\tau(\mathbf{r}) - \tau_W(\mathbf{r})$ はパウリの排他原理に起因する追加の運動エネルギー（Pauli運動エネルギー）を表している。

4.3 軌道重なり指標 $\alpha$ と化学環境の判別

meta-GGA汎関数では、上記の $\tau(\mathbf{r})$ と $\tau_W(\mathbf{r})$ を用いて、空間の各点における軌道の重なり合い（コヒーレンス）を特徴付ける軌道重なり指標（orbital overlap indicator） $\alpha(\mathbf{r})$ を以下のように定義する。 \(\alpha(\mathbf{r}) = \frac{\tau(\mathbf{r}) - \tau_W(\mathbf{r})}{\tau_{unif}(\mathbf{r})}\) ここで、分母の $\tau_{unif}(\mathbf{r})$ は、その点における局所密度 $\rho(\mathbf{r})$ と同じ密度を持つ均質電子ガス（一様電子ガス）の運動エネルギー密度（トーマス・フェルミ運動エネルギー密度）であり、以下のように定義される。 \(\tau_{unif}(\mathbf{r}) = \frac{3}{10} (3\pi^2)^{2/3} \rho(\mathbf{r})^{5/3}\)

軌道重なり指標 $\alpha(\mathbf{r})$ は、Weizsäcker限界の成立より常に $\alpha(\mathbf{r}) \ge 0$ を満たし、その値は局所的な化学結合の種類や物理的な状態（化学環境）を以下のように極めて精密に特徴付ける。

  1. $\alpha(\mathbf{r}) \approx 0$ （単一軌道支配領域、共有結合・孤立電子対）： $\tau \approx \tau_W$ となるのは、局所的な密度が単一の空間軌道（または互いにスピンが逆向きの2電子）によって形成されている領域である。
    • 具体例：共有結合の中間点（例：$H_2$, $N_2$ の結合中心）、原子の孤立電子対（lone pair）領域、あるいは希ガス原子の最外殻電子殻の境界。
    • 物理的意味：パウリ運動エネルギーがゼロであり、電子が局在化して強い化学結合や孤立電子対を形成している。
  2. $\alpha(\mathbf{r}) \approx 1$ （ゆっくり変化する密度領域、金属結合・均質流）： $\tau \approx \tau_{unif}$ かつ $\tau_W \ll \tau$ となる領域である。
    • 具体例：単純金属（Na, Alなど）の結晶内部、あるいは均質電子ガスに近い自由電子的な振る舞いを示す領域。
    • 物理的意味：電子軌道が空間的に広く非局在化し、多方向に重なり合って均一な電子の流れを作っている。
  3. $\alpha(\mathbf{r}) \gg 1$ （軌道の重なりが極めて弱い領域、非共有結合・分散力支配領域）： 電子密度が空間的に重なり合っていないが、波動関数のテールの重なりが存在する領域である。
    • 具体例：希ガス原子間の隙間、ファンデルワールス（vdW）力で結合した分子間の界面、または原子・分子の外側に広がる希薄なテール領域。
    • 物理的意味：電子密度 $\rho$ が極めて小さいため、分母の $\tau_{unif} \propto \rho^{5/3}$ が分子の $\tau - \tau_W$ よりもはるかに早くゼロに漸近するため、比である $\alpha$ が非常に大きくなる。

meta-GGA汎関数は、この局所的な $\alpha(\mathbf{r})$ の値をリアルタイムに参照することで、現在計算している空間点が「共有結合領域 ($\alpha \approx 0$)」なのか、「金属的な領域 ($\alpha \approx 1$)」なのか、あるいは「弱い相互作用の領域 ($\alpha \gg 1$)」なのかを自動的に判別し、それぞれの物理的状態に適した交換相関補正を動的に切り替えることができる。これがmeta-GGAがGGAを大きく凌駕する予測精度を持つ最大の理由である。

4.4 SCAN 汎関数の理論的設計と 17 の厳密な物理制約の完全詳報

2015年にJ. Sun、A. Ruzsinszky、J. P. Perdewによって提案されたSCAN（Strongly Constrained and Appropriately Normed）汎関数は、非経験的メタGGAの極致である。SCANは、半局所的な汎関数（密度、密度勾配、軌道運動エネルギー密度のみに依存する汎関数）が満たし得る、数学的および量子力学的に既知の17の厳密な物理的制約（exact constraints）をすべて満足するように数式モデルが数学的に構築されている。

以下に、SCAN汎関数が満足する17の厳密な物理的制約について、それぞれの物理的な定義と意義を詳細に列挙する。

I. 交換汎関数が満たす 6 つの厳密な制約

  1. 交換エネルギーの非正性（Negativity）： 任意の電子密度に対して、交換エネルギーは常に負またはゼロである。 \(E_x[\rho] \le 0\)
  2. スピンスケーリング関係式（Spin-scaling）： スピン偏極した系の交換エネルギーは、スピン非偏極系の交換エネルギーから厳密に変換可能である。 \(E_x[\rho_\uparrow, \rho_\downarrow] = \frac{1}{2} \left( E_x[2\rho_\uparrow] + E_x[2\rho_\downarrow] \right)\)
  3. 一様座標スケーリング限界（Uniform Coordinate Scaling）： 密度を一様に縮小・拡大させたとき（$\rho_\lambda(\mathbf{r}) = \lambda^3 \rho(\lambda \mathbf{r})$）、交換エネルギーは $\lambda$ に対して一次の均一性を持つ。 \(E_x[\rho_\lambda] = \lambda E_x[\rho]\)
  4. 緩やかに変化する限界における4次までの勾配拡張（Fourth-Order Gradient Expansion）： 空間変化が極めて小さいリミットにおいて、交換エネルギーは2次および4次の勾配拡張（GEA）を正しく回復する。
  5. 非一様スケーリングにおける大次元リミット（Non-Uniform Scaling in 1D）： 密度を1つの次元方向（例えば $z$ 方向）にのみ強く圧縮した場合（$\rho_\lambda(x, y, z) = \lambda \rho(x, y, \lambda z)$）、極限 $\lambda \to \infty$ において交換エネルギーは発散せず、有限な値にとどまる。 \(\lim_{\lambda \to \infty} E_x[\rho_\lambda] > -\infty\)
  6. 2電子系の交換に対する厳密な下限（Tight bound for two-electron densities）： 2電子系において、局所的な交換増強因子は $F_x(s) \le 1.174$ を満たし、一般系よりも厳しい下限を破らない。

II. 相関汎関数が満たす 6 つの厳密な制約

  1. 相関エネルギーの非陽性（Non-positivity）： 任意の電子密度に対して、相関による寄与は静電反発を緩和する効果であるため、常に負またはゼロである。 \(E_c[\rho] \le 0\)
  2. 緩やかに変化する限界における2次の勾配拡張（Second-Order Gradient Expansion）： $s \to 0$ において、相関エネルギーは半古典的勾配拡張の2次項を正しく再現する。
  3. 一様座標スケーリングの高密度極限（High-Density Limit of Uniform Scaling）： 一様に高密度化させた極限（$\lambda \to \infty$）において、相関エネルギーは有限な定数（スケーリングに依存しない一定値）に収束する。 \(\lim_{\lambda \to \infty} E_c[\rho_\lambda] = \text{const} > -\infty\)
  4. 一様座標スケーリングの低密度極限（Low-Density Limit of Uniform Scaling）： 一様かつ無限に希薄化させた極限（$\lambda \to 0$）において、相関エネルギーはスケーリングパラメータの $1/2$ 乗に比例して減少する。 \(E_c[\rho_\lambda] \propto \lambda^{1/2} \quad (\lambda \to 0)\)
  5. 1電子系の相関エネルギーの完全消失（Zero Correlation for One-Electron Densities）： 電子が1個しか存在しない系（スピン分極率 $\zeta = 1$、かつ単一占有）では、自分自身との動的相関は存在しないため、相関エネルギーは厳密にゼロにならなければならない。 \(E_c[\rho_\uparrow, 0] = 0\) SCANはこの条件を $\alpha \to 0$ の限界で完全に満足し、自己相関誤差（Self-Correlation Error）を本質的に排除している。
  6. 非一様スケーリングにおける大次元リミット（Non-Uniform Scaling for Correlation）： 1次元非一様スケーリング（$\rho_\lambda(x, y, z) = \lambda \rho(x, y, \lambda z)$）のもとで、$\lambda \to \infty$ において相関エネルギーは有限値に収束する。

III. 交換・相関の全体が満たす 5 つの厳密な制約

  1. サイズ延属性（Size Extensivity）： 空間的に無限に離れた2つの独立な系 $A$ と $B$ に対し、全交換相関エネルギーは個々の和に等しい。 \(E_{xc}[\rho_A + \rho_B] = E_{xc}[\rho_A] + E_{xc}[\rho_B]\)
  2. 一般のLieb-Oxford限界（General Lieb-Oxford Bound）： 任意の電子密度分布 $\rho(\mathbf{r})$ に対し、全交換相関エネルギーは以下の下限を下回らない。 \(E_{xc}[\rho] \ge -1.68 \int \rho(\mathbf{r})^{4/3} d\mathbf{r}\)
  3. 低密度極限における相対スピン分極依存性の弱化（Weak Spin-Polarization Dependence in Low-Density Limit）： 電子密度が極限まで希薄な場合（Wigner結晶相など）、交換相関エネルギーの相対スピン分極率 $\zeta$ に対する依存性は極めて微弱になる。
  4. 均質電子ガスの静的線形応答（Static Linear Response of Uniform Electron Gas）： 一様な密度の電子ガスに対する微小なポテンシャル摂動に対する応答（Lindhard誘電関数）の振る舞いを正しく再現する。
  5. 2電子系におけるLieb-Oxford限界の厳密充足（Lieb-Oxford Bound for Two-Electron Densities）： 2電子系において、一般系よりもさらに厳しい Lieb-Oxford bound 条件を満足する。

SCAN汎関数は、これらの17の物理的制約をすべて数式モデルに組み込むことで、経験的なパラメータフィッティングなしに極めて高い「普遍的記述力」を達成している。さらに、適切なノルム（Appropriate Norms）として、一様電子ガスや希ガス原子といった、半局所汎関数が厳密に記述すべき代表的な参照系（Normed systems）をいくつか選択し、そのエネルギーを再現するように関数の形状を決定している。

この結果、SCAN汎関数は結合エネルギーや分子・固体の幾何構造、さらにはファンデルワールス力に起因する長距離相互作用（中距離までの分散力効果）に至るまで、従来のGGAを大きく上回り、計算コストがはるかに高いハイブリッド汎関数に匹敵する精度を達成している。

5. ハイブリッド汎関数と断熱結合公式 (ACF) の詳細導出

5.1 断熱結合公式 (Adiabatic Connection Formula: ACF) の数学的導出

ハイブリッド汎関数は、Kohn-Sham密度汎関数理論（KS-DFT）の交換エネルギーの一部を、Hartree-Fock（HF）理論における「厳密な交換エネルギー（exact exchange）」で置き換える手法であり、ジェイコブの梯子の第4段に位置付けられる。このハイブリッド化の数学的および物理的正当性を保証するのが断熱結合公式（Adiabatic Connection Formula: ACF）である。

ACFは、電子密度 $\rho(\mathbf{r})$ を全プロセスで一定に保ったまま、電子間相互作用の強さを表す結合定数 $\lambda$ を $0$（非相互作用のKohn-Sham参照系）から $1$（実際の相互作用する実系）まで連続的に変化させる仮想的な断熱変化経路を考える。

この結合定数 $\lambda \in [0, 1]$ に対する補間ハミルトニアン $\hat{H}_\lambda$ は以下のように定義される。 \(\hat{H}_\lambda = \hat{T} + \hat{V}_{\text{ext}}^\lambda + \lambda \hat{V}_{ee}\) ここで、

  • $\hat{T} = -\frac{1}{2} \sum_{i=1}^N \nabla_i^2$ は運動エネルギー演算子である。
  • $\hat{V}{ee} = \sum{i < j} \frac{1}{ \mathbf{r}_i - \mathbf{r}_j }$ は電子間静電反発演算子である。
  • $\hat{V}{\text{ext}}^\lambda = \sum{i=1}^N v_{\text{ext}}^\lambda(\mathbf{r}_i)$ は、任意の $\lambda$ において基底状態の電子密度が常に実系の密度 $\rho(\mathbf{r})$ と一致するように動的に調整される局所的な一体系外部ポテンシャル演算子である。

$\hat{H}\lambda$ の基底状態波動関数を $\Psi\lambda$ とすると、基底状態の全エネルギー $E_\lambda$ は以下のように与えられる。 \(E_\lambda = \langle \Psi_\lambda | \hat{H}_\lambda | \Psi_\lambda \rangle\)

ここで、波動関数 $\Psi_\lambda$ に対するヘルマン-ファインマン（Hellmann-Feynman）の定理を適用する。基底状態波動関数の規格化条件 $\langle \Psi_\lambda \Psi_\lambda \rangle = 1$ より、エネルギー $E_\lambda$ の $\lambda$ に対する全微分は、ハミルトニアンの $\lambda$ に関する偏微分の期待値に等しくなる。  
$$\frac{d E_\lambda}{d \lambda} = \left\langle \Psi_\lambda \middle \frac{\partial \hat{H}_\lambda}{\partial \lambda} \middle \Psi_\lambda \right\rangle$$

ハミルトニアン $\hat{H}_\lambda$ の定義式を $\lambda$ で偏微分すると、運動エネルギー演算子 $\hat{T}$ は $\lambda$ に依存しないため消去され、以下を得る。 \(\frac{\partial \hat{H}_\lambda}{\partial \lambda} = \frac{\partial \hat{V}_{\text{ext}}^\lambda}{\partial \lambda} + \hat{V}_{ee}\) これを期待値の式に代入すると、以下のようになる。 \(\frac{d E_\lambda}{d \lambda} = \langle \Psi_\lambda | \hat{V}_{ee} | \Psi_\lambda \rangle + \left\langle \Psi_\lambda \middle| \frac{\partial \hat{V}_{\text{ext}}^\lambda}{\partial \lambda} \middle| \Psi_\lambda \right\rangle\)

一体系外部ポテンシャル演算子 $\hat{V}_{\text{ext}}^\lambda$ の期待値は、任意の $\lambda$ において電子密度 $\rho(\mathbf{r})$ が不変であるという拘束条件を用いることで、次のように空間積分として表すことができる。 \(\langle \Psi_\lambda | \hat{V}_{\text{ext}}^\lambda | \Psi_\lambda \rangle = \int \rho(\mathbf{r}) v_{\text{ext}}^\lambda(\mathbf{r}) d\mathbf{r}\) したがって、この項の $\lambda$ 偏微分の期待値は以下のようになる。 \(\left\langle \Psi_\lambda \middle| \frac{\partial \hat{V}_{\text{ext}}^\lambda}{\partial \lambda} \middle| \Psi_\lambda \right\rangle = \int \rho(\mathbf{r}) \frac{\partial v_{\text{ext}}^\lambda(\mathbf{r})}{\partial \lambda} d\mathbf{r}\) これより、エネルギー微分の式は以下のように書き換えられる。 \(\frac{d E_\lambda}{d \lambda} = \langle \Psi_\lambda | \hat{V}_{ee} | \Psi_\lambda \rangle + \int \rho(\mathbf{r}) \frac{\partial v_{\text{ext}}^\lambda(\mathbf{r})}{\partial \lambda} d\mathbf{r}\)

この式を $\lambda = 0$ から $\lambda = 1$ まで積分する。 \(E_{\lambda=1} - E_{\lambda=0} = \int_0^1 \frac{d E_\lambda}{d \lambda} d\lambda = \int_0^1 \langle \Psi_\lambda | \hat{V}_{ee} | \Psi_\lambda \rangle d\lambda + \int_0^1 d\lambda \int \rho(\mathbf{r}) \frac{\partial v_{\text{ext}}^\lambda(\mathbf{r})}{\partial \lambda} d\mathbf{r}\)

右辺第2項について、$\lambda$ に関する積分と空間積分 $\mathbf{r}$ の順序を入れ替えて評価する。 \(\int_0^1 d\lambda \int \rho(\mathbf{r}) \frac{\partial v_{\text{ext}}^\lambda(\mathbf{r})}{\partial \lambda} d\mathbf{r} = \int \rho(\mathbf{r}) \left( \int_0^1 \frac{\partial v_{\text{ext}}^\lambda(\mathbf{r})}{\partial \lambda} d\lambda \right) d\mathbf{r} = \int \rho(\mathbf{r}) \left[ v_{\text{ext}}^{\lambda=1}(\mathbf{r}) - v_{\text{ext}}^{\lambda=0}(\mathbf{r}) \right] d\mathbf{r}\) ここで、境界条件として：

  • $\lambda = 1$ では実系であるため、外部ポテンシャルは実際の核静電ポテンシャル $v_{\text{ext}}(\mathbf{r})$ に等しい：$v_{\text{ext}}^{\lambda=1}(\mathbf{r}) = v_{\text{ext}}(\mathbf{r})$
  • $\lambda = 0$ では非相互作用Kohn-Sham系であるため、外部ポテンシャルは有効ポテンシャル $v_{\text{eff}}(\mathbf{r})$ に等しい：$v_{\text{ext}}^{\lambda=0}(\mathbf{r}) = v_{\text{eff}}(\mathbf{r})$

したがって、 \(\int_0^1 d\lambda \int \rho(\mathbf{r}) \frac{\partial v_{\text{ext}}^\lambda(\mathbf{r})}{\partial \lambda} d\mathbf{r} = \int \rho(\mathbf{r}) v_{\text{ext}}(\mathbf{r}) d\mathbf{r} - \int \rho(\mathbf{r}) v_{\text{eff}}(\mathbf{r}) d\mathbf{r}\) となる。

次に、左辺のエネルギー差 $E_{\lambda=1} - E_{\lambda=0}$ を評価する。

  • $\lambda = 1$ の実系において、全エネルギーは以下のように定義される。 \(E_{\lambda=1} = \langle \Psi_1 | \hat{T} + \hat{V}_{\text{ext}} + \hat{V}_{ee} | \Psi_1 \rangle = T + V_{\text{ext}} + V_{ee}\) Kohn-Sham形式において、これは以下のように分割して表記される。 \(E_{\lambda=1} = T_s[\rho] + \int \rho(\mathbf{r}) v_{\text{ext}}(\mathbf{r}) d\mathbf{r} + E_H[\rho] + E_{xc}[\rho]\)
  • $\lambda = 0$ の非相互作用系において、全エネルギーは以下のように表される。 \(E_{\lambda=0} = \langle \Psi_0 | \hat{T} + \hat{V}_{\text{eff}} | \Psi_0 \rangle = T_s[\rho] + \int \rho(\mathbf{r}) v_{\text{eff}}(\mathbf{r}) d\mathbf{r}\)

これらを用いて、エネルギー差は次のように書き下せる。 \(E_{\lambda=1} - E_{\lambda=0} = \int \rho(\mathbf{r}) v_{\text{ext}}(\mathbf{r}) d\mathbf{r} - \int \rho(\mathbf{r}) v_{\text{eff}}(\mathbf{r}) d\mathbf{r} + E_H[\rho] + E_{xc}[\rho]\)

このエネルギー差の式を、積分の評価式に代入する。 \(\int \rho(\mathbf{r}) v_{\text{ext}}(\mathbf{r}) d\mathbf{r} - \int \rho(\mathbf{r}) v_{\text{eff}}(\mathbf{r}) d\mathbf{r} + E_H[\rho] + E_{xc}[\rho] = \int_0^1 \langle \Psi_\lambda | \hat{V}_{ee} | \Psi_\lambda \rangle d\lambda + \int \rho(\mathbf{r}) v_{\text{ext}}(\mathbf{r}) d\mathbf{r} - \int \rho(\mathbf{r}) v_{\text{eff}}(\mathbf{r}) d\mathbf{r}\) 両辺からポテンシャル積分項（$\int \rho v_{\text{ext}} d\mathbf{r}$ および $\int \rho v_{\text{eff}} d\mathbf{r}$）を消去すると、以下の極めて美しい関係式が得られる。 \(E_H[\rho] + E_{xc}[\rho] = \int_0^1 \langle \Psi_\lambda | \hat{V}_{ee} | \Psi_\lambda \rangle d\lambda\) これを整理することで、断熱結合公式（ACF）が導出される。 \(E_{xc}[\rho] = \int_0^1 \langle \Psi_\lambda | \hat{V}_{ee} | \Psi_\lambda \rangle d\lambda - E_H[\rho]\)

この公式の物理的な意味は極めて深い。相互作用エネルギーの期待値 $\langle \Psi_\lambda \hat{V}_{ee} \Psi_\lambda \rangle$ から古典的なHartreeエネルギー $E_H[\rho]$ を引いた非古典的な相互作用エネルギーを、結合定数 $\lambda$ について $0$ から $1$ まで積分（平均化）したものが、真の交換相関エネルギー $E_{xc}[\rho]$ に一致することを示している。この中には、実系における運動エネルギーの相関効果（$T[\rho] - T_s[\rho]$）が、断熱的に相互作用を導入する過程のポテンシャルエネルギー積分として自動的に包含されている。

5.2 Hartree-Fock交換（厳密交換）の積分形式

断熱結合公式における $\lambda = 0$ の極限では、電子間相互作用は完全に遮断されており、ハミルトニアンは独立粒子ハミルトニアンとなる。このとき、基底状態波動関数 $\Psi_0$ は非相互作用のKohn-Sham軌道 $\phi_i(\mathbf{r})$ から構成される単一の Slater 行列式で厳密に表される。

したがって、$\lambda = 0$ における非古典的相互作用の期待値は、軌道交換積分として以下のように厳密に計算可能である。 \(\langle \Psi_0 | \hat{V}_{ee} | \Psi_0 \rangle - E_H[\rho] = E_x^{HF}\) ここで、$E_x^{HF}$ は Hartree-Fock交換（厳密交換）エネルギーであり、スピン軌道表示（または閉殻系での空間軌道表示）を用いて以下の積分形式で記述される。

スピン分極を明示的に考慮した空間軌道表示では、スピン $\sigma \in {\uparrow, \downarrow}$ に対して以下のように書き表される。 \(E_x^{HF} = -\frac{1}{2} \sum_{\sigma \in \{\uparrow, \downarrow\}} \sum_{i,j}^{occ,\sigma} \iint \frac{\phi_{i\sigma}^*(\mathbf{r}) \phi_{j\sigma}(\mathbf{r}) \phi_{j\sigma}^*(\mathbf{r}') \phi_{i\sigma}(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r} d\mathbf{r}'\) この積分は、異なる空間点 $\mathbf{r}$ と $\mathbf{r}’$ における占有軌道の積同士の静電相互作用を表しており、非局所的（non-local）な特徴を持つ。自分自身との交換（$i=j$ の項）を含んでいるため、古典的なHartree反発に含まれる自己相互作用エネルギー（self-interaction energy）を数学的に正確に相殺する。

5.3 B3LYP (Becke, 3-parameter, Lee-Yang-Parr) 汎関数

断熱結合公式 $E_{xc} = \int_0^1 E_{xc,\lambda} d\lambda$ において、被積分関数 $E_{xc,\lambda} = \langle \Psi_\lambda | \hat{V}{ee} | \Psi\lambda \rangle - E_H[\rho]$ は、$\lambda$ の増加に伴って滑らかに変化する。ハイブリッド汎関数は、この $\lambda$ に関する積分を離散的な代表値の線形結合で数値積分近似することに相当する。

Beckeは、$\lambda=0$ （厳密交換 $E_x^{HF}$）と $\lambda=1$ （半局所的なDFT交換相関）の寄与を適切にブレンドすることで、化学反応や分子構造の記述が劇的に向上することを示した。その完成形が、パラメータフィッティングに基づく B3LYP 汎関数である。その定式化は以下の通りである。 \(E_{xc}^{B3LYP} = E_x^{LDA} + a_0 (E_x^{HF} - E_x^{LDA}) + a_x \Delta E_x^{B88} + E_c^{VWN} + a_c (E_c^{LYP} - E_c^{VWN})\) ここで、

  • $E_x^{LDA}$ は局所密度近似（ディラック交換）である。
  • $E_x^{HF}$ は非局所的なHartree-Fock厳密交換である。
  • $\Delta E_x^{B88} = E_x^{B88} - E_x^{LDA}$ は、Becke 88による交換エネルギーの勾配修正項である。
  • $E_c^{VWN}$ は、Vosko-Wilk-Nusair パラメータ化によるLDA相関汎関数である。
  • $E_c^{LYP}$ は、Lee-Yang-ParrによるGGA相関汎関数である。

3つのパラメータは、主族元素の原子化エネルギー, イオン化ポテンシャル, プロトン親和力などの実験データベース（G1セット）に対する二乗誤差が最小になるようにフィッティングされ、以下の値が決定された。 \(a_0 = 0.20, \quad a_x = 0.72, \quad a_c = 0.81\) B3LYPは、20%という適度な厳密交換の導入により、半局所汎関数に共通する自己相互作用誤差を大きく緩和し、特に有機分子の構造、反応熱、および活性化障壁の予測において圧倒的な標準として長年君臨している。

5.4 PBE0 汎関数

実験データへのフィッティングを排除し、純粋な物理的制約からハイブリッド汎関数を構築する試みとして開発されたのが PBE0 汎関数である。PBE0では、PBE交換相関汎関数をベースとして、厳密交換の混合比を理論的根拠に基づいて固定する。 \(E_{xc}^{PBE0} = \frac{1}{4} E_x^{HF} + \frac{3}{4} E_x^{PBE} + E_c^{PBE}\)

ここで、厳密交換の割合が $25\%$ （$1/4$）に固定されている理由は、結合定数 $\lambda$ に対する摂動論的考察に由来する。 ACFにおける被積分関数 $E_{xc,\lambda}$ を $\lambda = 0$ の周りでテイラー展開すると、以下のようになる。 \(E_{xc,\lambda} \approx E_{xc,0} + \left( \frac{d E_{xc,\lambda}}{d \lambda} \right)_{\lambda=0} \lambda + \dots\) これを積分すると、最低次の補正として以下の形式が得られる。 \(E_{xc} \approx E_{xc,0} + \frac{1}{2} \left( \frac{d E_{xc,\lambda}}{d \lambda} \right)_{\lambda=0} = \frac{1}{2} E_{xc,0} + \frac{1}{2} E_{xc,1}\) Perdewらは、多電子系における実際の $\lambda$ 依存性をさらに詳細に解析し、摂動展開の動的遮蔽効果を考慮すると、物理的なブレンド比は $1/4$ が最適値であることを数学的に示した。

PBE0は経験的なフィッティングパラメータを持たないにもかかわらず、分子の幾何構造や振動数だけでなく、固体の電子構造（特に半導体のバンドギャップの過小評価問題）を劇的に改善し、固体物理から量子化学まで広く信頼されている。

5.5 Heyd-Scuseria-Ernzerhof (HSE) 範囲分離型ハイブリッド汎関数

PBE0などの全領域（Full-range）ハイブリッド汎関数を金属や半導体などの固体（周期境界条件系）に適用する場合、非常に深刻な問題が発生する。実空間でのHartree-Fock交換積分のクーロン演算子 $1/r$ は長距離で極めて遅く減衰するため、計算収束に必要な実空間カットオフや $k$ 点サンプリング数が膨大になり、計算コストがGGAの100倍以上に跳ね上がる。また、金属中では長距離の静電反発は自由電子によって遮蔽（スクリーニング）されるため、長距離の裸の交換相互作用をそのまま取り入れることは物理的にも正しくない。

この問題を解決するために提案されたのが、HSE（Heyd-Scuseria-Ernzerhof）範囲分離型ハイブリッド汎関数である。HSEでは、2電子間のクーロン相互作用演算子 $1/r$ を、誤差関数 $\text{erf}$ および相補誤差関数 $\text{erfc}$ を用いて、短距離部（Short-Range: SR）と長距離部（Long-Range: LR）に数学的に分割する。 \(\frac{1}{r} = \underbrace{\frac{\text{erfc}(\omega r)}{r}}_{\text{SR}} + \underbrace{\frac{\text{erf}(\omega r)}{r}}_{\text{LR}}\) ここで、$\omega$ は範囲分離パラメータであり、実空間でのスクリーニング距離の逆数に対応する。

HSE汎関数の交換相関エネルギーは、この分割に基づいて以下のように定義される。 \(E_{xc}^{\text{HSE}} = a E_x^{HF,\text{SR}}(\omega) + (1-a) E_x^{PBE,\text{SR}}(\omega) + E_x^{PBE,\text{LR}}(\omega) + E_c^{PBE}\) ここで、

  • $a = 0.25$ は短距離部における厳密交換のブレンド比（25%）である。
  • $E_x^{HF,\text{SR}}(\omega)$ は、短距離ポテンシャル $\text{erfc}(\omega r)/r$ を用いて計算される短距離Hartree-Fock厳密交換エネルギーである。
  • $E_x^{PBE,\text{SR}}(\omega)$ および $E_x^{PBE,\text{LR}}(\omega)$ は、それぞれ短距離部および長距離部に対応するPBE交換エネルギーの寄与である。
  • 相関エネルギー $E_c^{PBE}$ は範囲分離を行わず、PBE相関をそのまま100%使用する。

スクリーニングパラメータ $\omega$ は、代表的な固体や分子の物理量（バンドギャップ、格子定数、形成熱）を最も良く再現するように最適化されており、代表的な HSE06 汎関数では以下のように設定されている。 \(\omega = 0.20 \text{ \AA}^{-1} \quad (\approx 0.11 \text{ a.u.}^{-1})\)

HSE06の最大の利点は、短距離部でのみ厳密交換を評価するため、実空間での交換相互作用が $\sim 1/\omega \approx 5 \text{ \AA}$ を超えると急速にゼロに減衰し、固体計算のコストをGGAの数倍程度に抑えることができる点である。同時に、固体の遮蔽効果を物理的に正しく記述できるため、半導体のバンドギャップ、格子定数、および欠陥エネルギーレベルを極めて高い精度で予測することが可能である。

6. 非局所相関・ダブルハイブリッド汎関数と分散力 (vdW) 相互作用

6.1 分散力（ファンデルワールス相互作用）の物理的起源と半局所汎関数の限界

分子間や層状物質（グラフェン、遷移金属ジカルコゲナイドなど）の層間に働く分散力（van der Waals: vdW 相互作用）は、電子密度の瞬間的な量子ゆらぎによって生じる一時的な誘起双極子間の長距離相関相互作用である。2つの離れた電気的に中性な系 $A$ と $B$ の間の距離を $R$ とすると、長距離極限での分散エネルギーは以下の漸近挙動を示す。 \(E_{\text{disp}}(R) \approx -\frac{C_6}{R^6} - \frac{C_8}{R^8} - \dots\)

これに対し、LDA、GGA、meta-GGAといった従来のすべての半局所（semi-local）汎関数は、空間の各点における局所的な電子密度とその勾配、および局所的な運動エネルギー密度のみを用いて交換相関エネルギーを表現する。 \(E_{xc}^{\text{semilocal}}[\rho] = \int f\big(\rho(\mathbf{r}), \nabla \rho(\mathbf{r}), \tau(\mathbf{r})\big) d\mathbf{r}\)

系 $A$ と系 $B$ が物理的に重なり合わないほど離れた距離 $R$ にあるとき、全電子密度は $\rho(\mathbf{r}) \approx \rho_A(\mathbf{r}) + \rho_B(\mathbf{r})$ となり、両者の密度の重なり（オーバーラップ）は距離に対して指数関数的に減衰する。 \(\rho_A(\mathbf{r}) \rho_B(\mathbf{r}) \propto e^{-\beta R}\) このため、半局所汎関数における交換相関エネルギーの相互作用項もまた、距離に対して指数関数的に減衰してしまう。 \(E_{xc}^{\text{semilocal}}(R) \propto e^{-\beta R}\)

これは、長距離でべき乗則（$-C_6/R^6$）に従って減衰する本質的な非局所分散力を、半局所汎関数が原理的に記述できないことを意味する。結果として、通常のGGAやmeta-GGAは、黒鉛の層間結合やDNA塩基対のスタッキング構造、自己組織化単分子膜などの非共有結合システムにおいて、引力を著しく過小評価（または完全に消失）し、誤った物理的予測を与える。

6.2 vdW-DF：非局所相関項と積分核 $\Phi$ の物理的意味

分散力をDFTの枠組みの中で第一原理に基づきシームレスに記述するため、2004年に M. Dion、H. Rydberg、E. Schröder、D. C. Langreth、B. I. Lundqvist らによって vdW-DF (van der Waals Density Functional) が提案された。

vdW-DFでは、全相関エネルギー $E_c[\rho]$ を以下のように局所的（あるいは半局所的）な相関項 $E_c^0[\rho]$ と、長距離の相関効果を担う非局所相関エネルギー項 $E_c^{nl}[\rho]$ に分割する。 \(E_c[\rho] = E_c^0[\rho] + E_c^{nl}[\rho]\) ここで、非局所相関エネルギー $E_c^{nl}[\rho]$ は、以下の2点空間積分（6次元積分）の形式で厳密に定式化される。 \(E_c^{nl}[\rho] = \frac{1}{2} \iint \rho(\mathbf{r}) \Phi(\mathbf{r}, \mathbf{r}') \rho(\mathbf{r}') d\mathbf{r} d\mathbf{r}'\)

ここで、$\Phi(\mathbf{r}, \mathbf{r}’)$ は非局所相関積分核（non-local correlation kernel）である。

  • 積分核 $\Phi$ の物理的意味： 積分核 $\Phi(\mathbf{r}, \mathbf{r}’)$ は、空間点 $\mathbf{r}$ と $\mathbf{r}’$ における電子密度のゆらぎが、分極ポテンシャルを介して互いに干渉し伝播するダイナミクスを表す。これは、誘電応答理論における乱雑位相近似（RPA）に基づき、電子ガスモデルの局所的応答性（プラズモン励起スペクトル）を媒介として導出される。 $\Phi$ は空間点間の距離 $d = |\mathbf{r} - \mathbf{r}’|$ に加え、各点での局所密度およびその勾配 $\big(\rho(\mathbf{r}), \nabla\rho(\mathbf{r}), \rho(\mathbf{r}’), \nabla\rho(\mathbf{r}’)\big)$ の陽な関数として定義される。

    距離 $d = |\mathbf{r} - \mathbf{r}’| \to \infty$ の大距離リミットにおいて、積分核 $\Phi$ は正確に以下の漸近挙動を示すように設計されている。 \(\lim_{|\mathbf{r} - \mathbf{r}'| \to \infty} \Phi(\mathbf{r}, \mathbf{r}') \propto -\frac{1}{|\mathbf{r} - \mathbf{r}'|^6}\) この漸近形を非局所相関エネルギーの式に代入すると、電子密度 $\rho(\mathbf{r})$ の積との積分により、経験的なパラメータフィッティングを行うことなく、厳密な分散力相互作用の長距離挙動 $-C_6/R^6$ が自動的かつ自然に回復する。 \(E_c^{nl} \xrightarrow{R \to \infty} -\frac{C_6}{R^6}\)

    vdW-DFは、原子ペアごとの分散力補正（DFT-D法など）のような経験的補正とは異なり、局所的な電子密度の変化に応じて分極率や分散係数 $C_6$ が自動的に変調されるため、より高精度で第一原理的な分散力の記述を可能にする。

6.3 ダブルハイブリッド（Double Hybrid）汎関数

ジェイコブの梯子の最上段である第5段に到達するアプローチの一つが、ダブルハイブリッド（Double Hybrid: DH）汎関数である。ダブルハイブリッド汎関数は、占有軌道の情報を利用するハイブリッド汎関数（第4段）をさらに拡張し、非占有軌道（仮想軌道: Virtual Orbitals）の情報をも記述に取り入れる。

具体的には、交換エネルギーに Hartree-Fock 厳密交換をブレンドするだけでなく、相関エネルギーの一部に2次多体摂動論（Møller-Plesset perturbation theory: MP2）による電子相関エネルギーをブレンドする。一般的なダブルハイブリッド汎関数の交換相関エネルギーは以下の形式で記述される。 \(E_{xc}^{\text{DH}} = (1 - a_x) E_x^{\text{DFT}} + a_x E_x^{HF} + (1 - a_c) E_c^{\text{DFT}} + a_c E_c^{\text{MP2}}\) ここで、

  • $E_x^{\text{DFT}}$ および $E_c^{\text{DFT}}$ は、ベースとなる半局所（GGAやmeta-GGA）交換および相関汎関数である。
  • $E_x^{HF}$ は、Kohn-Sham占有軌道を用いて評価されるHartree-Fock交換エネルギーである。
  • $E_c^{\text{MP2}}$ は、Kohn-Sham軌道（占有軌道 $\phi_i, \phi_j$ および非占有軌道 $\phi_a, \phi_b$）と、それぞれの軌道固有値 $\varepsilon$ を用いて計算されるMP2型相関エネルギーである。 \(E_c^{\text{MP2}} = -\frac{1}{4} \sum_{i,j}^{\text{occ}} \sum_{a,b}^{\text{vir}} \frac{\left| \langle \phi_i \phi_j || \phi_a \phi_b \rangle \right|^2}{\varepsilon_a + \varepsilon_b - \varepsilon_i - \varepsilon_j}\) ここで、$\langle \phi_i \phi_j || \phi_a \phi_b \rangle$ は、逆対称化された2電子クーロン積分（物理学者表記）である。 \(\langle \phi_i \phi_j || \phi_a \phi_b \rangle = \iint \frac{\phi_i^*(\mathbf{r}_1)\phi_j^*(\mathbf{r}_2)\big(\phi_a(\mathbf{r}_1)\phi_b(\mathbf{r}_2) - \phi_b(\mathbf{r}_1)\phi_a(\mathbf{r}_2)\big)}{|\mathbf{r}_1 - \mathbf{r}_2|} d\mathbf{r}_1 d\mathbf{r}_2\)

代表的なダブルハイブリッド汎関数である B2PLYP （Grimmeにより提案）では、パラメータは以下のように設定されている。 \(a_x = 0.53, \quad a_c = 0.27\) （ベースとして $E_x^{\text{DFT}} = E_x^{B88}$、$E_c^{\text{DFT}} = E_c^{LYP}$ を採用）。

  • ダブルハイブリッド汎関数の利点と課題： MP2相関エネルギー $E_c^{\text{MP2}}$ の計算において、非占有軌道への遷移励起エネルギー（分母の $\varepsilon_a + \varepsilon_b - \varepsilon_i - \varepsilon_j$）が組み込まれているため、長距離分散力を極めて高精度に自己矛盾なく記述でき、さらに自己相互作用誤差が大幅に削減される。その結果、分子の反応熱、結合解離エネルギー、水素結合、非共有結合相互作用において、量子化学計算における最高峰の精度を達成する。 しかし、非占有軌道に関する膨大な和を計算する必要があるため、計算コストは系のサイズ $N$ に対して $\mathcal{O}(N^5)$ でスケールし、標準的なDFTよりも大幅に重い。また、相関エネルギーの基底関数依存性が極めて強く、非常に大きな基底関数（Dunningの cc-pVTZ/cc-pVQZ 等）や基底関数完全系（CBS）への外挿を必要とする点に実用上の課題がある。

7. 代表的な交換相関汎関数の性能比較と適用限界

理論的な洗練度やジェイコブの梯子の階層の高さが、必ずしもすべての物理量に対して常に最高精度を保証するわけではない。実計算においては、計算コストと目標とする予測精度のバランスを考慮し、系に最適な汎関数を選択することが不可欠である。

7.1 バンドギャップの予測精度ベンチマーク

半導体や絶縁体の電子構造計算において、標準的な半局所DFTによるバンドギャップの極端な過小評価（underestimation）は、長年の深刻な課題であった。これは、局所的ポテンシャルが電子密度の微分不連続性を表現できないことに起因する。

以下に、窒化物半導体および多様な半導体・絶縁体（36種類の平均）におけるバンドギャップ予測の平均絶対誤差（MAE）および計算コスト比の比較を示す。

汎関数階層 代表的汎関数 平均絶対誤差 (MAE) 平均絶対誤差率 (MAPE) 計算コスト比 特徴と適用限界
第1段 (LDA) SVWN 1.22 eV 51.0% 1 構造記述は比較的安定しているが、バンドギャップは系統的に半減し、一部の半導体を金属と誤認する。
第2段 (GGA) PBE 1.17 eV 49.8% 1.5 ~ 2 LDAよりわずかに改善するが、依然として過小評価が極めて深刻。固体物性の標準的スクリーニングに利用される。
第3段 (meta-GGA) SCAN / r2SCAN 0.85 eV 35.2% 3 ~ 5 軌道運動エネルギー密度の導入により、GGAに比べバンドギャップが大幅に改善されるが、まだ十分ではない。
第4段 (Hybrid) HSE06 0.28 eV 11.1% 100 ~ 300 短距離部での厳密交換導入により、実験値に極めて近いバンドギャップを再現。半導体物性のデファクトスタンダード。
参考 (Semi-empirical) TB-mBJ 0.30 eV 12.1% 3 ~ 5 ポテンシャルを局所的に変調する経験的アプローチ。自己矛盾計算が不要で、GGA並みの超低コストでHSE06に迫るバンドギャップ精度を達成するが、全エネルギーの計算には使えない。

7.2 幾何構造と凝集エネルギーの予測精度

固体の格子定数や体積の記述において、LDAは軌道の重なりを過大評価するため 1% ~ 2% 程度の過小評価（過大結合：Overbinding）を引き起こす。一方、GGA（PBE）では逆に 1% ~ 2% 程度の過大評価（過小結合：Underbinding）を示す傾向がある。

第3段の meta-GGA である SCAN 汎関数は、多様な結晶構造において格子定数の平均絶対誤差を 0.5% 以内に抑え、極めて優れた構造予測能を示す。また、遷移金属酸化物の相安定性や、氷の水素結合ネットワークにおける多相安定相（氷の相図）の予測において、SCANは実験値を定性的・定量的に完全に再現することに成功しており、ハイブリッド汎関数に匹敵またはルーズする構造記述力を有している。

7.3 強相関系における局在化と対称性の破れ

遷移金属酸化物（NiO, CoOなど）や希土類元素を含む系などの強相関電子系（strongly correlated electron systems）では、 $d$ 軌道や $f$ 軌道に局在化した電子間の強いクーロン反発（オンサイト斥力 $U$）が支配的となる。 標準的なLDAやGGAは、自己相互作用誤差により電子密度を空間的に過度に非局在化（平滑化）させようとするため、本来モット絶縁体（Mott insulator）であるべき系を誤って金属と予測してしまう。

この課題に対し、実用的な回避策として、局在軌道にオンサイト斥力パラメータ $U$ を付加する DFT+U法 や、スピン対称性をあえて破る「対称性の破れた解（broken-symmetry solution）」を許容することが広く行われている。 SCAN汎関数は、強相関系における「軌道局在化の微弱な前駆現象」を、運動エネルギー密度 $\tau$ を介して敏感に検知し、DFT+Uのような手動パラメータなしに、銅酸化物超伝導体の母物質などの複雑な反強磁性磁気構造と絶縁体挙動を自律的に再現できることが確認されており、強相関系記述の限界を大きく押し広げている。

8. 結論と将来の展望

密度汎関数理論（DFT）における交換相関汎関数の発展の歴史は、電子多体問題という量子力学の極めて困難な壁に対する、理論家たちの論理的アプローチと絶え間ない挑戦の軌跡である。John Perdewが提唱した「ジェイコブの梯子」に沿って、最下段の局所密度近似（LDA）から、非局所相関や仮想軌道情報を内包する最上段（ダブルハイブリッドやRPA法）に至るまで、各階層はそれぞれ明確な物理的仮定と適用限界を伴いながら、物質科学の発展を支えてきた。

現代の交換相関汎関数の設計思想には、大きく分けて二つの対照的なアプローチが存在する。

  1. 非経験的アプローチ (Ab Initio / Constraint-Satisfaction)： PBEやSCANのように、実験値へのフィッティングを一切行わず、厳密な量子力学的制約条件（Lieb-Oxford限界、座標スケーリング等）を極限まで満足させる設計。予測の普遍性と信頼性が極めて高く、未知の極限環境や新物質開発において絶対的な威力を発揮する。
  2. 経験的アプローチ (Empirical / Database-Fitting)： B3LYPやMinnesota汎関数（M06など）、ダブルハイブリッド汎関数の多くに見られるように、特定の化学データベース（熱力学データ、反応障壁など）への適合度を最大化する設計。主族化学や一般的な有機化学反応において、極めて高い「化学的精度（sub-kcal/mol精度）」を提供する。

現場の計算科学においては、これら二つの潮流が互いの限界を補完し合う形で共存している。

今後の最大の展望として、機械学習（Machine Learning: ML）を用いた交換相関汎関数の構築が爆発的な勢いで進展している。これは、電子密度 $\rho(\mathbf{r})$ や各種局所記述子から交換相関エネルギー密度へのマッピングを、ニューラルネットワークなどの高表現力モデルに直接学習させる試みである。これにより、厳密な量子モンテカルロ（QMC）データやFCI（Full Configuration Interaction）の超高精度な密度情報を再現する「究極の汎関数」の実現が視野に入りつつある。

しかし、いかに高度な機械学習汎関数が登場したとしても、ディラック交換の $\rho^{1/3}$ 依存性や、断熱結合公式による厳密交換の定義、そして長距離における分散力の $-C_6/R^6$ 漸近挙動といった「物理の本質（厳密制約）」がその基盤であり続けることに変わりはない。DFTは、基礎理論の美しさと実用上の有用性の完璧なバランスを保ちながら、今後も未知なる極限環境やナノ材料、生命現象に至る物質科学の広大な地平を照らす灯火であり続けるであろう。

You May Also Enjoy