分子动力学模拟流程：从结构文件到轨迹分析

很多新手拿到一个蛋白结构，直接就扔进 GROMACS 开始跑模拟。结果跑了两天发现力场不匹配，或者体系没平衡好。嗯，我刚开始也是这样踩坑的。

分子动力学模拟看着步骤很多，但拆开来看，其实就是一条流水线。我把这套流程总结成七个核心阶段。你只要跟着走，基本不会出大问题。

第一步：结构准备与检查

不管你是从 PDB 数据库下载的，还是同源建模出来的，拿到手的结构第一件事——检查。我记得有一次，我从 PDB 下载了 1ake.pdb，直接跑 pdb2gmx 就报错。后来发现文件里居然有替代氨基酸残基，还有一个配体分子根本没在力场里。

我个人习惯，先做三件事：

1. 用可视化软件看一眼结构，有没有断裂、缺失侧链。
2. 检查有没有非标准残基或配体需要参数化。
3. 删除结晶水和非必要的配体。

然后是格式转换。GROMACS 虽然能读 PDB，但它更喜欢自家格式。

gmx pdb2gmx -f input.pdb -o processed.gro -p topol.top -ff amber99sb-ildn -water tip3p

这一步会生成三样东西：gro 坐标文件、top 拓扑文件、力场参数。别忘了，-ignh 参数可以忽略已有氢原子，让 GROMACS 自己加。否则可能出现氢原子类型不匹配的情况。

第二步：拓扑与力场构建

拓扑文件是整个模拟的「宪法」。它定义了原子类型、键合参数、电荷分布。说白了，你的蛋白质怎么动、能谱长什么样，全写在里面。

这里有一个坑——力场选择。千万别混用。你选了 amber99sb-ildn，那后续所有参数都要用同一套。我在一个项目里见过有人把 charmm27 的配体参数放进 amber 的体系里，结果能量直接发散。

# 在 topol.top 中引入配体参数
#include "ligand.itp"

[ molecules ]
; Compound        #mols
Protein_chain_A     1
LIG                 1
SOL             10000
Na+                 5
Cl-                 5

第三步：溶剂化与添加离子

蛋白在晶体里和在水里的行为完全不同。所以我们要把蛋白扔进一个「水盒子」里。

盒子形状怎么选？常见的有立方体、长方体、截断八面体。我个人偏爱截断八面体，因为它能用更少的水分子达到同等溶质溶剂比，跑起来快不少。

gmx editconf -f processed.gro -o box.gro -c -d 1.0 -bt dodecahedron

加水的命令很简单：

gmx solvate -cp box.gro -cs spc216.gro -o solv.gro -p topol.top

然后加离子。为什么要加？因为真实体系不是纯水。你需要中和电荷，还要模拟生理盐浓度（通常 0.15 M NaCl）。

gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral -conc 0.15

第四步：能量最小化

你在构建结构时，可能会引入原子重叠或键角不合理。直接跑 MD 的话，这些「坏接触」会让能量瞬间暴涨，模拟直接崩溃。

能量最小化的目的，就是找到一个局部能量最低点。常用两种方法：最陡下降法（steep）和共轭梯度法（cg）。

我习惯先用最陡下降法跑 5000 步。因为它在远离平衡点的时候收敛快。然后换共轭梯度法跑 5000 步收尾，精细调整。

gmx grompp -f em_steep.mdp -c solv_ions.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em

判断是否收敛：看输出的能量曲线。如果最后几百步的势能变化小于 1000 kJ/mol/nm，基本就合格了。我最关心的一个值是最大受力（Fmax），它应该小于 1000 kJ/mol/nm。如果还在 10000 以上，说明有严重的冲突，需要检查结构。

第五步：NVT 平衡

能量最小化只是热身。接下来要把体系从 0 K 加热到目标温度（通常是 300 K）。

NVT 系综下，粒子数、体积、温度是常数。我们控制温度，让分子获得动能。这个过程也叫做「加热」。温度耦合一般用 V-rescale 算法，稳定且不会引入太大扰动。

gmx grompp -f nvt.mdp -c em.gro -p topol.top -o nvt.tpr
gmx mdrun -v -deffnm nvt

跑完之后，务必检查温度曲线。有个快速验证方法：

gmx energy -f nvt.edr -o temperature.xvg
# 选择 Temperature

你看下温度是否在 300 K 附近震荡。如果偏差很大，比如一直停在 280 K，说明耦合参数或者初速度设置有问题。

第六步：NPT 平衡

NVT 阶段温度稳定了，但压力还不一定对。NPT 平衡要调节盒子的尺寸，让体系达到合适的密度。说白了，就是给体系「压缩」或「膨胀」到正确的压强。

这一步我最常犯的错是什么？忘记写入 compressibility。水在 300 K 下的压缩系数是 4.5e-5 bar^-1。你要是设成 0，盒子根本不变化，NPT 等于白跑。

; npt.mdp 关键参数
Pcoupl              = Parrinello-Rahman
Pcoupltype          = isotropic
tau_p               = 2.0
compressibility     = 4.5e-5
ref_p               = 1.0

跑完后检查密度：

gmx energy -f npt.edr -o density.xvg
# 选择 Density

水的密度在 300 K 下约 996 kg/m³。如果你看到 980 以下或 1020 以上，说明盒子尺寸不对或者力场参数有问题。

第七步：成品 MD 生产

平衡没问题了，终于可以跑成品模拟了。这一步就是把 NPT 最后一步的构象作为输入，接着跑更长的时间。

时间步长怎么选？最常见的是 2 fs。但如果你的体系里有氢原子运动很快，可以试试 4 fs 加虚拟氢原子约束。跑得快一倍，但前提是你得用 hydrogen mass repartitioning（HMR）技术。

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tpr
gmx mdrun -v -deffnm md_0_1 -nb gpu -pme gpu

注意这里我用了 -t 参数读入 NPT 的 checkpoint 文件。这样速度、温度、压力都是连续的，不会出现跳变。

轨迹后处理与分析

跑完 MD，你手里有一堆 .xtc 轨迹文件和 .edr 能量文件。怎么分析？

最基本的三板斧：

举个例子，RMSD 分析：

gmx rms -s md_0_1.tpr -f md_0_1.xtc -o rmsd.xvg
# 选择骨架为拟合组和计算组

RMSD 曲线如果一直在涨，没有平台，说明体系还没平衡。这时候你要回头检查 NPT 平衡是否足够长，或者是不是蛋白本身就不稳定。

我个人习惯，分析之前先做一步「修正」。轨迹里可能有周期性边界跳跃，不做修正的话，RMSD 会假性变大。

gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_center.xtc -center -pbc mol -ur compact

这一步把蛋白放到盒子中心，消除周期跳跃。你后面算任何结构相关量，都会准很多。