同样，由于 $N_o$ 极小，对角化步骤变得极为迅速。

通过将 CCSDt 基准态与上述活性 EOM 算符结合，论文定义了：

• DEA-EOMCCSDt(4p-2h){$N_u$}
• DIP-EOMCCSDt(4h-2p){$N_o$}

如果基准态采用完全的 CCSDT，而 EOM 采用活性限制，则分别称为 DEA-EOMCCSDT(4p-2h){$N_u$} 和 DIP-EOMCCSDT(4h-2p){$N_o$}。

2. 关键 Benchmark 体系与计算数据深度解析

论文对三个极具代表性的化学体系进行了详尽的高精度计算，以验证新方法的可靠性。以下是对这些 Benchmark 数据的深度剖析。

2.1 体系一：亚甲基（$\text{CH}_2$）的低躺激发态与单重态-三重态分裂

亚甲基 $\text{CH}_2$ 是量子化学中测试多参考态方法的经典“试金石”。它的基态是三重态 $X\,^3B_1$，而低躺的三个单重态分别是 $A\,^1A_1$、$B\,^1B_1$ 和 $C\,^1A_1$。其中 $A\,^1A_1$ 具有极其显著的多排态特征。

• 计算细节： 采用 TZ2P 基组。DEA 的计算参考态是闭壳层的 $(\text{CH}_2)^{2+}$ 离子（$N=6$ 电子体系），通过附着两个电子得到 $\text{CH}_2$。DIP 的计算参考态是闭壳层的 $(\text{CH}_2)^{2-}$ 离子（$N=10$ 电子体系），通过移去两个电子得到 $\text{CH}_2$。
• 活性空间选择： $N_u = 2$（对应于基准态的 $3a_1$ 和 $1b_1$ 轨道）。

表 2 给出了各个方法相对于 Full CI（全组态相互作用）的计算误差（单位为 kcal/mol）。我们将其核心数据提炼如下表：

数据深入分析与物理机制：

1. 高阶激发算符（4p-2h/4h-2p）的绝对必要性： 常规的 $3p-1h$ 截断方法（如 DEA-EOMCCSD(3p-1h) 和 DIP-EOMCCSD(3h-1p)）表现极差，其 MUE 分别高达 3.64 kcal/mol 和 4.63 kcal/mol。这主要是因为在描述强多参考态的 $A\,^1A_1$ 态时，体系包含明显的双激发贡献，这必须由 4p-2h（或 4h-2p）激发算符才能在数学上描述。一旦将截断级别提升到 4p-2h/4h-2p，DEA 和 DIP 误差立即降到了 0.5 kcal/mol 左右， NPE 也极大地减小。
2. 活性空间方法的完美复现： 在 DEA 计算中，活性空间版本 DEA-EOMCCSD(4p-2h){2}（MUE = 0.54 kcal/mol）和它的全空间父代方法 DEA-EOMCCSD(4p-2h)（MUE = 0.38 kcal/mol）精度几乎完全一致。而在 DIP 计算中，DIP-EOMCCSDt(4h-2p){2} 的 MUE 仅为 0.41 kcal/mol，其表现甚至在微弱程度上优于全空间 CCSDT(4h-2p) 父代方法（MUE = 0.58 kcal/mol），但其所需的计算时间只有父代方法的几个百分点。这证明了仅仅通过引入 2 个活性轨道限制，我们就可以锁住最重要的强相关作用。
3. 为什么三体集群 $T_3$ 在这里影响不大？ 从表中可以看出，加入完整的 $T_3$（CCSDT 级别）或活性三体（CCSDt 级别）对亚甲基单重态-三重态分裂能的改变只有大约 $0.1 - 0.2 \text{ kcal/mol}$。其背后的物理原因是：在这类等化学计量反应（同一分子的内部激发态能垒）中，我们比较的是相同电子数、相同电子排布区域的态。** Neutral 基底态中由 $T_3$ 带来的动力学相关能修正，在激发态与基态之间发生了极大的抵消（Error cancellation）**。然而，这一规律在涉及不同电子数的电离能计算中将会彻底失效。

2.2 体系二：三亚甲基甲烷（TMM）的单重态-三重态分裂能

三亚甲基甲烷（Trimethylenemethane, TMM）是另一个著名的非凯库勒型（Non-Kekulé）双自由基分子。它的基态是具有 $D_{3h}$ 对称性的三重态 $X\,^3A_2'$。而其最低的单重态会通过姜-泰勒效应（Jahn-Teller effect）发生几何畸变，降低对称性至 $C_{2v}$，表现为扭曲的开壳层单重态 $B\,^1A_1$。

• 计算细节： 采用 cc-pVDZ 基组。DEA 的计算参考态是闭壳层的 $(\text{TMM})^{2+}$ 二价阳离子，DIP 的参考态是闭壳层的 $(\text{TMM})^{2-}$ 二价阴离子。
• 活性空间选择： $N_o = N_u = 3$（对应于 TMM 分子的 3 个活性 $\pi$ 轨道：双重简并的 $1e''$ 和非简并的 $2a_2''$ 轨道）。

表 3 汇总了不同计算方法给出的绝热单重态-三重态能量差 $\Delta E_{\text{S-T}} = E(B\,^1A_1) - E(X\,^3A_2')$（单位：kcal/mol）：

数据深入分析：

1. 常规的 $3p-1h$ 和 $3h-1p$ 近似分别给出了 23.92 kcal/mol 和 21.24 kcal/mol 的分裂能，相比经零点能修正后的实验纯电子值（18.1 kcal/mol）存在 3 - 5 kcal/mol 的严重高估。这表明单激发和双激发截断完全无法平衡基态与具有强静态相关特征的开壳层单重态之间的能量关系。
2. 一旦在 DEA 或 DIP 中引入活性空间限制下的 4p-2h 或 4h-2p 激发，$\Delta E_{\text{S-T}}$ 瞬间收敛至 19.0 - 19.3 kcal/mol，与全空间父代方法高度一致，且与实验估计值（18.1 kcal/mol）的误差控制在了 1.0 kcal/mol（化学精度，Chemical Accuracy） 之内。
3. 与亚甲基类似，在 TMM 中，由 $T_3$ 带来的高阶三体相关效应对 $\Delta E_{\text{S-T}}$ 的影响微乎其微（小于 $0.15 \text{ kcal/mol}$），再次证明在等电子数体系间的分裂能计算中，活性 $t_3$ 算符在对角化能量层面发生了完美抵消。这也向计算工作者揭示了一个极具价值的工程经验：若目标仅仅是计算同分异构体分裂能或激发能，无需开启昂贵的三体集群基准态计算，采用经济的 active-space EOMCC(4p-2h/4h-2p) 即可保证精度。

2.3 体系三：23 种原子及分子的垂直双电离势（DIP）统计

当研究体系的电子数发生改变（例如电离过程：$N \rightarrow N-2$）时，物理图像将发生翻天覆地的变化。为了验证方法对动力学相关能不平衡的修复能力，作者对 23 种原子及小分子（包含 $\text{F}_2$, $\text{CO}$, $\text{N}_2$, $\text{H}_2\text{O}$ 等）进行了垂直双电离能计算，并以极高质量的 CIPSI（基于微扰选取的迭代组态相互作用）外推至 Full CI 极限值 作为金标准。

• 计算细节： 采用极其庞大的 aug-cc-pVTZ 弥散基组。
• 活性空间选择： 所有相关占据轨道均设为活性，即 $N_o = n_o$。

表 4 总结了 23 个体系在单重态和三重态双电离能计算中的统计误差（单位：eV）：

统计数据深度剖析与理论洞察：

1. 为什么在 DIP 中基准态 $T_3$ 变得极其关键？ 与前面两个体系不同，在计算双电离能时，DIP-EOMCCSD(3h-1p) 和 DIP-EOMCCSD(4h-2p) 的平均绝对误差分别为 0.59 eV 和 0.23 eV。即使引入了 4h-2p 激发，如果不提升基准态的 CC 级别，依然存在约 0.23 eV 的系统性低估（MSE = -0.23 eV）。这清晰地表明：在发生双电离时，由于 Neutral 基准态和双电离的 Dication 态之间存在巨大的动力学相关能不平衡，基准态中的三体集群（$T_3$）无法被抵消。要想消除这 0.23 eV 的系统性偏差，必须在 $N$ 电子中引入 $T_3$ 相关。
2. DIP-EOMCCSDt(4h-2p){$N_o$} 的惊艳表现： 当我们采用活性三体基准态方法 DIP-EOMCCSDt(4h-2p){$N_o$} 时，单重态的平均绝对误差（MUE）从 0.23 eV 直接骤降至 0.10 eV，最大误差（MaxE）从 0.67 eV 暴跌至 0.18 eV。这几乎完美地比肩了在数学和计算上极其昂贵的全空间 DIP-EOMCCSDT(4h-2p)（MUE = 0.05 eV, MaxE = 0.23 eV）。
3. 对传统微扰近似 [CCSD(T)(a)] 的超越： 传统的微扰三体方法 DIP-EOMCCSD(T)(a)(4h-2p) 虽有较好的平均精度（MUE = 0.07 eV），但其最大误差（MaxE）依然高达 0.43 eV。在面临具有强多参考态特征的体系时，由于微扰理论的分母极易趋近于零，微扰方法会发生严重的数值发散。而本文提出的活性空间非微扰变体 DIP-EOMCCSDt(4h-2p){$N_o$} 成功锁定了最大误差在 0.17 - 0.18 eV，表现出了极其卓越的数值鲁棒性（Robustness）。

3. 代码实现细节、数学缩合与复现指南

本研究所开发的高阶及活性空间 DEA/DIP-EOMCC 系列算法已被高效整合在著名开源量子化学计算包 GAMESS 中。以下是其底层实现的核心技术和用户复现指南。

3.1 核心技术：张量中间体因子化（Factorization）与 Table 1 剖析

在实现 $4p-2h$ 和 $4h-2p$ 投影方程时，如果直接根据公式（如方程 11-13）进行多重循环缩合，会导致极其恐怖的计算开销。为了实现高效的算法，作者引入了一系列精心构造的低维张量中间体（Intermediates）。

以论文中的 Table 1 为例，我们剖析其中两个最关键的中间体构建过程：

1. 中间体 $I_{mn}(\mu)$ 的构建：

这里的 $\bar{h}_{mn}^{ef}$ 是相似变换哈密顿量的双体分量，而 $r_{ef}(\mu)$ 是双电子结合算符的 $2p$ 振幅。通过提前缩合这两个指标，生成一个仅包含两个占据轨道指标的矩阵 $I_{mn}(\mu)$。在随后的三体收缩步骤中（如方程 12 中计算 $\frac{1}{12} I_{mn}(\mu) t_{abc}^{mnk}$），计算复杂度直接从 $\mathcal{O}(n_o^3 n_u^5)$ 降为了极轻量级的二次缩合，从而避免了高阶算符的直接相乘。

2. 中间体 $I^e_{abc}(\mu)$ 的构建：

该中间体包含了最耗时的 $4p-2h$ 振幅 $r_{abcf}^{mn}(\mu)$ 与基准态哈密顿量双体分量的收缩。在活性空间算法中，由于 $r_{abcf}^{mn}(\mu)$ 只有在某些指标属于活性空间时才不为零，程序在构建 $I^e_{abc}(\mu)$ 时，内部循环会自动跳过所有非活性指标。这使得该步骤的计算开销降低了数个数量级。

3.2 软件与开源链接

所有的公式推导通过 Piecuch 课题组自主开发的自动化公式推导与代码生成软件（Automated Formula Derivation and Implementation Software, 类似于 ALCHEMY 或 TCE）导出，从而保证了数十万行 Fortran/C 缩合代码的零错误率。

• 集成平台： GAMESS (General Atomic and Molecular Electronic Structure System) 软件包。
• 官方主页及源码下载： GAMESS Official Website
• 开发组链接： Piecuch Group at MSU

3.3 详尽复现指南与 GAMESS 输入文件编写

若想复现论文中亚甲基 $\text{CH}_2$ 的 DEA-EOMCCSDt(4p-2h){2} 计算，用户需要按照以下步骤准备输入文件和执行计算。

第一步：Neutral 基准态 RHF 计算与活性轨道空间分析

对于 $\text{CH}_2$ 的 DEA 计算，我们首先需要对 $(\text{CH}_2)^{2+}$ 阳离子进行 RHF 计算。计算完成后，在输出文件中寻找分子轨道（MO）的对称性和能量分布。我们需要确定：

1. 最高占据分子轨道（HOMO）和最低未占分子轨道（LUMO）。
2. 确定最活性的两个虚轨道：由于我们需要复现论文中的 $N_u = 2$ 活性空间，这两个活性虚轨道分别对应于基准态重构后的 $3a_1$ 和 $1b_1$ 轨道。

第二步：编写 GAMESS 计算输入文件（ch2_dea.inp）

在 GAMESS 中激活高阶及活性空间 EOMCC 计算，需要在 $CONTRL 和特定耦合集群卡中设置。以下是一个标准的复现输入模板：

 $CONTRL SCFTYP=RHF MULT=1 CCTYP=EOM-CC RUNTYP=ENERGY $END
 $SYSTEM MWORDS=500 $END
 $BASIS  GBASIS=TZ2P $END
 $GUESS  GUESS=HUCKEL $END

 $CCINP  
   ! 激活基准态的活性三体 CCSDt 计算
   CCTYP=CCSDt  
   
   ! 激活双电子结合方程运动耦合集群计算
   EOMTYP=DEA-EOM  
   
   ! 设定活性占据轨道和活性虚轨道数量
   NACTO=0  ! 对于 DEA-EOMCCSDt(4p-2h){2}，活性占据数设为0
   NACTU=2  ! 活性虚轨道数设为 2
   
   ! 设定具体的活性轨道指标（根据实际 RHF 轨道序号填写）
   ! 假设在 RHF 输出中，3a1 和 1b1 的轨道序号分别为 4 和 5
   MINACT=4
   MAXACT=5
   
   ! 计算的目标特征态数量
   NSTATE(1)=1,0,1,0  ! 分别计算不同对称性下的目标态个数
 $END

 $DATA
Methylene (CH2) 2+ Reference for DEA Calculation
C1
C     6.0     0.000000000     0.000000000    -0.122143000
H     1.0     0.000000000     0.865411000     0.723783000
H     1.0     0.000000000    -0.865411000     0.723783000
 $END

第三步：运行与输出结果提取

使用本地编译好的包含高阶 EOM-CC 模块的 GAMESS 可执行文件运行：

rungms ch2_dea.inp 00 1 > ch2_dea.log

运行结束后，在 ch2_dea.log 文件中搜索以下关键字提取关键数据：

• RESULTS OF CCSDt CALCULATION：获取基准态 $(\text{CH}_2)^{2+}$ 的总能量。
• DEA-EOMCCSDt(4p-2h) EXCITATION ENERGIES：获取各个目标特征态的垂直附着能（$\omega_{\mu}^{(+2)}$）。
• 通过将附着能进行差值，即可精确复现表 2 中扭曲开壳层单重态与三重态的分裂能。

4. 关键引用文献与局限性评论

4.1 关键引用文献

本研究建立在半个世纪以来耦合集群理论的演进之上，以下 5 篇文献构成了本项工作的基石：

1. Cížek, J. J. Chem. Phys. 45 (1966) 4256-4266.
   贡献： 首次将多体微扰及耦合集群（CC）理论引入化学领域，定义了 CC 理论的底层数学框架。 (本文引用 [1])
2. Stanton, J. F.; Bartlett, R. J. J. Chem. Phys. 98 (1993) 7029-7039.
   贡献： 形式化地建立了方程运动耦合集群（EOM-CC）理论，为激发态计算提供了非厄米对角化方案。 (本文引用 [3])
3. Piecuch, P. et al. J. Chem. Phys. 110 (1999) 6103-6122.
   贡献： 提出了活性空间耦合集群理论（CCSDt 和 CCSDtq），为在单参考框架下解决强多参考态问题开辟了高效的活性选择道路。 (本文引用 [35])
4. Shen, J.; Piecuch, P. J. Chem. Phys. 138 (2013) 194102.
   贡献： 首次设计并实现了 DEA-EOMCCSD(4p-2h) 和 DIP-EOMCCSD(4h-2p) 方法，论证了高阶激发在双电子过程中的关键作用。 (本文引用 [11])
5. Gururangan, K. et al. J. Chem. Phys. 162 (2025) 061101.
   贡献： 首次独立实现了包含完整三体集群的 full DIP-EOMCCSDT(4h-2p) 理论。 (本文引用 [16])

4.2 本工作局限性客观评论（计算化学家视角）

尽管 DEA/DIP-EOMCCSDt(4p-2h/4h-2p) 方法在精度和计算尺度的平衡上取得了令人瞩目的突破，但在实际的研究应用中，该方法仍存在一些亟待解决的局限性：

1. 活性轨道选择的“艺术性”与不确定性

活性空间方法（如 CCSDt 和 DEA-EOMCCSDt）的一个天然弱点在于，计算结果高度依赖于活性占据轨道 $N_o$ 和活性虚轨道 $N_u$ 的手动选择。在简单的双自由基（如 $\text{CH}_2$、TMM）中，依靠轨道对称性和化学直觉（如选择活性 $\pi$ 轨道）可以很容易挑出合理的活性空间。但在面对复杂过渡金属催化中心、多中心自由基或非对称大分子时，活性轨道的边界往往非常模糊。选择过小的活性空间会导致动力学相关能描述不全，而选择过大则会使计算量迅速失控。目前，该程序尚未集成自动化活性轨道选择算法（如基于自然轨道占有数或引力分析的自动选择）。

2. 双电离势计算中的基准态“双阴离子不稳定性”（Dianionic Instability）

在进行 DIP-EOMCC 计算时，为了利用对称性，其 underlying 闭壳层基准态往往是带负电荷的 $(N)$ 电子体系。例如，在计算中性分子 $M$ 的 DIP 时，通常需要对双阴离子 $M^{2-}$ 进行 RHF 和基准态 CC 计算。然而，在现实中，许多小分子或原子的双阴离子在气相中是极度不稳定的，它们倾向于自动自发脱附电子（Autoionization）。如果在计算中采用了包含大量弥散函数的基组（如 aug-cc-pVTZ），不稳定的基准态波函数极易发散，或者收敛到一个离域的、物理上不正确的连续态波函数，这会直接摧毁后续 EOM 计算的数值基础。虽然作者在文中指出，通过显式引入高阶 4h-2p 激发可以极大缓解这一问题，但在强弥散基组下的收敛性依然是一个潜在的雷区。

3. $N^6$ 计算尺度在大体系中的扩展瓶颈

虽然 active-space 近似成功将 full CCSDT(4h-2p) 的 $N^8$ 降为了 $N^6$，但这依然属于中等偏高的计算尺度。当分子体系的原子数超过 30 个、虚轨道数达到数百个时，$\mathcal{O}(N^6)$ 尺度的对角化步骤（尽管前因子很小）依然会带来沉重的 CPU 和内存带宽负担。该方法目前尚未与局部相关方法（Local Correlation Methods, 如 DLPNO 或 Domain-based Local Pair Natural Orbitals）结合，这限制了它向真正纳米级大分子体系的跨越。

5. 其他补充：应用场景、MRPT对比与张量收缩哲学

为了让读者能更全面、透彻地理解这项工作在量子化学版图中的定位，本节对相关背景进行技术拓展。

5.1 为什么 DIP-EOMCC 是计算俄歇光谱（Auger Spectroscopy）的黄金标准？

在俄歇衰变（Auger Decay）过程中，一个高能光子首先将分子的内层（Core）电子激发或电离，产生一个不稳定的核心空穴（Core-hole）。随后，一个外层价电子跌落入内层空穴，释放的能量将另一个价电子激发至连续态，最终产生一个双电离的二价阳离子（Dication）。整个俄歇谱的峰位和强度直接对应于分子的双电离能谱（DIP Spectrum）。

计算俄歇光谱极其困难，因为：

1. 末态是具有高能激发的双电离态，普通的单激发方法（如 CIS, TDDFT）完全无法描述。
2. 体系包含数以百计的可能末态，传统的 MRCI 方法由于其“非协同性”无法同时平衡描述所有这些高激发态。

DIP-EOMCCSDt(4h-2p) 几乎是为此量身定做的。 由于它的 EOM 算符 $R^{-2}$ 显式包含 4h-2p 激发，它可以极其完美、对称地描述所有可能的外层价空穴之间的排布。此外，由于 EOMCC 的非厄米对角化可以同时输出特征值矩阵的前几十甚至前几百个特征态，这使得在一大基组下一次性计算出高分辨率的俄歇能谱图成为了现实。

5.2 与多参考微扰理论（MRPT2 / CASPT2）的深度对比

许多研究人员在遇到双自由基等问题时，第一反应往往是采用 CASPT2（完全活性空间二级微扰理论）。那么，活性空间 DEA/DIP-EOMCC 相比 CASPT2 的技术优势在哪里？

从技术对比可以看出，DEA/DIP-EOMCC 实际上在理论完整性上显著超越了微扰多参考方法，是真正意义上的“单参考框架处理多参考问题”的终极方案。

5.3 现代耦合集群的张量收缩哲学：为什么“Connected”如此重要？

在公式 (9) 中，我们求解的是：

注意等式左侧下标的 $C$（Connected，已连接）。在多体物理和量子化学中，“已连接”是保证物理量大小协同（Size-extensibility）的唯一途径。根据著名的坎贝尔-贝克-豪斯多夫（Baker-Campbell-Hausdorff, BCH）展开，相似变换哈密顿量 $\bar{H}_N = e^{-T} H e^T$ 展开式可以严格截断在四体项：

这意味着它不包含无穷级数的非物理发散项。通过要求 $\bar{H}_{N, \text{open}}$ 与 EOM 算符 $R$ 的乘积必须全部以连通图（Connected diagrams）的形式收缩，我们在数学上斩断了所有“分离图（Disconnected diagrams）”的产生。这也是为什么不管分子的尺寸如何膨胀（例如，在 100 埃外的空间放一个氦原子），该方法计算出的 $\text{CH}_2$ 单重态-三重态分裂能都可以精确保持恒定不变，而传统的 CI 方法则会因为 disconnected 项的污染而发生灾难性的精度塌缩。

本项研究通过引入活性空间限制（CCSDt 和 active $R_4$），在完全不破坏这种 BCH BCH 严格连通图结构和大小协同性的前提下，将高阶激发从理论云端拉向了实际工程应用，是一项极具智慧的理论物理与计算化学杰作。