如何从RDKit选合适函数计算原子RMSD？结果差异解析

问题 #

我有两个存储在.xyz文件中的分子构象，原子顺序完全一致，需要计算原子RMSD。使用RDKit 2023.03.2（Python 3.10.12）的两个方法得到了不同结果：

from rdkit.Chem.rdmolfiles import MolFromXYZFile
from rdkit.Chem.AllChem import GetBestRMS, AlignMol

mol1 = MolFromXYZFile('mol1.xyz')
mol2 = MolFromXYZFile('mol2.xyz')

print(GetBestRMS(mol1, mol2))
# 输出：
>>> 1.4724067538388081

print(AlignMol(mol1, mol2))
# 输出：
>>> 1.6849895064063607

同时使用另一个工具计算，结果与AlignMol一致：

calculate_rmsd mol1.xyz mol2.xyz
>>> 1.68498950637163

请问哪个结果更具参考性？原因是什么？

输入的.xyz文件内容 #

mol1.xyz

53
mol1
C -1.2878 0.5763 2.5081
C -0.0940 0.1341 1.9360
C 1.1741 0.0622 2.6960
C 2.4687 0.3699 2.3374
C 3.5742 0.1850 3.3022
C 3.2088 -0.3028 4.6335
C 4.1333 -0.5219 5.6874
C 3.7003 -0.9868 6.9197
C 2.3568 -1.2459 7.1454
C 1.4287 -1.0423 6.1369
C 1.8659 -0.5789 4.9095
C -0.1188 -0.2426 0.5911
C -1.2905 -0.1922 -0.1478
C -2.4894 0.1923 0.4610
C -2.4687 0.5937 1.7888
O -3.6751 0.1816 -0.2061
O 0.8942 -0.4005 3.9555
O 4.7477 0.4131 3.0421
O 5.4703 -0.2990 5.5706
O 1.9479 -1.7025 8.3635
H -3.7342 -0.5414 -0.8759
H 0.7834 -0.5858 0.0965
H -1.3083 0.9226 3.5342
H 5.7190 0.0308 4.6726
H 4.4198 -1.1495 7.7154
H 0.9693 -1.8444 8.3853
H 0.3762 -1.2432 6.3074
C 3.2716 2.0790 0.6209
C 2.6150 2.5024 -0.7279
H 3.0594 2.8509 1.3592
C 5.1937 1.8062 -0.8108
C 3.4610 3.5968 -1.3656
H 1.6216 2.8999 -0.5251
C 4.9268 3.0978 -1.5712
H 4.6973 0.9805 -1.3186
H 3.0327 3.8733 -2.3278
H 5.6075 3.8633 -1.2015
O 2.6447 0.8584 1.0546
O 4.6471 1.9473 0.5163
O 2.4942 1.3966 -1.6177
H 2.0475 1.6913 -2.4495
O 3.4694 4.7317 -0.5048
H 4.0197 5.4470 -0.9094
O 5.1997 2.8790 -2.9520
H 6.1446 2.6081 -3.0614
C 6.6626 1.4789 -0.6410
H 7.0930 1.2618 -1.6264
H 6.7632 0.5731 -0.0304
O 7.4024 2.5363 -0.0454
H 8.3513 2.2694 0.0342
H -3.3856 0.9228 2.2665
O -1.2365 -0.5216 -1.4668
H -1.8719 -0.0059 -2.0191

mol2.xyz

53
mol2
C -1.0014821846 0.9626299172 2.5573785856
C 0.0496599188 0.3470275666 1.8824250702
C 1.3056173719 0.0052166453 2.5515440934
C 2.5351322319 0.0505765352 1.9720897612
C 3.7184917469 -0.3826176514 2.6867783559
C 3.5033977186 -0.7447874668 4.0680019924
C 4.5710689863 -1.1496900096 4.909962476
C 4.331012171 -1.4643148376 6.2323969408
C 3.0322050126 -1.3853039034 6.7309438348
C 1.9606405974 -0.9976519507 5.9352251045
C 2.2200803246 -0.6844866633 4.6126606673
C -0.0528746879 0.1344520608 0.507116763
C -1.1219184368 0.628333423 -0.1989442937
C -2.1750118806 1.258326425 0.4749709536
C -2.1156576788 1.3988437299 1.8542332231
O -3.2301270978 1.7355182031 -0.214346733
O 1.1556575067 -0.3281667805 3.8507217244
O 4.8279324073 -0.4347069578 2.1341908715
O 5.8133601226 -1.2251560697 4.434103765
O 2.8666906888 -1.7097569209 8.0310352514
H -3.0713245174 1.5796359161 -1.1554400192
H 0.7299831889 -0.3731461741 -0.0302850043
H -0.9388459439 1.1233243801 3.6248068287
H 5.7859501611 -0.9769864865 3.4829713595
H 5.1419673986 -1.7725754551 6.8757055943
H 1.9400722545 -1.62970029 8.2788987003
H 0.9493477107 -0.9458446957 6.3149706343
C 2.5307872509 1.7951110049 0.3561613886
C 2.3639651114 1.9388506917 -1.1814462985
H 1.6479610269 2.1607641485 0.879358051
C 4.8869466139 2.2418290229 0.3490523926
C 3.6720879214 2.3746048635 -1.8740440416
H 1.6401825681 2.7355051455 -1.3532583007
C 4.8654555682 1.7945686214 -1.1212738445
H 5.2822036245 1.4119392135 0.942689672
H 3.6552784387 1.9153847486 -2.8655390874
H 5.7929240643 2.1157036397 -1.6054301121
O 2.7255144868 0.4227796783 0.6686915653
O 3.5825033241 2.5760176524 0.833928119
O 1.8260857813 0.7731060515 -1.7765584379
H 2.5260815562 0.1043987317 -1.7561092302
O 3.7348130031 3.7633538109 -2.0506663046
H 4.3080557405 4.1705571602 -1.3821917325
O 4.7612665245 0.3796800179 -1.2361215032
H 4.7275141397 -0.0164450412 -0.3543935804
C 5.7429738789 3.4584223359 0.6600871543
H 6.7326676961 3.3376535301 0.2210362026
H 5.8537268537 3.5374523689 1.7445044248
O 5.1975857627 4.658562741 0.1413912577
H 4.3744067965 4.8217632085 0.6175036702
H -2.9382623189 1.8891314043 2.3568277443
O -1.2092892555 0.5529415687 -1.5558948255
H -0.3106612492 0.4735711915 -1.9168208442

答案 #

核心差异原因 #

两个RDKit函数的行为本质不同：

• GetBestRMS：忽略输入的原子顺序，自动寻找两个分子中拓扑等价的原子对（即化学环境相同的原子），通过重新配对这些原子计算能达到的最小RMSD值。即使原子顺序一致，它也会尝试拓扑重匹配，因此得到的是理论上的最小可能RMSD。
• AlignMol与第三方工具：严格按照输入的原子顺序配对，仅做刚性旋转和平移对齐，不改变原子对应关系，计算的是固定原子配对下的RMSD。

哪个结果更具参考性？ #

取决于你的计算需求：

1. 若需要固定原子顺序的构象差异（比如两个构象来自同一分子动力学轨迹、原子顺序完全一一对应），AlignMol和第三方工具的结果（1.68）更准确，它完全遵循给定的原子配对逻辑，符合常规构象轨迹分析需求。
2. 若需要两个分子的最优拓扑匹配差异（比如比较不同来源的同一分子，原子顺序可能存在差异，或需要找到最相似的原子配对），GetBestRMS的结果（1.47）是合理的，它代表了两个分子在拓扑等价原子配对下的最小构象差异。

验证说明 #

你的两个.xyz文件原子顺序一致，但GetBestRMS仍得到更小值，是因为分子中存在多个拓扑等价的原子（比如对称位置的氢原子、相同化学环境的氧原子），它通过重新配对这些原子，找到了能让RMSD更小的对齐方式。

内容的提问来源于stack exchange，提问作者vahahahart