RMSD 与 RMSF：MD 轨迹稳定性如何判断

RMSD（均方根偏差）与 RMSF（均方根涨落）是 MD 轨迹分析的第一组指标。它们计算简单，但误读的情况非常普遍——尤其是「RMSD 平了就是收敛了」这个说法。

两者的区别

计算：关键在于对齐

import MDAnalysis as mda
from MDAnalysis.analysis import rms, align
import numpy as np

u = mda.Universe("topology.pdb", "trajectory.dcd")

# ---- 蛋白 RMSD ----
R = rms.RMSD(u, u,
             select="backbone",              # 用于【对齐】的原子
             groupselections=["protein and name CA",
                              "resname LIG"],  # 额外计算 RMSD 的组
             ref_frame=0)
R.run()

# R.results.rmsd 的列：
#   [frame, time(ps), RMSD(对齐组), RMSD(组1), RMSD(组2), ...]
time = R.results.rmsd[:, 1] / 1000        # ns
rmsd_protein = R.results.rmsd[:, 2]
rmsd_ca = R.results.rmsd[:, 3]
rmsd_ligand = R.results.rmsd[:, 4]        # 【配体 RMSD，对齐蛋白后】

# 【最关键的一点】：
#   分析配体姿势稳定性时，必须先【对齐蛋白】，
#   然后计算【配体】的 RMSD
#
#   如果对齐配体自身，得到的是配体的内部构象变化，
#   而不是「配体是否在口袋中移动」—— 完全不同的问题
#
#   groupselections 参数正是为此设计：
#   用 backbone 对齐，但报告 LIG 的 RMSD

# ---- RMSF ----
# RMSF 计算前必须先对齐整条轨迹
average = align.AverageStructure(u, u, select="protein and name CA",
                                 ref_frame=0).run()
aligner = align.AlignTraj(u, average.results.universe,
                          select="protein and name CA",
                          in_memory=True).run()

ca = u.select_atoms("protein and name CA")
R_f = rms.RMSF(ca).run()

for res, rmsf in zip(ca.resids, R_f.results.rmsf):
    if rmsf > 2.0:
        print(f"残基 {res}: RMSF = {rmsf:.2f} Å  （柔性区）")

RMSD 的判读

「RMSD 平了 = 收敛了」是错的

# 这是 MD 分析中最普遍的误解
#
# RMSD 平稳只说明：
#   体系在【当前所处的状态】附近波动
#
# 它【不能】说明：
#   - 采样已经充分
#   - 已经找到了全局最稳定的状态
#   - 结果可以外推
#
# 反例：
#   体系陷在一个局部能量极小值中，
#   RMSD 会非常平稳 —— 但这正是采样不足的表现
#
# 更可靠的收敛判断：
#
# 1) 【多个独立副本给出一致的结果】
#    这是最重要的判据
#    不同初速度出发的模拟，结论是否相同？
#
# 2) 分块平均
#    把轨迹分成前后两半，各自算目标量
#    差异大 = 未收敛
#
# 3) 看具体的物理量而非只看 RMSD
#    - 关心的相互作用占有率是否稳定（见 109）
#    - 关心的距离/角度分布是否稳定
#    - 二级结构含量是否稳定
#
# 4) 自相关时间
#    估计关心量的去相关时间，
#    模拟长度应远大于它

def block_average_check(values, n_blocks=5):
    """分块检查收敛性"""
    blocks = np.array_split(np.asarray(values), n_blocks)
    means = [b.mean() for b in blocks]
    print("各块均值:", [f"{m:.2f}" for m in means])
    print(f"块间标准差: {np.std(means):.3f}")
    # 块间标准差远小于整体波动 = 较好的收敛迹象
    return np.std(means)

RMSF 的判读与应用

• 典型模式：二级结构区域（α 螺旋、β 折叠）RMSF 低（< 1 Å）；loop 与末端 RMSF 高（> 2 Å）；
• 与 B 因子对比：RMSF 应该与晶体学 B 因子有相关性——如果模拟中某个区域异常柔性而晶体中很刚性，可能是模拟有问题（或晶体接触限制了运动）；
• 结合位点的 RMSF：如果结合位点残基 RMSF 很高，说明口袋在波动——这对刚性对接是个警示，可能需要集合对接；
• 配体结合的影响：比较 apo 与 holo 模拟的 RMSF，能看出配体结合稳定了哪些区域——这有时能揭示变构效应；
• 末端要排除：链末端 RMSF 总是很高，这是几何必然，不代表功能相关的柔性。

配体姿势稳定性的完整判断

# 只看 RMSD 不够，应该综合多个维度
#
# 1) 配体 RMSD（对齐蛋白后）
#    < 2 Å 稳定
#
# 2) 关键相互作用的占有率（见 109）
#    共晶中的关键氢键，在模拟中占有率多少？
#    < 50% 说明该相互作用不稳定
#
# 3) 配体质心到口袋中心的距离
#    是否逐渐远离？
#
# 4) 埋藏表面积随时间的变化
#    减小 = 配体在往外走
#
# 5) 配体内部构象的变化
#    对齐配体自身算 RMSD
#    → 大幅变化说明初始构象是高能构象
#
# 6) 多副本一致性
#    3 个独立模拟都稳定 → 可信
#    1 个稳定 2 个漂走 → 不可信
#
# 综合判断示例：
def assess_pose_stability(u, ligand_sel="resname LIG"):
    from MDAnalysis.analysis import rms
    R = rms.RMSD(u, u, select="backbone",
                 groupselections=[ligand_sel], ref_frame=0).run()
    lig_rmsd = R.results.rmsd[:, 3]
    # 只看后半段（前半段可能还在平衡）
    tail = lig_rmsd[len(lig_rmsd)//2:]
    return {
        "mean": float(tail.mean()),
        "max": float(tail.max()),
        "drift": float(tail[-10:].mean() - tail[:10].mean()),
        "verdict": "稳定" if tail.mean() < 2.0 and tail.max() < 3.5 else "不稳定",
    }

延伸资源

• MD 入门：123《分子动力学 MD 入门》；MDAnalysis：129《MDAnalysis 分析轨迹》；OpenMM：128《OpenMM 跑 10–100 ns MD》；
• 相互作用占有率：109《ProLIF》；姿势验证：096《Binding Pose》。