📚 课程目录第4章分子积分计算与活性空间选择

4. 分子积分计算与活性空间选择

要把一个分子的电子结构问题交给量子计算机，第一步不是写量子线路，而是把化学问题“翻译”成量子比特能读懂的哈密顿量。这个翻译过程的核心，就是 Hartree-Fock 轨道、电子积分与活性空间选择。本节围绕“如何从分子结构出发，得到一个小而可算的二次量子化哈密顿量”展开，并用 PySCF 与 OpenFermion 完成一次可复现的小分子实战。

一、从原子坐标到 Hartree-Fock 轨道

量子化学计算的起点是原子核坐标、基组和电子数。以氢分子 H₂ 为例，若取键长 0.74 Å、基组 STO-3G，则每个氢原子提供 1 个 s 型基函数，共 2 个空间轨道。Hartree-Fock（HF）方法在 Born-Oppenheimer 近似下求解 Roothaan 方程：

F C = S C ε

其中 F 是 Fock 矩阵，S 是重叠矩阵，C 是分子轨道系数矩阵，ε 是轨道能。求解后得到一组正交归一的分子轨道，并按能量从低到高排列。对 H₂ 而言，两个电子填入能量最低的成键轨道 σ_g，反键轨道 σ_u 为空。HF 行列式给出了一个平均场近似下的参考态，但它忽略了电子间的瞬时关联，这正是后续需要量子计算或后 HF 方法修正的部分。

需要特别说明的是：HF 轨道不是唯一的。任何对占据轨道之间、虚轨道之间的酉变换都不改变 HF 行列式能量，但会改变积分值。因此在实际流程中，必须固定一套轨道表示，通常选择正则 HF 轨道，因为此时 Fock 矩阵对角，轨道能就是对角元，便于后续做 MP2、CCSD 或量子线路的参考态准备。

二、电子积分：哈密顿量的“原材料”

在二次量子化形式下，分子电子哈密顿量写作：

H = Σ_{pq} h_{pq} a†_p a_q
  + 1/2 Σ_{pqrs} h_{pqrs} a†_p a†_q a_r a_s + E_nuc

其中 h_pq 是单电子积分，包含动能与核吸引能；h_pqrs 是双电子积分，描述电子间库仑排斥；E_nuc 是核间排斥能。对 STO-3G 的 H₂，空间轨道只有 2 个，自旋轨道 4 个，双电子积分数量很少，可以手工列出。但对更大的分子，积分数量随基函数数目 N 以 N⁴ 增长，必须依赖程序自动计算。

在 PySCF 中，这些积分可以通过 mol.intor 或 scf 对象获得。例如：

from pyscf import gto, scf
mol = gto.M(atom='H 0 0 0; H 0 0 0.74', basis='sto-3g')
mf = scf.RHF(mol).run()
h1 = mf.mo_coeff.T @ mf.get_hcore() @ mf.mo_coeff
eri = mol.intor('int2e', aosym='s1')
h2 = np.einsum('pi,pqrs,qj,rk,sl->ijkl', mf.mo_coeff, eri,
               mf.mo_coeff, mf.mo_coeff, mf.mo_coeff)

这里 h1 是 MO 基下的单电子积分，h2 是 MO 基下的双电子积分。注意索引顺序与化学家记号、物理学家记号之间的差异，OpenFermion 与 PySCF 在自旋轨道排序上也有约定差别，工程实现时必须核对。

三、活性空间与冻结核心近似

全空间计算对量子硬件仍然过于昂贵。化学上有一个重要观察：内层核心电子几乎不参与成键与反应，强关联主要发生在费米能级附近的少数轨道上。于是引入冻结核心近似：把能量最低的核心轨道始终双占据，不纳入量子线路；同时把能量很高的虚轨道也冻结掉。剩下的占据轨道与部分虚轨道构成活性空间，记作 CAS(n, m)，即 n 个电子在 m 个空间轨道中做完全组态相互作用。

活性空间选择没有唯一答案，常用判据包括：

轨道能级：选取 HOMO/LUMO 附近若干轨道；
自然轨道占据数：偏离 0 或 2 较多的轨道应纳入活性空间；
化学直觉：断键、成键、孤对电子、d 轨道等关键轨道必须包含；
资源约束：活性空间每增加一个空间轨道，量子比特数增加 2，线路深度显著上升。

对 H₂/STO-3G，最小活性空间是 CAS(2,2)，即 2 个电子在 σ_g 与 σ_u 两个空间轨道中，对应 4 个自旋轨道、4 个量子比特。若冻结核心，则 H₂ 没有核心可冻；对 LiH 等分子，1s 核心可冻结，活性空间可取 CAS(2,3) 或 CAS(2,5)。

四、用 PySCF + OpenFermion 生成哈密顿量

下面给出一个完整可执行的最小流程。目标：对 H₂/STO-3G 生成 Jordan-Wigner 变换后的量子比特哈密顿量，并打印项数与常数项。

import numpy as np
from pyscf import gto, scf, ao2mo
from openfermion import InteractionOperator, get_fermion_operator
from openfermion.transforms import jordan_wigner
from openfermion.ops import QubitOperator

# 1. 构建分子并做 RHF
mol = gto.M(atom='H 0 0 0; H 0 0 0.74', basis='sto-3g', verbose=0)
mf = scf.RHF(mol).run()

# 2. 取 MO 基积分
n_orb = mol.nao_nr()
h1 = mf.mo_coeff.T @ mf.get_hcore() @ mf.mo_coeff
eri = ao2mo.restore(1, ao2mo.kernel(mol, mf.mo_coeff), n_orb)
h2 = np.einsum('pqrs->prqs', eri)  # 调整为 OpenFermion 约定

# 3. 构造 InteractionOperator（自旋轨道展开）
from openfermion.chem import MolecularData
from openfermion.transforms import get_fermion_operator
from openfermion.ops import InteractionOperator

# 这里用 OpenFermion 的 MolecularData 更直接
geometry = [('H', (0., 0., 0.)), ('H', (0., 0., 0.74))]
molecule = MolecularData(geometry, 'sto-3g', 1, 0)
molecule = run_pyscf(molecule, run_scf=True, run_mp2=False, run_cisd=False)
hamiltonian = molecule.get_molecular_hamiltonian()

# 4. Jordan-Wigner 变换
qubit_ham = jordan_wigner(hamiltonian)
print('项数:', len(qubit_ham.terms))
print('常数项:', qubit_ham.terms.get((), 0.0))

运行后，H₂/STO-3G 的量子比特哈密顿量通常包含约 15 项（含常数项），其中常数项对应核排斥能与 HF 能量中的部分贡献。若把活性空间扩大到 CAS(2,2) 以外的轨道，项数会迅速增加。工程上常用 get_fermion_operator 先得到费米子算符，再用 jordan_wigner 或 bravyi_kitaev 做映射，最后用 QubitOperator 的 compress 合并同类项。

五、结果判读与常见陷阱

生成哈密顿量后，应检查三件事：第一，常数项是否等于核排斥能加上冻结核心的贡献；第二，单电子项与双电子项的系数是否与手工估算量级一致；第三，量子比特数是否等于活性空间自旋轨道数。若发现项数异常多，通常是未做对称性约简或未冻结核心；若系数符号与预期相反，多半是索引顺序或自旋轨道排序约定不一致。

注意：PySCF 默认使用化学家记号 (pq|rs)，OpenFermion 使用物理学家记号，两者索引顺序不同。直接搬运积分而不做 einsum 重排，会导致哈密顿量错误，进而使后续 VQE 能量偏离真实值。

实用建议：先用 H₂/STO-3G 跑通全流程，确认 HF 能量、哈密顿量项数、VQE 基态能量与 FCI 一致，再迁移到 LiH、H₂O 等分子。活性空间选择可先用 CAS(n, m) 小窗口验证，再逐步扩大。

本节的核心链条是：分子结构 → HF 轨道 → 电子积分 → 活性空间与冻结核心 → 二次量子化哈密顿量 → 量子比特映射。掌握这条链条后，后续的 VQE、QPE 或量子线路设计才有可靠的输入。实际工程中，活性空间选择往往决定整个量子计算任务的可行性与精度，需要在化学直觉与硬件资源之间反复权衡。