高校までしか数学をやっていない人でも分かるハートリー・フォック法

投稿日: 投稿者:

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

計算化学プログラムで分子軌道を求める計算は何をやっているか知っておきたい。分子軌道と電子エネルギーを求める単純ヒュッケル法の次の段階として（DFTの前に）ハートリー・フォック法を理解するため、できるだけ全部の数字を書きながら計算を進めたいと思います。教科書は、

です。

以下、ウェブで閲覧できるハートリー・フォック法の解説、

ヘリウム原子の計算コード

計算の準備

高校数学の復習

指数関数の微分

※ 微分は \((f(x))’\) ではなく \(\frac{d f(x)}{dx}\) や\(\frac{d}{dx} f(x)\)と書きます。 \((f(x))'{}’\) は \(\frac{d^2 f(x)}{dx^2}\) と書きます。\[\frac{d}{dx}e^{kx} = k e^{kx} \]

指数関数の積分

\[\int_a^b e^{kx}dx = \left[\frac{1}{k} e^{kx} \right]_a^b = \frac{1}{k}\left(e^{kb} – e^{ka}\right)\]

関数の積の微分

\[ \frac{d}{dx} (f(x)g(x)) = \frac{d f(x)}{dx}g(x) + f(x)\frac{d g(x)}{dx} \]

部分積分

\[ \int \left(\frac{d}{dx}f(x)\right)g(x)\, dx = f(x)g(x) – \int f(x)\left(\frac{d}{dx}g(x)\right) dx \]

行列の足し算

\[ \begin{pmatrix} a & b \\ c & d \end{pmatrix} + \begin{pmatrix} p & q \\ r & s \end{pmatrix} = \begin{pmatrix} a+p & b+q \\ c+r & d+s \end{pmatrix} \]

行列の掛け算

\[ \begin{pmatrix} a & b \\ c & d \end{pmatrix} \times \begin{pmatrix} p & q \\ r & s \end{pmatrix} = \begin{pmatrix} ap+br & aq+bs \\ cp+dr & cq+ds \end{pmatrix} \]

ハートリー原子単位系

電子の位置を原子核からの距離\(r\)、極角\(\theta\)、方位角\(\phi\)で表わした時、水素原子の1s原子軌道は\[ \psi_{1s}(r,\theta,\phi) = \frac{1}{\sqrt{\pi a_{0}^3}}\exp\left(-\frac{r}{a_0}\right) \]です。球対称なので極角（\(\theta\)）と方位角（\(\phi\)）を含まない距離\(r\)だけの関数です。ここで、 \(a_0 = \frac{4\pi\varepsilon_0\hbar^2}{m_e e^2}\) はボーア半径（約0.529 Å）、\(m_e\)は電子の質量、\(\hbar \)はプランク定数\(h\)を\(2\pi\)で割った定数、\(e\)は電気素量、\(\varepsilon_0\)は真空の誘電率です。

原子単位系では、電荷の基本単位を\(e\)、質量の基本単位を\(m_e\)、作用の基本単位を\(\hbar\)、長さの基本単位を\(a_0\)、エネルギーの基本単位を\(E_h\)（\(\approx 4.36\times 10^{-18} \mathrm{J}\)）とします。誘電率の組立単位は\(\frac{e^2}{a_0 E_h} = 4\pi \varepsilon_0\) となります。したがって、原子単位系では1s原子軌道は\[ \psi_{1s}(r) = \frac{1}{\sqrt{\pi}}\exp(-r) \]という単純な式になりますし、水素原子核と電子の静電ポテンシャルエネルギー\(V(r) = -\frac{1}{4\pi \varepsilon_0}\frac{e^2}{r}\)も\(V(r) = -\frac{1}{r} \) と書けます。以下では距離もボーア半径 \(a_0\) を単位として表すため、例えば \(r=1\) は「原子核から \(1a_0\) 離れた位置」を意味します。

\(\exp(-r)\)は指数関数で、\(e^{-r}\)（ネイピア数\(e \approx 2.72\) の\(-r\)乗）の別の書き方です（ネイピア数と電気素量はどちらも\(e\)で表わすので紛らわしいのですが、原子単位では以後電気素量\(e\)は出てきません）。

エネルギー演算子

水素原子の1s軌道は、有名なシュレーディンガー方程式\[ \hat{H}\psi = E\psi \]を解いて得られました。方程式の意味は、関数\(\psi\) に\(\hat{H}\)という演算（掛け算や微分など）をすると、元の関数に係数\(E\)が掛かるような関数\(\psi\)を求めなさい、というものです。\(f(x) = e^{ax}\)という関数に微分という演算をすると\(f'(x)=a e^{ax}\)となって、元の関数に係数\(a\)が掛かりますね、みたいなことです。\(E\)はエネルギーであり、電子のエネルギーは、運動エネルギー（正の値）とポテンシャルエネルギー（負の値）からなります。\(\hat{H}\)は粒子の全エネルギー（運動エネルギーとポテンシャルエネルギーの和）を現わす演算子で、ハミルトニアンと呼ばれます。

\( \psi_{1s}(r) = \frac{1}{\sqrt{\pi}}e^{-r} \) に、\(\hat{H}\)という演算をして係数が掛かるか（\(\frac{1}{\sqrt{\pi}}e^{-r}\)がハミルトニアンの固有関数になっているかどうか）確かめてみます。

運動エネルギーを求めるには、\(x,y,z\)軸がある直交座標では関数を\(x,y,z\)でそれぞれ2回微分（偏微分）して足し合わせる計算をします。高校物理ではド・ブロイ波（物質波）を学びますが、粒子の運動量を\(p\)と置いたの一次元ド・ブロイ波式は\( f(x) = e^{i\frac{p}{\hbar}x}\)です（\(\hbar\)はプランク定数\(h\)を\(2\pi\)で割った定数）。

高校物理では\(運動量 p = 質量 × 速さ = mv\)、\(運動エネルギー T = 1/2 × 質量 × 速さの二乗 = \frac{1}{2}mv^2\)と習います。この2つの式から\(T = \frac{p^2}{2m}\) が導かれます。ド・ブロイ波の式を1回微分すると元の波に\(i\frac{p}{\hbar} \)が掛かり、2回微分すると\(-\frac{p^2}{\hbar^2}\)が掛かります。したがって、\(-\frac{\hbar^2}{2m}\frac{d^2}{dx^2}\)を作用させると、高校物理で習う古典力学の運動エネルギー\(\frac{p^2}{2m}\) が係数として現われます。原子単位系を使うので、以下の式には\(\hbar\)と\(m_e\)が含まれていません。

水素原子の1s原子軌道を\(x,y,z\)ではなく原子核からの距離\(r\)だけの関数で表わしているので、少し違う計算（演算）をします。演算子は

\[ \hat{T} = -\frac{1}{2}\left[\frac{1}{r^2}\frac{d}{dr}\left(r^2 \frac{d}{dr}\right)\right] \]

です。

\[ \hat{T}\psi_{1s}(r) = -\frac{1}{2}\left[\frac{1}{r^2}\frac{d}{dr}\left(r^2 \frac{d}{dr}\right)\right] \psi_{1s}(r) \]

という式で行うことは、\(\psi_{1s}(r)\)を微分して、\(r^2\)を掛けて、また微分した後に、\(\frac{1}{r^2}\)を掛けて、最後に\(-\frac{1}{2}\)を掛ける、です。実際に計算してみます。途中で、関数の積の微分 \((r^2 \times e^{-r})’ = (r^2)’ \times e^{-r} + r^2 \times (e^{-r})’ = (2r – r^2) e^{-r} \) を使います。

\[ -\frac{1}{2}\left[\frac{1}{r^2}\frac{d}{dr}\left(r^2 \frac{d}{dr}\right)\right] \frac{1}{\sqrt{\pi}}e^{-r} \]

\[ =  -\frac{1}{2}\left[\frac{1}{r^2}\frac{d}{dr}\right] r^2\left(-1\right)  \frac{1}{\sqrt{\pi}}e^{-r}\]

\[ = -\frac{1}{2}\left[\frac{1}{r^2}\right] (-1)(2r – r^2)   \frac{1}{\sqrt{\pi}}e^{-r} \]

\[ = \left(\frac{1}{r} – \frac{1}{2}\right)   \frac{1}{\sqrt{\pi}}e^{-r} \]

\(\frac{1}{r}\)という項が残っていますが、この項はポテンシャルエネルギー演算子によって打ち消されます。電子の全エネルギーを表わす演算子（ハミルトニアン）は

\[\hat{H} = \hat{T} + \hat{V}\]

で、水素原子の場合は\(\hat{V}= -\frac{1}{r}\)（原子核と電子のクーロン相互作用）なので、

\[ \hat{H}\psi_{1s}(r) = (\hat{T} + \hat{V})\psi_{1s}(r) = \left[\left(\frac{1}{r} – \frac{1}{2}\right) + \left(-\frac{1}{r}\right)\right] \psi_{1s}(r) = -\frac{1}{2} \psi_{1s}(r) \]

で、水素原子の基底状態の全エネルギー\(E =-0.5\;E_{\mathrm h}\) が係数として出てくることが確かめられました。

エネルギーの期待値

運動エネルギーの期待値\(\langle{}T\rangle\)とポテンシャルエネルギーの期待値\(\langle{}V\rangle\)を求めるには、\(\hat{H}\psi = E\psi\)の両辺に左から\(\psi^*\)（\(\psi\)の複素共役。例えば\(f(x) = \cos x + i\sin x\)の複素共役は\(f^*(x) = \cos x – i\sin x\)。水素原子の1s軌道は実関数なのでここでは\(\psi^{*} = \psi\)）を掛けて、全空間に渡って積分します。

任意の波動関数に対するエネルギー期待値は、
\[\langle{}H\rangle_\psi = \frac{\int\psi^* \hat{H}\psi d\tau}{\int|\psi|^2 d\tau} \]と定義されます。\(\psi\) が規格化されていれば分母は1です。さらに \(\psi\) がハミルトニアンの固有関数なら、\(\langle H\rangle=E\) になります。

分母は1として、積分を全部書くと、\[ \int_0^{2\pi}\!\!\!\int_0^{\pi}\!\!\!\int_0^{\infty} \psi^* \hat{H} \psi\ r^2 \sin\theta dr d\theta d\phi = E \int_0^{2\pi}\!\!\!\int_0^{\pi}\!\!\!\int_0^{\infty} \psi^* \psi\ r^2 \sin\theta dr d\theta d\phi. \]

\(E\)は単なる数値なので、積分記号の外に出せます。1s原子軌道は規格化されている（絶対値の二乗を積分すると1になるように係数が調整されている）ので、\(\int_0^{2\pi}\!\!\!\int_0^{\pi}\!\!\!\int_0^{\infty} \psi^* \psi\ r^2 \sin\theta dr d\theta d\phi =1\)です。したがって、

\[ E = \int_0^{2\pi}\!\!\!\int_0^{\pi}\!\!\!\int_0^{\infty} \psi^* \hat{H} \psi\ r^2 \sin\theta dr d\theta d\phi \]

\[ = \int_0^{2\pi}\!\!\!\int_0^{\pi}\!\!\!\int_0^{\infty} \psi^* \hat{T} \psi r^2 \sin\theta dr d\theta d\phi + \int_0^{2\pi}\!\!\!\int_0^{\pi}\!\!\!\int_0^{\infty} \psi^* \hat{V} \psi\cdot{} r^2 \sin\theta dr d\theta d\phi \]

運動エネルギーの期待値は、

\[ \langle{}T\rangle = \int_0^{2\pi}\!\!\!\int_0^{\pi}\!\!\!\int_0^{\infty} \frac{1}{\sqrt{\pi}}e^{-r} \left(\frac{1}{r} – \frac{1}{2}\right) \frac{1}{\sqrt{\pi}}e^{-r} \cdot{} r^2 \sin\theta dr d\theta d\phi \]

\[ = \frac{1}{\pi} \int_0^{2\pi}\!\!\!\int_0^{\pi}\!\!\!\int_0^{\infty} \left(r – \frac{r^2}{2}\right) e^{-2r} \sin\theta dr d\theta d\phi .\]

\(\theta\)と\(\phi\)の積分を先に済ませます。\(\int_0^{\pi} \sin\theta d\theta = [-\cos\theta]_0^{\pi} = -(-1 – 1) = 2 \)、\(\int_0^{2\pi} d\phi = [\phi]_0^{2\pi} = 2\pi \)。
\[ \langle{}T\rangle =4\int_0^{\infty} \left(r – \frac{r^2}{2}\right) e^{-2r} dr \]
部分積分の公式 \(\int f(r) g'(r) dr = f(r)g(r) – \int f'(r)g(r) dr \)を使います。

\[ \int_0^{\infty} r e^{-2r} dr = \left[r\left(–\frac{1}{2}\right)e^{-2r}\right]_0^{\infty} – \int_0^{\infty} \left(–\frac{1}{2}\right)e^{-2r} dr \]
\[= 0 – \left[\frac{1}{4} e^{-2r}\right]_0^{\infty} = \frac{1}{4} .\]
この積分には公式\(\int_0^{\infty} r^n e^{-ar} dr = \frac{n!}{a^{n+1}}\) があって、これを使うと、\[ \int_0^{\infty} r^2 e^{-2r} dr = \frac{2!}{2^3} = \frac{1}{4}.\]

したがって、\[ \langle{}T\rangle =4\int_0^{\infty} \left(r – \frac{r^2}{2}\right) e^{-2r} dr = 4\times \left(\frac{1}{4} – \frac{1}{2}\times \frac{1}{4}\right) = \frac{1}{2}\,E_{\mathrm h}.\]

次に、核引力エネルギーの期待値は、

\[ \langle{}V\rangle = 4\int_0^{\infty} e^{-r} \left(-\frac{1}{r}\right) e^{-r} \times r^2 dr = -4\int_0^{\infty} r e^{-2r} dr = -4 \times \frac{1}{4} = -1\,E_{\mathrm h}.\]

\(V\propto-1/r\) 型ポテンシャルに対する定常束縛状態では、ビリアル定理から、\[ \langle{}V\rangle = -2\langle{}T\rangle \] が成り立ちます。同じ \(1/r\) 型のポテンシャルをもつ重力系でもビリアル定理が成り立つので、重力的に十分平衡化された（ビリアル化した）銀河や銀河団の質量の推定にも利用できるそうです。ビリアル定理から推定した銀河団の総質量が光から見積られる物質の質量よりはるかに大きいことは、宇宙の暗黒物質（ダークマター）の存在を示す証拠の一つとされているそうです。

ヘリウム原子の計算

ヘリウム原子の波動関数

ヘリウム原子は電子が2つあります。厳密な空間波動関数は、電子1の原子核からの距離\(r_1\)、電子2の原子核からの距離\(r_2\)だけでなく、電子間距離\(r_{12}\)にも依存する\[ \psi(r_1,r_2,r_{12}) \]のような関数になります。しかし、水素原子のような閉じた形の厳密な解析解は知られていません。そこで、多電子波動関数を一電子軌道から作った1個のスレーター（Slater）行列式で近似します。それが何かは後から説明します。ヘリウム原子の基底状態では、2つの電子が同じ空間軌道を占有するため、その空間部分はまず\[ \Psi(\mathbf{r_1},\mathbf{r_2}) \approx \phi(\mathbf r_1)\phi(\mathbf r_2)\]のような関数の積で表わすことにします。

波動関数の反対称化

電子の位置に依存する空間軌道を\(\phi_{1s}(\mathbf r)\)と書きます。ここで、(\(\mathbf r\))は電子の位置ベクトルです。 上向きと下向きのスピン状態を表す関数を、それぞれ \(\alpha(\sigma)\)、\(\beta(\sigma)\) とします。空間軌道とスピン関数の積をスピン軌道と呼びます。位置とスピンの変数をまとめて \(x=(\mathbf r,\sigma)\) と書くと、今回用いる2つのスピン軌道は\[ \chi_{1s\alpha}(x)=\phi_{1s}(\mathbf r)\alpha(\sigma),
\qquad
\chi_{1s\beta}(x)=\phi_{1s}(\mathbf r)\beta(\sigma) \]です。この2つは空間部分が同じでも、スピン部分が異なるため、異なる一電子状態を表します。

まず、電子1に \(\alpha\)、電子2に \(\beta\) のスピン軌道を割り当てた積を考えます。

\[ \phi_{1s}(\mathbf r_1)\phi_{1s}(\mathbf r_2)\alpha(\sigma_1)\beta(\sigma_2). \]

しかし、電子は互いに区別できない同種のフェルミ粒子です。電子を表すラベル1と2を交換したときに、電子の全波動関数は符号が反転しなければいけません。そこで、ラベルを交換した（\(\phi_{1s}(\mathbf r_2)\phi_{1s}(\mathbf r_1)
\alpha(\sigma_2)\beta(\sigma_1)\)）という項も含めて考えます。ここで「ラベルを交換する」とは、電子を表すラベル1と2を、位置とスピンの両方で入れ換えることです。選択肢としては割り当てを入れ換えた項を元の式から足すか引くかです。2つの積を足した波動関数は、\[ \Psi_{+}(x_1,x_2) = \frac{1}{\sqrt{2}} \left\{\phi_{1s}(\mathbf r_1)\phi_{1s}(\mathbf r_2) \alpha(\sigma_1)\beta(\sigma_2)+\phi_{1s}(\mathbf r_2)\phi_{1s}(\mathbf r_1) \alpha(\sigma_2)\beta(\sigma_1)\right\} \]となります。ここで、\(x_i=(\mathbf r_i,\sigma_i)\) は、電子 \(i\) の位置とスピンの変数をまとめたものです。この波動関数は、電子1と電子2のラベルを入れ換えても、\[\Psi_{+}(x_2,x_1)=\Psi_{+}(x_1,x_2) \]となり、符号が変わりません。したがって、電子の波動関数に必要な反対称性を満たしません。

一方、2つの積を引いた波動関数は、\[\Psi_{-}(x_1,x_2) = \frac{1}{\sqrt{2}} \left\{\phi_{1s}(\mathbf r_1)\phi_{1s}(\mathbf r_2)\alpha(\sigma_1)\beta(\sigma_2) -\phi_{1s}(\mathbf r_2)\phi_{1s}(\mathbf r_1)\alpha(\sigma_2)\beta(\sigma_1)\right\} \]となります。この場合は、\[ \Psi_{-}(x_2,x_1)=-\Psi_{-}(x_1,x_2) \]となるため、必要な反対称性を満たしています。

この式に現われる（\(A\times B – C\times D\)）という形は\(2\times 2\)の行列式\[ \begin{vmatrix} A & C \\ D &B\end{vmatrix} \] と同じ形です。したがって上の波動関数は\[ \Psi_{-}(x_1,x_2) = \frac{1}{\sqrt{2}} \begin{vmatrix} \phi_{1s}(\mathbf r_1)\alpha(\sigma_1) & \phi_{1s}(\mathbf r_1)\beta(\sigma_1) \\ \phi_{1s}(\mathbf r_2)\alpha(\sigma_2) & \phi_{1s}(\mathbf r_2)\beta(\sigma_2) \end{vmatrix} \]と書けます。この行列式が「スレーター行列式」と呼ばれます。ヘリウム原子のような2電子系では、わざわざ行列式で書くとかえって式が大きく見えますが、電子数が増えると、反対称な波動関数を簡潔に表せるので非常に便利です。

パウリの排他原理との関係を確認するために、2つの電子を同じスピン軌道 \(\chi_a\) に入れることを考えます。この場合、反対称化した式は、\[ \frac{1}{\sqrt{2}} \left\{ \chi_a(x_1)\chi_a(x_2) – \chi_a(x_2)\chi_a(x_1) \right\} =0 \]となり、波動関数が消えてしまいます。つまり、2つの電子が同じスピン軌道を占める状態は作れません。これがパウリの排他原理に対応します。各スピン軌道が規格化され、互いに直交しているとき、係数 \(1/\sqrt{2}\) によって全波動関数も規格化されます。スピン軌道の定義を代入すると、共通の空間部分をくくり出せます。

α(σ1)β(σ2)−β(σ1)α(σ2)2⏟反対称なスピン部分\underbrace{ \frac{ \alpha(\sigma_1)\beta(\sigma_2) – \beta(\sigma_1)\alpha(\sigma_2) }{\sqrt{2}} }_{\text{反対称なスピン部分}}

このスピン部分は、全スピンが0のスピン一重項を表します。空間部分は電子の交換に対して対称ですが、スピン部分が反対称なので、全体は反対称になります。

ヘリウム原子のエネルギー演算子

ヘリウム原子は+2の電荷を持つ1つの原子核と2つの電子から構成されます。ヘリウム原子に対するエネルギー演算子、ハミルトニアンは、\[ \hat{H} = \left[-\frac{1}{2}\nabla_1^2 – \frac{2}{r_1}\right] + \left[-\frac{1}{2}\nabla_2^2 – \frac{2}{r_2}\right] + \frac{1}{r_{12}} \]となります。\({r_1}\)は電子1の原子核からの距離、\({r_2}\)は電子2の原子核からの距離、\({r_{12} = |\mathbf{r_1} – \mathbf{r_2} |}\)は電子1と電子2の距離です。\(-\frac{2}{r_1}\) の分子が2になっているのは、ヘリウム原子核の電荷が+2で核引力エネルギーが2倍になるからです。

このハミルトニアンは波動関数のスピン部分には作用しないので、エネルギーを計算するときには、規格化されたスピン部分を積分すると1になり、空間部分\(\phi_{1s}(r_1)\phi_{1s}(r_2)\) だけを使って計算することもできます。しかし、後で交換項がどこから現われるかを理解するために、ここではスピン部分を含めたスレーター行列式のまま計算を進めます。位置とスピンを合わせた変数を\[x=(\mathbf r,\sigma)\]とすると、ヘリウム原子の基底状態の波動関数は \[ \psi(x_1,x_2)= \frac{1}{\sqrt{2}}\begin{vmatrix} \chi_{1s\alpha}(x_1) & \chi_{1s\beta}(x_1)\\ \chi_{1s\alpha}(x_2) & \chi_{1s\beta}(x_2)\end{vmatrix} = \frac{1}{\sqrt{2}}(\chi_{1s\alpha}(x_1)\chi_{1s\beta}(x_2) – \chi_{1s\beta}(x_1)\chi_{1s\alpha}(x_2))\]のです。

また、スピン軌道について積分するときには、空間座標とスピン変数の両方について積分するので、\[
dx=d\tau\,d\sigma\]と書くことにします。

Slater（スレーター）型軌道（STO）を使ってエネルギー評価する

Heの1s空間軌道を\[ \phi_{1s}(r) = \sqrt{\frac{\zeta^3}{\pi}}e^{-\zeta r} \]とします。\(\zeta\)は後でエネルギーが最小になるように調整するためのパラメータです。

運動エネルギーの期待値

まず、電子1の運動エネルギーの期待値\[\langle{}T\rangle_1 = \iint\Psi^*(x_1,x_2)\left(-\frac{1}{2}\nabla_1^2\right)\Psi(x_1,x_2)\,dx_1dx_2\]を求めます。今回用いる軌道やスピン関数は実関数として扱えるので、以下では複素共役の記号 \(^*\) を省略します。

スレーター行列式を代入すると、\[ \langle{}T\rangle_1 = \iint \left(\frac{1}{\sqrt{2}}(\chi_{1s\alpha}(x_1)\chi_{1s\beta}(x_2) – \chi_{1s\beta}(x_1)\chi_{1s\alpha}(x_2))\right)\left(-\frac{1}{2}\nabla_1^2\right)\left(\frac{1}{\sqrt{2}}(\chi_{1s\alpha}(x_1)\chi_{1s\beta}(x_2) – \chi_{1s\beta}(x_1)\chi_{1s\alpha}(x_2))\right) dx_1 dx_2 \]

これを展開すると、4つの項が現われます。\[ \langle{}T\rangle_1 = \frac{1}{2}\iint \chi_{1s\alpha}(x_1)\chi_{1s\beta}(x_2)\left(-\frac{1}{2}\nabla_1^2\right)\chi_{1s\alpha}(x_1)\chi_{1s\beta}(x_2)dx_1 dx_2 \] \[+ \frac{1}{2}\iint \chi_{1s\beta}(x_1)\chi_{1s\alpha}(x_2)\left(-\frac{1}{2}\nabla_1^2\right)\chi_{1s\beta}(x_1)\chi_{1s\alpha}(x_2)dx_1 dx_2\]\[- \frac{1}{2}\iint \chi_{1s\alpha}(x_1)\chi_{1s\beta}(x_2)\left(-\frac{1}{2}\nabla_1^2\right)\chi_{1s\beta}(x_1)\chi_{1s\alpha}(x_2)dx_1 dx_2 \]\[- \frac{1}{2}\iint \chi_{1s\beta}(x_1)\chi_{1s\alpha}(x_2)\left(-\frac{1}{2}\nabla_1^2\right)\chi_{1s\alpha}(x_1)\chi_{1s\beta}(x_2)dx_1 dx_2 \]

真ん中の演算子の部分は電子1にだけ作用するので、電子2に関する部分は別に積分できて、試行1s軌道関数は規格化されているので、αスピン同士またはβスピン同士だと積分は1、αスピンとβスピンの組み合わせだと積分=0になります。
\[ \int \chi_{1s\alpha}(x_2)\chi_{1s\alpha}(x_2) dx_2= \int \chi_{1s\beta}(x_2)\chi_{1s\beta}(x_2) dx_2= 1\]\[ \int \chi_{1s\alpha}(x_2)\chi_{1s\beta}(x_2) dx_2= \int \chi_{1s\beta}(x_2)\chi_{1s\alpha}(x_2) dx_2= 0\]
となって第3項と第4項は消せます。残った部分は

\[ \langle{}T\rangle_1 = \frac{1}{2}\int \chi_{1s\alpha}(x_1) \left(-\frac{1}{2}\nabla_1^2\right) \chi_{1s\alpha}(x_1) dx_1 + \frac{1}{2}\int \chi_{1s\beta}(x_1) \left(-\frac{1}{2}\nabla_1^2\right) \chi_{1s\beta}(x_1) dx_1. \]

スピン部分の積分はどちらも1となり、第1項と第2項は同じ値になります。したがって、\[ \langle{}T\rangle_1 =\int \phi_{1s}(r_1) \left(-\frac{1}{2}\nabla_1^2\right) \phi_{1s}(\mathbf r_1) d\tau_1. \]
今回の 1s 軌道は球対称なので、\[ \nabla_1^2 = \frac{1}{r_1^2}\frac{d}{dr_1}\left(r_1^2 \frac{d}{dr_1}\right) \]を使うことができます。
したがって、
\[ \langle{}T\rangle_1= \int_0^{2\pi}\!\!\!\int_0^{\pi}\!\!\!\int_0^{\infty} \sqrt{\frac{\zeta^3}{\pi}}e^{-\zeta r_1} \left[-\frac{1}{2} \frac{1}{r_1^2}\frac{d}{dr_1}\left(r_1^2 \frac{d}{dr_1}\right) \right] \sqrt{\frac{\zeta^3}{\pi}}e^{-\zeta r_1} r_1^2 \sin\theta dr_1 d\theta d\phi \] \[= 4\zeta^3 \int_0^{\infty} \left(\zeta{}r_1 – \frac{\zeta^2 r_1^2}{2}\right) e^{-2\zeta{}r_1} dr_1 \]\[ = 4\zeta^3\left(\frac{\zeta}{4\zeta^2} – \frac{\zeta^2}{8\zeta^3}\right) = \frac{\zeta^2}{2}.\]

電子2についてもまったく同じ値なので、\[\langle{}T\rangle_2 = \frac{\zeta^2}{2}\]となり、ヘリウム原子全体の運動エネルギーは\[\langle{}T\rangle_1 + \langle{}T\rangle_2 = \zeta^2\]です。

核引力エネルギーの期待値

核引力エネルギー演算子\(-\frac{2}{r_1}\)も電子1にしか作用しない一電子演算子なので、運動エネルギーの計算と同様に、スピン部分の積分と電子2についての積分を先にしてからまとめると、以下の式になります。

\[ \langle{}V\rangle_1 = 4\zeta^3 \int_0^{\infty} e^{-\zeta r_1} \left(-\frac{2}{r_1}\right) e^{-\zeta r_1} \times r_1^2 dr_1 \]
\[ = -8\zeta^3 \int_0^{\infty} {r_1} e^{-2\zeta r_1} dr_1 = -8\zeta^3 \times \frac{1}{4\zeta^2} = -2\zeta \]

電子2についてもまったく同じ値なので、\[\langle{}V\rangle_2 = -2\zeta\]となり、ヘリウム原子全体の核引力エネルギーは\[\langle{}V\rangle_1 + \langle{}V\rangle_2 = -4\zeta\]です。

電子-電子反発エネルギーの期待値

電子間反発エネルギーの期待値を取ると、

\[\langle V_{ee}\rangle = \iint \left(\frac{1}{\sqrt{2}}(\chi_{1s\alpha}(x_1)\chi_{1s\beta}(x_2) – \chi_{1s\beta}(x_1)\chi_{1s\alpha}(x_2))\right)\frac{1}{r_{12}} \left(\frac{1}{\sqrt{2}}(\chi_{1s\alpha}(x_1)\chi_{1s\beta}(x_2) – \chi_{1s\beta}(x_1)\chi_{1s\alpha}(x_2))\right) dx_1 dx_2 \]
です。これを展開すると4つの項が現われます。
\[\langle V_{ee}\rangle = \frac{1}{2}\iint \chi_{1s\alpha}(x_1)\chi_{1s\beta}(x_2)\frac{1}{r_{12}}\chi_{1s\alpha}(x_1)\chi_{1s\beta}(x_2) dx_1dx_2 \]\[+ \frac{1}{2}\iint \chi_{1s\beta}(x_1)\chi_{1s\alpha}(x_2)\frac{1}{r_{12}}\chi_{1s\beta}(x_1)\chi_{1s\alpha}(x_2) dx_1dx_2 \]\[- \frac{1}{2}\iint \chi_{1s\alpha}(x_1)\chi_{1s\beta}(x_2)\frac{1}{r_{12}}\chi_{1s\beta}(x_1)\chi_{1s\alpha}(x_2) dx_1dx_2 \]\[- \frac{1}{2}\iint \chi_{1s\beta}(x_1)\chi_{1s\alpha}(x_2)\frac{1}{r_{12}}\chi_{1s\alpha}(x_1)\chi_{1s\beta}(x_2) dx_1dx_2 \]

最初の2項は、電子がそれぞれの軌道に分布していることによる通常のクーロン反発を表すクーロン項です。一方、後ろの2項は、スレーター行列式を反対称化したことによって生じた交換項です。

第3項と第4項では、スピン部分の積分が0\[ \int\alpha^*(\sigma)\beta(\sigma)\,d\sigma=0\]になります。したがって、第3項と第4項も0になります。

\[K_{\alpha\beta}=0.\]

一方、第一項と第二項では、スピン部分の積分が1\[ \int\alpha^*(\sigma)\alpha(\sigma)\,d\sigma=\int\beta^*(\sigma)\beta(\sigma)\,d\sigma=1\]になり、さらに電子1と電子2は同じ空間軌道 \(\phi_{1s}\) を占めているので、2つのクーロン項は等しくなります。そのため、\[\langle V_{ee}\rangle = \frac{1}{2}J_{12}+\frac{1}{2} J_{21} = J\]となり、
\[\langle V_{ee}\rangle = \iint\phi^*_{1s}(r_1)\phi^*_{1s}(r_2)\frac{1}{r_{12}} \phi_{1s}(r_1)\phi_{1s}(r_2) d\tau_1 d\tau_2 \]が得られます。

このように、スレーター行列式から出発すると一般には「\(クーロン項 – 交換項\)」という構造が自然に現われますが、ヘリウム原子の基底状態では2電子のスピンが異なるため、両電子間の交換項は0になります。ここから先は、空間部分だけを使って計算することができます。

電子2について先に積分する形に書き直すと、\[\langle V_{ee}\rangle = \int\phi^*_{1s}(\mathbf r_1)\left[\int\frac{{|\phi_{1s}(\mathbf r_2)|^2}}{r_{12}}d\tau_2\right] \phi_{1s}(\mathbf r_1) d\tau_1 . \]

真ん中の
\[ J(\mathbf r_1) = \int \frac{|\phi_{1s}(\mathbf r_2)|^2}{r_{12}} d\tau_2 \]
の部分は、\(r_1,r_2\)の関数を\(r_2\)について積分するので、\(r_1\)だけの関数になり、電子1が、試行軌道\(\phi_{1s}\)に分布している電子2から感じる平均的なクーロンポテンシャルと解釈できます。\(\sqrt{\frac{\zeta^3}{\pi}}e^{-\zeta r_2}\)を代入すると、
\[ J(\mathbf r_1) = \frac{\zeta^3}{\pi} \int_0^{2\pi}\!\!\!\int_0^{\pi}\!\!\!\int_0^{\infty} \frac{e^{-2\zeta r_2}}{r_{12}} r_2^2 \sin\theta dr_2 d\theta d\phi . \]
ここで\(r_1\)の方向を\(z\)軸に取れば、\(\theta\)はベクトル\(\mathbf r_1\)と\(\mathbf r_2\)がなす角であり、\[ r_{12} = \sqrt{r_1^2+r_2^2-2r_1r_2\cos\theta} \]となります。
この積分は高校数学でもがんばれば解くことができて（ワシントン大学の講義資料のpp.139-140を参照）、
\[ J(\mathbf r_1) = \frac{1-e^{-2\zeta r_1}}{r_1} – \zeta e^{-2\zeta r_1} \]
と、\(r_1\)だけの関数になりました。したがって、
\[\langle V_{ee}\rangle = \int_0^{2\pi}\!\!\!\int_0^{\pi}\!\!\!\int_0^{\infty} \sqrt{\frac{\zeta^3}{\pi}}e^{-\zeta r_1} \left(\frac{1-e^{-2\zeta r_1}}{r_1} – \zeta e^{-2\zeta r_1}\right) \sqrt{\frac{\zeta^3}{\pi}}e^{-\zeta r_1} r_1^2 \sin\theta dr_1 d\theta d\phi \]
の積分を計算をすると\[ \langle V_{ee}\rangle = \frac{5}{8}\zeta \]が得られます。

全エネルギーの期待値

全エネルギーの期待値は
\[ E = 2\langle{}T\rangle + 2 \langle{}V\rangle + \langle V_{ee}\rangle = \zeta^2 – 4\zeta + \frac{5}{8}\zeta = \zeta^2 – \frac{27}{8}\zeta\]
となります。この関数は上に開いた二次関数なので、最も低い\(E\)を与えるのは、\(E\)の傾きが0になる
\[ \frac{dE}{d\zeta} = 0\]
\(\zeta\)です。\(E\)を微分すると
\[ 2\zeta- \frac{27}{8}= 0 \]
\[\zeta = \frac{27}{16} \approx 1.6875.\]
この値を代入すると、

\[ E = -\frac{729}{256} \approx -2.848 E_{\mathrm h}\]

が得られます。高精度計算で得られたヘリウム原子のエネルギーが\(E \approx -2.9037 E_{\mathrm h}\) なので、比較的精度のよい値です。以上の方法は1s軌道を \[ \phi_{1s}(\mathbf r)= N(\zeta)e^{-\zeta r}\] という形に限定した、1パラメータのハートリー・フォック型変分計算でした。こで、\(N(\zeta)\) は規格化定数です。これに対して軌道関数 \(\phi\) の形そのものを変化させるハートリー・フォック法では、通常の微分の代わりに汎関数微分を使います。

ハートリー・フォック（Hartree–Fock）法

本来のハートリー・フォック法では、軌道を特定の関数形に限定せず、軌道\(\phi(\mathbf r)\)の形そのものを変化させて、エネルギーが停留する軌道を探します。このとき、エネルギー\(E\)を数値\(x\)の関数ではなく、関数\(\phi(\mathbf r)\)を入力すると数値を返すものになります。このような「関数を入力として数値を返すもの」を汎関数と呼び、\[ E[\phi] \] と書きます。


規格化条件とラグランジュの未定乗数

軌道を最適化する時には、軌道が規格化 \(\int |\phi(\mathbf r)|^2\,d\tau= 1\)（一電子軌道の存在確率を全空間で積分すると1になる）されていなければいけません。加えて、異なる軌道が直交していること\(\int \chi_i^*(x)\chi_j(x) dx = 0\)（異なる軌道\(i\neq j\)について軌道を掛け合わせて積分すると0になる）という規格直交条件を満たす必要があります。
この規格直交条件を保ちながらエネルギーを最適化するために、\(\varepsilon_{ij}\)というラグランジュの未定乗数を導入します。

ヘリウム原子の基底状態ではαスピン軌道とβスピン軌道の直交性が自動的に満たされるので、\[\int |\phi(\mathbf r)|^2\,d\tau= 1\]という1つの規格化条件だけを考えれば十分です。

「ラグランジュの未定乗数法」の説明をします。例えば \(f(x,y) = x + 2y\)を\(x^2+y^2 = 1\)（半径1の円）という制約を課して最大最小化する時は、\[ L(x,y,\lambda) = x + 2y -\lambda(x^2 + y^2 -1) \]という量を作って、
\[ \frac{\partial L(x,y,\lambda)}{\partial x}=0, \frac{\partial L(x,y,\lambda)}{\partial y}=0, \frac{\partial L(x,y,\lambda)}{\partial \lambda}=0 \] から得られる連立方程式を解くことで答えが得られます（最大値は\(\sqrt{5}\)、最小値は\(-\sqrt{5}\)）。

ヘリウム原子では、ラグランジュの未定乗数 \(\varepsilon\) を導入して、\[L[\phi] = E[\phi] – 2\varepsilon\left(\int|\phi(\mathbf r)|^2\,d\tau – 1 \right) \]という新しい汎関数を作ります。\(\varepsilon\)の係数2は、ヘリウム原子では同じ空間軌道を2電子が占有していることに対応させたものです。この係数を付けておくと、後で現れる \(\varepsilon\) を軌道エネルギーとしてそのまま解釈しやすくなります。

\(E[\phi]\)はエネルギー、\(\varepsilon\) の後ろの括弧の中は無次元なので、\(\varepsilon\) はエネルギーと同じ次元を持ちます。

「軌道を変分する」とはどういうことか

エネルギーを\(\phi^*\)で変分するのは、関数\(f(x)\)を最小化する時に変数\(x\)を\(x \longrightarrow x + \delta x\)と少しだけ変化させて、その時の \(f\)の変化を考えるのと考え方は同じです。ハートリー・フォック法では関数を少しだけ変化させます。
\[\phi(\mathbf r) \longrightarrow \phi(\mathbf r) + \delta\phi(\mathbf r)\]

エネルギーの変化を、この小さな変化量について展開すると、
\[ E[\phi+\delta\phi]-E[\phi] = \delta E+\text{二次以上の小さい項} \]
と書けます。\(\delta E\) は、\(\delta\phi\) または \(\delta\phi^*\) について一次の項を集めた「一次変分」です。

汎関数微分は、この一次変分を
\[ \delta E =\int \frac{\delta E}{\delta\phi(\mathbf r)}\,\delta\phi(\mathbf r)\,d\tau + \int \frac{\delta E}{\delta\phi^*(\mathbf r)}\,\delta\phi^*(\mathbf r)\,d\tau \]
と書いたときの係数に相当します。変分計算上は\(\phi\)と\(\phi^*\)を独立な変数として扱ったうえで、\(\phi\)は固定して、\(\phi^* \longrightarrow \phi^* + \delta\phi^*\) とすれば汎関数微分\(\frac{\delta L}{\delta\phi^*}\) を導くことができます。

ヘリウム原子のエネルギー汎関数

ヘリウム原子で実際に計算してみます。1つの電子にだけ作用する運動エネルギー演算子と核引力エネルギー演算子をまとめて
\[ \hat{h} = -\frac{1}{2}\nabla^2 – \frac{2}{r} \]
と書きます。ヘリウム原子では、同じ空間軌道\(\phi\)を2電子が占有するので、一電子部分のエネルギーは、\[ E_{1e} = 2\int \phi^*(\mathbf r)\hat{h}\phi(\mathbf r)\,d\tau \]と書けます。

電子間反発のエネルギーは以前求めた通り、\[E_{ee} = \iint\frac{|\phi(\mathbf r_1)|^2 |\phi(\mathbf r_2)|^2}{r_{12}} d\tau_1 d\tau_2\]です。

したがって、ヘリウム原子のエネルギー汎関数は\[ E[\phi] = 2\int \phi^*(\mathbf r)\hat{h}\phi(\mathbf r)\,d\tau +\iint\frac{|\phi(\mathbf r_1)|^2 |\phi(\mathbf r_2)|^2}{r_{12}} d\tau_1 d\tau_2 \]
です。これに規格化条件を加えた\[L = 2\int \phi^*(\mathbf r)\hat{h}\phi(\mathbf r)\,d\tau +\iint\frac{|\phi(\mathbf r_1)|^2 |\phi(\mathbf r_2)|^2}{r_{12}} d\tau_1 d\tau_2 – 2\varepsilon\left(\int|\phi(\mathbf r)|^2\,d\tau – 1 \right) \]を変分します。

一電子項

\(\phi^* \longrightarrow \phi^* + \delta\phi^*\) と変化させた\[ 2\int(\phi^* + \delta\phi^*)\hat{h}\phi d\tau \]を展開すると、\[ = 2\int\left[\phi^*\hat{h}\phi + \delta\phi^*\hat{h}\phi\right]d\tau . \]

一次変分として残るのは、\[\delta E_{1e} = 2\int\delta\phi^*\hat{h}\phi d\tau \]です。

電子間反発項

\[ (\phi^*(\mathbf r_1)+\delta\phi^*(\mathbf r_1))\phi(\mathbf r_1)(\phi^*(\mathbf r_2) +\delta\phi^*(\mathbf r_2))\phi(\mathbf r_2) \]の部分を先に展開すると、\[ |\phi_1|^2|\phi_2|^2 + (\delta\phi_1^*\phi_1)|\phi_2|^2 + |\phi_1|^2(\delta\phi_2^*\phi_2) + (\delta\phi_1^*\phi_1) (\delta\phi_2^*\phi_2).\]

したがって、一次の部分\(\delta E_{ee}\) は
\[ \delta E_{ee} = \iint \frac{\delta\phi_1^*\phi_1|\phi_2|^2}{r_{12}} d\tau_1d\tau_2 + \iint \frac{|\phi_1|^2\delta\phi_2^*\phi_2}{r_{12}} d\tau_1d\tau_2 .\]

電子1と電子2は同じ空間軌道\(\phi\)を占めています。また、\(r_{12}=r_{21}\) であり、積分変数 \(\mathbf r_1\) と \(\mathbf r_2\) は入れ換えることができます。したがって、
\[ \delta E_{ee} =2\iint \frac{\delta\phi_1^*\phi_1|\phi_2|^2}{r_{12}}d\tau_1 d\tau_2 \] \[ = 2\int\delta\phi^*(\mathbf r_1)\left[\int\frac{|\phi(\mathbf r_2)|^2}{r_{12}}d\tau_2\right]\phi(\mathbf r_1)d\tau_1. \]

ここで、\[ J(\mathbf r_1) = \int\frac{|\phi(\mathbf r_2)|^2}{r_{12}}d\tau_2 \]と定義します。これは、電子1が位置 \(\mathbf r_1\) にいるとき、軌道\(\phi\)に分布しているもう1つの電子から感じる平均的なクーロンポテンシャルです。

規格化条件

\[ – 2\varepsilon\left(\int (\phi^*+\delta\phi^*)\phi\, d\tau – 1\right)\] を展開すると、
\[ -2\varepsilon\left[\int \left(|\phi|^2 + \delta\phi^*\phi \right) d\tau – 1\right]. \] したがって、一次の項は \[ \delta L_{規格化} = -2\varepsilon\int \delta\phi^*\phi d\tau.\]

一電子項、電子間反発項、規格化条件をまとめる

\(\delta\phi^*\)についての一次の項、つまり第一変分をまとめると、
\[ \delta L =2\int\delta\phi^*\hat{h}\phi d\tau + 2\int\delta\phi^*(\mathbf r_1)\left[\int\frac{|\phi(\mathbf r_2)|^2}{r_{12}}d\tau_2\right]\phi(\mathbf r_1)d\tau_1\]\[ -2\varepsilon\int \delta\phi^*\phi d\tau. \]

普通の関数\(f(x)\) では、\(x\)をわずかに変化させる（\(x \longrightarrow x + \delta x\)）と一次の変化は\(\delta f = \frac{df}{dx}\delta x\)となり、\(\delta x\)に掛かっている係数が微分\(\frac{df}{dx}\)です。これと同じように\(\delta L =\int \delta\phi^*(\mathbf r)\frac{\delta L}{\delta \phi^*(\mathbf r)}\,d\tau\)と書いた時、\(\delta\phi^*\) に掛かっている係数が、汎関数微分です。したがって、\(\int\delta\phi^*\hat{h}\phi d\tau \)では、\(\delta\phi^*\)に掛かっている \(\hat{h}\phi\)を取り出して、\(\frac{\delta}{\delta \phi^*}\int\phi^*\hat{h}\phi\, d\tau = \hat{h}\phi\)となります。

\[ \frac{\delta L}{\delta\phi^*} = 2\left[\hat{h}\phi + \left(\int\frac{|\phi(\mathbf r_2)|^2}{r_{12}}d\tau_2\right)\phi -\varepsilon\phi\right] \] 停留条件\(\frac{\delta L}{\delta\phi^*} = 0 \) を課すと\(\frac{\delta L}{\delta\phi^*}\) 全体を2で割ることができて、さらに移項すると、 \[ \left( \hat{h}+ \int\frac{|\phi(\mathbf r_2)|^2}{r_{12}}d\tau_2\right)\phi = \varepsilon\phi \] が得られます。これは、関数\(\phi\)に \(\left( \hat{h}+ \int\frac{|\phi(\mathbf r_2)|^2}{r_{12}}d\tau_2\right)\) という演算をしたら、演算した結果が元の関数\(\phi\)の\(\varepsilon\)倍になるという固有値方程式になっています。普通と違うのは演算子の中に答えの\(\phi\)が含まれているので、矛盾がなくなるまで反復的に解かなければいけません。そして、規格化条件 \(\int \phi^*\phi d\tau = 1\) と連立させることで、規格化した固有関数を求めることになります。 \[\int\frac{|\phi(\mathbf r_2)|^2}{r_{12}}d\tau_2 = J(\mathbf r_1) \] と置き換えて、さらに\( \hat J\phi(\mathbf r_1) = J(\mathbf r_1)\phi(\mathbf r_1) \)と定義すれば、固有値方程式は \[ \left( \hat{h}+\hat{J} \right)\phi = \varepsilon\phi \] と表わされます。ここで得られた \((\hat h+\hat J_1)\phi_1=\epsilon_1\phi_1\) は、He基底状態で2電子が共有する占有空間軌道を求めるハートリー・フォック方程式です。Heのように占有空間軌道が1つだけの閉殻RHFでは、標準的なフォック演算子は \(\hat f=\hat h+2\hat J_1-\hat K_1\) です。占有軌道（被占軌道） \(\phi_1\) に対しては \(\hat K_1\phi_1=\hat J_1\phi_1\) なので、上の方程式に簡約されます。

ヘリウム原子のハートリー・フォック方程式を数値的に解く

ヘリウム原子の（非相対論的で原子核を固定した場合の）厳密な全エネルギーは \(E_{exact} = -2.9037\, E_{\mathrm h}\)です。ハートリー・フォック法では、多電子波動関数を単一のスレーター行列式で近似しているので、変分原理からHFエネルギーはこの値より高くなるはずです。しかし、HF限界における軌道\(\phi (r)\)がどのような単純な解析関数で表わせるかは、あらかじめ分かりません。少なくとも、軌道を単一の指数関数 \(\phi(r)\propto e^{-\zeta r}\)  に限定した変分計算では、\(\zeta=27/16\) のとき\(E =-2.8477\, E_{\mathrm h}\)が最低値でした。

そこで、距離\(r\)を\(0\,a_0\)から\(10 \,a_0\)まで\(0.01\)刻みに区切って、全部で1001個の格子点を使ってHF方程式を数値的に解くPythonコードをChatGPTに生成してもらいました。計算を簡単にするために動径波動関数\(\phi(r)\) そのものではなく、これに\(r\)を掛けた換算動径波動関数\(u(r) = r \phi(r)\)を最適化しているそうです。\(u(r)\)の一階微分（傾き）は「\((u_{i+1} -u_{i-1})/(2 \Delta r\))」と両隣りの値から求めて、\(u(r)\)の二階微分（曲率）は『右側の傾きと左側の傾きの差を、さらに \(\Delta r\) で割る』\( [(u_{i+1}-u_i) – (u_{i}-u_{i-1})]/\Delta r^2\) で求めます。初期値は\(e^{-\zeta r}\)（\(\zeta = \frac{27}{16}\)）から割り当てられますが、軌道の値は反復計算の間に更新されるので、最終的な軌道形は指数関数には制限されません。
計算の結果、\[\Delta r = 0.01\, a_0: E_{HF} \approx -2.86135\, E_{\mathrm h}\]となりました。さらに格子を細かくすると、\[\Delta r = 0.001\, a_0: E_{HF} \approx -2.86168\, E_{\mathrm h}\]と既知のハートリー・フォック限界\(E_{\mathrm{HF}}^{\mathrm{limit}} = -2.861679995612\,E_h\)とほぼ一致しました。\(E_{exact}\)と\(E_{\mathrm{HF}}^{\mathrm{limit}}\)の差が「電子相関エネルギー」\(E_{\mathrm{corr}}\)\[E_{\mathrm{corr}} = E_{\mathrm{exact}}- E_{\mathrm{HF}}^{\mathrm{limit}}\]
と定義されますが、それはまた別のお話。

ヘリウム原子のハートリー・フォック方程式をローターン方程式にして解く

数値的に格子点上で解くのではなく、やっぱり方程式を解いて軌道を求めたいところです。しかし、「軌道の形を変えて、、エネルギーが最小になる形を探す」と言われても、関数そのものをどう変化させて探せばいいのか、簡単には分かりません。そこで、実際の量子化学計算では、軌道をあらかじめ用意した複数の関数の足し合わせとして表現します。
\[ \phi =C_1 \chi_1 + C_2 \chi_2 \cdots \] つまり、「軌道の形を変える」という問題を「係数\(C_1, C_2, \cdots\)を変える」という問題に置き換えます。こうすると、もともとは関数を未知数とする微分方程式だったハートリー・フォック方程式を、係数\(C_1, C_2, \cdots\)を未知数とする連立方程式（行列の固有値問題）として表わすことができます。行列は以前は高校の数Cの内容に入っていたのですが、今はあまりやらないようです。この形にしたハートリー・フォック方程式をローターン（Roothaan）方程式またはローターン・ホール（Roothaan–Hall）方程式と呼びます。

足し合わせる関数の形としては、ガウス関数\[ g_i(r) = \left(\frac{2\alpha_i}{\pi}\right)^{3/4} e^{-\alpha_i r^2}\]がよく使われます。ガウス関数を使う大きな理由は、ガウス関数を使うことで異なる原子核を中心とする関数の積分（多中心積分）を解析的に効率よく計算できるためです。足し合わせるためにあらかじめ用意した関数のセットを基底関数系と呼びます。通常、それぞれの基底関数は原子核の位置を中心として配置されます。

この記事では、ジョン・ポープル（John Pople）による3-21G基底関数系を使用します。Heの3-21G基底では、s型基底関数を2本用います。1本目は2個の原始ガウス関数を縮約した関数、2本目は1個の原始ガウス関数からなる関数です。HF軌道は、この2本の基底関数の線形結合 \[ C_1\chi_1(r)+C_2\chi_2(r) \]として表わします。具体的に数字を書くと、\[= C_1 \left[ 0.885750e^{-13.6267r^2} + 1.070687e^{-1.99935r^2} \right] + C_2 \left[ 0.346979e^{-0.382993r^2} \right] \]となります。

原始ガウス関数は全部で3つですが、この計算で変化させるのは、それぞれの基底関数\(\chi_1,\chi_2\)をどれだけ混ぜるかを表わす2つの係数\(C_1,C_2\)です。2つ前の節では\[\phi(r)=Ne^{-\zeta r}\]の指数 \(\zeta\) を変化させてエネルギーが最低になる軌道を探しました。それに対して、3-21G基底を用いるローターン法では、ガウス関数の指数 \(\alpha_i\) や縮約係数はあらかじめ決めておき、分子軌道を作る係数 \(C_1,C_2\) を変化させます。

規格化

軌道 \(\phi_{1s}^{\mathrm{HF/3\text{-}21G}}(r)\) について、一電子の確率密度 \(|\phi_{1s}(r)|^2\) を全空間で積分すると1になるように規格化します。
\[ \int|\phi_{1s}|^2 d\tau = 1\]
\[\int (C_1\chi_1(r)+C_2\chi_2(r))^2\,d\tau = 1.\]

\((A + B)^2 = A^2 + 2AB+B^2\)なので、\[ C_1^2\int\chi_1^2d\tau + 2C_1C_2\int\chi_1\chi_2d\tau + C_2^2\int\chi_2^2d\tau =1. \]
ここで、基底関数同士の重なり積分を\(S_{ij} = \int \chi_i\chi_j d\tau \)と定義します。今回使っている基底関数\(\chi_1,\chi_2\)はそれぞれ規格化されているので、\[S_{11} =S_{22} = 1\]です。問題は2つの異なる基底関数の重なり
\[\int\chi_1(\mathbf r)\chi_2(\mathbf r)\,d\tau\]です。

\[ d\tau = r^2\sin\theta dr d\theta d\phi \]なので、\[ S_{12} = \int_0^{2\pi}\!\!\!\int_0^{\pi}\!\!\!\int_0^{\infty} \left[ 0.885750e^{-13.6267r^2} + 1.070687e^{-1.99935r^2} \right] \times \]\[ \left[ 0.346979e^{-0.382993r^2} \right] r^2\sin\theta dr d\theta d\phi \]

角度部分を先に積分すると、\[ S_{12} = 4\pi \left[ (0.885750 \times 0.346979) \int_0^{\infty} r^2 e^{-(13.6267+0.382993)r^2} dr + (1.070687 \times 0.346979) \int_0^{\infty} r^2 e^{-(1.99935+0.382993)r^2}dr\right]. \]

ここで積分公式\[ \int_0^{\infty} r^2 e^{-ar^2} dr = \frac{\sqrt{\pi}}{4a^{3/2}} \]
を使うと、\[ S_{12} = 4\pi(0.885750)(0.346979) \frac{\sqrt{\pi}} {4(14.009693)^{3/2}} + 4\pi(1.070687)(0.346979) \frac{\sqrt{\pi}} {4(2.382343)^{3/2}} \]
\[ = 0.03264 + 0.56258 \approx 0.5952. \]

数値を代入すると、\[ C_1^2 + 2(0.5952)C_1C_2 + C_2^2 =1. \]

これを行列の形で書くと\[ \mathbf C^{\mathrm T}\mathbf S\mathbf C =1\]です。 \(\mathbf C^{\mathrm T}\)は転置行列で、行列の行と列を引っくり返して作ります。具体的には、\[ \begin{pmatrix} C_1 & C_2\end{pmatrix} \begin{pmatrix} S_{11} & S_{12} \\ S_{21} & S_{22} \end{pmatrix} \begin{pmatrix} C_{1} \\ C_{2} \end{pmatrix}=1. \]すなわち、\[ \begin{pmatrix} C_1 & C_2\end{pmatrix} \begin{pmatrix} 1 & 0.5952 \\ 0.5952 & 1 \end{pmatrix} \begin{pmatrix} C_{1} \\ C_{2} \end{pmatrix} = 1 \]
となります。

この真ん中に挟まれた重なり行列 \(\mathbf S\) が、後でローターン方程式\[\mathbf F\mathbf C = \mathbf S\mathbf C\boldsymbol{\varepsilon}\]に現れます。

ローターン方程式を作る

ヘリウム原子のハートリー・フォック方程式は\[ \left( \hat{h}+ \int\frac{|\phi(\mathbf r_2)|^2}{r_{12}}d\tau_2\right)\phi = \varepsilon\phi \]でした。フォック演算子を使えば、 \[\hat{F}\phi = \varepsilon\phi \] と簡潔に書くことができます。ここで\[C_1\chi_1(r)+C_2\chi_2(r)\] をハートリー・フォック方程式に代入すると、 \[ \hat{F}(C_1\chi_1(r)+C_2\chi_2(r)) = \varepsilon(C_1\chi_1(r)+C_2\chi_2(r)) \]となります。この式の両辺に左から\(\chi_1(r)\)を掛けて全空間で積分します。\[ \int\chi_1(r)\hat{F} C_1\chi_1(r) d\tau+\int \chi_1(r)\hat{F} C_2\chi_2(r) d\tau = \int \chi_1(r)\varepsilon C_1\chi_1(r) d\tau+\int\chi_1(r)\varepsilon C_2\chi_2(r) d\tau.\] 係数 \(C_1,C_2\)と固有値\(\varepsilon\) は位置に依存しない数なので、積分の外に出すことができます。\[ C_1\int\chi_1(r)\hat{F}\chi_1(r) d\tau+C_2\int \chi_1(r)\hat{F}\chi_2(r)) d\tau = \varepsilon C_1 \int \chi_1(r)\chi_1(r) d\tau+\varepsilon C_2\int\chi_1(r)\chi_2(r)) d\tau.\]

少しややこしいのですが、これまで\(1,2\)という数字は電子1と電子2を指すラベルとして使ってきました。しかし、ここからは基底関数を指すラベルとして使っています。

ここで、\( F_{11} = \int\chi_1(r)\hat{F}\chi_1(r) d\tau \) をフォック行列要素、\( S_{11} = \int \chi_1(r)\chi_1(r) d\tau \)を重なり行列要素と定義します。すると、上の式は、\[ C_1 F_{11} + C_2 F_{12} = \varepsilon\left( C_1 S_{11} + C_2 S_{12} \right)\]と簡潔に書くことができます。前の節で既に計算したように、\(S_{11} = 1,\, S_{12} = 0.5952\)です。

同じように、今度はもとのハートリー・フォック方程式に左から \(\chi_2\) を掛けて積分します。途中を省略しますが、同様に、\[ C_1 F_{21} + C_2 F_{22} = \varepsilon\left( C_1 S_{21} + C_2 S_{22} \right)\] が得られました。これによって未知の関数\(\phi\)を求める微分方程式が、基底関数の係数 \(C_1\)と\(C_2\)を求める連立方程式に置き換わりました。

\[ C_1 F_{11} + C_2 F_{12} = \varepsilon\left( C_1 S_{11} + C_2 S_{12} \right)\]\[ C_1 F_{21} + C_2 F_{22} = \varepsilon\left( C_1 S_{21} + C_2 S_{22} \right)\]

連立一次方程式は行列の形で書くことができます。行列の形で書くと、\[ \begin{pmatrix} F_{11} & F_{12} \\ F_{21} & F_{22} \end{pmatrix} \begin{pmatrix} C_1 \\ C_2 \end{pmatrix} = \varepsilon \begin{pmatrix} S_{11} & S_{12} \\ S_{21} & S_{22} \end{pmatrix} \begin{pmatrix} C_1 \\ C_2 \end{pmatrix} \]となります。これが一つの軌道について書いたローターン方程式です。

次に、すべての項を左辺に移して \(=0\) の形にします。\[ C_1 (F_{11} – \varepsilon S_{11}) + C_2(F_{12} – \varepsilon S_{12}) = 0\]\[ C_1 (F_{21} – \varepsilon S_{21}) + C_2(F_{22} – \varepsilon S_{22}) = 0\]行列で書けば、 \[ \begin{pmatrix} F_{11} – \varepsilon S_{11} & F_{12} – \varepsilon S_{12} \\ F_{21} – \varepsilon S_{21} & F_{22} – \varepsilon S_{22} \end{pmatrix} \begin{pmatrix} C_1 \\ C_2 \end{pmatrix} = 0. \]

この連立方程式には\(C_1 = C_2 = 0\)という解が必ずあります。しかしそれでは、\(\psi = 0\)となって意味がありません。そこで、\(C_1,C_2\)が同時に0でない解を持つための条件を考えます。そのためには、係数行列の行列式が0でなければなりません。\[\begin{vmatrix}F_{11}-\varepsilon & F_{12}-0.5952\varepsilon \\ F_{21}-0.5952\varepsilon & F_{22}-\varepsilon \end{vmatrix}=[(F_{11}-\varepsilon) \times ( F_{22}-\varepsilon) ]-[(F_{12}-0.5952\varepsilon)\times ( F_{21}-0.5952\varepsilon)] = 0\]が条件です。

また、行列式が出てきました。以前出てきたスレーター行列式は、電子交換に対して反対称な多電子波動関数を簡潔に表すためのものでした。一方、ここで出てきた行列式は、連立方程式が「\(C_1=C_2= 0\)」以外の解を持つための条件です。同じ「行列式」という数学的な道具ですが、役割はまったく異なります。式という名前が付いていますが行列式は数値です。行列式は高校数学に登場しませんが、2×2の正方行列\(\begin{pmatrix}A & B \\ C & D \end{pmatrix}\)に対する行列式は\(A\times D – B\times C\)です。

ここで、適当な\(C_1,C_2\)を仮定してあらかじめ計算しておいて\( F_{11} , F_{12} ,F_{21} , F_{22}\)の数値を代入します。すると、上の式は\(\varepsilon\)に関する二次方程式になります。これを二次方程式の解の公式\[x = \frac{-b \pm\sqrt{b^2 – 4ac}}{2a}\]を使って解くと、2個の基底関数を使っているので、通常は2つの固有値\(\varepsilon_1, \varepsilon_2\)が得られます。さらに、それぞれの \(\varepsilon\) を連立方程式に戻せば、それぞれについて\(C_1\)と\(C_2\)の比\[\frac{C_1}{C_2}\]が決まります。最後に、規格化条件\[C_1^2+2(0.5952)C_1C_2+C_2^2=1\]を使うと、\(C_1\) と \(C_2\) の値そのものが決まります。

このような手順で2つの固有値\(\varepsilon_1,\varepsilon_2\)に対応して、\[\phi_{\mathrm{occupied}} = C_{11}\chi_1+C_{21}\chi_2\]\[\phi_{\mathrm{unoccupied}} = C_{12}\chi_1+C_{22}\chi_2\]という2つの軌道が得られます。ヘリウム原子の基底状態では、このうちエネルギーの低い方の軌道 \(\phi_{\mathrm{occupied}}\) に2電子を入れます。

ここまでの、\[\text{行列式}=0\]から固有値を求め、その固有値から係数を求めるという手順は、単純ヒュッケル法と基本的に同じです。

単純ヒュッケル法とハートリー・フォック法の大きな違いは、\(F_{ij}\)を計算するフォック演算子が
\[ \hat{F} = \hat{h} + \int\frac{|\phi(\mathbf r_2)|^2}{r_{12}}d\tau = \hat{h} + \int\frac{|C_1\chi_1(r_2)+C_2\chi_2(r_2)|^2}{r_{12}}d\tau \]と、解\(C_1,C_2\)を含んでいる点にあります。したがって、係数\(C_{11},C_{21}\)の値が更新される→\(F_{ij}\)の値も更新される→連立方程式を解いて係数\(C_{11},C_{21}\)の値を更新→\(F_{ij}\)の値を更新・・・、と\[(C_{11},C_{21}){\mathrm{新}}\approx(C_{11},C_{21}){\mathrm{前回}}\]となるまで反復計算しなければいけません。これが自己無撞着場（self-consistent field, SCF）計算です。

※ \(\phi_{\mathrm{unoccupied}}\)は非占有軌道（空軌道）ですが、HOMO/LUMOを議論する時に使われるハートリー・フォック法でいう正準仮想軌道（canonical virtual orbital）ではありません。ヘリウムの標準的な閉殻HFフォック演算子は\(\hat{f} = \hat{h} + 2\hat{J}_1 – \hat{K}_1\)です。占有軌道\(\phi_1\)に対しては\(\hat{J}\phi_1 = \hat{K}\phi_1\)なのに対して、別の軌道\(\phi_a\)に対しては一般に\(\hat{J}\phi_a \neq \hat{K}\phi_a\)です。そのため、正準仮想軌道を求めるには、収束した占有軌道から標準閉殻フォック行列 \(\hat h+2\hat J_1-\hat K_1\) を新たに組み立てます。

ローターン方程式を解く

SCF計算、1回目（\(C_1 = 1, C_2= 0\)）
フォック行列を計算する

ここでは、前節で導いたHeの占有空間軌道に対するHF方程式 \((\hat h+\hat J)\phi=\epsilon\phi\) を自己無撞着場（SCF）法で解きます。この節で用いる行列を \(F_{ij}\) と書きます。

\[ F_{ij}^{(n)} =\int\chi_i(\mathbf r_1) \left[ −\frac{1}{2}\nabla_1^2 –\frac{2}{r_1}+ \int\frac{| C_1^{(n)}\chi_1(\mathbf r_2)+C_2^{(n)}\chi_2(\mathbf r_2)|^2}{r_{12}}d\tau_2\right]\chi_j(\mathbf r_1) d\tau_1 \]

今回は、\(C_1 = 1, C_2= 0\)という極端な初期値を使ってみます。したがって、\[\phi^{(0)} = 1\times \chi_1(\mathbf r) + 0 \times \chi_2(\mathbf r) = 0.885750e^{−13.6267r^2}+1.070687e^{−1.99935r^2}\]
\[F_{ij}^{(0)}= \int\chi_i^*(\mathbf r_1) \left( −\frac{1}{2}\nabla_1^2 –\frac{2}{r_1}+ \int\frac{| \phi^{(0)}(\mathbf r_2)|^2}{r_{12}}d\tau_2\right)\chi_j(\mathbf r_1) d\tau_1\] です。

ガウス関数の積分公式を使って積分すると、運動エネルギー項、核引力項、電子間反発項の順に、
\[F_{11}^{(0)} = 3.91612-5.48968+1.80325=0.22969 \]
\[F_{12}^{(0)} = 0.57895-2.23529+0.87651=-0.77984\]
\[F_{21}^{(0)} = 0.57895-2.23529+0.87651=-0.77984 \]
\[F_{22}^{(0)} = 0.57449-1.97513+0.91701=-0.48363 \]
となります。したがって、1回目の行列は、\[\mathbf{F}^{(0)} = \begin{pmatrix} 0.22969 & -0.77984 \\ -0.77984 & -0.48363\end{pmatrix}\]となります。

解くべきローターン方程式は、\[ \begin{pmatrix} 0.22969 & -0.77984 \\ -0.77984 & -0.48363\end{pmatrix}\begin{pmatrix} C_1^{(1)}\\ C_2^{(1)} \end{pmatrix} =\begin{pmatrix}1 & 0.5952 \\ 0.5952 & 1\end{pmatrix}\begin{pmatrix}C_1^{(1)}\\ C_2^{(1)} \end{pmatrix}\varepsilon \]です。左辺にまとめて、行列式=0の条件を立てます。
\[ \begin{vmatrix} 0.22969 – \varepsilon & -0.77984 -0.5952\varepsilon\\ -0.77984-0.5952\varepsilon & -0.48363-\varepsilon\end{vmatrix}= 0\]
展開すると、
\[(0.22969 – \varepsilon)( -0.48363 – \varepsilon) – (-0.77984 – 0.5952\varepsilon)^2 = 0 \]\[ 0.64573696\varepsilon^2 -0.674381536\varepsilon -0.7192354003=0 \]となり、この\(\varepsilon\)の二次方程式を解くと、\[\varepsilon_1 \approx -0.65531\,E_h\]\[\varepsilon_2 \approx 1.69965\,E_h\]が求まります。エネルギーが低い方の\(\varepsilon_1 \approx -0.65531\,E_h\)を電子が占有するので、ローターン方程式に戻すと、\[ \begin{pmatrix} 0.885 & -0.38980 \\ -0.38980 & 0.17168\end{pmatrix}\begin{pmatrix} C_1^{(1)}\\ C_2^{(1)} \end{pmatrix} =0\]

新しい係数を求める

\(C_1^{(1)}とC_2^{(1)}\)の比は、\[\frac{C_2^{(1)}}{C_1^{(1)}}= 2.2705.\]
規格化条件に代入すると、\[\left(C_1^{(1)}\right)^2+2(0.5952)C_1^{(1)}(2.2705)C_1^{(1)}+\left(2.2705 C_1^{(1)}\right)^2=1\]\[C_1^{(1)} = \pm 0.33600.\]
正の値を取って、\[\left\{C_1^{(1)},C_2^{(1)}\right\} = \left\{0.33600,0.76289\right\} .\]

初期値として\(\left\{C_1^{(0)},C_2^{(0)}\right\} = \left\{1,0\right\}\)を入力したのに、新しい係数として\(\left\{C_1^{(1)},C_2^{(1)}\right\} = \left\{0.33600,0.76289\right\}\)が得られました。これを\(C_1,C_2\)の値が変化しなくなるまで繰替えします。

SCF計算の続き

続きはPythonに任せました。

iteration\(C_1\)\(C_2\)\(\varepsilon_1\) / \(E_h\)\(E\) / \(E_h\)max\(|\Delta C|\)
10.3359970.762871-0.655308-2.7764307.628708e-01
20.4851760.632086-0.970905-2.8325551.491784e-01
30.4515760.663054-0.888560-2.8355133.359954e-02
40.4593840.655939-0.907073-2.8356717.808075e-03
50.4575800.657587-0.902763-2.8356791.804063e-03
60.4579970.657206-0.903759-2.8356804.174173e-04
70.4579010.657294-0.903528-2.8356809.654941e-05
80.4579230.657274-0.903582-2.8356802.233373e-05
90.4579180.657279-0.903569-2.8356805.166129e-06
100.4579190.657278-0.903572-2.8356801.195009e-06
110.4579190.657278-0.903571-2.8356802.764245e-07
120.4579190.657278-0.903572-2.8356806.394137e-08
130.4579190.657278-0.903572-2.8356801.479066e-08
140.4579190.657278-0.903572-2.8356803.421314e-09
150.4579190.657278-0.903572-2.8356807.914042e-10
160.4579190.657278-0.903572-2.8356801.830645e-10
170.4579190.657278-0.903572-2.8356804.234574e-11

最終的に収束した値は、\[C_1 = 0.457919\]\[C_2 = 0.657278\]\[\varepsilon = -0.903572\, E_h\]\[E = -2.835680\, E_h\]でした。したがって、ヘリウム原子の一電子波動関数は、\[\phi_{\mathrm{3-21G}} = (0.457919)\left[ 0.885750e^{-13.6267r^2} + 1.070687e^{-1.99935r^2} \right] + (0.657278)\left[ 0.346979e^{-0.382993r^2} \right]\]

正準仮想軌道を求める
収束した係数から\(F = h + 2 J -K\)を作る

収束した係数を使って正準フォック行列\(\mathbf F\)を作ります。占有軌道の係数行列\(\begin{pmatrix}C_{11}\\C_{21}\end{pmatrix}\)に作用させると、\[\mathbf{J}\begin{pmatrix}C_{11}\\C_{21}\end{pmatrix} = \mathbf{K}\begin{pmatrix}C_{11}\\C_{21}\end{pmatrix}\]になって、打ち消されるので簡略化したフォック行列を使っていました。しかし、非占有軌道では、 \[\mathbf{J}\begin{pmatrix}C_{12}\\C_{22}\end{pmatrix} \neq \mathbf{K}\begin{pmatrix}C_{12}\\C_{22}\end{pmatrix}\] になるので、\(K\)の数値も計算して入れた正準フォック行列\(\mathbf F\)を作ります。

収束係数は、\[\begin{pmatrix}C_{11}\\C_{21}\end{pmatrix} = \begin{pmatrix}0.4579188741\\0.6572778604 \end{pmatrix} , \]重なり行列は\[S = \begin{pmatrix}1 & 0.5952159497 \\0.5952159497 & 1\end{pmatrix}.\]

正準HFの行列要素\[F_{11} = h_{11} + 2J_{11} – K_{11}\]を計算します。

\[h_{11} = \int \left(0.885750e^{-13.6267r^2} + 1.070687e^{-1.99935r^2}\right)\left[-\frac{1}{2}\nabla_1^2 – \frac{2}{r_1}\right]\left(0.885750e^{-13.6267r^2} + 1.070687e^{-1.99935r^2}\right)\,d\tau_1 \]\[ = -1.5735646753.\]

基底関数\[\chi_1 = 0.885750e^{-13.6267r^2} + 1.070687e^{-1.99935r^2} \]\[ \chi_2 = 0.346979e^{-0.382993r^2}\]同士の電子間反発積分をはじめに計算しておきます。

電子間反発積分を \[ (ij|kl)= \iint \chi_i(\mathbf r_1)\chi_j(\mathbf r_1) \frac{1}{r_{12}} \chi_k(\mathbf r_2)\chi_l(\mathbf r_2) \,d\tau_1d\tau_2\]と定義すると、\[(11|11)=1.8032545920,\]\[ (11|12)=0.8765068190,\]\[ (11|22)=0.9170082837,\]\[ (12|12)=0.4520772565. \]

クーロン項\(J_{11}\)は、\[\iint \chi_1(\mathbf r_1)^2\frac{1}{r_{12}}\phi(\mathbf r_2)^2\,d\tau_1 d\tau_2.\]

ここに \(\phi=C_{11}\chi_1+C_{21}\chi_2\) を代入して展開すると、\[J_{11} = C_{11}^2(11|11) +2C_{11}C_{21}(11|12) +C_{21}^2(11|22)\]です。数値を入れると、\[ J_{11} =(0.2096896952)(1.8032545920)+(0.6019598756)(0.8765068190)+(0.4320141858)(0.9170082837)\]\[=0.3781239058+0.5276219357+0.3961605870\]\[=1.3019064285\]

交換項\(K_{11}\)は、\[\iint \chi_1(\mathbf r_1)\phi(\mathbf r_1)\frac{1}{r_{12}}\phi(\mathbf r_2)\chi_1(\mathbf r_2)\,d\tau_1 d\tau_2.\]
 \(\phi=C_{11}\chi_1+C_{21}\chi_2\) を代入して展開すると、\[K_{11} =C_{11}^2(11|11) +C_{11}C_{21}(11|12) +C_{11}C_{21}(12|11) +C_{21}^2(12|12). \]\[ = C_{11}^2(11|11) +2C_{11}C_{21}(11|12) +C_{21}^2(12|12). \]
数値を代入すると、\[ K_{11} =(0.2096896952)(1.8032545920)+(0.6019598756)(0.8765068190)+(0.4320141858)(0.4520772565)\]\[ =0.3781239058+0.5276219357+0.1953037879\]\[ =1.1010496294. \]

これらを足し合わせると、\[F_{11} =h_{11}+2J_{11}-K_{11}=-1.5735646753 +2(1.3019064285) -1.1010496294\]\[ =-1.5735646753 +2.6038128570 -1.1010496294\]\[=-0.0708014476\ E_h.\]

残りも全部同様に計算すると、\[\mathbf{F} = \begin{pmatrix}-0.0708014476&-1.1180026983\\-1.1180026983&-0.4993641791\end{pmatrix}.\]

\[\begin{pmatrix}C_{11} & C_{12}\\ C_{21} & C_{22} \end{pmatrix}^{T}\mathbf{F}\begin{pmatrix}C_{11} & C_{12}\\ C_{21} &C_{22} \end{pmatrix}= \begin{pmatrix}C_{11} & C_{12} \\ C_{21} & C_{22} \end{pmatrix}^{T}\mathbf{S}\begin{pmatrix}C_{11} & C_{12}\\ C_{21} &C_{22} \end{pmatrix}\begin{pmatrix}\varepsilon_1 & 0 \\ 0 & \varepsilon_2\end{pmatrix} = \begin{pmatrix}\varepsilon_1 & 0 \\ 0 & \varepsilon_2\end{pmatrix} \]なので、\[\begin{pmatrix}0.4579188741 & 0.6572778604\\ C_{12} & C_{22} \end{pmatrix}
\begin{pmatrix}-0.0708014476&-1.1180026983\\-1.1180026983&-0.4993641791\end{pmatrix}\begin{pmatrix}0.4579188741 & C_{12}\\ 0.6572778604 &C_{22} \end{pmatrix}=\begin{pmatrix}\varepsilon_1 & 0 \\ 0 & \varepsilon_2\end{pmatrix} .\]
行列の掛け算をすると、\[ \begin{pmatrix}-0.9035715084 &-0.76725974\,C_{12} -0.84017556\,C_{22}  \\ -0.76725974\,C_{12} -0.84017556\,C_{22}  &-0.0708014476\,C_{12}^2 +2(-1.1180026983)\,C_{12}C_{22} +(-0.49936417)C_{22}^2 \end{pmatrix}= \begin{pmatrix}\varepsilon_1 & 0 \\ 0 & \varepsilon_2\end{pmatrix} .\]

非対角成分を0と置いた式\[-0.76725974\,C_{12} -0.84017556\,C_{22}  = 0\]から、\[ \frac{C_{22}}{C_{12}} = -\frac{0.76725974}{0.84017556}= -0.9132136.\]

規格化条件、\[C_{12}^2 + 2(0.5952159497)C_{12}C_{22}+ C_{22}^2 =1 \]
を解いて、\(C_{12} > 0\)の解を選ぶと、\[C_{12} = 1.1571404526,\]\[C_{22} = -1.0567163936 .\]

最後に\(\epsilon_2\)を計算すると、\[\varepsilon_2 = -0.0708014476C_{12}^2 + 2(-1.1180026983)C_{12}C_{22} + (-0.4993641791)C_{22}^2\]\[= 2.0817026437\ E_h\]

3-21G基底関数系を使った時のヘリウム原子のHF正準仮想軌道の波動関数は、\[\phi_{virtual} = (1.1571404526)\left[ 0.885750e^{-13.6267r^2} + 1.070687e^{-1.99935r^2} \right] + (-1.0567163936)\left[ 0.346979e^{-0.382993r^2} \right]\]です。図示すると以下の通りです。

この仮想軌道は占有軌道と直交する必要があるため、2本の正値のs型基底関数を反対符号で組み合わせることになり、1つの動径節（\(0.86796613\,a_0 \)）が生じます。そのため水素様2s軌道に似た形になりますが、ヘリウムの2s原子軌道そのものではありません。

(ここまで)

コメントを残す

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です

CAPTCHA


このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください。