构建任意阶内部坐标变换以改进对大振幅运动的研究

《The Journal of Physical Chemistry A》：Construction of Arbitrary-Order Internal Coordinate Transformations to Improve Studies of Large-Amplitude Motions

　　高分辨率图像
下载MS PowerPoint幻灯片

内部坐标及其导数是高效处理几何优化、高分辨率光谱模拟以及量子化学中势能面拟合的基础。现有的关于内部坐标导数构建的描述通常要么缺乏简洁性，要么缺乏通用性。在本文中，我们提供了一个简单的框架，用于评估任何阶数的内部坐标导数，并提供了一种自动获取相应逆变换的方法。通过进一步扩展到内部坐标系统之间的变换，这种方法为处理各种分子问题提供了一种完整且通用的方法。该构建的实用性通过研究内部坐标插值在异构化中的行为、量化有机金属系统中羰基伸缩与复杂伸缩坐标之间的耦合，以及分析机器学习得到的原子间势能在计算非谐频率时的性能得到了证明。这种方法在数值上高效且通用，并提供了参考实现。

特刊
作为《物理化学A杂志》的“John F. Stanton纪念特刊”的一部分发表。

引言
在量子化学中使用内部坐标与该领域本身一样悠久。虽然分子哈密顿量最初是在外部笛卡尔参考框架中构建的，但很快就描述了将其转换为固定体系内部坐标的适当变换，并在接下来的几十年中发展出了处理振动-旋转问题的方法。Wilson推广了构建和使用内部坐标的一般方法，随后针对不同类别的问题开发了其他形式主义。内部坐标为描述化学系统的基本物理提供了更直观的表示，并已被广泛用于理解复杂的化学和光谱问题。内部坐标方法也被长期用于开发更高效的势能面表示。最近，机器学习方法被应用于给定分子系统的势能面拟合，具有一定的可转移性。这些模型可能比密度泛函理论快数千倍，并通过自动微分提供解析导数。已经对这些模型的光谱适用性进行了研究，特别是使用这些方法计算的零点能量的准确性。作为系统特定机器学习势能面的替代方案，是拟合通用的机器学习原子间势能（MLIPs）。这些模型已被证明适用于能量研究等其他应用。对MLIP能量及其导数在动力学和光谱学中的性能的表征是一个开放的研究领域。内部坐标方法非常适合此类研究，可以访问模间耦合，冻结处理不佳的运动，并允许沿着化学相关的自由度进行高效采样。

在内部坐标系统中工作时，通常至少需要获得内部坐标相对于感兴趣的系统的相应笛卡尔坐标的一阶导数，因为这些导数用于构建分子动能算符，这是通过振动微扰理论或基组方法获得振动波函数所必需的。逆向导数，即系统的笛卡尔坐标相对于内部坐标的导数，用于将势能、偶极子和极化率表面转换为内部坐标表示，从而可以评估跃迁强度或进行更一般的势能面分析，其相关性远远超出了振动光谱或自由能计算的领域。关于最常用的内部坐标（例如键长、角度和二面角）相对于相应笛卡尔坐标的导数的表达式在文献中随处可见，更高阶的导数也在某些处理中出现。同样，数值微分也非常常用，以获得内部坐标相对于相应笛卡尔坐标的导数以及笛卡尔坐标相对于内部坐标的导数。然而，随着现代处理能力的提高以及用于线性和多线性代数运算的高效库的出现，拥有一个简单、高效且通用的方法来评估最常见的内部坐标的导数及其逆变换是非常有用的。

在这项工作中，我们详细介绍了这样的通用方法，并提供了具体的例子，说明如何使用内部坐标处理来理解复杂的化学现象。具体来说，这里描述的变换允许在内部坐标空间中进行改进的插值，应用于乙炔到乙烯基的异构化过程、环戊二烯基有机金属羰基复合物中低频复杂伸缩的处理及其与羰基伸缩的耦合，以及通过内部坐标振动微扰理论探究机器学习相互作用势的行为。

理论
为了使这项工作自包含，我们首先描述了坐标变换的一般理论，以及如何在坐标系统之间转换任何感兴趣的属性的导数的简单方法。从那里，我们用张量/多线性代数表达式为化学问题中最重要的坐标构建了简单的表达式，并描述了如何使用坐标变换来获得任何感兴趣的坐标的表达式。接下来，我们展示了如何隐式获得逆变换，特别注意如何仅使用内部坐标相对于笛卡尔结构的导数来获得嵌入在Eckart框架中的笛卡尔坐标的导数。最后，我们展示了如何在这个框架中表示更复杂的坐标和坐标系统之间的变换，同时保持可微性。

坐标变换理论
在构建内部坐标的哈密顿量表示或展开时，通常需要能够以系统化的方式评估属性导数在坐标系统之间的变换。通常将这些变换明确地写成可以在编程语言中评估的表达式；然而，也可以通用地表示这样的变换。符号的具体内容和相应的推导在S2节中给出，但为了方便起见，我们将引入符号来定义通过链式法则在坐标系统R和X之间变换标量量V的方法：
?RV = [R(1)ixV(1)x] = ?RX?XV
?RV = [R(2)ijxV(1)x] + [R(1)ixR(1)jyV(2)xy] = ?(2)RX?XV + ?RX(?(2)XV)?RXT(1)
其中?A(k)B表示B相对于坐标系统A的k阶导数张量。

使用这种符号，我们可以以紧凑的形式写出高达三阶的导数变换表达式：
?(3)RV = [R(3)ijkxV(1)x] + [R(2)ijxR(1)kyV(2)xy] + [R(1)ixR(1)jyR(1)kzV(3)xyz]
其中下标表达式(ijk + ikj + kij)表示确保对称性的相关置换集合，如S2节所述。应该注意的是，尽管这些表达式是逐元素书写的，但它们是使用一系列张量收缩和置换来评估的。

很容易看出这种构建的好处。加法比乘法更快，因此在没有硬件考虑的情况下，这种表达式用两次置换和两次加法替换了四次收缩，相对于简单的递归。此外，如后面将讨论的，这种方法通常比相应的数值微分快得多。通过考虑每个括号中的指数，我们可以进一步了解这个问题的结构。我们可以注意到，在每个括号内，R(k)项的指数之和等于导数的阶数，R(k)项的数量给出了相应V(n)的指数。最后，执行的置换是最小化的集合，以确保项保持完全对称，因为对于物理上相关的属性来说，导数运算的阶数并不重要。

因此，我们可以更一般地写出这个过程：
?(m)RV = ∑p∈P(m)[R(p1)...x1...R(pj)...xjV(j)x1...xj](T(p))
其中P(m)是m的整数分割集合，j是给定分割p的长度，T(p)是用于对称化的置换集合。我们在上面的表达式中保留了未收缩的索引。确定对称化置换是唯一排列理论的直接应用，在应用接触变换后广泛用于确定状态空间。这项工作的形式推导在S3节中给出。

坐标定义和展开
对于本研究中的问题，我们在键角-二面角坐标中展开量。对于给定的笛卡尔坐标X，这些量有以下标准定义：
rij = |xi?xj|
θijk = arctan2(sin(xi?xj, xk?xj), cos(xi?xj, xk?xj)
τijkl = arctan2(sin(nijk, njkl), cos(nijk, njkl)
其中
sin(a,b) = |?a×?b|
cos(a,b) = ?a·?b
nijk = xj?xi×xk?xj

应该注意的是，使用双参数反正切函数可以在角度接近0和π时提高数值稳定性，并且不需要额外的逻辑来获得正确符号的二面角。这些坐标在图1中显示。

图1. 本研究中考虑的不同类型的坐标。

首先，我们将这些量相对于相应原子坐标的分量进行微分。我们首先列出一系列来自标准多元微积分的属性：
?a|a| = a|a|
(a×b)γ = [ai?iγjbj]
?(k)Xarctan2(sinθ, cosθ)k = ?(k)Xsinθ??(k)Xcosθ
其中?是方向余弦或3D Levi-Cevita张量。这些表达式可以使用计算机代数系统验证，并且sin2θ + cos2θ = 1，与θ的值无关。

有了这些属性，我们现在可以使用链式法则将键长展开到一阶：
?Xrij = ?X(xi?xj)·(xi?xj|xi?xj|)
这本身并不值得注意，但它展示了我们将使用的一般结构。所有内部坐标量都将相对于最方便的坐标系统进行导出，然后使用上述重新展开程序来评估整个张量。

一个更值得注意的例子是键角的一阶导数，我们将其写为：
?Xθijk = ?X(sinθijk, cosθijk)·(cosθijk, sinθijk)
?Xcosθijk = ?X(a,b)·(b|b|, a|a|)
?Xsinθijk = ?X(a,b)·([ai?ixj][ai?ixjbj], [?ixjbj][ai?ixjbj])
a = xi?xjb = xk?xj
其中(a, b)表示向量a和b的连接。应该注意的是，?X(a, b)是稀疏的，对于高维系统应该使用稀疏矩阵方法来评估该张量积。

虽然这些导数的显式公式可能在低阶时提供几何洞察，但这种方法的强大之处在于它可以直接评估任意阶数，只要能够评估?a(k)|a|对于任何k。如方程7所示，向量范数的导数依赖于该距离的标量倒数，因此我们得到了一系列增长的导数：
?a|a| = a|a|
?(2)a|a| = I|a|?a?a|a|
?(3)a|a| = ?[I?a]ijk+ikj+kij|a|
?(3)a|a| = ?[I?a]ijk+ikj+kij|a|
我们可以看到，这与方程3中的表达式类别相同，只是用张量积代替了张量收缩，并且我们的导数展开限于二阶。再次使用了对称化置换集合T(p)。

有了这些，我们就有了构建上述任何形式内部坐标导数所需的所有组件，需要注意的是对于平面二面角，我们将arctan中的表达式替换为法线之间角度的正弦。这种替换的有效性可以通过数值微分来确认。

扩展这套基本集，通常使用键长、角度和二面角的一些潜在非线性组合作为大振幅运动的坐标。为了评估这些的导数，只需首先评估相对于初始键长、角度和二面角的导数，然后使用上述规则将这些导数转换到最终坐标系统。只要从初始内部坐标系统到最终坐标系统的变换是可微的，相对于笛卡尔坐标的导数就很容易评估。这些类变换的应用将在后面的章节中讨论。

逆坐标变换
我们描述了一种计算内部坐标相对于原子坐标的导数的程序。这就是评估分子动能算符所需要的一切。然而，为了用内部坐标重新表示像势能这样的标量量，也需要进行逆变换。为此，我们注意到给定坐标系Q和R之间存在一个可逆的变换，我们有?RQ?QR=I（13），因此对于任何k > 1，我们有?(k)R(?RQ?QR)=0（14）。设?Q(k)R = Q(k)和?R(k)Q = R(k)，通过展开上述关系并分离出R的最高阶导数项，我们得到一系列方程[R(2)ijxQ(1)xq]=?[R(1)ixR(1)jyQ(2)xyq][R(3)ijkxQ(1)xq]=?([R(2)ijxR(1)kyQ(2)xyq](ijk+ikj+kij)+[R(1)ixR(1)jyR(1)kzQ(3)xyzq])（15）。现在考虑沿着最后一个轴收缩，有(Q(1))?1 = R(1)，我们得到R(2)ijq=?[R(1)ixR(1)jyR(1)zqQ(2)xyz]R(3)ijkq=?([R(2)ijxR(1)kyR(1)zqQ(2)xyz](ijk+ikj+kij)+[R(1)ixR(1)jyR(1)kzR(1)aqQ(3)xyza])（16）。因此，给定Qk的展开，我们可以通过注意到逆变换展开中的高阶项必须消失来隐式获得逆变换的导数。这个过程与上面描述的重新展开过程直接类似，最简单的实现使用相同的程序。虽然笛卡尔坐标和内部坐标之间的变换是不可逆的，但在Eckart框架中的笛卡尔坐标之间的变换是可逆的（20,62）。因此，通过消除?X(k)R沿旋转或平移生成元的任何分量（63），可以使用相同的隐式反演过程，只是使用从?XR的伪逆得到的?RX。应该注意的是，如果提供了一组特定的笛卡尔位移及其相应的展开，也可以使用相同的过程来表示?XR，例如，系统的一个片段绕一个确定轴的旋转。

在不同分子的构象中，同一组原子可能有不同的物理直观坐标系。为了能够将属性从一组坐标平滑地转换到另一组坐标，需要能够在任何给定的内部坐标系之间进行转换。在这项工作中，例如，通过注意到在图1中，由rij、θijk和rjk定义的边-角-边三角形可以通过rik=√r2ij+r2jk?2rijrjkcosθijk转换为相应的边-边-边三角形（17）来做到这一点。其余的角度可以用类似的方法计算。最后，图1中的距离ril可以通过给定的坐标计算得出r2il=r2ij+r2jk+r2jk?2rijrjkcosθijk?2rjkrklcosθjkl?2rijrkl（18）。值得注意的是，可以使用上述框架直接从这些表达式计算导数张量?X(n)rik。然后可以通过基于图的直接方法自动确定初始坐标系R和目标坐标系S之间的直接转换，以确定所需的补全。

计算细节
电子结构计算在四个理论水平上进行。所有耦合簇计算，即CCSD(T)优化和单点能量，以及EOM-CCSD优化和垂直激发能量，都是使用Q-Chem 6.1软件包和cc-pVTZ基组进行的（71）。密度泛函理论计算，包括优化、Hessian和部分四次力场计算，在B3LYP和ωB97X-D3理论水平上使用cc-pVTZ基组进行，以及使用Martin及其同事的具有色散校正的PBE泛函（72）和def2-TZVP基组的优化和Hessian计算，是使用Gaussian 16软件包进行的（73）。能量、梯度和Hessian也是使用两种最近开发的机器学习原子间势（MLIPs）通过自动微分计算的，AIMNet2家族的aimnet2_wb97m_0模型是针对约107次在ωB97M-D3/def2-TZVP理论水平和基组上的DFT计算进行训练的（48），以及MACE-oMOL家族的额外大型模型是针对约108次在ωB97M/def2-TZVP理论水平上的DFT计算进行训练的（49）。每个MLIP的精确模型文件都包含在本工作的补充数据仓库中（74）。在本文的其余部分，这些模型将分别用AIMNet2和MACE表示。优化是使用拟牛顿步骤进行的，辅以基于单纯形的优化，使用Nelder–Mead无梯度优化算法，因为在执行约束优化时梯度对于拟牛顿优化器来说不够平滑。部分四次力场使用每个正常模式的步长0.5 a0进行评估，并在Hessian上计算3点中心差分。选择这个步长是为了与使用两种MLIP通过自动微分计算的第三导数保持一致，以及其在振动扰动理论的其他实现中的使用（75）。还使用通用的有限差分引擎比较了其他有限差分方案（5点和7点中心差分），但确定3点中心差分已经足够（76）。非谐频率是使用Psience Python软件包中可用的基于稀疏矩阵的振动扰动理论实现来评估的（33,77,78）。所有使用MLIP的计算都是在没有GPU的8核机器上进行的。

结果与讨论
在考虑直接应用内部坐标处理有趣的问题之前，讨论本工作中提供的内部坐标展开相对于数值微分的好处是重要的。我们不会将这种构建与其他直接解析实现的内部坐标进行比较，因为这种构建与此类方法的性能差异可能强烈依赖于实现的特定细节和所使用的系统。对于数值比较，我们将考虑蒽。为了获得运行二阶振动扰动理论计算所需的势能和Wilson G矩阵的导数，需要关于参考笛卡尔坐标的内部坐标的三阶导数。数值微分使用了笛卡尔坐标和内部坐标之间的向量化转换实现，并且有限差分使用了通用有限差分库中的5点中心差分模板，该库利用了导数张量的对称性。相对性能将取决于实现，但在获得关于相应笛卡尔坐标的内部坐标的导数时，对于一个有23个原子的系统，本文提供的方法大约快65倍，并且不需要代码来检测键角或二面角坐标是否发生了大的变化。获取Eckart嵌入笛卡尔坐标关于相应内部坐标的导数的逆问题具有更大的性能差异，第2.3节中描述的逆变换比直接数值微分快大约230倍，同时不需要额外的边缘情况处理。用于这些测试的笛卡尔坐标和Z矩阵可以在支持本工作的数据仓库中找到（74）。使用专门的或编译的实现可以实现更高效的数值导数，但这些时间代表的是这些方法的直接应用。

作为一个额外的好处，使用第2.3节中的过程，只需定义内部坐标关于笛卡尔坐标的导数，就可以获得到任何一组内部坐标的转换系数。这将在第4.2节中进一步讨论，但它允许在其他应用中构建“非局域化”的内部坐标（22,26,27）。

坐标系统之间的插值
一个经典的例子是异构化过程，其中最佳的内部坐标描述在不同的最小值之间变化。也许最简单且研究最充分的异构化是从乙炔到乙烯基的转化（79-85），如图2所示。
图2. 乙炔、弯曲的乙炔和乙烯基，显示了相关的坐标，以及指示异构化路径的箭头。为了方便起见，没有显示HCCH二面角。rCC、r1和θ1对两种构型都是共同的，但为了清晰起见，只显示了一次。

在类似乙炔的构型中，包括在S1表面上最稳定的弯曲乙炔构象，自然可以将系统视为两对CCH角及其相应的距离，以及定义系统平面的二面角。在类似乙烯基的构型中，自然坐标更类似于甲醛的坐标——一个HCH角、相应的CH距离、一个CC距离以及一个CCH角和HCCH二面角。这两个CCH角可以代替HCH角，但在两种情况下，关键的区别是一个CH距离被另一个取代了，因为氢从类似乙炔的状态转变为类似乙烯基的状态。

有多种可能的插值方法可以考虑，包括使用类似乙炔的坐标在最小值之间直接插值，对于类似乙烯基的坐标也是如此，第2.4节中讨论的两种插值形式，以及一种“直接”插值，它插值五组可能的CC和CH键长以及二面角以强制平面性。通过这些方法插值时，系统在基态和S1势面上的能量分别在图3中给出。使用AIMNet2和MACE机器学习原子间势得到的类似图S1的图表，其中系统用每种势重新优化并生成了相应的插值。

图3. (A) 在CCSD(T)/cc-pvTZ理论水平上，使用上述指数形式（蓝色，实线）插值（紫色，实线）和类似乙烯基的坐标（红色，虚线）的几何形状比较基态能量。(B) 在EOM-CCSD/cc-pvTZ理论水平上，使用类似乙炔的坐标（蓝色，实线）插值（绿色，点线）和类似乙烯基的坐标（红色，虚线）的几何形状比较S1能量。

显然，最直接的方法是基于距离矩阵的一部分直接插值，得到的过渡在物理意义上是最不明确的。在类似乙烯基的极限情况下，使用类似乙烯基的坐标比使用类似乙炔的坐标有明显的改进，如图3所示，指数和Sigmoid平滑插值的表现都相当好。在类似乙炔的极限情况下，情况稍微复杂一些。对于非常小的畸变，类似乙炔的坐标在畸变后会导致能量降低，但在障碍最大值之前这种情况就不再成立。在图3中，可以清楚地看到Sigmoid插值遵循这种行为，因为几何形状在任一最小值时都类似于类似乙炔或类似乙烯基的插值。有趣的是，指数插值在大约20,000 cm–1的范围内始终得到最低的能量结构。这种效应在两种MLIP理论水平上都可见。

这种差异的起源可以通过考虑相关的原子间距离来探索。考虑到图2，类似乙烯基和类似乙炔的坐标系统之间有两个共同的距离rCC和r1，以及θ1在两个系统中都存在，图2中未标记的相应第三个距离是隐式定义的。因此，在类似乙炔的插值中，r2可能会变化，而r3是线性插值的，而在类似乙烯基的插值中则相反。这在图4中可以清楚地看到。值得注意的是，类似乙烯基的插值（虚线）在氢原子旋转出来时导致r3（绿色）最初减小，而r2（蓝色）线性减小。对于类似乙炔的插值，r3线性增加，而r2最初增加然后减小。指数插值能够捕捉到r2的微小增加和r3的微小减少，因为它利用了两个最小值的信息，这种更大的灵活性使得它能够利用两个最小值的优点，从而实现能量更低的插值。然而，如果类似乙炔的坐标系统对类似乙炔的最小值更为适用，这种优势就会消失，此时S形插值的表现会更好。这正是在S1表面插值中观察到的情况（图3的B面板）。图4显示了乙炔到乙烯基乙烯转化过程中的关键原子间距离，图中标记了r2和r3。这些距离是根据第2.4节描述的指数插值（实线）、乙烯基乙烯插值（虚线）和类似乙炔的插值（点线）绘制的。相应的能量可以在图3中看到。

高分辨率图像下载MS PowerPoint幻灯片：有机金属配合物的广义内部坐标和振动耦合

通常情况下，某个感兴趣的坐标无法轻易地包含在一个非冗余的内部坐标集中。为了说明这一点，我们将考虑环戊二烯基钼三羰基([CpMo(CO)3]?1)阴离子。这个复合物的中性二聚体已经被研究过（86-88），而环戊二烯基金属复合物在无机合成中具有普遍的意义（89-93）。[CpMo(CO)3]?1的重要振动模式在图S2中展示，包括环戊二烯基的扭转和复合物的伸缩，这些振动模式改变了环戊二烯基部分与MoCO3基团之间的距离，但不改变它们的相对取向。这个坐标的理想化版本如图5所示。金属羰基化合物的光谱学已经被广泛研究，这些系统中的主要强度载体是羰基伸缩（94-97）。因此，研究环戊二烯基位置如何调节羰基伸缩频率作为这些运动之间耦合的指标可能是有意义的。

高分辨率图像下载MS PowerPoint幻灯片

我们首先考虑围绕扭转坐标进行的松弛扫描。所有关于这个复合物的电子结构计算都是使用Martin及其同事（72）提出的具有色散校正的PBE泛函，以及Kulik及其同事推荐的molybdenum原子的有效核势（98-100）进行的。这些电子结构计算使用了Gaussian 16（73）软件包。几何结构经过优化，扭转坐标以18°的步长扫描了72°的范围，每一步都允许其他坐标松弛。在这个范围内，扭转坐标的频率范围是2到8 cm–1，表明环戊二烯基部分可以自由旋转，对CO伸缩频率的影响可以忽略不计。虽然这样的低频坐标特别适合于内部坐标研究，但由于缺乏明确的耦合关系，因此无法得出确定的结论。因此，我们将考虑在这个研究中定义的复合物伸缩，即环戊二烯基碳原子的质心与钼原子之间的距离。这个坐标的平衡值为2.0645 ?，我们将考虑从1.8645 ?到2.3645 ?的位移。具体来说，在松弛过程中，在环戊二烯基的质心处放置了一个虚拟原子，并且Z矩阵被约束以确保环戊二烯基质心的位置保持不变。在这个扫描过程中，名义上的复合物伸缩模式从396 cm–1变化到148 cm–1，而反对称的CO伸缩从1844 cm–1变化到1831 cm–1，对称的CO伸缩从1940 cm–1变化到1951 cm–1。这在图6中绘制出来。这意味着这些坐标之间存在中等程度的耦合。

高分辨率图像下载MS PowerPoint幻灯片

然而，名义上的复合物伸缩也包含了MoCO弯曲和CO扭转，如图S2所示，这反过来会影响耦合。为了解释这一点，我们希望构建一个纯粹的复合物伸缩坐标和沿此扫描路径上每个点的局部CO伸缩，以便观察这些坐标之间的双线性耦合如何随着质心-Mo距离的变化而演变。这可以通过添加虚拟原子和进一步的约束来实现，但本文描述的方法提供了一种更直接的方法。我们将通过两种不同的方式来处理这个问题。首先，我们将利用系统的对称性。为了构建复合物伸缩模式，我们注意到这个坐标能够对称地增加或减少环戊二烯基中碳原子与钼原子之间的距离（CMo距离），同时保持CC和CH键长不变，也不改变环戊二烯基与MoCO3部分之间的任何欧拉角。因此，给定一组作为系统Z矩阵提供的内部坐标，我们可以通过添加所有的CMo、CC和CH距离来构建一组冗余的坐标。然后，感兴趣的坐标是以下对称组合：
S=1/√5(r1+r2+r3+r4+r5)（19）
其中ri是C–Mo距离之一。接着，我们找到这五个ri的四个线性组合，这些组合与这个坐标正交，以定义一组新的内部坐标{Si}。这组新的内部坐标{Si}是初始内部坐标{Ri}的线性变换，导数张量?X(n)S可以使用上述变换来计算，注意当n > 1时，?R(n)S为零。然而，通常这组坐标会过完备，因此需要将其简化为3N – 6个坐标。这是通过两个步骤完成的：首先以标准方式构建一组非局域化内部坐标（26）。接下来，使用最小二乘变换将这组3N – 6个坐标重新定位，使非局域化坐标变换矩阵尽可能接近单位矩阵。在这种构建中，我们确保方程19中定义的伸缩坐标和CO伸缩坐标不会混合。经过这些变换后，我们得到一组3N – 6个坐标Q及其导数?X(n)Q，可以用来评估CO伸缩与复合物伸缩坐标S之间的耦合。为了方便和通用性，在实施上述程序时，通过将C5v点群的对称操作嵌入子系统的框架中进行对称化（10），但通常任何包含S以及完整的CC和CH伸缩坐标的坐标集都可以使用。

另一种更直接的方法需要对感兴趣的坐标有更具体的了解。在这种方法中，我们注意到我们想要描述的坐标涉及沿连接环戊二烯基质心和钼原子轴线的原子位移。因此，我们可以找到系统内部坐标的线性组合，这些组合在线性阶上最接近这种模式。具体来说，我们取一组初始的内部坐标作为这个系统的Z矩阵，并找到导数?RX的线性组合，这些组合与感兴趣的位移向量最大程度地一致，使用最小二乘法。这种方法最好在质量加权坐标中应用。然后我们构建3N – 7个与这个坐标正交的内部坐标的线性组合，以获得一个完整的坐标系统。这种方法类似于用原始内部坐标表示平移-旋转内部坐标（101）。

这项分析的结果在表1中提供。比较复合物伸缩频率，我们注意到特定的坐标处理对相关局部模式频率有显著的影响，直接构建的方法一致地给出了更高的频率模式，特别是当Cp-Mo距离减小时，使用对称坐标得到的局部频率为596 cm–1，而使用局部坐标得到的局部频率为867 cm–1，相对于相应的名义Cp-Mo伸缩正常模式为397 cm–1。这种差异是由坐标处理的细微差异引起的，在应用纯缩放以确保局部模式与系统的正常模式等效时可以最小化这种差异（102）。缩放后的频率比较以及相应的双线性耦合元素在表2中提供。

表1. [CpMo(CO)3]?1中局部伸缩频率与正常模式频率的比较
a 对称直接 ΔS S COa COs SCOSCO–0.2397 1844 1940 596 1869 867 1869–0.1342 1839 1940 486 1863 687 1863 0.0288 1835 1941 393 1859 541 1859 0.1240 1833 1944 320 1858 427 1858 0.2193 1833 1948 248 1858 3191 1858 0.3148 1832 1951 1961 1858 229 1858

表2. [CpMo(CO)3]?1中局部伸缩频率和双线性耦合在缩放后的结果
a 对称直接 ΔS S COw SCOw–0.2443 1756–3047 31756 31–0.1372 1749–2939 51749–290.0312 1744–2632 51744–250.1257 1741–2226 61741 220.2205 1740–1920 81740–180.3156 1739–1515 11739 15

在缩放后，对称构建的复合物伸缩与直接构建的复合物伸缩之间的差异减小，趋势与名义模式一致。缩放后的CO伸缩的局部模式频率明显低于CO伸缩的正常模式频率。这是由于CO伸缩与相应的MoCO弯曲之间的强耦合造成的，当在CO伸缩和MoCO弯曲的子空间内进行对角化时，正常模式频率得以恢复。随着Cp-Mo距离的增加，CO伸缩与复合物伸缩坐标之间的耦合以及频率都会减小。这个结果并不特别令人惊讶，但不需要额外的计算或近似，类似的分析也可以使用上述方法高效地应用于更复杂的系统。这种方法还展示了不同的坐标系统如何导致对光谱可观测量的不同解释，以及如何缓解其中的一些差异。

理解使用内部坐标扰动理论训练的机器原子间势的性能

作为上述内部坐标变换实用性的最后一个应用，我们将分析机器学习的原子间势在使用近似简并振动扰动理论计算甲醇的非谐频率时的表现。甲醇是一个有用的测试系统，可以用来探究MLIP在处理非谐性方面的性能，因为它是一个简单的、共价结合的有机分子，包含一个孤立的非谐OH伸缩、一组近似简并的CH伸缩，并且有高质量的基准数据支持进一步研究，特别是扭转自由度中的振动-旋转耦合（103-109）。在表3中，我们提供了使用近似简并的二阶振动扰动理论和四次力场，在B3LYP和ωB97X-D3理论水平上计算的甲醇的OH伸缩、CH伸缩、HCH弯曲和HOCH扭转频率，以及AIMNet2机器学习的原子间势，模型的详细信息如上所述（48）。相应的模式向量在图S3中标记为模式1–7和12。选择B3LYP是因为文献中的先例（110），而AIMNet2势是根据ωB97M-D3理论水平的电子结构计算进行训练的，因此包括了这两个理论水平。计算使用了四次力场在笛卡尔位移坐标中（表3中标记为Cart.）和内部坐标展开（标记为Int.）。这些计算的笛卡尔坐标和Z矩阵以及非谐计算的完整输入和输出都包含在这项工作的补充数据存储库中（74）。在没有共振处理的情况下，坐标系统的选择不会影响最终频率和强度，如表S1所示（59,111），然而当考虑系统中的共振时，差异就会出现。首先，应当注意到ωB97X-D3和B3LYP在拉伸和弯曲频率上的差异在不同处理方法之间最多只有6厘米-1，而AIMNet2的结果差异可能达到70厘米-1。需要记住的是，AIMNet2并不是针对非谐数据进行训练的，因此这个结果仅仅表明需要更多的训练才能获得可靠的非谐结果。尽管如此，情况实际上更为复杂，这在扭转处理中表现得尤为明显。在谐波层面，三种方法都得到了超过200厘米-1的扭转频率，其中B3LYP和ωB97X-D3分别得到291厘米-1和298厘米-1的频率，而AIMNet2得到的频率为209厘米-1。一旦包括非谐校正，ωB97X-D3的频率降至约12厘米-1，这表明即使在AIMNet2训练所使用的理论水平上，该系统的这一模式也没有得到很好的处理。值得注意的是，相对于ωB97X-D3，内部坐标展开在AIMNet2水平上提供了更可靠的结果，因为扭转与其他坐标的耦合减少了。然而，无论在哪种理论水平上，结果似乎都不太可靠。

表3. 不同理论水平和坐标系统下甲醇的谐波和非谐频率比较
| 理论水平 | B3LYP/cc-PVTZ | ωB97X-D3/cc-pVTZ | AIMNet2/ωB97M-D3 |
|-----------------|--------------|-----------------|--------------|
| OH | 382 | 364 | 391 |
| CH3 | 108 | 298 | 299 |
| CH30 | 402 | 287 | 297 |
| CH2 | 994 | 281 | 281 |
| HCH | 150 | 146 | 146 |
| HCH14 | 98 | 145 | 145 |
| HCH147 | 97 | 144 | 144 |
| HCH148 | 95 | 145 | 145 |
| HCCH | 291 | 222 | 222 |
| | 206 | 276 | 276 |

a 使用几乎简并的振动微扰理论，结合哈密顿量在笛卡尔位移坐标（Cart.）和内部坐标（Int.）的展开，并在B3LPY/cc-pVTZ和ωB97X-D3/cc-pVTZ理论/基组水平上优化几何结构以及使用ωB97M-D3计算训练的AIMNet2机器学习势，计算得到的谐波（Harm.）和非谐频率（Anharmonic）频率（单位：厘米-1）。笛卡尔和内部结果之间的偏差是由微扰理论中的强耦合引起的，这些耦合是通过变分方法处理的。

为了解决扭转自由度处理的问题，我们将进一步使用内部坐标展开，并应用从反应路径正常模式分析中得出的技术，将扭转自由度从所有理论水平的展开中投影出去。操作上，这是通过将正常模式表示为内部坐标的展开，然后移除沿扭转自由度的分量来构建一组3N-7个扭转投影模式来实现的。然后使用这些模式来评估非谐频率，就像之前一样。为了更深入地了解MLIP与DFT理论水平之间的比较，我们将把这种分析扩展到平衡扭转角度之外的几何结构。为此，我们首先将扭转角度固定在偏离其平衡值某个程度，然后在剩余的自由度上进行优化，再次利用上述高效的内部坐标变换来转换笛卡尔坐标梯度。完成这一步后，我们应用之前的反应路径投影方法来获得没有沿扭转自由度分量的模式，并在这个子空间中进行振动微扰理论计算。类似的方法也被用来在鞍点处获得VPT2能量。表4提供了AIMNet2 MLIP在0°、30°和60°扭曲角度下的结果，所有理论水平的结果分别列在表S2和S3中。即使在扭曲的情况下，我们也看到笛卡尔坐标和内部坐标处理之间的差异仍然很小，因此对于所有后续结果，只报告内部坐标的结果。在没有扭转的情况下，笛卡尔和内部展开之间一致性的提高是笛卡尔表示中模式间耦合被高估的一个例子，也是如何规避这些高估的一个例子。

表4. 不同扭曲角度和坐标系统下甲醇的谐波和非谐频率比较
| 扭曲角度 | 0° | 30° | 60° |
|---------------|--------------|--------------|--------------|
| OH | 391 | 366 | 391 |
| CH3 | 135 | 296 | 297 |
| CH30 | 75 | 297 | 297 |
| CH30 | 308 | 297 | 297 |
| CH30 | 302 | 287 | 288 |
| CH30 | 303 | 288 | 288 |
| HCH | 150 | 146 | 146 |
| HCH14 | 98 | 145 | 145 |
| HCH147 | 97 | 144 | 144 |
| HCH148 | 95 | 145 | 145 |
| HCCH | 291 | 222 | 222 |

为了进行更全面的比较，我们将这种分析扩展到所有对称不同的扭曲范围。在图7中，反CH伸缩的未受扰动频率作为扭曲角度的函数绘制出来。未受扰动频率用于在不同理论水平之间进行更有用的比较。为了在MLIPs的空间中进行更丰富的比较，我们在分析中包括了MACE-oMOL系列势能。在性能方面，受限优化和部分四次力场评估使用MACE的速度大约是AIMNet2的10倍，但这是在没有GPU的8核机器上进行的。当可以使用GPU加速时，这两种方法的速度都可能会显著提高。我们清楚地看到，虽然两种DFT理论水平的这种频率都有所下降，但AIMNet2的伸缩频率却有所增加，且这种效应随着扭曲角度的增加而增大。相比之下，MACE的反CH伸缩频率有所下降，但其幅度比任何DFT理论水平都要大得多。所有理论水平的这种模式的特性都从平衡结构中的纯反CH伸缩平滑地转变为 eclipsed 结构中的反对称CH伸缩。

图7. 在不同扭曲角度下评估的甲醇反CH伸缩的非谐频率。

因此，评估不同DFT理论水平和MLIP处理之间哪些模式表现相似以及这些差异的起源就变得有趣了。我们将从扩展CH伸缩的讨论开始。在图8中，为不同的电子结构处理方式绘制了所有三种CH伸缩频率在不同扭曲角度下的变化。对于两种DFT理论水平，ν2和ν3在频率上区分明显，但使用AIMNet2时这些频率非常相似，而使用MACE时分裂被显著高估了。在谐波层面，所有三种方法总体上是一致的，如图S4所示，但通过评估受限模式空间中的非谐频率，可以洞察非谐频率趋势的起源。

图8. 不同扭曲角度和理论水平下三种CH伸缩频率的比较。能量、梯度和Hessian矩阵分别使用（A）B3LYP/cc-pVTZ理论/基组（B）ωB97X-D3/cc-pVTZ理论/基组（C）文中给出的特定模型的AIMNet2机器学习原子间势能（D）文中给出的特定模型的MACE-OMOL原子间势能进行评估。未标记的曲线是在完整的反应路径模式空间中评估的。标记表示在CH和OH伸缩子空间中评估的相应值。

在图8中带有标记的线条中，CH伸缩频率作为扭曲角度的函数绘制出来，同时使用仅包含OH和CH伸缩的模式空间进行微扰理论校正。通过这样做，我们能够在使用MLIPs评估非谐校正时消除伸缩与弯曲和其他低频自由度之间的耦合误差。值得注意的是，在不同理论水平上，ν3和ν4在形状上的定性变化方面都不受子空间限制的强烈影响，这表明对这些模式的校正主要来源于对角非谐性和与其他伸缩的耦合。相比之下，对于AIMNet2，当包括所有模式时ν2的频率会增加，但在任一子空间中频率都会降低。对于任何DFT理论水平或MACE，我们看到子空间限制对扭曲引起的频率变化影响相对较小，增加了大致恒定的频率增加，或相应的非谐性恒定减少。有趣的是，在伸缩子空间内评估的AIMNet2 CH伸缩频率在伸缩之间的分裂方面与DFT理论水平在定性上是一致的。这表明使用内部坐标变换将模式限制为伸缩和高频弯曲的线性组合可能提供了一种从机器学习原子间势能获得合理非谐校正的途径。

接下来讨论剩余的伸缩自由度，在图9中，OH伸缩频率作为扭曲角度的函数在理论水平和反应路径模式空间中绘制出来。在这里，令人惊讶的是，当在子空间中包括所有模式时，AIMNet2的非谐OH伸缩频率与B3LYP的OH伸缩频率在定性上非常吻合，而随着进一步限制，这种一致性会降低。对于MACE，也观察到了与非CH伸缩变化相同的高估非谐性的现象。

图9. 在不同扭曲角度和理论水平下评估的甲醇OH伸缩的非谐频率。标记表示在CH和OH伸缩子空间中评估的相应值。

与CH伸缩不同，即使在谐波水平上，OH伸缩也与DFT理论水平有显著差异，如图S5所示。在非谐水平上表面的一致性实际上是谐波频率的过度校正。当考虑OH伸缩的纯非谐贡献时，对于任何DFT理论水平，模式空间限制的效果是对非谐校正的恒定偏移，而对于AIMNet2势能，非谐校正随着子空间限制的增加而显著变化。这在图S7中得到了证明。特别是，这表明AIMNet2的结果将OH伸缩与系统的其他模式耦合得太强，尽管在将模式表示为内部坐标的展开时可以看到OH伸缩模式与系统的其他内部坐标没有混合。即使使用自动微分的第三导数，这种效应仍然存在，表明这不是数值微分的偶然现象。这可能与DFT和MLIP理论水平在谐波水平上的不一致有关。即使在限制到伸缩模式子空间的情况下，非谐校正之间的定性不一致仍然最小化，这表明使用具有三次导数（可能经过缩放的）的更高理论水平的谐波频率在MLIP水平上提供了有效的非谐校正途径。

总之，我们展示了如何实现和应用一种简单的组合方法来构建内部坐标导数。这种方法可以轻松与进一步的坐标变换结合使用，并提供了一个处理此类坐标变换的通用算法。以这种方式构建的导数在数值上是稳定的，并且比相应的数值微分计算在计算效率上更高，并且随着系统维度的增加而更有利于缩放。坐标系统之间的可微转换允许在内部坐标空间中进行更丰富的插值，而将复杂的笛卡尔空间运动表示为简单的、可能具有对称适应性的内部坐标函数，简化了非共价结合复合物中模式间振动耦合的探索。最后，高效构建高阶坐标变换导数使得可以在低频坐标上使用振动微扰理论进行反应路径形式主义的计算。这样的计算对现代机器学习原子间势能（MLIPs）的性能提供了强有力的测试，这些势能通常最多只能再现谐波频率。虽然使用MLIP沿扫描坐标获得的非谐校正可能与通过密度泛函理论获得的校正有显著差异，但通过将耦合模式的空间限制为仅包含高频模式，这种差异可能会被最小化。通过限制用于构建系统模式的坐标空间，这种差异可能会进一步减少。应当注意的是，本工作中详细的方法并没有完全解决内部坐标处理中坐标系统选择的核心问题。如前所述，冗余坐标处理是解决这一问题的一种可能方法。(26,29) 局部正常模式可能提供另一种解决方案，该方法可以通过内部坐标系统得到补充。(115?119) 高效地计算导数可以为给定配置下的问题选择最佳内部坐标提供途径，而这里提供的坐标系统转换功能可以连接局部最优的坐标系统。同样，本研究中对MLIPs和DFT的比较仅反映了两种最先进的机器学习原子间势能在当前时刻的性能。随着进一步训练的进行和非谐训练数据的纳入，这些比较结果需要更新。最后需要指出的是，本研究中提供的插值方法只是众多可能选择中的一种，利用能量信息或直接采样反应路径的插值方法可能会表现出更好的性能。在其他应用中，我们利用这种坐标构建方法进行了基于内部坐标的从头算分子动力学模拟，并研究了大型系统多维光谱中的光谱特征。已有使用张量运算的数值高效实现方法(76,78)，同时本研究中使用的所有代码和数据也存储在相应的仓库中(74)。上述方法可以轻松扩展到包括键长、角度和二面角之外的其他坐标类型，还实现了诸如垂直于平面的摆动和平面内的摇摆运动等各种隐式定义的坐标。以这种方式计算的导数和坐标可以集成到其他软件包中，特别是用于评估二阶以上的性质。(120?122) 通过使用GPU加速张量运算，还可以进一步提高计算效率。