综述：《RRKM理论实用指南：其简化多井版本AWATAR及辐射过程的主方程建模方法》

《Physical Chemistry Chemical Physics》：A practical guide to RRKM theory, its simplified multi-well version AWATAR and master equation modelling of radiative processes

　　在超高真空条件下，气相中的单分子反应和重排不能通过平衡热力学正确描述，因为缺乏热浴，并且对于小系统来说，辐射热化过程较慢。这就是Rice Ramsperger Kassel Marcus理论（通常称为RRKM理论）的领域，在该理论中，单分子反应速率常数由系统的内部能量E以及最小结构和过渡结构的能量可及的量子态决定。我们最近扩展了经典的RRKM理论，以描述在势能面上具有多个能量较低最小值的实体的解离过程，这种方法被称为AWATAR（All Wells And Transition structures Are Relevant）。这种方法用于对由辐射过程引起的解离进行主方程建模，例如黑体红外辐射解离（BIRD）或红外多光子解离（IRMPD）。在本教程综述中，我们首先介绍了RRKM和AWATAR计算的示例。辐射过程被实时建模，也可以通过将分子集合的种群收敛到稳定状态来建模。文中展示了BIRD和IRMPD实验的模拟示例。用于模拟的AWATAR 1.0程序也已提供。所有示例在补充信息（SI）中都有注释。

**1. 引言**

气相中的单分子异构化和解离反应在描述复杂过程（如燃烧或大气化学）中起着重要作用。1 通常也将双分子气相反应描述为单分子反应，因为它们的结果取决于碰撞复合物的命运，即是否回解离为反应物或发生辐射结合2,3 或产物形成。纯单分子反应所需的能量可以通过光子吸收或高能碰撞提供。在质谱仪中进行的气相离子反应性实验，在定义良好的条件下非常适合使用单分子统计速率理论进行定量分析。例如，使用引导离子束技术的碰撞诱导解离（CID）实验4，4 流管实验中速率系数的热依赖性建模5，黑体红外辐射解离（BIRD）6–10，生物分子离子解离11–14，真空紫外（VUV）区域的解离光电离途径15，具有潜在天体化学意义的有机自由基阳离子的单分子分解16，金属簇阳离子上的催化相关基本步骤17，气相中的对映选择性反应18，或大气化学中的潜在离子途径19。通过动理论同位素效应的建模10，或者通过Rice–Ramsperger–Kassel–Marcus（RRKM）理论可以获得基本见解16，该理论用于模拟技术相关的等离子体过程21或确定气相离子的内部能量分布22,23。在红外多光子解离（IRMPD）光谱学中，最近利用RRKM理论来评估水簇异构体的动力学稳定性24,25。用于RRKM计算的程序包括Eyringpy26–28，KiSThelP29，Multiwell30，Arkane31，TUMME32，MSMC33，MESS34，MESMER35，CHIMERA36和Zhu-Hase RRKM程序37。

气相离子化学实验通常伴随着对反应势能面的量子化学计算38–40，使用诸如ORCA41，Gaussian43，Molpro44，Q-Chem45，Turbomole46等强大的软件包。当使用这些软件包计算振动频率时，其输出可以用来计算热化学值，如焓ΔH和吉布斯自由能ΔG。在溶液相化学中，特别是在生物化学中，吉布斯自由能对于描述反应轮廓非常有用，因为它们包含了熵，并提供了中间体的化学平衡的实际情况。然而，使用ΔG意味着恒定的温度和压力，即系统与无限热浴的耦合，这在凝聚相或气相中的高碰撞率情况下是不会出现的。在低压力下的气相实验中，例如质谱仪的高真空或超高真空条件下，这些条件通常不满足，因此ΔG值是没有意义的。在这种情况下，描述反应势能面最实用的热化学量是零点校正能量，无需应用任何热或熵校正。然后可以在RRKM理论的框架内，根据这些零点校正能量和系统可及的量子态信息，推导出特定异构体的种群数量、单分子解离或异构化速率常数。本教程适用于具有相当量子化学知识的初学者，比如从事质谱学或红外多光子解离（IRMPD）光谱学研究的学生，他们通过密度泛函理论（DFT）计算来了解所研究离子的性质、可能的重排和解离途径。本教程综述的重点是将RRKM理论应用于更复杂的系统，如簇或生物分子，特别是在典型实验设置中的单分子分解。除了对单个反应步骤的标准RRKM计算外，我们还将指导读者分析多井势能面，这在研究生物分子或簇离子时经常遇到。为此，我们最近引入了All Wells And Transition structures Are Relevant（AWATAR）作为RRKM理论的简化多井扩展47,48。我们在Github/Zenodo49上提供了我们为BIRD实验的主方程建模开发的AWATAR代码。该代码还允许进行标准的RRKM和AWATAR计算，并分别计算反应物井的态密度和过渡结构的态之和作为内部能量的函数。它还提供了用于BIRD和IRMPD实验的实时模拟的参数。提供了示例输入文件以获得实践经验。我们希望使读者能够使用这些方法并根据自己的研究调整这些示例。

**2. 计算方法**

所有RRKM和AWATAR计算都是使用AWATAR1.0进行的，该程序作为可执行文件在Github/Zenodo上提供，同时提供源代码。49 快速入门指南介绍了程序的控制方法。所有参数和功能在用户手册中有详细描述。这两份文档也可在Youtube/Zenodo上找到。49

AWATAR1.0需要输入系统势能面上局部最小值（Wells）和过渡结构（TSs）的能量和振动频率计算。计算模型的选择，例如在DFT计算中的泛函和基组，由读者自行决定。然而，结果通常对相对能量和能垒高度非常敏感。在讨论AWATAR1.0的结果时，必须始终考虑量子化学计算的不确定性，例如通过使用不同的泛函和/或基组重复计算。该程序可以直接读取Gaussian43的输出文件。其他量子化学软件的输出需要转换为最小输入文件，这在手册中有说明。在AWATAR1.0以及本综述中，能量以cm?1为单位给出，参见表1了解转换为kJ mol?1、eV和Hartree的换算因子。系统势能面上全局最小值的零点校正能量设置为0 cm?1。然后，系统的总能量E在分子的内部自由度和势能面上更高稳态的势能之间进行分配。表1

能量换算因子

cm?1
kJ mol?1
eV
Hartree

1000
11.96
0.124
0.004556

8359
100
1.036
0.03809

8066
96.49
1.000
0.03675

21?947
262.6
2.721
0.1

**3. 单分子反应**

**3.1. RRKM速率常数**

**3.1.1. 异构化**

气相中单分子反应的描述始于系统势能面（PES）上的相关稳定点。图1显示了一个异构化反应（1）的PES作为示例。Ta(OH)4+ ? TaO(OH)2(H2O)+

(1)

为了形成水配体，一个质子从一个氢氧配体转移到另一个氢氧配体上。过渡结构TS1是反应路径上的最高点。与虚频相关的运动是质子在反应物和产物之间的移动。从顶点开始，这是一个一阶鞍点，在PES上，质子可以朝能量较低的方向移动，如箭头所示。补充信息（SI）中提供了一个动画（https://animation1.gif）。如果有足够的能量，即E > E0(TS1)，复合体将正向和反向进行异构化（1）。使用I1、I2和TS1的振动频率，我们可以通过RRKM理论计算异构化（1）在两个方向上的单分子速率常数，使用局部最小值i的态密度ρi(E)和能量为E0,y的过渡结构TS y的态之和Ny(E ? E0,y)，公式（2）。

文件example1.rrk提供了使用AWATAR1.0进行这些计算的相关输入。快速入门指南和用户手册详细说明了如何设置这样的文件，以及如何使用控制按钮和参数面板。每个井和过渡结构的面板包含振动频率及其红外强度，以及结构的旋转常数及其以Hartree为单位的能量。这些数据是使用Gaussian在B3LYP/def2-TZVPP理论水平上通过量子化学计算获得的。根据Gilbert和Smith的详细描述的标准近似50，不对称转子被近似为对称转子，旋转常数由程序自动调整。这导致一个具有量子数J的二维转子和一个具有量子数K的一维转子。由于角动量守恒，J是恒定的，只有在三维转子中才能与振动自由度交换能量。因此，一个不对称转子井具有（3N ? 6）个振动自由度和一个旋转自由度，总共（3N ? 5）个活跃自由度。有关AWATAR1.0输入和输出Example1的详细信息，请参见SI。图1

Ta(OH)4+通过TS1异构化为TaO(OH)2(H2O)+，随后失去H2O，这可以通过紧密或松弛的TS来表示，计算是在B3LYP/def2-TZVPP理论水平上进行的。参见SI中的TS1和紧密TS的动画GIF。密度和态之和是用Gilbert和Smith描述的Beyer–Swinehart算法51数值评估的。直接状态计数仅适用于较为简化的振动频率，这些频率必须是能量量子（BS-Quantum）的整数倍，BS-Quantum在AWATAR1.0的参数面板中定义。BS-Quantum的典型选择是1.0 cm?1，但对于小系统可能需要更大的值以避免非物理行为。这是因为在非常精确的能量水平下，可能在总和振动和旋转模式中的能量时无法完全匹配。较粗的粒度（较大的BS-Quantum）增加了与总能量的组合匹配数量。建议通过改变BS-Quantum来测试稳健性，并确保计算态密度和态之和的结果没有显著变化。有关BS-Quantum参数的更详细讨论，请参见SI。对于态密度ρi(E)，算法计算能量E可以在井的内部自由度上分配的方式数量，然后除以BS-Quantum得到ρi(E)的单位（cm?1）?1。态之和的计算使用类似的算法，但其初始化方式不同，算法会补偿不同的BS-Quantum值，因此结果在很大程度上独立于此参数。详见Gilbert和Smith的书籍。50

图2显示了反应（1）的RRKM速率常数kiy(E)，这些常数仅在E > E0,y时才有非零值。接近阈值时，这些值迅速增加，并在更高能量处逐渐趋于平稳。反向反应的速率常数比正向反应高几个数量级，表明离子主要占据全局最小结构I1。图2

**3.1.2. 更深入地理解态密度**

RRKM速率常数kiy还说明了态密度的物理意义。如果我们观察通过TSAB连接的两种异构体A和B之间的平衡状态（反应式3），我们可以根据A和B的态密度（方程式4）计算平衡常数KAB(E)，作为正向反应和反向反应速率常数的比值，因为TSAB的态之和相互抵消了。需要注意的是，在没有热平衡的情况下，我们仍然使用微正则系综，因此这里定义的KAB(E)不能直接与ΔG相关联。平衡常数根据定义是处于异构体B和A的系统概率pB(E)和pA(E)的比值。本质上，这描述了一个单分子的化学平衡情况，该分子可以从异构体A反应到B，然后再返回A，因为其内能E高于过渡结构，见图1。

（3）

找到分子处于异构体A的概率pA(E)由方程式（5）给出。

这意味着通过比较给定异构体的态密度，可以很容易地确定在某个能量下哪种异构体是最可能的。特别是对于较高的能量E，这比异构体之间的相对能量更有信息量。图3(a)显示了I1和I2的态密度，随着能量的增加而迅速增加，当E约为30,000 cm^-1时，I1的态密度达到了10^20 (cm^-1)^-1。直到E约为1200 cm^-1，ρI1(E)并不单调增加。这不影响这里的RRKM计算，因为异构化（1）的障碍接近20,000 cm^-1。然而，如果分析通过低能量过渡结构的异构化，可能会遇到类似的数值问题。图3(b)显示了使用方程式（5）计算的I2的概率。I1的概率p(I1)始终接近1，因此没有显示。

(a) I1和I2的态密度；(b) 随能量变化的处于I2状态的系统的概率。方程式（5）描述了两种异构体的情况，但可以推广到任何数量的异构体。找到分子或簇处于任何能量可及异构体i中的概率pi(E)由态密度ρi(E)与所有异构体态密度之和的比值给出（方程式（6）。我们使用这个方程来计算簇实验中较高能级异构体的分布，例如金硅烷52或氯化钠簇47。

3.1.3. 解离：紧密过渡结构与松散过渡结构

RRKM理论也非常适用于描述解离反应。在许多情况下，分子的解离是通过逐渐达到解离阈值来进行的，因此没有具有虚频的过渡结构。仍然有两种方式可以对解离进行建模：使用量子化学计算来表示过渡结构（紧密过渡结构），或者使用所谓的松散过渡结构。紧密过渡结构需要在参数面板中指定解离能量，并用合适的振动模式替换反应坐标。为此，可以借助GaussView、Molden或ChemCraft等工具来直观检查振动模式，并识别解离模式。通常，这是一种拉伸模式，其中分离的碎片以显著的振幅相互振荡。对于TaO(OH)2(H2O)+而言，这一点在图1和作为补充材料的animation2.gif中有所说明。在AWATAR1.0中，势能井的频率计算可以被理解为紧密过渡结构，解离模式将从TS的频率列表中移除。此外，必须调整TS能量以反映解离能量。有关示例2.rrk中实现的详细信息，请参见补充材料。紧密过渡结构模型需要识别一个真实的反应坐标。在Ta(OH)4+的解离过程中，分子离子首先必须重新排列为TaO(OH)2(H2O)+，这意味着我们必须分两步描述解离过程：首先重新排列，然后解离。此外，当相互作用减弱时，对应于两个碎片相对运动的振动自由度（如弯曲模式）会显著减弱。可以通过任意更改合适模式的频率来考虑这一点，但这会增加结果的不确定性。这些问题可以通过使用松散TS模型来解决。在这里，碎片被视为完全分离的，因此它们可以独立旋转。在非线性多原子碎片的一般情况下，紧密过渡结构模型中的五个振动自由度变为旋转自由度，见表2和第3.1.1节的讨论。这通常会导致态总数显著增加，从而导致更高的RRKM速率常数。AWATAR1.0中松散TS模型的实现遵循Rodgers、Ervin和Armentrout描述的程序。对于金属-配体复合物或像水簇这样的弱结合系统，通常更倾向于使用松散TS模型。图1显示了TaO(OH)2(H2O)+失去水分子的松散过渡状态。表2

3.2. AWATAR速率常数

在上面描述的紧密和松散过渡结构模型中，我们讨论了解离前的系统重排。我们已经看到，可以通过方程式（6）计算出在特定异构体中找到分子或簇的概率。Ta(OH)4+通过失去H2O的有效RRKM解离速率常数是从TaO(OH)2(H2O)+的单分子速率常数I2计算得出的，乘以在给定时间具有这种结构的离子的比例。更一般地说，当我们结合方程式（2）和方程式（6）来计算系统通过特定过渡结构y从特定势能井i进行单分子反应的有效速率keff,y(E)时，我们会立即看到特定势能井i的态密度相互抵消，剩下的是所有势能井的态密度之和（方程式（6）的分母。这是RRKM理论核心思想的一个扩展，即速率常数由过渡结构以上的可访问量子态数量和最小值处的量子态密度决定。在异构体之间的异构化比通过过渡结构的反应更快的情况下，态密度之和决定了反应的概率。

（7）

这说明所有能量可及的势能井在计算单分子速率常数时都是相关的。然而，类似的论点也适用于过渡结构：如果可以通过不同的过渡结构达到一个产物或一组产物，那么每一个过渡结构都会对单分子速率常数有所贡献。因此，我们引入了AWATAR速率常数（方程式（8），通过将所有导致一组产物的过渡结构y的态密度之和Ny(E ? E0,y)加上所有能量可及势能井i的态密度之和ρi(E)来计算速率常数kAWATAR(E)。一组产物例如是在光解实验中观察到的碎片离子，以及一个中性碎片，两者都可能由多种异构体表示。

AWATAR对于具有能量低级异构体的簇特别有利。它已被成功用于解释(NaCl)4Na+的“魔数”行为，以及在含有小水簇的BIRD实验中的主方程建模中。47,48,58

AWATAR速率常数的效果通过(NaCl)4Na+的解离示例3.rrk来说明。图5(a)显示了具有最低能级异构体的势能面。由于异构化过程很快（见补充材料中的图S1），这个系统应该使用AWATAR方法来获得更准确的结果。文件example3.rrk相当复杂，因为它包含了我们在量子化学计算中发现的(NaCl)4Na+的所有9种异构体，作为势能井i，见图5(b)和原始工作的详细信息。实验上，我们只观察到盐簇的化学计量解离，即(NaCl)n单位的损失，其中n = 1–4（反应式（9）。

(NaCl)4Na+ → (NaCl)4?nNa+ + (NaCl)n?n = 1–4

（9）

不仅是反应物簇以不同的异构体形式存在，碎片也可能以任何能量可及的异构体形式形成。图5(c)显示了与Na+一起描述反应（9）n = 4通道的(NaCl)4的异构体。每种碎片异构体的组合都对松散过渡结构模型中的态总和有贡献，因此必须考虑所有化学计量中性及带电碎片异构体的组合。

AWATAR特别适用于具有能量低级异构体的簇。它已被成功用于解释(NaCl)4Na+的“魔数”行为，以及在含有小水簇的BIRD实验中的主方程建模。47,48,58

AWATAR速率常数的效果通过(NaCl)4Na+的解离示例3.rrk来说明。势能面图5(a)显示了具有最低能级异构体的情况。由于异构化过程很快（见补充材料中的图S1），这个系统应该使用AWATAR方法来获得更准确的结果。文件example3.rrk相当复杂，因为它包含了我们在量子化学计算中发现的(NaCl)4Na+的所有9种异构体作为势能井i，见图5(b)和原始工作的详细信息。实验上，我们只观察到了盐簇的化学计量解离，即(NaCl)n单位的损失，其中n = 1–4（反应式（9）。

(a) (NaCl)4Na+的势能面(PES)，包括在CCSD//B3LYP-D3/aug-cc-pVDZ理论水平下计算的能量低级异构体C4a–C4d之间的重排和NaCl或Na+的损失过渡状态。(b) 通过量子化学计算得出的(NaCl)4Na+的异构体及其相对能量，计算在CCSD/DZ//B3LYP-D3理论水平上。标有*的能量是在CCSD/DZ//ωB97XD/aug-cc-pVDZ理论水平上计算的，因为这些异构体是B3LYP-D3/aug-cc-pVDZ水平上的鞍点。(c) 在AWATAR中用于描述Na+损失的松散过渡结构的中性碎片(NaCl)4的异构体。(d) 使用全局最小值对所有碎片通过松散过渡结构进行解离的RRKM速率常数kRRKM。(e) 根据反应式(9)对(NaCl)4Na+解离的AWATAR速率常数kAWATAR。图(a)–(c)改编自参考文献47，并在CC-BY许可下发布。? 2024 作者。在实验中，观察到了所有(NaCl)4–nNa+，n = 1–4的片段簇大小。因此，我们对不同片段簇大小的解离单分子速率常数感兴趣，并定义了四个AWATAR反应通道，用于从所有势能井解离到所有片段异构体组合。图5(e)显示了使用单势能井松散过渡模型获得的AWATAR解离速率常数，图5(d)显示了用于比较的单一势能井松散过渡模型得到的速率常数。AWATAR速率常数对于不同片段大小的曲线与图5(d)的质量不同：片段(NaCl)3Na+和(NaCl)4具有相对低能级的振动模式，产生了较高的态总和。因此，当n = 1时，AWATAR速率常数急剧增加至约3.3 eV；当n = 4时，速率常数的增加较为温和，约为3.9 eV。这表明所有能量可及的产物碎片异构体，即使它们构成不同的松散过渡结构，也都会对速率常数有所贡献，强调了“所有势能井和过渡结构都是相关的”(AWATAR)这一点。在这个例子中，AWATAR的速率常数解释了(NaCl)4Na+的“魔法行为”：Na+的丢失通道具有相对适度的障碍，但解离的速率常数仍然很小，因为这个通道的松弛过渡态的状态总和很小。

3.3. 使用单分子速率理论描述离子-分子反应

RRKM和AWATAR计算在分析低压下的气相离子-分子反应时非常有帮助。其原理是离子-分子反应A+ + B从形成碰撞复合体(AB)+*开始，反应(10)，其速率为kc。碰撞复合体(AB)+*的内能由A+ ? B的结合能以及反应物在碰撞前的内能和相对动能组成。(AB)+*要么通过kb发生逆解离，反应(11)，要么通过一个或多个过渡结构反应生成产物，其速率为kp，反应(12)。可以通过RRKM或AWATAR计算来计算描述碰撞复合体命运的单分子速率常数kb和kp。结合对碰撞速率kc的现实模型以及碰撞复合体以能量E形成的概率p(E)，就可以从第一原理预测离子-分子反应速率系数kabs，公式(13)，前提是碰撞复合体确实以碰撞模型预测的速率形成，并且动态效应不起重要作用。此外，Aristov和Armentrout等人讨论的角动量守恒效应并未包含在这种处理中。

3.3. 使用单分子速率理论描述离子-分子反应

RRKM和AWATAR计算在分析低压下的气相离子-分子反应时非常有用。其原理是离子-分子反应A+ + B从形成碰撞复合体(AB)+*开始，反应速率为kc。碰撞复合体(AB)+*的内能包括A+ ? B的结合能以及反应物在碰撞前的内能和相对动能。(AB)+*要么通过kb发生逆解离，反应(11)，要么通过一个或多个过渡结构反应生成产物，反应速率为kp。RRKM或AWATAR计算可以用来计算描述碰撞复合体命运的单分子速率常数kb和kp。结合对碰撞速率kc的现实模型以及碰撞复合体以能量E形成的概率p(E)，就可以从第一原理预测离子-分子反应速率系数kabs，公式(13)，前提是碰撞复合体确实以碰撞模型预测的速率形成，并且动态效应不起重要作用。此外，Aristov和Armentrout等人讨论的角动量守恒效应并未包含在这种处理中。

AWATAR中单分子速率常数kb(E)和kp(E)的计算是通过OCH+到N2O的质子转移反应为例进行演示的，反应(14)：HCO+ + N2O → CO + N2OH+

在example4.rrk文件中，所有能量较高的异构体都被忽略了，这是基于图6(a)中显示的势能面(PES)。反应的起点HCO+ + N2O由松弛的过渡结构TS 0表示，碰撞复合体由Well 0表示，产物CO + N2OH+由TS 1表示。因此，定义了两个反应通道：通道0代表从Well 0到TS 0的路径，即逆解离；通道1代表从Well 0到TS 1的路径，即产物形成。这些通道分别用于计算所有能量的kb和kp。注意，能量刻度的零点被设置为Well 0的能量。然而，如前所述，系统中的可用能量总是更高的，因为它包括了复合体的形成能、反应物的内能以及它们在形成碰撞复合体之前的相对动能。在这个计算中，没有考虑辐射关联。

在AWATAR中计算单分子速率常数kb(E)和kp(E)是以OCH+到N2O的质子转移反应为例进行的，反应式为HCO+ + N2O → CO + N2OH+。

在example4.rrk文件中，所有能量较高的异构体都被忽略了，这是基于图6(a)中的势能面(PES)。反应的起点HCO+ + N2O由松弛的过渡结构TS 0表示，碰撞复合体由Well 0表示，产物CO + N2OH+由TS 1表示。因此，定义了两个反应通道：通道0代表从Well 0到TS 0的路径，即逆解离；通道1代表从Well 0到TS 1的路径，即产物形成。这些通道分别用于计算所有能量的kb和kp。注意，能量刻度的零点被设置为Well 0的能量。然而，如前所述，系统中的可用能量总是更高的，因为它包括了复合体的形成能、反应物的内能以及它们在形成碰撞复合体之前的相对动能。在这个计算中，没有考虑辐射关联。

在AWATAR中计算单分子速率常数kb(E)和kp(E)是使用OCH+到N2O的质子转移反应作为示例的，反应式为HCO+ + N2O → CO + N2OH+。

在example4.rrk文件中，忽略了所有能量较高的异构体，这是根据图6(a)中的势能面(PES)进行的。反应的起点HCO+ + N2O由松弛的过渡结构TS 0表示，碰撞复合体由Well 0表示，产物CO + N2OH+由TS 1表示。因此，定义了两个反应通道：通道0代表从Well 0到TS 0的路径，即逆解离；通道1代表从Well 0到TS 1的路径，即产物形成。这些通道分别用于计算所有能量的kb和kp。请注意，能量刻度的零点被设置为Well 0的能量。然而，如前所述，系统中的可用能量总是更高的，因为它包括了复合体的形成能、反应物的内能以及它们在形成碰撞复合体之前的相对动能。在这个计算中，没有考虑辐射关联。

在AWATAR中计算单分子速率常数kb(E)和kp(E)是通过OCH+到N2O的质子转移反应作为示例的，反应式为HCO+ + N2O → CO + N2OH+。

在example4.rrk文件中，所有能量较高的异构体都被忽略了，这是基于图6(a)中的势能面(PES)进行的。反应的起点HCO+ + N2O由松弛的过渡结构TS 0表示，碰撞复合体由Well 0表示，产物CO + N2OH+由TS 1表示。因此，定义了两个反应通道：通道0代表从Well 0到TS 0的路径，即逆解离；通道1代表从Well 0到TS 1的路径，即产物形成。这些通道分别用于计算所有能量的kb和kp。请注意，能量刻度的零点被设置为Well 0的能量。然而，如前所述，系统中的可用能量总是更高的，因为它包括了复合体的形成能、反应物的内能以及它们在形成碰撞复合体之前的相对动能。在这个计算中，没有考虑辐射关联。

在AWATAR中计算单分子速率常数kb(E)和kp(E)是通过OCH+到N2O的质子转移反应作为示例的，反应式为HCO+ + N2O → CO + N2OH+。1秒后（如图8(b)中粉红色所示），粒子分布已经接近300 K的分布（用虚线表示），而10秒后，粒子分布完全达到平衡（用绿色表示）。文件example8.rrk模拟了一个最初处于冷却状态的粒子群的热膨胀过程，这些粒子群例如是在液氮冷却的离子源中产生的，初始温度为100 K。结果展示在图8(c)和(d)中。正如预期的那样，平均能量在几秒钟内上升，10秒后粒子分布完全平衡。然而，1秒后的粒子分布（用粉红色表示）并不类似于玻尔兹曼分布，因为它有一些局部最大值，最明显的是大约240 cm?1和980 cm?1处。为了理解这些特征，详细了解粒子分布的变化方式是有启发性的：在t = 0时，粒子分布是玻尔兹曼分布，如图8(d)中的蓝色所示。每个数据点代表一个能量区间。为了计算下一个时间步长内某个区间内粒子数的变化，需要知道总体的损失和增加量。这些损失是通过所有可能的跃迁到其他区间的吸收或发射率来计算的。这是利用量子化学频率计算提供的信息，并在AWATAR程序的井面板中总结的：振动模式的能量决定了哪些跃迁是可能的。这些跃迁的概率由爱因斯坦系数描述（用于自发发射），以及普朗克辐射定律描述（用于黑体辐射场中的受激发射和吸收）。特定振动模式的爱因斯坦系数是频率和量子化学计算得到的红外强度的函数。现在让我们回到问题，图8(d)中粉红色粒子分布的特征来源是什么：仔细观察“井0”面板中的振动模式可以看到，对于大多数振动模式，红外强度很小，但在680 cm?1到820 cm?1的范围内，某些振动模式的强度相对较高。这意味着，原始粒子群将优先被激发这些振动模式的光子所移动，移动量分别为684 cm?1、707 cm?1、773 cm?1或816 cm?1，从而形成了第二个峰值。高能量侧的肩部部分是由于吸收了多个光子，部分是由于吸收了能量更高的光子。后者可以通过将所有高能量振动模式的强度设为零，然后重新运行模拟来进行比较（结果未显示）。在图8(c)中，没有进行这样的修改，我们看到内部能量在6秒内达到平衡，10秒后粒子分布再次达到300 K的玻尔兹曼分布，如图8(d)中的绿色所示。文件example9.rrk描述了在吸收了一个hν = 3091.43 cm?1光子后的粒子群演变，这个能量对应于C–H对称伸缩振动。结果展示在图8(e)和(f)中。在这种情况下，粒子群在t = 0时被初始化为300 K的玻尔兹曼分布，然后由于吸收的光子能量而转移到更高的能量（如图8(f)中的蓝色所示）。如图8(e)所示，1秒后，一部分光子能量仍然存在于分子中，而大约3秒后，粒子群已经非常接近300 K。最后，文件example10.rrk模拟了恒定光子通量的影响，例如在IRMPD实验中由红外激光提供。在这个例子中，光子能量对应于COO反对称伸缩模式的频率（1694.44 cm?1）。光子吸收速率被设置为一个固定值。如果需要，可以为每个井设置单独的速率（见用户手册）。这里我们只有一个井，速率被设置为1 s?1。结果展示在图8(g)和(h)中。由于额外的光子通量，平均能量最初会增加。2秒内，达到了一个新的稳态，平均能量更高，因为发射率也增加了。10秒后的粒子分布（绿色）与1秒后的粒子分布（粉红色）没有显著差异。正如人们可能预期的那样，更高的光子吸收速率会导致更高的平均能量，如图S2中的SI所示。

4.2. 在黑体辐射场中具有解离现象的IRMPD实时分子动力学（MEM）

对于在非可忽略的占有量的内能下会解离的分子，情况会变得稍微复杂一些。除了光子的吸收和发射外，MEM还考虑了所有高于最低过渡结构能量的区间的解离。因此，现在需要定义反应通道和过渡结构，以便AWATAR能够计算相关的速率常数k(E)。解离会导致部分粒子群的损失，所以随着模拟时间的进行，母分子或团的强度会减少。由于解离也可以通过实验测量，例如使用质谱仪，因此可以将模拟的动力学与实验结果进行比较。通过将模拟的动力学与在不同温度下测量的实验动力学进行拟合，可以确定未知的AWATAR参数，如在IRMPD实验中的光子吸收速率。在这一章中，实时MEM用于模拟Ar2FeH+的IRMPD动力学，这些动力学已经通过参考文献65进行了实验测量和发表。该分子如图9(d)所示，只有一个稳定的异构体，并通过失去一个氩原子而衰变。失去一个氩原子所需的能量被计算为2423 cm?1.65。

图9展示了Ar2FeH+的MEM，重现了在80 K（a）、(c)、(e）和300 K（b）、(d）、(f）条件下的实验IRMPD动力学。在这两种情况下，初始团簇温度都设置为300 K，反映了离子源的温度。(a)和(b)显示了实验动力学（圆圈）和模拟（线条）；(c)和(d)显示了平均内能的演变；(e)和(f)显示了0秒时的粒子分布（无解离的玻尔兹曼分布），以及1秒和10秒模拟时间后的分布。文件example11.rrk模拟了实验场景，使用了一个光子能量为1855 cm?1的光子，这对应于Fe–H伸缩振动，以及80 K和300 K的环境温度。红外强度的缩放因子被设置为4.797，这样可以重现参考文献65中发表的实验BIRD动力学。通过将MEM模拟结果与300 K场景下的实验数据进行拟合，确定了光子吸收速率为0.34 s?1。团簇粒子群初始化为300 K，即离子源的温度，适用于两种环境温度。所得的动力学（线条）与实验动力学（圆圈）一起展示在图9(a)和(b)中。拟合几乎完美地再现了两种动力学，使用了一组拟合参数。为了更可靠地确定参数，需要在其他温度下进行额外的实验动力学测量，以验证所选参数集在更宽范围内的有效性。图9(c)和(d)显示了每次模拟的平均能量。最初，平均能量是相同的，因为两者都是从300 K的初始团簇温度开始的，这个温度是由激光蒸发源的温度决定的。在80 K的情景中，辐射冷却最初占主导地位，超过了红外激光的影响。在300 K的环境中，平均能量如IRMPD实验中所预期的那样增加，其中初始团簇温度与环境温度相匹配。图(e)和(f)显示了三个选定时间点的粒子分布。在80 K时，平均能量的减少也反映在粒子分布中。在300 K时，10秒模拟时间内的粒子分布非常相似。值得注意的是，在这两种情况下，玻尔兹曼分布的高能量尾部在过渡结构能量（2423 cm?1）以上都有所减少。这种快速下降反映了由于解离而导致的粒子损失。t = 0秒时的粒子分布是模拟开始时的玻尔兹曼分布。这个起始粒子分布还没有包括解离，因此在过渡结构能量以上没有减少。请注意，粒子分布并没有反映由于解离而导致的总体信号损失，它们都被归一化为100%。

4.3. 涉及多种异构体的BIRD实时MEM

AWATAR可以用来对具有多种异构体的系统进行MEM，前提是异构化发生的速度远快于解离。当包含多种异构体时，每个能量区间内的粒子分布会在每个时间步长后根据它们的状态密度在所有异构体之间重新分配。如果涉及多个过渡结构，解离可以通过所有能量上可访问的过渡结构发生。不同的通道定义允许在单个通道（AWATAR方法）中包含所有井和所有过渡结构，或者定义多个通道，每个通道将单个井与单个过渡结构关联起来。例如，在参考文献48中使用的，也可以通过为每个产物定义单独的通道来模拟解离成多种产物。在这种情况下，每个通道将所有井与导致相应产物的过渡结构连接起来（详见手册）。在AWATAR中进行包含多个最小值的MEM是直接的：除了加载所有井和过渡结构外，它的操作方式与单个井和单个过渡结构相同。文件example12.rrk设置为模拟在288 K环境温度下的双水合碳酸根自由基阴离子CO3(H2O)2?的BIRD动力学。图中包括了四个井和两个 loosen的TS，如图10(a)所示，同时提供了相对能量和简并度。简并度是通过将光学异构体的数量与自旋多重性相乘来确定的（详见用户手册）。在这种情况下，简并度为2，是由于井1的光学异构体数量为2。最低过渡结构的活化能被设置为3762 cm?1，这是根据参考文献64确定的。红外强度的缩放也来自于参考文献64中的一个参数。定义了一个AWATAR通道，包括了所有井和所有过渡结构。example12.rrk得到的动力学结果显示为线条，与参考文献64中报告的相应实验动力学数据（圆圈）一起展示在图10(b)中。值得注意的是，通过对活化能的拟合和红外强度的缩放，以及在不同温度下的多个数据集的Arrhenius图，模拟结果与实验结果非常吻合，详见参考文献64。图10展示了CO3(H2O)2?的BIRD动力学MEM模拟。(a)展示了带有相对能量（单位cm?1）和简并度（括号内）的井和过渡结构的优化结构；(b)将模拟数据（线条）与实验动力学数据（圆圈）进行比较；(c)模拟了不同初始团簇温度下的动力学，导致了诱导延迟。对于一个适当热化的粒子群，即当初始粒子群与环境温度下的稳态粒子群相匹配时，反应物的衰减预期是指数型的，在半对数图中表现为线性的，如图10(b)所示。这种行为确实在实验动力学中观察到，并且也被MEM模拟所再现，如图10(b)所示。然而，如果初始团簇温度低于环境温度，团簇在实验过程中会加热，导致一个称为诱导延迟的现象。最初，由于团簇中可用的能量较低，解离部分被抑制。随着时间的推移，当粒子群接近稳态时，预期的线性行为就会出现。这样的诱导延迟可以通过在参数窗口中设置粒子群的初始温度在AWATAR中模拟出来。图10(c)展示了初始温度为200 K和100 K时的模拟诱导延迟，分别用绿色和紫色表示。诱导延迟导致开始时的解离减少。然而，几秒钟后，所有三种情况下反应物衰减的梯度变得相同，这表明如果从半对数图的动力学线性部分中提取的话，解离速率常数最终与初始簇的温度无关。使用ATAR进行的MEM（分子统计力学）模拟可以帮助验证实验中是否真的满足了这些条件，或者可以指导实验者据此来规划实验。

4.4. 稳态MEM

除了模拟实时MEM之外，ATAR还可以模拟在环境条件保持不变的情况下，种群在一段时间后达到的稳态。稳态意味着种群的能量分布以及在有限时间步长Δt内形成的产物量是收敛的。文件example13.rrk提供了在两个环境温度288 K和250 K下进行的CO3(H2O)2?的稳态MEM计算。所有参数、井和TS（状态空间）都与前一章讨论的example12.rrk相同。稳态MEM参数设置得当，以使计算在相对较短的时间内完成（详见用户手册）。在这种模式下，有限的时间步长和收敛标准最初是粗略的，但在计算过程中会通过迭代进行细化。种群数据定期写入文件example13_pop.txt中，并且计算的状态会定期存储在检查点文件中。一旦种群在最小的时间步长内收敛，模拟就结束了。结果面板包含了所有温度和所有定义通道的MEM速率系数。这些速率系数可以绘制在Arrhenius图中，与实验数据的比较可以确定活化能和红外强度的缩放，具体描述参见参考文献46和48。图11显示了两种计算温度下的种群结果。与Ar2FeH+的情况一样，种群在最低过渡结构的能量以上急剧减少。总种群用黑色表示，而彩色线条表示各个井的贡献。有趣的是，全局最小值（橙色）并不是人口最多的井（绿色）。这是因为井2在E > 1000 cm?1时具有更高的态密度。因此，由于在模拟中包含了所有井，总种群向更高能量方向移动，平均解离速率比单井模型更高，单井模型用虚线表示以便对比。因此，在MEM中包含所有井可以改进活化能的测定，因为与单井模拟相比，它提供了更真实的实验表现。

在环境温度(a) T = 250 K和(b) T = 288 K下，使用ATAR方法和所有已知的井计算得到的CO3(H2O)?的收敛种群。请注意，人口最多的井（绿色）并非全局最小值（橙色）。虚线是单井方法下的种群。

5. 总结与展望

我们希望这里提供的示例能够说明单分子统计速率理论对从事气相反应研究的科学家们所能提供的附加价值。这里展示的示例来自典型的质谱设置。像TRAPed Ions的IRMPD这样的现代技术可以从实时的过程中获得定量模拟中受益。显然，了解特定分子或簇离子解离的时间尺度是非常有价值的。如果系统有量子化学计算，ATAR1.0可以相对轻松地提供这些信息。可以设想对该代码进行许多有益的扩展，例如包括非谐振子的量子态模型、除了辐射过程之外的碰撞能量交换，或在反应性碰撞中适当热化离子和中性粒子等。代码的任何更新都将参考本教程综述，并将在Github和Zenodo上发布。我们希望ATAR1.0有助于提高我们对分子尺度上辐射过程和反应过程相互作用的理解。

作者贡献

Magdalena Salzburger：概念化（支持），数据整理（主要负责），形式分析（主要负责），研究（主要负责），方法论（支持），软件（共同负责），验证（共同负责），可视化（主要负责），撰写——初稿（主要负责），撰写——审阅和编辑（共同负责）。Marc Reimann：形式分析（支持），资金获取（主要负责），研究（支持），软件（支持），撰写——审阅和编辑（主要负责）。Jessica C. Hartmann：形式分析（支持），研究（支持），可视化（支持），撰写——审阅和编辑（支持）。Milan On?ák：概念化（支持），资金获取（主要负责），研究（支持），方法论（共同负责），软件（支持），监督（共同负责），验证（共同负责），撰写——审阅和编辑（支持）。Martin K. Beyer：概念化（主要负责），数据整理（共同负责），形式分析（支持），资金获取（主要负责），研究（共同负责），方法论（主要负责），软件（主要负责），监督（主要负责），验证（共同负责），可视化（支持），撰写——初稿（主要负责），撰写——审阅和编辑（共同负责）。

利益冲突

没有需要声明的利益冲突。

数据可用性

示例文件可作为补充信息（SI）提供。更多的背景信息和注释也在SI中提供。补充信息包括MEM计算的详细描述；样本输入文件；关于BS-Quantum参数的讨论；对输入文件的评论。详见DOI: https://doi.org/10.1039/d6cp00705h。ATAR 1.0的源代码和可执行文件，连同快速入门指南、用户手册和样本输入文件，发布在Github和Zenodo上，网址为https://doi.org/10.5281/zenodo.19915221。

致谢

非常感谢德国研究基金会（Deutsche Forschungsgemeinschaft）的财政支持，授予编号为BE2505/3-1的资助（MKB）。本研究还得到了奥地利科学基金（Austrian Science Fund）的部分资助，授予编号为10.55776/W1259（MS, MKB）、10.55776/P35407（JCH, MKB）、10.55776/PAT3834224（MS, MKB）、10.55776/ESP662（MR）和10.55776/P35013（MO）的资助。所呈现的计算结果使用了因斯布鲁克大学（University of Innsbruck）的HPC基础设施LEO。作者感谢Michael Hütter为ATAR设置Github仓库，以及Gabriel Sch?pfer在讨论中提供的帮助。