本文是「第一性原理的微观计算模拟」系列的第 1 篇(共 11 篇),内容整理自同名书稿第 3 章,公式编号与原书一致。文中引用的文献依据原书参考文献清单整理并列于文末,按原书清单顺序从 1 开始编号。\nocite{*}
分子轨道理论是 20 世纪初由 F. Hund 和 R. S. Mulliken 发展起来的化学键理论,被广泛用于描述不同分子的结构和性质。早期的价键理论无法充分解释某些分子如何包含两个或多个等效键,以及其键序介于单键和双键之间的现象,而分子轨道理论比价键理论的解释能力要强大一些,它所描述的轨道可以准确反映所研究分子的几何形状。
玻恩-奥本海默近似
哈特里-福克(Hartree-Fock)方法,或者说大多数第一性原理计算方法的基础都是不含时薛定谔方程(见式(2.6))。从本质上看,这些计算方法可以视为对薛定谔方程所采取的不同的近似求解方法。设
$$ {\hat H}\Phi(\{\boldsymbol{r}_i\},\{\boldsymbol{R}_A\})={E}\Phi(\{\boldsymbol{r}_i\},\{\boldsymbol{R}_A\}) \tag{3.1} $$
式(3.1)中,$\hat H$ 是整个电子与原子核体系的哈密顿算符,作用于总波函数 $\Phi$;$E$ 是对应的总能量本征值。上帽表示这里的 $H$ 是作用于波函数的算符。
考虑到体系中的核运动的动能,电子的动能,核与核、核与电子、电子与电子之间的库仑相互作用,在国际单位制下哈密顿量可以表示为
$$ \begin{aligned} {\hat H}=-\sum_{A=1}^{M}\frac{\hbar^{2}}{2m_A}\boldsymbol{\nabla}_A^{2}-\sum_{i=1}^{N}\frac{\hbar^{2}}{2m_{\mathrm{e}}}\boldsymbol{\nabla}_i^{2}+\sum_{A=1}^{M}\sum_{B\gt A}^{M}\frac{Z_AZ_Be^{2}}{4\pi\varepsilon_0R_{AB}}\\ +\sum_{i=1}^{N}\sum_{j\gt i}^{N}\frac{e^{2}}{4\pi\varepsilon_0r_{ij}}-\sum_{i=1}^{N}\sum_{A=1}^{M}\frac{Z_Ae^{2}}{4\pi\varepsilon_0r_{iA}} \end{aligned} \tag{3.2} $$
图 3.1 粒子间相互作用力的示意图
式中:$A$、$B$ 分别为核的标号;$i$、$j$ 为电子的标号;$m_A$、$m_{\mathrm{e}}$ 分别为核子和电子的质量;$Z_A$、$Z_B$ 为各个核子所带的正电荷;$R_{AB}$、$r_{ij}$、$r_{iA}$ 分别为核与核、电子与电子、核与电子之间的距离;$\varepsilon_0$ 为真空介电常数,$\varepsilon_0=8.85419\times10^{-12}\ \mathrm{C^{2}\cdot J^{-1}\cdot m^{-1}}$;$e$ 为单位电荷,$e=1.6022\times10^{-19}\ \mathrm{C}$。该表达式的前两项分别是核子和电子的动能,后三项分别是核与核、电子与电子、核与电子间库仑相互作用,如图 3.1 所示。需要指出的是,求解多体薛定谔方程的最大困难在于,电子与电子间相互作用项的存在,使得薛定谔方程无法采用分离变量的方法求解。因此,如何引入适当的近似(平均场),将一个多体问题有效地转化为单体问题,是解决该问题的关键所在。Hartree-Fock 方法是实现以上目的的有效近似方法之一。后面提到的密度泛函,则是基于另外一种思路引入的。
哈密顿量中的系数,如 $\dfrac{\hbar^{2}}{2m_A}$、$\dfrac{Z_AZ_Be^{2}}{4\pi\varepsilon_0R_{AB}}$ 等,不仅使方程显得比较烦琐,而且在具体的数值计算中涉及很大的常数。由于计算机的数值模拟的精度有限,乘以或者除以很大的常数将大大降低计算的精度,因此,为了讨论和计算方便,人们通常采用原子单位制来重写方程,使之更易于处理。
表 3.1 给出了改写薛定谔方程时相关的几个量的国际单位制和原子单位制。改用原子单位制后,哈密顿量可简化为
$$ {\hat H}=-\sum_{A=1}^{M}\frac{1}{2M_A}\boldsymbol{\nabla}_A^{2}-\sum_{i=1}^{N}\frac{1}{2}\boldsymbol{\nabla}_i^{2}+\sum_{A=1}^{M}\sum_{B\gt A}^{M}\frac{Z_AZ_B}{R_{AB}}+\sum_{i=1}^{N}\sum_{j\gt i}^{N}\frac{1}{r_{ij}}-\sum_{i=1}^{N}\sum_{A=1}^{M}\frac{Z_A}{r_{iA}} \tag{3.3} $$
式中:$M_A=m_A/m_{\mathrm{e}}$。
表 3.1 国际单位制和原子单位制之间的对应关系
| 国际单位制 | 原子单位制 | |
|---|---|---|
| 质量 | 千克(kg) | 电子质量 $m_{\mathrm{e}}=9.1094\times10^{-31}\ \mathrm{kg}$ |
| 电荷 | 库仑(C) | 单位电荷 $e=1.6022\times10^{-19}\ \mathrm{C}$ |
| 角动量 | 千克二次方米每秒($\mathrm{kg\cdot m^{2}\cdot s^{-1}}$) | ${\hbar=1.0546\times10^{-34}}\ \mathrm{J\cdot s}$ |
| 介电常数 | 法每米($\mathrm{F\cdot m^{-1}}$) | $4\pi\varepsilon_0=1.1127\times10^{-10}\ \mathrm{F\cdot m^{-1}}$ |
| 长度 | 米(m) | 玻尔半径 $a_0=\dfrac{4\pi\varepsilon_0\hbar^{2}}{m_{\mathrm{e}}e^{2}}=5.2918\times10^{-11}\ \mathrm{m}$ |
| 能量 | 焦耳(J) | ${E_{\mathrm h}}=1\ \mathrm{Hartree}={\dfrac{m_{\mathrm{e}}e^{4}}{(4\pi\varepsilon_0)^{2}\hbar^{2}}}=4.3597\times10^{-18}\ \mathrm{J}$ |
由于 $\dfrac{Z_A}{r_{iA}}$ 项的存在,无法简单地对电子和核运动方程进行分离变量处理。将此薛定谔方程进一步简化的一个关键的近似方法是玻恩-奥本海默近似(Born-Oppenheimer approximation)。 由于核的质量远大于电子质量,因此缓慢的核运动方程和电子运动方程可以被有效地分开求解,而不会引入大的误差。从数学角度,我们可以显式地将总波函数 $\Phi(\{\boldsymbol{r}_i\},\{\boldsymbol{R}_A\})$ 中的核运动部分分离出来,有
$$ \Phi(\{\boldsymbol{r}_i\},\{\boldsymbol{R}_A\})=\phi(\{\boldsymbol{r}_i\};\{\boldsymbol{R}_A\})\chi(\{\boldsymbol{R}_A\}) \tag{3.4} $$
式中:$\phi(\{\boldsymbol{r}_i\};\{\boldsymbol{R}_A\})$ 代表在 $\{\boldsymbol{R}_A\}$ 构型下的电子运动波函数;$\chi(\{\boldsymbol{R}_A\})$ 代表相应的核运动波函数。假设电子运动波函数 $\phi(\{\boldsymbol{r}_i\};\{\boldsymbol{R}_A\})$ 满足:
$$ \underbrace{\left(-\sum_{i=1}^{N}\frac{1}{2}\boldsymbol{\nabla}_i^{2}+\sum_{i=1}^{N}\sum_{j\gt i}^{N}\frac{1}{r_{ij}}-\sum_{i=1}^{N}\sum_{A=1}^{M}\frac{Z_A}{r_{iA}}\right)}_{{\hat H}_{\mathrm{elec}}}\phi(\{\boldsymbol{r}_i\};\{\boldsymbol{R}_A\})={E}_{\mathrm{elec}}(\{\boldsymbol{R}_A\})\phi(\{\boldsymbol{r}_i\};\{\boldsymbol{R}_A\}) \tag{3.5} $$
式(3.5)中,$\hat H_{\mathrm{elec}}$ 是在给定核构型下作用于电子波函数的电子哈密顿算符,含电子动能、电子间排斥及电子与核的吸引,不含核动能与核间排斥;$E_{\mathrm{elec}}(\{\boldsymbol R_A\})$ 是该核构型下的电子能量本征值。
下面研究核运动波函数 $\chi(\{\boldsymbol{R}_A\})$ 需要满足什么样的条件,才能使得方程(3.1)成立。
$$ \begin{aligned} {\hat H}\Phi&=\left[-\sum_{A=1}^{M}\frac{\nabla_A^2}{2M_A} +\sum_{A\lt B}\frac{Z_AZ_B}{R_{AB}}+{\hat H}_{\mathrm{elec}}\right](\phi\chi)\\ &=\phi\left[-\sum_{A=1}^{M}\frac{\nabla_A^2}{2M_A} +\sum_{A\lt B}\frac{Z_AZ_B}{R_{AB}}+{E}_{\mathrm{elec}}(\boldsymbol R)\right]\chi\\ &\quad-\sum_{A=1}^{M}\frac{1}{2M_A} \left[2(\nabla_A\phi)\cdot(\nabla_A\chi)+\chi\nabla_A^2\phi\right], \qquad\Phi=\phi(\boldsymbol r;\boldsymbol R)\chi(\boldsymbol R). \end{aligned} \tag{3.6} $$
在绝热近似中,若核坐标对电子态的导数耦合相对于相关能隙足够小,可忽略式(3.6)的交叉导数项;二阶导数项也需另作估计,保留时可产生对角 Born–Oppenheimer 修正。电子态的归一化本身不能使这些项严格为零。在上述近似成立时,核运动波函数满足
$$ \left[-\sum_{A=1}^{M}\frac{1}{2M_A}\boldsymbol{\nabla}_A^{2}+\sum_{A=1}^{M}\sum_{B\gt A}^{M}\frac{Z_AZ_B}{R_{AB}}+{E}_{\mathrm{elec}}(\{\boldsymbol{R}_A\})\right]\chi(\{\boldsymbol{R}_A\})={E}\chi(\{\boldsymbol{R}_A\}) \tag{3.7} $$
则在玻恩-奥本海默近似下,可以得到体系的总波函数为
$$ \Phi(\{\boldsymbol{r}_i\},\{\boldsymbol{R}_A\})=\phi(\{\boldsymbol{r}_i\};\{\boldsymbol{R}_A\})\chi(\{\boldsymbol{R}_A\}) \tag{3.8} $$
由上面的推导可以看出,体系波函数中的电子自由度和核自由度可以被有效分离。在求解过程中,首先需要求得某个固定核构型下的电子基态,然后将电子能量的本征值(核构型的泛函)作为参数,来求解核运动的本征值问题。以下主要讨论电子自由度,也就是电子波函数 $\phi(\{\boldsymbol{r}_i\};\{\boldsymbol{R}_A\})$ 的求解问题。
平均场的概念
进行玻恩-奥本海默近似后,所需要求解的是固定核构型下电子的基态波函数(即电子运动波函数),它满足如下的薛定谔方程:
$$ \underbrace{\left(-\sum_{i=1}^{N}\frac{1}{2}\boldsymbol{\nabla}_i^{2}+\sum_{i=1}^{N}\sum_{j\gt i}^{N}\frac{1}{r_{ij}}-\sum_{i=1}^{N}\sum_{A=1}^{M}\frac{Z_A}{r_{iA}}\right)}_{{\hat H}_{\mathrm{elec}}}\phi(\{\boldsymbol{r}_i\};\{\boldsymbol{R}_A\})={E}_{\mathrm{elec}}(\{\boldsymbol{R}_A\})\phi(\{\boldsymbol{r}_i\};\{\boldsymbol{R}_A\}) $$
这仍然是一个相当具有挑战性的问题,这是因为电子哈密顿量 ${\hat H}_{\mathrm{elec}}$ 中的 $1/r_{ij}$ 项,使得我们无法用分离变量的办法求解上述方程。当然,在极端的情况下,也就是假设电子与电子之间的库仑相互作用 $\left(\displaystyle\sum_{i=1}^{N}\sum_{j\gt i}^{N}\frac{1}{r_{ij}}\right)$ 为零时,方程可以写成分离变量的形式进行求解(下面讨论原子核构型固定的情况,因此可省略波函数中的原子核坐标):
$$ \left[\sum_{i=1}^{N}\left(-\frac{1}{2}\boldsymbol{\nabla}_i^{2}+V_{\mathrm{ion}}\right)\right]\phi(\{\boldsymbol{r}_i\})={E}_{\mathrm{elec}}\phi(\{\boldsymbol{r}_i\}) \tag{3.9} $$
式中
$$ V_{\mathrm{ion}}(\boldsymbol r_i) =-\sum_{A=1}^{M}\frac{Z_A}{|\boldsymbol r_i-\boldsymbol R_A|} \tag{3.10} $$
但是在真实的物理体系中,电子与电子之间的相互作用是相当强的,至少和电子与核之间的相互作用在同一个数量级。实际上,为了达到分离变量的目的,也并不需要完全忽略电子间相互作用。这是因为我们可以用一个局域的势场来近似地描述其他电子所产生的作用,这个势场和由核产生的势场叠加所形成的“有效势”,就是独立电子空间运动所处的“平均场”(mean field)。因此在平均场近似下,方程(3.9)应当改写成
$$ \left[\sum_{i=1}^{N}\left(-\frac{1}{2}\boldsymbol{\nabla}_i^{2}+V_{\mathrm{eff}}\right)\right]\phi(\{\boldsymbol{r}_i\})={E}_{\mathrm{elec}}\phi(\{\boldsymbol{r}_i\}) \tag{3.11} $$
需要特别指出的是,虽然我们假设电子之间的运动是独立的(独立电子近似),但是这并不意味着求体系的总能等于各个电子能量的简单求和,也不意味着各电子空间分布概率完全不相关。这是因为和其他的微观粒子一样,电子也是不可区分的全同粒子,其波函数必须满足对称(玻色子)或者反对称(费米子)的量子力学要求。虽然没有经典的相互作用项,但是对波函数的交换对称性要求会对总能计算或者空间相对分布概率产生间接影响。总体上来说,交换对称性使相同量子态的玻色子呈现统计聚集,而反对称性使同自旋费米子呈现交换空穴;这不是由统计本身产生的相互作用力。费米子波函数的量子力学要求使同自旋费米子在空间上呈现统计反聚集。下面的例子说明这一点。
假设有两个自由粒子,分别处在自旋向上的动量本征态 $|\,\boldsymbol{k}_1\rangle$ 和 $|\,\boldsymbol{k}_2\rangle$,有
$$ \phi_1(x_1)=\frac{1}{(2\pi)^{3/2}}\mathrm{e}^{\mathrm{i}\boldsymbol{k}_1r_1}\alpha(s_1) $$
$$ \phi_2(x_2)=\frac{1}{(2\pi)^{3/2}}\mathrm{e}^{\mathrm{i}\boldsymbol{k}_2r_2}\alpha(s_2) $$
如果这两个粒子为费米子,则体系符合交换反对称性要求的总波函数可以写为
$$ \begin{aligned} \phi(x_1,x_2)&=\frac{1}{\sqrt{2}}\frac{1}{(2\pi)^{3}}\begin{vmatrix}\mathrm{e}^{\mathrm{i}\boldsymbol{k}_1r_1}\alpha(s_1)&\mathrm{e}^{\mathrm{i}\boldsymbol{k}_2r_1}\alpha(s_1)\\\mathrm{e}^{\mathrm{i}\boldsymbol{k}_1r_2}\alpha(s_2)&\mathrm{e}^{\mathrm{i}\boldsymbol{k}_2r_2}\alpha(s_2)\end{vmatrix}\\ &=\frac{1}{\sqrt{2}}\frac{1}{(2\pi)^{3}}\big[\mathrm{e}^{\mathrm{i}(\boldsymbol{KR}+\boldsymbol{kr})}-\mathrm{e}^{\mathrm{i}(\boldsymbol{KR}-\boldsymbol{kr})}\big]\alpha(s_1)\alpha(s_2)\\ &=\frac{\mathrm{i}\sqrt{2}}{(2\pi)^{3/2}}\sin(\boldsymbol{kr})\left[\frac{1}{(2\pi)^{3/2}}\mathrm{e}^{\mathrm{i}\boldsymbol{KR}}\right]\alpha(s_1)\alpha(s_2) \end{aligned} \tag{3.12} $$
式中:$\boldsymbol{K}=\boldsymbol{k}_1+\boldsymbol{k}_2$;$\boldsymbol{k}=\boldsymbol{k}_1-\boldsymbol{k}_2$;$\boldsymbol{R}=\dfrac{\boldsymbol{r}_1+\boldsymbol{r}_2}{2}$;$\boldsymbol{r}=\dfrac{\boldsymbol{r}_1-\boldsymbol{r}_2}{2}$;$\dfrac{1}{(2\pi)^{3/2}}\mathrm{e}^{\mathrm{i}\boldsymbol{KR}}$ 表示质心的运动,不影响两粒子的相对分布概率。对相对位矢方向做球面平均,可得到相对位置的角平均概率密度;它尚未乘上径向体积因子 $4\pi r^2\,\mathrm dr$,因此不是距离落在 $(r,r+\mathrm dr)$ 内的概率:
$$ P(r)=\frac{1}{4\pi}\iint\left|\frac{\mathrm{i}\sqrt{2}}{(2\pi)^{3/2}}\sin(kr\cos\theta)\right|^{2}\sin\theta\,\mathrm{d}\theta\,\mathrm{d}\phi=\frac{1}{(2\pi)^{3}}\left[1-\frac{\sin(2kr)}{2kr}\right] \tag{3.13} $$
可见,当 $r\to0$ 时,$P(r)\to0$;对于同自旋电子,反对称性使它们在同一点的联合密度为零。利用同样的推导过程,并将第二个电子的自旋方向改为 $\beta(s)$,可以证明 $P(r)\equiv\dfrac{1}{(2\pi)^{3}}$,表明在此无相互作用双电子例子中,异自旋电子没有交换空穴。
电子的空间轨道与自旋轨道
在平均场近似下,每个独立电子满足薛定谔方程
$$ \left(-\frac{1}{2}\boldsymbol{\nabla}_i^{2}+V_{\mathrm{eff}}\right)\phi_{n\sigma}(\boldsymbol{r}_i)=\varepsilon_{n\sigma}\phi_{n\sigma}(\boldsymbol{r}_i) \tag{3.14} $$
式中 $n$ 是单电子轨道指标,并不直接表示体系的第 $n$ 个激发态;$\varepsilon_{n\sigma}$ 为相应轨道本征值,$\sigma$ 标记自旋向上($\alpha$)或自旋向下($\beta$)。若平均场依赖自旋,有效势也须写成 $V_{\mathrm{eff}}^\sigma$。自旋轨道(spin orbital)可以分解为空间部分与自旋部分的直积,即
$$ \phi_{n\sigma}(\boldsymbol{r}_i)=\xi_n(\boldsymbol{r}_i)\sigma(s_i) \tag{3.15} $$
Hartree-Fock 方法
引入 Hartree-Fock 近似后,可以有效地把方程(3.5)等价地转化为 $N$ 个互相独立的可分离变量的方程,从而使得数值求解电子基态波函数成为可能。值得一提的是,类似于量子蒙特卡罗的求解方法则不需引入独立电子的概念,但是其计算量是可分离变量的 Hartree-Fock 方法计算量的上千倍甚至上亿倍。
首先考虑 Hartree-Fock 近似背后所代表的物理意义。电子是费米子的一种,因此其对波函数的要求首先是反对称性。如果体系中有 $N$ 个电子,一共有 $K$ 个可供占据的自旋轨道,则普遍来说,体系的基态(或者激发态)的波函数可以用自旋轨道所组成的反对称的 Slater 行列式进行展开:
$$ \begin{aligned} \phi\langle\boldsymbol{x}_1,\boldsymbol{x}_2,\cdots,\boldsymbol{x}_N\rangle={}&\frac{C_1}{\sqrt{N!}}\begin{vmatrix}\xi_i(\boldsymbol{x}_1)&\xi_j(\boldsymbol{x}_1)&\cdots&\xi_k(\boldsymbol{x}_1)\\\xi_i(\boldsymbol{x}_2)&\xi_j(\boldsymbol{x}_2)&\cdots&\xi_k(\boldsymbol{x}_2)\\\vdots&\vdots&&\vdots\\\xi_i(\boldsymbol{x}_N)&\xi_j(\boldsymbol{x}_N)&\cdots&\xi_k(\boldsymbol{x}_N)\end{vmatrix}+\frac{C_2}{\sqrt{N!}}\begin{vmatrix}\xi_{i'}(\boldsymbol{x}_1)&\xi_j(\boldsymbol{x}_1)&\cdots&\xi_k(\boldsymbol{x}_1)\\\xi_{i'}(\boldsymbol{x}_2)&\xi_j(\boldsymbol{x}_2)&\cdots&\xi_k(\boldsymbol{x}_2)\\\vdots&\vdots&&\vdots\\\xi_{i'}(\boldsymbol{x}_N)&\xi_j(\boldsymbol{x}_N)&\cdots&\xi_k(\boldsymbol{x}_N)\end{vmatrix}\\ &+\cdots+\frac{C'}{\sqrt{N!}}\begin{vmatrix}\xi_{i'}(\boldsymbol{x}_1)&\xi_{j'}(\boldsymbol{x}_1)&\cdots&\xi_{k'}(\boldsymbol{x}_1)\\\xi_{i'}(\boldsymbol{x}_2)&\xi_{j'}(\boldsymbol{x}_2)&\cdots&\xi_{k'}(\boldsymbol{x}_2)\\\vdots&\vdots&&\vdots\\\xi_{i'}(\boldsymbol{x}_N)&\xi_{j'}(\boldsymbol{x}_N)&\cdots&\xi_{k'}(\boldsymbol{x}_N)\end{vmatrix} \end{aligned} \tag{3.16} $$
方程(3.16)中的每一个行列式都称为一个组态。要精确地展开体系的波函数,有可能需要用到上千个组态。当然所有这些组态中,和体系基态最接近的应该是从 $K$ 个轨道中挑出 $N$ 个能量最低的自旋轨道所组成的行列式。Hartree-Fock 近似从本质上来说就是用 $N$ 个能量最低轨道所组成的单行列式来近似体系的真实波函数。我们用 $|\,\xi_i(1)\xi_j(2)\cdots\xi_k(N)\rangle_{\mathrm{S}}$ 来区别 Slater 形式波函数和普通的右矢 $|\,\xi_i(1)\xi_j(2)\cdots\xi_k(N)\rangle$,有
$$ \begin{aligned} &\phi(\boldsymbol{x}_1,\boldsymbol{x}_2,\cdots,\boldsymbol{x}_N)\simeq\phi^{0}(\boldsymbol{x}_1,\boldsymbol{x}_2,\cdots,\boldsymbol{x}_N)\\ &=|\,\xi_i(1)\xi_j(2)\cdots\xi_k(N)\,|\rangle_{\mathrm{S}}\\ &=(N!)^{-1/2}\sum_{n=1}^{N!}(-1)^{P_n}{\hat P_n^{(\mathrm{perm})}}\{\xi_i(1)\xi_j(2)\cdots\xi_k(N)\}=(N!)^{-1/2}\begin{vmatrix}\xi_i(\boldsymbol{x}_1)&\xi_j(\boldsymbol{x}_1)&\cdots&\xi_k(\boldsymbol{x}_1)\\\xi_i(\boldsymbol{x}_2)&\xi_j(\boldsymbol{x}_2)&\cdots&\xi_k(\boldsymbol{x}_2)\\\vdots&\vdots&&\vdots\\\xi_i(\boldsymbol{x}_N)&\xi_j(\boldsymbol{x}_N)&\cdots&\xi_k(\boldsymbol{x}_N)\end{vmatrix} \end{aligned} \tag{3.17} $$
式中:$n$ 枚举 $N!$ 种电子标号的置换;$\hat P_n^{(\mathrm{perm})}$ 按第 $n$ 种置换重新排列轨道乘积中的电子标号,$(-1)^{P_n}$ 则按置换的奇偶性取 $+1$ 或 $-1$,使求和结果在交换任意两个电子时变号。下面首先来推导当基态波函数为单个行列式的时候,体系的总能和单电子各个能级之间的关系。在给出普适的表达式之前,先来看一下双电子体系的例子。
归一化、反对称的双电子波基态函数在 Hartree-Fock 近似下可以通过最低占据的单电子轨道反对称化得到:
$$ \phi^{0}(\boldsymbol{x}_1,\boldsymbol{x}_2)=|\,\xi_i(1)\xi_j(2)\rangle_{\mathrm{S}}=\frac{1}{\sqrt{2}}\begin{vmatrix}\xi_i(\boldsymbol{x}_1)&\xi_j(\boldsymbol{x}_1)\\\xi_i(\boldsymbol{x}_2)&\xi_j(\boldsymbol{x}_2)\end{vmatrix}=\frac{1}{\sqrt{2}}\big[\xi_i(\boldsymbol{x}_1)\xi_j(\boldsymbol{x}_2)-\xi_i(\boldsymbol{x}_2)\xi_j(\boldsymbol{x}_1)\big] \tag{3.18} $$
同时,为了计算方便,将电子哈密顿量(原子单位制)分解成单电子部分和双电子部分,即
$$ \begin{aligned} {\hat H}_{\mathrm{elec}}&=-\sum_{i=1}^{N}\frac{1}{2}\boldsymbol{\nabla}_i^{2}+\sum_{i=1}^{N}\sum_{j\gt i}^{N}\frac{1}{r_{ij}}-\sum_{i=1}^{N}\sum_{A=1}^{M}\frac{Z_A}{r_{iA}}=\sum_{i=1}^{N}\left[-\frac{1}{2}\boldsymbol{\nabla}_i^{2}-\sum_{A=1}^{M}\frac{Z_A}{r_{iA}}\right]+\sum_{i=1}^{N}\sum_{j\gt i}^{N}\frac{1}{r_{ij}}\\ &=\sum_{i=1}^{N}h(i)+\sum_{i\lt j}v(i,j)={\hat O}_1+{\hat O}_2 \end{aligned} \tag{3.19} $$
这里 $\hat O_1=\sum_i h(i)$ 是一体算符的总和:$h(i)$ 只作用于第 $i$ 个电子,包含其动能和核的吸引作用;$\hat O_2=\sum_{i\lt j}v(i,j)$ 是二体算符的总和:$v(i,j)=1/r_{ij}$ 给出第 $i$、$j$ 个电子之间的库仑排斥。因此 $\hat H_{\mathrm{elec}}=\hat O_1+\hat O_2$。
接下来根据定义计算体系的基态能量,有
$$ E=\langle\phi^{0}(\boldsymbol{x}_1,\boldsymbol{x}_2)\,|\,{\hat H}_{\mathrm{elec}}\,|\,\phi^{0}(\boldsymbol{x}_1,\boldsymbol{x}_2)\rangle=\langle\phi^{0}(\boldsymbol{x}_1,\boldsymbol{x}_2)\,|\,{\hat O}_1+{\hat O}_2\,|\,\phi^{0}(\boldsymbol{x}_1,\boldsymbol{x}_2)\rangle \tag{3.20} $$
首先计算单电子部分 ${\hat O}_1$ 的贡献。由于自旋轨道之间的正交性,对于 ${\hat O}_1$,只有左矢和右矢完全等同时,积分才不等于零,于是有
$$ \langle\phi^{0}(\boldsymbol{x}_1,\boldsymbol{x}_2)\,|\,{\hat O}_1\,|\,\phi^{0}(\boldsymbol{x}_1,\boldsymbol{x}_2)\rangle $$
$$ \begin{aligned} &=\sum_{i=1}^{2}\langle\phi^{0}(\boldsymbol{x}_1,\boldsymbol{x}_2)\,|\,h(i)\,|\,\phi^{0}(\boldsymbol{x}_1,\boldsymbol{x}_2)\rangle_{\mathrm{S}}\\ &=2\langle\phi^{0}(\boldsymbol{x}_1,\boldsymbol{x}_2)\,|\,h(1)\,|\,\phi^{0}(\boldsymbol{x}_1,\boldsymbol{x}_2)\rangle\\ &=\langle\xi_i(\boldsymbol{x}_1)\xi_j(\boldsymbol{x}_2)-\xi_i(\boldsymbol{x}_2)\xi_j(\boldsymbol{x}_1)\,|\,h(1)\,|\,\xi_i(\boldsymbol{x}_1)\xi_j(\boldsymbol{x}_2)-\xi_i(\boldsymbol{x}_2)\xi_j(\boldsymbol{x}_1)\rangle\\ &=\langle\xi_i(\boldsymbol{x}_1)\,|\,h(1)\,|\,\xi_i(\boldsymbol{x}_1)\rangle+\langle\xi_j(\boldsymbol{x}_1)\,|\,h(1)\,|\,\xi_j(\boldsymbol{x}_1)\rangle\\ &=\sum_{i=1}^{2}\langle i\,|\,h\,|\,i\rangle \end{aligned} \tag{3.21} $$
以上双电子的情况非常容易推广到多电子。对于一个多电子的体系,在 Hartree-Fock 近似下,基态的电子波函数用 $N$ 个能量最低占据轨道的反对称波函数近似,而在此近似下,基态能量可以表示为单电子积分和双电子积分的形式,有
$$ \begin{aligned} E_0=\langle\phi^{0}\,|\,{\hat H}_{\mathrm{elec}}\,|\,\phi^{0}\rangle&=\sum_{i=1}^{N}\langle i\,|\,h\,|\,i\rangle+\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\big(\langle ij\,|\,ij\rangle-\langle ij\,|\,ji\rangle\big)\\ &=\sum_{i=1}^{N}\langle i\,|\,h\,|\,i\rangle+\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\langle ij\,\|\,ij\rangle \end{aligned} \tag{3.22} $$
式中
$$ \begin{aligned} \langle ij\,\|\,ij\rangle&=\langle\xi_i\xi_j\,|\,\xi_i\xi_j\rangle-\langle\xi_i\xi_j\,|\,\xi_j\xi_i\rangle\\ &=\int\mathrm{d}\boldsymbol{x}_1\,\mathrm{d}\boldsymbol{x}_2\,\xi_i^{*}(\boldsymbol{x}_1)\xi_j^{*}(\boldsymbol{x}_2)\frac{1}{r_{12}}\big[\xi_i(\boldsymbol{x}_1)\xi_j(\boldsymbol{x}_2)-\xi_j(\boldsymbol{x}_1)\xi_i(\boldsymbol{x}_2)\big] \end{aligned} \tag{3.23} $$
可以看到,在 Hartree-Fock 近似下,体系能量的表达式的物理意义非常明显。单电子项表达的是电子的动能项和电子与核之间的库仑相互作用。双电子项中的 $\langle\xi_i\xi_j\,|\,\xi_i\xi_j\rangle$ 可以根据电子密度的定义改写为 $\dfrac{\rho_i(x_1)\rho_j(x_2)}{r_{12}}$,表达电子之间的静电库仑斥能。双电子项中的 $\langle\xi_i\xi_j\,|\,\xi_j\xi_i\rangle$ 表达的则是相同自旋电子之间的交换作用,其源于 Slater 行列式的波函数中两个自旋相同的电子之间的交换关联作用。
需要注意的是,在部分化学书中,用到了另外一种不同但是等价的积分简写方式,即
$$ \begin{aligned} \langle ij|kl\rangle &=\iint\xi_i^*(x_1)\xi_j^*(x_2) \frac{\xi_k(x_1)\xi_l(x_2)}{r_{12}}\,\mathrm dx_1\,\mathrm dx_2\\ &=[ik|jl],\qquad r_{12}=|\boldsymbol r_1-\boldsymbol r_2|. \end{aligned} \tag{3.24} $$
Hartree-Fock 近似下的单电子自洽场方程
在 3.1.4 节中我们讨论了在 Hartree-Fock 近似下构造基态波函数的方法,以及体系的总能量和各个单电子轨道的单电子积分和双电子积分之间的关系。但是,如何构造、求解 1,2,…,$N$ 个被占据的单电子轨道方程呢?在这一节中,我们从 Hartree-Fock 的总能表达式出发,利用变分原理,推导 Hartree-Fock 近似下单电子轨道所满足的方程。
根据变分原理可知,任意归一化的试探波函数都满足
$$ \langle\hat{\Phi}\,|\,{\hat H}\,|\,\hat{\Phi}\rangle\geqslant{E}_0 \tag{3.25} $$
我们可以利用变分原理来有效地求解薛定谔方程的最佳近似解。也就是构造一系列含参数的试探波函数,然后通过变化参数,使得在这组参数下哈密顿量的期望值最小。这组参数所对应的波函数就是在相应子空间中,薛定谔方程的最佳近似解。
根据拉格朗日乘子法,需要在保持各个自旋轨道正交的情况下,变化各个轨道,使总能达到最小。也说是在 $\langle\xi_i\,|\,\xi_j\rangle=\delta_{ij}$ 的情况下,找到一组自旋轨道 $\xi$,使得
$$ \delta E_0=0 \tag{3.26} $$
其中
$$ E_0=\sum_{i=1}^{N}\langle i\,|\,h\,|\,i\rangle+\frac{1}{2}\sum_{i,j=1}^{N}\big(\langle ij\,|\,ij\rangle-\langle ij\,|\,ji\rangle\big) \tag{3.27} $$
在给出具体推导过程前,我们先给出最终结果,并对其物理意义进行分析。最终得到 Hartree-Fock 的单电子轨道满足的自洽场(self-consistent field,SCF)方程为
$$ \begin{aligned} &h(x_1)\xi_i(x_1)+\sum_{j\neq i}\left[\int\frac{\mathrm{d}x_2\,|\,\xi_j(x_2)\,|^{2}}{r_{12}}\right]\xi_i(x_1)-\sum_{j\neq i}\left[\int\frac{\mathrm{d}x_2\,\xi_j^{*}(x_2)\xi_i(x_2)}{r_{12}}\right]\xi_j(x_1)\\ &=\varepsilon_i\xi_i(x_1) \end{aligned} \tag{3.28} $$
式中,$\varepsilon_i$ 为第 $i$ 个自洽轨道的轨道能量。通常根据物理意义,引入库仑算符和交换算符,则
$$ {\hat J}_j(x_1)\xi_i(x_1)=\left(\int\mathrm{d}\boldsymbol{x}_2\,\xi_j^{*}(x_2)\frac{1}{r_{12}}\xi_j(x_2)\right)\xi_i(x_1) \tag{3.29} $$
$$ ({\hat K}_j\xi_i)(x_1) =\xi_j(x_1)\int\frac{\xi_j^*(x_2)\xi_i(x_2)}{r_{12}}\,\mathrm dx_2 \tag{3.30} $$
对于待作用的轨道 $\xi_i$,$\hat J_j$ 将它乘以由占据轨道 $\xi_j$ 的密度产生的平均库仑势,是局域乘法算符;$\hat K_j$ 则先对 $\xi_i$ 与 $\xi_j$ 在另一电子坐标上的乘积积分,再乘上 $\xi_j(x_1)$,因此是非局域交换算符。两个算符均由占据轨道 $\xi_j$ 决定。
利用这两个算符,Hartree-Fock 的单电子自洽场方程可以简洁地表示为
$$ \left(h(x_1)+\sum_{j\neq i}^{N}{\hat J}_j(x_1)-\sum_{j\neq i}^{N}{\hat K}_j(x_1)\right)\xi_i(x_1)=\varepsilon_i\xi_i(x_1) \tag{3.31} $$
上述方程中,对于不同的轨道,库仑算符和交换算符分别需要去掉 $j=i$ 的项,因此在形式表达上不方便。由于 $j=i$ 时库仑算符和交换算符相互抵消,还可以去除求和下标中 $j\neq i$ 的限制,即
$$ \left(h(x_1)+\sum_{j=1}^{N}{\hat J}_j(x_1)-\sum_{j=1}^{N}{\hat K}_j(x_1)\right)\xi_i(x_1)=\varepsilon_i\xi_i(x_1) \tag{3.32} $$
如果定义 Fock 算符为
$$ {\hat F}(x_1)=h(x_1)+\sum_{j=1}^{N}\big({\hat J}_j(x_1)-{\hat K}_j(x_1)\big) \tag{3.33} $$
$\hat F$ 是由单电子算符 $h$、平均库仑作用与交换作用组成的有效单电子算符。由于 $\hat J_j$ 和 $\hat K_j$ 依赖于占据轨道,$\hat F$ 也随所求轨道改变,故需自洽求解。
则正则 Hartree-Fock 方程有如下非常简洁的形式:
$$ {\hat F}\,|\,\xi_i(x_1)\rangle=\varepsilon_i\,|\,\xi_i(x_1)\rangle \tag{3.34} $$
下面我们通过变分原理给出正则 Hartree-Fock 方程的推导过程。记带正交约束的变分函数为 $L_{\mathrm{HF}}$。这里的变分函数是自旋轨道,也就是说,总能对正交的自旋轨道的变分为零:
$$ \xi_i\to\xi_i+\delta\xi_i\quad(i=1,2,\cdots,N) $$
$$ \delta{L_{\mathrm{HF}}}=\delta E_0-\delta\left[\sum_{i,j=1}^{N}\varepsilon_{ji}\big(\langle i\,|\,j\rangle-\delta_{ij}\big)\right]=\delta E_0-\sum_{i,j=1}^{N}\varepsilon_{ji}\,\delta\langle i\,|\,j\rangle=0 \tag{3.35} $$
根据 $E_0$ 的表达式
$$ E_0=\sum_{i=1}^{N}\langle i\,|\,h\,|\,i\rangle+\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\big(\langle ij\,|\,ij\rangle-\langle ij\,|\,ji\rangle\big) $$
容易得到
$$ \begin{aligned} \delta E_0&=\sum_{i=1}^{N}\delta\langle i\,|\,h\,|\,i\rangle+\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\big(\delta\langle ij\,|\,ij\rangle-\delta\langle ij\,|\,ji\rangle\big)\\ &=\sum_{i=1}^{N}\langle\delta\xi_i\,|\,h\,|\,\xi_i\rangle+\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\big(\langle\delta\xi_i\xi_j\,|\,\xi_i\xi_j\rangle+\langle\xi_i\delta\xi_j\,|\,\xi_i\xi_j\rangle\big)\\ &\quad-\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\big(\langle\delta\xi_i\xi_j\,|\,\xi_j\xi_i\rangle+\langle\xi_i\delta\xi_j\,|\,\xi_j\xi_i\rangle\big)+\mathrm{C.\,C.}\\ &=\sum_{i=1}^{N}\langle\delta\xi_i\,|\,h\,|\,\xi_i\rangle+\sum_{i=1}^{N}\sum_{j=1}^{N}\langle\delta\xi_i\xi_j\,|\,\xi_i\xi_j\rangle-\sum_{i=1}^{N}\sum_{j=1}^{N}\langle\delta\xi_i\xi_j\,|\,\xi_j\xi_i\rangle+\mathrm{C.\,C.} \end{aligned} \tag{3.36} $$
其中 C. C. 代表共轭项。
式(3.35)中第二项的变分为
$$ \sum_{i,j=1}^{N}\varepsilon_{ji}\,\delta\langle i\,|\,j\rangle=\sum_{i,j=1}^{N}\varepsilon_{ji}\langle\delta\xi_i\,|\,\xi_j\rangle+\mathrm{C.\,C.} \tag{3.37} $$
因此总能量对自旋轨道变分为零的条件等价转化为
$$ \begin{aligned} \delta{L_{\mathrm{HF}}}={}&\sum_{i=1}^{N}\langle\delta\xi_i\,|\,h\,|\,\xi_i\rangle+\sum_{i=1}^{N}\sum_{j=1}^{N}\langle\delta\xi_i\xi_j\,|\,\xi_i\xi_j\rangle-\sum_{i=1}^{N}\sum_{j=1}^{N}\langle\delta\xi_i\xi_j\,|\,\xi_j\xi_i\rangle\\ &-\sum_{i,j=1}^{N}\varepsilon_{ji}\langle\delta\xi_i\,|\,\xi_j\rangle+\mathrm{C.\,C.}=0 \end{aligned} \tag{3.38} $$
将方程(3.35)改写后,因为取极值的条件必须对于任意的 $\delta\xi_i$ 均成立,因此括号里的项必须为零,即
$$ \begin{aligned} \delta{L_{\mathrm{HF}}}&=\int\sum_{i=1}^{N}\mathrm{d}\boldsymbol{x}_1\,\delta\xi_i^{*}(x_1)\left[h(x_1)\xi_i(x_1)+\sum_{j=1}^{N}\big({\hat J}_j(x_1)-{\hat K}_j(x_1)\big)\xi_i(x_1)-\sum_{j=1}^{N}\varepsilon_{ji}\xi_j(x_1)\right]+\mathrm{C.\,C.}\\ &=0 \end{aligned} $$
得
$$ \left[h(x_1)+\sum_{j=1}^{N}\big({\hat J}_j(x_1)-{\hat K}_j(x_1)\big)\right]\xi_i(x_1)=\sum_{j=1}^{N}\varepsilon_{ji}\xi_j(x_1) \tag{3.39} $$
因此有
$$ {\hat F}(x_1)\,|\,\xi_i(x_1)\rangle=\sum_{j=1}^{N}\varepsilon_{ji}\,|\,\xi_j(x_1)\rangle \tag{3.40} $$
至此,我们得到了和 Hartree-Fock 方程等价的结果,但是与正则 Hartree-Fock 方程(式(3.34))相比较,还是略有不同。二者之间的差别可以通过幺正变换消除。首先考察当一组自旋轨道通过幺正变换成一组新的自旋轨道时,厄米算符(物理上的可观测量)所对应的期望值如何变化。有
$$ |\,\xi'_i(x_1)\xi'_j(x_2)\cdots\xi'_k(x_N)\rangle_{\mathrm{S}}=|\,\xi_i(x_1)\xi_j(x_2)\cdots\xi_k(x_N)\rangle_{\mathrm{S}}\cdot U $$
即
$$ \begin{vmatrix}\xi'_i(\boldsymbol{x}_1)&\xi'_j(\boldsymbol{x}_1)&\cdots&\xi'_k(\boldsymbol{x}_1)\\\xi'_i(\boldsymbol{x}_2)&\xi'_j(\boldsymbol{x}_2)&\cdots&\xi'_k(\boldsymbol{x}_2)\\\vdots&\vdots&&\vdots\\\xi'_i(\boldsymbol{x}_N)&\xi'_j(\boldsymbol{x}_N)&\cdots&\xi'_k(\boldsymbol{x}_N)\end{vmatrix}=\begin{vmatrix}\xi_i(\boldsymbol{x}_1)&\xi_j(\boldsymbol{x}_1)&\cdots&\xi_k(\boldsymbol{x}_1)\\\xi_i(\boldsymbol{x}_2)&\xi_j(\boldsymbol{x}_2)&\cdots&\xi_k(\boldsymbol{x}_2)\\\vdots&\vdots&&\vdots\\\xi_i(\boldsymbol{x}_N)&\xi_j(\boldsymbol{x}_N)&\cdots&\xi_k(\boldsymbol{x}_N)\end{vmatrix}\begin{vmatrix}U_{11}&U_{12}&\cdots&U_{1N}\\U_{21}&U_{22}&\cdots&U_{2N}\\\vdots&\vdots&&\vdots\\U_{N1}&U_{N2}&\cdots&U_{NN}\end{vmatrix} \tag{3.41} $$
由于幺正变换满足 $U^{*}\cdot U=1$,因此容易得出 $|\det(U)\,|^{2}=1$。也就是说 $\det(U)=\mathrm{e}^{\mathrm{i}\varphi}$。
我们可以得到,任意的厄米算符,包括总能量、动量等的期望值,在自旋轨道幺正变换下均保持不变。也就是说,自旋轨道的确定具有一定的任意性,给定的一组是 Hartree-Fock 方程解的自旋轨道,对其做幺正变换后得到的新的自旋轨道同样是方程的解。下面我们考察如何利用幺正变换,将方程(3.40)等价地转变成正则 Hartree-Fock 方程。
方程(3.40)左端包括库仑算符、交换算符。其中单电子算符并不依赖于自旋轨道。在占据轨道内部作幺正变换时,单个 $\hat J_i$ 或 $\hat K_i$ 一般会改变,但各自对全部占据轨道的求和保持不变。下式中的撇号表示用变换后的轨道构造的算符;这里先验证库仑算符之和,交换算符之和也可类似证明:
$$ \begin{aligned} \sum_i{\hat J}'_i(x_1) &=\sum_{i,j,k}U_{ji}^*U_{ki} \int\frac{\xi_j^*(x_2)\xi_k(x_2)}{r_{12}}\,\mathrm dx_2\\ &=\sum_{j,k}\delta_{jk}\int\frac{\xi_j^*(x_2)\xi_k(x_2)}{r_{12}}\,\mathrm dx_2 =\sum_j{\hat J}_j(x_1),\qquad U^\dagger U=I. \end{aligned} \tag{3.42} $$
因此可知,Fock 算符在自旋轨道的幺正变换下保持不变。进一步容易得到,拉格朗日乘子 $\varepsilon_{ji}$ 满足下式:
$$ \langle\xi_k(x_1)\,|\,{\hat F}(x_1)\,|\,\xi_i(x_1)\rangle=\sum_{j=1}^{N}\varepsilon_{ji}\langle\xi_k(x_1)\,|\,\xi_j(x_1)\rangle=\varepsilon_{ki} \tag{3.43} $$
因此在自旋轨道幺正变换下,有
$$ \begin{aligned} \varepsilon'_{ij}&=\int\mathrm{d}\boldsymbol{x}_1\,\xi'^{\,*}_i(x_1){\hat F}(x_1)\xi'_j(x_1)=\sum_{k,l}U_{ki}^{*}U_{lj}\int\mathrm{d}\boldsymbol{x}_1\,\xi_k^{*}(x_1){\hat F}(x_1)\xi_l(x_1)\\ &=\sum_{k,l}U_{ki}^{*}\varepsilon_{kl}U_{lj} \end{aligned} \tag{3.44} $$
解得
$$ \varepsilon'=U^{*}\varepsilon U $$
由此,可以看到,总是可以通过幺正变换得到一组自旋轨道,在此组自旋轨道下,$\varepsilon$ 成为一个对角矩阵。相应地,Hartree-Fock 方程退化为正则 Hartree-Fock 方程。
由上述讨论可知,如果我们以最低占据轨道构成的 Slater 行列式近似作为基态波函数,可以运用变分的方法得到一组自洽的 Hartree-Fock 方程。其中单电子的哈密顿量主要由三项构成:第一项 $h$ 是单电子算符,表达动能项和核的吸引作用项;第二项 ${\hat J}$ 是库仑斥能项,表示的是所有其他电子的密度分布对该电子的平均斥能,而并没有考虑电子与电子之间相斥对电子间关联函数的影响;第三项 ${\hat K}$ 是交换能项,在经典物理中没有对应,它所表现的是量子力学对费米子波函数的反对称性要求引起的一种关联作用。
Hartree-Fock 单电子波函数的讨论
本节中,我们将详细讨论用于构建 Hartree-Fock 基态波函数的 Slater 行列式的特点,以及处于此态的电子和独立电子之间的运动的差别,并且引入密度泛函中非常重要的费米空穴的概念。
在 3.1.5 节的讨论中,我们已经看到单电子的 Hartree-Fock 方程(式(3.31))中,Fock 算符有单电子项及双电子项。如果不考虑自旋,则可以相应地建立单电子密度分布函数 $\rho(\boldsymbol{r})$ 以及双电子密度分布函数 $\rho(\boldsymbol{r},\boldsymbol{r}')$:
$$ 1=\int|\,\Psi^{0}(\boldsymbol{r},\boldsymbol{r}_2,\boldsymbol{r}_3,\cdots,\boldsymbol{r}_N)\,|^{2}\,\mathrm{d}\boldsymbol{r}\,\mathrm{d}\boldsymbol{r}_2\cdots\mathrm{d}\boldsymbol{r}_N \tag{3.45} $$
$$ \rho(\boldsymbol{r})=N\int\cdots\int|\,\Psi^{0}(\boldsymbol{r},\boldsymbol{r}_2,\boldsymbol{r}_3,\cdots,\boldsymbol{r}_N)\,|^{2}\,\mathrm{d}\boldsymbol{r}_2\,\mathrm{d}\boldsymbol{r}_3\cdots\mathrm{d}\boldsymbol{r}_N \tag{3.46} $$
$$ \rho(\boldsymbol{r},\boldsymbol{r}')=N(N-1)\int\cdots\int|\,\Psi^{0}(\boldsymbol{r},\boldsymbol{r}',\boldsymbol{r}_3,\cdots,\boldsymbol{r}_N)\,|^{2}\,\mathrm{d}\boldsymbol{r}_3\,\mathrm{d}\boldsymbol{r}_4\cdots\mathrm{d}\boldsymbol{r}_N \tag{3.47} $$
式中:$\Psi^{0}$ 为体系的多体波函数;$\rho(\boldsymbol{r})$ 代表总数为 $N$ 的电子气中,在 $\boldsymbol{r}$ 的终点处发现电子的概率;系数 $N$ 表示体系中有 $N$ 个全同粒子,而每个粒子在空间 $\mathrm{d}\boldsymbol{r}$ 出现的概率相同;$\rho(\boldsymbol{r},\boldsymbol{r}')$ 表示在 $\boldsymbol{r}$ 的终点处发现一个电子的同时,在 $\boldsymbol{r}'$ 的终点处发现另一个电子的概率。式(3.47)前面的系数为 $N(N-1)$,和从 $N$ 个全同粒子中选取两个粒子放到空间两个位置的排列数相同。从上面的表达式还可以知道,电荷分布密度对全空间的积分等于电子总数,即
$$ \int\rho(\boldsymbol{r})\,\mathrm{d}\boldsymbol{r}=N \tag{3.48} $$
对于经典粒子,$\rho^{0}(\boldsymbol{r},\boldsymbol{r}')$ 为两个 $\rho(\boldsymbol{r})$ 的积,即
$$ \rho^{0}(\boldsymbol{r},\boldsymbol{r}')=\rho(\boldsymbol{r})\rho(\boldsymbol{r}') $$
而如果计入自能项(也就是电子不能和自己发生作用),则有
$$ \rho(\boldsymbol{r},\boldsymbol{r}')=\frac{N-1}{N}\rho(\boldsymbol{r})\rho(\boldsymbol{r}') $$
该式明显有别于式(3.47)。虽然经典粒子和独立电子都在势场中独立运动,但仍存在以下不同之处:① 电子遵循费米统计,基于泡利不相容原理,自旋相同的电子在空间上彼此疏离,若已有一个电子在 $\boldsymbol{r}$ 的终点,那么显然在 $\boldsymbol{r}'$ 的终点处发现另一个相同自旋态电子的概率比经典统计的要低;② 电子和电子之间由于带电,存在较强的库仑斥能,这种斥能同样会导致电子彼此疏离,因此,每个电子在它自身周围都存在一个低密度区,称为费米空穴或者交换关联空穴。在费米空穴处电子密度 $\rho_{\mathrm{xc}}(\boldsymbol{r},\boldsymbol{r}')$ 满足
$$ \rho(\boldsymbol{r},\boldsymbol{r}')=\rho(\boldsymbol{r})\rho(\boldsymbol{r}')+\rho(\boldsymbol{r})\rho_{\mathrm{xc}}(\boldsymbol{r},\boldsymbol{r}') \tag{3.49} $$
在 Hartree-Fock 近似中,$\Psi^{0}$ 近似为 $\phi^{0}$,即在轨道 $\xi_i$ 彼此正交的条件下,可以由方程(3.17)解析地给出 $\rho(\boldsymbol{r})$,$\rho(\boldsymbol{r},\boldsymbol{r}')$ 及 $\rho_{\mathrm{xc}}(\boldsymbol{r},\boldsymbol{r}')$。$\rho(\boldsymbol{r})$ 比较简单,有
$$ \rho(\boldsymbol{r})=\sum_{i}\xi_i^{*}(\boldsymbol{r})\xi_i(\boldsymbol{r}) \tag{3.50} $$
为求得 $\rho(\boldsymbol{r},\boldsymbol{r}')$,首先将 $\phi^{0}$ 按第 $i$ 列展开:
$$ \begin{aligned} \phi^{0}&=\frac{1}{\sqrt{N!}}\begin{vmatrix}\xi_1(\boldsymbol{r})&\xi_2(\boldsymbol{r})&\cdots&\xi_N(\boldsymbol{r})\\\xi_1(\boldsymbol{r}')&\xi_2(\boldsymbol{r}')&\cdots&\xi_N(\boldsymbol{r}')\\\vdots&\vdots&&\vdots\\\xi_1(\boldsymbol{r}_N)&\xi_2(\boldsymbol{r}_N)&\cdots&\xi_N(\boldsymbol{r}_N)\end{vmatrix}\\ &=\frac{1}{\sqrt{N!}}\sum_{i}\xi_i(\boldsymbol{r})(-1)^{i+1}\begin{vmatrix}\xi_1(\boldsymbol{r}')&\cdots&\xi_{i-1}(\boldsymbol{r}')&\xi_{i+1}(\boldsymbol{r}')&\cdots&\xi_N(\boldsymbol{r}')\\\xi_1(\boldsymbol{r}_3)&\cdots&\xi_{i-1}(\boldsymbol{r}_3)&\xi_{i+1}(\boldsymbol{r}_3)&\cdots&\xi_N(\boldsymbol{r}_3)\\\vdots&&\vdots&\vdots&&\vdots\\\xi_1(\boldsymbol{r}_N)&\cdots&\xi_{i-1}(\boldsymbol{r}_N)&\xi_{i+1}(\boldsymbol{r}_N)&\cdots&\xi_N(\boldsymbol{r}_N)\end{vmatrix} \end{aligned} \tag{3.51} $$
进一步将式(3.51)按第 $j$ 列展开:
$$ \begin{aligned} \phi^{0}={}&\frac{1}{\sqrt{N!}}\sum_{i,j\neq i}\xi_i(\boldsymbol{r})\xi_j(\boldsymbol{r}')(-1)^{C_{i,j}}\\ &\times\begin{vmatrix}\xi_1(\boldsymbol{r}_3)&\cdots&\xi_{i-1}(\boldsymbol{r}_3)&\xi_{i+1}(\boldsymbol{r}_3)&\cdots&\xi_{j-1}(\boldsymbol{r}_3)&\xi_{j+1}(\boldsymbol{r}_3)&\cdots&\xi_N(\boldsymbol{r}_3)\\\xi_1(\boldsymbol{r}_4)&\cdots&\xi_{i-1}(\boldsymbol{r}_4)&\xi_{i+1}(\boldsymbol{r}_4)&\cdots&\xi_{j-1}(\boldsymbol{r}_4)&\xi_{j+1}(\boldsymbol{r}_4)&\cdots&\xi_N(\boldsymbol{r}_4)\\\vdots&&\vdots&\vdots&&\vdots&\vdots&&\vdots\\\xi_1(\boldsymbol{r}_N)&\cdots&\xi_{i-1}(\boldsymbol{r}_N)&\xi_{i+1}(\boldsymbol{r}_N)&\cdots&\xi_{j-1}(\boldsymbol{r}_N)&\xi_{j+1}(\boldsymbol{r}_N)&\cdots&\xi_N(\boldsymbol{r}_N)\end{vmatrix} \end{aligned} \tag{3.52} $$
式中
$$ C_{i,j}=i+j-1+\frac{\mathrm{sgn}[i-j]+1}{2} \tag{3.53} $$
$\mathrm{sgn}[i-j]$ 为 $i-j$ 的符号函数。
将式(3.52)代入方程(3.47),由于 $\langle\xi_i\,|\,\xi_j\rangle=\delta_{ij}$,因此式(3.52)最后一行的 $N-2$ 阶行列式相乘,只有 $(N-2)!$ 个对角项(即包含相同 $\boldsymbol{r}_l$ 的 $\xi$ 项下标也相同)等于 1,而其他各项因为均包含至少一个形如 $\displaystyle\int\xi_k^{*}(\boldsymbol{r}_l)\xi_{m\neq k}^{*}(\boldsymbol{r}_l)\,\mathrm{d}\boldsymbol{r}_l$ 的项而等于零。因此有
$$ \begin{aligned} \rho_2(\boldsymbol r,\boldsymbol r') &=\rho(\boldsymbol r)\rho(\boldsymbol r') -\sum_\sigma\left|\gamma_\sigma(\boldsymbol r,\boldsymbol r')\right|^2,\\ \gamma_\sigma(\boldsymbol r,\boldsymbol r') &=\sum_i\xi_i^\sigma(\boldsymbol r)\xi_i^{\sigma*}(\boldsymbol r') \end{aligned} \tag{3.54} $$
利用 Slater 行列式与占据轨道的正交性,有序对密度可写成直接项减去同自旋交换项。 这里 $\rho_2$ 是有序电子对密度,$\gamma_\sigma$ 为自旋 $\sigma$ 的单粒子密度矩阵;在单 Slater 行列式近似下,空穴只有交换贡献。由方程(3.54)也可得
$$ \rho_{\mathrm{xc}}(\boldsymbol r,\boldsymbol r') =\frac{\rho_2(\boldsymbol r,\boldsymbol r')-\rho(\boldsymbol r)\rho(\boldsymbol r')}{\rho(\boldsymbol r)} =-\frac{\displaystyle\sum_\sigma\left|\sum_i\xi_i^\sigma(\boldsymbol r)\xi_i^{\sigma*}(\boldsymbol r')\right|^2}{\rho(\boldsymbol r)} \tag{3.55} $$
式(3.55)表明,$\rho_{\mathrm{xc}}(\boldsymbol{r},\boldsymbol{r}')$ 恒为负值,且有一个重要的性质:
$$ \int\rho_{\mathrm{xc}}(\boldsymbol r,\boldsymbol r')\,\mathrm d\boldsymbol r' =-\frac{\sum_{\sigma,i}|\xi_i^\sigma(\boldsymbol r)|^2}{\rho(\boldsymbol r)}=-1 \tag{3.56} $$
式(3.56)的物理意义非常明显:既然一个电子已经确定处于 $\boldsymbol{r}$ 终点位置,那么在所有其他 $\boldsymbol{r}'$ 终点位置所能找到的电子只有 $N-1$ 个,即一个电子不能同时存在于两处。在推导过程中,我们在 Hartree-Fock 近似下只考虑了自旋态相同的电子态,因此实际上只有交换效应而没有关联效应。所以式(3.56)中只有交换空穴密度 $\rho_{\mathrm{x}}$,而关联空穴密度 $\rho_{\mathrm{c}}\equiv0$。引入交换关联空穴 $\rho_{\mathrm{xc}}$,可以很方便地写出交换关联能的表达式为
$$ E_{\mathrm{xc}}^{\mathrm{HF}}=\frac{1}{2}\int\rho(\boldsymbol{r})\,\mathrm{d}\boldsymbol{r}\int\frac{\rho_{\mathrm{xc}}(\boldsymbol{r},\boldsymbol{r}')}{|\,\boldsymbol{r}-\boldsymbol{r}'\,|}\,\mathrm{d}\boldsymbol{r}' \tag{3.57} $$
如果将自旋变量 $\sigma$ 显式地表达出来,则方程(3.47)、方程(3.50)及方程(3.55)的形式略有变化:
$$ \begin{aligned} \rho(\boldsymbol{r},\sigma;\boldsymbol{r}',\sigma')={}&N(N-1)\times\sum_{\sigma_3,\sigma_4,\cdots,\sigma_N}\int|\,\Psi^{0}(\boldsymbol{r},\sigma;\boldsymbol{r}',\sigma';\boldsymbol{r}_3,\sigma_3;\cdots;\boldsymbol{r}_N,\sigma_N)\,|^{2}\,\mathrm{d}\boldsymbol{r}_3\,\mathrm{d}\boldsymbol{r}_4\cdots\mathrm{d}\boldsymbol{r}_N \end{aligned} \tag{3.58} $$
$$ \rho^{\sigma}(\boldsymbol{r})=\sum_{i}\xi_i^{\sigma*}(\boldsymbol{r})\xi_i^{\sigma}(\boldsymbol{r}) \tag{3.59} $$
$$ \rho_{\mathrm{xc}}(\boldsymbol r,\sigma;\boldsymbol r',\sigma') =-\delta_{\sigma\sigma'}\frac{\left|\sum_i\xi_i^{\sigma}(\boldsymbol r)\xi_i^{\sigma *}(\boldsymbol r')\right|^{2}}{\rho^\sigma(\boldsymbol r)} \tag{3.60} $$
上面不显含自旋情况的讨论在此仍然有效,这里不再详述。
根据上述讨论,可以引入双电子的对关联函数 $g(\boldsymbol{r},\sigma;\boldsymbol{r}',\sigma')$,且
$$ g(\boldsymbol{r},\sigma;\boldsymbol{r}',\sigma')=\frac{\rho(\boldsymbol{r},\sigma;\boldsymbol{r}',\sigma')}{\rho^{\sigma}(\boldsymbol{r})\rho^{\sigma'}(\boldsymbol{r}')}=1+\frac{\rho^{\sigma}(\boldsymbol{r})\rho_{\mathrm{xc}}(\boldsymbol{r},\sigma;\boldsymbol{r}',\sigma')}{\rho^{\sigma}(\boldsymbol{r})\rho^{\sigma'}(\boldsymbol{r}')} \tag{3.61} $$
在 Hartree-Fock 近似下,方程(3.55)或者方程(3.60)的分子显然就是单体密度矩阵 $n^{\sigma}(\boldsymbol{r},\boldsymbol{r}')$ 的二次方。因此式(3.61)可写为
$$ g(\boldsymbol{r},\sigma;\boldsymbol{r}',\sigma')=1-\delta_{\sigma\sigma'}\frac{|\,n^{\sigma}(\boldsymbol{r},\boldsymbol{r}')\,|^{2}}{\rho^{\sigma}(\boldsymbol{r})\rho^{\sigma'}(\boldsymbol{r}')} \tag{3.62} $$
从式(3.60)、式(3.62)可知,Hartree-Fock 近似只考虑了交换作用,而对另外的多体效应(如自旋相反波函数间的相互作用)未加考虑,这种影响通称为关联作用。
最后,给出显含自旋变量,但自旋非极化的情况下,对关联函数 $g_{\mathrm{x}}(\boldsymbol{r},\boldsymbol{r}')$ 的定义:
$$ g_{\mathrm{x}}(\boldsymbol{r},\boldsymbol{r}')=1-\frac{\displaystyle\sum_{\sigma}|\,n^{\sigma}(\boldsymbol{r},\boldsymbol{r}')\,|^{2}}{\rho(\boldsymbol{r})\rho(\boldsymbol{r}')} \tag{3.63} $$
在后文中我们将讨论特殊情况下 $g_{\mathrm{x}}(\boldsymbol{r};\boldsymbol{r}')$ 的一个解析解。
双电子对关联函数的引入,使得我们可以重新审视交换作用的物理意义。对于经典粒子,$g(\boldsymbol{r},\boldsymbol{r}')\equiv1$。因此,对关联函数对 1 的偏离反映了量子效应,而这种量子效应最终会体现在体系的能量表达式中。根据上面的讨论,可以知道对关联函数的大致行为。当 $\boldsymbol{r}'$ 趋于 $\boldsymbol{r}$ 时,$g(\boldsymbol{r},\boldsymbol{r}')$ 远小于 1(不考虑自旋时趋于 0,考虑自旋非极化时趋于 1/2);而当 $\boldsymbol{r}'$ 远离 $\boldsymbol{r}$ 时,$g(\boldsymbol{r},\boldsymbol{r}')$ 趋于 1。因此,若用经典电荷间的库仑相互作用描述电子与电子之间的相互作用,显然会严重高估排斥能,需要再引入一个吸引项以做校正 / 补偿,这就是交换能。所以,交换能并不反映某种新形式的粒子相互作用,而仅用于对被高估的库仑能进行修正。它是要求电子波函数反对称的必然结果。
闭壳层体系中的 Hartree-Fock 方程
正则 Hartree-Fock 方程(式(3.34))中的轨道是自旋轨道,分为空间部分和自旋部分。在实际求解体系中,由于自旋量子数是非经典的量子数,因此在数值求解中我们更关心如何得到自旋轨道中的空间部分。接下来介绍在闭壳层的情况下,如何将正则 Hartree-Fock 方程简化为空间轨道的一系列微分方程。闭壳层是指每一个占据的空间轨道都有自旋向上和自旋向下的电子配对填充的情况。
考虑一个含有偶数个电子的体系,用如下的方式对自旋轨道进行编号($i=1,2,\cdots,N/2$)
$$ \xi_{2i-1}(x)=\chi_i(r)\alpha(s)=\phi_{i\alpha} \tag{3.64} $$
$$ \xi_{2i}(x)=\chi_i(r)\beta(s)=\phi_{i\beta} \tag{3.65} $$
重新编号后,体系的基态波函数可以写为
$$ \phi^{0}_{\mathrm{RHF}}(\boldsymbol{x}_1,\boldsymbol{x}_2,\cdots,\boldsymbol{x}_N)=|\,\xi_i(1)\xi_j(2)\cdots\xi_k(N)\rangle_{\mathrm{S}}=|\,\phi_{1\alpha}(1)\phi_{1\beta}(2)\cdots\phi_{\frac{N}{2}\alpha}(N-1)\phi_{\frac{N}{2}\beta}(N)\rangle_{\mathrm{S}} \tag{3.66} $$
而双重求和可以表示为
$$ \sum_{i=1}^{N}\sum_{j=1}^{N}=\sum_{i\alpha}^{N/2}\sum_{j\alpha}^{N/2}+\sum_{i\alpha}^{N/2}\sum_{j\beta}^{N/2}+\sum_{i\beta}^{N/2}\sum_{j\alpha}^{N/2}+\sum_{i\beta}^{N/2}\sum_{j\beta}^{N/2} \tag{3.67} $$
另外,由于自旋态之间的正交性归一,可以得到
$$ \begin{gathered} \langle\phi_{i\alpha}\phi_{j\alpha}\,|\,\phi_{i\alpha}\phi_{j\alpha}\rangle=\langle\phi_{i\alpha}\phi_{j\beta}\,|\,\phi_{i\alpha}\phi_{j\beta}\rangle=\langle\phi_{i\beta}\phi_{j\alpha}\,|\,\phi_{i\beta}\phi_{j\alpha}\rangle=\langle\phi_{i\beta}\phi_{j\beta}\,|\,\phi_{i\beta}\phi_{j\beta}\rangle=\langle\chi_i\chi_j\,|\,\chi_i\chi_j\rangle\\ \langle\phi_{i\alpha}\phi_{j\alpha}\,|\,\phi_{j\alpha}\phi_{i\alpha}\rangle=\langle\phi_{i\beta}\phi_{j\beta}\,|\,\phi_{j\beta}\phi_{i\beta}\rangle=\langle\chi_i\chi_j\,|\,\chi_j\chi_i\rangle\\ \langle\phi_{i\alpha}\phi_{j\beta}\,|\,\phi_{j\beta}\phi_{i\alpha}\rangle=\langle\phi_{i\beta}\phi_{j\alpha}\,|\,\phi_{j\alpha}\phi_{i\beta}\rangle=0 \end{gathered} $$
因此 Hartree-Fock 基态能量表达式(3.22)可以化简为
$$ \begin{aligned} E_0&=\sum_{i=1}^{N}\langle\xi_i\,|\,h\,|\,\xi_i\rangle+\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\big(\langle\xi_i\xi_j\,|\,\xi_i\xi_j\rangle-\langle\xi_i\xi_j\,|\,\xi_j\xi_i\rangle\big)\\ &=2\sum_{i=1}^{N/2}\langle\chi_i\,|\,h\,|\,\chi_i\rangle+\sum_{i=1}^{N/2}\sum_{j=1}^{N/2}\big(2\langle\chi_i\chi_j\,|\,\chi_i\chi_j\rangle-\langle\chi_i\chi_j\,|\,\chi_j\chi_i\rangle\big)\\ &=2\sum_{i=1}^{N/2}h_{ii}+\sum_{i=1}^{N/2}\sum_{j=1}^{N/2}(2{J_{ij}}-{K_{ij}}) \end{aligned} \tag{3.68} $$
式中
$$ h_{ii}=\langle\chi_i\,|\,h\,|\,\chi_i\rangle=\int\mathrm{d}r\,\chi_i^{*}(r)h\chi_i(r) \tag{3.69} $$
$$ {J_{ij}}=\langle\chi_i\chi_j\,|\,\chi_i\chi_j\rangle=\int\mathrm{d}r_1\,\mathrm{d}r_2\,\chi_i^{*}(r_1)\chi_j^{*}(r_2)\frac{1}{r_{12}}\chi_i(r_1)\chi_j(r_2) \tag{3.70} $$
$$ {K_{ij}}=\langle\chi_i\chi_j\,|\,\chi_j\chi_i\rangle=\int\mathrm{d}r_1\,\mathrm{d}r_2\,\chi_i^{*}(r_1)\chi_j^{*}(r_2)\frac{1}{r_{12}}\chi_j(r_1)\chi_i(r_2) \tag{3.71} $$
这里的 $J_{ij}$ 和 $K_{ij}$ 是积分所得的标量,分别称为直接积分和交换积分;它们不是前文作用于轨道的算符 $\hat J_j$ 与 $\hat K_j$。
每一对自旋方向相反的电子贡献 ${J_{ij}}$ 的库仑斥能,而每一对自旋方向相同的电子贡献 ${J_{ij}}-{K_{ij}}$ 的库仑能和交换能。
同样,在求解单电子方程的时候,也需要将对自旋轨道的方程转成对空间轨道的微分方程。考察自旋 $\alpha$ 的自旋轨道
$$ \begin{aligned} &{\hat F}(x_1)\xi_i(x_1)=\varepsilon_i\xi_i(x_1)\\ \Rightarrow\quad&{\hat F}(x_1)\chi_i(r_1)\alpha(s)=\varepsilon_i\chi_i(r_1)\alpha(s)\\ \Rightarrow\quad&\int\mathrm{d}s\,\alpha^{*}(s){\hat F}(x_1)\chi_i(r_1)\alpha(s)=\int\mathrm{d}s\,\alpha^{*}(s)\varepsilon_i\chi_i(r_1)\alpha(s)\\ \Rightarrow\quad&\left\{\int\mathrm{d}s\,\alpha^{*}(s){\hat F}(x_1)\alpha(s)\right\}\chi_i(r_1)=\varepsilon_i\chi_i(r_1)\\ \Rightarrow\quad&{\hat F}(r_1)\chi_i(r_1)=\varepsilon_i\chi_i(r_1) \end{aligned} \tag{3.72} $$
通过上面的推导可知,Fock 算符在闭壳层情况下可以写为
$$ \begin{aligned} {\hat F}(r_1)&=\int\mathrm{d}s\,\alpha^{*}(s){\hat F}(x_1)\alpha(s)=\int\mathrm{d}s\,\alpha^{*}(s)\left[h(x_1)+\sum_{i=1}^{N}\big({\hat J}_i(x_1)-{\hat K}_i(x_1)\big)\right]\alpha(s)\\ &=h(r_1)+\sum_{i}^{N/2}\big[2{\hat J}_i(r_1)-{\hat K}_i(r_1)\big] \end{aligned} \tag{3.73} $$
式中
$$ {\hat J}_i(r_1)=\int\mathrm{d}r_2\,\chi_i^{*}(r_2)\frac{1}{r_{12}}\chi_i(r_2) \tag{3.74} $$
$$ ({\hat K}_i f)(\boldsymbol r_1)=\chi_i(\boldsymbol r_1)\int\!\mathrm d\boldsymbol r_2\,\frac{\chi_i^*(\boldsymbol r_2)f(\boldsymbol r_2)}{|\boldsymbol r_1-\boldsymbol r_2|} \tag{3.75} $$
式(3.74)给出空间轨道 $\chi_i$ 产生的库仑势,$\hat J_i$ 通过该势乘在待作用的空间函数上;式(3.75)明确给出 $\hat K_i$ 对任意空间函数 $f$ 的作用,交换项要对 $f$ 作积分,不能视为普通的局域势函数。
因此在闭壳层下,Hartree-Fock 方程为
$$ \left[h(\boldsymbol r_1)+\sum_{j=1}^{N/2}\bigl(2{\hat J}_j-{\hat K}_j\bigr)\right]\chi_i(\boldsymbol r_1)=\varepsilon_i\chi_i(\boldsymbol r_1) \tag{3.76} $$
开壳层体系中的 Hartree-Fock 方程
当体系中含有奇数个电子时,以及对于远离平衡态的解离过程等,由开壳层方法通常能够得到比闭壳层方法更加准确的结果。其处理方法和闭壳层的 Hartree-Fock 方法类似,但是我们不强制要求不同自旋方向的电子配对占据相同的空间轨道,而是允许同一个能级的不同自旋方向的电子占据不同的空间轨道。有:
$$ \Phi_{\mathrm{UHF}}=\bigl|\,\phi_{1\alpha}\alpha,\ldots,\phi_{N_\alpha\alpha}\alpha,\bar\phi_{1\beta}\beta,\ldots,\bar\phi_{N_\beta\beta}\beta\,\bigr\rangle_{\mathrm S},\qquad N=N_\alpha+N_\beta \tag{3.77} $$
在这种情况下,自旋方向相同的空间轨道互相正交,而自旋方向不同的空间轨道的交叠由矩阵 $\boldsymbol{S}$ 描述,其中 $S_{i\alpha,j\beta}=\langle\phi_{i\alpha}\,|\,\bar{\phi}_{j\beta}\rangle$。和闭壳层 Hartree-Fock 方程推导类似,可以通过对自旋积分将方程转化成只涉及空间轨道的微分方程。在开壳层 Hartree-Fock 方法中,将得到关于自旋 $\alpha$ 和 $\beta$ 的电子分立的两组关于空间轨道的方程:
$$ {\hat F}^\alpha\phi_{i\alpha}=\varepsilon_i^\alpha\phi_{i\alpha},\qquad i=1,\ldots,N_\alpha \tag{3.78} $$
$$ {\hat F}^\beta\bar\phi_{j\beta}=\varepsilon_j^\beta\bar\phi_{j\beta},\qquad j=1,\ldots,N_\beta \tag{3.79} $$
式中
$$ {\hat F}^{\alpha}(r_1)=h(r_1)+\sum_{i=1}^{N_\alpha}\big[{\hat J}_i^{\alpha}(r_1)-{\hat K}_i^{\alpha}(r_1)\big]+\sum_{i=1}^{N_\beta}{\hat J}_i^{\beta}(r_1) \tag{3.80} $$
$$ {\hat F}^{\beta}(r_1)=h(r_1)+\sum_{i=1}^{N_\alpha}{\hat J}_i^{\alpha}(r_1)+\sum_{i=1}^{N_\beta}\left[{\hat J}_i^{\beta}(r_1)-{\hat K}_i^{\beta}(r_1)\right] \tag{3.81} $$
上标 $\alpha$、$\beta$ 区分两个自旋通道的 Fock 算符。$\hat J_i^{\sigma}$ 是由自旋通道 $\sigma$ 中占据的第 $i$ 个空间轨道产生的库仑算符,对两种自旋的电子均有贡献;$\hat K_i^{\sigma}$ 是对应的交换算符,只进入同自旋通道的 Fock 算符。
Hartree-Fock 方程的矩阵表达
在实际求解 Hartree-Fock 自洽场方程的过程中,通常用已知的 $K$ 个基组对第 $i$ 个分子轨道的空间部分 $\chi_i(\boldsymbol{r})$ 进行展开:
$$ \chi_i(\boldsymbol{r})=\sum_{\mu=1}^{K}C_{\mu i}\zeta_\mu(\boldsymbol{r}) \tag{3.82} $$
式中:$\zeta_\mu(\boldsymbol{r})$ 中的下标 $\mu$ 用于同时标记不同原子中心和位于该中心的基组。
以闭壳层的 Hartree-Fock 方程为例,在上面的基组展开下,Hartree-Fock 方程可以写为
$$ \sum_{\nu}F_{\mu\nu}C_{\nu i}=\varepsilon_i\sum_{\nu}S_{\mu\nu}C_{\nu i} \tag{3.83} $$
式中
$$ S_{\mu\nu}=\int\mathrm{d}\boldsymbol{r}\,\zeta_\mu^{*}(\boldsymbol{r})\zeta_\nu(\boldsymbol{r}),\quad F_{\mu\nu}=\int\mathrm{d}\boldsymbol{r}\,\zeta_\mu^{*}(\boldsymbol{r}){\hat F}(\boldsymbol{r})\zeta_\nu(\boldsymbol{r}) $$
如果写成矩阵的形式,式(3.83)可以表示为
$$ \boldsymbol{FC}=\boldsymbol{SCE} \tag{3.84} $$
此方程称为 Roothaan 方程,其中
$$ \boldsymbol{C}=\begin{bmatrix}C_{11}&C_{12}&\cdots&C_{1K}\\C_{21}&C_{22}&\cdots&C_{2K}\\\vdots&\vdots&&\vdots\\C_{K1}&C_{K2}&\cdots&C_{KK}\end{bmatrix} \tag{3.85} $$
$$ \boldsymbol{E}=\begin{bmatrix}\varepsilon_1&&&\\&\varepsilon_2&&\\&&\ddots&\\&&&\varepsilon_K\end{bmatrix} \tag{3.86} $$
方程(3.84)中,$\boldsymbol{S}$、$\boldsymbol{C}$、$\boldsymbol{E}$ 三个矩阵都比较简单。$\boldsymbol{F}$ 矩阵的计算由于涉及双电子积分而比较复杂,因此给出 $\boldsymbol{F}$ 矩阵元比较详细的推导过程:
$$ \begin{aligned} F_{\mu\nu}&=\int\mathrm{d}\boldsymbol{r}\,\zeta_\mu^{*}(\boldsymbol{r}){\hat F}(\boldsymbol{r})\zeta_\nu(\boldsymbol{r})=\int\mathrm{d}\boldsymbol{r}\,\zeta_\mu^{*}(\boldsymbol{r})\left[h(\boldsymbol{r})+\sum_{i=1}^{N/2}\big(2{\hat J}_i(\boldsymbol{r})-{\hat K}_i(\boldsymbol{r})\big)\right]\zeta_\nu(\boldsymbol{r})\\ &=\int\mathrm{d}\boldsymbol{r}\,\zeta_\mu^{*}(\boldsymbol{r})h(\boldsymbol{r})\zeta_\nu(\boldsymbol{r})+\sum_{i=1}^{N/2}\int\mathrm{d}\boldsymbol{r}\,\zeta_\mu^{*}(\boldsymbol{r})\big(2{\hat J}_i(\boldsymbol{r})-{\hat K}_i(\boldsymbol{r})\big)\zeta_\nu(\boldsymbol{r})=H_{\mu\nu}+G_{\mu\nu} \end{aligned} \tag{3.87} $$
式中:$H_{\mu\nu}$ 为单电子积分;$G_{\mu\nu}$ 可以简化为单电子密度矩阵和双电子积分的乘积,即
$$ \begin{aligned} G_{\mu\nu}&=\sum_{\lambda,\sigma=1}^{K}P_{\lambda\sigma} \left[\langle\mu\sigma\,|\,\nu\lambda\rangle -\tfrac12\langle\mu\sigma\,|\,\lambda\nu\rangle\right],\\ \langle\mu\sigma\,|\,\nu\lambda\rangle &=\iint\!\mathrm d\boldsymbol r_1\,\mathrm d\boldsymbol r_2\, \frac{\zeta_\mu^*(\boldsymbol r_1)\zeta_\sigma^*(\boldsymbol r_2) \zeta_\nu(\boldsymbol r_1)\zeta_\lambda(\boldsymbol r_2)} {|\boldsymbol r_1-\boldsymbol r_2|} \end{aligned} \tag{3.88} $$
其中
$$ P_{\lambda\sigma}=2\sum_{i=1}^{N/2}C_{\lambda i}C_{\sigma i}^{*} \tag{3.89} $$
开壳层的 Hartree-Fock 方程在基组展开下可以类似地化为方程组进行求解,所对应的方程组称为 Pople-Nesbet 方程,在此不赘述,有兴趣的读者可以参考 Szabo 和 Ostlund 的 Modern Quantum Chemistry:Introduction to Advanced Electronic Structure Theory 一书。
Koopmans 定理
Koopmans 在 1933 年证明了下述定理\cite{koopmans1934zuordnung}:
Koopmans 定理 在 Hartree-Fock 近似下,一个占据(非占据)轨道 $\xi_k$ 的本征值 $\varepsilon_k$ 等于将一个电子从(向)该轨道移走(填充)且其他各轨道保持不变的情况下,Hartree-Fock 总能(见式(3.22))的变化 $E_0(N)-E_0(N-1)$。
该定理的证明比较简单,以从 $\xi_k$ 移走一个电子为例。将该电子从 $\xi_k$ 移走前后的所有占据态代入方程(3.27),可得
$$ E_0(N)-E_0(N-1)\,|_{\xi_k}=\langle\xi_k\,|\,h\,|\,\xi_k\rangle+\sum_{i=1}^{N}\big(\langle\xi_k\xi_i\,|\,\xi_k\xi_i\rangle-\langle\xi_k\xi_i\,|\,\xi_i\xi_k\rangle\big) \tag{3.90} $$
而将 $\langle\xi_k\,|$ 作用到 Hartree-Fock 自洽场方程(式(3.31)),可得
$$ \begin{aligned}\varepsilon_k&=\left\langle\xi_k\left|h+\sum_{j\neq k}^{N}({\hat J}_j-{\hat K}_j)\right|\xi_k\right\rangle\\&=\langle\xi_k|h|\xi_k\rangle+\sum_{j=1}^{N}\bigl(\langle\xi_k\xi_j|\xi_k\xi_j\rangle-\langle\xi_k\xi_j|\xi_j\xi_k\rangle\bigr)\end{aligned} \tag{3.91} $$
上面两个方程明显相等。至此,定理得证。
Koopmans 定理的重要意义在于它明确地给出了 Hartree-Fock 单电子自洽场方程本征能级的物理意义。需要强调的是,Hartree-Fock 方法本身是一种对多体体系并不十分准确的单电子近似。因此尽管 Koopmans 定理成立,但是 $\varepsilon_k$ 不能理解为分子轨道的真实本征能级。最为明显的一个例子就是 Hartree-Fock 近似没有包含关联效应,对自能修正的描述也不完全,因此其严重高估了最高占据态和最低非占据态的能量差,即能隙 $E_{\mathrm{g}}$。