RMSD 计算原理：坐标叠加、最小二乘法、旋转矩阵
聊 RMSD 之前，我想先问一个问题：两个结构的坐标，凭什么能直接相减？
你想想，如果两个蛋白骨架一个朝东、一个朝西，直接算原子距离那肯定大得离谱。实际测量的是“形状差异”，而不是“位置差异”。所以第一步——坐标叠加，就是把两个结构挪到同一个坐标系下，让它们尽可能重合。
坐标叠加：不只是平移
最早我入行时，以为叠加就是算个质心，平移过去完事。结果算出来的 RMSD 忽高忽低，直到前辈点醒我：旋转才是关键。
坐标叠加包含两步：
- 平移：把两个结构的质心对齐。
- 旋转：找到最优的旋转角度，让对应原子距离平方和最小。
平移简单，质心一减就行。但旋转就涉及我们下面要说的核心——最小二乘法。
最小二乘法：目标函数长什么样
我们有一组参考点集 X（比如实验结构坐标）和一组待对齐点集 Y（比如模型坐标）。想找一个旋转矩阵 R 和平移向量 t，使得：
min Σ || X_i - (R * Y_i + t) ||² (i = 1..N)说白了，就是让每个对应原子之间的距离平方和最小。为什么用平方？因为我做过一次绝对值最小化的尝试，结果优化时容易震荡，平方项能保证光滑可导，求解方便。
数学上，这个最小化问题有解析解。经典方法就是 Kabsch 算法，它基于奇异值分解（SVD）。
旋转矩阵求解：Kabsch 算法三步走
我个人习惯把 Kabsch 算法拆成三步，理解起来很顺：
- 去质心：构造去质心后的坐标矩阵 P 和 Q。
- 计算协方差矩阵：H = Pᵀ · Q。
- SVD 分解：H = U · Σ · Vᵀ，然后旋转矩阵 R = V · Uᵀ。
注意：如果行列式 det(R) < 0，要取 V 的最后一列取反来保证右手坐标系。
import numpy as np
def kabsch(P, Q):
# P, Q: Nx3 numpy arrays
centroid_P = np.mean(P, axis=0)
centroid_Q = np.mean(Q, axis=0)
P_centered = P - centroid_P
Q_centered = Q - centroid_Q
H = P_centered.T @ Q_centered
U, S, Vt = np.linalg.svd(H)
R = Vt.T @ U.T
if np.linalg.det(R) < 0:
Vt[-1, :] *= -1
R = Vt.T @ U.T
t = centroid_P - R @ centroid_Q
return R, t这段代码我反复用过几十次，看起来简洁，但 行列式校正那一步很容易忘。我第一次写时就漏了，结果算出来的 RMSD 偏大，还找半天 bug。
RMSD 最终计算
得到最优旋转矩阵 R 和平移 t 后，叠加好坐标，然后：
RMSD = sqrt( (1/N) * Σ || X_i - (R * Y_i + t) ||² )单位是 Å（埃）。我一般保留到小数点后四位，但可视化时只显示两位就够。
一张图理解整个流程
下面是我画的一个流程示意，从原始坐标到最终 RMSD。
这条路径走完，RMSD 就出来了。我个人的习惯是把 SVD 结果也打印出来看看奇异值大小，如果最后一个奇异值接近 0，说明两个点集几乎共面，旋转自由度受限，这时候需要小心处理。
几种旋转表示方式的对比
| 表示方式 | 优点 | 缺点 | 适用场景 |
|---|---|---|---|
| 欧拉角 | 直观，三个角度 | 万向锁，优化易局部极值 | 简单演示，非迭代 |
| 四元数 | 无万向锁，可迭代优化 | 不直观，需要归一化 | 分子动力学模拟 |
| 旋转矩阵 + SVD | 封闭解，不用迭代 | 必须保证正交性 | Kabsch / RMSD 标准方法 |
我比较推荐直接用 SVD 的方案，因为你不用操心初始值，一次到位。四元数也有它的好处，但在 RMSD 场景下 SVD 更成熟。
避坑指南
我曾经因为忘记滤除氢原子，导致 RMSD 计算特别慢。加氢后原子数翻倍，但大多结构比较随意，h 噪音太大。后来我养成了预处理先剔除氢原子的习惯。
另外，坐标单位要一致。有一次我从 PDB 读取坐标是 0.1 Å 的精度输入给代码，但另一个结构是 1.0 Å，结果算出来的 RMSD 直接差了 10 倍。真是低级错误。
小结
RMSD 计算核心其实就是三个东西：坐标叠加消除刚体运动，最小二乘法定义优化目标，旋转矩阵通过 SVD 直接求解。把这个流程刻在脑子里，后续画曲线、做分析就顺了。