Hartree-Fockの数値解法：Roothaan方程式の導出

---

1. Hartree-Fock方程式

Hartree-Fock方程式とは、\(N\)電子系の波動関数を近似的に求める方程式のことであり、それぞれの電子が独立の波動関数に従うという近似を使う。具体的には、軌道\(\{\phi_i\}\)のスレーター行列式 \[ \Phi_{\phi_1\phi_2\cdots\phi_N}(\b{r}_1,\b{r}_2,\cdots,\b{r}_N) = \frac{1}{\sqrt{N!}}\left|\begin{array}{cccc} \phi_1(\b{r}_1) & \phi_1(\b{r}_2) & \cdots & \phi_1(\b{r}_N) \\ \phi_2(\b{r}_1) & \phi_2(\b{r}_2) & \cdots & \phi_2(\b{r}_N) \\ \vdots & \vdots & \ddots & \vdots \\ \phi_N(\b{r}_1) & \phi_N(\b{r}_2) & \cdots & \phi_N(\b{r}_N) \\ \end{array} \right|\tag{1} \] を解の波動関数として仮定し、これがシュレディンガー方程式の真の解\(\psi\)にできるだけ近くなるような関数の組\(\{\phi_i\}\)を見つけ出すのが、Hartree-Fock 方程式である。このページでは、特に断らない限り\(\b{r}\)は空間座標とスピン座標をひっくるめた座標とする。方程式の形は以下のようであった。電荷\(Z_ne\)を持つ原子核が位置\(\b{R}_n\)に存在しているとし、1電子ハミルトニアンを \[\hat{H}^{(1)}(\b{r}) = -\frac{\hbar^2}{2m}\nabla^2 - \sum_n \frac{Z_n e^2}{4\pi\epsilon_0|\b{r}-\b{R}_n|} \tag{2}\] と書く。また、2電子ハミルトニアン (クーロン相互作用) を \[\hat{H}^{(2)}(\b{r},\b{r}') = \frac{e^2}{4\pi\epsilon_0|\b{r}-\b{r}'|} \tag{2}\] とする。Hartree-Fock近似の下で、古典的なクーロン相互作用ポテンシャルを表す\(U_C(\b{r})\)と、量子力学的なクーロン相互作用ポテンシャル (交換相互作用) を表すカーネル関数 \(V_E(\b{r},\b{r}')\)を \begin{align} U_C(\b{r}) &= \sum_{i=1}^N\int \phi_{i}^*(\b{r}')\hat{H}^{(2)}(\b{r},\b{r}')\phi_{i}(\b{r}') d\b{r}'\tag{3}\\ V_E(\b{r},\b{r}) &= \sum_{i=1}^N\phi_{i}^*(\b{r}')\hat{H}^{(2)}(\b{r},\b{r}')\phi_{i}(\b{r})\tag{4} \end{align} と定義すれば、Hartree-Fock方程式は、 \begin{align} \left(\hat{H}^{(1)}(\b{r}) + U_C(\b{r})\right) \phi_{k}(\b{r})-\int V_E(\b{r},\b{r}')\phi_{k}(\b{r}')d\b{r}' = E_k \phi_k(\b{r}) \tag{5} \end{align} である。これを解けば確かに1電子波動関数\(\{\phi_i\}\)が求まるのだが、これは\(k=1,2,\cdots,N\)のそれぞれについての方程式であり、「\(N\)本の非線形・偏微分・積分連立方程式」というとんでもない代物である。もちろん手では解けないので、今回は、これを数値的に解くときに利用されるRoothaan方程式を導出する。

---

---

2.基底関数を用意する

Hartree-Fock方程式を直接、ゼロから解くのはいくらなんでもしんどい。そこで\(\phi_i\)を適当な\(M\)個の基底関数\(f_j\)によって展開し、 \[\phi_i(\b{r}) = \sum_{j=1}^M c_{ji}f_j(\b{r})\tag{6}\] として、この係数\(c_{ji}\)を求める問題に帰着させてしまおう、というのがRoothaan方程式の考え方である。もちろん、基底関数\(f_j\)によって展開された\(\{\phi_i\}\)が、真のHartree-Fock方程式の解になっているとは限らない。どうしたって選んだ基底関数\(f_j\)だけでは手の届かない関数たちも存在するからだ。(無限個の基底関数を使うなら原理的には可能だが、だからといって計算機は無限個の関数は扱えない。) しかしながら、Hartree-Fock方程式自体がもともと変分法的に導出された方程式だったことから、ある仮定 (この場合基底関数の和として表せるという仮定) をしたとしても、Hartree-Fock方程式が近似的にでも成り立つ解が見つかれば、その解は真の解の良い近似になると考えられる。

なるべく少ない基底関数で良い近似を与えそうな\(f_i\)の選び方としては、例えば孤立原子の1s軌道や2p軌道なんかを使う方法があるだろう。\(f_i\)としてそのようなものを採用するとき、それは LCAO (Linear Combination of Atomic Orbitals) 法と呼ばれる。

---

3.係数の方程式に書き換える

Hartree-Fock 方程式 (5) に、(6) を代入して係数の方程式に書き換えよう。

まずクーロン相互作用項(3)は \begin{align} &U_C(\b{r}) \\ &= \sum_{i=1}^N\int \phi_{i}^*(\b{r}')\hat{H}^{(2)}(\b{r},\b{r}')\phi_{i}(\b{r}') d\b{r}'\\ &= \sum_{i=1}^N\int \left(\sum_{j=1}^M c_{ji}^*f_j^*(\b{r}')\right)\hat{H}^{(2)}(\b{r},\b{r}')\left(\sum_{k=1}^M c_{ki}f_k(\b{r}')\right)d\b{r}'\\ &= \sum_{i=1}^N\sum_{j=1}^M \sum_{k=1}^M c_{ji}^*c_{ki}\int f_j^*(\b{r}')\hat{H}^{(2)}(\b{r},\b{r}')f_k(\b{r}')d\b{r}' \end{align} となる。交換相互作用のカーネル (4) は、 \begin{align} &V_E(\b{r}) \\ &= \sum_{i=1}^N\sum_{j=1}^M \sum_{k=1}^M c_{ji}^*c_{ki}\left(f_j^*(\b{r}')\hat{H}^{(2)}(\b{r},\b{r}')f_k(\b{r})\right) \end{align} Hartree-Fock方程式 (5) に全て代入してしまうと、 \begin{align} &\left(\hat{H}^{(1)}(\b{r}) + U_C(\b{r})\right) \phi_{l}(\b{r})-\int V_E(\b{r},\b{r}')\phi_{l}(\b{r}')d\b{r}'\\ &~~~~~~~~~= E_l \phi_l(\b{r})\\\\ &\left(\hat{H}^{(1)}(\b{r}) + \sum_{i=1}^N\sum_{j=1}^M \sum_{k=1}^M c_{ji}^*c_{ki}\int f_j^*(\b{r}')\hat{H}^{(2)}(\b{r},\b{r}')f_k(\b{r}')d\b{r}'\right) \left(\sum_{m=1}^M c_{ml}f_m(\b{r})\right)\\ &~~~~~~~~~-\int \left(\sum_{i=1}^N\sum_{j=1}^M \sum_{k=1}^M c_{ji}^*c_{ki}\left(f_j^*(\b{r}')\hat{H}^{(2)}(\b{r},\b{r}')f_k(\b{r})\right)\right)\left(\sum_{m=1}^M c_{ml}f_m(\b{r}')\right)d\b{r}'= E_l \sum_{m=1}^M c_{ml}f_m(\b{r})\\\\ &\sum_{m=1}^Mc_{ml} \left(\hat{H}^{(1)}(\b{r})f_m(\b{r}) + \sum_{i=1}^N\sum_{j=1}^M \sum_{k=1}^M c_{ji}^*c_{ki}\int f_j^*(\b{r}')\hat{H}^{(2)}(\b{r},\b{r}')f_k(\b{r}')f_m(\b{r})d\b{r}'\right)\\ &~~~~~~~~~-\sum_{i=1}^N\sum_{j=1}^M \sum_{k=1}^M\sum_{m=1}^M c_{ji}^*c_{ki}c_{ml}\int f_j^*(\b{r}')\hat{H}^{(2)}(\b{r},\b{r}')f_m(\b{r}')f_k(\b{r})d\b{r}'= E_l \sum_{m=1}^M c_{ml}f_m(\b{r})\\ \end{align} さて、次に\(f_n^*(\b{r})\)をかけた後\(\b{r}\)で積分して、本当に係数だけの方程式にしてしまおう。そろそろ式が長くなりすぎたので、次のように記号を定義する。 \begin{align} S_{nm} &= \int f_n^*(\b{r})f_m(\b{r})d\b{r}\tag{7}\\ H_{nm} &= \int f_n^*(\b{r})\hat{H}^{(1)}(\b{r})f_m(\b{r})d\b{r}\tag{8}\\ U_{njmk} &= \int f_n^*(\b{r})f_j^*(\b{r}')\hat{H}^{(2)}(\b{r},\b{r}')f_m(\b{r})f_k(\b{r}')d\b{r}d\b{r}'\tag{9} \end{align} \(S_{nm}\)は重なり積分と呼ばれる。実際には、この積分は計算に使う基底関数を決めた時点で、数値的に (or 公式を使って) 評価してしまえば良い。 そうすると、Hartree-Fock方程式は \begin{align} &\sum_{m=1}^Mc_{ml} \left(H_{nm} + \sum_{i=1}^N\sum_{j=1}^M \sum_{k=1}^M c_{ji}^*c_{ki}U_{njkm}\right)\\ &~~~~~~~~~-\sum_{i=1}^N\sum_{j=1}^M \sum_{k=1}^M\sum_{m=1}^M c_{ji}^*c_{ki}c_{ml}U_{njmk} = E_k \sum_{m=1}^M c_{ml}S_{nm}\\\\ &\sum_{m=1}^Mc_{ml} \left(H_{nm} + \sum_{i=1}^N\sum_{j=1}^M \sum_{k=1}^M c_{ji}^*c_{ki}(U_{njkm}- U_{njmk})\right) = E_l \sum_{m=1}^M c_{ml}S_{nm} \end{align} となる。さらにきれいにするために、行列\(C = \{c_{ji}\}\)を定義しよう。まず、 \[\sum_{i=1}^N c_{ji}^*c_{ki} = (CC^\dagger)_{kj}\] である。さらに、\(C\)の関数としての行列 \[F_{nm}(C) = H_{nm} + \sum_{j=1}^M \sum_{k=1}^M (CC^\dagger)_{kj}(U_{njkm}- U_{njmk})\] \(F=\{F_{nm}\}\)を定義する。また、重なり積分も行列にしてしまって、\(S=\{S_{nm}\}\)、エネルギー\(E_l\)も\(E=\{\delta_{lj}E_j\}\)という対角行列にしてしまう。そうすると、Hartree-Fock方程式を係数行列\(C\)についての方程式に書き換えた、次の

Roothaan方程式

を得る。 \begin{align} F(C)C &= SCE \end{align}

---

4.まとめ

文字がたくさん出てきてわけが分からなくなりそうなので、最後にまとめよう。Hartree-Fock方程式を解くために、 \[\phi_i(\b{r}) = \sum_{j=1}^M c_{ji}f_j(\b{r})\tag{6}\] とおいて、方程式を書き換えた。結果は、 \begin{align} F(C)C &= SCE ~~~ (もしくは~F(C)\b{c}_l =E_lS\b{c}_l) \\ F_{nm}(C) &= H_{nm} + \sum_{j=1}^M \sum_{k=1}^M (CC^\dagger)_{kj}(U_{njkm}- U_{njmk})\\ S_{nm} &= \int f_n^*(\b{r})f_m(\b{r})d\b{r}\tag{7}\\ H_{nm} &= \int f_n^*(\b{r})\hat{H}^{(1)}(\b{r})f_m(\b{r})d\b{r}\tag{8}\\ U_{njmk} &= \int f_n^*(\b{r})f_j^*(\b{r}')\hat{H}^{(2)}(\b{r},\b{r}')f_m(\b{r})f_k(\b{r}')d\b{r}d\b{r}'\tag{9}\\ E_{lj} &= \delta_{lj}E_j \end{align} となった。\(N\)は電子数、\(M\)は基底関数の数である。\(C\)は\(M\times N\)行列、\(F(C),S\)は\(M\times M\)行列、\(E\)は\(N\times N\)行列だ。つまり、Roothaan方程式は\(M\times N\)の行列に関する方程式である。そして結局求めたかったのは、行列\(C,E\)の値である。

この方程式を実際に解く場合 のような手続きを踏む。\(F(C)\)が\(M\times M\)の行列なので、実際には固有ベクトルとして\(M\)個のベクトルが求まってしまうが、\(C\)にするときには、エネルギーの低い方から\(N\)個取ってやることで辻褄合わせをする。\(M-N\)の余った分は

仮想軌道 (virtual orbit)

と呼ばれる。

一般化固有値問題を解く際には、\(S\)は先に対角化しておいたほうが簡単だろう。つまり、 \[S = UDU^\dagger\] と対角化しておき、一般化固有値問題\(F(C)\b{c}=E_lS\b{c}\)を、通常の固有値問題 \[F'(C)\b{c}_l'=E_l\b{c}_l'\] \[\b{c}_l = D^{-1}\b{c}'_l\] \[F'(C) = U^\dagger F(C)U \] に帰着させてから解くわけだ。

一般化固有値問題を解く本当のところはどうやっているかしらない。こんなふうかな？という予想だ。 実際には、きっともっとかっこいいアルゴリズムが使われていることだろう。

ちなみに、スピン軌道がハミルトニアンに現れず、スピンを独立に考える場合は、 \[F_{nm}(C) = H_{nm} + \sum_{j=1}^M \sum_{k=1}^M (CC^\dagger)_{kj}(2U_{njkm}- U_{njmk}) \] となる。