分子动力学（MD）模拟的门槛很有意思：命令其实就那几条，跑通不难；难的是跑出来的结果可不可信。 这篇文章把我整理的一套标准流程记下来，从 PDB 结构一直到 RMSD/RMSF 分析，顺带列几个几乎每个新手都会踩的坑。
本文以 GROMACS 2020+ 版本为例，力场用 amber99sb-ildn，水模型用 TIP3P。不同版本个别参数名有差异，但主干流程一致。命令行细节以 GROMACS 官方手册 为准。
一、整体流程
一套完整的模拟大致分 8 步，前 4 步是体系搭建，第 5~7 步是平衡，最后一步才是真正的"生产跑"：
PDB 结构
↓ pdb2gmx 生成拓扑 + 坐标
↓ editconf 定义模拟盒子
↓ solvate 加水
↓ genion 加离子中和电荷
↓ 能量最小化 消除不合理接触
↓ NVT 平衡 稳定温度
↓ NPT 平衡 稳定压力/密度
↓ 正式 MD 采样
↓ 分析 RMSD / RMSF / Rg ...
新手最容易犯的错误就是跳过平衡直接跑 MD。 体系里原子坐标来自晶体结构，加水之后必然存在局部过近接触，不最小化直接跑，轻则 LINCS 报错，重则体系直接"炸开"。
二、体系搭建
1. 生成拓扑
gmx pdb2gmx -f protein.pdb -o conf.gro -water tip3p -ff amber99sb-ildn
这一步会让你选力场，并输出 topol.top（拓扑）和 conf.gro（坐标）。
常见报错与处理：
| 报错 | 原因 | 处理 |
|---|---|---|
Residue XXX not found in residue topology database |
含配体、离子或非标准残基 | 拆分出配体单独处理，或用 ACPYPE 等工具生成参数 |
Atom XXX is missing |
PDB 缺原子（尤其是氢） | 用 pdbfixer 补全，或 -ignh 让 GROMACS 自己加氢 |
| 二硫键没连上 | 未指定 | 加 -ss 让程序自动检测，或用 -cys 手动指定 |
PDB 文件里的水分子（HOH）要先删掉，否则会被当成残基处理而报错。
2. 定义盒子
gmx editconf -f conf.gro -o box.gro -c -d 1.0 -bt cubic
-d 1.0 意思是蛋白距盒壁至少 1.0 nm。这个数字不要省——盒子太小，蛋白会和自己的周期性镜像发生相互作用，模拟结果就废了。一般认为 1.0~1.5 nm 比较安全。
想省计算量可以用 -bt dodecahedron（菱形十二面体），体积比立方体小约 29%，对球状蛋白更划算。
3. 加水
gmx solvate -cp box.gro -cs spc216.gro -o solvated.gro -p topol.top
注意 -cs spc216.gro 虽然是 SPC 水盒子，但 GROMACS 会根据拓扑里的水类型做替换，配合 TIP3P 用也没问题。
4. 加离子
先准备一个 ions.mdp，然后：
gmx grompp -f ions.mdp -c solvated.gro -p topol.top -o ions.tpr
gmx genion -s ions.tpr -o neutral.gro -p topol.top -pname NA -nname CL -neutral
-neutral 会自动中和体系净电荷。执行时会让你选替换对象，选 SOL。
体系净电荷不为零会导致 PME（粒子网格 Ewald）计算出错，这一步不能跳。想模拟生理盐水浓度则用 -conc 0.15。
三、能量最小化
gmx grompp -f em.mdp -c neutral.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em
em.mdp 关键参数：
integrator = steep ; 最速下降法
emtol = 1000.0 ; 收敛判据 (kJ/mol/nm)
emstep = 0.01
nsteps = 50000
判据看势能是否收敛到 emtol 以下。如果提示未收敛，先把 nsteps 加大；还不行就说明体系里有严重冲突，得回头检查结构。
四、两步平衡
NVT（恒定粒子数、体积、温度）
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
gmx mdrun -deffnm nvt
关键参数：
integrator = md
dt = 0.002
nsteps = 50000 ; 100 ps
tcoupl = V-rescale ; 推荐，不要用 Berendsen
tc-grps = Protein Non-Protein
tau_t = 0.1
ref_t = 300
gen_vel = yes ; 生成初始速度
温度耦合建议用 V-rescale（Bussi-Parrinello 随机速度重标定）。老教程里的 Berendsen 恒温器不能产生正确的正则系综，GROMACS 官方现已不推荐用于生产模拟。
NPT（恒定粒子数、压力、温度）
gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -r nvt.gro -p topol.top -o npt.tpr
gmx mdrun -deffnm npt
压力耦合（Parrinello-Rahman 原始方法，或其改进版 C-rescale）：
pcoupl = C-rescale ; 或 Parrinello-Rahman
pcoupltype = Isotropic
tau_p = 2.0
ref_p = 1.0
compressibility = 4.5e-5
平衡是否成功，看密度是否稳定在 ~1000 kg/m³、温度是否稳定在设定值。
五、正式模拟
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md.tpr
gmx mdrun -deffnm md -ntmpi 1 -ntomp 8
有 GPU 的话加 -nb gpu，速度能上一个数量级。线程参数 -ntmpi（MPI 进程数）和 -ntomp（OpenMP 线程数）要根据机器核数调，不是越多越好。
六、分析
RMSD（整体构象漂移）
gmx rms -s md.tpr -f md.xtc -o rmsd.xvg -tu ns
选 Backbone 做最小二乘拟合组，再选 Backbone 计算 RMSD。曲线先升后平说明体系达到平衡——只有平衡之后的那段轨迹才应该用于分析。
RMSF（残基柔性）
gmx rmsf -s md.tpr -f md.xtc -o rmsf.xvg -res
高峰对应柔性区域（loop、端基），低谷常是二级结构或结合位点。
回转半径 Rg
gmx gyrate -s md.tpr -f md.xtc -o gyrate.xvg
Rg 突然变大往往意味着蛋白去折叠，是个很有用的"健康指标"。
画图
xmgrace rmsd.xvg
服务器无图形界面的话，导出数据用 Python/matplotlib 画更省事。
七、几个务实的提醒
- 先跑短、再跑长。 用 100 ps 试跑验证设置无误，再提交 100 ns 的长任务。否则跑了三天才发现参数错了。
- 轨迹文件很占空间。
xtc是压缩格式，别用trr；nstxout-compressed设为 5000（每 10 ps 存一帧）通常够用。 - PBC 问题。 分析前用
gmx trjconv -pbc mol -center处理轨迹，否则蛋白跨盒子边界时看起来像"裂成两半"。 - 结果要重复。 单次模拟说明不了什么，重要结论至少做 2~3 个独立重复。
小结
GROMACS 的门槛不在命令，在于对每一步为什么这么做有判断。同样的命令，参数选择不同，跑出来的东西可信度可以差一个数量级。把流程记熟之后，真正花时间的是设计模拟和分析结果——那才是计算化学本身。
愿你的体系都不炸，LINCS 都不报错。
参考文献
- Lindorff-Larsen, K., et al. (2010). Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins, 78(8), 1950–1958. https://doi.org/10.1002/prot.22711
- Jorgensen, W. L., et al. (1983). Comparison of simple potential functions for simulating liquid water. The Journal of Chemical Physics, 79(2), 926–935. https://doi.org/10.1063/1.445869
- Bussi, G., Donadio, D., & Parrinello, M. (2007). Canonical sampling through velocity rescaling. The Journal of Chemical Physics, 126(1), 014101. https://doi.org/10.1063/1.2408420
- Berendsen, H. J. C., et al. (1984). Molecular dynamics with coupling to an external bath. The Journal of Chemical Physics, 81(8), 3684–3690. https://doi.org/10.1063/1.448118
- Parrinello, M., & Rahman, A. (1981). Polymorphic transitions in single crystals: A new molecular dynamics method. Journal of Applied Physics, 52(12), 7182–7190. https://doi.org/10.1063/1.328693
- Darden, T., York, D., & Pedersen, L. (1993). Particle mesh Ewald: An N·log(N) method for Ewald sums in large systems. The Journal of Chemical Physics, 98(12), 10089–10092. https://doi.org/10.1063/1.464397
- GROMACS Development Team. GROMACS Manual. https://manual.gromacs.org/