5.2 RMSD 与 RMSF 判读

5.2 RMSD 与 RMSF 判读

本节摘要：RMSD（均方根偏差）回答"体系离参考结构整体走多远了"，RMSF（均方根涨落）回答"每个原子自己晃得多厉害"——前者是时间的函数、后者是原子的函数。本节从 Kabsch 叠合算法出发手写一个最小 RMSD 实现，验证 gmx rms 的口径；再讲三段式曲线判读、参考结构与叠合组的选择陷阱。这两个量是几乎所有分析报告的第一张图，也是平衡判据在结构层面的延伸。

5.1 修好了几何，本节问第一个问题：体系稳了吗。承上：3.6 的判据盯的是热力学量（温度、密度），它们稳了只说明宏观状态到位，结构层面的弛豫（聚合物链从初始构象松弛到平衡构象库）要靠本节的量来看。启下：RMSF 找出的高柔性部位，正是 5.3 局域结构分析该细看的区域，也是 6.2 自由能面里该重点采样的构象方向。

手写 RMSD：Kabsch 算法五步

RMSD 的定义：对叠合后的两组坐标，逐原子算距离平方、求平均、开根号：

RMSD(t) = √[ (1/N) Σ_i | r_i(t) − r_i^ref |² ]

定义里的"叠合"是全部技术含量所在：把 t 时刻的结构整体平移加旋转，与参考结构的偏差最小化（消除刚体运动的贡献，只留内部形变）。最优旋转由 Kabsch 算法给出，五步：两组坐标各自去质心；算 3×3 协方差矩阵 R = PᵀQ；对 R 做奇异值分解；修正反射（最后一列奇异向量乘以 det 符号）；最优旋转 = U·Vᵀ。Python 实现不到三十行：

import numpy as np

def kabsch_rmsd(P, Q):
    """P, Q: N×3 坐标数组（原子一一对应）。返回叠合后的 RMSD（nm）。"""
    Pc = P - P.mean(axis=0)          # 各自去质心
    Qc = Q - Q.mean(axis=0)
    R = Pc.T @ Qc                    # 协方差矩阵
    U, S, Vt = np.linalg.svd(R)
    d = np.sign(np.linalg.det(U @ Vt))   # 防反射：镜像不算旋转
    D = np.diag([1.0, 1.0, d])
    Rot = U @ D @ Vt
    P_fit = Pc @ Rot
    return np.sqrt(((P_fit - Qc) ** 2).sum() / len(P)), P_fit + Q.mean(axis=0)

# 自测：造一条"弯了一下的五珠链"，先与自身比（应为 0），再与微扰版比
rng = np.random.default_rng(42)
chain = np.cumsum(rng.normal(size=(5, 3)), axis=0) * 0.15
r0, _ = kabsch_rmsd(chain, chain.copy())
perturbed = chain + rng.normal(scale=0.02, size=chain.shape)
r1, _ = kabsch_rmsd(perturbed, chain)
print(f"自身 RMSD = {r0:.2e}（应为 0）")
print(f"扰动 RMSD = {r1:.4f} nm（量级应接近扰动幅度）")
​

最后两行就是自测口径：自身对比必须精确为零（算法对刚体运动不敏感的证明），扰动对比的量级必须与扰动幅度一致。把这个实现与 gmx rms 对同一条轨迹的结果对照，数值一致即口径验证通过——工具给你的每个数都值得这样抽验一次，此后才值得信任。

gmx 的实操与曲线判读

实操上两步走：先做叠合（把轨迹对齐到参考结构），再算 RMSD：

# 第一步：以骨架重原子为叠合组，把整条轨迹旋转平移对齐到 t0 结构
$ printf "Backbone\nBackbone\n" | gmx rms -s npt.tpr -f fix3.xtc -fit rot+trans \
        -o rmsd_bb.xvg
# 第二步：逐残基 RMSF（叠合同样重要，柔性分布才可比）
$ printf "PEO\n" | gmx rmsf -s npt.tpr -f fix3.xtc -fit -res -o rmsf.xvg
​

曲线判读看三段式。上升段（前几十 ns）：体系从初始构象弛豫走，RMSD 爬升——这段属于"还在平衡"，不进统计。平台段：RMSD 围绕某个均值涨落，幅度不随时间增长——这才是可采样的平衡区，平台"高度"本身没有绝对意义（它依赖参考结构与叠合组），有意义的是"是否平台"。二次爬升或台阶：平台后又上一个台阶，通常是构象转变（聚合物的链解缠结事件、脂膜的相变）或退化（聚集、塌缩）——前者是重大发现，后者要回 4.3 排查。

图 5-2 RMSD 曲线三段式判读

RMSF：柔性的地址簿

RMSF 把同一份信息换成"按原子索引"的视角：先让轨迹对参考结构叠合（去刚体运动），再对每个原子统计其位置的时间涨落。它回答的问题与 RMSD 严格不同——RMSD 是"整体走了多远"（时间轴），RMSF 是"谁在动"（原子轴）。聚合物的用法：RMSF 沿链号的曲线能直接读出端基比中段活泼（端基峰）、缠结点附近被压平（局部谷）；脂膜的用法：逐叶统计磷氮原子 RMSF，相变温度以上显著增大，是 7.2 判流动态的辅助指标。

RMSF 有一个必须警惕的解读陷阱：高 RMSF 不等于"重要"。端基与柔性侧链的 RMSF 天然大，但它们对力学性质的贡献未必大；反过来，RMSF 很小的部位（紧密堆积的核心区）哪怕动 0.01 nm 也可能是关键结构支撑。RMSF 只告诉你"谁在动、动多大"，"动了要紧吗"要靠 5.3 的结构量与 6.1 的自由能来回答。

演算：平台段有多"平"？用块平均量化

"看起来平了"同样要数字化。把平台段切成若干等长块，看块均值的标准误是否不再随块长增长——这是轨迹分析通用的收敛诊断（误差棒也靠它）：

import numpy as np

def block_error(y, nblocks=8):
    """把序列分成 nblocks 块，返回块均值的标准误（含块长相关的有效样本估计）。"""
    blocks = np.array_split(y, nblocks)
    means = np.array([b.mean() for b in blocks])
    return means.mean(), means.std(ddof=1) / np.sqrt(nblocks)

rng = np.random.default_rng(5)
n = 4000
# 模拟 RMSD 平台段：均值 0.31 nm 的强相关噪声（自相关时间约 200 帧）
noise = np.empty(n)
noise[0] = 0.31
for i in range(1, n):
    noise[i] = 0.995 * noise[i-1] + 0.31 * 0.005 + rng.normal(scale=0.004)

for nb in (2, 4, 8, 16):
    m, e = block_error(noise, nb)
    print(f"{nb:2d} 块: 均值 {m:.4f} nm, 标准误 {e:.4f}")

# 对比：把同样的数据当独立样本算误差，会低估几倍
print(f"天真标准差/根号N: {noise.std(ddof=1)/np.sqrt(n):.4f}（远小于块估计 → 低估证据）")
​

判读方法：块数从少到多，块均值的标准误应趋于稳定；若标准误随块数持续上升，说明相关时间长于当前块长，误差还没收敛——采样再加长。那个"天真误差"对比行就是审稿人抓包"误差棒过小"的计算依据。

两个高频疑问

RMSD 多小算"稳定"？ 这个问题本身是错的——RMSD 数值依赖参考结构、叠合组与体系柔度，聚合物熔体的平台 1 nm 也可能是健康的，蛋白核心区 0.15 nm 也可能藏着相变。正确的问法是"RMSD 是否进入有界平台、块平均误差是否收敛"（本节演算的两个量），数值高低只在与同类体系同口径对比时才有意义。

要不要去掉整体平动转动再算 RMSF？ 要，gmx rmsf 的 -fit 默认做。但注意一个组合陷阱：先对轨迹做过 -fit rot+trans 的 RMSD（5.2 的第一步），又让 rmsf 再叠合一次，两次叠合的参考不一致会让柔性分布整体平移。规矩是同一条分析链路里叠合参考与叠合组保持唯一，分析脚本开头注释写明。

易错点与本节速记

四个高频失误：叠合组选了全部原子（溶剂噪声淹没信号——叠合组选骨架或等效代表组，统计组另选）；参考结构选了 t=0 的刚最小化构象（平台值虚高且不可与文献比——文献惯例是用平衡段中某帧，且必须注明）；对未修 PBC 的轨迹直接算（5.1 的病，断裂分子让 RMSD 出现周期性尖刺）；对多分子体系问"体系的 RMSD"而其实想问"每条链的"（聚合物应逐链计算再统计分布，一条平均曲线会把解缠结事件抹平）。

速记：RMSD 是时间的函数、RMSF 是原子的函数；Kabsch 五步——去质心、协方差、SVD、防反射、旋转；曲线三段式——上升段别采样、平台段可采样、二次爬升要分辨是发现还是事故；叠合组与统计组分开选；高 RMSF 只说明会动，不说明重要。整体稳定性看完，下一节深入局域结构。