3. GROMACS简介：分子动力学模拟基础、常用命令、输入文件归纳总结
3.1 分子动力学模拟到底在干什么？
很多人第一次接触分子动力学（MD）时，容易把它想得太玄。其实说白了，就是让原子按照牛顿力学在计算机里“动起来”。你给定一个初始结构，给每个原子分配初速度，再定义好它们之间的作用力（力场），然后牛顿第二定律一步步积分——就这么简单。
我刚开始做MD的时候，总觉得这是“分子动画”，后来才意识到：只要力场够准、时间够长，模拟结果可以跟实验数据对上。我以前做过一个蛋白-配体结合自由能计算，模拟跑了两周，算出的ΔG跟实验值就差0.3 kcal/mol，那种感觉确实上瘾。
GROMACS是目前学术界用得最广的MD引擎之一。它最重要的特性是并行效率高、力场支持全面、而且完全免费。我个人的习惯是，每当新项目立项，先用GROMACS跑一轮快速验证，再说要不要上更复杂的工具。
3.2 GROMACS的典型工作流
一个标准的GROMACS模拟分为四步：
- 预处理 — 把PDB/结构文件转成GROMACS自己的格式，加上拓扑和力场参数。
- 能量最小化 — 消除初始结构中的原子冲突，避免模拟炸掉。
- 平衡 — 先NVT（控温），再NPT（控温控压），让系统稳定在目标条件。
- 成品模拟与产出分析 — 跑长时间的MD轨迹，然后算RMSD、Rg、氢键、自由能等等。
每个步骤在命令行里就是一条 gmx 命令加各种选项。别小看这些命令，参数写不对，跑一周可能全是废数据。我就曾经因为把压力耦合类型设错，整个NPT阶段跑偏了2小时才发现，最后只能重跑。
3.3 常用命令速查
GROMACS 主程序是 gmx，它的子命令特别多。我把最常用的几个列了个表，方便你快速检索：
| 命令 | 功能 | 典型用法示例 |
|---|---|---|
gmx pdb2gmx |
将PDB转为GROMACS格式并生成拓扑 | gmx pdb2gmx -f protein.pdb -o protein.gro -water spce |
gmx editconf |
修改盒子边界、转格式 | gmx editconf -f protein.gro -o boxed.gro -c -d 1.0 -bt cubic |
gmx solvate |
给体系加溶剂 | gmx solvate -cp boxed.gro -cs spc216.gro -o solv.gro -p topol.top |
gmx grompp |
将拓扑、坐标、参数打包成 .tpr 文件 | gmx grompp -f em.mdp -c solv.gro -p topol.top -o em.tpr |
gmx mdrun |
执行分子动力学运行（EM/MD均用） | gmx mdrun -v -deffnm em |
gmx energy |
提取能量、温度、压力等随时间变化 | gmx energy -f em.edr -o potential.xvg |
gmx traj |
处理轨迹文件（转换、截取、居中） | gmx traj -f md.xtc -s md.tpr -oxt filtered.xtc -b 1000 |
举个例子——能量最小化最常用的组合命令：
# 先准备 .tpr
gmx grompp -f em.mdp -c solv.gro -p topol.top -o em.tpr
# 再运行
gmx mdrun -v -deffnm em
# 检查能量是否收敛
gmx energy -f em.edr -o potential.xvg
💡 我的实战习惯： 每个步骤都用 -deffnm 指定输出文件名前缀，避免文件互相覆盖。另外 -v 加 verbose 输出，方便在日志里看收敛情况。
3.4 输入文件全景图
GROMACS 的输入文件有三种核心类型，我把它们整理成表格：
| 文件类型 | 后缀 | 内容 | 典型来源 |
|---|---|---|---|
| 结构/坐标文件 | .gro / .pdb |
原子坐标、盒子大小 | gmx pdb2gmx 或实验结构 |
| 拓扑文件 | .top / .itp |
原子类型、键、角、二面角、非键参数 | gmx pdb2gmx 自动生成 |
| 参数文件 | .mdp |
模拟控制参数（步长、温度、压力、约束等） | 用户手动编写 |
| 二进制结构+参数包 | .tpr |
经过 grompp 预处理后的完整输入 |
gmx grompp |
其中 .mdp 文件最容易被忽视，但也是最容易出问题的。我放一个典型的能量最小化 em.mdp 给大家参考：
integrator = steep
emtol = 1000.0
emstep = 0.01
nsteps = 5000
nstenergy = 10
cutoff-scheme = Verlet
ns_type = grid
rlist = 1.0
rcoulomb = 1.0
rvdw = 1.0
⚠️ 我曾经犯过的错： 把 integrator 写成了 md 而不是 steep，结果跑了一次“能量最小化”版的MD，步长太大直接崩了。后来查日志才发现。
3.5 新手常踩的坑
我从自己的经验出发，总结几个最常见的坑：
- 忘记加溶剂 — 在真空中跑蛋白，静电没屏蔽，体系直接发散。做完
pdb2gmx后一定记得solvate和genion。 - 拓扑与坐标不一致 — 你用
pdb2gmx生成的.gro原子顺序必须和.top完全对应，手改坐标文件会导致灾难。建议不要手动编辑.top里的原子列表。 - 平衡时间太短 — 我以前只给NVT跑了20 ps，结果NPT阶段体系压力一直震荡。后来发现控温约束还没稳定。一般至少100 ps。
- 并行线程设置不当 —
mdrun -ntmpi 4 -ntomp 2要根据物理核数来配，不然性能反而下降。建议用gmx mdrun -nt自动检测。
嗯，暂时讲到这里。GROMACS 的命令和文件体系虽然看起来多，但只要理解了一个典型工作流里的每一步做什么，剩下的就是反复练习。我刚开始的时候也是对着手册一条条敲，敲多了肌肉记忆自然就出来了。