Study Note：Hartree–Fock

之前在进行 QM 计算时，软件在做什么？里，简单介绍过 HF 和 SCF 的计算过程。这次把读 Szabo 与 Ostlund 的 Modern Quantum Chemistry 时留下的推导重新整理一遍。先从一个分子的能量由什么组成讲起，再看 HF 怎样处理电子之间的相互作用，最后把这些想法写成能在程序中求解的轨道方程。

本文先讨论一般的 HF，再转到闭壳层限制性 Hartree–Fock（restricted Hartree–Fock，RHF）。后半部分采用实基函数和实轨道系数，所有量均使用原子单位。需要的数学主要是积分、矩阵乘法、求导和拉格朗日乘子法。

1 分子能量与 HF 近似

1.1 BO 近似

一个分子包含原子核和电子。这些粒子运动时有动能，带电粒子之间又有静电作用：电子与核相互吸引，电子之间、原子核之间相互排斥。如果要同时描述所有粒子的运动和相互作用，问题会很复杂。

原子核远比电子重，通常运动得慢得多。因此，可以先固定一套原子核的位置，求电子在这套核构型下的状态；再由不同核构型下的能量，研究原子核如何运动。这就是 Born–Oppenheimer 近似（BO 近似）的基本思路。电子状态随所选核构型而改变，核构型固定以后，我们就可以先集中处理电子问题。

对于固定核构型下的一次单点能计算，可以先把能量分成两部分：一部分与电子有关，包括电子动能、电子与核的吸引能、电子之间的排斥能；另一部分是原子核之间的排斥能。后者只由核的电荷和核间距离决定，在这次电子计算中是一个常数。因此：

\[E_{\mathrm{tot}}=E_{\mathrm{elec}}+E_{\mathrm{nuc}}\]

式中：

• $E_{\mathrm{tot}}$：固定核构型下的总能量；下标 $\mathrm{tot}$ 表示总量。
• $E_{\mathrm{elec}}$：电子动能、电子与核的吸引能、电子之间的排斥能之和；下标 $\mathrm{elec}$ 表示电子部分。
• $E_{\mathrm{nuc}}$：所有原子核之间的排斥能；下标 $\mathrm{nuc}$ 表示原子核部分。

核间排斥能可以直接按每一对核之间的库仑排斥计算：

\[E_{\mathrm{nuc}}=\sum_{A<B}\frac{Z_AZ_B}{R_{AB}}\]

式中：

• $A,B$：两个原子核的编号。
• $Z_A,Z_B$：这两个原子核的电荷数，例如氢核为 $1$、碳核为 $6$。
• $R_{AB}$：第 $A$、$B$ 个原子核之间的距离。
• $\sum$：把指定范围内的各项相加；这里的 $A<B$ 使每一对不同原子核只算一次。

这样，只要求出电子能量，最后加上核间排斥能，就得到了这套核构型下的总能量。下面先看电子部分。

1.2 平均场近似

先单独看一个电子：它自身的运动贡献动能，它受到所有原子核的吸引又贡献势能。在原子核位置固定的前提下，这两部分的算符都只涉及一个电子的坐标，因此可以合称为单电子部分。

电子之间的排斥则要同时知道两个电子的位置：它们相距越近，排斥越强。这一部分要按电子对来计算，称为双电子部分。

单电子部分可以逐个电子处理，双电子排斥却把它们联系在了一起：一个电子的状态会受到其他电子的影响，而其他电子的状态也会随之改变。随着电子数增加，同时描述所有电子如何相互配合会变得很困难。

Hartree–Fock（HF）采用平均场近似（mean-field approximation）来简化这个问题。我们用一组单电子轨道描述电子状态，求某个轨道时，由其他占据轨道的电子分布共同给出它所受到的有效作用。这个有效作用包含平均的库仑排斥，以及由电子波函数反对称性产生的交换作用。这样，求解时就能转成一组相互依赖的单电子轨道方程。

这些轨道又决定了平均场本身，所以需要反复更新：先给出一组轨道，按它们的电子分布计算有效作用，再用这个作用求出新轨道，直到前后相互一致。这就是后文要写成程序的自洽场过程。

要把这个思路写成公式，先列出计算电子能量所用的哈密顿算符。将所有电子的单电子部分，以及所有电子对的排斥部分加起来：

\[\hat H_{\mathrm{elec}}=\sum_{k=1}^{N}\hat h(k)+\sum_{k<l}^{N}\frac{1}{r_{kl}}\]

式中：

• $\hat H_{\mathrm{elec}}$：电子哈密顿算符，用来计算电子体系的能量。
• $N$：电子总数。
• $k,l$：电子的编号，取值为 $1,\ldots,N$。
• $\hat h(k)$：只作用于第 $k$ 个电子的能量算符，包含动能和电子与核的吸引能。
• $r_{kl}$：第 $k$、$l$ 个电子之间的距离，$1/r_{kl}$ 给出这一对电子的排斥作用。
• $k<l$：让每一对不同电子只算一次。
• 字母上方的帽号（如 $\hat h$）：表示算符，也就是作用于波函数的一项数学操作。

将单电子部分也展开，就能看到前面说的动能和核吸引能：

\[\hat h(k)=-\frac{1}{2}\nabla_k^2-\sum_A\frac{Z_A}{r_{kA}}\]

式中：

• $\nabla_k^2$：对第 $k$ 个电子的三个空间坐标分别求二阶导数后相加；它描述波函数在空间中的弯曲程度，乘上 $-1/2$ 后给出动能算符。
• $r_{kA}$：第 $k$ 个电子到第 $A$ 个原子核的距离。
• $\sum_A$：对所有原子核求和，将这个电子受到的各个核的吸引作用加在一起。

上式第一项是动能算符，第二项是电子与原子核之间的吸引势能算符。这就写出了固定核构型下的非相对论电子哈密顿量。接下来，需要选取合适的波函数形式来实现 HF 的平均场描述。

1.3 Slater 行列式

要把上面的平均场想法落实成可以求解的数学形式，还需要决定怎样用单电子轨道构成整个体系的波函数。电子是费米子，交换任意两个电子的全部坐标，波函数必须变号。Slater 行列式恰好能满足这个要求。

用 $x=(\mathbf r,\omega)$ 表示空间与自旋坐标，用 $\eta_a(x)$ 表示自旋轨道，可以构造归一化的 Slater 行列式：

\[\Psi(x_1,\ldots,x_N)=\frac{1}{\sqrt{N!}}\det[\eta_a(x_k)]_{k,a=1}^{N}\]

式中：

• $\Psi$：整个 $N$ 电子体系的波函数，同时依赖所有电子的坐标。
• $x=(\mathbf r,\omega)$：一个电子的完整坐标，其中 $\mathbf r$ 是空间位置，$\omega$ 是用来区分电子自旋状态的变量。
• $x_k$：第 $k$ 个电子的完整坐标。
• $\eta_a(x)$：第 $a$ 个自旋轨道，即同时描述一个电子空间分布与自旋的单电子函数；$a$ 对这 $N$ 个占据自旋轨道编号。
• $\det[\eta_a(x_k)]$：由这些轨道值排成的矩阵的行列式；第 $k$ 行对应电子坐标，第 $a$ 列对应自旋轨道。
• $N!$：$1\times2\times\cdots\times N$；前面的 $1/\sqrt{N!}$ 保证在轨道正交归一时，整个波函数也归一化。

这里要求占据自旋轨道满足正交归一条件：

\[\int\eta_a^*(x)\eta_b(x)\,dx=\delta_{ab}\]

式中：

• $b$：另一个占据自旋轨道的编号，与 $a$ 的取值范围相同。
• 上标 $*$：取复共轭，即把复数的虚部变号；对实函数，它不改变函数值。
• $\int\cdots dx$：对电子的全部空间位置积分，并对自旋求和。
• $\delta_{ab}$：Kronecker delta，两个编号相同为 $1$，不同为 $0$。这里表示每个轨道自身的重叠为 $1$，不同轨道的重叠为 $0$。

交换两个电子的坐标相当于交换行列式的两行，行列式随之变号。如果两个占据自旋轨道相同，行列式的两列相同，整个波函数为零。

HF 的变分目标，就是在单个 Slater 行列式所能表示的波函数中寻找最低能量。 轨道本身是待优化的函数；采用有限基组以后，待优化的量就变成轨道展开系数。

现在，计算能量所需的哈密顿算符，以及 HF 采用的波函数形式都有了。接下来先算这个行列式的能量，看看电子之间的作用怎样变成库仑项和交换项，再用能量最小的要求确定轨道。

2 行列式的能量

单电子部分给出各轨道的动能和核吸引能；电子排斥则要对电子对求和。对 Slater 行列式计算排斥能时，反对称组合中的各个乘积会彼此交叉，因此结果中会同时出现库仑项和交换项。先从只有两个电子的情况展开，能直接看清这两项的来由。

两个正交归一的自旋轨道 $\eta_a$、$\eta_b$ 构成：

\[\Psi(1,2)=\frac{1}{\sqrt{2}}[\eta_a(1)\eta_b(2)-\eta_b(1)\eta_a(2)]\]

式中：

• 函数括号中的 $1,2$：分别是完整电子坐标 $x_1,x_2$ 的简写，例如 $\eta_a(1)$ 就是 $\eta_a(x_1)$。

为方便展开，将两个轨道乘积分别记为 $A$、$B$，将排斥算符记为 $g$。电子排斥能为：

\[\langle\Psi\vert g\vert\Psi\rangle=\frac{1}{2}[\langle A\vert g\vert A\rangle+\langle B\vert g\vert B\rangle-\langle A\vert g\vert B\rangle-\langle B\vert g\vert A\rangle]\]

式中：

• $A=\eta_a(1)\eta_b(2)$、$B=\eta_b(1)\eta_a(2)$：交换两个电子的轨道安排所得到的两个乘积函数。
• $g=1/r_{12}$：两个电子之间的库仑排斥算符，作用方式是乘以两电子距离的倒数。
• $r_{12}$：电子 $1$ 与电子 $2$ 的距离。
• $\langle A\vert g\vert B\rangle$：先让 $g$ 作用于 $B$，再乘以 $A$ 的复共轭，并对两个电子的全部坐标积分。
• $\langle\Psi\vert g\vert\Psi\rangle$：排斥算符在波函数 $\Psi$ 下的期望值，即这个态的平均电子排斥能。

前两项通过交换积分变量 $x_1,x_2$ 可知相等，定义为库仑积分：

\[J^{\mathrm{spin}}_{ab}=\iint\eta_a^*(x_1)\eta_a(x_1)\frac{1}{r_{12}}\eta_b^*(x_2)\eta_b(x_2)\,dx_1dx_2\]

式中：

• $J^{\mathrm{spin}}_{ab}$：自旋轨道 $a,b$ 的库仑积分，可以看作两个轨道对应的电子分布之间的静电排斥能。
• 上标 $\mathrm{spin}$：表示积分仍包含自旋，后面会将它转成空间轨道的积分。
• $\iint\cdots dx_1dx_2$：分别对两个电子的完整坐标积分。

后两项也相等，定义为交换积分：

\[K^{\mathrm{spin}}_{ab}=\iint\eta_a^*(x_1)\eta_b(x_1)\frac{1}{r_{12}}\eta_b^*(x_2)\eta_a(x_2)\,dx_1dx_2\]

式中：

• $K^{\mathrm{spin}}_{ab}$：自旋轨道 $a,b$ 的交换积分，是两个轨道乘积在反对称波函数中产生的交叉项所对应的积分。

所以：

\[\left\langle\Psi\middle\vert\frac{1}{r_{12}}\middle\vert\Psi\right\rangle=J^{\mathrm{spin}}_{ab}-K^{\mathrm{spin}}_{ab}\]

交换项的负号就来自波函数中的反对称组合。它是计算同一个电子排斥算符的期望值时自然出现的交叉项。

对于 $N$ 电子行列式，正交归一性使单电子部分化为各占据自旋轨道贡献的和；双电子部分则对所有占据轨道对求和：

\[E_{\mathrm{elec}}=\sum_{a=1}^{N}\langle\eta_a\vert\hat h\vert\eta_a\rangle+\frac{1}{2}\sum_{a,b=1}^{N}(J^{\mathrm{spin}}_{ab}-K^{\mathrm{spin}}_{ab})\]

式中：

• $\langle\eta_a\vert\hat h\vert\eta_a\rangle$：一个电子占据自旋轨道 $a$ 时的动能与核吸引能之和。
• $\sum_{a,b=1}^{N}$：让 $a,b$ 各自取遍全部 $N$ 个占据自旋轨道，再把所有轨道对的贡献相加。

这里的 $1/2$ 是因为双重求和把轨道对 $(a,b)$ 和 $(b,a)$ 都计入了。求和可以包含 $a=b$，因为这时 $J^{\mathrm{spin}}{aa}=K^{\mathrm{spin}}{aa}$，两项相消。

3 闭壳层 RHF 能量

前面的自旋轨道同时包含电子的空间与自旋信息。对于所有电子都成对占据空间轨道的闭壳层体系，可以先把自旋部分算完，将能量写成只含空间轨道的形式。这样也能看清楚双占据带来的系数是怎样出现的。

对于闭壳层 RHF，$N$ 为偶数，每个占据空间轨道 $\phi_i(\mathbf r)$ 分别与 $\alpha$、$\beta$ 自旋函数相乘，形成两个自旋轨道：

\[\eta_{i\alpha}(x)=\phi_i(\mathbf r)\alpha(\omega),\qquad\eta_{i\beta}(x)=\phi_i(\mathbf r)\beta(\omega)\]

式中：

• $\phi_i(\mathbf r)$：第 $i$ 个占据空间轨道，只描述电子的空间分布。
• $i,j$：占据空间轨道的编号，取值为 $1,\ldots,n_{\mathrm{occ}}$。
• $n_{\mathrm{occ}}=N/2$：占据空间轨道数；下标 $\mathrm{occ}$ 是 occupied 的缩写，表示被电子占据。
• $\alpha(\omega),\beta(\omega)$：两个相互正交的自旋函数，分别对应选定方向上自旋投影的两个取值，通常称为自旋向上和向下。
• $\eta_{i\alpha},\eta_{i\beta}$：共用同一个空间轨道 $\phi_i$，但自旋不同的两个自旋轨道。

自旋函数满足：

\[\langle\alpha\vert\alpha\rangle=\langle\beta\vert\beta\rangle=1,\qquad\langle\alpha\vert\beta\rangle=0\]

式中：

• $\langle\alpha\vert\beta\rangle$：两个自旋函数的重叠，即将前者取复共轭后与后者相乘，再对自旋求和；同一种自旋的重叠为 $1$，两种不同自旋的重叠为 $0$。

单电子算符不作用于自旋，每个空间轨道贡献两次相同的单电子能量，因此：

\[E_{\mathrm{1e}}=2\sum_{i=1}^{n_{\mathrm{occ}}}h_{ii},\qquad h_{ii}=\langle\phi_i\vert\hat h\vert\phi_i\rangle\]

式中：

• $E_{\mathrm{1e}}$：总电子能量中的单电子部分，即所有电子的动能与核吸引能之和。
• $h_{ii}$：一个电子占据空间轨道 $i$ 时的单电子能量。
• 前面的 $2$：每个占据空间轨道中有两个电子。

用 $J_{ij}$、$K_{ij}$ 表示对自旋积分以后得到的空间轨道库仑、交换积分。对应的积分形式在下面给出。

双电子部分可以按自旋逐项数出来。给定一对空间轨道 $i,j$：

• 库仑项有 $\alpha\alpha$、$\alpha\beta$、$\beta\alpha$、$\beta\beta$ 四种组合，自旋积分都为 $1$，共得到 $4J_{ij}$。
• 交换项含有两个自旋函数之间的重叠。只有 $\alpha\alpha$、$\beta\beta$ 两种组合保留下来，共得到 $2K_{ij}$。

再乘原来双重求和前面的 $1/2$，就得到 $2J_{ij}-K_{ij}$。于是：

\[\boxed{E_{\mathrm{elec}}=2\sum_i^{\mathrm{occ}}h_{ii}+\sum_{ij}^{\mathrm{occ}}(2J_{ij}-K_{ij})}\]

式中：

• $J_{ij}$：占据空间轨道 $i,j$ 的库仑积分，描述两份轨道电子分布之间的静电排斥。
• $K_{ij}$：占据空间轨道 $i,j$ 的交换积分，描述同自旋电子因波函数反对称性而产生的交换贡献。
• $\sum_i^{\mathrm{occ}}$、$\sum_{ij}^{\mathrm{occ}}$：分别对全部占据空间轨道、全部占据空间轨道对求和。

后文所有未特别标注的轨道与积分均为空间形式。

采用化学家记号：

\[(pq\vert rs)=\iint\phi_p(\mathbf r_1)\phi_q(\mathbf r_1)\frac{1}{r_{12}}\phi_r(\mathbf r_2)\phi_s(\mathbf r_2)\,d\mathbf r_1d\mathbf r_2\]

式中：

• $p,q,r,s$：任意空间轨道的编号，可以是占据轨道，也可以是未被占据的虚轨道。
• $\phi_p,\phi_q,\phi_r,\phi_s$：这四个编号对应的空间轨道。
• $\mathbf r_1,\mathbf r_2$：两个电子的空间位置，$d\mathbf r_1d\mathbf r_2$ 表示对这两个三维位置分别积分。
• $(pq\vert rs)$：双电子积分的化学家记号；竖线左侧的两个轨道取电子 $1$ 的坐标，右侧的两个轨道取电子 $2$ 的坐标。

对于这里采用的实轨道：

\[J_{ij}=(ii\vert jj),\qquad K_{ij}=(ij\vert ji)\]

例如，一个空间轨道被两个自旋相反的电子占据时，$J_{11}=K_{11}$，因而电子排斥能是 $2J_{11}-K_{11}=J_{11}$，正好是这一对电子之间的一次库仑排斥。

4 AO 展开与密度矩阵

现在已经知道，给定一组轨道，就能用积分计算它们构成的行列式能量。要继续优化轨道，程序还需要一种具体的轨道表示方法。一个直接的办法是选取一组已知函数，用它们的线性组合表示未知轨道；这样，待求的函数就变成了一组待求的系数。

4.1 轨道展开

取 $M$ 个已知的基函数 ${\chi_\mu}$，将分子轨道展开为：

\[\phi_p(\mathbf r)=\sum_{\mu=1}^{M}C_{\mu p}\chi_\mu(\mathbf r)\]

式中：

• $M$：选取的基函数总数。
• $\chi_\mu(\mathbf r)$：第 $\mu$ 个已知的基函数，用来组合出待求的分子轨道。
• $\mu,\nu,\lambda,\sigma$：基函数编号，后文的 $\kappa$ 也表示基函数编号，取值均为 $1,\ldots,M$。
• $C_{\mu p}$：第 $\mu$ 个基函数在第 $p$ 个分子轨道中的展开系数，也就是线性组合时乘在该基函数前的权重。
• $C$：由全部展开系数组成的矩阵，每一列对应一个分子轨道。

在通常的分子计算中，这些函数以原子为中心，因此习惯上称为原子轨道（atomic orbital，AO）基函数；由它们组合出的轨道称为分子轨道（molecular orbital，MO）。

MO 满足正交归一条件：

\[\langle\phi_p\vert\phi_q\rangle=\delta_{pq}\]

将展开式代入：

\[\langle\phi_p\vert\phi_q\rangle=\sum_{\mu\nu}C_{\mu p}C_{\nu q}\langle\chi_\mu\vert\chi_\nu\rangle=\sum_{\mu\nu}C_{\mu p}S_{\mu\nu}C_{\nu q}\]

其中重叠矩阵为：

\[S_{\mu\nu}=\langle\chi_\mu\vert\chi_\nu\rangle\]

式中：

• $S_{\mu\nu}$：第 $\mu$、$\nu$ 个基函数的重叠积分，即将两个基函数相乘后对空间积分，衡量它们在空间中的重叠。
• $S$：将所有这些重叠积分排成的 $M\times M$ 矩阵，称为重叠矩阵。
• $\delta_{pq}$：前面同一类 Kronecker delta，用于空间轨道编号；$p=q$ 时为 $1$，否则为 $0$。

右侧正是矩阵 $C^TSC$ 的第 $(p,q)$ 个元素，所以：

\[\boxed{C^TSC=I}\]

式中：

• 上标 $T$：矩阵转置，即把行与列互换。
• $I$：单位矩阵，对角元素为 $1$，其余元素为 $0$；此处表示各 MO 之间的重叠组成单位矩阵。
• 上标 $\dagger$：共轭转置，即先取复共轭再转置；如果使用复系数，正交归一条件相应写为 $C^\dagger SC=I$。

AO 基函数之间通常存在重叠，$S$ 因而进入归一化条件。本文采用实系数，所以后文使用转置 $T$。

4.2 密度矩阵

有了轨道，还可以问一个更直观的问题：电子在空间中怎样分布？每个归一化轨道的模平方给出一个电子的概率密度，闭壳层的每个占据空间轨道中又有两个电子。因此，把各轨道的贡献加起来，电子数密度为：

\[\rho(\mathbf r)=2\sum_i^{\mathrm{occ}}\phi_i(\mathbf r)^2\]

式中：

• $\rho(\mathbf r)$：位置 $\mathbf r$ 处的电子数密度，表示该位置附近单位体积内的平均电子数。
• $\phi_i(\mathbf r)^2$：实空间轨道 $i$ 对应的单电子概率密度；前面的 $2$ 计入这个轨道中的两个电子。

将每个轨道的 AO 展开代入：

\[\rho(\mathbf r)=2\sum_i^{\mathrm{occ}}\sum_{\mu\nu}C_{\mu i}C_{\nu i}\chi_\mu(\mathbf r)\chi_\nu(\mathbf r)\]

交换求和顺序，把只依赖系数的部分放在一起：

\[\rho(\mathbf r)=\sum_{\mu\nu}\left(2\sum_i^{\mathrm{occ}}C_{\mu i}C_{\nu i}\right)\chi_\mu(\mathbf r)\chi_\nu(\mathbf r)\]

括号中的系数在后面的能量计算里还会反复出现。把这些成对的轨道系数组合收集起来，定义包含双占据因子 $2$ 的密度矩阵：

\[\boxed{P_{\mu\nu}=2\sum_i^{\mathrm{occ}}C_{\mu i}C_{\nu i}},\qquad P=2C_{\mathrm{occ}}C_{\mathrm{occ}}^T\]

式中：

• $P_{\mu\nu}$：密度矩阵的第 $(\mu,\nu)$ 个元素，是电子密度展开后乘在基函数对 $\chi_\mu\chi_\nu$ 前面的系数。
• $P$：包含所有 $P_{\mu\nu}$ 的密度矩阵，用一组数记录占据轨道产生的电子分布。
• $C_{\mathrm{occ}}$：只保留占据轨道各列的系数矩阵，维数为 $M\times n_{\mathrm{occ}}$。

于是 $P$ 是一个 $M\times M$ 的实对称矩阵，并且：

\[\rho(\mathbf r)=\sum_{\mu\nu}P_{\mu\nu}\chi_\mu(\mathbf r)\chi_\nu(\mathbf r)\]

对密度积分应当得到电子数：

\[\int\rho(\mathbf r)\,d\mathbf r=\sum_{\mu\nu}P_{\mu\nu}S_{\mu\nu}=\operatorname{Tr}(PS)=2\operatorname{Tr}(C_{\mathrm{occ}}^TSC_{\mathrm{occ}})=N\]

式中：

• $\operatorname{Tr}$：矩阵的迹，即把所有对角元素相加。
• $\operatorname{Tr}(PS)$：在这组有重叠的 AO 基函数下，由密度矩阵算出的电子总数。
• $\sum_{\mu\nu}$：让 $\mu,\nu$ 各自遍历全部 $M$ 个基函数后求和；后文带三个或四个 AO 下标的求和也采用这一约定。

同样，从占据轨道的正交归一性可以得到：

\[PSP=4C_{\mathrm{occ}}(C_{\mathrm{occ}}^TSC_{\mathrm{occ}})C_{\mathrm{occ}}^T=2P\]

这两个关系都很适合拿来检查程序中的密度矩阵。

5 从密度矩阵推导 RHF 能量

前面已经把电子密度中的占据轨道求和收集成了密度矩阵。接下来对能量作同样的展开：先将每个 MO 换成 AO 的线性组合，再找出其中的密度矩阵。这样就能用固定的 AO 积分和待优化的密度矩阵来计算能量。

5.1 单电子能量

定义 AO 单电子积分：

\[h_{\mu\nu}=\langle\chi_\mu\vert\hat h\vert\chi_\nu\rangle\]

式中：

• $h_{\mu\nu}$：单电子算符在第 $\mu$、$\nu$ 个基函数之间的矩阵元，即让算符作用于后一个基函数，再与前一个基函数相乘积分。
• $h$：由全部 $h_{\mu\nu}$ 组成的单电子积分矩阵，包含动能与核吸引两部分，也称芯哈密顿矩阵。

将轨道展开代入 $E_{\mathrm{1e}}$：

\[E_{\mathrm{1e}}=2\sum_i^{\mathrm{occ}}\sum_{\mu\nu}C_{\mu i}C_{\nu i}h_{\mu\nu}=\sum_{\mu\nu}\left(2\sum_i^{\mathrm{occ}}C_{\mu i}C_{\nu i}\right)h_{\mu\nu}\]

括号中正好是 $P_{\mu\nu}$，因此：

\[\boxed{E_{\mathrm{1e}}=\sum_{\mu\nu}P_{\mu\nu}h_{\mu\nu}}\]

5.2 双电子能量

定义 AO 双电子积分：

\[(\mu\nu\vert\lambda\sigma)=\iint\chi_\mu(\mathbf r_1)\chi_\nu(\mathbf r_1)\frac{1}{r_{12}}\chi_\lambda(\mathbf r_2)\chi_\sigma(\mathbf r_2)\,d\mathbf r_1d\mathbf r_2\]

式中：

• $(\mu\nu\vert\lambda\sigma)$：用四个 AO 基函数计算的双电子积分，记录两对基函数通过电子排斥算符产生的相互作用。
• $\chi_\mu\chi_\nu$、$\chi_\lambda\chi_\sigma$：分别取电子 $1$、电子 $2$ 坐标的一对基函数乘积；四个编号各自独立遍历全部基函数。

实基函数的乘法可交换，两个电子的积分变量也可互换，所以有：

\[(\mu\nu\vert\lambda\sigma)=(\nu\mu\vert\lambda\sigma)=(\mu\nu\vert\sigma\lambda)=(\lambda\sigma\vert\mu\nu)\]

先展开库仑积分：

\[J_{ij}=(ii\vert jj)=\sum_{\mu\nu\lambda\sigma}C_{\mu i}C_{\nu i}C_{\lambda j}C_{\sigma j}(\mu\nu\vert\lambda\sigma)\]

再展开交换积分。为了最后也能凑出 $P_{\mu\nu}P_{\lambda\sigma}$，这里把轨道 $i$ 的系数下标写成 $\mu,\nu$，轨道 $j$ 的系数下标写成 $\lambda,\sigma$：

\[K_{ij}=(ij\vert ji)=\sum_{\mu\nu\lambda\sigma}C_{\mu i}C_{\nu i}C_{\lambda j}C_{\sigma j}(\mu\lambda\vert\sigma\nu)=\sum_{\mu\nu\lambda\sigma}C_{\mu i}C_{\nu i}C_{\lambda j}C_{\sigma j}(\mu\lambda\vert\nu\sigma)\]

最后一个等号使用了同一电子坐标下两个实基函数可交换的性质。与库仑项相比，交换项中基函数在两个电子坐标上的配对变了。

代入 $E_{\mathrm{2e}}=\sum_{ij}^{\mathrm{occ}}(2J_{ij}-K_{ij})$，交换求和顺序：

\[E_{\mathrm{2e}}=\sum_{\mu\nu\lambda\sigma}\left(\sum_i^{\mathrm{occ}}C_{\mu i}C_{\nu i}\right)\left(\sum_j^{\mathrm{occ}}C_{\lambda j}C_{\sigma j}\right)[2(\mu\nu\vert\lambda\sigma)-(\mu\lambda\vert\nu\sigma)]\]

式中：

• $E_{\mathrm{2e}}$：总电子能量中的双电子部分，即电子与电子之间的排斥能在行列式波函数下的期望值，包含库仑与交换两部分。

由于密度矩阵的定义中包含 $2$，上面每个占据轨道求和都等于相应密度矩阵元素的一半：

\[\sum_i^{\mathrm{occ}}C_{\mu i}C_{\nu i}=\frac{1}{2}P_{\mu\nu}\]

两个求和因而共同带来 $1/4$：

\[E_{\mathrm{2e}}=\frac{1}{4}\sum_{\mu\nu\lambda\sigma}P_{\mu\nu}P_{\lambda\sigma}[2(\mu\nu\vert\lambda\sigma)-(\mu\lambda\vert\nu\sigma)]\]

将括号里的 $2$ 提出来：

\[\boxed{E_{\mathrm{2e}}=\frac{1}{2}\sum_{\mu\nu\lambda\sigma}P_{\mu\nu}P_{\lambda\sigma}\left[(\mu\nu\vert\lambda\sigma)-\frac{1}{2}(\mu\lambda\vert\nu\sigma)\right]}\]

所以，闭壳层 RHF 的电子能量可以完全写为：

\[E_{\mathrm{elec}}[P]=\sum_{\mu\nu}P_{\mu\nu}h_{\mu\nu}+\frac{1}{2}\sum_{\mu\nu\lambda\sigma}P_{\mu\nu}P_{\lambda\sigma}\left[(\mu\nu\vert\lambda\sigma)-\frac{1}{2}(\mu\lambda\vert\nu\sigma)\right]\]

式中：

• $E_{\mathrm{elec}}[P]$：由密度矩阵 $P$ 算出的电子能量；方括号 $[P]$ 表示这个量依赖整个密度矩阵。

给定核构型和基组以后，$h$ 与双电子积分都已经确定，接下来需要优化的是 $P$ 所对应的占据轨道。

6 Fock 矩阵的推导

能量写出来以后，下一步是寻找让它降低的轨道变化。先从密度矩阵入手：让它发生一点微小变化，看看能量会怎样变化。把能量的一阶变化按各个密度矩阵元素的变化量展开，前面的系数就会组成 Fock 矩阵。

为了看清求导时的系数，暂时把积分组合记为：

\[A_{\mu\nu,\lambda\sigma}=(\mu\nu\vert\lambda\sigma)-\frac{1}{2}(\mu\lambda\vert\nu\sigma)\]

式中：

• $A_{\mu\nu,\lambda\sigma}$：将一个库仑积分与半个交换积分合在一起的简写，便于写出后面的能量变化。

由实基函数双电子积分的对称性可知：

\[A_{\mu\nu,\lambda\sigma}=A_{\lambda\sigma,\mu\nu}\]

双电子能量包含两个 $P$，取一阶变化时两者都要变化：

\[\delta E_{\mathrm{2e}}=\frac{1}{2}\sum_{\mu\nu\lambda\sigma}\delta P_{\mu\nu}P_{\lambda\sigma}A_{\mu\nu,\lambda\sigma}+\frac{1}{2}\sum_{\mu\nu\lambda\sigma}P_{\mu\nu}\delta P_{\lambda\sigma}A_{\mu\nu,\lambda\sigma}\]

式中：

• $\delta P_{\mu\nu}$：密度矩阵元素发生的一点微小变化；$\delta P_{\lambda\sigma}$ 含义相同，只是元素编号不同。
• $\delta E_{\mathrm{2e}}$：由这些微小变化引起的双电子能量的一阶变化，只保留与变化量成正比的部分。
• 单独写在一个量前面的 $\delta$：表示对这个量取一阶变化；带两个编号的 $\delta_{ab}$ 则仍是前面定义的 Kronecker delta。

在第二项中交换哑指标 $(\mu,\nu)$ 与 $(\lambda,\sigma)$，利用 $A$ 的上述对称性，两项相同。因此：

\[\delta E_{\mathrm{2e}}=\sum_{\mu\nu}\left(\sum_{\lambda\sigma}P_{\lambda\sigma}A_{\mu\nu,\lambda\sigma}\right)\delta P_{\mu\nu}\]

加上单电子能量的变化：

\[\delta E_{\mathrm{elec}}=\sum_{\mu\nu}\left[h_{\mu\nu}+\sum_{\lambda\sigma}P_{\lambda\sigma}A_{\mu\nu,\lambda\sigma}\right]\delta P_{\mu\nu}\]

式中：

• $\delta E_{\mathrm{elec}}$：密度矩阵发生微小变化时，总电子能量相应的一阶变化。

将 $\delta P_{\mu\nu}$ 前的系数定义为 Fock 矩阵元素：

\[\boxed{F_{\mu\nu}=h_{\mu\nu}+\sum_{\lambda\sigma}P_{\lambda\sigma}\left[(\mu\nu\vert\lambda\sigma)-\frac{1}{2}(\mu\lambda\vert\nu\sigma)\right]}\]

式中：

• $F_{\mu\nu}$：电子能量对密度矩阵元素 $P_{\mu\nu}$ 的一阶变化系数。
• $F$：由这些元素组成的 Fock 矩阵，描述一个电子的动能、核吸引以及电子之间的平均库仑与交换作用。

于是可以简写为：

\[\delta E_{\mathrm{elec}}=\sum_{\mu\nu}F_{\mu\nu}\delta P_{\mu\nu}\]

这个推导说明，双电子能量前面的 $1/2$ 在求一阶变化时，被两个 $P$ 各贡献的一项抵消了。

接下来要在轨道正交归一的约束下继续求驻点。$P$ 由占据轨道构成，满足 $\operatorname{Tr}(PS)=N$ 和 $PSP=2P$ 等关系，因此实际的变分自由度仍取为轨道系数 $C_{\mathrm{occ}}$。

7 Roothaan 方程的推导

密度矩阵来自占据轨道，所以实际要调整的是轨道展开系数。调整时还要保证每个轨道归一化、不同轨道相互正交。下面先把能量对密度矩阵的变化换成对轨道系数的变化，再把这些约束一起纳入求导。

7.1 约束变分

由 $P_{\mu\nu}=2\sum_i^{\mathrm{occ}}C_{\mu i}C_{\nu i}$ 可得：

\[\frac{\partial P_{\mu\nu}}{\partial C_{\kappa i}}=2(\delta_{\mu\kappa}C_{\nu i}+\delta_{\nu\kappa}C_{\mu i})\]

式中：

• $C_{\kappa i}$：占据轨道 $i$ 在基函数 $\kappa$ 上的展开系数，是这里要改变的一个变量。
• $\partial P_{\mu\nu}/\partial C_{\kappa i}$：偏导数，表示只改变这个展开系数、其余系数暂时不变时，密度矩阵元素随它变化的快慢。
• $\delta_{\mu\kappa}$、$\delta_{\nu\kappa}$：对应两个 AO 编号的 Kronecker delta，用来选出被求导的那个系数。

利用链式法则以及上一节的一阶变化关系：

\[\frac{\partial E_{\mathrm{elec}}}{\partial C_{\kappa i}}=2\sum_\nu F_{\kappa\nu}C_{\nu i}+2\sum_\mu F_{\mu\kappa}C_{\mu i}\]

实轨道下 $F$ 为实对称矩阵，两项相同：

\[\frac{\partial E_{\mathrm{elec}}}{\partial C_{\kappa i}}=4\sum_\nu F_{\kappa\nu}C_{\nu i}\]

占据轨道的约束为：

\[(C_{\mathrm{occ}}^TSC_{\mathrm{occ}})_{ij}=\delta_{ij}\]

引入实对称的拉格朗日乘子矩阵 $\varepsilon$，构造：

\[\mathcal L=E_{\mathrm{elec}}-2\sum_{ij}^{\mathrm{occ}}\varepsilon_{ij}[(C_{\mathrm{occ}}^TSC_{\mathrm{occ}})_{ij}-\delta_{ij}]\]

式中：

• $\mathcal L$：把能量与正交归一约束合在一起的辅助函数，称为拉格朗日函数。
• $\varepsilon_{ij}$：约束对应的拉格朗日乘子，即给第 $(i,j)$ 项约束分配的一个待定系数；它们共同组成实对称矩阵 $\varepsilon$。
• $(C_{\mathrm{occ}}^TSC_{\mathrm{occ}})_{ij}$：占据轨道 $i,j$ 的重叠；矩阵括号外的下标 $(i,j)$ 表示取第 $i$ 行、第 $j$ 列的元素。
• $\delta_{ij}$：轨道 $i,j$ 应达到的重叠值，同一轨道为 $1$，不同轨道为 $0$。
• 方括号中的差：轨道重叠偏离正交归一要求的程度，满足约束时为 $0$。

约束项前的 $2$ 是乘子的定义约定，便于让最终的 $\varepsilon$ 直接对应通常的轨道能量。

约束项求导时，左右两个 $C_{\mathrm{occ}}$ 各贡献一项。利用 $S$ 与 $\varepsilon$ 的对称性，驻点条件为：

\[\frac{\partial\mathcal L}{\partial C_{\kappa i}}=4\sum_\nu F_{\kappa\nu}C_{\nu i}-4\sum_\nu\sum_j^{\mathrm{occ}}S_{\kappa\nu}C_{\nu j}\varepsilon_{ji}=0\]

这里第二项的 $\nu$ 遍历全部 AO 基函数，$j$ 遍历占据 MO。两边约去 $4$：

\[\sum_\nu F_{\kappa\nu}C_{\nu i}=\sum_\nu\sum_j^{\mathrm{occ}}S_{\kappa\nu}C_{\nu j}\varepsilon_{ji}\]

写成矩阵形式：

\[FC_{\mathrm{occ}}=SC_{\mathrm{occ}}\varepsilon\]

7.2 正则轨道

目前的 $\varepsilon$ 是一个实对称矩阵，可以通过正交矩阵 $U$ 对角化：

\[U^T\varepsilon U=\varepsilon_{\mathrm{diag}},\qquad\widetilde C_{\mathrm{occ}}=C_{\mathrm{occ}}U\]

式中：

• $U$：将原有占据轨道重新线性组合的正交矩阵，满足 $U^TU=UU^T=I$，因此变换后轨道仍正交归一。
• $\varepsilon_{\mathrm{diag}}$：由乘子矩阵的本征值组成的对角矩阵；下标 $\mathrm{diag}$ 表示只有对角位置可以非零。
• $\widetilde C_{\mathrm{occ}}$：重新组合后的占据轨道系数矩阵；上方的波浪号表示经过这次轨道变换。

在占据轨道之间作这样的正交变换，密度矩阵不变：

\[\widetilde P=2C_{\mathrm{occ}}UU^TC_{\mathrm{occ}}^T=P\]

式中：

• $\widetilde P$：由变换后的占据轨道算出的密度矩阵，它与变换前的 $P$ 相同。

所以电子密度、能量和由 $P$ 构造的 Fock 矩阵都不变。在变换后的轨道下：

\[F\widetilde C_{\mathrm{occ}}=S\widetilde C_{\mathrm{occ}}\varepsilon_{\mathrm{diag}}\]

这种使乘子矩阵对角化的轨道称为正则轨道。给定 $P$ 后，可以求解整个 $M$ 维 AO 空间中的广义本征值问题，得到占据轨道与虚轨道，统一写成：

\[\boxed{FC=SC\varepsilon}\]

式中：

• $C$：包含全部 $M$ 个正则空间轨道的系数矩阵，其中部分列为占据轨道，其余列为虚轨道。
• $\varepsilon_p$：第 $p$ 个正则轨道的轨道能量，也就是该轨道对应的广义本征值。
• $\varepsilon=\operatorname{diag}(\varepsilon_1,\ldots,\varepsilon_M)$：按各个轨道能量排成的对角矩阵；$\operatorname{diag}$ 表示把括号内的数依次放在矩阵对角位置，其余位置填 $0$。

这就是通常使用的 Hartree–Fock–Roothaan 方程。Roothaan 的原始论文讨论了将分子轨道展开为已知基函数后得到的自洽矩阵方程。

8 Fock 算符与 HF 方程

开头说过，HF 希望让每个电子在由占据轨道决定的有效作用下运动。上面的变分已经给出了轨道方程和 Fock 矩阵，现在可以把它们重新写成作用于轨道函数的算符，看看这个有效作用具体包含什么。将 $F_{\mu\nu}$ 写成算符在 AO 基组中的矩阵元：

\[F_{\mu\nu}=\langle\chi_\mu\vert\hat f\vert\chi_\nu\rangle\]

式中：

• $\hat f$：Fock 算符，是把上面的有效单电子能量关系写成作用于轨道函数的形式；$F$ 则是它在 AO 基函数中的矩阵表示。

对于闭壳层 RHF，这个 Fock 算符为：

\[\hat f=\hat h+\sum_j^{\mathrm{occ}}(2\hat J_j-\hat K_j)\]

式中：

• $\hat J_j$：占据空间轨道 $j$ 产生的库仑算符，用这个轨道的电子分布所产生的静电势作用于另一个轨道。
• $\hat K_j$：占据空间轨道 $j$ 对应的交换算符，计算由轨道交换产生的那部分作用。
• $2\hat J_j-\hat K_j$：轨道 $j$ 双占据时对 Fock 算符的贡献，库仑作用计两次，同自旋交换作用计一次。

库仑算符作用在任意轨道 $\phi_p$ 上时：

\[\hat J_j\phi_p(\mathbf r_1)=\left[\int\frac{\phi_j(\mathbf r_2)^2}{r_{12}}\,d\mathbf r_2\right]\phi_p(\mathbf r_1)\]

方括号中是轨道 $j$ 的电子分布产生的库仑势，它在 $\mathbf r_1$ 处乘上待作用的轨道。

交换算符的作用为：

\[\hat K_j\phi_p(\mathbf r_1)=\phi_j(\mathbf r_1)\int\frac{\phi_j(\mathbf r_2)\phi_p(\mathbf r_2)}{r_{12}}\,d\mathbf r_2\]

这里 $\phi_p$ 出现在积分内部，结果依赖它在其他位置的取值，因此交换算符是非局域的。算符中的 $2J-K$，与前面自旋积分得到的计数一致。

把 $P_{\lambda\sigma}=2\sum_j^{\mathrm{occ}}C_{\lambda j}C_{\sigma j}$ 代入 Fock 矩阵，就能直接验证 $F_{\mu\nu}=\langle\chi_\mu\vert\hat f\vert\chi_\nu\rangle$。例如交换部分的系数 $-P/2$ 恰好给出对占据空间轨道的 $-\sum_j\hat K_j$。

在不限制轨道展开空间时，正则 HF 方程写作：

\[\hat f\phi_p=\varepsilon_p\phi_p\]

将轨道近似展开为 $\phi_p=\sum_\nu C_{\nu p}\chi_\nu$，要求方程残差与每个基函数 $\chi_\mu$ 正交，即：

\[\langle\chi_\mu\vert(\hat f-\varepsilon_p)\vert\phi_p\rangle=0\]

代入展开式后得到：

\[\sum_\nu\langle\chi_\mu\vert\hat f\vert\chi_\nu\rangle C_{\nu p}=\varepsilon_p\sum_\nu\langle\chi_\mu\vert\chi_\nu\rangle C_{\nu p}\]

这就回到了 $FC=SC\varepsilon$。有限基组求解的是 HF 方程在所选基函数空间内的投影；扩大基组可以改进轨道的表示，单行列式的变分限制则仍然保留。

9 基组正交化

普通的矩阵本征值问题可以直接通过对角化求解，而这里的轨道方程右侧还有重叠矩阵。可以先将原有基函数重新组合成一组正交归一的函数，使重叠矩阵变成单位矩阵，再求轨道，最后把系数换回原来的 AO 表示。

对于线性无关的基函数，$S$ 是实对称正定矩阵。先对它作本征分解：

\[S=UsU^T\]

式中：

• $U$：由重叠矩阵 $S$ 的归一化本征向量按列组成的正交矩阵，维数为 $M\times M$。
• $s$：由 $S$ 的本征值组成的对角矩阵；本征值表示沿相应本征向量方向的缩放倍数，$S$ 正定时这些值全部为正。

定义对称正交化矩阵：

\[X=S^{-1/2}=Us^{-1/2}U^T\]

式中：

• $X$：对称正交化矩阵，用来把原先有重叠的基函数组合成一组正交归一的基函数。
• $s^{-1/2}$：对角矩阵，对 $s$ 的每个对角元素取平方根后再取倒数。
• $S^{-1/2}$：通过上述本征分解定义的矩阵逆平方根，作用是补偿重叠矩阵在各本征方向上的缩放。

于是：

\[X^TSX=I\]

令 $C=XC’$，将它代入 $FC=SC\varepsilon$，再左乘 $X^T$：

\[X^TFXC'=X^TSXC'\varepsilon=C'\varepsilon\]

式中：

• $C’$：轨道在正交化基函数中的展开系数矩阵；撇号 $’$ 表示换到了这组基函数下，$C=XC’$ 则将它变回原来的 AO 表示。

因此只需对角化：

\[F'C'=C'\varepsilon,\qquad F'=X^TFX\]

式中：

• $F’$：Fock 算符在正交化基函数中的矩阵；与 $F$ 描述的是同一个算符，只是采用了不同的基函数。

得到 $C’$ 后，通过 $C=XC’$ 变回原来的 AO 表示。此时：

\[C^TSC=C'^TX^TSXC'=C'^TC'=I\]

如果基函数接近线性相关，$S$ 会有很小的本征值，$s^{-1/2}$ 会放大数值误差，需要先处理相应的线性相关问题。下面的教学代码假设 $S$ 正定且条件良好。

10 SCF 迭代

给定密度矩阵以后，现在已经知道怎样构造 Fock 矩阵、求出轨道，并由占据轨道算出新的密度矩阵。这也回到了开头平均场的相互依赖关系：用来求轨道的有效作用，正是由这些轨道自身决定的。因此，要让输入与输出相互一致，需要重复下面的过程：

\[P^{(k)}\;\longrightarrow\;F[P^{(k)}]\;\longrightarrow\;C^{(k+1)}\;\longrightarrow\;P^{(k+1)}\]

式中：

• 上标 $(k)$、$(k+1)$：第 $k$ 轮与下一轮 SCF 迭代的编号，括号中的数表示迭代次数。
• $P^{(k)}$：第 $k$ 轮使用的密度矩阵。
• $F[P^{(k)}]$：由这轮密度矩阵构造的 Fock 矩阵；方括号表示它依赖输入的密度矩阵。
• $C^{(k+1)}$：对这轮 Fock 矩阵求解轨道方程后得到的新轨道系数。
• $P^{(k+1)}$：由这些新轨道的占据部分构造的新密度矩阵。

于是可以从一个初猜出发反复迭代，直到输入密度与输出密度一致，这就是自洽场（self-consistent field，SCF）。对于通常按能量从低到高占据的闭壳层基态计算，每轮取最低的 $N/2$ 个空间轨道构造新密度。

一种简单的启动方式是先求解 $hC=SC\varepsilon$，用芯哈密顿矩阵的轨道构造初始密度。程序中也可以从临时矩阵 $P^{(0)}=0$ 开始，此时第一次构造的 $F$ 就是 $h$；对角化后即得到具有正确电子数的轨道密度。

10.1 能量计算

记 $G[P]=F[P]-h$。前面的双电子能量可以写为：

\[E_{\mathrm{2e}}=\frac{1}{2}\sum_{\mu\nu}P_{\mu\nu}G_{\mu\nu}[P]\]

式中：

• $G[P]=F[P]-h$：Fock 矩阵中的电子间相互作用部分，去掉了单电子的动能和核吸引项。
• $G_{\mu\nu}[P]$：上述矩阵的第 $(\mu,\nu)$ 个元素，同样由密度矩阵 $P$ 决定。

加上单电子部分：

\[\boxed{E_{\mathrm{elec}}[P]=\frac{1}{2}\sum_{\mu\nu}P_{\mu\nu}\{h_{\mu\nu}+F_{\mu\nu}[P]\}}\]

这里 $P$ 和 $F[P]$ 必须来自同一个密度。如果在某一轮已经更新到 $P^{(k+1)}$，就应当用它重新构造 $F[P^{(k+1)}]$，再按上式计算该密度对应的电子能量。

在自洽的正则轨道下，$\varepsilon_i=h_{ii}+\sum_j^{\mathrm{occ}}(2J_{ij}-K_{ij})$。对占据轨道能量求和得到：

\[2\sum_i^{\mathrm{occ}}\varepsilon_i=E_{\mathrm{1e}}+2E_{\mathrm{2e}}\]

因此电子能量也可以写为：

\[E_{\mathrm{elec}}=2\sum_i^{\mathrm{occ}}\varepsilon_i-E_{\mathrm{2e}}\]

轨道能量求和把电子对的相互作用算了两次，需要减去一次。以只有一个双占据空间轨道的二电子体系为例，$\varepsilon_1=h_{11}+J_{11}$，而 $E_{\mathrm{elec}}=2h_{11}+J_{11}$，两者的区别可以直接看出来。

10.2 程序实现

下面保留显式求和，便于对照公式。假定已经有重叠矩阵 S、单电子积分 h、完整的双电子积分数组 eri、电子数 nelec 和核排斥能 E_nuc，并且 eri[mu, nu, lam, sig] 按化学家记号存储。

import numpy as np


def build_fock(P, h, eri):
    M = h.shape[0]
    F = h.copy()
    for mu in range(M):
        for nu in range(M):
            for lam in range(M):
                for sig in range(M):
                    F[mu, nu] += P[lam, sig] * (
                        eri[mu, nu, lam, sig]
                        - 0.5 * eri[mu, lam, nu, sig]
                    )
    return F


def rhf(S, h, eri, nelec, E_nuc, max_cycle=100):
    M = h.shape[0]
    if nelec % 2 != 0 or not 0 < nelec <= 2 * M:
        raise ValueError("RHF requires an even electron count in 1..2M")
    nocc = nelec // 2

    s, U = np.linalg.eigh(S)
    if s.min() <= 1e-10:
        raise ValueError("Handle linear dependence before orthogonalizing")
    X = (U * s**-0.5) @ U.T

    P = np.zeros_like(h)
    E_old = None
    for cycle in range(1, max_cycle + 1):
        F = build_fock(P, h, eri)
        eps, C_prime = np.linalg.eigh(X.T @ F @ X)
        C = X @ C_prime
        C_occ = C[:, :nocc]
        P_new = 2 * C_occ @ C_occ.T

        F_new = build_fock(P_new, h, eri)
        E_elec = 0.5 * np.sum(P_new * (h + F_new))
        E_total = E_elec + E_nuc
        dP = np.linalg.norm(P_new - P)
        residual = F_new @ P_new @ S - S @ P_new @ F_new
        rnorm = np.linalg.norm(residual)
        dE = np.inf if E_old is None else abs(E_total - E_old)
        print(f"{cycle:3d}  E = {E_total:.12f}  "
              f"dE = {dE:.3e}  dP = {dP:.3e}  R = {rnorm:.3e}")

        if dE < 1e-10 and dP < 1e-8 and rnorm < 1e-8:
            return E_total, P_new
        P, E_old = P_new, E_total

    raise RuntimeError("SCF did not converge")

其中最核心的几行分别对应：

• C = X @ C_prime：把正交基中的轨道系数变回 AO 表示。
• P_new = 2 * C_occ @ C_occ.T：由占据轨道构造闭壳层密度矩阵。
• build_fock(...)：将密度矩阵与库仑、交换积分收缩，得到 $F[P]$。
• 0.5 * np.sum(P_new * (h + F_new))：计算同一个密度对应的电子能量。

除能量变化和密度变化外，代码还检查了自洽残差：

\[R=FPS-SPF\]

式中：

• $R$：自洽残差矩阵，用来检查当前密度矩阵与它产生的 Fock 矩阵是否满足轨道方程要求的相容关系；自洽时各元素应接近 $0$。
• 此处的 $F$：由当前 $P$ 构造的 $F[P]$，与能量公式采用相同的密度。

在自洽解处，将 $P=2C_{\mathrm{occ}}C_{\mathrm{occ}}^T$ 和占据轨道的 Roothaan 方程代入，就能得到 $R=0$。每轮由正交归一轨道构造的 $P$ 还应满足 $\operatorname{Tr}(PS)=N$ 与 $PSP=2P$。

这段程序用于对应推导，四重循环构造 Fock 矩阵的代价随基函数数目按 $M^4$ 增长。直接迭代也可能出现振荡或不收敛，实际程序通常还会使用阻尼、DIIS 等收敛技术。SCF 达到上述驻点条件以后，是否为所需的最低能量解，还涉及初猜、占据选择以及轨道稳定性。

参考资料

本文按 AO 展开、密度矩阵、能量和轨道变分的顺序整理。与 Szabo–Ostlund 对照阅读时，可以先看第 2 章的行列式矩阵元，再看第 3.2 节的 HF 变分推导和第 3.4 节的闭壳层 Roothaan 方程。

• Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory. Dover Publications, 1996. 第 2、3 章。
• 《现代量子化学》中文翻译项目：本文整理时对照了其中的单行列式能量、闭壳层自旋积分与 Roothaan 方程推导。
• Roothaan, C. C. J. New Developments in Molecular Orbital Theory. Rev. Mod. Phys. 1951, 23, 69–89。