第3章 分子动力学模拟基本流程与文件类型
聊起GROMACS，很多人一上来就盯着跑模拟的命令。其实，理解整个流程和文件流转，才是真正高效工作的关键。
我这几年带过不少新手，最常见的坑就是：折腾半天跑出结果，才发现输入文件就错了。白白浪费计算资源。
所以这一章，我们先画清楚骨架——模拟到底分几步，每步产生什么文件。说白了，你得知道自己在做什么。
3.1 标准模拟工作流
一次完整的分子动力学模拟，大致能分成五个阶段。不是每个项目都必须走完所有阶段，但标准流程是这样的：
- 模型构建与准备 —— 拿到初始结构（PDB），加氢、生成拓扑
- 能量最小化 —— 消除不合理接触，释放空间位阻
- 平衡模拟 —— NVT + NPT，让体系达到目标温度和压力
- 生产模拟（MD Run） —— 正式采样，记录轨迹
- 结果分析 —— 计算均方根偏差、径向分布函数、自由能等
每个阶段都会产生特定类型的文件。我自己的习惯是，给每个项目建一个清晰的目录结构，把不同阶段的文件分开存放。不然到最后，你连哪个是平衡后的构型都分不清。
3.2 各阶段的核心文件清单
下表列出了GROMACS在每个阶段涉及的主要文件。嗯，这张表我建议你收藏一下，实际干活时用得到。
| 阶段 | 输入文件 | 输出文件 | 说明 |
|---|---|---|---|
| 模型构建 | protein.pdb, forcefield.itp | protein.gro, topol.top, posre.itp | 生成坐标和拓扑 |
| 能量最小化 | protein.gro, topol.top, em.mdp | em.gro, em.edr, em.log, em.trr | 检查势能是否收敛 |
| NVT/NPT平衡 | em.gro, topol.top, nvt.mdp | nvt.gro, nvt.edr, nvt.log, nvt.trr | 稳定温度/压力 |
| 生产模拟 | npt.gro, topol.top, md.mdp | md.gro, md.edr, md.log, md.xtc | 主要采样轨迹 |
| 结果分析 | md.xtc, md.gro, topol.top | rmsd.xvg, rdf.xvg, 等 | 各类分析输出 |
需要注意：有些文件是累进的——比如 topol.top 从头用到尾，而轨迹文件 .trr 和 .xtc 会在不同环节产生不同版本，千万别混用。
核心要点：同一分子体系在不同阶段会生成多套“.gro”文件（如 em.gro、nvt.gro、npt.gro、md.gro）。它们代表不同平衡阶段的原子坐标，文件中原子顺序必须保持一致，否则后续分析会出错。
3.3 关键文件类型详解
3.3.1 坐标文件：.gro 与 .pdb
.gro 是GROMACS原生坐标格式，紧凑高效。它不含元素信息，只有原子名、残基名、坐标和盒子尺寸。
MD simulation, t= 0.00000 step= 0
50SOL OW 1 0.126 1.624 1.481
50SOL HW1 2 0.187 1.661 1.543
50SOL HW2 3 0.133 1.536 1.517
1.82461 1.82461 1.82461.pdb 格式大家更熟悉，但用它做大体系轨迹文件时，体积会非常大。所以我建议：能用 .gro 就用 .gro，最终分析阶段才可能用到 .pdb 格式的输出。
3.3.2 拓扑文件：.top 与 .itp
.top 是体系的全局拓扑，定义了原子类型、键合参数、非键参数、以及包含哪些分子。
; Include forcefield parameters
#include "amber99sb-ildn.ff/forcefield.itp"
[ molecules ]
SOL 12000
NA 10
CL 12.itp 是单个片段（如蛋白质、配体、溶剂）的拓扑文件。有一个细节我要提醒你：.itp 里的原子编号要从1开始，并在 .top 中通过 [ molecules ] 指定每个片段的拷贝数。
我曾经 在配体拓扑里忘了加 [ position_restraints ] 段，能量最小化时配体直接飞出了盒子。后来我养成习惯：每次修改拓扑后，先跑个极短的最小化试试。
3.3.3 参数文件：.mdp
所有模拟控制参数（温度、压力、步长、输出频率等）都放在 .mdp 文件里。GROMACS不支持在命令行设置这些参数，全得靠 .mdp 传递。
title = Protein in water
define = -DPOSRES
; Run parameters
integrator = md
nsteps = 500000
dt = 0.002
; Output control
nstxout = 1000
nstvout = 1000
nstenergy = 1000
nstlog = 1000
; Temperature coupling
tcoupl = v-rescale
tc-grps = Protein Non-Protein
tau-t = 0.1 0.1
ref-t = 300 300这里有个坑：不同GROMACS版本支持的 .mdp 参数不完全一样。如果你用5.x版本的 .mdp 去跑2023版的模拟，有些参数会直接报错，有些会被静默忽略——后一种更危险。
3.3.4 轨迹文件：.trr 与 .xtc
.trr 是完整精度轨迹，包含坐标、速度、力。体积大，但能用于精确重启动。
.xtc 是压缩轨迹，只保留坐标，丢失速度信息。体积小，适合长时间采样。
我个人习惯是：平衡阶段用 .trr 记录全部信息，生产阶段只存 .xtc 以节省磁盘。
3.3.5 能量文件：.edr 与 .log
.edr 是二进制能量文件，包含势能、动能、温度、压力、密度等时间序列。用 gmx energy 提取。
.log 是文本日志，记录每步的收敛信息和警告。跑完模拟先看一眼日志末尾，有没有提示“LINCS warning”或“DOMAIN DECOMPOSITION error”。
一个小技巧：用 gmx check -f md.xtc -f2 md.edr 可以快速核对轨迹帧数与能量步数是否一致。如果不一致，分析时肯定会出问题。
3.4 流程与文件类型总览
下面这张SVG图，把整个过程串了起来。你盯着它看两分钟，应该就能在脑子里建立清晰的文件流转图景。
3.5 文件命名的个人经验
文件名的混乱程度，直接决定你回溯问题时需要花多少时间。我个人建议一套命名规范：
- 使用三位数字前缀：
001_protein.gro、002_em.gro、003_nvt.gro…… - 不同模拟条件用后缀区分：
md_300K.xtc、md_310K.xtc、md_325K.xtc - 永远保留原始的
topol.top版本号，比如topol_v2.top
用工具链跑自动作业时，gmx会默认生成很多同名文件（比如 topol.top），手工验证时很容易搞混。这一点你千万要当心。
再次提醒：不要在生产模拟结束后，又回过头修改拓扑并重新运行，却不更新分析脚本。我见过有人花了三天时间分析，最后发现跑的是一套旧轨迹——顶级浪费。