企业官网建设流程全解析

在量子化学基准测试里，铁硫簇（iron–sulfur clusters）属于典型的“硬骨头”：金属中心多、未配对电子多、自旋态密集，传统的单参考方法很容易翻车。最近几年，SQD、QSCI 这类基于量子子空间构造的量子-经典混合算法开始在强关联体系上崭露头角，但随之而来的一个新问题被很多人忽略——你以为算出来的是目标自旋态，实际上可能已经混入了大量自旋污染。本文会从铁硫簇的量子化学难点出发，讲清楚 SQD/QSCI 基准测试里为什么要专门做 spin audit（自旋审计），并且用一个可运行的最小模型完整演示审计流程。

适合两类读者：一类是做量子化学计算、想尝试量子算法但没头绪的研究生和工程师；另一类是做量子计算算法、需要找真实化学体系做 benchmark 的开发者。读完你会掌握 S² 期望值怎么算、SQD 风格采样子空间对角化怎么做、以及如何用 spin audit 快速定位“看起来能量很低但波函数已经不对”的问题。

1. 背景：铁硫簇为什么需要“自旋审计”

1.1 铁硫簇：量子化学的“压力测试”

铁硫簇是生物体内一类非常古老且重要的金属辅因子，常见的结构包括 [2Fe-2S]、[4Fe-4S] 等。它们在光合作用电子传递链、固氮酶、线粒体呼吸链里都承担着电子转移和催化功能。从化学本质上看，铁硫簇的核心特征有两个：一是铁原子之间存在较强的交换耦合，二是 d 电子数多，活性空间里多个轨道能量接近，电子关联效应非常强。

这意味着什么呢？简单说，铁硫簇的基态往往不是单参考态能描述的。DFT 的 Kohn-Sham 单行列式波函数可能出现严重的自旋污染，CASSCF 如果活性空间选得不够又会漏掉关键的动态关联，CCSD(T) 这类高精度方法在大体系上又贵得离谱。所以铁硫簇成为了量子化学方法发展的“压力测试”：你真的有本事，就把 [4Fe-4S] 的多重态能量差算准。

1.2 从经典计算到量子-经典混合算法

面对铁硫簇这种强关联体系，量子计算提供了一条看起来很自然的路径：把 Hamiltonian 编码到量子比特上，用浅层量子电路制备一个包含多参考特征的试探态，再通过测量得到的关键组态，交给经典计算机做最后的对角化处理。

SQD（Sample-based Quantum Diagonalization，基于采样的量子对角化）和 QSCI 就属于这一类思路。它们的共同点是：

1. 用量子态制备电路生成一个参考态，通常是一个包含多个重要 Slater 行列式的叠加态。
2. 对这个量子态做测量采样，得到一系列 bitstring，也就是占据数向量。
3. 用这些 bitstring 张成一个有效子空间，子空间大小远小于完整的 Fock 空间。
4. 在经典计算机上构造 Hamiltonian 在该子空间中的矩阵并求解本征问题，得到能量和波函数。

这类方法的优势在于，量子电路只需要制备一个“足够好但不完美”的试探态，剩下的纠错和精确化交给经典对角化完成。相比完全依赖量子硬件的变分特征值求解器（VQE），它对噪声的容忍度更高，也更能处理强关联体系。

但这类方法也有一个隐患：采样子空间可能没有精确保持体系的自旋对称性，最终解出的“基态”虽然能量很低，却可能是一个自旋混合态。对于铁硫簇这种自旋态密集、多重态能量差很小的体系，这会直接导致错误的多重态排序和错误的磁耦合常数。

1.3 什么是自旋审计

在量子化学里，自旋算符的平方 S² 的本征值是 S(S+1)，其中 S 是总自旋量子数。比如单重态 S=0，S² 期望值为 0；三重态 S=1，S² 期望值为 2；五重态 S=2，S² 期望值为 6。如果一个波函数的 S² 期望值偏离这些整数本征值，就说明它混合了不同自旋的成分，这就是自旋污染（spin contamination）。

Spin audit，就是我在这类量子化学基准测试里加的一道检查流程：不只看能量，还要检查最终波函数的 S²、Sz 等自旋相关量是否落在目标自旋态的合理范围内。它的目标是回答几个问题：

• 算出来的基态是目标自旋态吗？
• 波函数是否混入了其他自旋成分？
• 如果改变采样数量、电路深度或截断阈值，自旋纯度是否稳定？
• 能量差（如单重态-三重态能隙）是否建立在自旋一致的波函数之上？

如果是传统 CASCI 或 CASSCF，自旋对称性通常由算法本身保证；但在 SQD/QSCI 这类“采样 + 子空间截断”的方法里，自旋对称性并不会自动成立。采样数量不够、电路表达能力不足、噪声扰动，都可能导致有效子空间缺少某些关键的自旋补偿项。所以，把 spin audit 写进基准测试流程，不是可选项，而是必选项。

2. 环境准备与版本说明

2.1 运行环境

本文的演示代码主要依赖 Python 和 NumPy，不需要真实量子硬件，也不需要安装重量级的量子化学软件包。核心思路是把计算基下的 Hamiltonian 矩阵和自旋算符矩阵显式构造出来，然后演示“采样 → 子空间对角化 → 自旋审计”的完整链路。

推荐环境如下：

• 操作系统：Windows 10/11、Ubuntu 20.04 或 macOS 均可。
• Python：3.9 及以上版本。
• 依赖库：numpy、scipy。
• 可选扩展：pyscf、openfermion、qiskit（迁移到真实化学体系时使用）。

版本需要根据你的项目实际情况调整。本文示例以常见环境为例，重点演示算法思路，而不是绑定某个特定版本。

2.2 项目结构

建议用一个干净的目录组织代码和输出：

spin-audit-sqd/ ├── README.md # 项目说明与运行方式 ├── requirements.txt # 依赖清单 ├── main.py # 主脚本：构造模型、采样、自旋审计 └── outputs/ └── audit_report.log # 审计日志

requirements.txt可以写成这样：

numpy>=1.24 scipy>=1.10

安装依赖：

pip install -r requirements.txt

如果你的环境里已经有 PySCF、OpenFermion 等量子化学工具，后续可以替换成真实铁硫簇模型，但核心的 spin audit 函数完全不需要改动。

3. 核心原理拆解

3.1 SQD 的思路：采样构造子空间

SQD 的核心思想可以概括成一句话：不要在完整 Fock 空间里做对角化，而是用测量采样找到少数重要的组态，再在它们张成的子空间里对角化。

具体流程通常包含这几步：

1. 制备一个参数化量子态 |ψ(θ)⟩，这个态可以是某个浅层电路作用在 Hartree-Fock 态上的结果。
2. 测量 |ψ(θ)⟩，得到 M 个 bitstring，比如 0101、1001 等。每一个 bitstring 对应一个占据数向量，也就是一个 Slater 行列式（计算基矢态）。
3. 对 bitstring 做去重和筛选，去掉重复项和能量上明显无关的组态，得到一个子空间基底 {|φ_i⟩}。
4. 构造 Hamiltonian 在该子空间中的矩阵 H_ij = ⟨φ_i|H|φ_j⟩。
5. 经典对角化 H_ij，得到近似基态能量和波函数。

这种方法之所以有效，是因为强关联体系的基态虽然复杂，但通常集中分布在少数重要组态张成的空间里。只要量子态制备电路能让这些重要组态有可观的概率被采样到，子空间对角化就能恢复出高精度结果。

3.2 QSCI 与同族方法的共同点

QSCI 在文献中常与 SQD 并列讨论，它们本质上都是“用量子态采样确定子空间，用经典对角化求能量”的混合策略。不同实现可能在某些细节上有差异，比如如何生成试探态、如何选择 bitstring、如何在含噪声环境下做纠错，但核心框架是一致的。

所以在讨论 spin audit 时，可以把它们当作同一族方法看待。正因为这种“采样 + 截断”的框架，自旋对称性才可能被破坏，需要额外的检查。

3.3 自旋算符与 S² 期望值计算

在二次量子化框架下，自旋算符可以写成产生和湮灭算符的组合。假设分子轨道按 α、β 自旋区分，那么：

• Sz = (Nα - Nβ) / 2，其中 Nα 是 α 电子数算符，Nβ 是 β 电子数算符。
• S⁺ = ∑_i c_{iα}† c_{iβ}，S⁻ = ∑_i c_{iβ}† c_{iα}。
• S² = Sz² + (S⁺S⁻ + S⁻S⁺) / 2。

对于任意一个量子态 |Ψ⟩，只要能得到 S² 矩阵，就可以计算：

S²_expected = ⟨Ψ|S²|Ψ⟩

然后与目标自旋态的理论值比较。以两电子体系为例，如果目标是单重态，S² 应接近 0；如果目标是三重态，S² 应接近 2。如果计算结果明显偏离，就必须谨慎解读能量。

3.4 Spin audit 的检查清单

在实际基准测试中，我建议至少检查以下内容：

• 基态的 S² 期望值是否接近目标自旋态的理论值。
• Sz 是否守恒，即总磁量子数是否符合设定。
• 不同采样数量下 S² 是否稳定。
• 子空间截断阈值改变后，自旋态是否发生跳变。
• 单重态与三重态等不同多重态，是否在同等自旋纯度条件下比较能量。

这几点构成了一个最小可用的 spin audit 协议。下面用代码完整演示。

4. 完整实战案例：SQD 风格采样与自旋审计

4.1 模型选择与轨道约定

为了在一个示例里同时做到“可运行”和“接近铁硫簇物理”，这里使用两轨道 Hubbard 模型。它有两个空间轨道，每个轨道可容纳 α、β 两个自旋电子，一共 4 个自旋轨道，Fock 空间大小为 2⁴ = 16 个基态。模型 Hamiltonian 为：

H = -t ∑_σ (c_{1σ}† c_{2σ} + c_{2σ}† c_{1σ}) + U ∑_i n_{iα} n_{iβ}

其中 t 表示轨道间的跃迁积分，U 表示同一轨道内两个电子占据时的库仑排斥能。当 U 比较大时，体系倾向于形成局域自旋，表现出反铁磁耦合特征，这和铁硫簇中金属中心间的超交换相互作用有类似之处。实际铁硫簇计算只需要把这里的 H 矩阵替换成由 PySCF 等程序生成的 ab initio Hamiltonian，自旋审计代码完全复用。

轨道顺序约定如下：

0 号自旋轨道：空间轨道 1 的 α 电子 1 号自旋轨道：空间轨道 1 的 β 电子 2 号自旋轨道：空间轨道 2 的 α 电子 3 号自旋轨道：空间轨道 2 的 β 电子

4.2 构造费米子算符与 Hamiltonian 矩阵

下面的代码实现了费米子产生、湮灭算符在 Fock 态上的作用，并构造 Hamiltonian 矩阵。注意费米子反对易符号由“目标自旋轨道之前占据数之和”决定。

# 文件路径：main.py import numpy as np # 自旋轨道编号：0=1α, 1=1β, 2=2α, 3=2β def spinorb_idx(orb, spin): """orb 取 0/1，spin 取 0=α/1=β，返回自旋轨道编号""" return 2 * orb + spin def occupied_before(bitstring, idx): """返回 bitstring 中 idx 之前占据数之和，用于确定费米子反对易符号""" mask = (1 << idx) - 1 return bin(bitstring & mask).count("1") def create(bitstring, idx): """产生算符 c_idx† 作用在 Fock 态 bitstring 上""" if (bitstring >> idx) & 1: return None sign = (-1) ** occupied_before(bitstring, idx) return bitstring | (1 << idx), sign def annihilate(bitstring, idx): """湮灭算符 c_idx 作用在 Fock 态 bitstring 上""" if not ((bitstring >> idx) & 1): return None sign = (-1) ** occupied_before(bitstring, idx) return bitstring & ~(1 << idx), sign def apply_operator(ops, bitstring): """按顺序作用一串费米子算符。 ops 中 '+' 表示产生，'-' 表示湮灭，例如 [('-', 2), ('+', 0)] 表示 c_0† c_2 """ current = bitstring sign = 1.0 for op, idx in ops: if op == '+': res = create(current, idx) else: res = annihilate(current, idx) if res is None: return None current, s = res sign *= s return current, sign def build_H(t=0.5, U=3.0): """构造两轨道 Hubbard 模型在 16 维 Fock 空间中的 Hamiltonian 矩阵""" n = 16 H = np.zeros((n, n)) # 跃迁项：-t (c_{1σ}† c_{2σ} + c_{2σ}† c_{1σ}) for spin in range(2): i1 = spinorb_idx(0, spin) i2 = spinorb_idx(1, spin) for nu in range(n): # -t c_{1σ}† c_{2σ} res = apply_operator([('-', i2), ('+', i1)], nu) if res is not None: mu, phase = res H[mu, nu] += -t * phase # -t c_{2σ}† c_{1σ} res = apply_operator([('-', i1), ('+', i2)], nu) if res is not None: mu, phase = res H[mu, nu] += -t * phase # Hubbard U：同一空间轨道上 α、β 同时占据时增加能量 U for orb in range(2): ia = spinorb_idx(orb, 0) ib = spinorb_idx(orb, 1) for nu in range(n): if ((nu >> ia) & 1) and ((nu >> ib) & 1): H[nu, nu] += U return (H + H.T) / 2

这里有一个容易忽略的细节：create和annihilate返回的符号不能省。如果直接在计算基下用 0/1 占据数拼接 Hamiltonian，忽略费米子反对易符号，得到的矩阵在物理上是错的。本文的做法是把产生湮灭算符显式作用到每个 Fock 态上，由occupied_before自动处理符号，适合小体系模型验证。

4.3 构造 S² 矩阵

自旋算符矩阵同样在 16 维 Fock 空间构造。先构造 Sz、S⁺、S⁻，再组合成 S²。

def build_sz(): """构造 Sz 矩阵：Sz = (Nα - Nβ) / 2""" n =