本文是「第一性原理的微观计算模拟」系列的第 2 篇（共 11 篇），内容整理自同名书稿第 3 章，公式编号与原书一致。文中引用的文献依据原书参考文献清单整理并列于文末，按原书清单顺序从 1 开始编号。\nocite{*}
均匀电子气模型
在本章开头讨论 Hartree-Fock 方法时我们已经提到，相互作用电子气由于多体效应而出现了一些非经典的能量项，比如交换能与关联能，与之相应的电子气分布也不同于独立电子气系统。对于绝大部分体系，都无法求得这些能量项的解析表达式。但是对于极个别情况，如凝胶模型（jellium model）下的均匀电子气，体系的各项能量可以解析地或者比较精确地给出。这无疑有助于我们对交换关联项的理解，因此有必要对其进行详细介绍。在凝胶模型下，电子气在空间中均匀分布，且嵌入同样在空间中均匀分布的正电荷背景。
引入电子密度参数（Wigner–Seitz 半径） $r_{\mathrm{s}}^{0}$，其物理意义为：按密度 $\rho_0$ 均匀分布的电子，平均每个电子所占据的球体的半径。为了讨论方便，我们将 $r_{\mathrm{s}}^{0}$ 写为 $r_{\mathrm{s}}a_0$，其中 $a_0$ 为玻尔半径。显而易见，有关系式
$$ \frac{4\pi}{3}r_{\mathrm{s}}^{3}=\frac{1}{\rho_0a_0^{3}} \tag{3.92} $$
利用 2.4.1.3 小节中的结果，可知在不考虑自旋的情况下，有
$$ \rho_0=\frac{2}{(2\pi)^{3}}\int f(E(\boldsymbol{k}))\,\mathrm{d}\boldsymbol{k}=\frac{1}{\pi^{2}}\int_{0}^{k_{\mathrm{F}}}k^{2}\,\mathrm{d}k=\frac{k_{\mathrm{F}}^{3}}{3\pi^{2}} \tag{3.93} $$
式中：$k_{\mathrm{F}}$ 是费米波矢的大小。由式（3.92）及式（3.93）可得
$$ k_{\mathrm{F}}a_0=(3\pi^{2}\rho_0)^{1/3}a_0=\left(\frac{9\pi}{4}\right)^{1/3}\left(\frac{4\pi\rho_0a_0^{3}}{3}\right)^{1/3}=\left(\frac{9\pi}{4}\right)^{1/3}\frac{1}{r_{\mathrm{s}}}=\frac{1.9192}{r_{\mathrm{s}}} \tag{3.94} $$
根据表 3.1，有
$$ \frac{\hbar^{2}}{m_{\mathrm{e}}a_0^{2}}={E_{\mathrm h}}=1\,\mathrm{Hartree} \tag{3.95} $$
也就是说，在原子单位制下，取 $\hbar=m_{\mathrm{e}}=4\pi\varepsilon=1$，坐标单位取 $a_0$，能量单位为 Hartree，$1\,\mathrm{Hartree}=27.2\ \mathrm{eV}$。
在下面的讨论中，我们将详细、定量地推导该模型下体系能量的各项贡献。
库仑能
体系总的库仑能分别由电子与电子间库仑相互作用、正电荷背景间库仑相互作用以及电子-正电荷背景间库仑相互作用贡献。因为电子与正电荷背景均在空间中均匀分布，即
$$ \rho^{-}(\boldsymbol{r})=\rho^{+}(\boldsymbol{r})\equiv\rho_0=\frac{N}{V} \tag{3.96} $$
所以体系的库仑能为
$$ U_{\mathrm{Col}}=U_{\mathrm{ee}}+U_{\mathrm{II}}+U_{\mathrm{eI}}=e^{2}\left(\frac{N}{V}\right)^{2}\iint\left(\frac{1}{2}\frac{1}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}+\frac{1}{2}\frac{1}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}-\frac{1}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\right)\mathrm{d}\boldsymbol{r}\,\mathrm{d}\boldsymbol{r}'=0 \tag{3.97} $$
因此均匀分布的经典直接库仑项相互抵消；电子交换和关联贡献仍需另外计算。
动能与交换能
因为均匀电荷的经典直接库仑项与背景项相互抵消，因此，在均匀电子气模型的 Hartree–Fock 近似下，单电子的本征态可以用平面波 $|\,\boldsymbol{k}_i\rangle=\varOmega^{-1/2}\mathrm{e}^{\mathrm{i}\boldsymbol{k}_i\cdot\boldsymbol{r}}$ 表示。同时，体系的基态波函数可表示为 Slater 行列式，其中 $\boldsymbol{k}$ 的取值充满半径为 $k_{\mathrm{F}}$ 的费米球。不考虑自旋极化，也即每个 $|\,\boldsymbol{k}\rangle$ 态上占据两个电子，可具体写出该多体基态波函数 $\phi^{0}$：
$$ \phi^{0}=(N!)^{-1/2}\begin{vmatrix}\langle\boldsymbol{r}_1\,|\,\boldsymbol{k}_1\rangle\uparrow&\langle\boldsymbol{r}_2\,|\,\boldsymbol{k}_1\rangle\uparrow&\cdots&\langle\boldsymbol{r}_N\,|\,\boldsymbol{k}_1\rangle\uparrow\\\langle\boldsymbol{r}_1\,|\,\boldsymbol{k}_1\rangle\downarrow&\langle\boldsymbol{r}_2\,|\,\boldsymbol{k}_1\rangle\downarrow&\cdots&\langle\boldsymbol{r}_N\,|\,\boldsymbol{k}_1\rangle\downarrow\\\langle\boldsymbol{r}_1\,|\,\boldsymbol{k}_2\rangle\uparrow&\langle\boldsymbol{r}_2\,|\,\boldsymbol{k}_2\rangle\uparrow&\cdots&\langle\boldsymbol{r}_N\,|\,\boldsymbol{k}_2\rangle\uparrow\\\vdots&\vdots&&\vdots\\\langle\boldsymbol{r}_1\,|\,\boldsymbol{k}_{N/2}\rangle\downarrow&\langle\boldsymbol{r}_2\,|\,\boldsymbol{k}_{N/2}\rangle\downarrow&\cdots&\langle\boldsymbol{r}_N\,|\,\boldsymbol{k}_{N/2}\rangle\downarrow\end{vmatrix} \tag{3.98} $$
因为库仑能为零，所以正则 Hartree-Fock 方程（见式（3.34））中的 Fock 算符仅有动能算符及交换算符。动能算符的形式比较简单，而交换算符的普遍形式已经由方程（3.30）给出。在平面波基下，方程（3.34）可改写为（原子单位制下）
$$ \begin{aligned} {\hat F}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{r}}&=\left(-\frac{\boldsymbol{\nabla}^{2}}{2}-\sum_{j}{\hat K}_j\right)\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{r}}=\frac{k^{2}}{2}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{r}}-\frac{1}{\varOmega}\sum_{\boldsymbol{k}'}^{(\mathrm{occ})}\mathrm{e}^{\mathrm{i}\boldsymbol{k}'\cdot\boldsymbol{r}}\int\mathrm{e}^{-\mathrm{i}\boldsymbol{k}'\cdot\boldsymbol{r}'}\frac{1}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{r}'}\,\mathrm{d}\boldsymbol{r}'\\ &=\frac{k^{2}}{2}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{r}}-\frac{1}{\varOmega}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{r}}\sum_{\boldsymbol{k}'}^{(\mathrm{occ})}\int\frac{\mathrm{e}^{-\mathrm{i}(\boldsymbol{k}'-\boldsymbol{k})\cdot(\boldsymbol{r}-\boldsymbol{r}')}}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}'=\left(\frac{k^{2}}{2}-\frac{1}{\varOmega}\sum_{k'\lt k_{\mathrm{F}}}\frac{4\pi}{|\,\boldsymbol{k}-\boldsymbol{k}'\,|^{2}}\right)\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{r}}\\ &=\varepsilon_k\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{r}} \end{aligned} \tag{3.99} $$
上述计算过程中最后一步采用了 $1/|\,\boldsymbol{r}-\boldsymbol{r}'\,|$ 的傅里叶变换。其中本征值的交换能部分可以转化为费米球内的积分：
$$ \begin{aligned} \frac{1}{\varOmega}\sum_{k'\lt k_{\mathrm{F}}}\frac{4\pi}{|\,\boldsymbol{k}-\boldsymbol{k}'\,|^{2}}&=\frac{4\pi}{(2\pi)^{3}}\iiint\frac{k'^{2}\sin\theta\,\mathrm{d}k'\,\mathrm{d}\theta\,\mathrm{d}\varphi}{k^{2}-2kk'\cos\theta+k'^{2}}=\frac{1}{\pi k}\int_{0}^{k_{\mathrm{F}}}k'\ln\left|\frac{k+k'}{k-k'}\right|\mathrm{d}k'\\ &=\frac{k_{\mathrm{F}}}{\pi}F\left(\frac{k}{k_{\mathrm{F}}}\right) \end{aligned} \tag{3.100} $$
式中
$$ F(x)=1+\frac{1-x^{2}}{2x}\ln\left|\frac{1+x}{1-x}\right| \tag{3.101} $$
方程（3.99）至方程（3.101）表明，$|\,\boldsymbol{k}\rangle$ 确实是均匀电子气系统的 Fock 算符的本征函数，相应的本征值为
$$ \varepsilon(k)=\frac{k^{2}}{2}-\frac{k_{\mathrm{F}}}{\pi}F\left(\frac{k}{k_{\mathrm{F}}}\right) \tag{3.102} $$
图 3.2 给出了均匀电子气考虑和未考虑交换作用的约化色散曲线。可以看到，计入交换能会高估导带的带宽。而且在费米动量为 $k_{\mathrm{F}}$ 处，$\mathrm{d}\varepsilon_{\mathrm{HF}}/\mathrm{d}k$ 发散，从而导致此处的态密度为零，这显然与实际情况不符。这些都是利用 Hartree-Fock 方法处理均匀电子气的局限。
图 3.2 约化非相互作用均匀电子气及相互作用均匀电子气在不同电子密度下的 Hartree-Fock 能 $\varepsilon_{\mathrm{HF}}/\varepsilon_0$（$\varepsilon_0=\hbar^{2}k_{\mathrm{F}}^{2}/(2m)$）
由方程（3.102）出发，为了得到基态下每个电子的平均 Hartree-Fock 能量，需要对费米球内所有的态求和，并乘以 2（因为自旋简并度），而交换能部分还应再乘以 1/2（因为求和导致每个电子对交换能的贡献计入了两次）。因此
$$ E_0^{\mathrm{HF}}=2\sum_{k\lt k_{\mathrm{F}}}\frac{k^{2}}{2}-\sum_{k\lt k_{\mathrm{F}}}\frac{k_{\mathrm{F}}}{\pi}F\left(\frac{k}{k_{\mathrm{F}}}\right)=\frac{4\pi\varOmega}{8\pi^{3}}\int_{0}^{k_{\mathrm{F}}}k^{2}\left[k^{2}-\frac{k_{\mathrm{F}}}{\pi}F\left(\frac{k}{k_{\mathrm{F}}}\right)\right]\mathrm{d}k \tag{3.103} $$
容易求出，积分的第一项动能为 $\dfrac{8\pi\varOmega k_{\mathrm{F}}^{3}}{3(2\pi)^{3}}\times\dfrac{3}{5}\dfrac{k_{\mathrm{F}}^{2}}{2}$。
利用不定积分公式\cite{grosso2000solid}，有
$$ \int x(1-x^{2})\ln\frac{1+x}{1-x}\,\mathrm{d}x=\frac{x}{2}-\frac{x^{3}}{6}-\frac{1}{4}(1-x^{2})^{2}\ln\frac{1+x}{1-x} \tag{3.104} $$
可以得到积分的第二项交换能为 $\dfrac{8\pi\varOmega k_{\mathrm{F}}^{3}}{3(2\pi)^{3}}\times\dfrac{3}{4}\dfrac{k_{\mathrm{F}}}{\pi}$。
又因为电子总数 $N$ 为
$$ N=\frac{8\pi\varOmega k_{\mathrm{F}}^{3}}{3(2\pi)^{3}} \tag{3.105} $$
则可得
$$ \bar{E}_0=\frac{E_0^{\mathrm{HF}}}{N}=\frac{3}{5}\frac{k_{\mathrm{F}}^{2}}{2}-\frac{3}{4}\frac{k_{\mathrm{F}}}{\pi} \tag{3.106} $$
比较方程（3.105）和方程（3.106），可以看到，在均匀电子气系统中，可以认为每个电子的交换能正比于系统电子密度的 $\rho^{1/3}$；交换能的体密度正比于 $\rho^{4/3}$。如果认为交换能密度只与局域的电子密度有关，则可以近似认为
$$ V_{\mathrm{x}}(\boldsymbol{r})\propto\rho(\boldsymbol{r})^{1/3} \tag{3.107} $$
特别需要说明的是，交换作用只涉及自旋相同的电子态，因此方程（3.106）的第二项——交换能密度实际上应显含自旋指标\cite{martin2004electronic}：
$$ \varepsilon_{\mathrm{x}}^{\sigma}=-\frac{3}{4}\frac{k_{\mathrm{F}}^{\sigma}}{\pi}=-\frac{3}{4}\left(\frac{6\rho^{\sigma}}{\pi}\right)^{1/3} \tag{3.108} $$
对于自旋非极化情况，显然有
$$ \varepsilon_{\mathrm{x}}^{\uparrow}=\varepsilon_{\mathrm{x}}^{\downarrow}=-\frac{3}{4}\left(\frac{3\rho^{\mathrm{tot}}}{\pi}\right)^{1/3} \tag{3.109} $$
从本节开始的讨论可知，费米波矢 $\boldsymbol{k}_{\mathrm{F}}$ 可以表示为电子密度参数 $r_{\mathrm{s}}$ 的函数，因此平均动能与交换能也可以表示为 $r_{\mathrm{s}}$ 的函数。利用式（3.94）及式（3.95），可以将式（3.106）表示为
$$ \bar E_0=\left[\frac{3}{10}(k_{\mathrm F}a_0)^2-\frac{3}{4\pi}(k_{\mathrm F}a_0)\right]{E_{\mathrm h}} =\left(\frac{1.1050}{r_{\mathrm s}^{2}}-\frac{0.4581}{r_{\mathrm s}}\right){E_{\mathrm h}} \tag{3.110} $$
这里 ${E_{\mathrm h}}=1\,\mathrm{Hartree}$，两侧均为能量。
关联能
利用 Hartree-Fock 理论无法计入所有的多体效应，习惯上将除去交换作用以外所有其他的多体效应称为关联作用（correlation）。即使对于均匀电子气模型，平均每个电子的关联能 $E_{\mathrm{c}}$ 也很难精确求得。这个问题的定量解决是由 Gellmann 与 Brueckner 在 1957 年完成的\cite{gellmann1957correlation}。Gellmann 和 Brueckner 在高电子密度极限（$r_{\mathrm{s}}\to0$）下利用微扰将 $E_{\mathrm{c}}$ 展开，找出各阶的发散项，将这些项转为子级数各项积分的求和，从而求得 $E_{\mathrm{c}}$。定量的计算需要用到多体理论，这大大超出了本书的讨论范围，因此这里只给出非自旋极化三维电子气的高密度渐近结构（每电子能量单位为 Hartree）\cite{gellmann1957correlation,mahan1981many,macke1950wechselwirkungen,pines1953collective,onsager1966integrals}：
$$ E_{\mathrm c}(r_{\mathrm s})=\frac{1-\ln2}{\pi^2}\ln r_{\mathrm s} +C_{\mathrm{GB}}+O(r_{\mathrm s}\ln r_{\mathrm s}),\qquad r_{\mathrm s}\to0 \tag{3.111} $$
其中 $C_{\mathrm{GB}}$ 为高密度展开的常数项，不能由未经核对的分项数值相加得到。实际 LDA 计算常用 Perdew–Zunger 对 Ceperley–Alder 数据的参数化；其高密度分支为
$$ \varepsilon_{\mathrm c}^{\mathrm{PZ,high}}(r_{\mathrm s}) =0.0311\ln r_{\mathrm s}-0.048 +0.0020r_{\mathrm s}\ln r_{\mathrm s}-0.0116r_{\mathrm s} \quad (r_{\mathrm s}\le1) \tag{3.112} $$
式（3.112）是 Perdew–Zunger 参数化的高密度分支；它保持了 Gell-Mann–Brueckner 对数项的系数，常数 $-0.048$ 则是拟合参数，不应与原始渐近常数混同。若电子密度不符合上述极限，一般采用 Wigner 公式\cite{wigner1934interaction}：
$$ E_{\mathrm{c}}=-\frac{0.44}{r_{\mathrm{s}}+7.8} \tag{3.113} $$
图 3.3 相互作用均匀电子气的关联能 $E_{\mathrm{c}}(r_{\mathrm{s}})$
图 3.3 给出了由 Gellmann-Brueckner 公式（3.112）以及 Wigner 公式（3.113）得出的均匀电子气关联能 $E_{\mathrm{c}}(r_{\mathrm{s}})$，图中实线是作为基准的精确结果（见 3.2.5 节），标识“Ceperley-Alder”的结果取自 Perdew-Wang 对原始计算的拟合结果（见 3.2.5 节）。在高密度情况下（$r_{\mathrm{s}}\leqslant1$），Gellmann-Brueckner 公式非常精确。但是当 $r_{\mathrm{s}}\geqslant2$ 时，Wigner 公式更加准确。而当 $r_{\mathrm{s}}\gt 5$ 时，根据 Gellmann-Brueckner 公式计算出的关联能明显错误，这表明需要加入更高阶的项（如 $r_{\mathrm{s}}$ 以及 $r_{\mathrm{s}}\ln r_{\mathrm{s}}$ 等）加以修正。
电子对关联函数
Hartree-Fock 近似下，电子对关联函数 $g(\boldsymbol{r},\sigma;\boldsymbol{r}',\sigma')$ 由方程（3.61）或者方程（3.62）给出。$\tilde{\rho}(r)$ 可计算如下：
$$ \begin{aligned} \tilde{\rho}(r)&=\frac{2}{V}\sum_{\boldsymbol{k}}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{r}}=\frac{2}{8\pi^{3}}\int\mathrm{d}\boldsymbol{k}\,\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{r}}f[E(\boldsymbol{k})]\\ &=\frac{1}{4\pi^{3}}\sum_{l}\int_{0}^{k_{\mathrm{F}}}k^{2}\,\mathrm{d}k(2l+1)\mathrm{i}^{l}\mathrm{j}_l(kr)\cdot\int_{0}^{\pi}\mathrm{P}_0(\cos\theta)\mathrm{P}_l(\cos\theta)\sin\theta\,\mathrm{d}\theta\int_{0}^{2\pi}\mathrm{d}\varphi\\ &=\frac{1}{\pi^{2}}\int_{0}^{k_{\mathrm{F}}}k^{2}\mathrm{j}_0(kr)\,\mathrm{d}k\\ &=\frac{1}{\pi^{2}r^{3}}\big[\sin(rk_{\mathrm{F}})-rk_{\mathrm{F}}\cos(rk_{\mathrm{F}})\big] \end{aligned} \tag{3.114} $$
设体系自旋非极化，$\tilde{\rho}^{\uparrow}(r)=\tilde{\rho}^{\downarrow}(r)=\tilde{\rho}(r)/2$，而 $\rho_0$ 由式（3.93）给出，将其代入式（3.62），即得
$$ g_{\mathrm{x}}(r)=1-\frac{9}{2}\left[\frac{\sin rk_{\mathrm{F}}-rk_{\mathrm{F}}\cos(rk_{\mathrm{F}})}{(rk_{\mathrm{F}})^{3}}\right]^{2} \tag{3.115} $$
图 3.4 给出了 $g_{\mathrm{x}}(r)$ 的具体形状。
图 3.4 相互作用的均匀电子气在不同电子密度下的对关联函数 $g_{\mathrm{x}}(r)$
Hartree-Fock 方程的数值求解和基组选取
对于一个实际分子体系的轨道，我们通常用类似于原子轨道的基组对波函数进行展开。Slater 形式的原子轨道是最自然的选择。普遍地，中心在 $\boldsymbol{r}_A$ 终点处、衰减指数为 $\kappa$ 的 Slater 函数可以表示为
$$ S_n(\kappa,A)=(2\kappa)^{n+\frac{1}{2}}\big[(2n)!\big]^{-\frac{1}{2}}r^{n-1}\mathrm{e}^{-\kappa r} \tag{3.116} $$
式中：$A$ 代表空间中的原子核的坐标 $(A_x,A_y,A_z)$。
由于实际计算过程涉及多中心积分，且涉及 Slater 原子轨道的多中心积分，计算量庞大，人们发展了高斯函数作为波函数的展开基组以简化计算。目前许多常用的计算软件如 Gaussian16 等都支持高斯基组。高斯基组受欢迎的根本原因在于两中心的高斯函数的积分可以简化为单中心的高斯积分，此性质递推使用，可以极大简化多电子体系的 Hartree-Fock 方程中的三中心、四中心积分的计算。下面对高斯基组做一个简单介绍。我们用 $G(\alpha,A)$ 表示中心在 $\boldsymbol{r}_A$ 终点 $A$ 处、衰减指数为 $\alpha$ 的未归一化高斯函数，用 $\tilde{G}(\alpha,A)$ 表示归一化的高斯函数，则它们的定义式分别为
$$ G(\alpha,\boldsymbol{r}-\boldsymbol{r}_A)=\mathrm{e}^{-\alpha|\boldsymbol{r}-\boldsymbol{r}_A|^{2}} \tag{3.117} $$
$$ \tilde{G}(\alpha,\boldsymbol{r}-\boldsymbol{r}_A)=\left(\frac{2\alpha}{\pi}\right)^{3/4}\mathrm{e}^{-\alpha|\boldsymbol{r}-\boldsymbol{r}_A|^{2}} \tag{3.118} $$
而广义高斯函数可以写成如下形式：
$$ G(\alpha,\boldsymbol{r}-\boldsymbol{r}_A,l,m,n)=(x-x_A)^{l}(y-y_A)^{m}(z-z_A)^{n}\mathrm{e}^{-\alpha|\boldsymbol{r}-\boldsymbol{r}_A|^{2}} \tag{3.119} $$
$$ \tilde{G}(\alpha,\boldsymbol{r}-\boldsymbol{r}_A,l,m,n)=N(x-x_A)^{l}(y-y_A)^{m}(z-z_A)^{n}\mathrm{e}^{-\alpha|\boldsymbol{r}-\boldsymbol{r}_A|^{2}} \tag{3.120} $$
式中：$x_A$、$y_A$、$z_A$ 分别为中心 $A$ 的笛卡尔坐标。归一化常数 $N$ 由下式得到：
$$ N=\left(\frac{2\alpha}{\pi}\right)^{3/4}\left[\frac{(4\alpha)^{l+m+n}}{(2l-1)!!(2m-1)!!(2n-1)!!}\right]^{1/2} \tag{3.121} $$
式中：!! 代表双阶乘，$(2l-1)!!=(2l-1)(2l-3)\cdots(3)(1)$。
在实际分子轨道计算中，由于高斯函数与原子轨道（更接近 Slater 函数）相差较远，因此通常用一组高斯函数来线性拟合 Slater 基组，这样的高斯函数的集合称为编缩高斯基组。通常用 $K$ 个高斯函数就记为 STO-$K$G，比如量化计算中最常用的最小基组 STO-3G 代表用三个高斯函数来线性展开一个近似为 Slater 形式的函数。人们通常用 $G(\alpha,\boldsymbol{r}-\boldsymbol{r}_A,l=0,m=0,n=0)$ 来拟合 s 电子的 Slater 轨道，用 $G(\alpha,\boldsymbol{r}-\boldsymbol{r}_A,l=1,m=0,n=0)$ 来拟合 $p_x$ 电子的 Slater 轨道，用 $G(\alpha,\boldsymbol{r}-\boldsymbol{r}_A,l=1,m=1,n=0)$ 等展开 $d$ 电子的 Slater 轨道。接下来我们说明如何用最小二乘法确定 Slater 轨道的高斯展开系数。以 1s 的 Slater 函数为例：
$$ S_{1s}(\kappa=1.0,\boldsymbol{r}-\boldsymbol{r}_A)=\left(\frac{1}{\pi}\right)^{1/2}\mathrm{e}^{-|\boldsymbol{r}-\boldsymbol{r}_A|} \tag{3.122} $$
其中，取衰减系数 $\kappa=1.0$。
我们需要将式（3.122）展开为高斯函数的线性叠加，并且利用非线性的最小二乘法确定各个高斯函数的衰减系数 $\alpha_i$ 及高斯函数前的展开系数 $c_i$。
$$ S_{1s}(\kappa,\boldsymbol{r}-\boldsymbol{r}_A)\approx\sum_{i=1}^{K}c_iG(\alpha_i,\boldsymbol{r}-\boldsymbol{r}_A) \tag{3.123} $$
优化系数后得到 STO-1G、STO-2G、STO-3G 的结果分别如下：
$$ \begin{aligned} \text{STO-1G：}S_{1s}(\kappa=1.0,\boldsymbol{r}-\boldsymbol{r}_A)&=\tilde{G}(0.270950,\boldsymbol{r}-\boldsymbol{r}_A)\\ \text{STO-2G：}S_{1s}(\kappa=1.0,\boldsymbol{r}-\boldsymbol{r}_A)&=0.678914\times\tilde{G}(0.151623,\boldsymbol{r}-\boldsymbol{r}_A)\\ &\quad+0.430129\times\tilde{G}(0.851819,\boldsymbol{r}-\boldsymbol{r}_A)\\ \text{STO-3G：}S_{1s}(\kappa=1.0,\boldsymbol{r}-\boldsymbol{r}_A)&=0.444635\times\tilde{G}(0.109818,\boldsymbol{r}-\boldsymbol{r}_A)\\ &\quad+0.535328\times\tilde{G}(0.405771,\boldsymbol{r}-\boldsymbol{r}_A)\\ &\quad+0.154329\times\tilde{G}(2.22766,\boldsymbol{r}-\boldsymbol{r}_A) \end{aligned} $$
图 3.5 给出了 Slater 函数分别用 1、2、3 个高斯函数拟合结果的对比。可以看到，随着用于拟合的高斯基组的不断增大，编缩的高斯基组也越来越接近 Slater 轨道。但是，同时也可以看到，Slater 函数在原子核所在的空间坐标处导数不连续（这是由库仑势在距离等于零时的发散引起的），而高斯函数的线性组合则在原点处导数平滑连续。
将基组展开成高斯函数后，可以利用高斯函数的约化法则简化双中心的积分计算。只需将 $|\,\boldsymbol{r}-\boldsymbol{r}'\,|^{2}$ 用 $r^{2}+r'^{2}-2\boldsymbol{r}\cdot\boldsymbol{r}'$ 展开。易证明以下等式成立：
图 3.5 用 STO-1G、STO-2G、STO-3G 的高斯基组分别拟合 Slater 函数
$$ \begin{aligned} &e^{-\alpha|\boldsymbol r-\boldsymbol r_A|^2}e^{-\beta|\boldsymbol r-\boldsymbol r_B|^2}\\ &=e^{-\frac{\alpha\beta}{\alpha+\beta}|\boldsymbol r_A-\boldsymbol r_B|^2} e^{-(\alpha+\beta)|\boldsymbol r-\frac{\alpha\boldsymbol r_A+\beta\boldsymbol r_B}{\alpha+\beta}|^2} \end{aligned} \tag{3.124} $$
如图 3.6 所示，以 $A$ 点和 $B$ 点为中心的两个高斯函数的乘积可以约化成一个与 $r_{AB}$ 相关的常数与一个以 $AB$ 连线上的重心 $C$ 为中心的高斯函数的乘积。
图 3.6 高斯函数双中心积分约化示意
如果不考虑高斯函数的归一化，上述规律可以表示为
$$ G(\alpha,\boldsymbol{r}-\boldsymbol{r}_A)G(\beta,\boldsymbol{r}-\boldsymbol{r}_B)=K\cdot G(\alpha+\beta,\boldsymbol{r}-\boldsymbol{r}_C) \tag{3.125} $$
如果进一步考虑归一化系数，则可以表示为
$$ \tilde{G}(\alpha,\boldsymbol{r}-\boldsymbol{r}_A)\tilde{G}(\beta,\boldsymbol{r}-\boldsymbol{r}_B)=\tilde{K}\cdot\tilde{G}(\alpha+\beta,\boldsymbol{r}-\boldsymbol{r}_C) \tag{3.126} $$
其中
$$ K=\exp\left(-\frac{\alpha\beta}{\alpha+\beta}r_{AB}^{2}\right) \tag{3.127} $$
$$ \tilde{K}=\left[\frac{2\alpha\beta}{\pi(\alpha+\beta)}\right]^{3/4}\exp\left(-\frac{\alpha\beta}{\alpha+\beta}r_{AB}^{2}\right) \tag{3.128} $$
$$ \boldsymbol{r}_C=\frac{\alpha}{\alpha+\beta}\boldsymbol{r}_A+\frac{\beta}{\alpha+\beta}\boldsymbol{r}_B \tag{3.129} $$
上述结果容易推广到广义高斯函数的双中心积分。如果忽略高斯函数前面的归一化系数，其约化规律如下
$$ \begin{aligned} &G(\alpha,\boldsymbol{r}-\boldsymbol{r}_A,l,m,n)G(\beta,\boldsymbol{r}-\boldsymbol{r}_B,l',m',n')\\ &=[(x-x_C)+R_x]^{l}[(x-x_C)+R'_x]^{l'}[(y-y_C)+R_y]^{m}[(y-y_C)+R'_y]^{m'}\\ &\quad\cdot[(z-z_C)+R_z]^{n}[(z-z_C)+R'_z]^{n'}\cdot K\cdot G(\alpha+\beta,\boldsymbol{r}-\boldsymbol{r}_C) \end{aligned} \tag{3.130} $$
其中
$$ R=\frac{\beta}{\alpha+\beta}(\boldsymbol{r}_B-\boldsymbol{r}_A) \tag{3.131} $$
$$ R'=\frac{\alpha}{\alpha+\beta}(\boldsymbol{r}_A-\boldsymbol{r}_B) \tag{3.132} $$
通过递推运用高斯函数的约化规则，可以将多中心的积分转变成单中心的积分。Hartree-Fock 方程矩阵表达中的积分项可以有解析表达式。现以 1s 型的高斯函数为例说明计算简化步骤。
假设我们感兴趣的是分别位于 $\boldsymbol{r}_A$ 和 $\boldsymbol{r}_B$ 终点处的两个高斯基组 $\zeta_\mu$ 和 $\zeta_\nu$ 之间的矩阵元
$$ \zeta_\mu=\tilde{G}^{*}(\alpha,\boldsymbol{r}-\boldsymbol{r}_A) \tag{3.133} $$
$$ \zeta_\nu=\tilde{G}^{*}(\beta,\boldsymbol{r}-\boldsymbol{r}_B) \tag{3.134} $$
重叠积分 $S_{\mu\nu}$ 的计算如下：
$$ \begin{aligned} S_{\mu\nu}&=\int\tilde{G}^{*}(\alpha,\boldsymbol{r}-\boldsymbol{r}_A)\tilde{G}(\beta,\boldsymbol{r}-\boldsymbol{r}_B)\,\mathrm{d}\boldsymbol{r}\\ &=\left(\frac{2\sqrt{\alpha\beta}}{\pi}\right)^{3/2}K\int_{-\infty}^{\infty}G(\alpha+\beta,\boldsymbol{r}-\boldsymbol{r}_C)\,\mathrm{d}\boldsymbol{r}\\ &=\left(\frac{2\sqrt{\alpha\beta}}{\alpha+\beta}\right)^{3/2}\exp\left(-\frac{\alpha\beta}{\alpha+\beta}r_{AB}^{2}\right) \end{aligned} \tag{3.135} $$
式（3.135）的推导中用到了定积分
$$ \int_{\mathbb R^3}\exp\!\left[-(\alpha+\beta)|\boldsymbol r|^2\right]\,\mathrm d^3\boldsymbol r =\left(\frac{\pi}{\alpha+\beta}\right)^{3/2} \tag{3.136} $$
动能项矩阵元 $T_{\mu\nu}$ 可以写为
$$ \begin{aligned} T_{\mu\nu} &=\int\!\mathrm d^3\boldsymbol r\;\tilde G(\alpha,\boldsymbol r-\boldsymbol r_A) \left(-\tfrac12\nabla^2\right)\tilde G(\beta,\boldsymbol r-\boldsymbol r_B)\\ &=\int\!\mathrm d^3\boldsymbol r\;\tilde G(\alpha,\boldsymbol r-\boldsymbol r_A) \left[3\beta-2\beta^2|\boldsymbol r-\boldsymbol r_B|^2\right] \tilde G(\beta,\boldsymbol r-\boldsymbol r_B)\\ &=S_{\mu\nu}\left[\frac{3\alpha\beta}{\alpha+\beta} -\frac{2\alpha^2\beta^2}{(\alpha+\beta)^2}r_{AB}^2\right]. \end{aligned} $$
$$ =\left(\frac{2\sqrt{\alpha\beta}}{\alpha+\beta}\right)^{3/2}\cdot\left[\frac{3\alpha\beta}{\alpha+\beta}-\frac{2\alpha^{2}\beta^{2}}{(\alpha+\beta)^{2}}r_{AB}^{2}\right]\exp\left(-\frac{\alpha\beta}{\alpha+\beta}r_{AB}^{2}\right) \tag{3.137} $$
式（3.137）的推导中用到了定积分
$$ \int_{-\infty}^{\infty}x^{2}\exp[-(\alpha+\beta)x^{2}]\,\mathrm{d}x=\frac{1}{2(\alpha+\beta)}\sqrt{\frac{\pi}{\alpha+\beta}}=\frac{1}{2(\alpha+\beta)}\int_{-\infty}^{\infty}\exp[-(\alpha+\beta)x^{2}]\,\mathrm{d}x $$
电子-核库仑吸引矩阵元可由高斯乘积定理以及径向积分直接求得。记 $p=\alpha+\beta$、$\boldsymbol r_C=(\alpha\boldsymbol r_A+\beta\boldsymbol r_B)/p$，并用 $N_\alpha=(2\alpha/\pi)^{3/4}$ 表示归一化常数。积分恒等式为
$$ \int_{\mathbb R^3}\frac{e^{-p|\boldsymbol r-\boldsymbol r_C|^2}}{|\boldsymbol r-\boldsymbol r_M|}\,\mathrm d^3\boldsymbol r =\frac{2\pi}{p}F_0\!\left(p|\boldsymbol r_C-\boldsymbol r_M|^2\right),\qquad p\gt 0 \tag{3.138} $$
其中零阶 Boys 函数定义为
$$ F_0(x)=\int_0^1 e^{-xt^2}\,\mathrm dt,\qquad x\geqslant0 \tag{3.139} $$
因此，代入式（3.124）即可得到
$$ \begin{aligned} V_{\mu\nu}^{M} &=-Z_M N_\alpha N_\beta \int\!\mathrm d^3\boldsymbol r\; \frac{e^{-\alpha|\boldsymbol r-\boldsymbol r_A|^2-\beta|\boldsymbol r-\boldsymbol r_B|^2}} {|\boldsymbol r-\boldsymbol r_M|}\\ &=-Z_M N_\alpha N_\beta\frac{2\pi}{\alpha+\beta} \exp\!\left[-\frac{\alpha\beta}{\alpha+\beta}r_{AB}^2\right] F_0\!\left((\alpha+\beta)|\boldsymbol r_C-\boldsymbol r_M|^2\right)\\ &=-\frac{2^{5/2}(\alpha\beta)^{3/4}}{\sqrt\pi(\alpha+\beta)}Z_M \exp\!\left[-\frac{\alpha\beta}{\alpha+\beta}r_{AB}^2\right] F_0\!\left((\alpha+\beta)|\boldsymbol r_C-\boldsymbol r_M|^2\right) \end{aligned} \tag{3.140} $$
在实际计算过程中，$F_0$ 可以很容易通过程序包自带的误差函数得到：
$$ F_0(x)=\begin{cases}\dfrac{1}{2}\sqrt{\dfrac{\pi}{x}}\,\mathrm{erf}\sqrt{x},&x\gt 0,\\1,&x=0.\end{cases} \tag{3.141} $$
方程（3.88）所呈现的双电子积分的约化思路与电子-核库仑吸引矩阵元的简化过程类似。如图 3.7 所示，我们首先可以将分别位于点 $A$、$C$ 处的高斯函数乘积约化成位于 $A$、$C$ 连线中心的点 $M$ 处的高斯函数，将分别位于点 $B$、$D$ 处的高斯函数乘积约化成位于 $B$、$D$ 连线中心的点 $N$ 处的高斯函数。分别位于点 $M$、$N$ 处的高斯函数则可以类似地利用计算电子-核库仑吸引矩阵元时的傅里叶变换技巧进行简化。由于和计算电子-核库仑吸引矩阵元的相似度较高，我们省略了具体的推导过程，直接给出该积分最终的表达式：
$$ \begin{aligned} \langle AB\,|\,CD\rangle&=\iint\mathrm{d}\boldsymbol{r}_1\,\mathrm{d}\boldsymbol{r}_2\,\tilde{G}^{*}(\alpha,\boldsymbol{r}_1-\boldsymbol{r}_A)\tilde{G}^{*}(\beta,\boldsymbol{r}_2-\boldsymbol{r}_B)\frac{1}{|\,\boldsymbol{r}_1-\boldsymbol{r}_2\,|}\tilde{G}(\gamma,\boldsymbol{r}_1-\boldsymbol{r}_C)\tilde{G}(\delta,\boldsymbol{r}_2-\boldsymbol{r}_D)\\ &=\frac{16(\alpha\beta\gamma\delta)^{3/4}}{(\alpha+\gamma)(\beta+\delta)\sqrt{\pi(\alpha+\beta+\gamma+\delta)}}\exp\left(-\frac{\alpha\gamma}{\alpha+\gamma}\,|\,\boldsymbol{r}_A-\boldsymbol{r}_C\,|^{2}\right.\\ &\quad\left.-\frac{\beta\delta}{\beta+\delta}\,|\,\boldsymbol{r}_B-\boldsymbol{r}_D\,|^{2}\right)\cdot F_0\left[\frac{(\alpha+\gamma)(\beta+\delta)}{\alpha+\beta+\gamma+\delta}\,|\,\boldsymbol{r}_M-\boldsymbol{r}_N\,|^{2}\right] \end{aligned} \tag{3.142} $$
图 3.7 高斯函数双电子积分约化示意
可见高斯函数的引入大大简化了 Hartree-Fock 方程求解过程中的其他交换能、库仑斥能、交换能等积分项，用高斯积分的解析表达式结合约化法则替代耗时的数值积分过程，可以加速求解的过程。
$X_\alpha$ 方法和超越 Hartree-Fock 近似
现代的分子轨道计算理论很大程度上是以 Hartree-Fock 方法为基础的。在这里我们给出与 Hartree-Fock 方法联系非常紧密的两种方法。第一种是 $X_\alpha$ 方法，从实用角度看可以认为它是对 Hartree-Fock 方法的一种简化，从另外一个角度看，也可以认为它是密度泛函理论的前身。在 3.1.12 节中我们可以看到，Hartree–Fock 非局域交换算符的计算涉及多中心的双电子积分，这是最为耗时的一部分，因此人们借鉴均匀电子气的结果（见式（3.108）），将交换作用近似表示为局域交换势（原子单位制）
$$ \mu_{\mathrm{x}}^{X_\alpha}(\boldsymbol r)=-3\alpha\left(\frac{3\rho(\boldsymbol r)}{8\pi}\right)^{1/3} \tag{3.143} $$
式中：$\alpha$ 是一个可调参数。实际计算中取 $\alpha=0.7$ 可以得到较好的结果。
另外，在现代量子化学计算中，许多高精度方法是在 Hartree-Fock 方法的基础上发展起来的，统称为超越 Hartree-Fock 近似的方法。其中最主要的一种为组态相互作用（configuration interaction，CI）方法。在这里，我们简要地介绍一下该方法。回顾 Hartree-Fock 方法的基本假设，可以看到其核心是用单行列式的波函数近似体系的真实基态，其本征解给出基态能量的上限估计。显然，单行列式的波函数不足以构成能展开一个多电子体系的完备基组，更精确的近似是行列式波函数的线性组合。假设有 $2K$ 个单电子轨道，从中挑出 $N$ 个轨道组成行列式的可能性为 $\mathrm{C}_N^{2K}$。可见，即使是一个较小的体系，将其波函数用所有的行列式波函数展开，计算量也是相当大的。因此，在实际计算中，人们通常取几阶较低的近似。如果用于展开的行列式中，$N$ 个轨道的组成和 Hartree-Fock 中的单行列式分别仅相差一个轨道，则称为单激发的组态相互作用；如果各个行列式与 Hartree-Fock 行列式相差两个轨道，则称该组态相互作用为双激发的组态相互作用。
我们从 Hartree-Fock 基态波函数出发，有
$$ |\,\phi\rangle=\phi(\boldsymbol{x}_1,\boldsymbol{x}_2,\cdots,\boldsymbol{x}_N)=|\,\xi_i(1)\xi_j(2)\cdots\xi_l(a)\xi_m(b)\cdots\xi_k(N)\rangle_{\mathrm{S}} \tag{3.144} $$
用如下的符号代表第 $a$ 个电子从 $\xi_l$ 轨道到 $\xi_p$ 轨道的激发：
$$ |\,\phi_l^{p}\rangle=|\,\xi_i(1)\xi_j(2)\cdots\xi_p(a)\xi_m(b)\cdots\xi_k(N)\rangle_{\mathrm{S}} \tag{3.145} $$
$|\,\phi_{lm}^{pq}\rangle$ 代表第 $a$ 个电子从 $\xi_l$ 轨道到 $\xi_p$ 轨道激发的同时，第 $b$ 个电子从 $\xi_m$ 轨道被激发到 $\xi_q$ 轨道，即
$$ |\,\phi_{lm}^{pq}\rangle=|\,\xi_i(1)\xi_j(2)\cdots\xi_p(a)\xi_q(b)\cdots\xi_k(N)\rangle_{\mathrm{S}} \tag{3.146} $$
如果用一组单电子波函数 $\{\xi_i\}$ 作为完备基组，则符合反对称性质的体系真实波函数可以如下行列式的形式完备展开（其中 $l\lt m$，$p\lt q$，该限制条件保证不重复对各组态计数）：
$$ \psi=c_0\,|\,\phi\rangle+\sum_{l,p}c_l^{p}\,|\,\phi_l^{p}\rangle+\sum_{l\lt m,\,p\lt q}c_{lm}^{pq}\,|\,\phi_{lm}^{pq}\rangle+\cdots \tag{3.147} $$
因此，从理论上来讲，只要有足够的计算能力，我们就可以穷举所有的组态，并计算哈密顿量在各个组态下的矩阵元，则对角化哈密顿量得到的最低能级就是体系的基态能量。在变分过程中，由于引入了更多的自由度，通过 CI 方法计算出来的基态能量要低于用 Hartree-Fock 方法得到的基态能量，两者之间的差值通常被定义为关联能。
现以氢分子 H2 的电子结构为例，阐述组态相互作用的基本计算过程及其与 Hartree-Fock 方法的区别。为简洁起见，仅考虑以氢原子的 1s 轨道作为基组，并将 H2 的分子波函数近似为原子轨道的线性组合。对角化单电子哈密顿量后，可以得到四个单电子自旋轨道：
$$ \begin{cases} \xi_1=\phi_{1s}^{+}\alpha(s)\\ \xi_2=\phi_{1s}^{+}\beta(s)\\ \xi_3=\phi_{1s}^{-}\alpha(s)\\ \xi_4=\phi_{1s}^{-}\beta(s) \end{cases} \tag{3.148} $$
式中：$\xi_1$ 和 $\xi_2$ 为简并的成键轨道；$\xi_3$ 和 $\xi_4$ 为简并的反键轨道。Hartree-Fock 近似下的基态轨道为 $|\,\phi\rangle=|\,\xi_1\xi_2\rangle_{\mathrm{S}}$。在组态相互作用计算中，我们需要从四个自旋轨道中挑出两个构成 Slater 行列式，作为电子的一个组态，则这样的取法共有 $\mathrm{C}_2^{4}=6$ 种，如图 3.8 所示，它们分别是：
（1）基态组态 $|\,\xi_1\xi_2\rangle_{\mathrm{S}}$；
（2）单激发组态 $|\,\xi_3\xi_2\rangle_{\mathrm{S}}$、$|\,\xi_4\xi_2\rangle_{\mathrm{S}}$、$|\,\xi_1\xi_3\rangle_{\mathrm{S}}$、$|\,\xi_1\xi_4\rangle_{\mathrm{S}}$；
（3）双激发组态 $|\,\xi_3\xi_4\rangle_{\mathrm{S}}$。
图 3.8 H2 分子中考虑 1s 轨道时所有可能的组态
和 Hartree-Fock 基态波函数类似，CI 基态波函数也具有偶宇称，即在空间反演操作下该基态波函数保持不变。因此，其基态波函数可以表示成基态组态和双激发组态的线性叠加：
$$ |\,\phi_{\mathrm{CI}}\rangle=c_0\,|\,\xi_1\xi_2\rangle_{\mathrm{S}}+c_{12}^{34}\,|\,\xi_3\xi_4\rangle_{\mathrm{S}} \tag{3.149} $$
不可约的二阶哈密顿矩阵可以写为
$$ \begin{aligned} \boldsymbol{H}&=\begin{bmatrix}\langle\phi|{\hat H}|\phi\rangle&\langle\phi|{\hat H}|\phi_{12}^{34}\rangle\\\langle\phi_{12}^{34}|{\hat H}|\phi\rangle&\langle\phi_{12}^{34}|{\hat H}|\phi_{12}^{34}\rangle\end{bmatrix}\\ &=\begin{bmatrix}\langle1|h|1\rangle+\langle2|h|2\rangle+\langle12|12\rangle-\langle12|21\rangle&\langle12|34\rangle-\langle12|43\rangle\\\langle34|12\rangle-\langle34|21\rangle&\langle3|h|3\rangle+\langle4|h|4\rangle+\langle34|34\rangle-\langle34|43\rangle\end{bmatrix} \end{aligned} $$
将自旋部分积分后，哈密顿矩阵可以进一步简化为
$$ \boldsymbol{H}=\begin{bmatrix}2\langle\phi_{1s}^{+}|h|\phi_{1s}^{+}\rangle+\langle\phi_{1s}^{+}\phi_{1s}^{+}|\phi_{1s}^{+}\phi_{1s}^{+}\rangle&\langle\phi_{1s}^{+}\phi_{1s}^{-}|\phi_{1s}^{+}\phi_{1s}^{-}\rangle\\\langle\phi_{1s}^{-}\phi_{1s}^{+}|\phi_{1s}^{-}\phi_{1s}^{+}\rangle&2\langle\phi_{1s}^{-}|h|\phi_{1s}^{-}\rangle+\langle\phi_{1s}^{-}\phi_{1s}^{-}|\phi_{1s}^{-}\phi_{1s}^{-}\rangle\end{bmatrix} $$
对角化此哈密顿矩阵，即可得到仅考虑 1s 轨道时的基态能量和基态波函数。值得指出的是，矩阵元 $H_{11}$ 即为 Hartree-Fock 方法的基态能量。因此，由 CI 方法得到的基态能量是考虑激发组态时对系统总能的一个修正。