基于半经验方法的激发态伞形采样

《JACS Au》：Umbrella Sampling for Excited States Using a Semiempirical Method

　　

伞形采样（US）被广泛用于探索具有高能量壁垒的基态反应路径以及识别稳定的分子构象。研究人员将伞形采样扩展至激发态，利用电子态之间的能量差作为集体变量（CV），不仅可以定位交叉缝，还能系统性地探索通向锥形交叉（CI）的激发态弛豫路径。基于这一方法，研究人员提出了研究激发态弛豫机制的完整工作流程。作为案例研究，研究人员将该方法应用于推-拉式芪衍生物4-(N,N-二甲氨基)-4'-硝基芪（DANS），使用基于正交化模型2（OM2）方法的半经验多参考组态相互作用（MRCI）模拟光异构化。结合势能面（PES）扫描、非绝热动力学模拟和激发态伞形采样，研究人员识别并表征了反式和顺式异构体的五个不同锥形交叉。伞形采样通过确定将这些锥形交叉与激发态极小值分离的自由能壁垒，评估了这些锥形交叉的可及性。此外，它提供了沿弛豫路径的自由能和熵演化 insight，突出了系统接近锥形交叉时的关键热力学特征。

光化学反应自20世纪初Ciamician的开创性实验以来一直吸引着研究人员的关注。（1,2）对光化学反应机制的早期解释基于非交叉规则，该规则认为相同对称性的绝热电子态不能相交，而是产生回避交叉——这一原理严格来说仅对双原子分子有效。（3,4）在多原子分子中，绝热态之间的交叉，即锥形交叉（CI），可以发生，使得不同绝热电子态之间能够进行超快、无辐射跃迁。（5-7）在这些点上，由于电子波函数随核构型快速变化，玻恩-奥本海默近似失效。CI在光化学中起着关键作用，使分子能够访问通过单纯热活化无法达到核构型，例如在光异构化反应中。（8-11）此外，CI提供了一种无辐射的、通常更快的替代途径回到基态，防止长寿命、高反应活性激发态的积累，从而确保例如核酸的光稳定性。（12,13）因此，理解CI的结构、能量学和反应路径至关重要。鉴于其核心作用，数十年来方法开发致力于CI的识别和表征，更多细节总结在第2节中。然而，大多数已建立的CI定位方法都是针对气相中的孤立分子（或隐式溶剂模型中的分子）设计的，对于大型生色团或嵌入在显式环境中的情况，如溶剂、生物大分子、蛋白质或表面，这些方法并不有效。随着分子系统尺寸的增长，单个优化结构（无论是全局还是局部极小值）对完整平衡构象系综的代表性越来越差，（14,15）使得广泛的构象采样对于准确捕捉热力学和动力学行为至关重要。这一局限性同样适用于电子基态和激发态，因此也延伸到CI附近区域，其中多种几何构型可能都对系统的弛豫有贡献。在这种高维图景中，势能面的崎岖性使得单个优化的CI变得无关紧要。

本研究引入了一种激发态伞形采样方法，克服了基于优化和单几何描述方法的局限性，允许沿选定反应坐标进行全面的统计采样。该框架特别适用于表征交叉区域，在这些系统中，单个或少量优化的最低能量锥形交叉（MECIs）或最低能量交叉点（MECPs）不足以捕捉控制非绝热行为的相关核构型。通过在激发态势能面上沿选定的集体坐标进行受控采样，激发态伞形采样能够高效地探索近简并区域的构象系综。这种系综视角特别适用于大型生色团或存在显式环境的情况，其中非绝热弛豫由构型分布而非孤立的交叉点控制。该方法已在开源的SHARC（包含任意耦合的表面跳变）分子动力学套件中实现。（16）研究人员提供了计算方案的详细描述，将其确立为研究近简并区域激发态动力学的通用且稳健的方法。为证明其实用价值，研究人员以4-(N,N-二甲氨基)-4'-硝基芪（DANS，见图1）为例，逐步展示了应用过程。DANS是一种推-拉式光开关，（17-19）可通过荧光辐射衰变到电子基态，或通过导致光异构化的CI进行无辐射衰变。尽管DANS在多种溶剂中已有广泛的实验表征（20-28），并在众多应用中具有重要意义，（29-33）但其弛豫机制仍未完全解决。为此，研究人员采用伞形采样，以能量差为反应坐标，研究从激发态极小值到与电子基态交叉缝的弛豫路径。模拟使用基于正交化校正方法2（OM2）（34-36）的半经验多参考组态相互作用（MRCI），在MNDO代码（37）中实现并与SHARC接口。（16,38）OM2/MRCI提供了描述顺-反异构化所需的多组态框架，以及采样大型、柔性光活性系统相空间所需的效率。

本研究首先介绍了寻找CI的现有计算方法，将其分为三个概念性"时代"：第一时代（约1980年代至2000年代初）聚焦于优化MECI，即CI缝上能量最低的点；第二时代（约2000年代中期至2010年代中期）发展识别多个局部MECI的方法；第三时代（约2010年代中期至今）采用增强采样策略随机探索CI缝附近的构象空间。梯度投影法、拉格朗日乘数法和罚函数法是三种主要的MECI优化技术。随着计算工具和算法的发展，非绝热激发态动力学变得可行，提供了机制洞察以及不同弛豫通道被访问的时间信息。

研究的核心方法论是激发态伞形采样。为构建势能均值力（PMF）和自由能曲线，需要沿所选反应坐标或集体变量采样构型。由于垒阻碍探索，研究人员采用伞形采样方法，通过引入偏置势能来聚焦特定区域的采样。在重叠窗口中进行模拟后，使用加权直方图分析方法（WHAM）重建无偏的自由能曲线。本研究选择相交态之间的能量差作为集体变量，这一选择与Marcus理论相关，该理论以能量差ΔE_IJ作为反应坐标来框架化电子转移过程。该方法已在SHARC分子动力学程序套件中实现。

电子结构计算采用OM2/MRCI方法，基于OM2半经验哈密顿量。OM2方法属于忽略双原子微分重叠（NDDO）族的半经验电子哈密顿方法，通过正交化校正改进了Pauli交换。分子轨道计算使用标准或浮动占据分子轨道自洽场技术，随后在第二步使用MRCI计算电子态。SHARC接口可提取电子能量、梯度、非绝热耦合矢量和偶极矩，并计算波函数重叠以传播电子波函数。

作为案例研究，研究人员研究了DANS在气相中的光化学行为。DANS的反式和顺式异构体的可能弛豫路径通过最小能量扫描确定。对于反式DANS，主要弛豫路径对应于序列1→3→4→5，即异构化路径：弛豫到S₁极小值（S₁^t_min）后，越过0.33 eV的壁垒到达S₁/S₀^iso MECI，导致异构化。该MECI被称为"垂直"几何构型，因为碳桥处的HCCH二面角约为90°。对于顺式DANS，类似的异构化路径也存在，且壁垒稍低（0.28 eV）。此外，研究人员还发现了三条替代路径，涉及硝基基团的锥形化和环化，这些新路径是通过非绝热动力学模拟和激发态伞形采样发现的。

非绝热动力学模拟揭示了不同的弛豫时间尺度：反式异构体的弛豫时间常数为21.7±2.6 ps，而顺式异构体仅为1.9±0.3 ps，约快一个数量级。这归因于顺式异构体从激发态极小值到CI的能垒更低。在反式异构体中，61%的布居在20 ps后到达基态，其中42%异构化为顺式；在顺式异构体中，96%的轨迹在5 ps内到达基态，其中不到一半转化为反式。

激发态伞形采样提供了静态扫描无法获得的定量信息。对于反式异构体，自由能曲线显示在能量差约2.3 eV处有~8.2 kcal/mol（0.36 eV）的壁垒，与静态计算结果（0.33 eV）吻合。在顺式异构体中，异构化和环化路径显示出明显不同的自由能行为：异构化路径在大能量差处呈现浅势阱，而环化路径则显示更明显的能量增加。特别值得注意的是，当系统接近CI时，自由能急剧上升，同时熵显著下降，这是由于在CI附近系统被限制在3N-8个自由度，分支平面矢量（g和h）的简并条件严重约束了核几何构型，限制了可及的构象空间。

研究还验证了能量差作为集体变量相对于几何坐标（如二面角）的优势。当使用HCCH二面角作为CV时，采样往往只覆盖特定区域，可能遗漏其他重要的CI；而能量差CV能更精确地引导采样到达实际的交叉区域。

本研究的主要技术方法包括：（1）半经验OM2/MRCI电子结构方法，使用MNDO代码实现，采用四电子四轨道的激发空间（HOMO-1、HOMO、LUMO、LUMO+1），包含六个参考组态及所有单双激发；（2）SHARC非绝热动力学程序套件，用于表面跳变模拟和伞形采样实现；（3）Wigner采样生成初始条件；（4）加权直方图分析方法（WHAM）重建自由能曲线；（5）基于Dietschreit等人开发的公式计算自由能、内能和熵，引入de Broglie热波长修正以消除非线性坐标变换带来的畸变；（6）nudged elastic band计算和过渡态优化；（7）image-dependent pair potential方法进行路径插值。

研究结果部分详细展示了以下发现：静态势能面扫描揭示了反式和顺式DANS的五条不同弛豫路径，包括已知的异构化路径和新发现的硝基锥化及环化路径。非绝热动力学结果显示了不同的时间尺度和量子产率，并识别了多个CI构型族。激发态伞形采样的热力学剖面揭示了异构化路径中第二极小值的存在，以及接近CI时自由能的急升和熵的骤降。配置分布分析证明了能量差CV相对于几何CV的优越性，能够全面采样多种CI变体。

讨论部分强调了激发态伞形采样的优势，特别是在处理大型系统和显式环境时的重要性。研究人员指出，对于小型气相分子，静态势能面计算可能提供合理的近似，但对于更大的系统，采样方法对于正确考虑熵贡献和识别最相关的弛豫路径至关重要。

研究结论总结道：本工作提出了研究光弛豫过程的工作流程，以光开关DANS为例。应用OM2/MRCI方法进行电子结构计算，首先计算静态势能面扫描，随后进行基于表面跳变的非绝热动力学模拟，最后通过新实现的激发态伞形采样探索从S₁极小值到各种CI的弛豫路径。后者在静态势能面扫描不足且需要适当采样弛豫路径的情况下提供了深入认识。

应用于气相DANS的完整方案确定了五个CI：可从反式和顺式异构体到达的异构化CI；反式异构体的NO₂锥化CI；以及顺式异构体的三个环化CI（两个含NO₂扭转，一个不含）。对于反式异构体，研究证实了实验研究（22,23,25）所提示的与NO₂基团振动模式相关的CI的存在，尽管该路径的势能面计算壁垒高于异构化路径。对于顺式异构体，发现了与母体化合物芪一致的环化路径（178,181,182,187-191），以及两个新的（此前未报道的）伴随NO₂扭转的环化路径，这些路径被证明是竞争的弛豫通道。反式异构体的弛豫显著慢于顺式异构体（21.7 ps vs 1.9 ps），这可归因于顺式情况下分离S<sub>1极小值与CI的能量壁垒更小。

使用能量差作为CV的激发态伞形采样，结合文献（118）中描述的自由能和熵形式，提供了沿弛豫路径性质的详细认识。计算的自由能曲线也揭示了异构化路径的独特特征，突出显示了它与其他弛豫路径的区别，并提供了相应熵曲线的认识。对于两种异构体的异构化路径，越过能垒后S1态中出现第二极小值。另一个显著特征是，当系统接近CI时，自由能迅速增加伴随熵的降低，表明在核经典处理的模拟中，电子态之间的跃迁将在小但非零的能量差处发生（76）。

本工作展示了激发态伞形采样在探索弛豫路径中的实用性，特别是使用能量差而非几何参数作为CV时。能量差为跟踪弛豫过程提供了自然且直接的手段，允许沿反应路径而非整个构象空间进行定向采样。此外，它提供了对大小能量差区域的采样控制。伞形采样可能特别适合非绝热增强采样（112），因为像元动力学这类方法施加态特定的偏置势能，使其不太适合系统在多个电子态之间演化的情况。

对于自由度较少的小型气相分子，静态势能计算通常能提供对自由能曲线的合理近似，如DANS中所观察到的。相比之下，更大的系统——特别是在存在环境效应的情况下——展现出更为复杂的弛豫图景，具有众多的激发态极小值、能垒和MECI。在这些情况下，基于采样的方法如伞形采样变得至关重要，以正确考虑熵贡献并识别最相关的弛豫路径。

类似于其在溶剂化和生物系统基态采样研究中的既定作用（199-203），将伞形采样扩展到激发态过程，为阐明复杂环境中的光化学机制提供了有力策略，其中环境相互作用可以将路径重定向至CI并最终决定弛豫结果。在这方面，本研究提出的方法特别适用于QM/MM框架内的实现，从而能够对环境效应进行显式描述。基于先前对无定形二氧化硅表面上DANS的研究（164,204），未来的工作将致力于研究玻璃表面上的动力学。</sub>