できるだけ簡単にHartree-Fock法を解説してみる（C++17とJuliaのコード付き）

はじめに

この記事は，Physics Advent Calendar 2017 25日目の記事です。この記事では，できる限り簡単（？）に，Hartree-Fock（ハートリー・フォック）法を解説していきます。

量子力学における多体問題（量子多体問題）とHartree-Fock法

Physics Advent Calendar 2017 3日目の，adhara_mathphysさんによる非相対論的水素原子における力学的対称代数の記事にあるように，水素原子は高度な力学的対称性を有するため，水素原子に対する（非相対論的）Schrödinger方程式は何通りもの方法で，解析的に解くことができます（参考サイト[1]）。
しかし一般に，ヘリウム原子などの多電子原子，あるいは分子においては，この対称性が破れることによって可積分系ではなくなるため，これらの系に対するSchrödinger方程式を解析的に解くことは不可能です。これは，量子多体問題と呼ばれています（参考サイト[2]）。そこで，これらの系に対するSchrödinger方程式を解く（エネルギー固有状態を求める）ためには，何らかの近似を行わなければなりませんが，Hartree-Fock法はそのような近似の一つです。
なおこの記事では，（特殊）相対論効果，ラムシフト，超微細構造などについては無視します（これらを扱っていると，それだけで記事がいくつも必要になります）。

Hartree-Fock法のあらまし

Hartree-Fock法（またはHartree-Fock近似）とは，多電子原子，または分子に対するSchrödinger方程式を解くために，V. A. Fockによって考案された方法であり，その際に現れる方程式をHartree-Fock方程式と呼びます（参考サイト[3]）。
この方法は，N個の電子に対して適切に反対称化した行列積の波動関数を作り，相互作用を全て取り入れたHamiltonianに対して，全エネルギーを極小にする単一の行列積を探す，というものです。と言っても何のことか分からないと思いますが，要は，多体問題を一体問題に帰着させる方法です。つまり，（量子）多体問題が解けないなら，何らかの近似で，それを（比較的解きやすい）一体問題に置き換えちゃえばいいじゃない，ということです。

ヘリウム原子へのHartree-Fock法の「考え方」の適用

ヘリウム原子は，原子核1個と電子2個を有する系であり，最も単純な多電子原子です。Hartree-Fock法について一般的な議論をする前に，ヘリウム原子についてHartree-Fock法の「考え方」を適用することは有用でしょう。また，その際に現れる方程式の解き方についても，以下で示したいと思います。

Born-Oppenheimer近似

多電子原子や分子に対して，Hartree-Fock法を適用する際にまず重要なことは，（原子）核の取り扱いです。通常のHartree-Fock法においては，電子が核に相対的に運動している間は，核が「静止」していると見なし，電子と核の運動を分離します。つまり，（電子の運動について考えている間は）核の運動については一切考えません。これをBorn-Oppenheimer（ボルン・オッペンハイマー）近似と呼びます。断熱近似という用語もあり，しばしば，Born-Oppenheimer近似と同義として用いられますが，厳密には両者は異なります（参考サイト[4]）。なお，水素原子の場合には，電子と核の運動は厳密に分離できます（参考サイト[5]）。

ヘリウム原子に対するSchrödinger方程式

それでは早速，ヘリウム原子に対するSchrödinger方程式を書き下してみましょう。ヘリウム原子に対するSchrödinger方程式は，Born-Oppenheimer近似を用いると以下のように書けます（電子の質量$ m $とディラック定数$ \hbar $を1とする，Hartree原子単位系（参考サイト[6]）を用いていることに注意してください）。

\left[  -\dfrac{1}{2}\nabla_{1}^{2}-\dfrac{1}{2}\nabla_{2}^{2}-\dfrac
{2}{r_{1}}-\dfrac{2}{r_{2}}+\dfrac{1}{\left\vert \vec{r}_{1}-\vec{r}%
_{2}\right\vert }\right]  \psi\left(  \vec{x}_{1},\vec{x}_{2}\right)
=E\psi\left(  \vec{x}_{1},\vec{x}_{2}\right) \qquad(1)

ただし，座標$ \vec{x}_{1} $および$ \vec{x}_{2} $ は，空間位置座標$ \vec{r}_{i} $及びスピン座標$ s_{i} $をひとまとめにした$ \vec{x}_{i}=\left( \vec{r}_{i},s_{i}\right) $です。またここで，左辺括弧の中の第一項は電子1の運動エネルギーポテンシャル，同第二項は電子2の運動エネルギーポテンシャル，同第三項は核と電子1間のCoulombポテンシャル，同第四項は核と電子2間のCoulombポテンシャル，同第五項は電子1と電子2間のCoulombポテンシャル（この項が問題）です。
左辺括弧の中の第五項があるせいで，$ (1) $式を変数分離して解析的に解くのは不可能なのですが，ひとまずこの項を無視して，波動関数$ \psi\left( \vec{x}_{1},\vec{x}_{2}\right) $が以下のように変数分離できると仮定してみましょう。

\psi\left(  \vec{x}_{1},\vec{x}_{2}\right)  =\phi_{1}\left(  \vec{x}%
_{1}\right)  \phi_{2}\left(  \vec{x}_{2}\right) \qquad(2)

$ (2) $式はHartree積と呼ばれますが，これは，第一近似としては悪くありません（実際，この方法はHartree近似（参考サイト[7]）として知られています）。しかし$ (2) $式は，電子はフェルミオンであるので，波動関数は2個の座標$ \vec{x}_{1} $および$ \vec{x}_{2} $に対して反対称でなければならない（$ \vec{x}_{1} $と$ \vec{x}_{2} $を入れ替えたら，波動関数の符号が逆転しなければならない），という要請を満たしていません。$ (2) $式の右辺を，この要請を満たすように改良するのは簡単で，それは以下のように書けるでしょう。

\psi\left(  \vec{x}_{1},\vec{x}_{2}\right)  =\dfrac{1}{\sqrt{2}}\left[
\phi_{1}\left(  \vec{x}_{1}\right)  \phi_{2}\left(  \vec{x}_{2}\right)
-\phi_{2}\left(  \vec{x}_{1}\right)  \phi_{1}\left(  \vec{x}_{2}\right)
\right] \qquad(3)

$ (3) $式は，$ \vec{x}_{1} $と$ \vec{x}_{2} $を入れ替えると，波動関数の符号が逆転します（係数$ 1/\sqrt{2} $は，波動関数の規格化に必要です）。従って確かに，2個の座標$ \vec{x}_{1} $および$ \vec{x}_{2} $に対して反対称になっています。あるいは，ヘリウム原子では，軌道$ \phi $を，2個の電子が占有していることを考慮すれば，

\psi\left(  \vec{x}_{1},\vec{x}_{2}\right)  =\phi\left(  \vec{r}_{1}\right)
\phi\left(  \vec{r}_{2}\right)  \dfrac{1}{\sqrt{2}}\left[  \alpha\left(
s_{1}\right)  \beta\left(  s_{2}\right)  -\alpha\left(  s_{2}\right)
\beta\left(  s_{1}\right)  \right] \qquad(4)

と書くことも可能でしょう。ここで，$ \alpha\left( s\right) $は上向きスピン，$ \beta\left( s\right) $は下向きスピンの波動関数です。さて，$ (4) $式の（近似）波動関数を，$ (1) $式に代入してみましょう（なお，これは変分法であり，$ (4) $式の波動関数は試行関数に対応します。以下では，$ (4) $式の波動関数を試行関数と表記します）。ここで，ヘリウム原子に対するSchrödinger方程式のHamiltonian

\hat{H}=-\dfrac{1}{2}\nabla_{1}^{2}-\dfrac{1}{2}\nabla_{2}^{2}-\dfrac
{2}{r_{1}}-\dfrac{2}{r_{2}}+\dfrac{1}{\left\vert \vec{r}_{1}-\vec{r}%
_{2}\right\vert } \qquad(5)

は，スピンには全く作用しないので，スピンに依存する部分はSchrödinger方程式の両辺で落ちることを用いると，試行関数を$ (1) $式に代入した結果の式は，以下のようになります。

\left[  -\dfrac{1}{2}\nabla_{1}^{2}-\dfrac{1}{2}\nabla_{2}^{2}-\dfrac
{2}{r_{1}}-\dfrac{2}{r_{2}}+\dfrac{1}{\left\vert \vec{r}_{1}-\vec{r}%
_{2}\right\vert }\right]  \phi\left(  \vec{r}_{1}\right)  \phi\left(  \vec
{r}_{2}\right)  =E\phi\left(  \vec{r}_{1}\right)  \phi\left(  \vec{r}%
_{2}\right) \qquad(6)

より簡単な式を得るため，$ (6) $式の両辺に左から$ \phi^{\ast}\left( \vec{r}_{2}\right) $を乗じて，$ \vec{r}_{2} $について積分し，$ \vec{r}_{2} $に関する依存性を取り除くと，以下の式が得られます。

\left[  -\dfrac{1}{2}\nabla_{1}^{2}-\dfrac{2}{r_{1}}+%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}_{2}\left\vert \phi\left(  \vec{r}_{2}\right)  \right\vert ^{2}%
\dfrac{1}{\left\vert \vec{r}_{1}-\vec{r}_{2}\right\vert }\right]  \phi\left(
\vec{r}_{1}\right)  =E^{\prime}\phi\left(  \vec{r}_{1}\right) \qquad(7)

ここで，積分して定数になる項（すなわち$ \vec{r}_{1} $に依存しない項）は$ E^{\prime} $に含めました。すなわち，

E^{\prime}=E-%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}_{2}\phi^{\ast}\left(  \vec{r}_{2}\right)  \left(  -\dfrac{1}{2}%
\nabla_{2}^{2}-\dfrac{2}{r_{2}}\right)  \phi\left(  \vec{r}_{2}\right) \qquad(8)

です（ただし，$ \phi $は1に規格化されていることを用いました）。ここで重要なことは，$ (7) $式のSchrödinger方程式が，自己無撞着問題の形をしていることです。すなわち，$ \phi $はこのSchrödinger方程式の解ですが，このSchrödinger方程式はまた，左辺括弧の中の第三項が存在するために，$ \phi $自身にも依存しています。従って，このSchrödinger方程式は非線形偏微分方程式となっており，（すでに大胆な近似を行ったにもかかわらず）解析的には解けません。従って，何らかの方法で数値的に解かなければなりません。

ヘリウム原子に対するSchrödinger方程式の数値的解法と基底による解法

$ (7) $式のSchrödinger方程式は，差分法，あるいは有限要素法などで，完全に数値的に解くことも可能です。ここで，$ (7) $式は，ポテンシャルが球対称であるため，動径方向の常微分方程式と角度方向の偏微分方程式に厳密に分離できることを用いれば，$ (7) $式は，単なる常微分方程式に還元され，数値的に解くことがより容易になります（角度方向の偏微分方程式の解は，球面調和関数として解析的に得られます）。ヘリウム原子において，（動径方向の）常微分方程式を完全に数値的に解いたコードとしては，例えば拙作のschracがあります。
しかし一般に，Hartree-Fock方程式を完全に数値的に解くことは不可能です（参考サイト[3]，参考文献[1]）。従って，Hartree-Fock方程式を解く際は，専ら，基底で展開して行列方程式に変換して解くという方法が用いられます（この方法によって，行列方程式に変換されたHartree-Fock方程式は，Roothaan（ローターン）方程式と呼ばれます）。ここで$ (7) $式を，基底で展開して行列方程式に変換して解くことは，Roothaan方程式の一般的な議論をする前の導入として有用でしょう。従って以下では，この方法について説明します。

基底の導入

ヘリウム原子のことは一旦忘れ，基底での展開によって微分方程式が行列方程式に変換できることを，簡単な例を挙げて具体的に説明したいと思います。まず，波動関数$ \phi\left( \vec{r}\right) $が，以下のように2個の基底関数の線形結合で表せる（2個の基底で展開できる）と仮定してみましょう。

\phi\left(  \vec{r}\right)  =C_{1}\chi_{1}\left(  \vec{r}\right)  +C_{2}%
\chi_{2}\left(  \vec{r}\right) \qquad(9)

この$ \phi\left( \vec{r}\right) $を，Schrödinger方程式に代入すると，

\hat{H}\left[  C_{1}\chi_{1}\left(  \vec{r}\right)  +C_{2}\chi_{2}\left(  \vec
{r}\right)  \right]  =E\left[  C_{1}\chi_{1}\left(  \vec{r}\right)  +C_{2}%
\chi_{2}\left(  \vec{r}\right)  \right] \qquad(10)

となります。さらに，左から$ \phi^{\ast}\left( \vec{r}\right) $を乗じると，

\left[  C_{1}^{\ast}\chi_{1}^{\ast}\left(  \vec{r}\right)  +C_{2}^{\ast}%
\chi_{2}^{\ast}\left(  \vec{r}\right)  \right]  \hat{H}\left[  C_{1}\chi_{1}\left(
\vec{r}\right)  +C_{2}\chi_{2}\left(  \vec{r}\right)  \right]  =\left[
C_{1}^{\ast}\chi_{1}^{\ast}\left(  \vec{r}\right)  +C_{2}^{\ast}\chi_{2}%
^{\ast}\left(  \vec{r}\right)  \right]  E\left[  C_{1}\chi_{1}\left(  \vec
{r}\right)  +C_{2}\chi_{2}\left(  \vec{r}\right)  \right] \qquad(11)

となります。何も難しいことはありません。ここで，$ \chi_{1}\left( \vec{r}\right) $と，$ \chi_{2}\left( \vec{r}\right) $を実関数に限定すると，$ C_{1} $と$ C_{2} $も実数となり，従って$ (11) $式は，

\left[  C_{1}\chi_{1}\left(  \vec{r}\right)  +C_{2}\chi_{2}\left(  \vec
{r}\right)  \right]  \hat{H}\left[  C_{1}\chi_{1}\left(  \vec{r}\right)  +C_{2}%
\chi_{2}\left(  \vec{r}\right)  \right]  =\left[  C_{1}\chi_{1}\left(  \vec
{r}\right)  +C_{2}\chi_{2}\left(  \vec{r}\right)  \right]  E\left[  C_{1}%
\chi_{1}\left(  \vec{r}\right)  +C_{2}\chi_{2}\left(  \vec{r}\right)  \right] \qquad(12)

となります。これを一旦展開します。

C_{1}\chi_{1}\left(  \vec{r}\right)  \hat{H}C_{1}\chi_{1}\left(  \vec{r}\right)
+C_{1}\chi_{1}\left(  \vec{r}\right)  \hat{H}C_{2}\chi_{2}\left(  \vec{r}\right)
+C_{2}\chi_{2}\left(  \vec{r}\right)  \hat{H}C_{1}\chi_{1}\left(  \vec{r}\right)
+C_{2}\chi_{2}\left(  \vec{r}\right)  \hat{H}C_{2}\chi_{2}\left(  \vec{r}\right)= \\ 
C_{1}\chi_{1}\left(  \vec{r}\right)  EC_{1}\chi_{1}\left(  \vec{r}\right)
+C_{1}\chi_{1}\left(  \vec{r}\right)  EC_{2}\chi_{2}\left(  \vec{r}\right)
+C_{2}\chi_{2}\left(  \vec{r}\right)  EC_{1}\chi_{1}\left(  \vec{r}\right)
+C_{2}\chi_{2}\left(  \vec{r}\right)  EC_{2}\chi_{2}\left(  \vec{r}\right) \qquad(13)

さらにこれを全空間で積分すると，

%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\left[  C_{1}\chi_{1}\left(  \vec{r}\right)  HC_{1}\chi_{1}\left(
\vec{r}\right)  +C_{1}\chi_{1}\left(  \vec{r}\right)  HC_{2}\chi_{2}\left(
\vec{r}\right)  +C_{2}\chi_{2}\left(  \vec{r}\right)  HC_{1}\chi_{1}\left(
\vec{r}\right)  +C_{2}\chi_{2}\left(  \vec{r}\right)  HC_{2}\chi_{2}\left(
\vec{r}\right)  \right]  =%
\\ %TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\left[  C_{1}\chi_{1}\left(  \vec{r}\right)  EC_{1}\chi_{1}\left(
\vec{r}\right)  +C_{1}\chi_{1}\left(  \vec{r}\right)  EC_{2}\chi_{2}\left(
\vec{r}\right)  +C_{2}\chi_{2}\left(  \vec{r}\right)  EC_{1}\chi_{1}\left(
\vec{r}\right)  +C_{2}\chi_{2}\left(  \vec{r}\right)  EC_{2}\chi_{2}\left(
\vec{r}\right)  \right] \qquad(14)

となります。ごちゃごちゃしていますが，難しいことは何もしていません。ここで，$ (14) $式を以下のような「連立方程式」に書き直します。

\left\{
\begin{array}{ll}
C_{1}%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\left[  \chi_{1}\left(  \vec{r}\right)  H\chi_{1}\left(  \vec
{r}\right)  C_{1}+\chi_{1}\left(  \vec{r}\right)  H\chi_{2}\left(  \vec
{r}\right)  C_{2}\right]  =C_{1}E%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\left[  \chi_{1}\left(  \vec{r}\right)  \chi_{1}\left(  \vec
{r}\right)  C_{1}+\chi_{1}\left(  \vec{r}\right)  \chi_{2}\left(  \vec
{r}\right)  C_{2}\right] \\
C_{2}%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\left[  \chi_{2}\left(  \vec{r}\right)  H\chi_{1}\left(  \vec
{r}\right)  C_{1}+\chi_{2}\left(  \vec{r}\right)  H\chi_{2}\left(  \vec
{r}\right)  C_{2}\right]  =C_{2}E%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\left[  \chi_{2}\left(  \vec{r}\right)  \chi_{1}\left(  \vec
{r}\right)  C_{1}+\chi_{2}\left(  \vec{r}\right)  \chi_{2}\left(  \vec
{r}\right)  C_{2}\right]
\end{array}
\right. \qquad(15)

$ (15) $式が成り立つのなら，$ (14) $式が成立することは自明です。ここで，この連立方程式を$ C_{1} $と$ C_{2} $でそれぞれ除し，さらに行列形式で書き直すと，以下の式が得られます。

\left(
\begin{array}
[c]{cc}%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{1}\left(  \vec{r}\right)  H\chi_{1}\left(  \vec{r}\right)   &
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{1}\left(  \vec{r}\right)  H\chi_{2}\left(  \vec{r}\right)  \\%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{2}\left(  \vec{r}\right)  H\chi_{1}\left(  \vec{r}\right)   &
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{2}\left(  \vec{r}\right)  H\chi_{2}\left(  \vec{r}\right)
\end{array}
\right)  \left(
\begin{array}
[c]{c}%
C_{1}\\
C_{2}%
\end{array}
\right)  =E\left(
\begin{array}
[c]{cc}%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{1}\left(  \vec{r}\right)  \chi_{1}\left(  \vec{r}\right)   &
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{1}\left(  \vec{r}\right)  \chi_{2}\left(  \vec{r}\right)  \\%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{2}\left(  \vec{r}\right)  \chi_{1}\left(  \vec{r}\right)   &
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{2}\left(  \vec{r}\right)  \chi_{2}\left(  \vec{r}\right)
\end{array}
\right)  \left(
\begin{array}
[c]{c}%
C_{1}\\
C_{2}%
\end{array}
\right) \qquad(16)

ここで，

\mathbf{H}=\left(
\begin{array}
[c]{cc}%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{1}\left(  \vec{r}\right)  H\chi_{1}\left(  \vec{r}\right)   &
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{1}\left(  \vec{r}\right)  H\chi_{2}\left(  \vec{r}\right)  \\%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{2}\left(  \vec{r}\right)  H\chi_{1}\left(  \vec{r}\right)   &
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{2}\left(  \vec{r}\right)  H\chi_{2}\left(  \vec{r}\right)
\end{array}
\right) \quad(17a) \\
\mathbf{S}=\left(
\begin{array}
[c]{cc}%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{1}\left(  \vec{r}\right)  \chi_{1}\left(  \vec{r}\right)   &
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{1}\left(  \vec{r}\right)  \chi_{2}\left(  \vec{r}\right)  \\%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{2}\left(  \vec{r}\right)  \chi_{1}\left(  \vec{r}\right)   &
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{2}\left(  \vec{r}\right)  \chi_{2}\left(  \vec{r}\right)
\end{array}
\right) \quad(17b) \\
\vec{C}=\left(
\begin{array}
[c]{c}%
C_{1}\\
C_{2}%
\end{array}
\right) \quad(17c)

とおくと（$ \mathbf{S} $は重なり行列と呼ばれます），

\mathbf{H}\vec{C}=E\mathbf{S}\vec{C} \qquad(18)

なる行列方程式が得られます。$ (18) $式は一般化固有値問題となっており，$ E $が（エネルギー）固有値で，$ \vec{C} $が固有ベクトルです（なお，基底が正規直交系であれば，$ \mathbf{S} $は単位行列となり，$ (18) $式は通常の固有値問題になります）。さてこれで，当初の目的である，「微分方程式を行列方程式に変換する」ことが達成できました。つまり，波動関数$ \phi\left( \vec{r}\right) $を，2個の基底で展開すれば，微分方程式から行列方程式への変換が可能だということが分かったわけです。
ここで一般の場合，つまり波動関数を任意の個数（$ n $個）の基底で展開した場合を考えます。この場合，波動関数$ \phi\left( \vec{r}\right) $は以下のように書けるでしょう。

\phi\left(  \vec{r}\right)  =%
%TCIMACRO{\dsum \limits_{p=1}^{n}}%
%BeginExpansion
{\displaystyle\sum\limits_{p=1}^{n}}
%EndExpansion
C_{p}\chi_{p}\left(  \vec{r}\right) \quad(19)

この場合も，全く同様の手順によって$ (18) $式が得られます（ただし，行列$ \mathbf{H} $および重なり行列$ \mathbf{S} $が，$ n $次正方行列となることと，固有ベクトル$ \vec{C} $が$ n $次元ベクトルとなることに注意してください）。

基底関数の種類

前節で，波動関数を任意の個数の基底で展開することで，微分方程式を行列方程式に変換できることを見ました。しかし，肝心の基底関数には一体どのような形の関数を選べば良いのでしょうか。結論から言うと，Hartree-Fock法で最も多用されている基底関数は，ガウス関数（$ e^{-ar^{2}} $）であり，これはGTO（Gaussian Type Orbital，ガウス型軌道）と呼ばれています（これが多用されている理由は後述します）。また，$ e^{-ar} $型の関数も用いられており，これはSTO（Slater Type Orbital，スレーター型軌道）と呼ばれます。他には，Bessel関数や，B-splineなどが基底関数として用いられることもあります。今は（暗黙に）非周期的かつ有限の系を考えているのですが，平面波（平たく言えば三角関数）を基底関数とすることも可能のようです（参考文献[3]）。

(7)式の基底による展開

ヘリウム原子の話に戻って，$ (7) $式の$ \phi\left( \vec{r}_{1}\right) $と$ \phi\left( \vec{r}_{2}\right) $を，以下のようにn個の基底で展開してみましょう。

\phi\left(  \vec{r}_{1}\right)  =%
%TCIMACRO{\dsum \limits_{q=1}^{n}}%
%BeginExpansion
{\displaystyle\sum\limits_{q=1}^{n}}
%EndExpansion
C_{q}\chi_{q}\left(  \vec{r}_{1}\right) \quad(20) \\
\phi\left(  \vec{r}_{2}\right)  =%
%TCIMACRO{\dsum \limits_{u=1}^{n}}%
%BeginExpansion
{\displaystyle\sum\limits_{u=1}^{n}}
%EndExpansion
C_{u}\chi_{u}\left(  \vec{r}_{2}\right) \quad(21)

$ (7) $式に，$ (20) $式および$ (21) $式を代入すると，以下の式が得られます。

\left[  -\dfrac{1}{2}\nabla_{1}^{2}-\dfrac{2}{r_{1}}+%
%TCIMACRO{\dsum \limits_{r,s=1}^{n}}%
%BeginExpansion
{\displaystyle\sum\limits_{r,s=1}^{n}}
%EndExpansion
C_{r}C_{s}%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}_{2}\chi_{r}\left(  \vec{r}_{2}\right)  \chi_{s}\left(  \vec{r}%
_{2}\right)  \dfrac{1}{\left\vert \vec{r}_{1}-\vec{r}_{2}\right\vert }\right]
%
%TCIMACRO{\dsum \limits_{q=1}^{n}}%
%BeginExpansion
{\displaystyle\sum\limits_{q=1}^{n}}
%EndExpansion
C_{q}\chi_{q}\left(  \vec{r}_{1}\right)  =E^{\prime}%
%TCIMACRO{\dsum \limits_{q=1}^{n}}%
%BeginExpansion
{\displaystyle\sum\limits_{q=1}^{n}}
%EndExpansion
C_{q}\chi_{q}\left(  \vec{r}_{1}\right) \quad(22)

さらに$ (22) $式に，左から$ \chi_{p}\left( \vec{r}_{1}\right) $を乗じ，$ \vec{r}_{1} $について全空間で積分すると，以下の式が得られます。

%
%TCIMACRO{\dsum \limits_{p,q=1}^{n}}%
%BeginExpansion
{\displaystyle\sum\limits_{p,q=1}^{n}}
%EndExpansion
\left(  h_{pq}+%
%TCIMACRO{\dsum \limits_{r,s=1}^{n}}%
%BeginExpansion
{\displaystyle\sum\limits_{r,s=1}^{n}}
%EndExpansion
C_{r}C_{s}Q_{prqs}\right)  C_{q}=E^{\prime}%
%TCIMACRO{\dsum \limits_{p,q=1}^{n}}%
%BeginExpansion
{\displaystyle\sum\limits_{p,q=1}^{n}}
%EndExpansion
S_{pq}C_{q} \quad(23)

ただし，

h_{pq}=%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}_{1}\chi_{p}\left(  \vec{r}_{1}\right)  \left(  -\dfrac{1}{2}%
\nabla_{1}^{2}-\dfrac{2}{r_{1}}\right)  \chi_{q}\left(  \vec{r}_{1}\right) \quad(24a) \\
Q_{prqs}=%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}_{1}d\vec{r}_{2}\chi_{p}\left(  \vec{r}_{1}\right)  \chi_{r}\left(
\vec{r}_{2}\right)  \dfrac{1}{\left\vert \vec{r}_{1}-\vec{r}_{2}\right\vert
}\chi_{q}\left(  \vec{r}_{1}\right)  \chi_{s}\left(  \vec{r}_{2}\right) \quad(24b) \\
S_{pq}=%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}_{1}\chi_{p}\left(  \vec{r}_{1}\right)  \chi_{q}\left(  \vec{r}%
_{1}\right) \quad(24c)

です。なお，$ (24a) $式は1電子積分，$ (24b) $式は2電子積分，$ (24c) $式は重なり積分と呼ばれます。さらに，

F_{pq}=h_{pq}+%
%TCIMACRO{\dsum \limits_{r,s=1}^{n}}%
%BeginExpansion
{\displaystyle\sum\limits_{r,s=1}^{n}}
%EndExpansion
C_{r}C_{s}Q_{prqs} \quad(25)

とおくと（$ F_{pq} $を行列と見なした$ \mathbf{F} $は，Fock行列と呼ばれます），結局$ (7) $式は，

\mathbf{F}\vec{C}=E^{\prime}\mathbf{S}\vec{C} \quad(26)

なる，行列方程式（一般化固有値問題）に帰着できます（やや複雑ですが，やっていることは「基底の導入」の節で説明したことと全く同じです）。ここまでくれば，あとは具体的に積分を計算して，行列要素を埋め，一般化固有値問題を数値的に解くだけです。ここで，「具体的に積分を計算して，行列要素を埋める」ためには，具体的な基底関数の形を決めなければなりませんが，水素原子の1s軌道の波動関数が$ e^{-ar} $の形，つまりSTOの形をしていたことを思い出せば，ヘリウム原子の波動関数も，基底にSTOを用いて展開できると考えるのが自然でしょう。
が，基底にSTOを用いるのには問題があります。基底にSTOを使うと，一般的には2電子積分が解析的に求められないのです（参考文献[4]によると，ヘリウム原子の場合は例外で，解析的に求められるらしいです）。一方，基底にGTOを使うと，2電子積分を解析的に求めることが可能になります。これが，Hartree-Fock法で基底としてGTOが多用されている理由です。

GTOで1電子積分，2電子積分および重なり積分を求める

基底にGTOを用いると，基底状態については，$ l = 0 $の関数（s関数）だけが必要であり，これは，原子核に中心がある（そのためそこを原点に置く）ので，$ \chi_{p}\left( \vec{r}_{1}\right) $は，以下のようになります。

\chi_{p}\left(  \vec{r}_{1}\right)  =e^{-\alpha_{p}r^{2}} \quad(27)

$ (27) $式を用いると，ヘリウム原子における，1電子積分$ h_{pq} $と2電子積分$ Q_{prqs} $，および重なり積分$ \mathbf{S} $の要素$ S_{pq} $は，以下のように求められます（計算は省略します）。

h_{pq}=3\dfrac{\alpha_{p}\alpha_{q}\pi^{\frac{3}{2}}}{\left(  \alpha
_{p}+\alpha_{q}\right)  ^{\frac{5}{2}}}-\dfrac{4\pi}{\alpha_{p}+\alpha_{q}} \quad(28) \\
Q_{prqs}=\dfrac{2\pi^{\frac{5}{2}}}{\left(  \alpha_{p}+\alpha_{q}\right)
\left(  \alpha_{r}+\alpha_{s}\right)  \sqrt{\alpha_{p}+\alpha_{q}+\alpha
_{r}+\alpha_{s}}} \quad(29) \\
S_{pq}=\left(  \dfrac{\pi}{\alpha_{p}+\alpha_{q}}\right)  ^{\frac{3}{2}} \quad(30)

ここで，具体的な$ \alpha_{p} $の値についてですが，参考文献[2]に，$ n = 4 $の場合（$ \phi\left( \vec{r}_{1}\right) $と$ \phi\left( \vec{r}_{2}\right) $を4つのGTOで展開した場合）の値が記載されているので，それをそのまま引用します（なお一般に$ \alpha_{p} $の値は，非線形の最適化問題を解くことによって求められますが，ここではその方法については立ち入りません）。その値は，

\begin{array}{ll}
\alpha_{1}=0.297104 \\
\alpha_{2}=1.236745 \\
\alpha_{3}=5.749982 \\
\alpha_{4}=38.216677
\end{array} \qquad(31)

です。

自己無撞着問題の解法

前節で，行列要素を埋めるための情報が全て得られたので，後は$ (26) $式の一般化固有値問題を数値的に解くだけ…と言いたいところですが，そうは問屋が卸しません。$ (25) $式の第二項を見てください。この項の$ C_{r} $と$ C_{s} $は，$ (26) $式の一般化固有値問題の解である，固有ベクトル$ \vec{C} $の要素です。つまり，問題の「答え」が，問題を解くために必要なのです。これは，「ヘリウム原子に対するSchrödinger方程式」の節で述べたように，$ (7) $式が自己無撞着問題だったからです。じゃあどうすんねん，と突っ込みたくなりますが，一応解決法はあります。どうするかというと，以下のような手続きを踏みます。

1. 固有ベクトル$ \vec{C} $の要素$ C_{p} $の初期値を適当に選びます。これは例えば，全部$ 0 $とか，あるいは全部$ 1 $とかで構いません（今はヘリウム原子という単純な系を考えているので，たいていの初期値から正しい解に収束するはずです）。また，1電子積分$ h_{pq} $と2電子積分$ Q_{prqs} $，および重なり積分$ \mathbf{S} $の要素$ S_{pq} $も計算しておきます。
2. この$ \vec{C} $から，$ (25) $式によって，Fock行列$ \mathbf{F} $を計算します。
3. $ (26) $式の一般化固有値問題を数値的に解きます。
4. エネルギー

E=2%
%TCIMACRO{\dsum \limits_{p,q=1}^{n}}%
%BeginExpansion
{\displaystyle\sum\limits_{p,q=1}^{n}}
%EndExpansion
C_{p}C_{q}h_{pq}+%
%TCIMACRO{\dsum \limits_{p,q,r,s=1}^{n}}%
%BeginExpansion
{\displaystyle\sum\limits_{p,q,r,s=1}^{n}}
%EndExpansion
C_{p}C_{q}C_{r}C_{s}Q_{prqs}=E^{\prime}+%
%TCIMACRO{\dsum \limits_{p,q=1}^{n}}%
%BeginExpansion
{\displaystyle\sum\limits_{p,q=1}^{n}}
%EndExpansion
C_{p}C_{q}h_{pq} \qquad(32)

を求めます（$ E $と$ E^{\prime} $の関係については，$ (8) $式を思い出してください。ここで，

h_{pq}=%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}_{1}\chi_{p}\left(  \vec{r}_{1}\right)  \left(  -\dfrac{1}{2}%
\nabla_{1}^{2}-\dfrac{2}{r_{1}}\right)  \chi_{q}\left(  \vec{r}_{1}\right)  =%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}_{2}\chi_{p}\left(  \vec{r}_{2}\right)  \left(  -\dfrac{1}{2}%
\nabla_{2}^{2}-\dfrac{2}{r_{2}}\right)  \chi_{q}\left(  \vec{r}_{2}\right) \qquad(33)

は，系の対称性から自明です）。
5. ステップ3で求まった固有ベクトル$ \vec{C} $を用いて，$ (25) $式からFock行列$ \mathbf{F} $を計算します。
6. ステップ3～5を繰り返します。

ここで，このままでは無限ループになってしまうので，どこかで計算を打ち切らなくてはなりません。通常は，前回計算したエネルギー$ E_{old} $と，今回計算したエネルギー$ E_{new} $の差の絶対値$ \left\vert E_{new}-E_{old}\right\vert $が，ある閾値未満になったところで計算を打ち切ります。なお，この計算法を，SCF（Self-Consistent Field，自己無撞着場）計算と呼びます。

ヘリウム原子のエネルギーを求めるC++のプログラム

上述の方法によって，ヘリウム原子の（基底状態の）エネルギーを求める，C++17のプログラムは以下のようになります（多次元配列のためにBoost C++ Librariesを，また一般化固有値問題の解を求めるためにEigenをそれぞれ使っていますので，注意してください）。なお，このプログラムの完全なソースコードは，GitHubのこちらにあります。

/*! \file helium_hf.cpp
    \brief Hartree-Fock法で、ヘリウム原子のエネルギーを計算する
    Copyright © 2017 @dc1394 All Rights Reserved.
    (but this is originally adapted by Paolo Giannozzi for helium_hf_gauss.c from http://www.fisica.uniud.it/~giannozz/Corsi/MQ/Software/C/helium_hf_gauss.c )
    This software is released under the BSD 2-Clause License.
*/

#include <cmath>                                // for std::pow, std::sqrt
#include <cstdint>                              // for std::int32_t
#include <iostream>                             // for std::cerr, std::cin, std::cout
#include <optional>                             // for std::make_optional, std::optional, std::nullopt
#include <valarray>                             // for std::valarray
#include <boost/assert.hpp>                     // for BOOST_ASSERT
#include <boost/format.hpp>                     // for boost::format
#include <boost/math/constants/constants.hpp>   // for boost::math::constants::pi
#include <boost/multi_array.hpp>                // for boost::multi_array
#include <Eigen/Core>                           // for Eigen::MatrixXd, Eigen::VectorXd
#include <Eigen/Eigenvalues>                    // for Eigen::GeneralizedSelfAdjointEigenSolver

namespace {
    //! A global variable (constant expression).
    /*!
        バッファサイズの上限
    */
    static auto constexpr MAXBUFSIZE = 32;

    //! A global variable (constant expression).
    /*!
        SCF計算のループの上限
    */
    static auto constexpr MAXITER = 1000;

    //! A global variable (constant expression).
    /*!
        SCF計算のループから抜ける際のエネルギーの差の閾値
    */
    static auto constexpr SCFTHRESHOLD = 1.0E-15;

    //! A global function.
    /*!
        SCF計算を行う
        \return SCF計算が正常に終了した場合はエネルギーを、しなかった場合はstd::nulloptを返す
    */
    std::optional<double> do_scfloop();

    //! A global function.
    /*!
        nalpha個のGTOによるヘリウム原子のエネルギーを計算する
        \param c 固有ベクトルC
        \param ep 一般化固有値問題のエネルギー固有値E'
        \param h 1電子積分
        \return ヘリウム原子のエネルギー
    */
    double getenergy(Eigen::VectorXd const & c, double ep, boost::multi_array<double, 2> const & h);

    //! A global function.
    /*!
        使用するGTOの数をユーザに入力させる
        \return 使用するGTOの数
    */
    std::int32_t input_nalpha();

    //! A global function.
    /*!
        GTOの肩の係数が格納された配列を生成する
        \param nalpha 使用するGTOの個数
        \return GTOの肩の係数が格納されたstd::vector
    */
    std::valarray<double> make_alpha(std::int32_t nalpha);

    //! A global function.
    /*!
        全ての要素が、引数で指定された値で埋められたnalpha次元ベクトルを生成する
        \param nalpha 使用するGTOの個数
        \param val 要素を埋める値
        \return 引数で指定された値で埋められたベクトル (Eigen::VectorXd)
    */
    Eigen::VectorXd make_c(std::int32_t nalpha, double val);

    //! A global function.
    /*!
        nalphaの数で、固有ベクトル、1電子積分および2電子積分からFock行列を生成する
        \param c 固有ベクトルC
        \param h 1電子積分hpq
        \param q 2電子積分Qprqs
        \return Fock行列 (Eigen::MatrixXd)
    */
    Eigen::MatrixXd make_fockmatrix(Eigen::VectorXd const & c, boost::multi_array<double, 2> const & h, boost::multi_array<double, 4> const & q);

    //! A global function.
    /*!
        1電子積分が格納された、nalpha×nalphaの2次元配列を生成する
        \param alpha GTOの肩の係数が格納されたstd::vector
        \return 1電子積分が格納された2次元配列 (boost::multi_array)
    */
    boost::multi_array<double, 2> make_oneelectroninteg(std::valarray<double> const & alpha);

    //! A global function.
    /*!
        nalpha次正方行列の重なり行列を生成する
        \param alpha GTOの肩の係数が格納されたstd::vector
        \return 重なり行列 (Eigen::MatrixXd)
    */
    Eigen::MatrixXd make_overlapmatrix(std::valarray<double> const & alpha);

    //! A global function.
    /*!
        2電子積分が格納されたnalpha×nalpha×nalpha×nalphaの4次元配列を生成する
        \param alpha GTOの肩の係数が格納されたstd::vector
        \return 2電子積分が格納された4次元配列 (boost::multi_array)
    */
    boost::multi_array<double, 4> make_twoelectroninteg(std::valarray<double> const & alpha);
}

int main()
{
    if (auto const res(do_scfloop()); res) {
        std::cout << boost::format("SCF計算が収束しました: energy = %.14f (Hartree)") % (*res) << std::endl;

        return 0;
    }
    else {
        std::cerr << "SCF計算が収束しませんでした" << std::endl;

        return -1;
    }

}

namespace {
    std::optional<double> do_scfloop()
    {
        // 使用するGTOの数を入力
        auto const nalpha(input_nalpha());

        // GTOの肩の係数が格納された配列を生成
        auto alpha = make_alpha(nalpha);

        // 1電子積分が格納された2次元配列を生成
        auto const h(make_oneelectroninteg(alpha));

        // 2電子積分が格納された4次元配列を生成
        auto const q(make_twoelectroninteg(alpha));

        // 重なり行列を生成
        auto const s(make_overlapmatrix(alpha));

        // 全て0.0で初期化された固有ベクトルを生成
        auto c(make_c(nalpha, 0.0));

        // 新しく計算されたエネルギー
        auto enew = 0.0;

        // SCFループ
        for (auto iter = 1; iter < MAXITER; iter++) {
            // Fock行列を生成
            auto const f(make_fockmatrix(c, h, q));

            // 一般化固有値問題を解く
            Eigen::GeneralizedSelfAdjointEigenSolver<Eigen::MatrixXd> es(f, s);

            // E'を取得
            auto const ep = es.eigenvalues()[0];

            // 固有ベクトルを取得
            c = es.eigenvectors().col(0);

            // 前回のSCF計算のエネルギーを保管
            auto const eold = enew;

            // 今回のSCF計算のエネルギーを計算する
            enew = getenergy(c, ep, h);

            std::cout << boost::format("Iteration # %2d: HF eigenvalue = %.14f, energy = %.14f\n") % iter % ep % enew;

            // SCF計算が収束したかどうか
            if (std::fabs(enew - eold) < SCFTHRESHOLD) {
                // 収束したのでそのエネルギーを返す
                return std::make_optional(enew);
            }
        }

        // SCF計算が収束しなかった
        return std::nullopt;
    }

    double getenergy(Eigen::VectorXd const & c, double ep, boost::multi_array<double, 2> const & h)
    {
        auto const nalpha = static_cast<std::int32_t>(c.size());
        auto e = ep;
        for (auto p = 0; p < nalpha; p++) {
            for (auto q = 0; q < nalpha; q++) {
                // E = E' + Cp * Cq * hpq
                e += c[p] * c[q] * h[p][q];
            }
        }

        return e;
    }

    std::int32_t input_nalpha()
    {
        std::int32_t nalpha;

        while (true) {
            std::cout << "使用するGTOの個数を入力してください (3, 4 or 6): ";
            std::cin >> nalpha;

            if (!std::cin.fail() && (nalpha == 3 || nalpha == 4 || nalpha == 6)) {
                break;
            }

            std::cin.clear();
            std::cin.ignore(MAXBUFSIZE, '\n');
        }

        return nalpha;
    }

    std::valarray<double> make_alpha(std::int32_t nalpha)
    {
        switch (nalpha) {
        case 3:
            return { 0.31364978999999998, 1.1589229999999999, 6.3624213899999997 };
            break;

        case 4:
            return { 0.297104, 1.236745, 5.749982, 38.2166777 };
            break;

        case 6:
            return { 0.18595935599999999, 0.45151632200000003, 1.1627151630000001, 3.384639924, 12.09819836, 65.984568240000002 };
            break;

        default:
            BOOST_ASSERT(!"switch文のdefaultに来てしまった！");
            return std::valarray<double>();
            break;
        }
    }

    Eigen::VectorXd make_c(std::int32_t nalpha, double val)
    {
        Eigen::VectorXd c(nalpha);

        // 固有ベクトルCの要素を全てvalで初期化
        for (auto i = 0; i < nalpha; i++) {
            c[i] = val;
        }

        return c;
    }

    Eigen::MatrixXd make_fockmatrix(Eigen::VectorXd const & c, boost::multi_array<double, 2> const & h, boost::multi_array<double, 4> const & q)
    {
        auto const nalpha = static_cast<std::int32_t>(c.size());
        Eigen::MatrixXd f = Eigen::MatrixXd::Zero(nalpha, nalpha);

        for (auto p = 0; p < nalpha; p++) {
            for (auto qi = 0; qi < nalpha; qi++) {
                // Fpq = hpq + ΣCr * Cs * Qprqs
                f(p, qi) = h[p][qi];

                for (auto r = 0; r < nalpha; r++) {
                    for (auto s = 0; s < nalpha; s++) {
                        f(p, qi) += c[r] * c[s] * q[p][r][qi][s];
                    }
                }
            }
        }

        return f;
    }

    boost::multi_array<double, 2> make_oneelectroninteg(std::valarray<double> const & alpha)
    {
        using namespace boost::math::constants;

        auto const nalpha = static_cast<std::int32_t>(alpha.size());
        boost::multi_array<double, 2> h(boost::extents[nalpha][nalpha]);

        for (auto p = 0; p < nalpha; p++) {
            for (auto q = 0; q < nalpha; q++) {
                // αp + αq
                auto const appaq = alpha[p] + alpha[q];

                // hpq = 3αpαqπ^1.5 / (αp + αq)^2.5 - 4π / (αp + αq)
                h[p][q] = 3.0 * alpha[p] * alpha[q] * std::pow((pi<double>() / appaq), 1.5) / appaq -
                          4.0 * pi<double>() / appaq;
            }
        }

        return h;
    }

    Eigen::MatrixXd make_overlapmatrix(std::valarray<double> const & alpha)
    {
        using namespace boost::math::constants;

        auto const nalpha = static_cast<std::int32_t>(alpha.size());
        Eigen::MatrixXd s = Eigen::MatrixXd::Zero(nalpha, nalpha);

        for (auto p = 0; p < nalpha; p++) {
            for (auto q = 0; q < nalpha; q++) {
                // Spq = (π / (αp + αq))^1.5
                s(p, q) = std::pow((pi<double>() / (alpha[p] + alpha[q])), 1.5);
            }
        }

        return s;
    }

    boost::multi_array<double, 4> make_twoelectroninteg(std::valarray<double> const & alpha)
    {
        using namespace boost::math::constants;

        auto const nalpha = static_cast<std::int32_t>(alpha.size());
        boost::multi_array<double, 4> q(boost::extents[nalpha][nalpha][nalpha][nalpha]);

        for (auto p = 0; p < nalpha; p++) {
            for (auto qi = 0; qi < nalpha; qi++) {
                for (auto r = 0; r < nalpha; r++) {
                    for (auto s = 0; s < nalpha; s++) {
                        // Qprqs = 2π^2.5 / [(αp + αq)(αr + αs)√(αp + αq + αr + αs)]
                        q[p][r][qi][s] = 2.0 * std::pow(pi<double>(), 2.5) /
                            ((alpha[p] + alpha[qi]) * (alpha[r] + alpha[s]) *
                            std::sqrt(alpha[p] + alpha[qi] + alpha[r] + alpha[s]));
                    }
                }
            }
        }

        return q;
    }
}

なお，このプログラムは，参考サイト[8]のhelium_hf_gauss.cを参考にさせて頂きました。また，$ n = 3 $および$ n = 6 $の場合の$ \alpha_{p} $の値は，参考サイト[9]の，basis-sto3g.cppおよびbasis-sto6g.cppから引用させて頂きました。

ヘリウム原子のエネルギーを求めるJuliaのプログラム

上述の方法によって，ヘリウム原子の（基底状態の）エネルギーを求める，Juliaのプログラムは以下のようになります。なお，このプログラムの完全なソースコードは，GitHubのこちらにあります。

using Match
using Printf

const DR = 0.01
const MAXR = 30.0
const MAXBUFSIZE = 32
const MAXITER = 1000
const SCFTHRESHOLD = 1.0E-15

function do_scfloop(verbose=true)
    # 使用するGTOの数を入力
    nalpha = input_nalpha()

    # GTOの肩の係数が格納された配列を生成
    alpha = make_alpha(nalpha)

    # 1電子積分が格納された2次元配列を生成
    h = make_oneelectroninteg(alpha)

    # 2電子積分が格納された4次元配列を生成
    q = make_twoelectroninteg(alpha)

    # 重なり行列を生成
    s = make_overlapmatrix(alpha)

    # 全て0.0で初期化された固有ベクトルを生成
    c = make_c(nalpha, 0.0)

    # 新しく計算されたエネルギー
    enew = 0.0

    # SCFループ
    f = Symmetric(similar(h))
    stmp = similar(s)

    for iter = 1:MAXITER
        # Fock行列を生成
        make_fockmatrix!(f, c, h, q)

        # 一般化固有値問題を解く
        stmp .= s
        eigenval, eigenvec = eigen!(f, stmp)

        # E'を取得
        ep = eigenval[1]

        # 固有ベクトルを取得
        c .= @view(eigenvec[:,1])

        # 固有ベクトルを正規化
        normalize!(alpha, c)

        # 前回のSCF計算のエネルギーを保管
        eold = enew

        # 今回のSCF計算のエネルギーを計算する
        enew = getenergy(c, ep, h)

        verbose && @printf "Iteration # %2d: HF eigenvalue = %.14f, energy = %.14f\n" iter ep enew

        # SCF計算が収束したかどうか
        if abs(enew - eold) < SCFTHRESHOLD
            # 収束したのでα, c, エネルギーを返す
            return alpha, c, enew
        end
    end     

    # SCF計算が収束しなかった
    return nothing
end

function getenergy(c, ep, h)
    nalpha = length(c)
    e = ep

    @inbounds @simd for q = 1:nalpha
        for p = 1:nalpha
            # E = E' + Cp * Cq * hpq
            e += c[p] * c[q] * h[p, q]
        end
    end

    return e
end

function input_nalpha()
    nalpha = nothing

    while true
        @printf "使用するGTOの個数を入力してください (3, 4 or 6): "
        s = readline()

        nalpha = tryparse(Int64, s)
        if nalpha != nothing
            @match nalpha begin
                3 => break
                4 => break
                6 => break
                _ => continue
            end
        end
    end

    return nalpha
end

function make_alpha(nalpha)
    return @match nalpha begin
                3 => [0.31364978999999998, 1.1589229999999999, 6.3624213899999997]
                4 => [0.297104, 1.236745, 5.749982, 38.2166777]
                6 => [0.18595935599999999, 0.45151632200000003, 1.1627151630000001, 3.384639924, 12.09819836, 65.984568240000002]
                _ => []
           end
end

function make_c(nalpha, val)
    return fill(val, nalpha)
end

function make_fockmatrix!(f, c, h, q)
    nalpha = length(c)
    f .= 0

    @inbounds for qi = 1:nalpha
        for p = 1:qi
            # Fpq = hpq + ΣCr * Cs * Qprqs
            f.data[p, qi] = h[p, qi]

            for s = 1:nalpha
                for r = 1:nalpha
                    f.data[p, qi] += c[r] * c[s] * q[p, r, qi, s]
                end
            end
        end
    end
end

function make_oneelectroninteg(alpha)
    nalpha = length(alpha)
    h = Symmetric(zeros(nalpha, nalpha))

    @inbounds for q = 1:nalpha
        for p = 1:q
            # αp + αq
            appaq = alpha[p] + alpha[q]

            # hpq = 3αpαqπ^1.5 / (αp + αq)^2.5 - 4π / (αp + αq)
            h.data[p, q] = 3.0 * alpha[p] * alpha[q] * ((pi / appaq) ^ 1.5) / appaq -
                    4.0 * pi / appaq
        end
    end

    return h
end

function make_overlapmatrix(alpha)
    nalpha = length(alpha)
    s = Symmetric(zeros(nalpha, nalpha))

    @inbounds for q = 1:nalpha
        for p = 1:q
            # Spq = (π / (αp + αq))^1.5
            s.data[p, q] = (pi / (alpha[p] + alpha[q])) ^ 1.5
        end
    end

    return s
end

function make_twoelectroninteg(alpha)
    nalpha = length(alpha)
    q = Array{Float64}(undef, nalpha, nalpha, nalpha, nalpha)

    @inbounds for p = 1:nalpha
        for qi = 1:nalpha
            for r = 1:nalpha
                for s = 1:nalpha
                    # Qprqs = 2π^2.5 / [(αp + αq)(αr + αs)√(αp + αq + αr + αs)]
                    q[p, r, qi, s] = 2.0 * (pi ^ 2.5) /
                        ((alpha[p] + alpha[qi]) * (alpha[r] + alpha[s]) *
                        sqrt(alpha[p] + alpha[qi] + alpha[r] + alpha[s]))
                end
            end
        end
    end

    return q
end

function normalize!(alpha, c)
    nalpha = length(c)

    sum = 0.0
    @inbounds @simd for i = 1:nalpha
        for j = 1:nalpha
            sum += c[i] * c[j] / (4.0 * (alpha[i] + alpha[j])) * (pi / (alpha[i] + alpha[j])) ^ 0.5
        end
    end

    c ./= sqrt(4.0 * pi * sum)
end

プログラムの実行結果

前節で紹介したプログラムを実行した結果は，以下の画像のようになります。

$ n = 4 $で，固有ベクトル$ \vec{C} $の要素$ C_{p} $の初期値がすべて0の場合，$ \left\vert E_{new}-E_{old}\right\vert $が，閾値（今の場合$ 10^{-15} $）未満になるまでに，29回のループを要していることが分かります。

プログラムによるヘリウム原子のエネルギーの計算結果の結果と考察

ここで，このプログラムによる$ n = 3 $，$ n = 4 $，$ n = 6 $のときのエネルギーの計算結果を表に示しますと，以下のようになります（ついでに，Hartree-Fock極限（この値は参考文献[5]より引用）と，厳密な値（この値は参考文献[6]より引用）も併せて，以下の表に示します）。

n	計算結果 (Hartree)
3	-2.8162
4	-2.8552
6	-2.8600
Hartree-Fock極限	−2.8617
厳密	−2.9037

上の表のように，$ n $が大きくなるほど，Hartree-Fock極限の値に近づくことが分かります。ここで，Hartree-Fock極限というのは，おおざっぱに言うと，$ n $を無限に大きくしていったときに収束する値のことです（正確には，Hartree-Fock極限に収束させるためには，ノードが異なる基底関数も必要です）。また，厳密な値とHartree-Fock極限との差のエネルギー（-0.042 (Hartree)）のことを相関エネルギーと呼び，これは$ (4) $式で，$ \psi\left( \vec{x}_{1},\vec{x}_{2}\right) $を，$ \phi\left( \vec{r}_{1}\right) \phi\left( \vec{r}_{2}\right) $(×スピン部分)と変数分離してしまったことに起因する，誤差のエネルギーです。
ここまで，ヘリウム原子のSchrödinger方程式を，Hartree-Fock法の「考え方」を用いて解く方法を解説してきました。以下では，一般の多電子原子や分子に対するSchrödinger方程式を解くために，この議論を一般化することにします。

Hartree-Fock方程式

N電子系（電子がN個ある系）における，非相関の反対称な波動関数は，以下のような式で表されることが知られています（非相関というのは，多体問題の効果を無視していることを意味しています）。

\psi\left(  \vec{x}_{1,\ldots,}\vec{x}_{N}\right)  =\dfrac{1}{\sqrt{N!}%
}\left\vert
\begin{array}
[c]{cccc}%
\phi_{1}\left(  \vec{x}_{1}\right)   & \phi_{2}\left(  \vec{x}_{1}\right)   &
\cdots & \phi_{N}\left(  \vec{x}_{1}\right)  \\
\phi_{1}\left(  \vec{x}_{2}\right)   & \phi_{2}\left(  \vec{x}_{2}\right)   &
\cdots & \phi_{N}\left(  \vec{x}_{2}\right)  \\
\vdots & \vdots & \ddots & \vdots\\
\phi_{1}\left(  \vec{x}_{N}\right)   & \phi_{2}\left(  \vec{x}_{N}\right)   &
\cdots & \phi_{N}\left(  \vec{x}_{N}\right)
\end{array}
\right\vert \quad(34)

$ (34) $式は，1929年にJ. C. Slaterによって導かれたので（参考サイト[10]），Slater行列式と呼ばれています。ここで，$ \phi_{i}\left( \vec{x}_{j}\right) =\psi_{i}^{\sigma}\left( \vec{r}%
_{j}\right) \alpha_{\sigma}\left( s_{j}\right) $は1粒子「スピン軌道」であり，空間位置座標$ \vec{r}_{j} $の関数$ \psi_{i}^{\sigma}\left( \vec{r}_{j}\right) $とスピン座標$ s_{j} $の関数$ \alpha_{\sigma}\left( s_{j}\right) $の積です（ここで，$ \phi_{i} $の添字$ i $は軌道の状態とスピンの状態 $ \sigma $を表しており，$ j $は電子の番号です）。ここで本来は，RHF（Restricted Hartree-Fock，制限ハートリー・フォック）法と，UHF（Unrestricted Hartree-Fock，非制限ハートリー・フォック）法の細かい説明が必要なのですが，ここでは省略し，以下ではRHF法についてのみ述べます。
$ (34) $式を使うと，（核間反発エネルギーを除いた）エネルギー固有値$ E $は，以下の式で与えられます。

E=%
%TCIMACRO{\dsum \limits_{i,\sigma}}%
%BeginExpansion
{\displaystyle\sum\limits_{i,\sigma}}
%EndExpansion%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\psi_{i}^{\sigma\ast}\left(  \vec{r}\right)  \left[  -\dfrac{1}{2}\nabla^{2}-%
%TCIMACRO{\dsum \limits_{n}}%
%BeginExpansion
{\displaystyle\sum\limits_{n}}
%EndExpansion
\dfrac{Z_{n}}{\left\vert \vec{r}-\vec{R}_{n}\right\vert }\right]  \psi
_{i}^{\sigma}\left(  \vec{r}\right)  + \\
\dfrac{1}{2}%
%TCIMACRO{\dsum \limits_{i,j,\sigma,\sigma^{\prime}}}%
%BeginExpansion
{\displaystyle\sum\limits_{i,j,\sigma,\sigma^{\prime}}}
%EndExpansion%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}d\vec{r}^{\prime}\psi_{i}^{\sigma\ast}\left(  \vec{r}\right)  \psi
_{j}^{\sigma^{\prime}\ast}\left(  \vec{r}^{\prime}\right)  \dfrac
{1}{\left\vert \vec{r}-\vec{r}^{\prime}\right\vert }\psi_{i}^{\sigma}\left(
\vec{r}\right)  \psi_{j}^{\sigma^{\prime}}\left(  \vec{r}^{\prime}\right) \\
-\dfrac{1}{2}%
%TCIMACRO{\dsum \limits_{i,j,\sigma}}%
%BeginExpansion
{\displaystyle\sum\limits_{i,j,\sigma}}
%EndExpansion%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}d\vec{r}^{\prime}\psi_{i}^{\sigma\ast}\left(  \vec{r}\right)  \psi
_{j}^{\sigma\ast}\left(  \vec{r}^{\prime}\right)  \dfrac{1}{\left\vert \vec
{r}-\vec{r}^{\prime}\right\vert }\psi_{j}^{\sigma}\left(  \vec{r}\right)
\psi_{i}^{\sigma}\left(  \vec{r}^{\prime}\right) \quad(35)

ただし，添字$ n $ に対する総和は，系に存在する核の数にわたってとることにします。ここで，

h\left(  \vec{r}\right)  =-\dfrac{1}{2}\nabla^{2}-%
%TCIMACRO{\dsum \limits_{n}}%
%BeginExpansion
{\displaystyle\sum\limits_{n}}
%EndExpansion
\dfrac{Z_{n}}{\left\vert \vec{r}-\vec{R}_{n}\right\vert } \quad(36a) \\
J\left(  \vec{r}\right)  \psi\left(  \vec{r}\right)  =%
%TCIMACRO{\dsum \limits_{l}}%
%BeginExpansion
{\displaystyle\sum\limits_{l}}
%EndExpansion%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}^{\prime}\psi_{l}^{\ast}\left(  \vec{r}^{\prime}\right)  \psi
_{l}\left(  \vec{r}^{\prime}\right)  \dfrac{1}{\left\vert \vec{r}-\vec
{r}^{\prime}\right\vert }\psi\left(  \vec{r}\right) \quad(36b) \\
K\left(  \vec{r}\right)  \psi\left(  \vec{r}\right)  =%
%TCIMACRO{\dsum \limits_{l}}%
%BeginExpansion
{\displaystyle\sum\limits_{l}}
%EndExpansion%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}^{\prime}\psi_{l}^{\ast}\left(  \vec{r}^{\prime}\right)  \psi\left(
\vec{r}^{\prime}\right)  \dfrac{1}{\left\vert \vec{r}-\vec{r}^{\prime
}\right\vert }\psi_{l}\left(  \vec{r}\right) \quad(36c)

とおきます。ただし，添字$ l $に対する総和は，占有軌道にわたってとることとします。また，$ (36b) $式の$ J\left( \vec{r}\right) $はCoulomb演算子，$ (36c) $式の$ K\left( \vec{r}\right) $は交換演算子と呼ばれます。$ h\left( \vec{r}\right) $，Coulomb演算子および交換演算子を用いて$ (35) $式を書き直すと，

E=2%
%TCIMACRO{\dsum \limits_{k}}%
%BeginExpansion
{\displaystyle\sum\limits_{k}}
%EndExpansion%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\psi_{k}^{\ast}\left(  \vec{r}\right)  h\left(  \vec{r}\right)
\psi_{k}\left(  \vec{r}\right)  +%
%TCIMACRO{\dsum \limits_{k}}%
%BeginExpansion
{\displaystyle\sum\limits_{k}}
%EndExpansion
\left(  2%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\psi_{k}^{\ast}\left(  \vec{r}\right)  J\left(  \vec{r}\right)
\psi_{k}\left(  \vec{r}\right)  -%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\psi_{k}^{\ast}\left(  \vec{r}\right)  K\left(  \vec{r}\right)
\psi_{k}\left(  \vec{r}\right)  \right) \quad(37)

と多少見通しが良く（？）なります。ここで，添字$ k $に対する総和は，$ l $と同様に，占有軌道にわたってとります。また$ (37) $式から，以下の方程式が得られます（導出は省略します）。

F\left(  \vec{r}\right)  \psi_{k}\left(  \vec{r}\right)  =\epsilon_{k}%
\psi_{k}\left(  \vec{r}\right)　\quad(38)

$ (38) $式が，Hartree-Fock方程式です。ここで，$ F\left( \vec{r}\right) $は，Fock演算子と呼ばれており，

F\left(  \vec{r}\right)  =h\left(  \vec{r}\right)  +2J\left(  \vec{r}\right)
-K\left(  \vec{r}\right) \quad(39)

です。

Roothaan方程式

すでに述べたように，$ (38) $式のHartree-Fock方程式を完全に数値的に解くことは，一般的には不可能です。従って，有限個の基底を用いて行列方程式に変換して解く必要があります。$ (20) $式および $ (21) $式と同様に，空間軌道関数$ \psi_{k}\left( \vec{r}\right) $は，$ m $個の基底関数$ \chi_{p} $の線形結合で展開して，

\psi_{k}\left(  \vec{r}\right)  =%
%TCIMACRO{\dsum \limits_{p=1}^{m}}%
%BeginExpansion
{\displaystyle\sum\limits_{p=1}^{m}}
%EndExpansion
\chi_{pk}\left(  \vec{r}\right) \quad(40)

と表すことができます。このとき$ (38) $式は，$ (26) $式と同様に，

\mathbf{F}\vec{C}_{k}=\epsilon_{k}\mathbf{S}\vec{C}_{k} \quad(41)

なる行列方程式（一般化固有値問題）になります。$ (41) $式は，Roothaan方程式と呼ばれています。ここで，Fock行列$ \mathbf{F} $の要素$ F_{pq} $は，

F_{pq}=h_{pq}+%
%TCIMACRO{\dsum \limits_{k}}%
%BeginExpansion
{\displaystyle\sum\limits_{k}}
%EndExpansion%
%TCIMACRO{\dsum \limits_{r,s=1}^{m}}%
%BeginExpansion
{\displaystyle\sum\limits_{r,s=1}^{m}}
%EndExpansion
C_{rk}^{\ast}C_{sk}\left(  2\left\langle pr\right\vert g\left\vert
qs\right\rangle -\left\langle pr\right\vert g\left\vert sq\right\rangle
\right) \quad(42)

で与えられます。ただし，

h_{pq}=%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{p}^{\ast}\left(  \vec{r}\right)  \left[  -\dfrac{1}{2}\nabla
^{2}-%
%TCIMACRO{\dsum \limits_{n}}%
%BeginExpansion
{\displaystyle\sum\limits_{n}}
%EndExpansion
\dfrac{Z_{n}}{\left\vert \vec{r}-\vec{R}_{n}\right\vert }\right]  \chi
_{q}\left(  \vec{r}\right) \quad(43a) \\
\left\langle pr\right\vert g\left\vert qs\right\rangle =%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}_{1}d\vec{r}_{2}\chi_{p}^{\ast}\left(  \vec{r}_{1}\right)  \chi
_{r}^{\ast}\left(  \vec{r}_{2}\right)  \dfrac{1}{\left\vert \vec{r}_{1}%
-\vec{r}_{2}\right\vert }\chi_{q}\left(  \vec{r}_{1}\right)  \chi_{s}\left(
\vec{r}_{2}\right) \quad(43b) \\
\left\langle pr\right\vert g\left\vert sq\right\rangle =%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}_{1}d\vec{r}_{2}\chi_{p}^{\ast}\left(  \vec{r}_{1}\right)  \chi
_{r}^{\ast}\left(  \vec{r}_{2}\right)  \dfrac{1}{\left\vert \vec{r}_{1}%
-\vec{r}_{2}\right\vert }\chi_{s}\left(  \vec{r}_{1}\right)  \chi_{q}\left(
\vec{r}_{2}\right) \quad(43c)

です。ただし，添字$ k $は，軌道$ \psi_{k} $を表し，また添字$ p $, $ q $, $ r $, $ s $は基底を表します。ここで，$ (43a) $を1電子積分，$ (43b) $をCoulomb積分， $ (43c) $を交換積分と呼びます。また，$ (43b) $と$ (43c) $をあわせて，2電子積分と呼びます。また，$ \mathbf{S} $は重なり行列であり，その要素$ S_{pq} $は，

S_{pq}=%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{p}^{\ast}\left(  \vec{r}\right)  \chi_{q}\left(  \vec{r}\right) \quad(44)

で定義されます。系の（核間反発エネルギーを除いた）全エネルギー$ E $については，密度行列$ \mathbf{P} $と呼ばれるものを導入すると便利であり，その要素$ P_{pq} $は，

P_{pq}\mathbf{=2}%
%TCIMACRO{\dsum \limits_{k}}%
%BeginExpansion
{\displaystyle\sum\limits_{k}}
%EndExpansion
C_{pk}C_{qk}^{\ast} \quad(45)

で定義されます。ここで，$ (45) $式を使うと，フォック行列$ \mathbf{F} $の要素$ F_{pq} $は，

F_{pq}=h_{pq}+\dfrac{1}{2}%
%TCIMACRO{\dsum \limits_{r,s=1}^{m}}%
%BeginExpansion
{\displaystyle\sum\limits_{r,s=1}^{m}}
%EndExpansion
P_{rs}\left(  2\left\langle pr\right\vert g\left\vert qs\right\rangle
-\left\langle pr\right\vert g\left\vert sq\right\rangle \right) \quad(46)

と書けます。従って，系の全エネルギー$ E $は，

E=%
%TCIMACRO{\dsum \limits_{p,q=1}^{m}}%
%BeginExpansion
{\displaystyle\sum\limits_{p,q=1}^{m}}
%EndExpansion
P_{pq}h_{pq}+\dfrac{1}{2}%
%TCIMACRO{\dsum \limits_{p,q,r,s=1}^{m}}%
%BeginExpansion
{\displaystyle\sum\limits_{p,q,r,s=1}^{m}}
%EndExpansion
P_{pq}P_{rs}\left(  \left\langle pr\right\vert g\left\vert qs\right\rangle
-\dfrac{1}{2}\left\langle pr\right\vert g\left\vert sq\right\rangle \right) \quad (47)

となります。

ヘリウム原子におけるSlater行列式，Hartree-Fock方程式およびRoothaan方程式

最後に，ヘリウム原子において，Slater行列式，Hartree-Fock方程式およびRoothaan方程式を考え，これらがそれぞれ $ (3) $式，$ (7) $式および$ (26) $式と等しくなることを確認してみましょう。

ヘリウム原子におけるSlater行列式とHartree-Fock方程式

まず，ヘリウム原子（2原子系）に対するSlater行列式は，

\psi\left(  \vec{x}_{1},\vec{x}_{2}\right)  =\dfrac{1}{\sqrt{2!}}\left\vert
\begin{array}
[c]{cc}%
\phi_{1}\left(  \vec{x}_{1}\right)  & \phi_{2}\left(  \vec{x}_{1}\right) \\
\phi_{1}\left(  \vec{x}_{2}\right)  & \phi_{2}\left(  \vec{x}_{2}\right)
\end{array}
\right\vert =\dfrac{1}{\sqrt{2}}\left[  \phi_{1}\left(  \vec{x}_{1}\right)
\phi_{2}\left(  \vec{x}_{2}\right)  -\phi_{1}\left(  \vec{x}_{2}\right)
\phi_{2}\left(  \vec{x}_{1}\right)  \right] \quad(48)

となり，$ (48) $式は，確かに$ (3) $式と全く同じ式となっています。$ (36a) $式については，ヘリウム原子の場合，核は1個しか存在せず($ n=1 $)，かつ核の電荷は$ 2e $($ Z_{1}=2 $)なので，

h\left(  \vec{r}\right)  =-\dfrac{1}{2}\nabla^{2}-\dfrac{2}{\left\vert
\vec{r}\right\vert }=-\dfrac{1}{2}\nabla^{2}-\dfrac{2}{r} \quad(49)

となります。ここで，$ \vec{R}_{1} $を原点にとりました($ \vec{R}_{1}=\vec{0} $)。また，Coulomb演算子$ J\left( \vec{r}\right) $は，占有軌道が1個なので，総和をとる必要がなく，従って添字$ l $は不要であることに注意すると，

J\left(  \vec{r}\right)  \psi\left(  \vec{r}\right)  =%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}^{\prime}\psi^{\ast}\left(  \vec{r}^{\prime}\right)  \psi\left(
\vec{r}^{\prime}\right)  \dfrac{1}{\left\vert \vec{r}-\vec{r}^{\prime
}\right\vert }\psi\left(  \vec{r}\right) \quad(50)

です。よって，

J\left(  \vec{r}\right)  =%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}^{\prime}\psi^{\ast}\left(  \vec{r}^{\prime}\right)  \psi\left(
\vec{r}^{\prime}\right)  \dfrac{1}{\left\vert \vec{r}-\vec{r}^{\prime
}\right\vert }=%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}^{\prime}\left\vert \psi\left(  \vec{r}^{\prime}\right)  \right\vert
^{2}\dfrac{1}{\left\vert \vec{r}-\vec{r}^{\prime}\right\vert } \quad(51)

となります。同様に，$ l $が不要であることに注意すると，交換演算子$ K\left( \vec{r}\right) $は，

K\left(  \vec{r}\right)  \psi\left(  \vec{r}\right)  =%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}^{\prime}\psi^{\ast}\left(  \vec{r}^{\prime}\right)  \psi\left(
\vec{r}^{\prime}\right)  \dfrac{1}{\left\vert \vec{r}-\vec{r}^{\prime
}\right\vert }\psi\left(  \vec{r}\right) \quad(52)

より，

K\left(  \vec{r}\right)  =%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}^{\prime}\psi^{\ast}\left(  \vec{r}^{\prime}\right)  \psi\left(
\vec{r}^{\prime}\right)  \dfrac{1}{\left\vert \vec{r}-\vec{r}^{\prime
}\right\vert }=%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}^{\prime}\left\vert \psi\left(  \vec{r}^{\prime}\right)  \right\vert
^{2}\dfrac{1}{\left\vert \vec{r}-\vec{r}^{\prime}\right\vert } \quad(53)

となります。よって，Fock演算子$ F\left( \vec{r}\right) $は，

F\left(  \vec{r}\right)  =-\dfrac{1}{2}\nabla^{2}-\dfrac{2}{r}+2%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}^{\prime}\left\vert \psi\left(  \vec{r}^{\prime}\right)  \right\vert
^{2}\dfrac{1}{\left\vert \vec{r}-\vec{r}^{\prime}\right\vert }-%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}^{\prime}\left\vert \psi\left(  \vec{r}^{\prime}\right)  \right\vert
^{2}\dfrac{1}{\left\vert \vec{r}-\vec{r}^{\prime}\right\vert } \\
=-\dfrac{1}{2}\nabla^{2}-\dfrac{2}{r}+%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}^{\prime}\left\vert \psi\left(  \vec{r}^{\prime}\right)  \right\vert
^{2}\dfrac{1}{\left\vert \vec{r}-\vec{r}^{\prime}\right\vert } \quad(54)

となります。ここで，添字$ k $は不要であることを用いると，$ (38) $式のHartree-Fock方程式は結局，

\left[  -\dfrac{1}{2}\nabla^{2}-\dfrac{2}{r}+%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}^{\prime}\left\vert \psi\left(  \vec{r}^{\prime}\right)  \right\vert
^{2}\dfrac{1}{\left\vert \vec{r}-\vec{r}^{\prime}\right\vert }\right]
\psi\left(  \vec{r}\right)  =\epsilon\psi\left(  \vec{r}\right) \quad(55)

となりますが，$ (55) $式は，確かに$ (7) $式と全く同じ式となっています（$ \vec{r} $を$ \vec{r}_{1} $に，$ \vec{r}^{\prime} $を$ \vec{r}_{2} $に，また$ \epsilon $を$ E^{\prime} $に適宜読み替えてください）。つまり$ (7) $式は，$ (38) $式の特別な場合であることが，これで確認できたことになります。

ヘリウム原子におけるRoothaan方程式

次に，ヘリウム原子におけるRoothaan方程式を考えてみましょう。ここで，基底関数$ \chi_{p}\left( \vec{r}\right) $は，実関数に限定するものとします。1電子積分$ h_{pq} $は$ (49) $式から容易に導くことができ，

h_{pq}=%
%TCIMACRO{\dint }%
%BeginExpansion
{\displaystyle\int}
%EndExpansion
d\vec{r}\chi_{p}\left(  \vec{r}\right)  \left(  -\dfrac{1}{2}\nabla^{2}%
-\dfrac{2}{r}\right)  \chi_{q}\left(  \vec{r}\right) \quad(56)

となりますが，$ (56) $式は，$ (24a) $式そのものです（$ \vec{r} $を$ \vec{r}_{1} $に適宜読み替えてください）。次に2電子積分ですが，ヘリウム原子の場合，Coulomb積分$ \left\langle pr\right\vert g\left\vert qs\right\rangle $と，交換積分$ \left\langle pr\right\vert g\left\vert sq\right\rangle $は等しくなります（これは，$ (29) $式で，添字$ q $と$ s $を入れ替えても，結果が変わらないことから理解できます）。従って，ヘリウム原子におけるFock行列の要素$ F_{pq} $は，

F_{pq}=h_{pq}+%
%TCIMACRO{\dsum \limits_{r,s=1}^{m}}%
%BeginExpansion
{\displaystyle\sum\limits_{r,s=1}^{m}}
%EndExpansion
C_{r}C_{s}\left\langle pr\right\vert g\left\vert qs\right\rangle \quad(57)

となります（ただし，添字$ k $は不要であることと，$ C_{r}^{\ast} $が実数になることを用いました）。$ (57) $式は，$ (25) $式に他なりません（$ \left\langle pr\right\vert g\left\vert qs\right\rangle $を$ Q_{prqs} $に適宜読み替えてください）。さらに，基底関数$ \chi_{p}\left( \vec{r}\right) $を，実関数に限定したので，$ (44) $式は，$ (24c) $式と全く等しくなります。よって，$ (26) $式が，$ (41) $式の特別な場合であることが，これで確認できました。
さらに，添字$ k $が不要であることに注意すると，密度行列$ \mathbf{P} $の要素$ P_{pq} $は，

P_{pq}=2C_{p}C_{q} \quad(58)

となります。ここで，$ (58) $式を$ (47) $式に代入すると，系の全エネルギー$ E $は，

E=2%
%TCIMACRO{\dsum \limits_{p,q=1}^{m}}%
%BeginExpansion
{\displaystyle\sum\limits_{p,q=1}^{m}}
%EndExpansion
C_{p}C_{q}h_{pq}+%
%TCIMACRO{\dsum \limits_{p,q,r,s=1}^{m}}%
%BeginExpansion
{\displaystyle\sum\limits_{p,q,r,s=1}^{m}}
%EndExpansion
C_{p}C_{q}C_{r}C_{s}\left\langle pr\right\vert g\left\vert qs\right\rangle \quad(59)

となりますが，$ (59) $式は$ (32) $式そのものです。

まとめ

1. Hartree-Fock法の考え方を，ヘリウム原子を例にとって説明しました。また，ヘリウム原子に対する試行関数の，基底による展開を例にとって，Roothaan方程式の考え方を説明しました。
2. 一般的なHartree-Fock法について，RHF法に限り説明しました。Slater行列式から，Hartree-Fock方程式が導かれることを示しました（実際の導出はしていませんが…）。
3. Hartree-Fock方程式の軌道$ \psi_{k}\left( \vec{r}\right) $を基底で展開することにより，Roothaan方程式が導かれることを説明しました。
4. ヘリウム原子におけるSlater行列式，Hartree-Fock方程式およびRoothaan方程式を導出し，これらが1.で導いた方程式と一致することを確認しました。

この記事の続きをできるだけ簡単にHartree-Fock法を解説してみる part2と題して書きました。よろしければご覧ください。

参考サイト

[1] 水素様原子のエネルギースペクトル解法（その１）〜 Schrödingerによるラゲール陪多項式を用いた波動方程式解法 〜
[2] 多体問題 (量子論) - Wikipedia
[3] ハートリー＝フォック方程式 - Wikipedia
[4] 23. Born−Oppenheimer近似と断熱近似 - 「yam＠広島大」物理化学Monographシリーズ
[5] 量子力学ホームページ (2007 年度前期) ノート(§4-§6)
[6] 原子単位系 - Wikipedia
[7] ハートリー近似 - Wikipedia
[8] Numerical Methods in Quantum Mechanics
[9] HFCXX - Hartree-Fock C++ code
[10] Slater determinant - Wikipedia

参考文献

[1] A.ザボ, N.S.オストランド『新しい量子化学―電子構造の理論入門〈上〉』東京大学出版会（1987）
[2] J.M.ティッセン 『計算物理学』シュプリンガー・フェアラーク東京（2003）
[3] J. M. Perez-Jorda, Phys. Rev. E. 90, 053307 (2014).
[4] J. T. Zung, J. Chem. Phys. 41, 2888 (1964).
[5] K. Szalewicz and H. J. Monkhorst, J. Chem. Phys. 75, 5785 (1981).
[6] G. Tanner, et al., Rev. Mod. Phys. 72, 497 (2000).
[7] R.M.マーチン 『物質の電子状態 上』シュプリンガー・ジャパン株式会社（2010）