序言

针对不同的体系选择不同的计算级别(泛函/基组)做结构相关的计算是任何一个应用计算化学工作者都一直在学习、讨论、测试的问题。与此同时几乎每年都有一些方法被提出来，对于一些体系来说使用新的方法跟过去的方法相比看结构相似度是同时有理论价值和使用价值的。

重叠算法的原理介绍

基础中的基础:SVD分解

在介绍算法前首先需要介绍一下SVD分解的概念，任何一个实矩阵AA(不必是方的)都能被分解为三个矩阵的积:

A=USVTA = USV^T

其中UU和VV矩阵都是正交矩阵，正交矩阵的意思是矩阵和自己的转置互逆，即有:

VVT=I,UUT=IVV^T = I,\quad UU^T = I

而其中SS矩阵是一个方对角阵。

分子最大重叠问题的数学描述和Kabsch算法:

如果拿到两个相当类似的分子，首先要做的要做的是对两个分子进行最大重叠。我们拥有的东西是原子序号相互对应的两组xyz坐标。两个结构的变换从几何角度看无非是平动和转动。我们就是要针对平动和转动进行两块处理。

平动的处理:

对平动的处理是相当容易的，只需要让前后两个分子重心相同。常见的处理是把两个分子的坐标从绝对坐标变成以分子中心(重心)为原点的相对坐标，这样就最大程度上消灭了平动的问题。

转动的处理:

这时候我们得到了两组变换后的xyz坐标，最常见的想法是想利用最经典的3D旋转矩阵来表示，我们用RMSD作为目标函数。

Root-Mean-Square Deviation (RMSD) 的定义:

RMSD=1N∑i=1Ndi2RMSD = \sqrt{\frac{1}{N}\sum_{i=1}^{N}d_i^2}

目标函数最小化，跟上面的RMSD非常像，因为我们的结构是一一对应的，所以把上面的根号和系数1/N1/N去掉也无妨，下面这个式子最小一定是RMSD最小。其中RR是转动矩阵，PiP_i是首结构第ii号原子的xyz坐标，QiQ_i是尾结构第ii号原子的xyz坐标:

E=∑i=1N∥RPi−Qi∥2=∑i=1N(PiTRTRPi−2QiTRPi+QiTQi)E = \sum_{i=1}^{N}\left\|RP_i-Q_i\right\|^2 = \sum_{i=1}^{N}\left(P_i^TR^TRP_i - 2Q_i^TRP_i + Q_i^TQ_i\right)

其中因为旋转矩阵本身是一个正交矩阵，而PP和QQ本身又是给定的定值，所以完全可以忽略不管。同时基于定义的非负性，我们只需要让下面这项最大就可以了:

∑i=1NQiTRPi\sum_{i=1}^{N}Q_i^TRP_i

当然机器学习爱好者在这里估计是毫不犹豫的开始反向传播优化，得到一力降十会的效果。这是3D旋转矩阵RR有三个参数，估计计算可以让机器学习爱好者吃些苦头。最后肯定能得到一样的结果。但这其实是一个不需要机器学习的问题。

注意到，上面的求和可以被写成积的形式，Tr⁡\operatorname{Tr}里面的乘法是可交换的:

∑i=1NQiTRPi=Tr⁡(QTRP)=Tr⁡(RPQT)\sum_{i=1}^{N}Q_i^TRP_i = \operatorname{Tr}(Q^TRP) = \operatorname{Tr}(RPQ^T)

这时候对PQTPQ^T做SVD分解，再次利用Tr⁡\operatorname{Tr}里面矩阵可交换的性质，注意到V/R/UV/R/U都是正交矩阵，它们的乘积必然也是一个正交矩阵。可以记做AA:

Tr⁡(RPQT)=Tr⁡(RUSVT)=Tr⁡(VTRUS)=Tr⁡(AS)\operatorname{Tr}(RPQ^T) = \operatorname{Tr}(RUSV^T) = \operatorname{Tr}(V^TRUS) = \operatorname{Tr}(AS)

SS是一个对角矩阵，AA是一个对角矩阵，最好的情况就是AA为单位矩阵，这样可以让这个Tr⁡\operatorname{Tr}最大。

VTRU=I→R=VUTV^TRU = I \rightarrow R = VU^T

我们直接得到了旋转矩阵RR，至于得到了旋转矩阵之后，剩下的所有工作都变成了体力劳动。上述的所有过程又被称作Kabsch算法1，只利用一个非常高效的SVD分解就得到了最大匹配的旋转矩阵，把一个逆向计算的问题变成一个直接计算的问题。

实例演示

南方科技大学刘柳老师课题组于今年年初在Nature上发表了铂宾催化炔烃三聚化的反应2进行了详细的表征，当然也包括DFT计算反应路径。这篇文章采取了非常甚至有些过度稳健的M06-2X-D3/def2-SVP级别在气相条件下进行了几何优化，然后单点能采取的是M06-2X-D3/def2-TZVP SMD(Benzene)的级别，均用Gaussian软件完成。算出了和实验一致的结果。但如果让我用ORCA重复这篇文章的计算部分，我并不太想用M06-2x。两条理由:

  • M06-2X对格点比较敏感，要得到稳健的结果需要至少用Grid3级别的格点，比默认Grid2格点要多不少计算量。
  • M06-2X是一个杂化泛函，在ORCA里面不如纯泛函跑得快。

基于以上两点问题，容易想到B97-3c这类复合方法，但随之而来的一个问题是，B97-3c类方法能否产生类似的结构？这时候正是用RMSD验证B97-3c能否在一定程度上代替M06-2X-D3/def2-SVP做几何优化的时候。拿本文的1a结构来说，使用B97-3c的话RMSD仅为0.11 A˚0.11\ \text{\AA}，如果在重心匹配的时候压根不考虑HH原子的话，那么RMSD会进一步降低到0.086 A˚0.086\ \text{\AA}，这在100多个原子的体系中完全可以认为是一样的。所以至少在这个阶段，用B97-3c完全OK。但是值得注意的是仅仅是这一个例子还不足以说明替代性，而是要对整条反应势能面有相当的理解和Benchmark才可以，这里纯属抛砖引玉。

Footnotes

  1. Kabsch W. A solution for the best rotation to relate two sets of vectors. Acta Crystallographica Section A. 1976;32(5):922-923. doi:10.1107/S0567739476001873 ↩

  2. Zhang X, Liu L L. Aluminium redox catalysis enables cyclotrimerization of alkynes. Nature. 2026;650:353-360. doi:10.1038/s41586-025-09941-9 ↩