求助:分子重组能较小但 RMSD 值偏大，这种情况是否合理？

gonglijing

想向各位请教一个问题：我计算了分子的重组能，同时通过 VMD 软件计算 RMSD 以表征两个结构间的几何偏差，得到具体数据如下：空穴重组能 = 0.185，对应 RMSD=1.1764；电子重组能 = 0.085，对应 RMSD=4.3345。对应空穴重组能和电子重组能的两个结构的 RMSD 几何偏差图分别见下面左图和右图。

目前发现重组能数值偏小，但对应的 RMSD 值却显著偏大，即两结构间的几何偏差较大，想请教各位该结果是否合理？若存在问题，该从哪些方面排查解决（跑了两边计算，数据提取已核对过，没问题）？恳请指点，谢谢！

wal

感觉没对齐

gonglijing

wal 发表于 2026-1-31 21:36
感觉没对齐

您好，点击align了，难道还有其他对其的方式吗？

snljty2

gonglijing 发表于 2026-1-31 22:15
您好，点击align了，难道还有其他对其的方式吗？

VMD的align对齐是考虑原子顺序打乱的，有可能会收敛到一些局部极小值。如果你两个分子的原子顺序完全一致，可以手写一个对齐，比如这个Python小脚本（需要先用VMD把两个结构存成.xyz格式，然后执行下面命令）：

1. python xxx.py structure1.xyz structure2.xyz

复制代码

这里面没考虑只选择一部分原子对齐的功能，需要的话可以自己改一下。

1. #!/usr/bin/env python3
   
2. 
   
3. import os
   
4. import sys
   
5. import numpy as np
   
6. from numpy.typing import NDArray
   
7. 
   
8. def K_cross_mat(vec: NDArray[np.double]) -> NDArray[np.double]: # cross matrix of a 3d vector
   
9.     return np.array([[0., - vec[2], vec[1]], [vec[2], 0., - vec[0]], [- vec[1], vec[0], 0.]], dtype=np.double)
   
10. 
    
11. def L_quat_mat(q: NDArray[np.double]) -> NDArray[np.double]: # left-multiply matrix of a quaternion
    
12.     return np.block([[q[0], - q[1:4][np.newaxis, :]], [q[1:4][:, np.newaxis], q[0] * np.identity(3, dtype=np.double) + \
    
13.         K_cross_mat(q[1:4])]])
    
14. 
    
15. def R_quat_mat(q: NDArray[np.double]) -> NDArray[np.double]: # right-multiply matrix of a quaternion
    
16.     return np.block([[q[0], - q[1:4][np.newaxis, :]], [q[1:4][:, np.newaxis], q[0] * np.identity(3, dtype=np.double) - \
    
17.         K_cross_mat(q[1:4])]])
    
18. 
    
19. def quat_rot_mat(q: NDArray[np.double]) -> NDArray[np.double]: # 4x4 rotate matrix of a quaternion
    
20.     return L_quat_mat(q) @ R_quat_mat(q).T
    
21. 
    
22. def get_rotation_matrix(coord1: NDArray[np.double], coord2: NDArray[np.double], allow_reflection: bool=False):
    
23.     assert coord1.shape == coord2.shape and coord1.shape[1] == 3
    
24.     H = coord1.T @ coord2
    
25.     U, Sigma, Vt = np.linalg.svd(H)
    
26.     R = Vt.T @ U.T
    
27.     if not allow_reflection:
    
28.         # only half of O(3) is SO(3)
    
29.         if np.linalg.det(R) < 0.:
    
30.             Vt[-1, :] = - Vt[-1, :]
    
31.             R = Vt.T @ U.T
    
32.     return R
    
33. 
    
34. def get_rot_quaternion(coord1: NDArray[np.double], coord2: NDArray[np.double]) -> NDArray[np.double]:
    
35.     assert coord1.shape == coord2.shape and coord1.shape[1] == 3
    
36.     H = coord1.T @ coord2
    
37. 
    
38.     Delta = np.array([H[1, 2] - H[2, 1], H[2, 0] - H[0, 2], H[0, 1] - H[1, 0]], dtype=np.double)
    
39.     tr_H = np.trace(H)
    
40. 
    
41.     N = np.block([[tr_H, Delta[np.newaxis, :]], [Delta[:, np.newaxis], H + H.T - tr_H * np.identity(3, dtype=np.double)]])
    
42. 
    
43.     eig_val, eig_vec = np.linalg.eigh(N)
    
44. 
    
45.     q = eig_vec[:, -1].copy()
    
46.     if q[0] < 0.: q = - q # SU(2) has a 2->1 mapping to SO(3)
    
47. 
    
48.     return q
    
49. 
    
50. def rotate_by_rotation_matrix(coord1: NDArray[np.double], coord2: NDArray[np.double], allow_reflection: bool=False):
    
51.     assert coord1.shape == coord2.shape and coord1.shape[1] == 3
    
52.     coord1 -= coord1.mean(axis=0)
    
53.     coord2 -= coord2.mean(axis=0)
    
54.     R = get_rotation_matrix(coord1, coord2, allow_reflection=allow_reflection)
    
55.     coord1 @= R.T
    
56. 
    
57. def rotate_by_quaternion(coord1: NDArray[np.double], coord2: NDArray[np.double]):
    
58.     assert coord1.shape == coord2.shape and coord1.shape[1] == 3
    
59.     coord1 -= coord1.mean(axis=0)
    
60.     coord2 -= coord2.mean(axis=0)
    
61.     q = get_rot_quaternion(coord1, coord2)
    
62.     coord1 = ((quat_rot_mat(q) @ (np.hstack([np.zeros((len(coord1), 1), dtype=np.double), coord1])).T).T)[:, 1:]
    
63. 
    
64. def align_molecules(ifilename1: str, ifilename2: str, imethod: int=2, allow_reflection: bool=False):
    
65.     # 1 for Kabsch, 2 for quaternion
    
66.     if imethod == 1:
    
67.         align_molecules_Kabsch(ifilename1, ifilename2, allow_reflection=allow_reflection)
    
68.     elif imethod == 2:
    
69.         if allow_reflection:
    
70.             raise ValueError("Quaternion does not allow reflection.")
    
71.         align_moleculs_quaternion(ifilename1, ifilename2)
    
72.     else:
    
73.         raise ValueError(f"Unsupported imethod: {imethod:d}")
    
74. 
    
75. def align_molecules_Kabsch(ifilename1: str, ifilename2: str, allow_reflection: bool=False):
    
76.     ifilename1_split = os.path.splitext(ifilename1)
    
77.     ifilename2_split = os.path.splitext(ifilename2)
    
78.     assert ifilename1_split[1] == ".xyz" and ifilename2_split[1] == ".xyz"
    
79.     elements = np.loadtxt(ifilename1, unpack=True, usecols=0, skiprows=2, dtype=np.dtype("<U2"))
    
80.     coord1 = np.loadtxt(ifilename1, usecols=(1, 2, 3), skiprows=2, dtype=np.double)
    
81.     coord2 = np.loadtxt(ifilename2, usecols=(1, 2, 3), skiprows=2, dtype=np.double)
    
82.     rotate_by_rotation_matrix(coord1, coord2, allow_reflection=allow_reflection)
    
83.     ofilename1 = ifilename1_split[0] + "_aligned.xyz"
    
84.     ofilename2 = ifilename2_split[0] + "_aligned.xyz"
    
85.     with open(ofilename1, "w") as ofile1, open(ofilename2, "w") as ofile2:
    
86.         print(len(elements), file=ofile1)
    
87.         print(len(elements), file=ofile2)
    
88.         print(f"aligned from {ifilename1:s}", file=ofile1)
    
89.         print(f"aligned from {ifilename2:s}", file=ofile2)
    
90.         for iatom in range(len(elements)):
    
91.             print(f" {elements[iatom]:<2s}" + ("    {:13.8f}" * 3).format(* coord1[iatom]), file=ofile1)
    
92.             print(f" {elements[iatom]:<2s}" + ("    {:13.8f}" * 3).format(* coord2[iatom]), file=ofile2)
    
93. 
    
94. def align_molecules_quaternion(ifilename1: str, ifilename2: str):
    
95.     ifilename1_split = os.path.splitext(ifilename1)
    
96.     ifilename2_split = os.path.splitext(ifilename2)
    
97.     assert ifilename1_split[1] == ".xyz" and ifilename2_split[1] == ".xyz"
    
98.     elements = np.loadtxt(ifilename1, unpack=True, usecols=0, skiprows=2, dtype=np.dtype("<U2"))
    
99.     coord1 = np.loadtxt(ifilename1, usecols=(1, 2, 3), skiprows=2, dtype=np.double)
    
100.     coord2 = np.loadtxt(ifilename2, usecols=(1, 2, 3), skiprows=2, dtype=np.double)
     
101.     rotate_by_quaternion(coord1, coord2)
     
102.     ofilename1 = ifilename1_split[0] + "_aligned.xyz"
     
103.     ofilename2 = ifilename2_split[0] + "_aligned.xyz"
     
104.     with open(ofilename1, "w") as ofile1, open(ofilename2, "w") as ofile2:
     
105.         print(len(elements), file=ofile1)
     
106.         print(len(elements), file=ofile2)
     
107.         print(f"aligned from {ifilename1:s}", file=ofile1)
     
108.         print(f"aligned from {ifilename2:s}", file=ofile2)
     
109.         for iatom in range(len(elements)):
     
110.             print(f" {elements[iatom]:<2s}" + ("    {:13.8f}" * 3).format(* coord1[iatom]), file=ofile1)
     
111.             print(f" {elements[iatom]:<2s}" + ("    {:13.8f}" * 3).format(* coord2[iatom]), file=ofile2)
     
112. 
     
113. if __name__ == "__main__":
     
114.     ifilename1 = sys.argv[1]
     
115.     ifilename2 = sys.argv[2]
     
116.     # ifilename1 = "1.xyz"
     
117.     # ifilename2 = "2.xyz"
     
118.     # align_molecules(ifilename1, ifilename2, imethod=1, allow_reflection=True)
     
119.     align_molecules(ifilename1, ifilename2)
     
120.    
     

复制代码

gonglijing

谢谢，是没对齐，已按您们建议调整好了，再次感谢！