高校数学で理解したいコーン・シャム密度汎関数法 - 香川大学農学部ケミカルバイオロジー研究室

本記事の作成にあたり、説明・数式の確認と修正にChatGPTを活用しました。内容の最終的な責任は著者にあります。

本記事では高校数学を足掛かりに計算の流れを追います。汎関数微分などの必要事項は補足し、LDAの交換・相関汎関数の厳密な導出は扱いません。それでも、できるだけ簡単な系について、最後まで計算してみたいと思います。

主な参考書は

常田貴夫「密度汎関数法の基礎」2012年、講談社サイエンティフィク.

です。

以下、ハートリー原子単位系を使用します。記事中の電子密度 \(\rho(\mathbf r)\) は電子の数密度であり、空間全体で積分すると電子数 \(N\) になります。

ハートリー・フォック法の勉強過程は別記事にまとめています。高校数学の微分積分、行列の復習と、変分、汎関数微分、ラグランジュの未定乗数法など高校数学を超える部分の説明も別記事に書いてあります。

理論的背景の要点は、記事の最後にまとめています。

ベリリウム（Be）原子の計算

ハートリー・フォック・ローターン・ホール法による計算

スレーター行列式から正準ハートリー・フォック方程式までの復習

原子番号4のベリリウム原子は核電荷が4、電子数が4です。規格直交している軌道Aと軌道Bを逆スピンの2電子がそれぞれ占有した閉殻電子配置の波動関数を単一のスレーター行列式で表すと、

\[\Psi_{\mathrm{HF}}(x_1,x_2,x_3,x_4) = \frac{1}{\sqrt{24}}\begin{vmatrix}\psi_{A\alpha}(x_1) & \psi_{A\beta}(x_1) &\psi_{B\alpha}(x_1) & \psi_{B\beta}(x_1) \\
\psi_{A\alpha}(x_2) & \psi_{A\beta}(x_2) &\psi_{B\alpha}(x_2) & \psi_{B\beta}(x_2) \\
\psi_{A\alpha}(x_3) & \psi_{A\beta}(x_3) &\psi_{B\alpha}(x_3) & \psi_{B\beta}(x_3) \\ \psi_{A\alpha}(x_4) & \psi_{A\beta}(x_4) &\psi_{B\alpha}(x_4) & \psi_{B\beta}(x_4)
\end{vmatrix}\] となります。ここで、\(x\)は空間座標とスピンを合わせたスピン座標\(x=(r,\sigma)\)、\(\int dx = \Sigma_\sigma\int d^3 r\)。

ハミルトニアン（エネルギー演算子）は、\[\hat{H} = \left[-\frac{1}{2}\nabla^2_1 – \frac{4}{r_1} \right] + \left[-\frac{1}{2}\nabla^2_2 – \frac{4}{r_2} \right] + \left[-\frac{1}{2}\nabla^2_3 – \frac{4}{r_3} \right] + \left[-\frac{1}{2}\nabla^2_4 – \frac{4}{r_4} \right] + \frac{1}{r_{12}} + \frac{1}{r_{13}} + \frac{1}{r_{14}}+ \frac{1}{r_{23}}+ \frac{1}{r_{24}}+ \frac{1}{r_{34}}\]です。高校数学で習うシグマ記号を使うと、\[ \hat{H}= -\frac{1}{2}\sum_{i=1}^{4}\nabla_i^2
-\sum_{i=1}^{4}\frac{4}{r_i}+\sum_{1\le i<j\le4}\frac{1}{r_{ij}} \]と書けます。ここで、\(r_i=|\mathbf r_i|\) は原子核と電子 \(i\) の距離、\(r_{ij}=|\mathbf r_i-\mathbf r_j|\) は電子間距離です。

スレーター行列式を展開すると24個の項が現れます。エネルギー期待値を取ると\[E = \iiiint \Psi^*_{\mathrm{HF}}(x_1,x_2,x_3,x_4) \hat{H} \Psi_{\mathrm{HF}}(x_1,x_2,x_3,x_4) dx_1 dx_2 dx_3 dx_4\] となって、\(24 \times 24 =576 \)の組み合わせ掛ける\(14\)個の演算子（運動エネルギー4項、核引力4項、電子間反発6項の計14項）で合計\(8,064\)個の積分を計算しないといけないような気がしますが、一電子だけに作用する演算子、スピンの組み合わせで消える組み合わせ、積分結果（エネルギー期待値）が一緒なのでまとめられる項などがあって実際に計算する数は少ないです（スレーター・コンドン則）。

ベリリウム原子の電子ハミルトニアンの期待値は、\[\begin{aligned}E = &2\left[\int \phi^*_{A}\left(-\frac{1}{2}\nabla^2\right)\phi_{A} d\tau \right] + 2\left[\int \phi^*_{B}\left(-\frac{1}{2}\nabla^2\right)\phi_{B} d\tau \right] \\&- 8\int \frac{\phi^*_{A}\phi_{A}}{r}d\tau – 8\int \frac{\phi^*_{B}\phi_{B}}{r}d\tau \\&+ \iint \frac{|\phi_{A}(\mathbf r_1)|^2 | \phi_{A}(\mathbf r_2)|^2 }{r_{12}} d\tau_1 d\tau_2 \\&+ 4\iint \frac{|\phi_{A}(\mathbf r_1)|^2 | \phi_{B}(\mathbf r_2)|^2 }{r_{12}} d\tau_1 d\tau_2 \\&+ \iint \frac{|\phi_{B}(\mathbf r_1)|^2 | \phi_{B}(\mathbf r_2)|^2 }{r_{12}} d\tau_1 d\tau_2 \\&- 2\iint \frac{\phi_{A}^*(\mathbf r_1)\phi_{B}(\mathbf r_1) \phi_{B}^*(\mathbf r_2)\phi_{A}(\mathbf r_2)}{r_{12}} d\tau_1 d\tau_2 . \end{aligned}\] ここで、\(d\tau=d^3r\)、\(r_{12}=|\mathbf r_1-\mathbf r_2|\) です。したがって、与えられた軌道でこの期待値を評価するために必要なのは、運動エネルギー積分2種類、核引力積分2種類、クーロン積分3種類、交換積分1種類の、計8種類の空間積分です。

次に、軌道を最適化するために規格直交条件を課してエネルギーを変分します。ラグランジュの未定乗数法を使って、\[L = E – 2\lambda_A\left(\int \phi^*_{A}\phi_{A} d\tau – 1 \right) – 2\lambda_B\left(\int \phi^*_{B}\phi_{B} d\tau – 1 \right) -2\lambda_{AB}\left(\int \phi^*_{A}\phi_{B} d\tau\right) -2\lambda_{AB}^*\left(\int \phi^*_{B}\phi_{A} d\tau\right)\]を定義します。この\(L\)を軌道\(\phi^*_{A}\)と\(\phi^*_{B}\)で変分する（\(\phi^* \to \phi^* + \delta\phi^*\)と変化させて、\(\delta\phi^*\)について一次の項を取り出す）と汎関数微分（一次変分における\(\delta\phi^*\)の係数）は、\[\begin{aligned}\frac{\delta L}{\delta \phi^*_{A}} = &2\left(-\frac{1}{2}\nabla^2- \frac{4}{r}\right)\phi_{A} + 2\left(\int\frac{\phi_{A}^*(\mathbf r_2)\phi_{A}(\mathbf r_2)}{r_{12}}d\tau_2 \right)\phi_{A} + 4\left(\int\frac{\phi_{B}^*(\mathbf r_2)\phi_{B}(\mathbf r_2)}{r_{12}}d\tau_2 \right)\phi_{A} – \\ &2\left(\int\frac{\phi_{B}^*(\mathbf r_2)\phi_{A}(\mathbf r_2)}{r_{12}}d\tau_2 \right)\phi_{B} – 2\lambda_A \phi_{A} – 2\lambda_{AB} \phi_{B} \end{aligned}\] です。この汎関数微分をゼロと置いたのがハートリー・フォック方程式です。

後ろから三番目の項に\(\phi_{B}\)が掛かっていて何だか変に感じますが、「任意の関数\(f\)に作用させると、空間全体の\(f\)の値を、軌道\(\phi_{B}\)と距離に応じて重み付けして足し合わせ（積分）、その結果に、その場所の\(\phi_{B}\)を掛ける」演算子を交換演算子\(\hat{K}_{B}\)として定義すると、\[\left(\int\frac{\phi_{B}^*(\mathbf r_2)\phi_{A}(\mathbf r_2)}{r_{12}}d\tau_2 \right)\phi_{B} = \hat{K}_{B}\phi_{A}\] と書き直せて変数として\(\phi_{A}\)が現れます。同様に、「\(\phi_{A}\)に作用させると、空間全体の密度\(|\phi_B(\mathbf r_2)|^2\)が作るクーロンポテンシャルを、位置\(\mathbf r_1\)で評価し、\(\phi_A(\mathbf r_1)\)に掛ける」演算子をクーロン演算子\(\hat{J}_{B}\)として定義すると、\[\left(\int\frac{\phi_{B}^*(\mathbf r_2)\phi_{B}(\mathbf r_2)}{r_{12}}d\tau_2 \right)\phi_{A} = \hat{J}_{B}\phi_{A} \]と書き直せます。さらに一電子に作用する部分を\[\left(-\frac{1}{2}\nabla^2- \frac{4}{r}\right)\phi_{A} = \hat{h}\phi_{A}\]と書き直すと、全体は、\[\left(2\hat{h} + 2\hat{J}_{A} + 4\hat{J}_{B} -2\hat{K}_{B} -2\lambda_A \right)\phi_{A} -2\lambda_{AB}\phi_{B} = 0.\]全体を2で割って、\(\lambda\)の項を右辺に移すと、\[ \left(\hat{h} + \hat{J}_{A} + 2\hat{J}_{B} -\hat{K}_{B}\right)\phi_{A} = \lambda_A \phi_{A} + \lambda_{AB}\phi_{B}. \] さらに、\((\hat{J}_{A} – \hat{K}_{A})\phi_{A} = 0\)を利用して、両軌道に共通のフォック演算子の形に書き直すと、\[ \left(\hat{h} + 2\hat{J}_{A} – \hat{K}_{A} + 2\hat{J}_{B} -\hat{K}_{B}\right)\phi_{A} = \lambda_A \phi_{A} + \lambda_{AB}\phi_{B}. \]

同様に、\(\phi^*_{B}\)で変分した汎関数微分を0と置くと、\[ \left(\hat{h} + 2\hat{J}_{A} – \hat{K}_{A} + 2\hat{J}_{B} -\hat{K}_{B}\right)\phi_{B} = \lambda_B \phi_{B} + \lambda_{AB}^*\phi_{A} \]です。

演算子の部分をフォック演算子の記号で\[\hat{F} = \hat{h} + 2\hat{J}_{A} – \hat{K}_{A} + 2\hat{J}_{B} -\hat{K}_{B}\]と置き換えると、この2つの連立方程式は行列の形式で\[\hat{F} \begin{pmatrix}\phi_A & \phi_B\end{pmatrix} = \begin{pmatrix}\phi_A & \phi_B\end{pmatrix} \begin{pmatrix}\lambda_A & \lambda^*_{AB} \\ \lambda_{AB} & \lambda_B\end{pmatrix}\]と書けます。これがハートリー・フォック方程式です。\(\phi_{A}\)と\(\phi_{B}\)は規格直交（各軌道の絶対値の二乗の積分が1かつ一方の複素共役を掛けて積分したら0）していますが、未定乗数行列\(\begin{pmatrix}\lambda_A & \lambda^*_{AB} \\ \lambda_{AB} & \lambda_B\end{pmatrix}\)の非対角要素（右上と左下）は一般に\(0\)ではありません。そこでこれらが\(0\)になるように、軌道\(\phi_A\) と \(\phi_B\)を混ぜ合わせて\(\phi_{1s}\) と \(\phi_{2s}\)を作ると、正準ハートリー・フォック方程式\[\begin{aligned}\hat{F}\phi_{1s} &= \phi_{1s}\varepsilon_{1s}, \\ \hat{F}\phi_{2s} &= \phi_{2s}\varepsilon_{2s} \end{aligned}\] が得られます（球対称な閉殻Beの解を対象として、2つの占有正準軌道を1s、2sとここでは呼んでいます）。ここで、\(\hat{F} = \hat{h} + 2\hat{J}_{1s} – \hat{K}_{1s} + 2\hat{J}_{2s} -\hat{K}_{2s}\)です。また、\(\varepsilon_\pm = \frac{\lambda_A + \lambda_B}{2} \pm \sqrt{\left(\frac{\lambda_A – \lambda_B}{2}\right)^2 + |\lambda_{AB}|^2} \)という関係があります（二次方程式の解の公式の変形）。

ローターン・ホール方程式（3-21G基底関数系）

ベリリウム原子の3-21G基底関数はs型関数3個とp型関数6個の計9個です。
\[\begin{aligned}
\chi_1 = \chi_{s,1} &= 1.13361044\,e^{-71.8876r^2} +1.54675504\,e^{-10.7289r^2} +0.90270029\,e^{-2.22205r^2},\\
\chi_2 = \chi_{s,2} &=-0.36440247\,e^{-1.29548r^2}+0.32575156\,e^{-0.268881r^2}, \\
\chi_3 = \chi_{s,3} &=0.10453356\,e^{-0.0773501r^2}, \\
\chi_4 = \chi_{p_x,1} &=x\left(0.40412129\,e^{-1.29548r^2}+0.24356707\,e^{-0.268881r^2}\right),\\
\chi_5 = \chi_{p_y,1}&=y\left(0.40412129\,e^{-1.29548r^2}+0.24356707\,e^{-0.268881r^2}\right),\\
\chi_6 = \chi_{p_z,1}&=z\left(0.40412129\,e^{-1.29548r^2}+0.24356707\,e^{-0.268881r^2}\right),\\
\chi_7 = \chi_{p_x,2}&=0.05814551\,x\,e^{-0.0773501r^2},\\
\chi_8 = \chi_{p_y,2}&=0.05814551\,y\,e^{-0.0773501r^2},\\
\chi_9 = \chi_{p_z,2}&=0.05814551\,z\,e^{-0.0773501r^2}. \end{aligned} \]
SPシェルのs型とp型の指数は共通にしてあって、これによってプログラムでの計算が高速になるそうです。ハートリー・フォック方程式の\(\phi\)に\(\sum_{k=1}^9 C_k \chi_k = C_1 \chi_{1} + C_2 \chi_{2} + C_3 \chi_{3} + C_4 \chi_{4} + C_5 \chi_{5} + C_6 \chi_{6} + C_7 \chi_{7} + C_8 \chi_{8} + C_9 \chi_{9} \) を代入します。\[\hat{F} \sum_{k=1}^9 C_k \chi_k = \left(\sum_{k=1}^9 C_k \chi_k \right)\varepsilon_1 .\] この式の両辺に左から\(\chi_{1}\)を掛けて全空間で積分します。\(\chi_1\)と\(\chi_k\)で\(\hat{F}\)を挟んだ積分結果の数値を\(F_{1k}\)、\(\chi_1\)と\(\chi_k\)を掛け合わせた積分結果を\(S_{1k}\)とすると、\[\sum_{k=1}^{9} C_k F_{1k} = \left(\sum_{k=1}^{9} C_k S_{1k} \right)\varepsilon_1 \]が得られます。両辺とも9つの項からなります。これを\(\chi_9\)まで繰り返すと、9元連立方程式が得られます。行列形式に書き直すと、\[\begin{pmatrix}
F_{11} & F_{12} & \dots & F_{19} \\
F_{21} & F_{22} & \dots & F_{29} \\
\vdots & \vdots & \ddots & \vdots \\
F_{91} & F_{92} & \dots & F_{99}
\end{pmatrix}\begin{pmatrix}
C_{11} & C_{12} & \dots & C_{19} \\
C_{21} & C_{22} & \dots & C_{29} \\
\vdots & \vdots & \ddots & \vdots \\
C_{91} & C_{92} & \dots & C_{99}
\end{pmatrix} = \begin{pmatrix}
S_{11} & S_{12} & \dots & S_{19} \\
S_{21} & S_{22} & \dots & S_{29} \\
\vdots & \vdots & \ddots & \vdots \\
S_{91} & S_{92} & \dots & S_{99}
\end{pmatrix} \begin{pmatrix}
C_{11} & C_{12} & \dots & C_{19} \\
C_{21} & C_{22} & \dots & C_{29} \\
\vdots & \vdots & \ddots & \vdots \\
C_{91} & C_{92} & \dots & C_{99}
\end{pmatrix} \begin{pmatrix}
\varepsilon_{1} &0 & \dots & 0 \\
0& \varepsilon_{2} & \dots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \dots & \varepsilon_{9}
\end{pmatrix}\]が得られます。

1) 2種類の占有軌道について規格直交条件を満たすように係数を決めて、\(F_{ij}\)をすべて計算して数値を埋める。
2) 一般化固有値問題 \(\mathbf{FC}=\mathbf{SC}\varepsilon\) を解いて、軌道エネルギーと係数行列を求める。
3) 更新された占有軌道（下から2つ）の係数から、\(F_{ij}\)をすべて計算し直す。
4) 2) と 3) をエネルギーが収束するまで繰り返す。

ことで、HF軌道が求められます。2個の基底関数を使うヘリウム原子の時は計算しないといけない2電子積分の数はわずか6個でしたが、ベリリウム原子ではs型の3個だけ使う場合21個、p型も含めて9個を使う場合添字の置換による重複だけを除くと、最大1,035種類の2電子積分があります。ただし、球対称性があって、2px, 2py, 2pzは同じ動径部分を持ち、向きだけが異なるので、回転で対応する積分値を共用できます。したがって、実際は少ない数の動径積分からフォック行列を構築できます。
次に、s型関数とp型関数で挟まれたフォック行列要素について考えてみます。今回のs型占有軌道から作るフォック演算子は空間反転\((x,y,z) \to (-x,-y,-z)\)に対して不変であり、関数の偶奇性を保存します。したがって、偶関数のs型基底と奇関数のp型基底の間の行列要素はゼロになります。つまり、\(\hat{F}\chi_{s}\)は偶関数のまま、\(\hat{F}\chi_{p}\)は奇関数のままであり、偶関数と奇関数を掛けて積分するフォック行列要素は、\(F_{ps} = \int \chi^*_{p}(\mathbf r)(\hat{F}\chi_{s})(\mathbf r) d^3 r = 0\)、 \(F_{sp} = \int \chi^*_{s}(\mathbf r)(\hat{F}\chi_{p})(\mathbf r) d^3 r = 0\)です。したがって、9×9のフォック行列は、\[\mathbf{F} = \begin{pmatrix}
F_{11} & F_{12} & F_{13} & 0 & \dots & 0 \\
F_{21} & F_{22} & F_{23} & 0 & \dots & 0 \\ F_{31} & F_{32} & F_{33} & 0 &\dots & 0 \\ 0 & 0 & 0 & F_{44} &\dots & F_{49} \\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & F_{94}& \dots &F_{99}
\end{pmatrix}\]となり、占有軌道を求めるだけなら、球対称な閉殻Beの占有軌道がs型なので3×3行列だけで済みます。重なり行列\(\mathbf S\)についても、同じ偶奇性によって\(S_{sp} = S_{ps} = 0\)です。以下ではs型関数の部分だけで計算を進めます。

**図1．Be原子の3-21G基底関数系に含まれる3個のs型基底関数。** 原子核からの距離 (r) に対する \(\chi_1(r)\)、\(\chi_2(r)\)、\(\chi_3(r)\) を示す。各基底関数は全空間で規格化されているが、互いに直交しているわけではない。横軸の \(a_0\) はボーア半径を表す。

試行関数として、\[\phi_{A}^{(0)} = \chi_{1}\,\,\,(C_1=1, C_2 = C_3 = 0)\]と\[\phi_{B}^{(0)} = -0.141583535\chi_{1} + 1.009973216 \chi_{2} \,\,\,(C_3 = 0)\]を使います。

\(\chi_1\)と\(\chi_2\)の重なりは、\[\begin{aligned}S_{12} &= \int \chi_1\chi_2 d\tau\\ &= \int_0^{2\pi}\!\!\!\int_0^{\pi}\!\!\!\int_0^{\infty} \left(1.13361044\,e^{-71.8876r^2} +1.54675504\,e^{-10.7289r^2} +0.90270029\,e^{-2.22205r^2}\right) \left(-0.36440247\,e^{-1.29548r^2}+0.32575156\,e^{-0.268881r^2}\right) r^2 \sin\theta\, dr\, d\theta\, d\phi. \end{aligned}\] ガウス関数の積分公式を使って計算すると、\(S_{12} = 0.140185435\)です。

\[\phi_A = \chi_1,\,\,\,\,\,\phi_B = a\chi_1 + b\chi_2\]の係数\(a\)と\(b\)を決めたいので、直交条件から、\[\int \phi_A \phi_B d\tau = \int \chi_1 a\chi_1 d\tau + \int \chi_1 b\chi_2 d\tau = a\int \chi_1^2 d\tau + b \int \chi_1\chi_2 d\tau = 0\]なので、\[a + S_{12} b = 0.\]規格化条件から、\[\int \phi_B^2 d\tau = \int (a\chi_1 + b\chi_2)^2 d\tau = a^2\int \chi_1^2 d\tau + 2ab\int\chi_1\chi_2 d\tau + b^2\int\chi_2^2 d\tau =1. \]したがって、\[a^2 + 2S_{12}ab + b^2= 1.\] \(a = -S_{12}b\)を代入して、\(b > 0\)の解を選べば、\[b = \frac{1}{\sqrt{1-(S_{12})^2}}= 1.0099732161,\,\,\,\,a = -\frac{S_{12}}{\sqrt{1-(S_{12})^2}} = -0.1415835346.\]

1回目のフォック行列\(\mathbf{F}_s^{(0)}\)を上の初期軌道から計算すると、\[\mathbf{F}_s^{(0)} =\begin{pmatrix}
-4.371097127&-0.848118824&-0.903680504\\
-0.848118824&-0.232386963&-0.391305865\\
-0.903680504&-0.391305865&-0.332779004
\end{pmatrix} ,\] 重なり行列は\[\mathbf{S}_s=
\begin{pmatrix}
1&0.140185435&0.177121893\\
0.140185435&1&0.798850479\\
0.177121893&0.798850479&1
\end{pmatrix}\]です。したがって、1サイクル目に解くべきローターン・ホール方程式は、\[\begin{pmatrix}
-4.371097127&-0.848118824&-0.903680504\\
-0.848118824&-0.232386963&-0.391305865\\
-0.903680504&-0.391305865&-0.332779004
\end{pmatrix}\begin{pmatrix}C_{11}&C_{12}&C_{13}\\
C_{21}&C_{22}&C_{23}\\
C_{31}&C_{32}&C_{33}
\end{pmatrix} =\begin{pmatrix}
1&0.140185435&0.177121893\\
0.140185435&1&0.798850479\\
0.177121893&0.798850479&1
\end{pmatrix} \begin{pmatrix}
C_{11}&C_{12}&C_{13}\\
C_{21}&C_{22}&C_{23}\\
C_{31}&C_{32}&C_{33}
\end{pmatrix}\begin{pmatrix}
\varepsilon_1&0&0\\
0&\varepsilon_2&0\\
0&0&\varepsilon_3
\end{pmatrix}\]です。

この式を左辺 = 0にまとめて、「係数\(C\)がすべて0」以外の解を持つための条件として、行列式 = 0の条件で\(\varepsilon\)の値を求めれば、固有値\(\varepsilon\)と係数\(C\)を求めることができます。ヘリウムの時は二次方程式の解の公式で解けました。今回は三次方程式なので難しそうです。そこで、行列の道具を使って解くのですが、内部的には、\(\phi_A\)と\(\phi_B\)を混ぜて\(\phi_{1s}\)と\(\phi_{2s}\)を作ったように、\(\chi_{1},\chi_{2},\chi_{3}\)を混ぜて、重なり行列が単位行列\(\begin{pmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0& 0& 1\end{pmatrix}\)になるように基底を変換してから、フォック行列を変換して、新しいフォック行列の固有ベクトルと固有値を求めて、最後に元の基底に戻した係数を得る、という手順で計算するのが普通のやり方だそうです。

import numpy as np

np.set_printoptions(precision=9, suppress=True)

F = np.array([
    [-4.371097127, -0.848118824, -0.903680504],
    [-0.848118824, -0.232386963, -0.391305865],
    [-0.903680504, -0.391305865, -0.332779004],
])

S = np.array([
    [1.0,         0.140185435, 0.177121893],
    [0.140185435, 1.0,         0.798850479],
    [0.177121893, 0.798850479, 1.0        ],
])

# 重なり行列を対角化する（固有値と固有ベクトルを求める）
s, U = np.linalg.eigh(S)

print("固有値", s)
print("固有ベクトル", U)

# 新しい基底を規格直交させる変換
X = U @ np.diag(1.0 / np.sqrt(s)) @ U.T

# 新しい基底でフォック行列を表す
F_prime = X.T @ F @ X

print("新しいフォック行列\n", F_prime)

# 通常の固有値問題を解く
epsilon, C_prime = np.linalg.eigh(F_prime)

# 元の基底での軌道係数に戻す
C = X @ C_prime

print("係数行列\n",C)
print("固有値\n",epsilon)

得られた係数行列\(\mathbf{C}^{(1)}\)は、 \[\begin{pmatrix}-0.992573578 & -0.217113238 & -0.007268534\\-0.076483319 & 0.386531599 & -1.615062078\\ 0.027650271 & 0.678389502 & 1.528445305\end{pmatrix}\]です。各固有ベクトルの符号は独立に選べます。固有値（軌道エネルギー）が低い方2つを新しい占有軌道\[\phi_{A}^{(1)} =0.992573578 \chi_{1} + 0.0764833207 \chi_{2} -0.0276502727 \chi_{3},\]\[\phi_{B}^{(1)} =-0.2171132376 \chi_{1} +0.3865315999 \chi_{2}+0.6783895011 \chi_{3}\]として、フォック行列を更新して、以後収束するまで繰り返したのが以下の表です。全エネルギー差 \(10^{-10}E_h\) 未満、直交基底表示での密度行列変化のRMS \(10^{-8}\) 未満、SCF残差 \(10^{-8}E_h\) 未満を収束条件としました。

ITERATION	\(C_{11}\)	\(C_{21}\)	\(C_{31}\)	\(C_{12}\)	\(C_{22}\)	\(C_{32}\)	\(E\, /\, E_h\)	\(\varepsilon_1\,/\, E_h\)	\(\varepsilon_2\,/\, E_h\)
0	1	0	0	-0.1415835	1.00997322	0	-14.161712
1	0.99257358	0.07648332	-0.0276503	-0.2171132	0.3865316	0.6783895	-14.473945	-4.3855413	-0.1905787
2	0.99276829	0.07678131	-0.0288572	-0.2161103	0.25925336	0.79569118	-14.486358	-4.6185539	-0.289476
3	0.99280038	0.07652332	-0.0287721	-0.2157989	0.23498271	0.81734664	-14.486804	-4.675113	-0.302876
4	0.99280794	0.07644606	-0.0287369	-0.2157295	0.23040386	0.82140719	-14.48682	-4.685887	-0.3051334
5	0.9928095	0.07642921	-0.0287288	-0.2157157	0.22954139	0.82217116	-14.48682	-4.6879141	-0.3055482
6	0.99280981	0.07642586	-0.0287271	-0.2157131	0.22937897	0.822315	-14.48682	-4.6882955	-0.305626
7	0.99280986	0.07642521	-0.0287268	-0.2157126	0.22934838	0.82234209	-14.48682	-4.6883673	-0.3056406
8	0.99280988	0.07642509	-0.0287268	-0.2157125	0.22934262	0.82234719	-14.48682	-4.6883808	-0.3056434
9	0.99280988	0.07642507	-0.0287268	-0.2157125	0.22934153	0.82234815	-14.48682	-4.6883834	-0.3056439
10	0.99280988	0.07642506	-0.0287267	-0.2157125	0.22934133	0.82234833	-14.48682	-4.6883839	-0.305644
11	0.99280988	0.07642506	-0.0287267	-0.2157125	0.22934129	0.82234837	-14.48682	-4.6883839	-0.305644
12	0.99280988	0.07642506	-0.0287267	-0.2157125	0.22934128	0.82234837	-14.48682	-4.688384	-0.305644

最終的に収束した全エネルギーと空間占有軌道は、\[\begin{aligned}E_{\mathrm{HF/3-21G}} &= -14.48682\,E_h \\ \phi_{1s} &= 0.99280988\chi_1 + 0.07642506\chi_2 – 0.0287267\chi_3 \\ \phi_{2s} &= -0.2157125 \chi_1 + 0.22934128\chi_2 + 0.82234837\chi_3 \end{aligned}\]でした。

ベリリウムのハートリー・フォック限界は\(E_{\mathrm{HF}}^{\mathrm{limit}} = -14.573023\,E_h\)、厳密な非相対論的エネルギーは\(E_{\mathrm{exact}} = −14.66735649269\,E_h\) とのことです。

コーン・シャム・ローターン法による計算

コーン・シャム密度汎関数法の理論は後ろの方で説明しています。ここでは、Be原子の4個の電子が2つの空間軌道を占有する閉殻系についてコーン・シャム参照系を導入して、コーン・シャム方程式を導きます。

電子密度は2つのコーン・シャム空間軌道\(\phi_A\)と\(\phi_B\)を使って、\(\rho = 2(\phi^*_A\phi_A + \phi^*_B\phi_B)\)と表わされます。参照系の運動エネルギーはこのコーン・シャム軌道から計算して、全エネルギー汎関数を\[E[\rho] = 2\left(\int \phi^*_A \left(-\frac{1}{2}\nabla^2\right)\phi_A d\tau + \int \phi^*_B \left(-\frac{1}{2}\nabla^2\right) \phi_B d\tau\right) – \int \frac{4}{r}\rho(\mathbf r) d\mathbf{r} + \frac{1}{2}\iint \frac{\rho(\mathbf r)\rho(\mathbf r’)}{|\mathbf{r}-\mathbf{r}’|} d\mathbf{r}\, d\mathbf{r}’ + E_{xc}[\rho] \]と表わします。ここからは、ハートリー・フォック方程式の導出と手順は同じです。軌道の規格直交条件を課すため、ラグランジュ関数\[\mathcal{L} = E[\rho] – 2\lambda_A\left(\int\phi^*_A \phi_A d\tau – 1\right) – 2\lambda_B\left(\int\phi^*_B \phi_B d\tau – 1\right) – 2\lambda_{AB}\left(\int\phi^*_A \phi_B d\tau \right) – 2\lambda_{AB}^*\left(\int\phi^*_B \phi_A d\tau \right)\]を作ります。軌道の複素共役\(\phi^*\)についてそれぞれ変分して、汎関数微分を0と置きます。\[\begin{aligned}\frac{\delta\mathcal{L}}{\delta\phi^*_A} &= 2\left(-\frac{1}{2}\nabla^2\phi_A(\mathbf{r}) \right) – 2\left(\frac{4}{r}\phi_A(\mathbf{r})\right) +2\int\frac{\rho(\mathbf r’)}{|\mathbf{r} – \mathbf{r}’|}d\mathbf{r}’ \phi_A(\mathbf{r}) + 2v_{xc}(\mathbf{r})\phi_A(\mathbf{r}) – 2\lambda_A\phi_A(\mathbf{r}) – 2\lambda_{AB}\phi_B(\mathbf{r}) = 0, \\ \frac{\delta\mathcal{L}}{\delta\phi^*_B} &= 2\left(-\frac{1}{2}\nabla^2\phi_B(\mathbf{r}) \right) – 2\left(\frac{4}{r}\phi_B(\mathbf{r})\right) +2\int\frac{\rho(\mathbf r’)}{|\mathbf{r} – \mathbf{r}’|}d\mathbf{r}’ \phi_B(\mathbf{r}) + 2v_{xc}(\mathbf{r})\phi_B(\mathbf{r}) – 2\lambda_B\phi_B(\mathbf{r}) – 2\lambda_{AB}^*\phi_A(\mathbf{r}) = 0.\end{aligned}\] \(\rho = 2(\phi^*_A\phi_A + \phi^*_B\phi_B)\)なので、電子密度\(\rho\)を\(\phi^*\)で変分すると係数2が付きます。さらに、有効ポテンシャルを、\[v_{\mathrm{eff}}(\mathbf r) =-\frac{4}{r} +\int\frac{\rho(\mathbf r’)}{|\mathbf r-\mathbf r’|}\,d\mathbf r’+v_{xc}(\mathbf r)\]と定義すると、2つの方程式は\[
\begin{aligned}
\left(-\frac{1}{2}\nabla^2+v_{\mathrm{eff}}\right)\phi_A
&=\lambda_A\phi_A+\lambda_{AB}\phi_B,\\
\left(-\frac{1}{2}\nabla^2+v_{\mathrm{eff}}\right)\phi_B
&=\lambda_{AB}^*\phi_A+\lambda_B\phi_B
\end{aligned}\]と簡単に表せます。

\(\lambda_{AB} = \lambda_{AB}^* = 0\)になるように\(\phi_A\)と\(\phi_B\)を混合して新しく軌道\(\phi_{1s}\)と\(\phi_{2s}\)を作ると、正準コーン・シャム方程式\[
\begin{aligned}
\left(-\frac{1}{2}\nabla^2+v_{\mathrm{eff}}\right)\phi_{1s}
&=\varepsilon_1\phi_{1s},\\
\left(-\frac{1}{2}\nabla^2+v_{\mathrm{eff}}\right)\phi_{2s}
&=\varepsilon_2\phi_{2s}\end{aligned}\]になりました。ここで数学的に保証されているのは、方程式を正準形式に変換できることであって、具体的な\(\phi\)や\(\varepsilon\)を求めているわけではありません。

コーン・シャム法でも、ハートリー・フォック・ローターン法と同様に、基底関数を用いて行列の固有値問題を解きます。コーン・シャム・ローターン方程式を作るためには、軌道\(\phi\)を基底関数を足し合わせたもの（線形結合）で置き換えます。\[\begin{aligned}\left(-\frac{1}{2}\nabla^2+v_{\mathrm{eff}}\right)(C_{11} \chi_1 + C_{21} \chi_2 + C_{31}\chi_3)
&=\varepsilon_1(C_{11} \chi_1 + C_{21} \chi_2 + C_{31}\chi_3)\\\left(-\frac{1}{2}\nabla^2+v_{\mathrm{eff}}\right)(C_{12} \chi_1 + C_{22} \chi_2 + C_{32}\chi_3)
&=\varepsilon_2(C_{12} \chi_1 + C_{22} \chi_2 + C_{32}\chi_3).\end{aligned}
\] 次に、左から\(\chi_1\)を掛けて積分する計算を\(\chi_2\)、\(\chi_3\)についても行うと、ハートリー・フォック・ローターン方程式と同じように行列形式に変換されます。ただし、フォック行列\(\mathbf{F}\)の代わりに、コーン・シャム行列\(\mathbf{H}^{\mathrm{KS}}\)を使っています。\[\mathbf H^{\mathrm{KS}}\mathbf C
=\mathbf S\mathbf C\boldsymbol\varepsilon\]

コーン・シャム法の一電子演算子を確認します。演算子は\[\hat{h}_{\mathrm{KS}} = -\frac{1}{2}\nabla^2 + v_{\mathrm{ext}}(\mathbf r) + v_{\mathrm{H}}[\rho](\mathbf r) + v_{\mathrm{x}}[\rho](\mathbf{r}) + v_{\mathrm{c}}[\rho](\mathbf{r})\]です。ここで、\(v_{\mathrm{ext}}\) は原子核による外部ポテンシャル、\(v_{\mathrm H}\) は電子密度が作る静電ポテンシャル（ハートリーポテンシャル）、\(v_{\mathrm x}\) と \(v_{\mathrm c}\) は、それぞれ交換ポテンシャルと相関ポテンシャルです。これらのポテンシャルは位置の関数であり、軌道に掛け算することで作用します。

Be原子について、交換エネルギー汎関数としてLDA（ディラック交換）、相関エネルギー汎関数としてVWN5を使うと、球対称な電子密度 \(\rho(r)\) に対するコーン・シャム演算子は次のように表されます。\[\begin{aligned}
\hat h_{\mathrm{KS}}
={}&-\frac12\nabla^2-\frac{4}{r}\\
&+4\pi\left[
\frac{1}{r}\int_0^r\rho(t)t^2\,dt
+\int_r^\infty\rho(t)t\,dt
\right]\\
&-\left(\frac{3}{\pi}\right)^{1/3}\rho(r)^{1/3}
+v_c^{\mathrm{VWN5}}(\rho(r)).
\end{aligned}.\]ここで、第2行はハートリーポテンシャル、第3行の第1項はLDAによる交換ポテンシャルです。VWN5による相関ポテンシャルは複雑な関数形になっているのでここでは省略します。

LDA交換エネルギー汎関数の式は、\[E_{x}[\rho] = -\frac{3}{4}\left(\frac{3}{\pi}\right)^{\frac{1}{3}} \int \rho(\mathbf r)^{\frac{4}{3}}d\mathbf{r}\]です。3分の4乗は「4乗の3乗根 \(\sqrt[3]{x^4}\)」を表わします。この交換汎関数の汎関数微分を導くために、電子密度を\(\rho(\mathbf{r}) \to \rho(\mathbf{r}) + \delta\rho(\mathbf{r})\)と少しだけ変化させます。\(x^{\frac{4}{3}}\) を微分すると、\[\frac{d}{dx}x^{\frac{4}{3}}=\frac{4}{3}x^{\frac{1}{3}}\]です。したがって、\(x\)を微小に変化させた時の一次変分（二次以上の変化は無視）は\((x+\delta x)^{4/3}\approx x^{4/3}+\frac43x^{1/3}\delta x\)と近似できます。電子密度に当てはめると、電子密度を微小に変化させた時の、LDA交換エネルギーの変化は、\[\begin{aligned}
\delta E_x&=E_x[\rho+\delta\rho]-E_x[\rho]\\
&\approx-\left(\frac{3}{\pi}\right)^{\frac{1}{3}}\int\rho(\mathbf r)^{\frac{1}{3}}\delta\rho(\mathbf r)d\mathbf r
\end{aligned}\]となります。\(\delta\rho(\mathbf r)\) の係数が汎関数微分なので、\[v_x(\mathbf r) = \frac{\delta E_x[\rho]}{\delta\rho(\mathbf r)}=-\left(\frac{3}{\pi}\right)^{\frac{1}{3}}\rho(\mathbf r)^{\frac{1}{3}}\]となります。

なぜこれが交換ポテンシャルなのでしょうか。コーン・シャム方程式を導く時には、エネルギー汎関数を軌道で変分すると、\[\delta\rho(\mathbf r) = 2\phi_A(\mathbf r)\delta\phi^*_A(\mathbf r),\]\[\delta E_x = \int v_x(\mathbf r)\delta\rho(\mathbf r)\,d^3 r = \int 2v_x(\mathbf r)\phi_A(\mathbf r)\delta\phi^*_A(\mathbf r)\,d^3 r\]\[\frac{\delta E_x}{\delta\phi^*_A(\mathbf r)} = 2v_x(\mathbf r)\phi_A(\mathbf r)\]
という交換ポテンシャルを軌道に掛けた形が現れます。これは一電子に作用するポテンシャルエネルギーと同じ形なので、交換ポテンシャルと呼ばれます。

次に、コーン・シャムの一電子演算子を実関数\(\chi_1\)で挟んで積分したコーン・シャム行列要素\(H_{11}^{\mathrm{KS,(0)}}\)を計算します。

コーン・シャム行列要素\(H_{11}^{\mathrm{KS,(0)}}\)の計算

1) 初期密度を作る

KS-DFTにおいて電子密度はコーン・シャム軌道から作られます。電子密度は\(\rho = 2|\phi_A|^2 + 2|\phi_B|^2\)なので、初期KS軌道を\[\begin{aligned}
\phi_A^{(0)}(r)&=\chi_1(r),\\
\phi_B^{(0)}(r)&=-0.141583535\,\chi_1(r)
+1.009973216\,\chi_2(r)
\end{aligned}\]とすると、\[\begin{aligned}\rho^{(0)}(r)
&=2\chi_1(r)^2
+2\left[-0.141583535\,\chi_1(r)
+1.009973216\,\chi_2(r)\right]^2\\ &= 2.040091795 \chi_1(r)^2 + 2.040091794 \chi_2(r)^2 – 0.5719823127063942 \chi_1(r) \chi_2(r) \end{aligned}\]です。ここで、\[\begin{aligned} \chi_1(r)^2 =& (1.13361044)^2\,e^{-2\times 71.8876r^2} +(1.54675504)^2\,e^{-2\times 10.7289r^2} +(0.90270029)^2\,e^{-2\times 2.22205r^2}\\ &+ 2 (1.13361044)(1.54675504)\,e^{-(71.8876+10.7289)r^2} + 2 (1.13361044)(0.90270029)\,e^{-(71.8876+2.22205)r^2} \\ &+ 2 (1.54675504)(0.90270029)\,e^{-(10.7289+2.22205)r^2}, \\ \chi_2(r)^2 =& (0.36440247)^2\,e^{-2\times1.29548r^2}+(0.32575156)^2\,e^{-2\times 0.268881r^2}-2(0.36440247)(0.32575156)\,e^{-(1.29548+0.268881)r^2}, \\ \chi_1(r)\chi_2(r) =& (1.13361044)(-0.36440247)\,e^{-(71.8876+1.29548)r^2} + (1.13361044)(0.32575156)\,e^{-(71.8876+0.268881)r^2} \\ &+ (1.54675504)(-0.36440247)\,e^{-(10.7289+1.29548)r^2} + (1.54675504)(0.32575156)\,e^{-(10.7289+0.268881)r^2} \\ &+ (0.90270029)(-0.36440247)\,e^{-(2.22205+1.29548)r^2} + (0.90270029)(0.32575156)\,e^{-(2.22205+0.268881)r^2} \end{aligned} \]です。
このように、初期電子密度は15個のガウス関数の和として表せます。

**図2．Be原子の初期電子密度と動径電子数分布。** 左図は初期電子密度 \(\rho^{(0)}(r)\) を対数目盛で、右図は動径電子数分布 \(4\pi r^2\rho^{(0)}(r)\) を示す。実線は合計、破線と一点鎖線はそれぞれ初期軌道 \(\phi_A^{(0)}\)、\(\phi_B^{(0)}\) に由来する寄与を表す。右図の分布を半径について積分すると、その範囲の球殻内に含まれる電子数が得られる。これはSCF計算を開始するための初期分布であり、収束後の電子密度ではない。\(a_0\) はボーア半径を表す。

2) 運動エネルギーと核引力エネルギー

この部分は電子密度に依存しません。\[
T_{11}
=\int\chi_1\left(-\frac12\nabla^2\right)\chi_1\,d^3r,
\]
\[
V_{11}^{\mathrm{ext}}
=\int\chi_1\left(-\frac4r\right)\chi_1\,d^3r.
\]

\(\chi_1\)は3つの原始ガウス関数からなる\(\chi_1(\mathbf r) = A e^{-\alpha r^2} + B e^{-\beta r^2} + C e^{-\gamma r^2}\)ので、どちらも展開すると\(3×3=9\)項ですが、対称性により独立な積分は6種類です。

運動エネルギーの行列要素の計算

\[\begin{aligned}\langle A e^{-\alpha r^2}|\hat{T}| B e^{-\beta r^2}\rangle &= AB \int_0^{2\pi}\!\!\!\int_0^{\pi}\!\!\!\int_0^{\infty} e^{-\alpha r^2}\left(-\frac{1}{2}\right) \left(\frac{d^2}{dr^2} + \frac{2}{r}\frac{d}{dr}\right) e^{-\beta r^2} \cdot{} r^2 \sin\theta dr d\theta d\phi \\ &= 4\pi AB \int_0^{\infty} e^{-\alpha r^2} \left(-\frac{1}{2}\right) \left((-2\beta+4\beta^2 r^2) +\frac{2}{r}(-2\beta r )\right) e^{-\beta r^2}r^2\,dr \\ & = 4\pi AB\left\{
\int_0^{\infty}(3\beta – 2\beta^2 r^2) e^{-(\alpha + \beta) r^2} r^2\,dr \right\} \\ &= 4\pi AB\left\{ 3\beta\int_0^{\infty} r^2 e^{-(\alpha + \beta) r^2} dr – 2\beta^2 \int_0^{\infty} r^4 e^{-(\alpha + \beta) r^2} dr \right\} \end{aligned} \]

積分 \(\int_0^{\infty} r^n e^{-a r^2} dr\) （\(n\)が0以上の整数で\(a > 0\)の時）は、部分積分を繰り返すことで求められます。公式は、\(n\)が偶数（\(n = 2m\)）の時、\[I_{2m} = \frac{(2m)!}{2^{2m+1}m!}\frac{\sqrt{\pi}}{a^{m+1/2}},\]\(n\)が奇数（\(n = 2m+1\)）の時、 \[I_{2m+1} = \frac{m!}{2a^{m+1}}\]です。

今回必要なのは、\[
I_2(\alpha+\beta)
=\frac{\sqrt{\pi}}{4(\alpha+\beta)^{3/2}},
\qquad
I_4(\alpha+\beta)
=\frac{3\sqrt{\pi}}{8(\alpha+\beta)^{5/2}}
\]なので、運動エネルギーの行列要素は、\[
\begin{aligned}
\langle Ae^{-\alpha r^2}|\hat T
|Be^{-\beta r^2}\rangle
&=4\pi AB
\left[
\frac{3\beta\sqrt{\pi}}{4(\alpha+\beta)^{3/2}}
-\frac{6\beta^2\sqrt{\pi}}{8(\alpha+\beta)^{5/2}}
\right]\\
&=\frac{3\pi^{3/2}AB
\left[\beta(\alpha+\beta)-\beta^2\right]}
{(\alpha+\beta)^{5/2}}\\
&=
\frac{3\pi^{3/2}AB\,\alpha\beta}
{(\alpha+\beta)^{5/2}}
\end{aligned}
\]となります。

核引力エネルギーの行列要素の計算

\[\begin{aligned}\langle A e^{-\alpha r^2}|\hat{V}_{\mathrm{ext}}| B e^{-\beta r^2}\rangle &= AB \int_0^{2\pi}\!\!\!\int_0^{\pi}\!\!\!\int_0^{\infty} e^{-\alpha r^2}\left(-\frac{4}{r}\right) e^{-\beta r^2} \cdot{} r^2 \sin\theta dr d\theta d\phi \\ & = -16\pi AB\int_0^{\infty} r e^{-(\alpha + \beta) r^2} dr \\ &= -16\pi AB\int_0^{\infty} -\frac{1}{2(\alpha + \beta)} \left(\frac{d}{dr}e^{-(\alpha + \beta) r^2} \right) dr \\ &= \frac{8\pi AB}{\alpha + \beta}\left[ e^{-(\alpha + \beta) r^2} \right]_0^{\infty} \\ &= -\frac{8\pi AB}{\alpha + \beta}\end{aligned} \]

途中で \(\frac{d}{dr}e^{-(\alpha +\beta) r^2 } = -2(\alpha + \beta) r e^{-(\alpha +\beta) r^2} \) を使っています。

全ての積分を計算して足し合わせると、\[
T_{11}=6.905049356,\qquad
V_{11}^{\mathrm{ext}}=-14.815214657
\]
なので、
\[
h_{11}=T_{11}+V_{11}^{\mathrm{ext}}
=-7.910165301
\]
です。

3) ハートリー項

全電子密度\(\rho^{(0)}\)が作るハートリーポテンシャル\(v_{\mathrm{H}}\)を、基底関数\(\chi_1\)で期待値として評価した行列要素です。

\[
V_{11}^{\mathrm H}
=\iint
\frac{\chi_1(\mathbf r)^2\rho^{(0)}(\mathbf r’)}
{|\mathbf r-\mathbf r’|}
\,d^3r\,d^3r’.
\]
今回は、\[V_{11}^{\mathrm H} = 5.877695401\]です。

3-21基底関数系のs型部分には縮約（contracted）基底関数が3つ（\(\chi_1, \chi_2, \chi_3\)）しかないので、4つの縮約関数が全て異なる場合はありません。ただし、原始ガウス関数まで考えると、最初に4つの異なるガウス関数が含まれる積分を計算しておく必要があります（例えば、\(\chi_1\)の1つ目の原始関数、\(\chi_1\)の2つ目の原始関数、\(\chi_2\)の1つ目の原始関数、\(\chi_3\)の組み合わせ）。

同じ原子上の4つの原始s型ガウス関数の電子間反発積分では、\[\begin{aligned} (ij|kl) &= \iint \frac{N_i e^{-\alpha_i r_1^2} N_j e^{-\alpha_j r_1^2} N_k e^{-\alpha_k r_2^2} N_l e^{-\alpha_l r_2^2}}{|\mathbf{r}_1 – \mathbf{r}_2|} d\mathbf{r}_1 d\mathbf{r}_2 \\ &=N_i N_j N_k N_l \iint \frac{e^{-(\alpha_i + \alpha_j) r_1^2} e^{-(\alpha_k + \alpha_l) r_2^2}}{|\mathbf{r}_1 – \mathbf{r}_2|} d\mathbf{r}_1 d\mathbf{r}_2 \\ |\mathbf{r}_1 – \mathbf{r}_2| = \sqrt{r_1^2 + r_2^2 – 2r_1r_2\cos\theta} &\text{と置いて角度部分を先に積分すると、}\\ &= 16\pi^2 N_i N_j N_k N_l \left[ \int_0^\infty
r_1e^{-(\alpha_i + \alpha_j) r_1^2}
\left(\int_0^{r_1}r_2^2e^{-(\alpha_k + \alpha_l) r_2^2}dr_2\right)dr_1 \right. \\ &\,\,\,\left. +\int_0^\infty
r_2e^{-(\alpha_k + \alpha_l) r_2^2}
\left(\int_0^{r_2}r_1^2e^{-(\alpha_i + \alpha_j) r_1^2}dr_1\right)dr_2 \right]. \end{aligned}\] 積分を進めると最終的には、\[(ij|kl)=N_iN_jN_kN_l
\frac{2\pi^{5/2}}
{(\alpha_i+\alpha_j)(\alpha_k+\alpha_l)
\sqrt{\alpha_i+\alpha_j+\alpha_k+\alpha_l}}\]となります。

4) LDA交換項

ハートリー・フォック法では交換項を計算するために、基底関数の数に依存して、四つの基底関数を含む二電子積分（一般の分子では四中心積分）を計算する必要があります。一方、今回使うLDAでは交換項を電子密度から直接評価できるので、HF型の交換積分が必要ありません。

各半径で密度を求め、交換ポテンシャル\[ v_x(r)=-\left(\frac3\pi\right)^{1/3}\rho^{(0)}(r)^{1/3} \]を評価します。今回の電子密度では、密度の分数乗が積分に入るのでガウス積分の公式が使えず、数値積分する必要があります。このため、計算化学プログラムでDFT計算をする時には、座標のグリッドを指定するオプションがあります。\[
V_{11}^{x}
=-\left(\frac{3}{\pi}\right)^{1/3}
\int\chi_1^*(\mathbf r)\,
\rho(\mathbf r)^{1/3}\,
\chi_1(\mathbf r)\,d^3r
\]を数値積分すると、\[{V_{11}^{x}=-1.338674832}\]です。

5) VWN5相関項

\[ E_c^{\rm VWN5}[\rho]=\int\rho(\mathbf r)\varepsilon_c^{\rm VWN5}(\rho(\mathbf r))\,d^3r,
\qquad v_c(\rho)=\frac{d[\rho\varepsilon_c(\rho)]}{d\rho} \]

\[V_{11}^{c}
=4\pi\int_0^\infty
r^2\chi_1(r)^2v_c(r)\,dr\]を数値積分すると、\[V_{11}^{c}=-0.085261629\]です。

6) 5つを足し合わせる

\[H_{11}^{\mathrm{KS}}
=T_{11}+V_{11}^{\mathrm{ext}}
+V_{11}^{\mathrm H}+V_{11}^{x}+V_{11}^{c}\]なので、\[H_{11}^{\mathrm{KS,(0)}} = -3.456406361\]となります。

自己無撞着場（SCF）計算

同様に行列要素をすべて計算すると、\[H_s^{\mathrm{KS},(0)}
=
\begin{pmatrix}
-3.456406361&-0.765056014&-0.770064201\\
-0.765056014&-0.085156251&-0.263514314\\
-0.770064201&-0.263514314&-0.246804663
\end{pmatrix}\]となります。HF法の時のフォック行列と似たような値です。1サイクル目に解くべきコーン・シャム・ローターン方程式は\[\begin{pmatrix}
-3.456406361&-0.765056014&-0.770064201\\
-0.765056014&-0.085156251&-0.263514314\\
-0.770064201&-0.263514314&-0.246804663
\end{pmatrix}\begin{pmatrix}C_{11}&C_{12}&C_{13}\\
C_{21}&C_{22}&C_{23}\\
C_{31}&C_{32}&C_{33}
\end{pmatrix} =\begin{pmatrix}
1&0.140185435&0.177121893\\
0.140185435&1&0.798850479\\
0.177121893&0.798850479&1
\end{pmatrix} \begin{pmatrix}
C_{11}&C_{12}&C_{13}\\
C_{21}&C_{22}&C_{23}\\
C_{31}&C_{32}&C_{33}
\end{pmatrix}\begin{pmatrix}
\varepsilon_1&0&0\\
0&\varepsilon_2&0\\
0&0&\varepsilon_3
\end{pmatrix}\]です。HF法の時と同じPythonコードで解いてみると、係数行列\[\mathbf{C}^{(1)}=\begin{pmatrix}-0.988139323 &-0.233032528 &0.0408634 \\ -0.108402074 & 0.154002943 &-1.6517306 \\ 0.037472168 & 0.887127481 & 1.417294102\end{pmatrix}\]が得られました。各固有ベクトルの符号は独立に選べます。以下では、固有値（軌道エネルギー）が低い方2つを新しいKS占有軌道\[\phi_{A}^{(1)} =0.988139323 \chi_{1} + 0.108402074 \chi_{2} -0.037472168 \chi_{3},\]\[\phi_{B}^{(1)} =-0.233032528 \chi_{1} +0.154002943 \chi_{2} +0.887127481 \chi_{3}\]として、新しい電子密度を作って、コーン・シャム行列を更新して、以後エネルギーが収束するまで繰り返します。この初期値と球対称・閉殻の制約の下では、以下の解に収束しました。

全エネルギーは、\[E_{\mathrm{KS}}=2\langle\phi_A|\hat{ h}|\phi_A\rangle+2\langle\phi_B|\hat{h}|\phi_B\rangle+E_{\mathrm H}[\rho]+E_{\mathrm x}[\rho]+E_{\mathrm c}[\rho].\]

Iteration	\(E\,/\,E_h\)	\(C_{11}\)	\(C_{21}\)	\(C_{31}\)	\(C_{12}\)	\(C_{22}\)	\(C_{32}\)	\(\varepsilon_{1}\,/\,E_h\)	\(\varepsilon_{2}\,/\,E_h\)
0	-14.037457	1	0	0	-0.1415835	1.00997322	0
1	-14.354198	0.98813932	0.10840208	-0.0374722	-0.2330325	0.15400294	0.88712748	-3.4809809	-0.0826515
2	-14.357123	0.98952234	0.09953175	-0.035148	-0.229482	0.24107613	0.81129465	-3.8380542	-0.2150877
3	-14.357197	0.98934164	0.10071279	-0.0354598	-0.2299793	0.22711445	0.82365246	-3.7940845	-0.2012269
4	-14.357199	0.98936947	0.10053049	-0.0354111	-0.2299061	0.22934949	0.82167885	-3.8010036	-0.2034579
5	-14.357199	0.98936502	0.10055964	-0.0354189	-0.2299178	0.22899218	0.82199449	-3.7998968	-0.2031011
6	-14.357199	0.98936573	0.10055498	-0.0354177	-0.2299159	0.2290493	0.82194403	-3.8000737	-0.2031582
7	-14.357199	0.98936562	0.10055572	-0.0354179	-0.2299162	0.22904017	0.8219521	-3.8000454	-0.2031491
8	-14.357199	0.98936564	0.1005556	-0.0354178	-0.2299162	0.22904163	0.82195081	-3.8000499	-0.2031505
9	-14.357199	0.98936564	0.10055562	-0.0354178	-0.2299162	0.2290414	0.82195101	-3.8000492	-0.2031503
10	-14.357199	0.98936564	0.10055562	-0.0354178	-0.2299162	0.22904143	0.82195098	-3.8000493	-0.2031503
11	-14.357199	0.98936564	0.10055562	-0.0354178	-0.2299162	0.22904143	0.82195099	-3.8000493	-0.2031503

Be原子の収束したLDA（Dirac交換＋VWN5相関）／3-21G のエネルギー内訳は、以下のとおりです。単位はすべてhartree（E_h）です。

成分	記号	エネルギー
KS非相互作用系の運動エネルギー	\(T_s\)	\(+14.4229115261\)
核引力エネルギー	\(V_{\mathrm{ne}}\)	\(-33.3982887803\)
ハートリーエネルギー	\(J\)	\(+7.1480571442\)
交換エネルギー	\(E_x\)	\(-2.3050452405\)
相関エネルギー	\(E_c\)	\(-0.2248339062\)
全エネルギー	\(E\)	\(-14.3571992567\)

厳密なKS理論では、相関エネルギーには相関運動エネルギー\(T_c\)が含まれています。定義は相互作用系とKS参照系の運動エネルギーの差\(T_c=T−T_s\)です。ここではVWN5相関汎関数の一様座標スケーリングから、近似的な\(T_c\)と\(U_c\)を評価しました（運動エネルギーは座標の2階微分、クーロン反発エネルギーは距離の逆数を含むので、波動関数を一様に座標スケーリング \(\psi_\gamma(\mathbf r_1,\dots,\mathbf r_N)=\gamma^{3N/2}\psi(\gamma\mathbf r_1,\dots,\gamma\mathbf r_N)\) とスケーリングすると、運動エネルギーの期待値は \(\gamma^2\)倍、電子間反発エネルギーの期待値は \(\gamma\) 倍になります。厳密な相関汎関数ではこのスケーリング則と変分原理から、相関エネルギーについて\(\left.\frac{dE_{\mathrm c}[\rho_\gamma]}{d\gamma}\right|_{\gamma=1}
=2T_{\mathrm c}+U_{\mathrm c}\)が成り立ちます）。以下の値はVWN5の枠内での分解であり、Be原子の厳密な運動エネルギーや電子間反発エネルギーを直接求めたものではありません。\[
E_c=T_c+U_c
=\underbrace{0.1381078096}_{T_c}
+\underbrace{-0.3629417158}_{U_c}
=-0.2248339062\]です。

比較

CCSD(T)は、結合クラスター法と呼ばれるpost-HF法の高精度な計算レベルです。DFT汎関数の検証や開発に用いる高精度な参照計算の一つです。ただし、厳密解ではなく、強い静的相関や基底関数の不足による誤差があり得ます。

CCSD(T)列は、HFの各成分をそのまま残し、全エネルギーの差をまとめて相関補正として示した帳簿上の分解です。上4行はCCSD(T)の密度・密度行列から評価した各エネルギー成分ではありません。

成分	HF/3-21G	LDA(VWN5)/3-21G	HF/cc-pCVQZ-F12	CCSD(T)/cc-pCVQZ-F12 <HF基準>
運動エネルギー\(T\, \text{or}\, T_s\)	+14.478061	+14.422912	+14.572967	+14.572967
核引力エネルギー\(V_{\mathrm{ne}}\)	−33.463645	−33.398289	−33.635078	−33.635078
ハートリーエネルギー\(J\)	+7.166001	+7.148057	+7.156002	+7.156002
交換エネルギー\(E_x^{\mathrm{HF}}\, \text{or}\, E_x^{\mathrm{LDA}}\)	−2.667237	−2.305045	−2.666903	−2.666903
相関エネルギー\(E_c^{\mathrm{HF}基準} \,\text{or}\, E_c^{\mathrm{VWN5}}\)	0	−0.224834	0	−0.090533
全エネルギー	−14.486820	−14.357199	−14.573013	−14.663545

i) HF/cc-pCVQZ-F12とCCSD(T)/cc-pCVQZ-F12レベルはORCA 6.1.1で計算した。ii) ハートリー・フォック限界は-14.573023 E_h、厳密な非相対論的エネルギーは−14.66735649269 E_h。iii) CCSD(T)-F12/cc-pCVQZ-F12レベルでの全エネルギーは−14.666660946795 E_h。iV) CCSD(T)/cc-pXVQZ-F12のHF基準の相関エネルギーは、CCSD(T)の緩和応答密度に基づく分解によって、運動エネルギー成分\(\Delta T = +0.084797991 E_{h}\) とポテンシャルエネルギー成分\(\Delta V = -0.175330519 E_h\)に分けられる。

今回、LDA交換の絶対値はHF交換より小さく、VWN5相関の絶対値はCCSD(T)とHFのエネルギー差より大きいという結果でした。ただし、密度と定義が異なるため、各成分の誤差を直接比較はできません。

**図3．Be原子の占有1s・2s軌道の比較。** HF/3-21G、KS-LDA（Dirac交換＋VWN5相関）/3-21G、およびHF/cc-pCVQZ-F12による軌道を示す。上段は動径関数 \(R(r)\)、下段は \(u(r)=rR(r)\) で、左列が1s軌道、右列が2s軌道である。本文のs型空間軌道との関係は \(\phi(\mathbf r)=R(r)/\sqrt{4\pi}\) であり、\(\int_0^\infty |R(r)|^2r^2dr=\int_0^\infty |u(r)|^2dr=1\) と規格化される。軌道全体の符号は任意である。

**図4．Be原子の電子密度と動径電子数分布の比較。** HF/3-21G、KS-LDA（Dirac交換＋VWN5相関）/3-21G、およびCCSD(T)/cc-pCVQZ-F12による結果を示す。CCSD(T)には緩和密度（relaxed density。外部摂動に対するエネルギーの応答から求め、軌道の応答も含めた密度）を使用した。左図は電子の数密度 \(\rho(r)\) を対数目盛で、右図は動径電子数分布 \(4\pi r^2\rho(r)\) を示す。\(4\pi r^2\rho(r)\,dr\) は半径 \(r\) から \(r+dr\) の球殻内に含まれる電子数を表し、全空間で積分するとBe原子の電子数4になる。

**図5．CCSD(T)を基準としたBe原子の電子密度の差。** CCSD(T)/cc-pCVQZ-F12の緩和密度を基準として、HF/3-21GおよびKS-LDA（Dirac交換＋VWN5相関）/3-21Gとの差を比較した。左図は \(\Delta\rho(r)=\rho_{\mathrm{各方法}}(r)-\rho_{\mathrm{CCSD(T)}}(r)\)、右図は動径電子数分布の差 \(4\pi r^2\Delta\rho(r)\) を示す。正の値は基準より電子密度が大きく、負の値は小さいことを表す。計算方法と基底関数系の両方が異なるため、この差には両者の影響が含まれる。

コーン・シャム密度汎関数法の理論的背景の勉強

系の波動関数\(\Psi\)の絶対値の二乗 \(|\Psi|^2\) から電子密度\(\rho\)を求めることができます。電子密度は、ほかの電子の座標を積分し、スピンについて和を取り、電子数 \(N\) を掛けて求めます。波動関数理論（wave function theory, WFT）では、外部ポテンシャル\(v_{\mathrm{ext}}\)（原子核の位置と原子番号）と電子の数\(N\)からハミルトニアンが決まるので、\[(N,v_{\mathrm{ext}}) \rightarrow \hat{H} \rightarrow \Psi_0 \rightarrow \rho_0\]です。

一方、DFTの理論的基礎となるのは、基底状態について逆向きの対応も成立するという事実です。以下では、基底状態が縮退していない（全体の位相の違いを除いて基底状態が一意に決まる）場合を考えます。基底状態の電子密度から、電子の数が分かり、電子密度に対応する原子核の位置と原子の種類（外部ポテンシャル）が（加法定数 \(+ C\) を除いて）一意に決まるので、ハミルトニアンが決まって、波動関数も決まって、（エネルギーの零点を固定すれば）エネルギーも決まることが証明されています。\[\rho_0 \rightarrow (N,v_{\mathrm{ext}}) \rightarrow \hat{H} \rightarrow \Psi_0 \rightarrow E_0 .\] さらに、エネルギーを電子密度の汎関数として表すことができ、変分原理から、エネルギー（汎関数）を最小化することで基底状態の電子密度を求めることができます。ここで、基底状態密度から外部ポテンシャルが一意に決まるという話と、ある系の基底状態密度を変分によって求める議論を区別しておきます。後者では、原子核の配置などが定める外部ポテンシャルと電子数を固定して試行密度を変化させます。

ポイントを箇条書きで整理します。

原子や分子の系の非相対論的シュレーディンガー（Schrödinger）方程式を解析的に解き、波動関数を求めることは一般的に困難です。
波動関数を単一のスレーター行列式に制限して、エネルギーを最小化することで得られるのがハートリー・フォック方程式です。シュレーディンガー方程式の厳密解を与えるわけではありません。
DFTでは、基底状態のエネルギーを電子密度の汎関数として表わせます。これは、電子密度という関数を入力するとエネルギーという数値が返るということです。ただし、表わせることが分かっているだけで、厳密な汎関数の具体的な形は一般に分かっていません。
特に、運動エネルギーや交換・相関エネルギーを、電子密度から直接計算することが難題です。運動エネルギーの誤差は、結合エネルギーなどの小さなエネルギー差の評価に大きく影響します。
そこで、軌道から運動エネルギーと電子密度を計算して、残りを電子密度の汎関数で扱う枠組みとして、コーン・シャム（Kohn–Sham）法が考案されました。
コーン・シャム方程式を基底関数で展開して、方程式を行列表示したものが、コーン・シャム・ローターン（Kohn–Sham–Roothaan）方程式です。
電子密度による変分原理は、レヴィ（Levy）の制限付き探索によって、より一般的な密度を扱える形に定式化されました。ホーヘンバーグ・コーン (Hohenberg–Kohn) の第2定理はこの制限付き探索に含まれます。電子密度と外部ポテンシャル（原子核の位置と原子の種類）の関係を説明するための重要な定理としてホーヘンバーグ・コーンの第1定理があります。

以下では原子核を固定して、核間反発を除いた電子系のエネルギーを考えます。

数学的な証明上、DFTにおける変分計算は、「外部ポテンシャル（原子核）を考えずに、与えられた電子密度に対して運動エネルギーと電子間反発エネルギーを最小化する反対称化された波動関数を探索する」過程と「電子密度を変化させて全エネルギーを最小化する」過程に分けられます。前者の探索で得られる最小値を、密度の汎関数として定義したものが「普遍汎関数（universal functional）」\[F[\rho]=\min_{\Psi\to\rho} \langle\Psi|\hat T+\hat V_{ee}|\Psi\rangle\]です。汎関数なので、電子密度関数を入力するとエネルギーを返しますが、内部的には「与えられた電子密度に対して運動エネルギーと電子間反発エネルギーを最小化する反対称化された波動関数を探索」しています。3種類の電子密度\(\rho_A,\rho_B,\rho_C\)を考え、この仮想例では、この3種類が許される密度のすべてであるとします。それぞれの電子密度を生成する波動関数は無数に存在しますが、各密度を固定して、\(\langle\Psi|\hat T+\hat V_{ee}|\Psi\rangle\)を最小にする\(\Psi\)を探索しておきます。そして、外部ポテンシャルも含めて全エネルギーを評価します。\[\begin{aligned}&E[\rho_A] = F[\rho_A] + \int \rho_A(\mathbf r)v_{\mathrm{ext}}(\mathbf r) d\mathbf{r} = -0.4\,\,\,（探索された波動関数は\Psi_{A}）\\ &E[\rho_B] = F[\rho_B] + \int \rho_B(\mathbf r)v_{\mathrm{ext}}(\mathbf r) d\mathbf{r} = -0.5\,\,\,（探索された波動関数は\Psi_{B}）\\&E[\rho_C] = F[\rho_C] + \int \rho_C(\mathbf r)v_{\mathrm{ext}}(\mathbf r) d\mathbf{r} = -0.2\,\,\,（探索された波動関数は\Psi_{C}）\end{aligned}\] この3つの中では\(\rho_B\)が全エネルギーの最小値を与える密度\(\rho_0\)なので、\(\Psi_{B}\)が真の基底状態波動関数\(\Psi_{0}\)です。これは数学的な理論上の話で、実際のDFT計算手順とは違います。

理論上、基底状態波動関数\(\Psi_{0}\)は既に探索済みなので、普遍汎関数の中の運動エネルギー期待値はこの波動関数から定義できます。\[T[\rho_0] = \int\cdots\int \Psi^*_{0} \hat{T} \Psi_{0} dx_1\cdots dx_N = \langle\Psi_0|\hat{T}|\Psi_0\rangle\] 次に、電子間反発エネルギー\(\ \int\cdots\int \Psi^*_{0} \hat{V}_{ee} \Psi_{0} dx_1\cdots dx_N = \langle\Psi_0|\hat{V}_{ee}|\Psi_0\rangle\)のうち、古典的な電子間反発エネルギー（ハートリー項）は波動関数ではなく電子密度だけから定義できます。\[J[\rho_0] = \frac{1}{2}\iint\frac{\rho_0(\mathbf r)\rho_0(\mathbf r’)}{|\mathbf{r} – \mathbf{r’}|}d^3 r d^3 r’ \] そして、普遍汎関数の残りの電子間相互作用の非古典的部分を\(W_{nc}\)と置きます。\[W_{nc}[\rho_0] = F[\rho_0] – T[\rho_0] – J[\rho_0] \] 上の期待値の定義から非古典的部分を直接評価するには基底状態波動関数\(\Psi_{0}\)が必要です。さらに、核引力エネルギーを記号で置き換えておきます。\[ E_{\mathrm{ext}}[\rho_0] = \int \rho_0(\mathbf r)v_{\mathrm{ext}}(\mathbf r) d^3 r \] したがって、全エネルギーを各要素に分割すると以下の式となります。\[E[\rho_0] = T[\rho_0] + E_{\mathrm{ext}}[\rho_0] + J[\rho_0] + W_{nc}[\rho_0] \]

運動エネルギー \(T[\rho_0]\) は理論上、波動関数\(\Psi_{0}\)から計算され、電子密度から直接計算するよい近似汎関数がありませんでした。そこで、コーンとシャムが導入したのが、実際の系と同じ基底状態電子密度を持つ仮想的な「相互作用のない参照系 (non-interacting reference system)」です。この「参照系」のそれぞれの電子は、有効ポテンシャル\(v_{s}\)を感じながら運動していると考えます。「相互作用のない」とは、この参照系のハミルトニアンに、実際の電子系のハミルトニアンに含まれる二体演算子（例えばヘリウム原子では電子間反発演算子は\(\frac{1}{r_{12}}\)）などが含まれない、という意味です。「厳密な交換・相関汎関数（上の\(W_{nc}[\rho]\)ではなく後で定義される\(E_{xc}[\rho]\)）を使って、適切な参照系が存在する」という条件の下では、基底状態の解は実際の系の電子密度を再現します。

文献によって「補助系（auxiliary system）」という言い方と「参照系（reference system）」という言い方（あるいは両方 “auxiliary, noninteracting, “Kohn-Sham” reference system“）が使われています。ニュアンスとしては、「補助系」が計算上の仮想的な系であることを強調した言い方で、「参照系」は実在の系とのエネルギー差を定義する基準の系であることを強調した言い方でしょうか。

参照系の一電子軌道をコーン・シャム軌道と呼びます。ここでは、これらの軌道から単一のスレーター行列式 \(\Phi_s\) を作れる場合を考えます。この行列式を使って、参照系の運動エネルギー \(T_s\) を計算します。相互作用のない参照系（コーン・シャム参照系）と実際の系では運動エネルギーが一般に異なります。その差を\[\Delta T = T – T_s\]と定義します。この差を他の部分で補正しなくてはいけません。

ここで、電子間相互作用の非古典的エネルギー\(W_{nc}[\rho_0]\)を交換項と相関の電子間反発成分に分けます。\[W_{nc}[\rho_0] = E_x[\rho_0] + U_c[\rho_0].\]

\(E_x[\rho_0] \) は、KS軌道から作ったスレーター行列式 \(\Phi_s\) を使って計算した電子間反発エネルギーの期待値とハートリーエネルギー\(J[\rho_0]\)の差として定義されます。\[E_x[\rho_0] = \langle\Phi_s|\hat V_{ee}|\Phi_s\rangle – J[\rho_0]. \] この交換項によって、ハートリー項に含まれる各占有スピン軌道の自己相互作用が相殺されます（この相殺は上記の厳密交換について成り立ちますが、LDAなどの近似交換汎関数では、一般に自己相互作用の誤差が残ります。）。

\(U_c[\rho_0]\)は残りの電子間反発の補正エネルギーで、実際の系の波動関数 \(\Psi_0\) を使った電子間反発エネルギーの期待値とKS軌道から作ったスレーター行列式 \(\Phi_s\) を使って計算した電子間反発エネルギーの期待値の差\[U_c[\rho_0]=\langle\Psi_0|\hat V_{ee}|\Psi_0\rangle-\langle\Phi_s|\hat V_{ee}|\Phi_s\rangle\]と定義されます。さらに、運動エネルギーの差も含めて、相関エネルギーを定義します。\[E_c[\rho_0] = \underbrace{\left(T[\rho_0]-T_s[\rho_0]\right)}_{\text{運動エネルギーの差}} + \underbrace{U_c[\rho_0]}_{\text{電子間反発エネルギーの差}}\]実際にできるかはともかく運動エネルギーの差を電子密度から計算する汎関数\(\Delta T[\rho_0]\)が理論上は存在します。ここで定義した \(E_c\) が、KS参照系を基準とする相関エネルギーです。KS行列式とHF行列式は一般に異なるため、WFTで用いられる電子相関エネルギーの定義 \(E_{\mathrm{exact}}-E_{\mathrm{HF}}\) とは、一般には一致しません。

以上から、電子系の全エネルギーは、\[\begin{align}E[\rho_0] &= T_s + E_{ext}+ J + E_{x} + E_c \\
&= T_s + E_{ext}+ J + E_{x} + \left(T-T_s\right) + (V_{ee} – J – E_{x})\\
&= T + E_{\mathrm{ext}} + V_{ee}.\end{align}\]となります。このように、運動エネルギーの差\(T – T_s\)と電子間相互作用の非古典的な効果を交換相関汎関数にまとめることで、元の電子系のエネルギーを厳密に書き直せます。

汎関数は基底状態以外の試行密度に対しても定義されています。\[E_{xc}[\rho] = F[\rho] – T_s[\rho] – J[\rho]\]\[E[\rho] = T_s[\rho]+ E_{\mathrm{ext}}[\rho] + J[\rho] + E_{xc}[\rho]\]

電子数を固定する条件 \(\int\rho(\mathbf r)\,d\mathbf r=N\) のもとで\(E[\rho]\)を変分すると、微分可能な場合、基底状態密度では\[\left.\frac{\delta E[\rho]}{\delta\rho(\mathbf r)}\right|_{\rho=\rho_0}=\mu\]となります。展開すると、\[\left.\frac{\delta T_s[\rho]}{\delta\rho(\mathbf r)}\right|_{\rho_0}+v_{\mathrm{ext}}(\mathbf r)+v_{\mathrm H}[\rho_0](\mathbf r)+v_{\mathrm{xc}}[\rho_0](\mathbf r)=\mu,\qquad v_{\mathrm{xc}}[\rho](\mathbf r)=\frac{\delta E_{\mathrm{xc}}[\rho]}{\delta\rho(\mathbf r)}.\] ここで、\(E_{\mathrm{xc}}=E_x+E_c\)、 \(\mu\) は電子数を固定するためのラグランジュ未定乗数です。

問題は、何度も「厳密な交換・相関汎関数を使えば」と書いている通りで、厳密な汎関数の具体的な形は一般に分かっていない点にあります。この記事では、基本的な近似である局所密度近似（LDA、各点の密度に、その密度の均一電子ガスの交換相関エネルギーを当てはめる近似）に基づく交換・相関汎関数を使いました。これらは、背景に一様な正電荷があって、その中を多数の電子が運動していると考える均一電子ガス（homogeneous electron gas）という理想的な系から導かれています。