来源论文: https://arxiv.org/abs/2607.03395v2 生成时间: Jul 08, 2026 10:45

0. 执行摘要

吡嗪（Pyrazine, $C_4H_4N_2$）分子的超快光弛豫过程是分子物理、光化学与理论化学界公认的非绝热动力学基准（Benchmark）体系。其激发态势能面上密集的锥形交叉（Conical Intersections, CIs）使得光激发后的非辐射衰减在几十飞秒内发生，这一现象的精确模拟对电子结构理论和量子动力学方法提出了极高要求的双重挑战。

本研究利用新开发的多态相似性约束耦合簇（Multistate Similarity Constrained Coupled Cluster, Multistate SCCSD）理论，克服了标准耦合簇（CCSD）在锥形交叉等准简正区域出现非物理复数解与势能面断裂的宿疾；同时，结合从头算多态产卵方法（Ab Initio Multiple Spawning, AIMS），在包含全部24个内部振动自由度（全维势能面）的基础上，实现了首次完全基于高精度耦合簇理论的吡嗪光弛豫全飞秒动力学模拟。

通过对模拟轨迹进行直接光谱学可观测量的重构，本工作计算了：

1. 时间分辨光电子能谱（Time-Resolved Photoelectron Spectroscopy, TR-PES）
2. 氮和碳 K 边的时间分辨 X 射线吸收光谱（Time-Resolved X-ray Absorption Spectroscopy, TR-XAS）

结果不仅在定量上与实验测量结果高度契合（如预测的亮态 $^1B_{2u}$ 衰减寿命为 11 fs，实验值为 $13 \pm 3$ fs），而且从根本上澄清了吡嗪是否存在“激发态显著复苏”以及 $^1B_{3u}$ 与 $^1A_u$ 两状态之间量子相干相消（Beats）的微观图像。这一工作标志着高精度电子相关效应理论与全维非绝热量子核动力学的大尺度融合取得了实质性突破。

1. 核心科学问题、理论基础与方法细节

1.1 吡嗪超快光弛豫的科学争议与多态机制

吡嗪在接受紫外光照射后，会迅速跃迁到其亮态 $^1B_{2u}(\pi\pi^*)$。由于该状态与邻近的暗态 $^1B_{3u}(n\pi^*)$ 以及 $^1A_u(n\pi^*)$ 存在强烈的非绝热和振动耦合（Vibronic Coupling），分子会在极其短暂的时间内发生非辐射内转换（Internal Conversion, IC）。

长期以来，围绕吡嗪的衰减机制存在两个主要学术争论：

1. 两态机制 vs. 三态机制：早期受限于计算能力和低维简化模型，许多研究倾向于简单的 $^1B_{2u} \rightarrow ^1B_{3u}$ 两态内转换图像。然而，近年来高精度势能面计算表明，介于亮态和暗态之间的暗态 $^1A_u$ 在弛豫路径中扮演了不可忽视的角色，即必须引入三态机制。
2. 亮态 $^1B_{2u}$ 的复苏（Revival）与量子相干（Beats）：许多理论动力学模拟预测，亮态 $^1B_{2u}$ 在衰减后会在 80 fs 左右出现显著的种群复苏。然而，最先进的超快光电子成像和 TR-PES 实验仅观测到了极微弱甚至无法确证的复苏信号。此外，$^1B_{3u}$ 与 $^1A_u$ 之间是否存在周期为 30-60 fs 的相干相消振荡，其在实验光谱上的指征也一直模糊不清。

要想解决这些争议，必须在动力学计算中满足两个条件：一是在全维（24个简正模式）空间中演化核波包以防降维伪影，二是利用高精度多态电子结构方法实时计算受电子相关效应强烈影响的非绝热势能面与非绝热耦合矢量。

1.2 高级电子结构理论：耦合簇理论 (CCSD) 及其非绝热局限性

对于处理弱相关或单参考主导的基态与激发态，耦合簇单双激发方法（CCSD）或方程式耦合簇方法（EOM-CCSD）由于其大小一致性（Size-Extensivity）和高精度的相关能恢复能力，是极佳的选择。然而，当直接将标准耦合簇方法应用于分子动力学，特别是非绝热动力学时，会遭遇致命的方法学失效（Methodological Failure）：

1. 非厄米哈密顿量非物理行为：耦合簇理论基于相似变换（Similarity Transformation）方法： $$\bar{H} = e^{-T} H e^T$$ 由于变换算符 $e^{-T}$ 与 $e^T$ 不互为伴随，相似变换后的有效哈密顿量 $\bar{H}$ 是**非厄米（Non-Hermitian）**的。在大多数平衡态附近，这不会产生严重后果，但在势能交叉区域，非厄米哈密顿量会出现复共轭特征值对（Complex Conjugate Eigenvalue Pairs），导致激发态能量变为复数，势能面在此断裂。
2. 势能交叉点（Conical Intersection）描述错误：在非厄米框架下，两个状态交叉处的数学拓扑不再是单点式的“锥形交叉（Conical Intersection）”，而是一个非物理的“异常点（Exceptional Point, EP）”，在该点上不仅特征值简并，连特征矢量也会发生合并。这导致在 EP 附近，激发态受扰动极度敏感，核梯度和非绝热耦合矢量（Non-Adiabatic Coupling Vectors, NACVs）发生数学发散，使非绝热动力学轨迹无法连续积分。

1.3 相似性约束耦合簇理论 (SCCSD) 及其多态扩展

为了克服标准耦合簇理论在锥形交叉附近的失效，Kjønstad 和 Koch 等人此前开发了相似性约束耦合簇方法（SCCSD）。其核心思想是通过在耦合簇方程中显式引入拉格朗日乘子，强制约束相互交叉或准简正的两个状态（例如状态 $i$ 和 $j$）相互正交，从而保证即便在非厄米算符下也具有实数能量：

通过强制正交化约束，相似变换哈密顿量的特征值在交叉区域退化为真实的单点锥形交叉，避免了复数根和异常点的出现，NACV 表现得平滑且物理意义明确。然而，传统的 SCCSD 只能描述两个状态的简正。在吡嗪这种涉及 $^1B_{2u}$、$^1B_{3u}$ 和 $^1A_u$ 三个电子态超快内转换的过程中，传统的两态约束极易因为第三个状态的逼近而发生失效。

本工作发展了自适应多态相似性约束耦合簇（Multistate SCCSD）方法。其自适应控制方案如下： 在每个核动力学步，首先在标准 CCSD 级别下计算所有感兴趣激发态的能量。计算所有激发态两两之间的能量差：

定义一个能量阈值 $\tau$（本工作设为 0.01 eV）。若对于某些状态对有 $\Delta E_{ij} < \tau$，说明系统正在穿过锥形交叉区，此时系统自适应地将计算切换为约束这些状态相互正交的 SCCSD 水平；而在能量间距大于 $\tau$ 的宽阔区域，计算则平滑回归到高效率的标准 CCSD。通过这种方案，多态动力学模拟能够无缝穿越多个势能面相交区域，且整个势能面上能量与梯度的连续性得到了完美的维持。

拉格朗日函数可表示为：

其中 $\mathcal{C}$ 代表被约束的准简正状态对集合，$\lambda_{ij}$ 是对应的拉格朗日乘子。通过求解平稳点方程，可以得到自适应约束下的激发态能量、受约束梯度（受力矢量）和相应的非绝热耦合矢量。

1.4 非绝热动力学模拟：从头算多态产卵 (AIMS) 理论

动力学模拟采用了基于全量子核动力学框架的从头算多态产卵（Ab Initio Multiple Spawning, AIMS）方法。AIMS 将核波包展开为一组运动的、多维冷冻高斯基函数（Trajectory Basis Functions, TBFs）的线性组合：

其中 $I$ 表示电子态指数，$\chi_{Ij}$ 是在第 $I$ 个电子态上运行的第 $j$ 个高斯核波包。$\mathbf{\xi}_{Ij}(t)$ 包含了高斯波包的中心位置、动量及相位。$C_{Ij}(t)$ 是波包展开系数，其含时演化由 Schrödnger 方程控制：

其中 $\mathbf{S}$ 和 $\mathbf{H}$ 分别为 TBFs 之间的重叠矩阵和有效哈密顿量矩阵，哈密顿量矩阵元包含了非绝热耦合项：

这里 $\mathbf{d}_{IK} = \langle \Phi_I | \nabla_R | \Phi_K \rangle$ 即非绝热导数耦合（Derivative Coupling）矢量。当某一个 TBF 靠近两个激发态的锥形交叉区，且它们之间的非绝热耦合强度大于预设阈值时，AIMS 算法会执行“产卵（Spawning）”操作，即在目标激发态上创生（Spawn）一个新的 TBF。随着动力学演化，新旧波包发生干涉、分裂，从而在全维空间下物理、连续地重现非绝热内转换过程。

1.5 动力学-电子结构自适应接口与非绝热耦合评价

为了将 AIMS 与多态 CCSD/SCCSD 结合，本研究开发了高效的自适应数据传输接口。该接口的核心在于如何精准且廉价地在标准 CCSD 与多态 SCCSD 之间进行切换，并计算非绝热耦合。接口的执行序列逻辑如下：

                    [ 动力学程序 (FMS90) 发送几何构型 R ]
                                     │
                       [ 计算多态激发态 CCSD 能量 E_i ]
                                     │
                     ┌───────────────┴───────────────┐
                     │                               │
             [ 所有 ΔE_ij > τ ]             [ 存在某个 ΔE_ij < τ ]
                     │                               │
         [ 采用标准 CCSD 计算 ]             [ 采用多态 SCCSD 约束 ]
         [ 计算无约束梯度和耦合 ]             [ 计算约束梯度和正则化耦合 ]
                     │                               │
                     └───────────────┬───────────────┘
                                     │
                      [ 返回能量、梯度与耦合至 FMS90 ]
                                     │
                     [ FMS90 进行核波包传播与产卵决策 ]

接口在计算受约束状态与非受约束状态之间的非绝热耦合（导数耦合矢量 $\mathbf{d}_{kl}$）时，采用了双正交（Biorthonormal）耦合簇形式。对 Lagrangian 函数进行变分求导所得的耦合算符形式为：

为了防止多态动力学中导数耦合不满足总和规则（Sum Rule）从而引入非物理不连续，算法直接忽略了二次响应中违反总和规则的微扰项（Neglected terms with minor contributions to ensure translational invariance），从而极大地增强了积分的数值稳定性。

2. 关键 Benchmark 体系、计算所得数据与性能数据

2.1 静态临界点与势能面扫描

为确保动力学轨线的可靠性，首先在 CCSD/cc-pVDZ 级别下对吡嗪基态（$S_0$）及两个最低激发态（$S_1$ 和 $S_2$）的关键静态几何特征进行了表征。计算所得的临界点几何构型及振动频率列于下表：

表 1：基态（$S_0$）平衡几何构型（CCSD/cc-pVDZ, 长度单位：a.u.）

同时，计算了 $S_0$ 的简正模式振动频率（共30个模式，单位：$\text{cm}^{-1}$），其中最低频模式为 $360\text{ cm}^{-1}$，最高频模式在 $3227\text{ cm}^{-1}$。这些频率直接用于构建 $0\text{ K}$ 下的 Wigner 波包分布函数以抽取动力学初始构型。

针对 $S_1$ 激发态，表征了两个局域极小值点：

1. $^1B_{3u}$ 局域极小值（见图 4B 中的 ‘Min’）：其平衡几何构型对称性为 $D_{2h}$。最低振动频率为 $183\text{ cm}^{-1}$，无虚频，证实其为真正的稳定驻点（Table S5）。
2. $^1A_u$ 过渡态（见图 4B 中的 ‘TS’）：其最低振动频率包含一个虚频 $565i\text{ cm}^{-1}$（Table S8），该虚频模式（$B_{3g}$）强烈地将系统引导至 $^1B_{3u}$ 极小值。两态之间的过渡垒能量仅为 $0.06\text{ eV}$（见图 4C 势能面等高线图）。

通过对简正坐标进行一维能量扫描（Potential Energy Scans, 见图 4F 和图 S27-S31），发现：

• 对称化坐标 $\nu_{8a}$（$A_g$ 模式）：沿着该坐标运动时，$^1B_{3u}$ 和 $^1A_u$ 的能量曲线出现显著的交叉行为。这直接印证了 $\nu_{8a}$ 对系统在不同激发态布居分配中所起的调制作用。
• 非绝热交叉面（MECI）：通过多态 SCCSD 定位了 $S_2(^1B_{2u})/S_3(^1A_u)$ 最小能量锥形交叉点（MECI, Table S10）以及 $S_1(^1A_u)/S_2(^1B_{3u})$ 锥形交叉点（MECI, Table S12），它们的势能面能隙均在 $0.001\text{ eV}$ 以下，证实了算法描述锥交叉的优越性。

2.2 时间分辨光电子能谱 (TR-PES) 模拟与定量对比

利用动力学模拟产生的 445 个轨迹基函数（TBFs），计算了每个几何位置上的 Dyson 轨道跃迁强度，从而重构了时间分辨光电子能谱（TR-PES）。光谱经过了 $0.73\text{ eV}$ 的整体平移以对齐基态实验能阶，并应用了 $13\text{ fs}$ 时间半峰宽（FWHM）和 $0.3\text{ eV}$ 能量半峰宽的双维高斯拓宽。这一模拟重构直接与经典实验（Ref. 30）进行了定量对标，如图 1 所示。

关键能谱窗口的光谱指征提取与物理寿命分析：

1. 亮态 $^1B_{2u}$ 的超快衰减寿命（$6.9 - 7.4\text{ eV}$ 探测窗口）：

   • 理论模拟：该窗口对应的光谱带极其短暂（见图 1B）。对其积分强度时间曲线（图 1E 中的实线）进行单指数拟合，提取出的衰减特征寿命为： $$\tau_{theory} = 11\text{ fs}$$ 如果采用更宽的能谱积分范围 $6.9 - 7.9\text{ eV}$（包纳整条 $^1B_{2u}$ 带，图 S2D），衰减寿命增加至 $24\text{ fs}$。
   • 实验对比：实验测得的衰减时间在这一窗口中为： $$\tau_{exp} = 13 \pm 3\text{ fs} \quad (\text{Ref. 30})$$ 而在其他实验重构中分别报告为 $22 \pm 3\text{ fs}$ (Ref. 29) 和 $23 \pm 4\text{ fs}$ (Ref. 31)。高精度理论模拟与实验结果在误差范围内取得了高度的定量契合。

2. 暗态布居的相干振荡（$5.0 - 5.2\text{ eV}$ 与 $6.3 - 6.5\text{ eV}$ 探测窗口）：

   • 在低于 $6.8\text{ eV}$ 的结合能区域，能谱由一个宽阔、长寿命的特征带主导。这是由跃迁到 $^1B_{3u}$ 和 $^1A_u$ 状态产生的 Dyson 信号高度重叠形成的（图 1D）。
   • 通过精确积分两能段（图 1D），可以发现，在 $6.3 - 6.5\text{ eV}$ 的积分信号（红色实线）与 $^1B_{3u}$ 的动力学极大值时间点高度同步（垂直红色虚线，约 $75\text{ fs}$ 与 $135\text{ fs}$）；而在 $5.0 - 5.2\text{ eV}$ 的信号（绿色实线）振荡趋势则与 $^1A_u$ 的布居极大值时间同步（垂直绿色虚线，约 $55\text{ fs}$ 与 $95\text{ fs}$）。这直接证实了即便在 Dyson 轨道重叠严重的超快光谱中，利用选择性能量探测窗口提取单个激发态布居演化的可行性。

3. 亮态复苏的否定：

   • 图 1E（及图 S2D）的积分曲线在 $50\text{ fs}$ 之后保持平坦，未见任何明显的“二次抬升”。这以极具说服力的高精度模拟结果宣告：吡嗪激发态并不存在传统低维模型中预测的亮态在 $80\text{ fs}$ 之后的强复苏（Revival）行为，与实验谱线的观测结论高度一致。

2.3 氮/碳 K 边时间分辨 X 射线吸收光谱 (TR-XAS)

为了展示多态模拟方法在更先进的超快光源（如 X 射线自由电子激光，XFEL）实验设计中的指导价值，本工作计算了吡嗪在 $N$-edge（氮原核 K 边）和 $C$-edge（碳原核 K 边）的 TR-XAS 信号。

• 在 $N$-edge 探测（$393 - 397\text{ eV}$ 范围，见图 2）：

  • 亮态 $^1B_{2u}$ 的光谱峰主要分布在 $393.5\text{ eV}$ 和 $395.1\text{ eV}$；暗态 $^1B_{3u}$ 和 $^1A_u$ 的主峰分别位于 $394.8\text{ eV}$ 和 $396.5\text{ eV}$。由于在这一能量范围内核心激发态（Core-excited states）具有显著的单激发特征，计算中使用了基于 CVS 技术的 CC3 高阶方法恢复了核心空穴（Core hole）相关能。
  • 通过对理论能谱进行不同窗口的积分提取（图 2D-E），清晰复现了实验中观测到的强周期振荡（量子 Beats）。在 $394.3 - 395.1\text{ eV}$ 探测窗口中，$^1B_{3u}$ 的振荡峰值位于约 $74\text{ fs}$ 和 $138\text{ fs}$，对应周期的相干时间约为 $64\text{ fs}$。而在 $396.0 - 396.9\text{ eV}$ 探测窗口中，$^1A_u$ 在 $55\text{ fs}$ 和 $90\text{ fs}$ 显示了动力学峰，这为三态弛豫机制提供了不容置疑的光谱学铁证。

• 在 $C$-edge 探测（$279 - 283\text{ eV}$ 范围，见图 3）：

  • 碳边光谱主要由 $^1B_{2u}$（在 $280\text{ eV}$ 附近产生极亮的光谱带）在超快尺度下的消逝主导。通过重构 $279.5 - 280.3\text{ eV}$ 能段的光谱积分（图 3D），其超快消逝表现为 $\tau \approx 27\text{ fs}$ 的特征寿命，并在约 $40\text{ fs}$ 处表现出一个短暂的平台阶梯（Shoulder）。这一结论得到了实验（Ref. 21, 图 3E）的定性支撑，证明了碳边 XAS 能谱对亮态衰减动力学极度敏感。

2.4 动力学布居演化与核波包运动机制

通过对多态核动力学波包进行透彻地分析，提取了最本质的微观动力学图景（见图 4A 中的二阶 diabatic 种群演化）：

• 亮态衰减：在光激发瞬间，$100\%$ 种群布居于 $^1B_{2u}(\pi\pi^*)$（蓝色曲线）。由于强烈的振动耦合作用，在最初的 $25\text{ fs}$ 内，该状态发生雪崩式衰减，并在约 $50\text{ fs}$ 时其布居数已降至 $10\%$ 以下，伴随着单指数寿命约 $34\text{ fs}$（在波包真实演化层面）。这一种群向 $^1B_{3u}$（红色曲线）与 $^1A_u$（绿色曲线）分流。
• 量子 Beats 的本质：在 $50\text{ fs}$ 以后，亮态种群宣告枯竭，此后体系主要局限于暗态 $^1B_{3u}$ 和 $^1A_u$。这两个状态的二阶种群呈现出了极为规则的反相干（Out-of-phase）振荡特征（Beats）：当 $^1B_{3u}$ 种群攀升至极大值（约 $60\%$）时，$^1A_u$ 的种群降至极小（约 $20\%$），反之亦然。振荡周期在 $30\text{ fs}$ 到 $56\text{ fs}$ 之间，这正是由于体系核波包在第一激发态（$S_1$）的绝热势能面上，在 $^1B_{3u}$ 局域极小值（Min）与 $^1A_u$ 过渡态（TS）之间的简正模式（主要是 $\nu_{8a}$）路径上来回强力振荡导致的。这一量子运动机制在图 4D-E 的平均位移曲线中得到了无可争议的印证。

3. 代码实现细节、复现指南与开源生态

3.1 开源软件包核心：$e^T$ 程序与 FMS90

进行本项研究所使用的计算软件生态主要由两个开源的高性能计算库构建：

1. $e^T$ 程序（$e^T$ program）：作为电子结构的核心引擎。这是一个用 Fortran 编写的、专注于大规模含时与激发态耦合簇计算的开源程序。本研究所用的多态相似性约束耦合簇（Multistate SCCSD）方法和 CVS-CC3 核心激发能方法已被开发并集成于该代码的开发版本中。
   • 开源仓库：$e^T$ Program GitHub Repository
2. FMS90（Full Multiple Spawning 90）：作为量子核动力学演化的主干，控制 TBFs 高斯波包的经典运动、量子干涉和产卵机制。
   • 开源仓库：FMS90 相关代码与 AIMS 核心算法可参考 Martínez 组开发的分支 FMS90 Source

3.2 接口实现细节与输入文件模板设计

电子结构程序 $e^T$ 与动力学程序 FMS90 之间的耦合是通过一个基于 TCP/IP 套接字（Socket）或文件读写的自适应接口实现的。接口程序读取来自 FMS90 发送的高斯基函数中心几何坐标，在运行过程中动态做出自适应约束（SCCSD）判断。

在执行非绝热动力学复现时，首先要在 $e^T$ 的输入文件中指定相似约束阈值。一个典型的 $e^T$ 运行输入配置模板如下：

# eT program input file - Multistate SCCSD energy and gradients
# Adaptive threshold calculation

# 系统分子声明
molecule
  charge = 0
  multiplicity = 1
  geometry = {
    N   5.415992742468  -2.693688539865   0.000000000000
    N   5.415992990121   2.693688540299   0.000000000000
    C   7.556648910809  -1.327086026150   0.000000000000
    C   7.556649032784   1.327085828427   0.000000000000
    C   3.275336654687  -1.327085804547   0.000000000000
    C   3.275336777956   1.327086002459   0.000000000000
    H   9.340593105527  -2.378688274083   0.000000000000
    H   9.340593324553   2.378687912089   0.000000000000
    H   1.491392456703  -2.378688031330   0.000000000000
    H   1.491392676646   2.378688393153   0.000000000000
  }
end molecule

# 基组定义
basis
  default = cc-pVDZ
end basis

# 电子结构计算控制
method
  ccsd
end method

# 激发态与非绝热设置
cc_excited
  n_states = 3                 # 计算最低三个 singlet 激发态
  multistate_sccsd = true      # 启用自适应多态相似性约束耦合簇
  sccsd_threshold = 0.000368   # 准简正阈值，单位为 Hartree (对应 0.01 eV)
  gradients = true             # 启用受力梯度计算
  nacv = true                  # 计算非绝热导数耦合矢量
end cc_excited

在 AIMS 动力学输入文件 fms90.in 中，控制波包产生和演化的参数配置应当包含以下关键参数：

&dynamics_control
  N_states = 3                   ! 电子激发态数目
  Initial_state = 3              ! 初始激发亮态 (S3, 在 FC 处对应 1B2u)
  Timestep = 20.0                ! 默认演化积分步长 (a.u.，约 0.48 fs)
  Coupling_timestep = 5.0        ! 强耦合区域积分步长缩减至 5.0 a.u. (约 0.12 fs)
  Max_time = 8000.0              ! 总模拟时长 (a.u.，约 200 fs)
  Spawning_threshold = 20.0      ! 产卵耦合矢量模长阈值 (a.u.)
  Min_spawn_population = 0.05    ! 触发产卵的最低分子轨道种群分数 (5%)
  Biorthogonal_coupling = .true. ! 采用双正交耦合簇核准方案
/

3.3 动力学复现指南与计算资源开销评估

动力学复现步骤指引：

1. Wigner 采样：在 $S_0$ 势能面（CCSD/cc-pVDZ 级别）下运行解析频率计算（Table S2）。根据所得的 30 个简正模式，在 0 K 分布下抽取 100 个初始微观构型。根据计算的 UV-Vis 吸光能谱（图 S26），筛选出跃迁能落在 $4.56 - 5.13\text{ eV}$ 窗口内的 20 个物理初始构型。
2. 波包演化：以这些构型作为初始条件，初始波包放置于亮激发态（$S_2$ 或 $S_3$）。使用 FMS90 结合 $e^T$ 作为计算后端，自适应调用多态 CCSD/SCCSD。若在积分步出现不收敛（通常在离锥交叉极近处），接口会调用原子微位移辅助收敛算法（S11），即沿简正模式反方向微调构型使 DIIS 重新收敛，进而返回交叉点（Convergence stabilizer via small displacement steps）。
3. 数据提取与拓宽：将 20 个初始动力学演化所产生的所有 TBFs 的含时布居数和跃迁极化进行加权。使用 incoherent approximation（不相干近似），调用自编 Python 脚本读取 et_output 中的 Dyson 轨道极化和激发能，应用二维 Gaussian 函数进行能谱拓宽，重构出含时 TR-PES 与 TR-XAS 谱图。

计算资源开销评估：

• 单步计算开销：在常规 24 核（如 Intel Xeon Gold 6248R @ 3.0GHz）计算节点上，一个构型点的 CCSD/cc-pVDZ 的激发能、梯度和 NACV 联合计算约需 3 - 5 分钟。在激活 SCCSD 约束的情况下，由于要解拉格朗日平稳方程，单步计算时间增加约 30% - 50%。
• 总计算吞吐量：本研究所模拟的总轨迹时间为 200 fs。由于产卵机制，动力学过程中 TBF 数量从初始的 20 个膨胀到最终的 445 个。这导致总共演化了数十万步电子结构计算。这是一个计算密集的项目，累积消耗的计算资源约为 150,000 核时（CPU core-hours）。该方法目前难以直接应用于更大的多环芳烃分子，但作为小分子的精准 Benchmark 是极其卓越的。

4. 局限性深度剖析与关键文献

4.1 方法学与计算局限性批判

尽管本工作在吡嗪非绝热动力学和光谱重构上取得了前所未有的定量精度，作为一门严谨的学术成果，它依然存在以下几处方法学局限：

1. 双基组与核心激发态能级处理的脱节（Core-Valence Separation, CVS 的局限）： 对于价键激发态的动力学演化，使用的是较小的 cc-pVDZ 基组。这一基组在描述价激发态时能给出合理的几何构型，但是，在计算 TR-XAS 这种高精度核心激发态光谱时，cc-pVDZ 会引入显著的轨道截断误差，导致理论激发能相比实验整体高出数十 eV。虽然本工作通过应用整体能移（$N$-edge 整体平移 $-1.58\text{ eV}$，见图 2；$C$-edge 整体平移 $-1.66\text{ eV}$，见图 3）在光谱特征上取得了完美符合，但未能实现从头算上的绝对能量精准匹配。在未来的动力学中，有必要采用大基组（如 aug-cc-pVTZ）并混合核激发基组以达到自洽。
2. 单双激发近似（CCSD）对高阶多电子激发的缺失： 在 K 边吸收光谱的高能端（396-397 eV 探测区之外，见图 2C 实验图），实验观测到了数个高度重叠的谱峰。由于这些激发本质上是双电子激发主导（Double-excitation dominated core-excited states），而在 CCSD 级别下，双激发能的分辨率非常粗糙（通常高估 $1 - 2\text{ eV}$）。为了在理论上捕捉这些谱峰，必须在核心激发态计算中引入具有高阶相关性的三激发（CCSDT）或甚至四激发理论。本工作中作者尝试利用 CCSDT（图 S13）进行了修正，虽有改善但仍未彻底根治，证明更高精度的多参考（Multi-reference）耦合簇或者 CCSDTQ 才是描述核心高阶激发的最理想方案。
3. 经典核高斯近似的半经典局限： 尽管 AIMS 基于量子产卵，但高斯基函数的宽、高在演化中被硬性约束为“冷冻高斯（Frozen Gaussian）”，无法动态拓宽和变形。这意味着长期演化后（> 150 fs），量子隧穿和极端非谐性引起的波包严重畸变无法被 AIMS 完美重现，这或许是 100-125 fs 处模拟振荡周期与实验不一致的深层力学原因。

4.2 核心参考文献

为了完整把握本领域的演化脉络，强烈推荐精读以下奠基性文献：

1. 多态量子动力学基准（吡嗪光弛豫）：
   • Sala, M., Lasorne, B., Gatti, F. & Guérin, S. The role of the low-lying dark $n\pi^*$ states in the photophysics of pyrazine: a quantum dynamics study. Phys. Chem. Chem. Phys. 16, 15957–15967 (2014).
     (该文献系统阐释了三态光弛豫机制在低维模型下的基础力学图像。)
2. 相似性约束耦合簇（SCCSD）理论的创立：
   • Kjønstad, E. F. & Koch, H. Resolving the notorious case of conical intersections for coupled cluster dynamics. J. Phys. Chem. Lett. 8, 4801–4807 (2017).
     (这篇奠基性论文首次提出了通过拉格朗日乘子约束突破非厄米哈密顿量锥交叉断裂限制的方法。)
3. 从头算多态产卵（AIMS）方法论：
   • Ben-Nun, M. & Martínez, T. J. Ab initio quantum molecular dynamics. Adv. Chem. Phys. 121, 439–512 (2002).
     (Martínez 教授撰写的经典综述，全面讲解了 FMS90 和 AIMS 动力学中轨迹基函数的物理推导。)
4. 时间分辨光电子谱（TR-PES）测量实验：
   • Karashima, S., Humeniuk, A. & Suzuki, T. Vibrational motions in ultrafast electronic relaxation of pyrazine. J. Am. Chem. Soc. 146, 11067–11071 (2024).
     ( Suzuki 课题组提供的极高信噪比超快 PES 实验光谱，本工作与其进行了直接的定量对比。)

5. 补充理论推导与技术细节

5.1 绝热到非绝热(Diabatization)变换的数学推导

在动力学分析（图 4A）和光谱拆解中，必须将沿绝热势能面（Adiabatic states $S_1, S_2, S_3$）演化的波包投影到具有化学意义的非绝热态（Diabatic states $^1B_{2u}, ^1B_{3u}, ^1A_u$）上。本工作发展了一种利用激发态**电子跃迁强度（Transition Strengths）**进行二阶变换的准正交（Orthogonal Transformation）数学方案：

在 Franck-Condon（FC）参考几何点，由于空间对称性，$S_1$ 态纯粹由 $^1B_{3u}$ 主导，$S_2$ 由 $^1A_u$ 主导，而 $S_3$ 由 $^1B_{2u}$ 主导。定义一个在任意位移构型 $\mathbf{R}$ 处的跃迁强度矩阵 $\mathbf{F}(\mathbf{R})$，其每一列对应一个绝热态的跃迁极化强度的平方根：

在 FC 构型处，该矩阵由于跃迁偶极禁戒性，具有极其简洁的形式：

这对应于 $S_1(^1B_{3u})$ 仅具有沿 $x$-方向的跃迁，而 $S_3(^1B_{2u})$ 仅具有沿 $y$-方向的跃迁。对于任意动力学几何点 $\mathbf{R}$，利用此对称性构建如下正交变换矩阵 $\mathbf{U}(\mathbf{R})$：

1. 首先提取第二行元素： $$\mathbf{V}_{B2u} = \mathbf{F}_{2, :}(\mathbf{R}) = \begin{pmatrix} F_{21} & F_{22} & F_{23} \end{pmatrix}$$ 将其进行归一化，得到 $^1B_{2u}$ 在绝热态空间中的混合系数： $$\mathbf{u}_1 = \frac{\mathbf{V}_{B2u}}{\|\mathbf{V}_{B2u}\|}$$
2. 类似地，提取第一行元素作为 $^1B_{3u}$ 的基础投影向量： $$\mathbf{V}_{B3u} = \mathbf{F}_{1, :}(\mathbf{R}) = \begin{pmatrix} F_{11} & F_{12} & F_{13} \end{pmatrix}$$
3. 为了确保变换正交化，使用施密特（Gram-Schmidt）正交化方法将 $\mathbf{V}_{B3u}$ 投影，使其严格垂直于 $\mathbf{u}_1$： $$\mathbf{V}_{B3u}' = \mathbf{V}_{B3u} - (\mathbf{V}_{B3u} \cdot \mathbf{u}_1)\mathbf{u}_1$$ 归一化得到第二个转换向量： $$\mathbf{u}_2 = \frac{\mathbf{V}_{B3u}'}{\|\mathbf{V}_{B3u}'\|}$$
4. 最后，通过向量叉乘（Cross Product）自动得到第三个正交方向，即暗态 $^1A_u$ 的投影向量： $$\mathbf{u}_3 = \mathbf{u}_1 \times \mathbf{u}_2$$

最终的正交变换矩阵即为 $\mathbf{U}(\mathbf{R}) = \begin{pmatrix} \mathbf{u}_1 \\ \mathbf{u}_2 \\ \mathbf{u}_3 \end{pmatrix}$。通过这一矩阵，可以在每个时间步将含时绝热密度矩阵（Populations in adiabatic representation）转换为精确的非绝热态布居数：

这一优雅的数学推导，不仅比传统的波函数重叠二阶变换更易编码，且计算开销可忽略不计。

5.2 简正振动模式位移投射分析

非绝热动力学轨迹的复杂性在于 24 维势能面内多重振动运动的交织。为厘清决定光弛豫的核心模式，对 AIMS 的平均核波包位移在 Franck-Condon 简正模式上进行了几何投射。简正分析表明：

• 对称化呼吸模式 $\nu_{6a}$（$A_g$ 模式，频率 $606\text{ cm}^{-1}$）： 在激发后的前 $40\text{ fs}$ 内，核波包沿着该振动方向出现剧烈的极化位移（见图 4E 橙色实线，位移值由 0 陡峭跌落至 -1.0）。由于该模式直接调节了亮态 $^1B_{2u}$ 和两个暗态势能面的能级差，正是这一快速的单向振动运动将核波包迅速推向了 $S_2/S_1$ 锥交叉窗口，导致了亮态的急剧消逝。在 $50\text{ fs}$ 以后，该模式位移围绕其激发态新极小值（$-0.7$ a.u.）进行微弱地弛豫和耗散振荡。
• 对称化伸缩模式 $\nu_{8a}$（$A_g$ 模式，频率 $1670\text{ cm}^{-1}$）： 与 $\nu_{6a}$ 在前期占统治地位不同，$\nu_{8a}$ 的平均位移在全时域（$0 - 200\text{ fs}$）呈现了完美的正弦规律性振荡（见图 4D 蓝色实线），振荡周期正好对应 $34\text{ fs}$。更加令人惊叹的是，$\nu_{8a}$ 位移的极大值时间点与 $^1A_u$ 的种群极大值（垂直绿色虚线）高度吻合，极小值时间点则对应 $^1B_{3u}$ 的极大值（垂直红色虚线）。这一完美的匹配向人们证实：$S_1$ 势能面上的绝热演化本质上是核波包在 $\nu_{8a}$ 简正模式上进行的大幅度、不简谐来回穿梭振动。这一简正振动正是调制暗态量子相干 Beats 的分子弹簧。

5.3 核心激发态模拟中的 CVS-CC3 理论

X 射线吸收光谱（XAS）涉及分子最内层的 $1s$ 轨道电子向未占据轨道跃迁。在理论上，核心激发态的能量极高，极易与外层价激发态连续区产生强烈的共振混合。为了解决激发态计算中高阶态根缺失的严重难题，本工作引入了核心-价键分离（Core-Valence Separation, CVS）近似。

CVS 近似基于投影算符技术，将相似哈密顿量矩阵中不涉及核心空穴（Core hole）激发的矩阵元强制清零：

其中投影算符 $P_{CVS}$ 仅保留那些在算符激发指数中包含核心轨道指标（例如核心轨道 $1s$ 上的退激或激发）的振幅算符。在求解 $e^T$ 方程组时，CVS 技术能够完全隔离出纯净的核心激发态，杜绝了其与高能价激发态的混合，从而允许在 $S_1, S_2, S_3$ 的动力学几何轨迹上精确计算 K 边核心激发谱线能级。

为了恢复高精度的核心轨道弛豫（Core relaxation）和电子相关能，在核心态使用了具有近似三激发的 CC3 理论，其激发能量精度相比常规 EOM-CCSD 提升了数个数量级。这使得本研究重构的 TR-XAS 光谱无论在理论严谨性还是在物理真实性上，均树立了非绝热动力学领域的全新行业标杆。