哈特里-福克方程

哈特里-福克方程（Hartree–Fock equation），又稱為HF方程，是一個應用變分法計算多電子系統波函數的方程式，是量子物理、凝聚態物理學、量子化學中最重要的方程之一。HF方程形式上是單電子本徵方程，求得的本徵態是單電子波函數，即分子軌道。以HF方程為核心的數值計算方法稱為「哈特里-福克方法」（Hartree–Fock method）。

基於分子軌道理論的所有量子化學計算方法都是以HF方法為基礎的。鑑於分子軌道理論在現代量子化學中的廣泛應用，HF方程被視為現代量子化學的基石。

## 哈特里-福克近似

哈特里-福克近似也稱為分子軌道近似或單行列式近似，認為多電子體系波函數可以由體系分子軌道波函數構造的單個斯萊特行列式表示：

$$\Psi=\frac{1}{\sqrt{N!}}\begin{vmatrix} \chi_1(1)&\chi_2(1) & \cdots &\chi_N(1) \\ \chi_1(2)&\chi_2(2) & \cdots &\chi_N(2) \\ \vdots&\vdots&\ddots&\vdots\\ \chi_1(N)&\chi_2(N) & \cdots &\chi_N(N) \\ \end{vmatrix}=|\chi_1\chi_2\cdots\chi_N\rangle$$

這裡$N$為電子數，$\chi_i(x)$是能量最低的$N$個分子軌道波函數。而分子軌道波函數是由體系中所有原子軌道波函數經過線性組合構成的。

$$\chi_a(x)=\sum_i C_{ai}\phi_i(x)$$

通常原子軌道波函數$\phi_i(x)$由計算時指定的基組給出。

在此近似下，由變分法，體系的多電子薛丁格方程可化為單電子方程

$$f_{(1)}\chi_{a}(1)=\varepsilon_{a}\chi_{a}(1)$$

其中$f_{(1)}$為福克算符，$\chi_{a(1)}$為體系的分子軌道，$\varepsilon_{a}$是分子軌道$\chi_{a(1)}$的軌道能。這就是哈特里-福克方程。求解該方程，即求解使方程成立的一組分子軌道；這組分子軌道使得在HF近似下的體系總能量最低。

## 發展簡史

1927年，物理學家瓦爾特·海特勒和弗里茨·倫敦完成了氫氣分子的量子力學計算之後，開啟了量子化學的時代。從此，人們便開始嘗試使用量子力學理論來解釋化學物質結構和化學現象。

為了解決多電子體系薛丁格方程近似求解的問題，量子化學家道格拉斯·哈特里在1928年提出了哈特里假設，將每個電子看做是在其他所有電子構成的平均勢場中運動的粒子，並且首先提出了迭代法的思路。哈特里根據他的假設，將體系電子哈密頓算子分解為若干個單電子哈密頓算子的簡單代數和，每個單電子哈密頓算子中只包含一個電子的坐標，因而體系多電子波函數可以表示為單電子波函數的簡單乘積，這就是哈特里方程。但是由於哈特里沒有考慮電子波函數的反對稱要求，他的哈特里方程實際上是非常不成功的。

1930年，哈特里的學生弗拉基米爾·福克和約翰·斯萊特分別提出了考慮泡利原理的自洽場迭代方程和單行列式型多電子體系波函數，這就是今天的哈特里-福克方程。但是由於計算上的困難，HF方程誕生後整整沉寂了二十年。在1950年，量子化學家克萊門斯·羅特漢想到將分子軌道用原子軌道的線性組合來近似展開，而得到了閉殼層結構的羅特漢方程。

1953年，美國的R.帕里瑟、R.帕爾和英國的約翰·波普爾花費兩年時間使用手搖計算器分別獨立地實現了對氮氣分子的RHF自洽場計算，這是人類首次通過求解HF方程獲得對化學結構的量子力學解釋，也是量子化學計算方法第一次實際完成。在第一次成功之後，伴隨著電腦技術的迅猛發展，HF方程與量子化學一道獲得長足發展，在HF方程的基礎上，人們發展出了高級量子化學計算方法，使得計算精度進一步提高，通過對HF方程電子積分的簡化和參數化，人們大大縮減了量子化學的計算量，使得對超過1000個原子的中等大小分子的計算成為可能。

## 方程的導出

哈特里-福克方程源於對多電子體系電子波函數的變分法處理。在玻恩-奧本海默近似條件下，一個多電子體系的電子運動與能量可以與原子核的運動和能量相互分離，這樣利用電子哈密頓算符和多電子波函數便可以計算體系的電子能量。其能量的表達式為:

$$E_0=\langle\Psi_0|H_{ele}|\Psi_0\rangle$$

式中$E_0$表示體系基態電子能量；$H_{ele}$表示體系的電子哈密頓算符。

$$H_{ele}=\sum_i^N \frac{-1}{2}\nabla_i^2 - \sum_{i=1}^N\sum_{a=1}^M \frac{Z_{a}}{R_{ia}}+\sum_{i=1}^N\sum_{j>i}^N \frac{1}{r_{ij}}$$

根據作用方式，可以將$H_{ele}$分解為兩部分

$$H_{ele}=O_1+O_2$$

其中$O_1$為單電子算符

$$O_1=\sum_i^N \frac{-1}{2}\nabla_i^2 - \sum_{i=1}^N\sum_{a=1}^M \frac{Z_{a}}{R_{ia}}=\sum_i^N -\left(\frac{1}{2}\nabla_i^2 +\sum_{a=1}^M \frac{Z_{a}}{R_{ia}} \right)=\sum_i^N h_i$$

描述單個電子的動能和原子核吸引勢能；而$O_2$為二電子算符

$$O_2=\sum_{i=1}^N\sum_{j>i}^N \frac{1}{r_{ij}}$$

描述電子間相互作用。

$\Psi_0$代表基態多電子波函數，是體系單電子分子軌道波函數為基函數組建的斯萊特行列式。構建$\Psi_0$的各個分子軌道相互之間是正交歸一的，即約束條件為

$$\langle\chi_a|\chi_b\rangle=\delta_{ab}$$

考慮此約束條件，應用拉格朗日乘數法對函數$L=E_0[\{\chi_i\}]-\sum_{a,b}^N \varepsilon_{ab}\left(\langle a|b\rangle-\delta_{ab}\right)$變分求極值。式中$\varepsilon_{ab}$是拉格朗日待定因子，$\langle a|b\rangle$是$\langle\chi_a|\chi_b\rangle$的簡寫。

變分法的處理過程如下：令

$$\delta L=\delta E_0[\{\chi_i\}]-\sum_{a,b}^N\varepsilon_{ab}\delta[\chi_a|\chi_b]=0$$

其中

$$\begin{align}\delta E_0=&\sum_a^N [\delta \chi_a|h|\chi_a]+\sum_a^N [\chi_a|h|\delta \chi_a]\\ &+\frac{1}{2}\sum_{a,b}^N \left([\delta\chi_a\chi_a|\chi_b\chi_b]+[\chi_a\delta\chi_a|\chi_b\chi_b]+[\chi_a\chi_a|\delta\chi_b\chi_b]+[\chi_a\chi_a|\chi_b\delta\chi_b]\right)\\ &+\frac{1}{2}\sum_{a,b}^N \left([\delta\chi_a\chi_b|\chi_b\chi_a]+[\chi_a\delta\chi_b|\chi_b\chi_a]+[\chi_a\chi_b|\delta\chi_b\chi_a]+[\chi_a\chi_b|\chi_b\delta\chi_a]\right)\end{align}$$

這裡定義單電子積分記號

$$[\chi_a|\chi_b]=\langle\chi_a|\chi_b\rangle=\int dx\chi_a^*(x)\chi_b(x)$$

以及二電子積分記號

$$[\chi_a\chi_b|\chi_c\chi_d]=\int dx_1dx_2\chi_a^*(x_1)\chi_b(x_1)r_{12}^{-1}\chi_c^*(x_2)\chi_d(x_2)$$

考慮到流動坐標的不可區分性，可以簡化為：

$$\begin{align}\delta E_0=&\sum_a^N [\delta \chi_a|h|\chi_a]+\sum_{a,b}^N \left([\delta\chi_a\chi_a|\chi_b\chi_b]-[\delta\chi_a\chi_b|\chi_b\chi_a]\right)\\ &+\left(\sum_a^N [\delta \chi_a|h|\chi_a]+\sum_{a,b}^N \left([\delta\chi_a\chi_a|\chi_b\chi_b]-[\delta\chi_a\chi_b|\chi_b\chi_a]\right)\right)^*\end{align}$$

同理，$\delta L$中的$\sum_{a,b}^N\varepsilon_{ba}\delta[a|b]$項有：

$$\sum_{a,b}^N\varepsilon_{ba}\delta[a|b]=\sum_{a,b}^N\varepsilon_{ba}[\delta\chi_a|\chi_b]+\left(\sum_{a,b}^N\varepsilon_{ba}[\delta\chi_a|\chi_b]\right)^*$$

將兩項相加，$\delta L$表示為：

$$\begin{align}\delta L=&\sum_a^N [\delta \chi_a|h|\chi_a]+\sum_{a,b}^N \left([\delta\chi_a\chi_a|\chi_b\chi_b]-[\delta\chi_a\chi_b|\chi_b\chi_a]\right)+\sum_{a,b}^N\varepsilon_{ba}[\delta\chi_a|\chi_b]\\ &+\left(\sum_a^N [\delta \chi_a|h|\chi_a]+\sum_{a,b}^N \left([\delta\chi_a\chi_a|\chi_b\chi_b]-[\delta\chi_a\chi_b|\chi_b\chi_a]\right)+\sum_{a,b}^N\varepsilon_{ba}[\delta\chi_a|\chi_b]\right)^*\\ =&0\end{align}$$

若$L$函數處於極值點，則變量$\left\{\chi_i\right\}$向各個方向的微小變化都應該有$\delta L=0$。可以取$\left\{\delta \chi_i\right\}=\left\{i\delta \chi_i\right\}$沿虛軸變分，則在$\delta L$表達式中，第一項前會產生一個$i$的係數，第二項復共軛會產生一個$-i$係數：

$$\begin{align}\delta L=&i\sum_a^N [\delta \chi_a|h|\chi_a]+\sum_{a,b}^N \left([\delta\chi_a\chi_a|\chi_b\chi_b]-[\delta\chi_a\chi_b|\chi_b\chi_a]\right)+\sum_{a,b}^N\varepsilon_{ba}[\delta\chi_a|\chi_b]\\ &-i\left(\sum_a^N [\delta \chi_a|h|\chi_a]+\sum_{a,b}^N \left([\delta\chi_a\chi_a|\chi_b\chi_b]-[\delta\chi_a\chi_b|\chi_b\chi_a]\right)+\sum_{a,b}^N\varepsilon_{ba}[\delta\chi_a|\chi_b]\right)^*\\ =&0\end{align}$$

消去虛數單位，兩式相加，可以消去表達式中的復共軛項：

$$\delta L=\sum_a^N [\delta \chi_a|h|\chi_a]+\sum_{a,b}^N \left([\delta\chi_a\chi_a|\chi_b\chi_b]-[\delta\chi_a\chi_b|\chi_b\chi_a]\right)+\sum_{a,b}^N\varepsilon_{ba}[\delta\chi_a|\chi_b]=0$$

在引入庫侖算符和交換算符的概念之後，上述表達式可以改寫為：

$$\sum_a^N \int dx_1 \delta \chi_{a}^*(1)\left[h_1\chi_{a}(1)+\sum_b^N\left(J_{b(1)}-K_{b(1)}\right)\chi_{a}(1)-\sum_b^n\varepsilon_{ba}\chi_{b}(1)\right]=0$$

對任意$\left\{\delta\chi_i\right\}$上述等式均應成立，因而必須有:

$$h_1\chi_{a}(1)+\sum_b^N\left(J_{b(1)}-K_{b(1)}\right)\chi_{a}(1)-\sum_b^n\varepsilon_{ba}\chi_{b}(1)=0$$

整理得到：

$$\left[h_{(1)}+\sum_b^N\left(J_{b(1)}-K_{b(1)}\right)\right]\chi_{a}(1)=\sum_b^N\varepsilon_{ba}\chi_{b}(1)$$

定義福克算符

$$f_{(1)}=h_{(1)}+\sum_b^N\left(J_{b(1)}-K_{b(1)}\right)$$

方程可以表達為

$$f_{(1)}\chi_{a}(1)=\sum_b^N\varepsilon_{ba}\chi_{b}(1)$$

即哈特里-福克方程。為了求解，通過對分子軌道波函數進行酉變換處理，使得軌道能矩陣$\varepsilon_{ba}$對角化，將一般的不可解的HF方程轉化為正則哈特里-福克方程：

$$f_{(1)}\chi_{a}(1)=\epsilon_{a}\chi_{a}(1)$$

此方程形式上為本徵方程，但是福克算符中的庫侖算符和交換算符都與分子軌道有關，因此只能夠通過自洽迭代的方法近似求解，即哈特里-福克自洽場（HF-SCF）方法。HF-SCF方法是組態相互作用方法、多體微擾理論、半經驗量子化學計算等現代量子化學計算方法的基礎。

## 方程的解與自洽場方法

正則哈特里-福克方程雖然具有簡單的本徵方程形式，但福克算子中的庫侖算子和交換算子中含有所有$\chi_{a}$的表達式，因而方程的實際形式非常複雜，無法求得精確的解析解，只能使用迭代法求解，即量子化學中的自洽場方法。

在實際操作中，人們會首先給定基組，將本徵方程轉化為矩陣方程。這一轉變使分子軌道可以表示為$\chi_a(x)=\sum_i C_{ai}\phi_i(x)$，可以用一個變化的$K$維向量$\mathbb{C_a}$來代表分子軌道，$K$為基組中基函數的數量。同理，可以將福克算符轉化為福克矩陣。然而由於基函數不一定正交，存在重疊積分

$$S_{ij}=\int dx\phi_i^*(x)\phi_j(x)$$

最終HF方程的形式轉化為廣義本徵方程

$$\mathbb{F}\mathbb{C}=\mathbb{\varepsilon}\mathbb{S}\mathbb{C}$$

此方程稱為哈特里-福克-羅特漢方程，或羅特漢方程。這樣，相對複雜的本徵方程就轉化為只需要進行簡單代數計算就可以求解的矩陣本徵方程，而原方程中複雜的積分則在上述轉化過程中一次性完成了。

求解時首先需要給定一個分子軌道的初始猜測$\mathbb{C}$，用這套分子軌道計算福克矩陣$\mathbb{F}$，再求解哈特里-福克-羅特漢方程，得到一組新的分子軌道$\mathbb{C}$和相應的本徵值矩陣$\mathbb{\varepsilon}$，再使用新軌道$\mathbb{C}$計算新的福克矩陣$\mathbb{F}$......重複上述過程，直到$\mathbb{C}$的各個矩陣元數值不再有明顯的變化，稱作自洽迭代收斂，就得到了HF方程的解。

得到收斂的$\mathbb{C}$矩陣後，將這些係數代入$\chi_a=\sum_i C_{ia}\phi_i$與基函數結合，便獲得了最終的分子軌道波函數以及體系電子總能量等各種性質。

## 方程的應用

HF方程在量子化學中有著廣泛的應用，所有分子軌道理論的量子化學計算都是以HF方程為基礎的。

- 組態相互作用方法（CI）：在CI方法中，通過HF方程解得的一系列分子軌道用於構建多電子基函數集，在構建了多電子基函數集後再通過變分法處理獲得CI能量的最低點，因而進行CI計算必須首先完成HF方程的求解。
- 多體微擾理論方法（MPn）：MPn計算將體系所有福克算符的代數和定義為哈特里-福克哈密頓算符，將電子間相互作用所產生的能量項看作是對哈特里-福克哈密頓算符的微擾，經過微擾處理後可以獲得體系能量的近似值。進行多體微擾計算也需要首先進行HF方程的求解，以獲得需要的分子軌道波函數形式和分子軌道能量。
- 半經驗量子化學計算：半經驗量子化學計算是對HF方程求解過程的簡化。在HF方程的求解中，絕大部分計算量都分布在由正則HF方程向矩陣本徵方程形式轉變的過程中，如果將這一過程中大量的電子積分用經驗數值代替，便可以極大地縮短HF方法的計算時間。為此，針對不同的研究體系，量子化學家開發了不同的積分經驗常數，與之相應地產生了各種半經驗量子化學計算方法。本質上講，半經驗計算仍然是通過自洽場方法求解HF方程的過程。

## 外部連結

- A Numerical Hartree-Fock Program