计算化学公社
标题: 在VMD中计算RMSD衡量两个结构间的差异以及叠合两个结构 [打印本页]
作者Author: sobereva 时间: 2015-5-11 01:04
标题: 在VMD中计算RMSD衡量两个结构间的差异以及叠合两个结构
在VMD中计算RMSD衡量两个结构间的差异以及叠合两个结构
Calculating RMSD in VMD to measure the difference between two structures and to superimpose two structures
First release: 2015-May-11 Last update: 2020-Feb-25
1 原理
RMSD（Root mean square displacement/deviation）搞分子模拟的人都很熟悉，但是搞量子化学的人可能好多都不知道。有人问怎么衡量激发态结构相对于基态的改变，虽然可以挨个比较键长键角的变化，但是若想用一个值衡量整体的结构改变程度，计算RMSD是最恰当也是最省事的，这里就简单说一下怎么算。对任意两个结构间都可以计算RMSD，比如可以衡量两种极小点构型的差异、电子激发/电离/外加电场前后几何结构的整体改变程度、形成复合物后单体结构的变化等等。
RMSD定义为
(, 下载次数 Times of downloads: 161)
其中i循环所有原子，x_i和x_i'分别是第i个原子在第一个结构和第二个结构中的x坐标，y、z类似。
计算两个结构间的RMSD之前必须先进行叠合（align），相当于令一个结构平移和旋转来最小化两个结构之间的RMSD，不做这一步的话计算出的RMSD是不能衡量两个结构内部的差异的。
能计算RMSD、做叠合的程序不计其数，下面用免费的VMD程序（http://www.ks.uiuc.edu/Research/vmd/）的1.9.3版来演示，而且还将体现如何通过叠加图直观地展现两个结构之间的差异。
2 实例：中性和阳离子状态下的丙烷
本例演示如何计算丙烷在中性状态下和在+1阳离子状态下优化的结构之间的RMSD值，并且绘制出叠合后的图。
首先把中性的丙烷用量子化学程序优化，把优化后的结构保存成VMD能认的一种格式，比如pdb、mol2、xyz等。这里我们用很常用的xyz格式，见《谈谈记录化学体系结构的xyz文件》（http://sobereva.com/477）。我们创建一个名为neutral.xyz的文本文件，把优化后的中性丙烷坐标拷进去，并在第一行写上原子数，第二行是注释行（内容随便写）。具体内容如下
11
Neutral
6 0.000000 1.277193 -0.259856
1 0.884678 1.321887 -0.907445
1 -0.884678 1.321887 -0.907445
1 0.000000 2.176504 0.367166
6 0.000000 0.000000 0.586502
1 0.877607 0.000000 1.247354
1 -0.877607 0.000000 1.247354
6 0.000000 -1.277193 -0.259856
1 0.000000 -2.176504 0.367166
1 -0.884678 -1.321887 -0.907445
1 0.884678 -1.321887 -0.907445
类似地也优化丙烷阳离子，把结构存为cation.xyz，内容如下
11
Ionized
6 1.249671 -0.343789 0.000081
1 1.315464 -0.947246 0.909976
1 1.314553 -0.947793 -0.909563
1 2.180924 0.290109 -0.000844
6 0.192643 0.673245 0.000031
1 0.022018 1.237089 0.919487
1 0.022792 1.237187 -0.919539
6 -1.442724 -0.286303 -0.000053
1 -2.107101 0.577835 -0.000705
1 -1.372956 -0.853783 -0.925484
1 -1.373233 -0.852316 0.926316
PS：如果你用的是Gaussian，把优化任务最后一帧的结构转化为xyz格式最简单的方法是用Multiwfn（http://sobereva.com/multiwfn）。把Multiwfn的settings.ini文件里的iloadGaugeom设为1，然后用Multiwfn载入Gaussian优化任务的输出文件，选择主功能100的子功能2，就可以看见导出xyz文件的选项，选择即可。
启动VMD，把neutral.xyz和cation.xyz都拖到VMD Main窗口里载入。然后选择Extensions - Analysis - RMSD Calculator，把左上角的文本框的内容改为all，取消选择Backbond only复选框，然后点击Align，在图形窗口你会发现两个结构已经叠合了，重叠度较高。然后再点RMSD按钮，在RMSD Tool窗口下方显示的Total RMSD就是两个结构间的RMSD值。当前窗口看到的情况如下所示，可见RMSD是0.093埃。
(, 下载次数 Times of downloads: 182)
做叠合和计算RMSD的时候，我们可以只对指定部分进行。比如我们不输入all，而是输入比如serial 2 to 4 8，那么点Align按钮时只会根据2、3、4、8号原子来做align，点击RMSD按钮时也只会计算这些原子间的RMSD。如果你对VMD的选择语句不熟悉，建议参看《VMD里原子选择语句的语法和例子》（http://sobereva.com/504）。
VMD还可以把两帧结构用不同颜色绘制出来便于直观比较。做法是选Graphics-Representation，在Selected Molecule里先选择对应中性的体系，把Coloring Method设为Color ID，在旁边选择13 Mauve，并把Thickness设为4使细线显示得更粗。然后在Selected Molecule里选择对应阴离子的体系，以类似操作把颜色设为蓝色并加粗细线。最后在文本窗口输入color Display Background white命令把背景改为白色，此时在图形窗口会看到下图，可见结构差异体现得十分清晰直观
(, 下载次数 Times of downloads: 168)
3 实例：不同的取代苯体系
有时候我们想对两个主体结构相同，但某些部分不同的两个分子的局部区域计算RMSD来衡量局部结构的差异。本例就用苯酚以及间硝基苯腈为例进行说明，二者在B3LYP/6-31G*优化后的结构如下
(, 下载次数 Times of downloads: 164)
二者的差异在于苯环上连的基团。本例我们要让两个分子的苯环部分相互叠合（只需考虑碳原子）并计算RMSD。为了能够正确叠合，一定要注意被叠合的部分要满足两个条件：(1)原子数相同 (2)原子顺序一致。由上图可见，这两个体系中的苯环上的碳都是按连接顺序排列的，所以可以直接做叠合和计算RMSD。如果其中一个体系的苯环上的碳是比如1,4,2,6,5,3这样交错排列的就没法直接做叠合了，而需要在叠合之前自行编辑结构文件，调整原子顺序使得条件(2)被满足。
还是先载入这两个结构到VMD里并进入RMSD Calculator，左上角的文本框输入serial 1 to 6（对应两个分子的苯环上的碳原子序号），取消Backbone only复选框，然后点Align。按理说此时两个分子的苯环应该被叠合了、之后再点击RMSD就完事了，然而当前的情况比较特殊，图形窗口里目前看到如下情况
(, 下载次数 Times of downloads: 165)
显然还没有真正叠合好。这是VMD用的算法的缺陷，类似于优化落到了局部极小点而非全局最小点的意味。解决办法就是手动旋转其中一个分子，令其苯环部分的各个原子与另一个分子的苯环部分相应序号的原子整体比较接近，类似于提供一个更好的初猜。调节结构具体的做法是选Mouse - Move - Molecule，之后按住alt键并用鼠标左键拖动某原子就可以平移整个体系，按住alt+shift并用鼠标左键拖动某原子可以旋转整个分子，按住alt+shift并用鼠标右键拖动某原子可以令整个分子沿屏幕旋转。等调节得差不多了，再次按Align按钮，就会发现两个分子的苯环部分确实很好地叠合在了一起，然后再点击RMSD按钮，此时情况如下所示。
(, 下载次数 Times of downloads: 147)
可见两个分子的苯环碳原子间的RMSD仅为0.013埃，体现出苯环非常刚，取代基对这部分结构的影响可以忽略不计。
作者Author: 赵云跳槽 时间: 2016-3-22 10:56
每次运算完一次RMSD计算后，都需要重启一次VMD，不然报错，怎样清楚上一次作业的痕迹？
作者Author: jackie 时间: 2017-6-20 16:23
老师我有两个问题：
1.在这里讲的是两帧构型的的计算，那如果我想做更多个，比如我在.xyz文件中保存100个分子，根据上述描述进行同样动作的操作，VMD可不可以实现对于多个构型RMSD值的计算？
2.老师您提到，在计算RMSD值之前要先对分子进行align，也就是通过平移翻转尽可能和第一帧构型接近，那对于这一部分程序内部是通过什么算法实现的，有没有计算RMSD的align的代码，我想学习一下？
谢谢老师~
作者Author: KiritsuguPapa 时间: 2017-6-20 22:40
关于第二个问题 可以参考https://github.com/charnley/rmsd 里面也有算法的介绍
第一个问题 如果你会python的话 调用上面的那个程序也很方便
作者Author: jackie 时间: 2017-6-21 10:26
So cool!!!非常感谢
作者Author: sobereva 时间: 2017-6-21 13:48
RMSD Trajectory Tool里，点了align，所有帧（xyz里每个构型）都会相对于第一帧（默认情况下如此）做align，然后可以直接把所有帧的RMSD都算出来并导出
作者Author: jackie 时间: 2017-6-23 10:41
老师 为什么我把好多帧的蛋白构型文件用VMD打开通过Graphics-Representation-Drawing method选项没法变成cartoon构型？文件已上传
作者Author: sobereva 时间: 2017-6-23 15:50
上传大文件前一定要压缩以节约论坛空间、加快下载速度
xyz文件里没有残基、原子名信息，不足以令VMD通过stride程序判断二级结构
作者Author: 陈丢丢 时间: 2018-8-6 16:23
sobereva老师，打开VMD后，Extensions-Analysis目标下面有RMSD calculator，RMSD Trajectory Tool，RMSD visualizer Tool三个选项，这三个选项计算RMSD有什么区别吗？
作者Author: sobereva 时间: 2018-8-7 04:12
陈丢丢 发表于 2018-8-6 16:23
sobereva老师，打开VMD后，Extensions-Analysis目标下面有RMSD calculator，RMSD Trajectory Tool，RMSD vi ...
原理都一样，只不过具体用处不同。一般只用得到RMSD Trajectory Tool
作者Author: 陈丢丢 时间: 2018-8-7 09:53
好的，谢谢sobereva老师！
作者Author: sobereva 时间: 2020-2-25 07:54
大幅重写了本文，使用法更为灵活，并且加入了对原子数不同的两个体系做叠合并计算RMSD的说明
作者Author: 小土豆 时间: 2020-5-19 11:55
你好老师，我想问一下，我有一个目标单体和一个已知三聚体结构，想通过叠加做出目标的三聚体结构，怎么操作呢。
作者Author: sobereva 时间: 2020-5-21 08:01
我不知道你说的“已知三聚体结构”具体指什么，里面的分子和当前单体有什么差别
如果里面的分子没法和当前单体相叠合，不可能靠叠加得到当前单体的三聚体结构
作者Author: Aristotler 时间: 2020-6-23 16:15
老师，我想引用一下RMSD计算公式，但也没找到您写的这个，有参考可以提供一下吗？
作者Author: sobereva 时间: 2020-6-24 16:35
不用引用。都是众所周知的东西
作者Author: wangxm 时间: 2021-12-23 16:33
请问如何将重叠之后的导出为图片格式呢？
作者Author: sobereva 时间: 2021-12-23 17:06
看此文里关于保存图像的说明
用Multiwfn+VMD做RDG分析时的一些要点和常见问题
http://sobereva.com/291（http://bbs.keinsci.com/thread-1206-1-1.html）
作者Author: wangxm 时间: 2021-12-23 17:21
谢谢老师
作者Author: 一条君 时间: 2022-3-25 16:37
老师，请问md中一般大于多少原子数或者分子量后需要看rmsd来判断是否平衡呢，最近在模拟低聚物，大约1500分子量，遇到了这方面疑问，谢谢老师
作者Author: sobereva 时间: 2022-3-25 19:37
跟原子数没关系，取决于体系
不同体系判断平衡方式不同，此文说得很充分了
谈谈怎么判断分子动力学模拟是否达到了平衡
http://sobereva.com/627（http://bbs.keinsci.com/thread-27122-1-1.html）
作者Author: smutao 时间: 2022-3-26 10:46
在pymol中 可以使用 pair_fit 命令 轻松实现任意两个分子的叠合以及rmsd计算
大致的语法为 pair_fit a_1,b_1,a_2,b_2
https://pymolwiki.org/index.php/Pair_fit
作者Author: qcn1211 时间: 2023-5-29 10:37
请问用VMD 计算出的RMSD值为多少才说明这两个结构相近？
作者Author: sobereva 时间: 2023-5-29 10:40
没一刀切的标准。结合叠合的图像凭直觉判断最靠谱
作者Author: zjxitcc 时间: 2023-5-29 10:46
说说你的体系多少个原子，你比较的是所有原子还是部分原子
作者Author: qcn1211 时间: 2023-7-20 17:32
请问sob老师，用molclus调动高斯计算筛选低能构象，结合VMD align, RMSD值为多少时说明结构相近？进而需要isostat重新归簇。
作者Author: sobereva 时间: 2023-7-20 19:22
用isostat默认的结构偏离标准就够了，没必要结合VMD
作者Author: qcn1211 时间: 2023-7-20 20:12
谢谢sob老师的指点
作者Author: granvia 时间: 2023-10-2 10:04
求教一下：如果两个结构的原子排列顺序不一致，有没有不手动而完全自动化计算RMSD的算法？
作者Author: sobereva 时间: 2023-10-3 14:12
可以通过距离矩阵对比差异
molclus的isostat比较结构差异的时候就是基于距离矩阵对比的
作者Author: granvia 时间: 2023-10-4 12:22
是不是比较距离矩阵的特征值就够了？
作者Author: sobereva 时间: 2023-10-4 14:31
这是一个办法，也可以把距离值按从大到小排序成为矢量，比较两个矢量的差异，molclus的isostat比较结构相似性的时候就是这么做的
作者Author: granvia 时间: 2023-10-4 16:54
sobereva 发表于 2023-10-4 14:31
这是一个办法，也可以把距离值按从大到小排序成为矢量，比较两个矢量的差异，molclus的isostat比较结构相 ...
明白啦，谢谢！
作者Author: Uus/pMeC6H4-/キ 时间: 2025-4-8 08:14
本帖最后由 Uus/pMeC6H4-/キ 于 2025-4-8 08:23 编辑
昨天我发的这帖就用此文的方法，简单展示了某个质子转移反应的反应物、过渡态、产物的叠合情况。
做些笔记：在VMD中启用File - Log Tcl Commands to Console后，可见从Extensions - Analysis打开RMSD Calculator窗口对应于命令行menu rmsd on指令，打开RMSD Trajectory Tool窗口对应于命令行menu rmsdtt on指令，打开RMSD Visualizer Tool窗口对应于命令行menu rmsdvt on指令。在窗口内操作时并不会在命令行显示对应指令，不过定义几个工具的原始tcl脚本可以在VMD目录下plugins/noarch/tcl/rmsd1.0或rmsdtt3.0或rmsdvt1.0文件夹内找到。实际上VMD自带的doc/ug.pdf用户指南的12.4节专门介绍了RMSD相关的tcl脚本编写技巧，核心指令就是measure fit $sel1 $sel2返回4x4矩阵，使得选区$sel1的原子改变坐标后与$sel2间的RMSD最小，以及measure rmsd $sel1 $sel2计算RMSD本身。measure fit有一个选项order可以指定作为参比的$sel2里的原子index列表。VMD的RMSD算法来自Kabsch (Acta Cryst. (1978) A34, 827-828)。
有一点不足之处是，这个方法得到的是视角（具体说就是http://bbs.keinsci.com/thread-17814-1-1.html涉及的若干矩阵）不变的前提下原子坐标的变换矩阵，但实际上有时可能需要的是体系的视角变换矩阵。比如对两个结构相似但朝向不同的分子用波函数信息计算实空间函数后将cube文件载入VMD，希望将原子位置和等值面同时在视觉上叠合，这就不能只做原子坐标的矩阵变换而不管格点数据了。不知道除了在产生波函数信息前就把分子对齐以外，有没有后处理上的解决方法。
作者Author: zhouyulin 时间: 2025-5-2 23:50
老师，我想问一下优化结构对比结果后的RMSD怎么保存下来
作者Author: sobereva 时间: 2025-5-3 00:01
如果你是指保存RMSD数值，复制粘贴就完了
如果你是指保存align后的结构，file - save coordinate
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |