Study note：量子化学方法

笔者是初学者，理解可能有误，以下内容请慎重相信。

1 Hartree-Fock方法 (HF)

Hartree-Fock方法用一个Slater行列式近似多电子波函数。行列式的反对称性保证了泡利不相容原理，并自然产生交换作用；电子之间其余的瞬时相关并未被纳入。各个轨道又依赖于其他电子形成的平均场，因此需要通过自洽场（Self-Consistent Field, SCF）迭代求解。

以下以偶数电子闭壳层体系的限制性Hartree-Fock方法（Restricted Hartree-Fock, RHF）为例。每个占据空间轨道容纳一对自旋相反的电子。

1.1 单行列式近似

HF波函数写成一组自旋轨道构成的Slater行列式：

\[\Psi_{\mathrm{HF}} = \frac{1}{\sqrt{N!}} \det\left[\psi_i(\mathbf{x}_j)\right],\]

其中，$N$为电子数，$\mathbf{x}$同时包含空间坐标和自旋坐标。RHF中两个自旋相反的电子共享同一个空间轨道$\phi_i$。

1.2 用原子轨道展开分子轨道

实际计算只能使用有限组基函数。将分子轨道写成原子轨道（Atomic Orbital, AO）基函数的线性组合：

\[\phi_i=\sum_\mu C_{\mu i}\chi_\mu,\]

其中，$\chi_\mu$为AO基函数，$C_{\mu i}$为分子轨道系数。AO之间通常并不正交，因此定义重叠矩阵

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

分子轨道满足正交归一条件：

\[\langle\phi_i|\phi_j\rangle=\delta_{ij},\]

代入AO展开后可得

\[C^TSC=I.\]

1.3 密度矩阵

对闭壳层体系，电子密度为

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

将分子轨道展开到AO基组中：

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

其中AO密度矩阵定义为

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

系数2来自空间轨道的双占据。密度矩阵还满足电子数检查：

\[\operatorname{Tr}(PS)=N_{\mathrm e}.\]

1.4 单电子积分、双电子积分与HF能量

单电子哈密顿算符包含电子动能和电子-核吸引：

\[\hat h = -\frac{1}{2}\nabla^2 - \sum_A\frac{Z_A}{r_A},\]

相应的单电子积分为

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

电子-电子排斥由双电子积分表示：

\[(\mu\nu|\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.\]

采用上述密度矩阵约定，RHF电子能量可以写成

\[E_{\mathrm{elec}} = \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|\lambda\sigma) - \frac{1}{2}(\mu\lambda|\nu\sigma) \right].\]

方括号中的第一项为库仑项，描述电子密度之间的经典排斥；第二项为交换项，来自Slater行列式的反对称性。体系总能量还需要加上核间排斥能：

\[E_{\mathrm{HF}}=E_{\mathrm{elec}}+E_{\mathrm{nuc}}.\]

1.5 Fock矩阵

对电子能量关于密度矩阵求导，可以得到Fock矩阵：

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

记电子-电子部分为$G(P)$，则

\[F=h+G(P).\]

电子能量也可以简写为

\[E_{\mathrm{elec}} = \frac{1}{2} \sum_{\mu\nu}P_{\mu\nu}\left(h_{\mu\nu}+F_{\mu\nu}\right).\]

1.6 Hartree-Fock-Roothaan方程

HF方法在$C^TSC=I$的约束下最小化电子能量。引入拉格朗日乘子并令能量对轨道系数的一阶导数为零，可得

\[\sum_\nu F_{\mu\nu}C_{\nu i} = \sum_\nu S_{\mu\nu}C_{\nu i}\varepsilon_i.\]

写成矩阵形式为

\[FC=SC\varepsilon.\]

这就是Hartree-Fock-Roothaan方程。它是广义本征值问题：$\varepsilon$给出轨道能量，$C$给出分子轨道系数。由于$F$依赖$P$，而$P$又由$C$构造，该方程不能一次求完，只能迭代至自洽。

1.7 SCF迭代

一次RHF计算可以按以下步骤进行：

1. 读入重叠矩阵$S$、单电子积分$h$、双电子积分$(\mu\nu

   \lambda\sigma)$、核间排斥能$E_{\mathrm{nuc}}$和电子数，并给出初始密度矩阵$P$。

2. 根据当前密度矩阵构造Fock矩阵$F(P)$。
3. 求解$FC=SC\varepsilon$，并按照轨道能量从低到高排列分子轨道。
4. 取最低的$n_{\mathrm{occ}}=N_{\mathrm e}/2$个轨道作为占据轨道，构造新的密度矩阵$P_{\mathrm{new}}$。
5. 用$P_{\mathrm{new}}$计算电子能量和总能量，并检查$\operatorname{Tr}(P_{\mathrm{new}}S)$是否等于电子数。
6. 比较相邻两次迭代的能量和密度矩阵；若尚未收敛，则令$P\leftarrow P_{\mathrm{new}}$并继续迭代。

若初始密度矩阵取零矩阵，第一次迭代有$F=h$，相当于从core Hamiltonian guess开始。

下面是一段与上述公式逐项对应的最小RHF实现。为了便于观察公式，Fock矩阵和密度矩阵都使用显式循环构造：

import numpy as np
from scipy.linalg import eigh


def build_fock(h, eri, P):
    """F_mn = h_mn + sum_ls P_ls [(mn|ls) - 1/2 (ml|ns)]."""
    K = h.shape[0]
    F = np.zeros_like(h)

    for mu in range(K):
        for nu in range(K):
            value = h[mu, nu]
            for lam in range(K):
                for sig in range(K):
                    coulomb = eri[mu, nu, lam, sig]
                    exchange = eri[mu, lam, nu, sig]
                    value += P[lam, sig] * (coulomb - 0.5 * exchange)
            F[mu, nu] = value

    return F


def build_density_matrix(C_occ):
    """P_mn = 2 sum_i^occ C_mi C_ni."""
    K, nocc = C_occ.shape
    P = np.zeros((K, K))

    for mu in range(K):
        for nu in range(K):
            for i in range(nocc):
                P[mu, nu] += 2.0 * C_occ[mu, i] * C_occ[nu, i]

    return P


def compute_electronic_energy(P, h, F):
    """E_elec = 1/2 sum_mn P_mn (h_mn + F_mn)."""
    return 0.5 * np.sum(P * (h + F))


nocc = nelec // 2
P = np.zeros_like(h)
E_old = None

for cycle in range(50):
    # F(P) -> C -> P_new
    F = build_fock(h, eri, P)
    eps, C = eigh(F, S)
    C_occ = C[:, :nocc]
    P_new = build_density_matrix(C_occ)

    # 能量必须使用与P_new对应的Fock矩阵
    F_new = build_fock(h, eri, P_new)
    E_elec = compute_electronic_energy(P_new, h, F_new)
    E_total = E_elec + E_nuc

    print(f"cycle = {cycle + 1}")
    print("Tr(PS) =", np.trace(P_new @ S))
    print("E_total =", E_total)

    if E_old is not None:
        dE = abs(E_total - E_old)
        dP = np.linalg.norm(P_new - P)
        print("dE =", dE, "dP =", dP)

        if dE < 1e-10 and dP < 1e-8:
            print("SCF converged")
            break

    P = P_new
    E_old = E_total

HF通过交换项处理了同自旋电子的交换效应，但仍然缺少单行列式之外的电子相关。因此，HF通常既是一个独立的近似方法，也是MP、CI和CC等post-HF方法的参考态。

2 post HF方法

HF不擅长描述电子相关作用，因此后来基于HF波函数又开发了许多方法。

---

微扰理论 (Møller-Plesset Perturbation Theory, MP)

Møller-Plesset微扰理论的目标是通过对电子波函数和能量的逐阶修正来描述电子相关效应，核心思想为将哈密顿算符中难以求解的电子相关部分提取出来，作为微扰项处理：

\[\hat{H} = \hat{H_0} + \lambda\hat{H'}\]

其中：

• $\hat{H_0}$为HF哈密顿量

• $\hat{H’}$为微扰项，它可以表示为：

  \[\hat{H}' = \hat{H}_{\text{ee}} - \hat{V}_{\text{HF}},\]

  其中：

  • $\hat{H}_{\text{ee}}$ 是电子之间的精确相互作用；
  • $\hat{V}_{\text{HF}}$ 是 HF 的平均场近似下的电子相互作用。

使用微扰参数 $\lambda$ 将总能量和波函数展开为级数：

\[E = E_0 + \lambda E_1 + \lambda^2 E_2 + \lambda^3 E_3 + \cdots,\] \[|\Psi\rangle = |\Psi_0\rangle + \lambda |\Psi_1\rangle + \lambda^2 |\Psi_2\rangle + \cdots.\]

其中：

• 零阶项$E_0$为HF能量，即 $\langle\Psi_0|\hat{H}_0|\Psi_0\rangle$
• 由于HF波函数使用了变分法，一阶项$E_1$为0。

• 二阶项$E_2$可以写为：

  \[E_2 = \sum_{i \neq 0} \frac{|\langle \Psi_0 | \hat{H}' | \Psi_i \rangle|^2}{E_0 - E_i}\]

  在实际计算中，二阶能量校正的表达式可以写为：

  \[E_2 = \sum_{ij}^{\text{occupied}} \sum_{ab}^{\text{virtual}} \frac{| \langle ij | ab \rangle - \langle ij | ba \rangle |^2}{\epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b}\]

  其中：

  • $i, j$：占据轨道（HF 波函数中被填满的轨道）。
  • $a, b$：虚轨道（HF 波函数中未被填满的轨道）。
  • $\epsilon_i, \epsilon_j, \epsilon_a, \epsilon_b$：HF 轨道能量。
  • $\langle ij | ab \rangle$：双电子积分，表示两个电子在占据轨道 $i, j$ 和虚轨道 $a, b$ 之间的相互作用。

最常用的微扰方法为MP2，更高阶的方法由于形式复杂，计算代价较大，并不常用。

---

组态相互作用 (Configuration Interaction, CI)

CI是一种多电子波函数近似方法，核心思想是通过线性组合多个单电子波函数来表示体系的真实波函数，以避免使用平均场近似处理电子相关问题。

就像用多个共振式共同描述体系的真实结构一样，体系的真实波函数也可以被表示为多个Slater行列式的线性组合：

\[\Psi_{\text{CI}} = c_0 \Phi_0 + \sum_i c_i \Phi_i + \sum_{i<j} c_{ij} \Phi_{ij} + \cdots\]

其中：

• $\Phi_0$ 是参考态波函数，通常为HF波函数。
• $\Phi_i, \Phi_{ij}, \dots$ 是通过单电子、双电子等激发生成的激发组态。
• $c_0, c_i, c_{ij}, \dots$ 是各个组态的系数。

CI的目标是通过变分原理最小化总能量，因此需要构造哈密顿矩阵 $H_{\text{CI}}$，矩阵元定义为： \(H_{ij} = \langle \Phi_i | \hat{H} | \Phi_j \rangle,\) 其中：$\hat{H}$ 是全电子哈密顿算符；$\Phi_i$ 和 $\Phi_j$ 是不同组态的 Slater 行列式。

随后通过对CI哈密顿矩阵的对角化，可以得到体系的基态和激发态波函数与能量。然而，完全考虑所有组态时，计算复杂度将达到$O(N!)$，计算量随体系增大而爆炸式增长，仅适用于极小体系。因此，实际使用的CI多数为截断CI，如仅包含一阶项的CIS、考虑到二阶项的CISD等。

---

耦合簇 (Coupled Cluster, CC)

为了避免显示计算CI的行列式系数，提出了Coupled Cluster方法。该方法的核心思想是通过指数形式的簇算符$e^{\hat{T}}$描述体系激发，并通过指数函数的泰勒展开将高阶激发纳入考虑：

\(|\Psi_{\text{CC}}\rangle = e^{\hat{T}} |\Phi_0\rangle,\) 其中：

• $\Phi_0$ 是参考态波函数，通常为HF波函数。

• $\hat{T}$ 是激发算符，定义为： \(\hat{T} = \hat{T}_1 + \hat{T}_2 + \hat{T}_3 + \dots,\)

  $\hat{T}_1, \hat{T}_2$ 等算符用于描述不同阶的电子跃迁。比如：

  \[\hat{T}_1 = \sum_{ia} t_i^a a_a^\dagger a_i\] \[\quad \hat{T}_2 = \frac{1}{4} \sum_{ijab} t_{ij}^{ab} a_a^\dagger a_b^\dagger a_j a_i,\]

  其中：

  • $t_i^a$ 和 $t_{ij}^{ab}$ 为激发振幅，可以调整激发态的组态系数。
  • $a_a^\dagger, a_i$ 为创建和湮灭算符，分别对应占据和虚轨道。

体系总能量为：

\[E_{\text{CC}} = \langle \Phi_0 | \hat{H} e^{\hat{T}} | \Phi_0 \rangle.\]

将 Schrödinger 方程投影到不同的组态空间，得到一系列关于激发振幅 $t_i^a, t_{ij}^{ab}$ 等的非线性方程组：

\[\langle \Phi_i^a | \hat{H} e^{\hat{T}} | \Phi_0 \rangle = 0, \quad \langle \Phi_{ij}^{ab} | \hat{H} e^{\hat{T}} | \Phi_0 \rangle = 0, \dots\]

最后通过数值迭代求解这些振幅方程，求得正确的振幅值后，该等式将会成立，就可以得到体系能量了。

CC的优点是避免了显示计算组态系数，通过指数的泰勒展开和激发振幅来构造了等效的表示方法，只要考虑到二阶项就能得到很高精度的结果。缺点是需要解非线性方程，可能带来数值收敛性问题。

---

3 多参考组态自洽场 (Multi-Configurational Self-Consistent Field, MCSCF)

上述方法着眼于解决动态相关问题，然而除高阶CI和高阶CC以外，这些方法对静态相关的考虑比较欠缺。因此当HF波函数对体系描述存在定性错误时，上述方法除高阶CI和高阶CC外均会失效。MCSCF方法是一种针对多电子体系的量子化学方法，用于精确描述静态相关效应，核心思想是直接优化多电子波函数，从而更好地描述静态相关效应。

---

完整活性空间自洽场(Complete Active Space Self-Consistent Field, CASSCF)

CASSCF方法是MCSCF的经典实现之一，广泛用于研究分子基态和低激发态的电子结构。其核心思想是定义一个完整活性空间(Complete Active Space, CAS)，在活性空间内进行FCI处理，考虑所有电子组态：

\[\Psi_{\text{CASSCF}} = \sum_i c_i \Phi_i,\]

其中：

• $\Phi_i$ 是活性空间内不同电子组态的Slater行列式；
• $c_i$ 是每个组态的线性组合系数。

体系能量为：

\[E_{\text{CASSCF}} = \langle \Psi_{\text{CASSCF}} | \hat{H} | \Psi_{\text{CASSCF}} \rangle + V_\text{HF}\]

其中：$V_\text{HF}$是活性空间外的电子在平均场近似下的能量。

交替优化$c_i$和$E_{\text{CASSCF}}$，直到二者均收敛，就能得到CASSCF级别的轨道和能量。

---

完整活性空间二阶微扰理论(Complete Active Space Perturbation Theory 2nd Order, CASPT2) 与N电子价态二阶微扰理论(N-Electron Valence Perturbation Theory 2nd Order, NEVPT2)

CASSCF对活性空间外的电子仍然使用平均场近似，因此对动态相关描述欠缺。CASPT2和NEVPT2是基于CASSCF波函数的二阶微扰理论方法，核心思想是通过MP2描述动态相关，从而得到更准确的能量。按照微扰理论，以CASSCF波函数为参考波函数将哈密顿量分解为零级哈密顿量和微扰哈密顿量：

\[\hat{H} = \hat{H}_0 + \hat{H}'.\]

其中，动态相关的贡献通过二阶微扰能量校正计算得到： \(E^{(2)} = \sum_{k \neq 0} \frac{\langle \Psi_0 | \hat{H}' | \Psi_k \rangle^2}{E_0 - E_k}.\)

CASPT2和NVEPT2的区别在于NEVPT2对波函数和零阶哈密顿量的定义更加严格，所有修正能量满足正定性条件，结果往往更加准确。然而，NEVPT2的形式更复杂，计算效率不如CASPT2。

MRCI、MRCC

类似地，以CASSCF波函数作为参考态，按照前述原理进行CI和CC计算，可以同时精确描述动态相关和静态相关效应。缺点是计算量极大，只能研究极小体系。

4 密度泛函理论(Density Functional Theory, DFT)

DFT是一种基于电子密度而非波函数的量子力学方法，起源于第一性原理，在量子化学方面同样表现出色。其核心思想是通过电子密度$\rho(\mathbf{r})$而非复杂的多电子波函数来描述多电子体系。

a. Hohenberg-Kohn定理

在外势 $v(\mathbf{r})$ 给定的条件下，体系的基态电子密度 $\rho(\mathbf{r})$ 可唯一确定体系的外势 $v(\mathbf{r})$ 和哈密顿量 $\hat{H}$，进而确定体系的基态波函数 $\Psi_0$ 和能量 $E_0$，即：$\rho(\mathbf{r}) \rightarrow v(\mathbf{r}) \rightarrow \Psi_0 \rightarrow E_0$。在HK定理下，3N维的多电子波函数可以被简化为3维的电子密度泛函，大大降低了计算复杂度。

b. Kohn-Sham方程

虽然经过简化，直接处理多电子体系的相互作用密度仍然非常复杂。Kohn和Sham提出了将问题转化为一组单电子问题的方法。首先，将真实的相互作用体系映射到一个无相互作用的虚拟体系，通过引入单电子轨道 $\phi_i(\mathbf{r})$ 表达电子密度：

\[\rho(\mathbf{r}) = \sum_i |\phi_i(\mathbf{r})|^2.\]

然后通过以下Kohn-Sham方程计算单电子轨道：

\[\left[ -\frac{\hbar^2}{2m} \nabla^2 + v_{\text{eff}}(\mathbf{r}) \right] \phi_i(\mathbf{r}) = \epsilon_i \phi_i(\mathbf{r}),\]

其中$v_{\text{eff}}(\mathbf{r})$ 是有效外势，包括库仑势和交换-相关势。库伦势由电子密度产生的电子-电子相互作用，可以精确求解；而交换-相关势是DFT最重要的部分，描述量子力学效应和电子相关。精确的交换-相关泛函形式未知，实际常用各种近似方法构造交换-相关泛函：

• 局域密度近似(Local Density Approximation, LDA)假设在均匀电子气模型下，交换-相关能仅依赖于局域的电子密度。LDA泛函的优点是形式简单，计算高效，缺点是对密度梯度大的体系(非金属体系)描述不够准确。例：SVWN

• 广义梯度近似(Generalized Gradient Approximation, GGA)虑了密度梯度对交换-相关能的修正，可以分别考虑交换部分和相关部分，提供更精确的结果。例：
  • 交换泛函：B88 (Becke 1988)等。
  • 相关泛函：LYP (Lee-Yang-Parr)等。
  • 交换-相关泛函：BLYP、PBE等。
• 元广义梯度近似(meta-Generalized Gradient Approximation, meta-GGA)引入了电子动能密度或Laplacian等额外信息，可以更加细致地刻画复杂的电子相互作用。例：M06-L，TPSS等。

c. 杂化(Hybrid)

为了改进泛函性能，可以在交换相关项$E_{\text{XC}}$中引入其他理论方法。

• HF通过单电子波函数的反对称化Slater行列式反映电子间排斥的非局域性质，能够提供比DFT更准确的交换能量。使用一部分HF交换项代替DFT交换项，可以得到杂化泛函。例：B3LYP等。
• post-HF方法，如MP2等对电子相关描述较好，在杂化泛函的基础上引入MP2相关项，可以得到双杂化泛函。例：B2PLYP等。

构造出泛函后，可以像HF的自洽场方法一样迭代求解KS方程，直到电子密度收敛，就能得到DFT级别的能量。

5. 对比

方法	静态相关	动态相关	备注
HF	无法考虑	无法考虑
DFT	可以考虑	描述较好	静态相关强的问题应使用纯泛函
单参考CI	几乎没有考虑	描述较好	高阶CI一定程度上能解决静态相关问题
微扰	几乎没有考虑	可以考虑
单参考CC	几乎没有考虑	描述较好	高阶CC一定程度上能解决静态相关问题
CASSCF	描述较好	几乎没有考虑
NEVPT2、CASPT2	描述较好	描述较好
MRCC、MRCI	描述较好	描述较好

(http://bbs.keinsci.com/thread-1177-1-1.html)