企业官网建设流程全解析

做蛋白-配体分子动力学模拟这些年，Gromacs 一直是我用得最顺手的工具，没有之一。很多刚上手的人拿着一个蛋白加配体的复合物结构，直接丢进pdb2gmx，结果报错一片，误以为是软件装坏了。这套流程看着简单，真正跑通的人都知道，坑基本集中在配体拓扑生成、系统构建顺序和结果分析这三块。这篇文章我把从结构准备、力场选择、建盒子、平衡模拟到 RMSD/RMSF/氢键/结合能分析的完整流程完整走一遍，命令、参数都是可以直接抄的。适合刚接触 Gromacs 的初学者，也适合想系统梳理一遍蛋白-配体模拟流程的同行参考。

1. 项目整体设计与方案选型思路

1.1 全原子显式溶剂方案怎么定

拿到一个蛋白-配体复合物体系，第一步不是急着敲命令，而是想清楚用哪种模拟方案。蛋白-配体相互作用的核心是氢键、疏水接触、水桥和静电相互作用，这些细节只有在全原子显式溶剂模型里才能保真。隐式溶剂虽然计算快，但没法描述水分子介导的相互作用，配体结合袋里的水分子被排挤、置换这些关键事件直接看不见。粗粒化更适合膜蛋白大尺度聚集这类问题，研究结合模式早期筛选也不太合适。所以常规项目直接选全原子显式溶剂，代价是计算量大，但结果可靠。

盒子形状和大小也是有讲究的。正交盒子简单，但同样的截断距离下，同体积的立方形盒子要比十二面体盒子多装 20% 左右的水，白白浪费算力。常规做法是选择十二面体盒子，蛋白边缘到盒子边界留出至少 1.0 纳米的距离，保守一点我给 1.2 到 1.5 纳米。边界留太薄，配体在周期边界附近有可能撞上自己的镜像，PME 静电算出来的结果也容易失真，这种错误会直接体现在结合能分析上。

温度和压力条件，生化体系一般定为 300 K 和 1 bar，模拟人体内环境可以用 310 K。这里要留一个心眼：温度耦合算法、压强耦合算法在不同阶段用的策略不一样，后面第 4 章我会把参数细节展开。

1.2 蛋白和配体力场怎么搭配

力场是整个模拟的底层参数来源，这个选择直接决定了模拟结果的物理意义。蛋白链用 AMBER99SB-ILDN，配体用 GAFF2，配体电荷用 AM1-BCC，这是我现在最常用的组合。原因很简单，AMBER 系的蛋白参数和 GAFF2 的配体参数都遵循相似的原子类型定义和电荷处理习惯，接口顺滑，配合 Gromacs 自带的 pdb2gmx 和 acpype 这套工具链，不容易出参数打架的问题。

近两年 CHARMM36m 搭配 CGenFF 也很多人用，优势是蛋白和配体参数出自同一套参数体系，处理糖基化、脂修饰这类特殊残基时更友好，CGenFF 服务器可以直接给配体生成参数。缺点是需要联网提交，而且 CGenFF 对复杂杂环药物的电荷分配偶尔会给你警告。下表是我常用的三个方案对比：

我的建议是，没有特殊需求就直接照着 A 方案走，后面整篇教程的命令也统一用 A 方案演示。

1.3 全套模拟流程怎么安排

整个流程可以拆成六个环节：结构准备、配体参数生成、系统构建、能量最小化与平衡、成品模拟、结果分析。我习惯在项目目录下建子目录，每个环节一层，这样中途换力场、换条件时不会把中间文件搞乱。

project/ 01_struct/ # 初始结构和配体 02_prep/ # 拓扑文件和 itp 参数 03_build/ # 建盒子、溶剂化、加离子 04_eq/ # 能量最小化、NVT、NPT 05_md/ # 成品模拟 06_analysis/ # 轨迹处理与分析

这套结构看着土，但效果特别好。Gromacs 的中间文件命名混乱是出了名的，有了固定目录，三个月后回头翻数据还能知道当时到底跑了什么条件。

2. 前期准备：抓结构、配体参数与拓扑生成

2.1 环境与工具链

Gromacs 的安装方式我不展开太多，直接说最省心的路线。有 conda 就用 conda-forge 的包，一条命令装完，CUDA 版本也能自动匹配。需要编译优化的同学再考虑手动编译，手动编译的要点是 CPU 指令集和 GPU 算力得对上，不然跑起来性能差距很大。

conda create -n gmx python=3.10 conda activate gmx conda install -c conda-forge gromacs openbabel ambertools

这里注意，配体参数生成依赖 AmberTools 里的antechamber和parmchk2，很多人漏装。装完检查版本：

gmx --version antechamber --version

实际工作中我建议再装一个 PyMOL 或 VMD，用来可视化检查结构和轨迹。命令行能跑是一回事，结构有没有问题，眼睛看一遍比任何脚本都靠谱。

2.2 从 RCSB 获取复合物结构并做检查

从 RCSB 下载 PDB 结构时，我一般会先看一眼生物学组装（biological assembly）还是不对称单元（asymmetric unit）。做复合物模拟，最好用包含完整配体结合界面的那条链，不要拿一个晶格堆积产生的截断状态去模拟，否则配体可能直接被切掉一半接触面。

拿到 PDB 后在 PyMOL 里做三件事：删掉水分子、检查配体残基、确认有没有缺失残基。命令很简单：

remove solvent select LIG, resn LIG

这里要盯一下配体的残基名，有的 PDB 里小分子配体叫LIG，有的是INH、STI之类的专用代码。后面所有拓扑文件、索引文件里的组名都要和这里的名字保持一致。

还有一点容易被忽略，PDB 里如果存在多个构象的残基（ALT LOC），或者有磷酸化等修饰残基，pdb2gmx 很容易报错或者生成错误的拓扑。稳妥做法是用 PyMOL 清理成只保留一个构象，必要时把修饰残基单独处理。

2.3 配体三维结构与 GAFF2 参数生成

配体参数是整个流程里最容易翻车的地方。分三步走：三维结构、电荷计算、拓扑生成。

第一步，拿到配体的 SMILES 式，用 OpenBabel 生成三维结构。建议在原始报告或 PubChem 上找规范的 SMILES，不要在 PDB 里直接拿配体的初始坐标去跑参数，因为晶体结构里的氢原子和质子化状态不一定对。

obabel "CC(=O)Oc1ccccc1C(=O)O" -O ligand.sdf --gen3d obabel ligand.sdf -O ligand.mol2 -p 7.4

-p 7.4这一步是按生理 pH 分配氢原子和质子化状态。对于含羧基、氨基的配体，这一步很关键，电荷算错后面整个体系都会出问题。严格一些可以用 pKa 工具重新评估，但大部分常规配体 OpenBabel 的分配够用。

第二步，用antechamber计算 AM1-BCC 电荷：

antechamber -i ligand.mol2 -fi mol2 -o ligand_charge.mol2 -fo mol2 -c bcc -nc 0 -at gaff2 parmchk2 -i ligand_charge.mol2 -f mol2 -o ligand.frcmod

-nc 0是配体净电荷，根据实际质子化状态填。算完看输出文件里的Total charges，不是整数说明氢分配可能有问题，需要回到上一步检查。

第三步，用acpype生成 Gromacs 拓扑：

acpype -i ligand_charge.mol2 -o gmx -a gaff2

生成的文件在ligand_charge.amb2gmx/目录下，里面有LIG.itp、LIG.gro和POSRE_LIG.itp。把LIG.itp拷到工作目录时，记得先打开看一眼，确认里面的原子名称和顺序与LIG.gro完全一致，这一步不一致会导致后面合并系统时坐标错位。

3. 系统构建实操：从 pdb2gmx 到中性体系

3.1 pdb2gmx 生成蛋白拓扑

pdb2gmx 本身不认识配体，它会默认把非标准残基当作错误处理。所以第一步是把配体从复合物中分离出来，只留蛋白去生成拓扑。

gmx pdb2gmx -f protein_only.pdb -o protein_processed.gro -ignh -water spce -ff amber99sb-ildn

-ignh表示忽略输入结构里的所有氢原子，由 pdb2gmx 按力场规则重新加氢。这一步能规避很多氢原子命名不匹配的问题。如果蛋白有多条链，且链间有稳定的相互作用，可以加-merge no保持链分离；如果想让链间形成链间二硫键，也可以手动确认-ss参数。

生成后目录里会出现topol.top、protein_processed.gro和若干posre*.itp文件。先别急着往下走，打开topol.top看一眼：

[ molecules ] ; Compound #mols Protein_chain_A 1

确认蛋白链数和残基数正常。这时候还没有配体，因为 pdb2gmx 完全不认识它。

下一步是把配体坐标和拓扑合并到蛋白系统里。我有一个比较笨但可靠的做法：手动拼接。先把protein_processed.gro复制一份，然后在它末尾粘贴LIG.gro除文件头之外的坐标行，再把总原子数改成两者之和。例如蛋白是 3800 个原子，配体是 28 个原子，合并后第一行写上 3828。这一步也可以用脚本实现，但手动做一次反而能帮你理解 Gromacs 坐标系的结构。

然后在topol.top里做两处修改：一是在蛋白 itp 引入后面加一行#include "LIG.itp"，二是在[ molecules ]段加一行LIG 1。顺序一定要和 gro 里的原子顺序对应：蛋白坐标在前，配体坐标在后，那就先写 Protein_chain_A，再写 LIG。顺序错了，grompp 会直接报错，或者在更坏的情况下静默错位，后面结果全废。

3.2 建盒子和溶剂化

有了完整的 complex 拓扑和坐标，可以建盒子了：

gmx editconf -f complex.gro -o complex_box.gro -bt dodecahedron -d 1.2

-d 1.2的意思是蛋白到盒子边缘的最小距离是 1.2 纳米。用白话说，就是给体系周围留出一圈 1.2 纳米厚的水层。盒子建太大浪费水分子数量，太小会出现周期影像干扰，这个参数我调过很多次，1.2 是性价比比较高的值。

溶剂化用solvate命令，它会自动往盒子里填充水分子：

gmx solvate -cp complex_box.gro -cs spc216.gro -o complex_solv.gro -p topol.top

-cs spc216.gro是预先平好的 216 个 SPC/E 水分子构型文件，Gromacs 直接内置。如果你在 pdb2gmx 里选了 TIP3P，这里要先跑一个gmx solvate -cs tip3p.gro之类的命令生成对应构型，我为了省事，从 pdb2gmx 到 solvate 都用 spce，力场也选 amber99sb-ildn，两边的水模型参数是配套的。

solvate 运行完，topol.top会自动追加SOL的行数。这时候再打开看一眼，确认 SOL 的分子数不是 0。如果 SOL 是 0，最常见原因是 gro 文件没有和 topol 同步，回头核对[ molecules ]顺序。

3.3 加离子中和体系

蛋白质和配体在水溶液里带净电荷，直接模拟会因为同号电荷排斥导致体系爆炸，加离子的目的就是中和净电荷，同时维持生理离子强度。

先生成 tpr 文件，这里用的 EM.mdp 是最简化的参数，只要让 grompp 能处理坐标和拓扑就行。

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

执行genion时，程序会提示你选择哪一组要替换成离子，输入SOL组的编号即可。-pname NA -nname CL指定正负离子种类为钠和氯离子；-neutral要求体系总的净电荷为零；-conc 0.15顺便加上 150 mM 的盐浓度，模拟生理环境。

需要注意，如果 genion 报错说电荷无法中和，多半是配体电荷算错了，回到 2.3 节检查配体质子化状态。另一个坑是加完离子后 topol.top 应该在[ molecules ]的末尾多出 NA 和 CL 行，如果没看到，说明 genion 用错了 tpr。

4. 三步模拟：能量最小化、NVT 与 NPT 平衡

4.1 能量最小化：为什么必须做，怎么算算完

很多人刚建好体系就急着重心跑长模拟，结果第一步就出现原子重叠导致的数值爆炸。能量最小化的本质是找局部势能极小值，把初始结构里不合理的近距离原子对推开。这一步不追求采样，只求体系达到一个稳定的起点。

EM.mdp 我一般写成这样：

integrator = steep emtol = 1000 emstep = 0.01 nsteps = 50000 cutoff-scheme = Verlet coulombtype = PME rcoulomb = 1.0 rvdw = 1.0

最陡下降法（steep）对付初始重叠最稳，emtol设 1000 的意思是当最大受力降到 1000 kJ/mol/nm 以下就认为收敛。正常情况几百步就收住，如果 5 万步还没结束，说明初始结构冲突太严重，先回 PyMOL 里看是不是配体和残基嵌进去了。

跑和检查：

gmx grompp -f EM.mdp -c complex_neutral.gro -p topol.top -o em.tpr gmx mdrun -deffnm em -v gmx energy -f em.edr -o potential.xvg

用gmx energy输入Potential这个词，画出的势能曲线应该在末尾趋于平缓。如果曲线还在持续大幅下降，就再加大 nsteps 重新跑一次 EM，直到势能不再明显变化再进平衡。

4.2 NVT 加温：温度耦合怎么选

能量最小化完成后，体系处于 0 K 的势能极小，需要逐渐加热到目标温度。NVT 阶段用位置限制把蛋白重原子固定住，只让水、离子、配体和蛋白的氢原子运动，这样体系能在不偏离初始结构的前提下慢慢升温。

NVT.mdp 核心参数：

integrator = md dt = 0.002 nsteps = 50000 define = -DPOSRES tcoupl = v-rescale tc-grps = Protein_LIG Water_and_ions tau_t = 0.1 0.1 ref_t = 300 300 pcoupl = no constraints = h-bonds constraint-algorithm = LINCS

tc-grps把蛋白和配体作为一组，水和离子作为另一组，各自独立耦合。温度耦合算法我强烈建议用 v-rescale 而不是 Berendsen，v-rescale 能在保持合理温度涨落的同时避免热浴过度干扰体系动力学。

运行命令：

gmx grompp -f NVT.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr gmx mdrun -deffnm nvt -v gmx energy -f nvt.edr -o temperature.xvg

输入Temperature看温度曲线，正常情况下应该在 300 K 附近波动，偏差在正负 5 K 以内。如果温度持续偏高且不收敛，检查-r是否给了正确的参考结构，位置限制文件的原子选择有没有覆盖全部重原子。

4.3 NPT 平衡：压强耦合与体系密度

NVT 平衡后体系温度已经正确，但盒子的体积和压强还没有达到目标。NPT 阶段在继续限制蛋白重原子的同时，打开压浴，让盒子体积自动调整到与 1 bar 压强匹配。

NPT.mdp 相比 NVT 需要加这几项：

pcoupl = C-rescale pcoupltype = isotropic tau_p = 2.0 ref_p = 1.0 compressibility = 4.5e-5

Gromacs 2021 之后我常用 C-rescale 代替 Berendsen，C-rescale 对压强控制更稳定，不容易出现 Berendsen 那种在极端体系里压强涨落被低估的问题。compressibility设为水的等温压缩系数，这个值直接照抄 4.5e-5 1/bar 即可。

运行：

gmx grompp -f NPT.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr gmx mdrun -deffnm npt -v gmx energy -f npt.edr -o pressure.xvg -o density.xvg

分别看 Pressure 和 Density。水的密度应该接近 1000 kg/m^3，如果密度结果差出几十个单位，检查盒子大小和配体拓扑，很可能是配体体积参数明显不对导致溶剂壳层异常。压强瞬时波动有几百 bar 是正常现象，看整体平均在 0 附近就行。

平衡阶段总时长我一般给 100 ps，也就是 NVT 和 NPT 各 5 万步，对大多体系足够。如果体系的配体很大或者溶剂化层很特别，可以延长到 500 ps，这个灵活调整就好。

5. 成品 MD 运行与轨迹预处理

5.1 成品模拟参数与算力评估

NPT 平衡完成后，去掉位置限制，开始正式的成品模拟。MD.mdp 的关键点是压强耦合改用 Parrinello-Rahman，温度耦合继续保持 v-rescale，积分器和步长不变。

integrator = md dt = 0.002 nsteps = 50000000 tcoupl = v-rescale tc-grps = Protein_LIG Water_and_ions tau_t = 0.1 0.1 ref_t = 300 300 pcoupl = Parrinello-Rahman pcoupltype = isotropic tau_p = 2.0 ref_p = 1.0 compressibility = 4.5e-5 constraints = h-bonds constraint-algorithm = LINCS nstxout-compressed = 5000

nsteps = 50000000和dt = 0.002 ps算下来是 100 ns，对蛋白-配体体系来说，100 ns 已经可以观察到配体的稳定结合模式，如果要算结合自由能，建议跑到 200 ns 以上取平衡段。nstxout-compressed = 5000表示每 10 ps 存一帧轨迹，这样 100 ns 的 xtc 文件不会大到没法处理。

运行前先做一个小测试，对 nsteps 按 1000 步试跑一次，用gmx mdrun -deffnm md_test加-v输出速度，估算一天能跑多少 ns，再决定正式跑多长。

正式运行：

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

并行参数的设置，如果机器是单机多核，直接gmx mdrun -deffnm md -nt 16就行，不需要纠结 MPI 和 OpenMP 的层数。有 GPU 就把-nb gpu或-bonded gpu打开，速度提升明显。

5.2 轨迹预处理：一定要做的 pbc 处理

Gromacs 模拟用的是周期性边界条件，轨迹里的原子坐标在跨越盒子边界时会发生“跳跃”。如果直接拿原始轨迹算 RMSD 或距离，你会得到一条完全错误的暴涨曲线。

正确顺序是先让分子通过周期性边界恢复完整，再做质心对中。实际命令是两步合并成一条：

gmx trjconv -s md.tpr -f md.xtc -o md_center.xtc -pbc whole -center

第一次会让你选对中使用的组，选择Protein_LIG（蛋白和配体整体），确保整个结合复合物在盒子中心；第二次让你选输出组，选System，这样水分子也保留在轨迹里。如果之后只分析蛋白-配体相互作用，也可以选Protein_LIG输出，文件更小。

我见过不少人跳过这步直接分析，结果 RMSD 上窜到几十纳米，其实不是蛋白真的漂了，只是原子在边界复制箱之间跳来跳去。这个预处理一定要养成肌肉记忆。

5.3 结合能快速评估：MM/PBSA 可以做但别贪快

如果想对配体和蛋白的结合强度做初步评估，可以用 MM/PBSA 方法。Gromacs 本身不带这个工具，我用的是gmx_MMPBSA，它底层调用 AMBER 的 MMPBSA.py，但能直接吃 Gromacs 的 tpr 和 xtc 文件。

conda install -c conda-forge gmx_MMPBSA

准备一个输入文件：

&general sys_name="complex", startframe=1, endframe=100, interval=10, / &gb igb=2, saltcon=0.150, /

运行：

gmx_MMPBSA -O -i mmpbsa.in -cs md.tpr -ct md_center.xtc -cg Protein_LIG LIG

注意，MM/PBSA 的结果非常依赖轨迹段选取和参数模型，不同文章里的绝对值没有可比性，它更适合用来在同一条轨迹内比较不同突变体，或者评估同一配体不同结合姿态的差异。这个工具玩熟了可以出一篇单独的文章，这里点到为止。

6. 结果分析：把轨迹变成结论

6.1 RMSD：先看体系稳不稳

RMSD 是整个分析里最基础也是最重要的指标，它回答一个核心问题：模拟过程中蛋白骨架有没有跑偏，配体有没有稳定在结合口袋里。

先建立索引文件：

gmx make_ndx -f md.tpr -o index.ndx

在交互界面里输入keep 1并回车后，可以用a Protein之类命令查看分组编号。实际使用时 Gromacs 输出的编号可能不同，大家按自己体系的实际编号操作。

蛋白骨架 RMSD：

gmx rms -s md.tpr -f md_center.xtc -o rmsd_backbone.xvg -tu ns -n index.ndx

选择用蛋白骨架做最小二乘拟合的参考组，再选择要计算的组。一般选 Backbone 或 C-alpha。看输出的曲线，蛋白质 RMSD 在 1 到 3 埃的范围内波动的体系算稳定；如果 RMSD 一直攀升，说明蛋白还在缓慢漂移，要么模拟时长不够，要么初始结构本身不够稳定。

配体 RMSD 要小心一个陷阱：如果直接用原始轨迹去算，配体在结合口袋里转动会导致 RMSD 虚高。正确的做法是先用蛋白骨架对轨迹做最小二乘对齐，再算配体 RMSD，这样考察的是配体相对于结合口袋的相对位置变化。

gmx rms -s md.tpr -f md_center.xtc -o rmsd_lig.xvg -tu ns -n index.ndx

参考组选 Backbone，计算组选 LIG。配体 RMSD 在零点几到 1 纳米内波动，说明结合模式稳定；如果配体 RMSD 持续上升并伴有跳跃，可能在模拟过程中发生了解结合或姿态翻转，这时候要回到轨迹动画确认。

6.2 RMSF：看哪些残基在动

RMSF 看的是每个残基在整个轨迹中相对平均位置的波动幅度，用来找蛋白里柔性大的区域，比如 loop 区和配体结合位点附近的动态变化。

gmx rmsf -s md.tpr -f md_center.xtc -o rmsf_CA.xvg -n index.ndx -p

-p参数输出每个残基的平均 RMSF。建议选择 C-alpha 原子组来计算，结果画出来之后，高值区域通常对应表面 loop；结合位点附近的残基 RMSF 低，说明口袋在模拟过程中是稳定的，这是个很关键的质量指标。

RMSF 分析需要结合结构来解释，单纯看峰值没有意义。比如某个 loop 是已知的底物入口，RMSF 高说明这个区域在配体存在时仍有明显柔性，这一信息在描述蛋白构象变化和结合机制时经常用得上。

6.3 回转半径、氢键与溶剂可及表面积

回转半径（Rg）反映蛋白整体压实程度。蛋白越紧凑，Rg 越小；蛋白膨胀或部分解折叠，Rg 会上升。

gmx gyrate -s md.tpr -f md_center.xtc -o gyrate.xvg

正常的蛋白-配体复合物在平衡段 Rg 应该在平均值附近小幅波动。如果 Rg 有单调上升趋势，要警惕蛋白在逐渐解折叠，或者配体结合导致了不稳定的构象变化。

蛋白与配体之间的氢键是判断结合稳定性的直接证据。计算前先做一个包含 LIG 和 Protein 的索引，然后运行：

gmx hbond -s md.tpr -f md_center.xtc -n index.ndx -num hbond_prot_lig.xvg

选择参与氢键计算的两组，一般选 Protein 和 LIG。输出文件会显示氢键数量随时间的变化。如果整个模拟过程中持续存在一到三个氢键，说明配体与蛋白在结合口袋里有稳定的极性相互作用。氢键数的波动大，说明结合是动态的，可能在多个亚态之间切换。

溶剂可及表面积（SASA）用来评估蛋白整体或结合袋区域的溶剂暴露情况。结合口袋的 SASA 如果明显低于同尺寸表面残基的平均值，说明口袋被配体占据得比较踏实。

gmx sasa -s md.tpr -f md_center.xtc -o sasa.xvg -surface Protein -output Protein

6.4 配体-蛋白关键相互作用距离追踪

RMSD 和氢键是全局指标，很多细节需要精确到原子对。比如配体的某个羧基氧和蛋白的某个精氨酸胍基之间形成了盐桥，这个距离随时间的变化直接决定这个相互作用能不能稳定存在。

用gmx distance可以追踪任意两个原子或两团原子质心之间的距离：

gmx distance -s md.tpr -f md_center.xtc -select 'com of group LIG plus com of group Protein' -oall dist.xvg

select语法里，com of group会计算一个组的质心，plus表示计算两质心之间的距离。更细致的做法是直接在索引文件里定义指定的原子对，比如配体某原子的编号和蛋白某残基侧链某原子的编号，这样距离曲线非常明确。

处理结果时设定一个判据：氢键的距离阈值一般是 0.35 纳米，盐桥是 0.4 纳米左右，疏水接触看的是 0.5 纳米内的原子对数量。距离曲线稳定在阈值以下，说明该相互作用在整个模拟中基本维持，可以在论文里作为关键相互作用证据。

7. 常见问题与排查技巧实录

7.1 拓扑与原子报错

这一块是刚上手时的重灾区。我把这几年遇到最多的报错整理成一个速查表：

配体拓扑合并完成后，我强烈建议用一条软验证命令查看系统是否正常：

gmx grompp -f EM.mdp -c complex_neutral.gro -p topol.top -o test.tpr

如果 grompp 能通过，拓扑基本没问题。很多人在合并阶段出错，其实是在 gro 的总原子数上少加或多加了一行，这种错误往往要到跑 NVT 才会爆炸，排查起来更费劲。

7.2 模拟发散与温度压力异常

模拟跑到一半突然报 NaN 是让人最头疼的问题。第一步看输出里最后几行的能量信息，如果是 LJ 能量突然变成天文数字，说明有原子在模拟中撞到了一起。常见原因有：时间步长太大、温度耦合组定义有问题、某个约束没加成功。

优先排查方案：把 dt 从 0.002 降到 0.001，同时确认constraints = h-bonds已经开启；检查 NVT 阶段有没有真正加上-DPOSRES位置限制；再确认温度耦合组tc-grps的命名与实际索引文件里的组名一致。

温度异常过热且下不来，大概率是位置限制没作用。NPT 阶段出现了压强严重漂移，检查compressibility是否写成了水的值，这一点经常有人写错导致盒子体积崩溃。

压力曲线瞬时波动几百 bar 是正常的，不需要干预。真正需要关注的是压力平均值持续偏正或偏负超过几十 bar，这时候检查ref_p和tau_p，顺便看一眼盒子密度是否正常。

7.3 结果分析里的那些“坑”

跑完模拟只是完成了一半工作，分析阶段也有很多陷阱。第一个坑是 pbc 处理绕过，直接算距离和 RMSD，出来的曲线自己都解释不通。第二个坑是用配体 RMSD 时没有用蛋白骨架做对齐，结果把结合口袋内的旋转误判断为去结合。

第三个坑是氢键分析里组选择搞错。gmx hbond 要求两组原子互为供体和受体，如果选择的组里包含了太多无关残基，氢键数会虚高，建议在 make_ndx 里单独定义配体组和结合口袋残基组，再跑分析。

第四个坑是拿 MM/PBSA 的绝对值不同体系横向比较，这个我前面强调过。不同轨迹长度、不同帧间隔、不同 igb 模型算出来的结合能有很大差异，只能作相对比较。

踩过这些坑之后，我现在接手一个新体系的固定动作是：先用最小化流程把结构、拓扑跑通，再跑一段 1 纳米的测试模拟确认温度和压力稳定，然后才正式排产。这套流程看着多跑几步，实际省下的排查时间远比多花的算力值。配体拓扑合并、pbc 预处理、RMSD 对齐这三点，只要在第一次跑就按照正确流程来做，后面基本不会遇到让人熬夜的疑难杂症。