3. 轨迹分析基础：玩转五张数据画像
做分子动力学模拟，最怕的是什么？
跑了几百纳秒，拿到轨迹文件，却发现不同批次结果对不上。
这时候，轨迹分析就是你的第一道防线。
我个人习惯，拿到轨迹后先做五件事：RMSD、RMSF、回旋半径、氢键、二级结构。
这五张图看下来，模拟有没有出问题，心里就有底了。
3.1 RMSD：批次间的“基准尺”
RMSD，全称 Root Mean Square Deviation。
说白了，就是看某一帧结构相对于参考结构，原子位置偏移了多少。
数学上长这样：
RMSD(t) = sqrt( (1/N) * Σ|r_i(t) - r_i_ref|² )其中 N 是原子数，r_i(t) 是第 i 个原子在 t 时刻的坐标，r_i_ref 是参考坐标。
在批次检测中，RMSD 是最直观的指标。
两个批次如果体系一致，RMSD 的时间演化曲线应该形态相似、收敛值接近。
如果一个批次始终比另一个高 0.2 nm，那肯定有问题。
| 指标 | 目标 | 批次差异特征 |
|---|---|---|
| RMSD 均值 | 所有批次 < 0.15 nm（对蛋白而言） | 偏差 > 0.05 nm 需警惕 |
| RMSD 波动幅度 | 稳定后波动 < 0.02 nm | 某些批次持续震荡 ≠ 正常构象采样 |
| 平衡时间 | 所有批次在 2 ns 内进入平台 | 某批次迟迟不收敛 → 初始化可能出错 |
实用代码（MDAnalysis）：
import MDAnalysis as mda
from MDAnalysis.analysis import rms
u = mda.Universe('topol.tpr', 'traj.xtc')
ref = mda.Universe('topol.tpr')
R = rms.RMSD(u, ref, select='backbone')
R.run()
# R.rmsd 是一个数组，第 2 列是 RMSD 值
print(R.rmsd[:, 2])3.2 RMSF：揪出那些“不老实”的区域
RMSF 和 RMSD 是孪生兄弟。
RMSD 看的是时间轴上的整体漂移，RMSF 看的是每个残基在整个轨迹中的平均波动幅度。
公式几乎是翻版：
RMSF(i) = sqrt( (1/T) * Σ|x_i(t) - x_i_mean|² )这里是对时间平均，非参考结构。
我在项目中遇到过这么个情况：两个批次整体 RMSD 都很漂亮，0.12 nm 左右。
但一画 RMSF，发现第一批次的 loop 区波动是第二批次的两倍。
查到最后，发现第一批次里某个二硫键没加上。这个问题单看 RMSD 根本看不出来。
- 如果你发现某个残基 RMSF > 0.3 nm，标注出来。
- 如果这种高波动残基在批次间不一致，那模拟可能有问题。
避坑经验：我推荐用 Cα 原子算 RMSF，而不是所有重原子。
侧链波动太大，容易掩盖骨架上的真实差异。
3.3 回旋半径：判断蛋白是“胖了”还是“瘦了”
回旋半径（Rg）衡量的是蛋白的紧凑程度。
值越大，结构越松散；值越小，折叠越紧。
怎么算？
Rg² = (Σm_i * |r_i - r_cm|²) / Σm_ir_cm 是质心，m_i 是原子质量。
我觉得 Rg 在批次差异检测里是个隐形的信号灯。
某个批次的 Rg 突然多了 0.1 nm，并且持续不回弹，多半是模拟跑飞了——
常见原因包括：力场参数不匹配、温度耦合出问题、或者某个约束没加上。
你可以用这个代码快速算：
import MDAnalysis as mda
from MDAnalysis.analysis import rdf
u = mda.Universe('topol.tpr', 'traj.xtc')
protein = u.select_atoms('protein')
rg_values = []
for ts in u.trajectory:
rg = protein.radius_of_gyration()
rg_values.append(rg)
print(f"平均 Rg: {np.mean(rg_values):.3f} nm")3.4 氢键数目：稳定性的微观晴雨表
氢键是维持蛋白二级结构和配体结合的核心力量。
我说句实话，在批次检测里，氢键数目变化往往比 RMSD 更敏感。
判断氢键有三个几何标准：
- 供体-受体距离 ≤ 0.35 nm
- 角度（供体-氢-受体）≥ 120° 或 150°（看定义）
- 持续时间：≥ 10% 的轨迹存在才算稳定氢键
| 场景 | 正常波动 | 需要怀疑 |
|---|---|---|
| 球状蛋白内部 | 80-120 个 | 某批次少 15 个以上 |
| 蛋白-配体界面 | 3-8 个 | 某个批次少 2 个以上 |
| 溶剂暴露的 α-螺旋 | 每个螺旋转角 1 个（i→i+4） | 断裂成 310-螺旋 → 可能力场不对 |
警告：不同工具对氢键的判定标准不一样。
当你比较批次时，请确保用同一个分析工具、同一套参数。
我见过有人用 GROMACS 算 Batch A、用 VMD 算 Batch B，结果差了 30 个——其实只是阈值没对齐。
3.5 二级结构变化：看“骨架”的稳定性
二级结构是蛋白的骨架。
α-螺旋、β-折叠、转角这些东西变了，通常意味着功能区有构象变化。
我常用的工具有两种：
- DSSP（Dictionary of Secondary Structure of Proteins）——算法重，但准确。
- STRIDE——速度更快，对大体系友好。
在 Python 里用 MDAnalysis 调 DSSP：
from MDAnalysis.analysis import dssp
u = mda.Universe('topol.tpr', 'traj.xtc')
D = dssp.DSSP(u)
D.run()
# D.results['structure'] 包含每条链每帧的二级结构代码
# 你可以统计每种结构的占比变化批次检测看什么？
- 关键功能区的二级结构在批次间是否一致。
- α-螺旋含量的平均值偏差 ≤ 5%。
- 如果发现某个批次比别的批次多了一段 β-折叠，要小心——可能是错误折叠。
整体流程一张图
下面这张图，是我自己归纳的轨迹分析流程。
从原始轨迹到五张数据画像，再到批次判断，一目了然。
实战小结
好，咱们做个总结。
- RMSD：看整体漂移，找批次间的收敛差异。
- RMSF：揪出高波动区域，排查局部异常。
- Rg：判断蛋白是否“松散”或“塌缩”。
- 氢键：微观稳定的晴雨表，比 RMSD 更敏感。
- 二级结构：检查骨架有没有发生不该有的变化。
别偷懒，每次做完模拟，先把这五张图画出来。
养成这个习惯，你能省下大量排查批次差异的时间。