使用UniMoVib进行GSVA广义子系统振动分析

1. 引言

简谐振动分析(harmonic vibrational analysis) 是量子化学计算中一项常用的技术手段。一方面，这种振动分析可以给出红外、拉曼等振动光谱。在最常用的高斯程序中，指定关键词freq则可以进行简谐振动分析。如果加上谐振频率校正因子，计算得到的振动频率可以更接近实验测量频率值。相关内容可参见《红外光谱的理论计算》一文。因此简谐振动分析可以以一个较低的计算代价得到质量还不错的振动光谱，但如果需要更精确的振动光谱，则需要考虑非谐振效应。另一个方面，简谐振动分析可以帮助我们确定结构优化过的体系在势能面上驻点（stationary point）的性质。假如振动分析得到的振动频率都是正值，那么此时体系位于能量局部极小点处。假如分析得到的振动频率有一个是虚频而其余都是正值，说明此时体系在一个鞍点上（即反应过渡态）。

对于一个含有N0个原子的孤立分子（比如水分子），它的振动模式是比较容易理解的，一共有3N0−6个（直线型分子为3N0−5）。但是，假如我们把这个孤立的分子跟其他分子放到一起成为一个稳定的复合物（比如水二聚体），这个时候我们将得到3N1−6 个振动模式，其中每个振动模式将不再是简单的局限于那个孤立分子的3N0−6个振动模式，而是会牵涉到其他地方的原子参与其中。就拿孤立的水分子和水二聚体来说，水分子有如下的三个振动模式(N0=3, 3N0−6=3)：(1) HOH夹角弯曲振动，(2)OH键对称伸缩，以及(3) OH键非对称伸缩；对水二聚体来说，它一共有12个振动模式(N1=6)，假如我们想在二聚体体系中找到孤立水分子对应的3个振动模式，将会非常困难，因为此时多数振动模式同时牵涉到的是两个水分子而非只局限于其中一个。为了更具体地说明这一点，我们在ωB97XD/6-311++G(d,p)这一级别下计算了优化过的水二聚体的振动模式，以其中第7和8个振动模式为例（见下图），两个水分子的HOH夹角弯曲振动产生了严重的混合，无法再将这两个振动模式归属到其中任意一个水分子了。

从上面的例子我们看到：假如我们想要在一个复合物体系的结构中得到其中一个子系统结构对应于的相同子系统结构在孤立体系下的简正模式（normal mode）时，直接计算该复合物的振动模式将不再可行。以水二聚体的例子来讲，我们想要得到其中一个水分子对应于孤立水分子的振动模式的时候，直接计算水二聚体的振动模式不再合适。

2. GSVA—广义子系统振动分析

为了解决上一节提到的问题，笔者于2018年开发了广义子系统振动分析方法（Generalized Subsystem Vibrational Analysis, GSVA）。该方法的大致思路为：通过一些矩阵变换从复合物结构（即full system，完整体系）中提取出子系统（即subsystem）片段的有效Hessian矩阵（effective Hessian matrix），然后我们对这一子系统专属的Hessian矩阵进行简谐振动分析，得到的子系统的简正模式（称为内禀片段振动，intrinsic fragmental vibration），则可以和相同子系统在孤立体系下的振动直接比较了。

当然，文献中有报道过其他类似的可以获得子系统振动模式的方法，GSVA有什么特别之处呢？在笔者的文章J. Chem. Theory Comput. 2018, 14, 2558一文中指出，GSVA方法最大的优点在于其计算的子系统的有效Hessian矩阵可以保留完整体系势能面在子系统内任意一个内坐标方向的曲率 (curvature)，该曲率也对应了Konkoli-Cremer局域振动模式理论下的局域振动力常数或称柔性力常数（相关综述请见doi: 10.1002/wcms.1480以及sob老师的博文http://sobereva.com/364）。因此，GSVA相较于其他方法完整地保留了子系统Hessian矩阵应有物理意义。

最近笔者对GSVA方法的计算方法进行了简化，原GSVA方法需要用户手动构建子系统的非冗余内坐标集合，新方法将不再需要这一繁杂的步骤，相关的文章请见 Theor. Chem. Acc. 2021, 140, 31。与此同时，笔者将GSVA的新计算方法写入了UniMoVib程序

https://github.com/zorkzou/UniMoVib

这样我们将可以很方便地进行GSVA分析了。 3. 准备工作

进行GSVA分析需要先安装UniMoVib程序，在笔者之前的《使用UniMoVib+PyVibMS显示其他量化程序振动分析结果》一文中，介绍了UniMoVib程序的安装和基本使用方法。如果后续需要显示子系统振动模式，则还应该安装好PyVibMS插件，相关的教程请见《使用PyVibMS可视化分子和固体中的振动模式》一文。

4. 实例

在这一节中，笔者将对三种不同的GSVA方法应用场景进行介绍：(1) 处于能量极小点体系的GSVA分析；(2) 处于过渡态体系的GSVA分析；(3) 柔性扫描路径下的GSVA分析。

本文涉及的实例文件中，Hessian数据文件均为高斯计算得到的fchk文件。fchk文件可以在普通freq计算得到chk文件后，通过formchk程序（https://gaussian.com/formchk/）转换得到。当然，其他量子化学程序得到的Hessian数据也可以用来做GSVA分析，只要是UniMoVib程序支持的就都可以。

4.1 水二聚体的GSVA分析以及水单体子系统内禀片段振动的PyVibMS可视化

本例涉及的文件可以在

https://github.com/smutao/UniMoVib/tree/master/test/Gaussian09/water-dimer

找到。

在ωB97XD/6-311++G(d,p)级别下优化得到了水二聚体的结构（如图），其中1-3号原子为氢键受体水，4-6号原子为氢键供体水。

假如我们想知道受体水（1-3号原子）的内禀片段振动，则需要这样编写UniMoVib的输入文件：

a test job
 
 $contrl
   qcprog="gaussian"
   ifgsva=.true.
 $end
 
 $gsva
   subsystem="1,2,3"
 $end
 
 $qcdata
   fchk="0000-w2-lm0.fchk"
 $end

其中最关键的是ifgsva关键词需要设置为.true.。在相应的UniMoVib输出文件中，定位到文件最末的 “Generalized Subsystem Vibrational Analysis (GSVA)”，给出了子系统所包含的原子以及这3个原子构成的子系统（即受体水）的3个内禀片段振动模式的频率等信息。

假如我们想知道这3个内禀片段振动模式分别对应什么样的振动，则可以用PyVibMS插件来显示。此时我们只需要在UniMoVib输入文件中添加 ifpyvibms关键词即可：

a test job
 
 $contrl
   qcprog="gaussian"
   ifgsva=.true.
   ifpyvibms=.true.
 $end
 
 $gsva
   subsystem="1,2,3"
 $end
 
 $qcdata
   fchk="0000-w2-lm0.fchk"
 $end

在运行UniMoVib之后，除了主输出文件外，程序还会给出额外的3个文件：(1) system.xyz，为水二聚体结构的xyz文件；(2) pyvibms_mode_full.txt，为水二聚体振动模式的mode文件；(3) pyvibms_mode_subsys.txt，为受体水振动模式的mode文件。假如我们想显示受体水的3个内禀片段振动，只需要在打开PyVibMS插件后先载入system.xyz坐标文件，再载入pyvibms_mode_subsys.txt mode文件即可（见下图）。

通过PyVibMS可视化我们得知，受体水的第2、3个内禀片段振动分别为OH键的对称、非对称伸缩，和孤立水分子的情况一致。

假如我们想得到供体水分子的内禀片段振动，只需要在运行UniMoVib的时候将subsystem关键词的值设置为 "4-6"即可。

4.2 SARS-CoV-2主蛋白酶与抑制剂的反应过渡态的GSVA分析

本例涉及的文件可以在

https://github.com/smutao/UniMoVib/tree/master/test/Gaussian09/gsva-paper-examples/transition-state-gsva

找到。

在本例中，我们计算了SARS-CoV-2病毒的主蛋白酶Mpro与酮酰胺(ketoamide)抑制剂发生反应时的过渡态，结构如下图

该反应的大致机理为蛋白酶的Cys残基（以CH3SH简化模型表示）和抑制剂分子发生了质子转移（S→O），同时S原子和63号C原子成键。

众所周知，反应过渡态结构有且只有一个虚频振动。该体系反应过渡态的虚频频率为−725 cm−1。假如我们对反应中心区域进行GSVA分析，会得到什么结果呢？

假如我们只选取质子转移的部分(S85, H65, O64)来做GSVA分析，结果得到内禀片段振动中并没有虚频振动的存在。假如我们再将C原子考虑进去(S85, H65, O64, C63)，此时GSVA得到的内禀片段振动中便出现了1个虚频振动，频率为−969 cm−1。这说明只有将参与断键、成键的原子都包含进去时，才能得到虚频振动。

假如我们以这4个原子为中心，逐步扩展子系统的区域，则会发现相应内禀片段振动中的那个虚频振动会很快向着完整体系的虚频振动频率收敛（见下表）。

由上表可知，当子系统的原子数目增加到7个时，虚频振动以及可以重复出完整体系96%的虚频振动频率了。此时，选定的7个原子可以理解为决定了该反应机理最重要的部分。

因此，GSVA在分析反应过渡态时，可以帮助我们确定反应中心区域不同原子对反应机理的贡献程度。假如选定的一个子系统已经可以贡献90%甚至更高，那么我们已经可以断定改变这个子系统之外的其他原子对反应能垒的影响微乎其微。

4.3 水二聚体柔性扫描路径下的GSVA分析

本例涉及的文件可以在

https://github.com/smutao/UniMoVib/tree/master/test/Gaussian09/water-dimer-relaxed-scan

找到。

在本例中，我们把例子 4.1 中的水二聚体的两个氧原子（1、4号原子）往左右两侧拉开，同时对其他原子进行弛豫。此时，我们仍然可以对供、受体水分子分别进行GSVA分析，因为对于这两个水分子来说，它们内部的受力均为0（是稳定的）。

下图展示了拉伸O-O长度过程中的能量变化（黑色实线，左侧Y轴）。

同时，我们还对扫描路径中其中4个点进行了GSVA计算，得到的彩色实线部分（右侧Y轴）是供体水分子(4-6)的内禀片段振动，彩色虚线部分是受体水(1-3)的内禀片段振动，最右侧的彩色小箭头是孤立水分子的3个振动频率。绿色、蓝色和红色分别代表了水分子的(1) HOH夹角弯曲振动，(2) OH键对称伸缩，以及(3) OH键非对称伸缩。

从这个例子中，我们不难看出：随着OO距离的逐渐增大，该“二聚体”中的两个水分子的行为越来越趋近于孤立水分子的振动。

柔性扫描是量子化学中一种常用的计算方法，假如我们对两个组分之间的距离、夹角等参数进行了柔性扫描，那么对于单个组分，我们仍旧可以进行GSVA分析。但是需要注意的是，直接对柔性扫描的完整体系（即两个组分的集合）进行振动分析是不合理的，因为此时整个体系的梯度不为0，且还可能出现虚频。

5. 结语

本文介绍了笔者开发的广义子系统振动分析GSVA方法在UniMoVib程序中的使用，以及利用PyVibMS插件对子系统的内禀片段振动进行可视化。GSVA方法是目前唯一一种保留了Hessian物理意义的子系统振动分析方法。该方法不仅适用于能量极小点体系和一阶鞍点反应过渡态体系，同时还适用于柔性扫描路径下的情况。

如果GSVA方法在你的科研中有所帮助，请引用我们的论文：

Tao, Y., Zou, W., Nanayakkara, S. et al. A revised formulation of the generalized subsystem vibrational analysis (GSVA). Theor. Chem. Acc. 140, 31 (2021).

https://doi.org/10.1007/s00214-021-02727-y

下文的引用可选：

Tao, Y., et al. Recovering Intrinsic Fragmental Vibrations Using the Generalized Subsystem Vibrational Analysis. J. Chem. Theory Comput. 14, 2558 (2018).

https://pubs.acs.org/doi/10.1021/acs.jctc.7b01171