使用Python numpy计算原子坐标RMSD结果有误，如何修正公式？

核心错误点 #

• 维度平均逻辑不符合标准RMSD公式：原子坐标数组通常形状为(N, 3)（N为原子数，3对应x/y/z三个空间维度），你直接全局调用mean()等价于将3个空间维度也纳入了平均统计，分母为N*3，而标准RMSD的分母仅为原子数N，二者数值完全不一致。
• 若你需要计算的是分子结构领域的常规RMSD，还缺少*刚性对齐（Kabsch对齐）*步骤：两套坐标没有经过平移、旋转做最优叠加，直接计算的坐标偏差不具备结构比较的意义。

修正代码 #

场景1：计算未对齐的原始坐标RMSD #

仅修正公式实现错误，不需要做结构叠加：

import numpy as np

def cal_rmsd_numpy(coord_1, coord_2):
    # 输入要求：coord_1、coord_2均为形状(N,3)的numpy数组，N为原子数
    diff = coord_1 - coord_2
    # 先对每个原子的三个维度差平方求和，再对原子数取平均后开根号
    rmsd = np.sqrt((diff ** 2).sum(axis=1).mean())
    return rmsd

rmsd = cal_rmsd_numpy(coord_1, coord_2)
print(rmsd)

场景2：计算标准结构RMSD（先做刚性对齐） #

适用于分子结构比较的常规场景，叠加对齐后再计算RMSD：

import numpy as np

def cal_aligned_rmsd(coord_1, coord_2):
    # 平移到质心重合
    centroid1 = coord_1.mean(axis=0)
    centroid2 = coord_2.mean(axis=0)
    coord1_center = coord_1 - centroid1
    coord2_center = coord_2 - centroid2

    # Kabsch算法求解最优旋转矩阵
    H = coord1_center.T @ coord2_center
    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
    coord1_aligned = coord1_center @ R

    # 计算对齐后的RMSD
    diff = coord1_aligned - coord2_center
    rmsd = np.sqrt((diff ** 2).sum(axis=1).mean())
    return rmsd

rmsd = cal_aligned_rmsd(coord_1, coord_2)
print(rmsd)

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