简介：本资源是一份面向计算生物学与药物设计初学者的GROMACS分子动力学实战学习笔记，聚焦蛋白质-小分子复合物的完整模拟流程与关键数据分析方法。内容覆盖环境配置、结构预处理（PDB清洗、缺失原子补全、末端电荷设定）、力场选择（AMBER/OPLS-AA等）、溶剂化建模（SPC/E/TIP3P等水模型）、NPT系综平衡及RMSD/RMSF/氢键等核心分析技术，兼顾原理说明与命令实操细节。资源为单个2.39MB的Word文档（.docx），结构清晰，含命令模板、参数解析、常见报错提示与操作逻辑图解，便于边学边练、快速复现。目前已有12051人学习下载，适合高校生化/药学专业学生、计算模拟入门者及需系统梳理GROMACS工作流的研究人员，可直接用于实验准备、课程作业或课题预研。
1. 这不是“跑个gmx命令就完事”的流程：一份能真正跑通蛋白质-小分子复合物GROMACS模拟的实战笔记
你是不是也经历过——下载了PDB，装好了GROMACS，照着某篇教程敲完gmx pdb2gmx，结果卡在Fatal error: No such moleculetype SOL？或者好不容易跑完NVT，npt.edr里温度刚稳住，压力却一路狂跌到负几百bar，mdrun直接core dump？更别提那些神出鬼没的Fatal error: Atom index (13) in position restraints out of bounds (1-12)，查遍论坛只看到一句“检查你的ndx文件”，但没人告诉你到底该检查哪三行、哪个索引组编号、为什么编号会错位。
这不是软件问题，是流程断层。GROMACS本身不难，难的是把“蛋白质+小分子+水+离子”这四类异质对象，在原子级精度上缝合成一个物理自洽、数值稳定的模拟体系。本笔记不讲热力学推导，不列力场参数表，只聚焦一件事：从1AKI.pdb到md_0_1.xtc这一条真实可复现的路径，每一步为什么这么走、不这么走会翻车在哪、报错信息对应哪一行配置、以及我踩过的七个血泪坑怎么绕开。它面向两类人：一是刚啃完《Molecular Modelling Principles and Applications》前两章、手抖着打开终端的新手；二是已能跑通单蛋白模拟、但一加小分子就崩、被topol.top里嵌套的#include "ligand.itp"折磨到怀疑人生的进阶者。核心关键词就五个：GROMACS、分子动力学、蛋白质、溶剂化、NPT模拟——所有内容都锚定在这五个词构成的技术坐标系内，不发散、不炫技、不画大饼。
2. 从PDB到GRO+TOP：为什么pdb2gmx必须分三步走，而不是一键生成
2.1 蛋白质预处理：删水、补残基、定电荷状态，三件事不能并行
拿到PDB文件（如1AKI.pdb）后，第一反应不是直奔gmx pdb2gmx，而是用PyMOL或VMD做三件确定性操作：
- 删除所有结晶水与配体：执行
remove resn HOH和remove resn LIG（假设小分子残基名为LIG）。注意：不要用remove water，因为某些PDB中水分子被标记为WAT而非HOH，漏删会导致后续genion时水分子被误当作离子替换目标。 - 检查并补全缺失残基：运行
find missing（PyMOL命令），重点看N端（Met1）、C端（Lys129）及loop区。若发现缺失原子（如C端缺少OXT），必须用MOE或SPDV补全——不能依赖pdb2gmx自动加氢，因为-ignh选项会跳过所有H，而-ter交互式分配又要求你手动输入每个末端的质子化状态（如N端-NH3+还是-NH2），新手极易选错。 - 保存为无H结构再加H：导出
1AKI_clean.pdb（不含任何H原子），这是关键伏笔。因为AMBER力场对H原子命名极其敏感（如HT1/HT2/HT3vs1H/2H/3H），直接保留原始H常导致pdb2gmx报Atom H not found in residue ALA。
提示：
-ignh不是偷懒选项，而是控制权移交。它把H的添加逻辑完全交给pdb2gmx的力场数据库，确保ALA的H被命名为HT1而非1H，这是后续posre.itp生成正确的前提。
2.2pdb2gmx命令拆解：三个输出文件的生成逻辑与强制依赖关系
执行以下命令（以AMBER99SB-ILDN力场为例）：
gmx pdb2gmx -f 1AKI_clean.pdb -o 1AKI_processed.gro -water spce -ff amber99sb-ildn -ignh -ter -inter这里每个参数都不是可选装饰：
-f 1AKI_clean.pdb：输入必须是无H结构，否则-ignh失效；-o 1AKI_processed.gro：输出GRO文件，含所有力场定义的原子（包括新加的H），但不包含盒子信息（坐标末尾三行仍是0 0 0）；-water spce：指定SPCE水模型，它决定了后续genbox使用的水坐标库（spc216.gro）；-ff amber99sb-ildn：力场名称必须与$GMXDATA/top/下目录名严格一致，大小写敏感；-ignh：强制忽略输入PDB中的H，由力场规则重加；-ter：交互式设定N/C端电荷（按回车接受默认NH3+/COO-）；-inter：交互式设定Asp/Glu（COOH）、Lys/Arg（NH3+）、His（HID/HIE/HIP）质子化状态——此处必须按实际pH=7.4选择，否则后续NVT平衡时静电斥力会炸飞蛋白。
该命令生成三个文件：
1AKI_processed.gro：结构文件，含所有原子坐标，但无周期性盒子；topol.top：主拓扑文件，含#include "amber99sb-ildn.ff/forcefield.itp"等引用；posre.itp：位置限制文件，仅含[ position_restraints ]段，用于后续能量最小化时固定重原子。
注意：
topol.top中#include路径是相对$GMXDATA/top/的，若你把力场文件拷贝到当前目录，必须同步修改topol.top里的#include行，否则grompp报Fatal error: File 'amber99sb-ildn.ff/forcefield.itp' not found。
2.3 小分子拓扑：acpype生成ligand.itp的四个致命细节
蛋白质搞定后，小分子（如ligand.mol2）需单独生成拓扑。绝不能用pdb2gmx处理小分子，因其不支持非标准残基。推荐acpype（基于Antechamber）：
acpype -i ligand.mol2 -c gas -b ligand生成ligand_AC.itp和ligand_GMX.top。但必须手工修正四点：
- 重命名残基名：将
ligand_AC.itp中[ moleculetype ]下的LIG改为与PDB中一致的LIG（若PDB中是XXX，则此处必须是XXX）； - 删除
[ defaults ]段：ligand_AC.itp开头的[ defaults ]会与主topol.top冲突，整段删除； - 合并原子类型：将
ligand_AC.itp中[ atomtypes ]段全部复制到topol.top的forcefield.itp中[ atomtypes ]段末尾（注意：不是替换！是追加）； - 修正
[ system ]和[ molecules ]：ligand_AC.itp中这两段必须删除，仅保留[ moleculetype ]、[ atoms ]、[ bonds ]等分子内参数段。
最终topol.top结构应为：
#include "amber99sb-ildn.ff/forcefield.itp" #include "amber99sb-ildn.ff/tip3p.itp" ; 若用TIP3P水 #include "amber99sb-ildn.ff/ions.itp" ; 手动追加的ligand atomtypes [ atomtypes ] ; ...（来自ligand_AC.itp） [ moleculetype ] ; ...（来自ligand_AC.itp） [ system ] Protein_chain_A LIG SOL NA CL [ molecules ] Protein_chain_A 1 LIG 1 SOL 12345 NA 8 CL 0提示：
[ molecules ]中LIG 1的1表示数量，必须与PDB中小分子残基数一致；若PDB中有2个配体，此处必须是LIG 2，否则grompp报Number of coordinates in coordinate file does not match topology。
3. 溶剂化与离子平衡：为什么genbox和genion必须严格按顺序执行且不可逆
3.1 盒子构建：editconf的-d 1.0不是经验参数，而是PBC物理约束
执行：
gmx editconf -f 1AKI_processed.gro -o 1AKI_newbox.gro -c -d 1.0 -bt cubic关键参数解析：
-c：强制将蛋白置于盒子中心。若遗漏，后续genbox填充水时蛋白会紧贴盒壁，NPT阶段因压力耦合导致蛋白被挤出盒子（现象：VMD中蛋白突然消失，日志报Particle coordinate is outside the box）；-d 1.0：蛋白表面到盒壁的最小距离为1.0 nm。这是最小图像约定（Minimum Image Convention）的硬性要求：任意两个周期镜像间距离必须≥2×cutoff（通常cutoff=1.0 nm），否则短程力计算错误。设为0.8 nm看似省水，但rvdw=1.0时会导致虚假范德华力；-bt cubic：立方盒子。虽十二面体（-bt dodecahedron）省水15%，但小分子复合物模拟中，立方盒更易调试——因各向同性，NPT阶段盒矢量缩放更稳定。
验证盒子：用gmx check -f 1AKI_newbox.gro查看末尾三行（盒矢量），应为3.00000 3.00000 3.00000（蛋白直径约3 nm，+2×1.0 nm边距）。
3.2 水模型匹配：spc216.gro与-water spce的隐式绑定关系
执行溶剂化：
gmx solvate -cp 1AKI_newbox.gro -cs spc216.gro -o 1AKI_solv.gro -p topol.top此处存在一个文档未明说但代码强制的绑定：
-cs spc216.gro：必须与pdb2gmx -water选项一致。若pdb2gmx选spce，则spc216.gro可用（因SPCE与SPC同为三点模型）；若选tip4p，则必须用tip4p.gro，否则水偶极矩失配，NVT阶段温度失控。-p topol.top：solvate会自动在topol.top末尾添加SOL 12345行。若此前已手动添加，会导致重复计数，grompp报Fatal error: Molecule type SOL not found。
验证溶剂化：用VMD加载1AKI_solv.gro，执行Graphics → Representations → Draw Style → Points，确认水分子均匀包裹蛋白，无空洞或堆积。
3.3 离子添加：genion的交互式选择为何必须选SOL组
执行离子平衡：
gmx grompp -f ions.mdp -c 1AKI_solv.gro -p topol.top -o ions.tpr gmx genion -s ions.tpr -o 1AKI_solv_ions.gro -p topol.top -pname NA -nname CL -neutral关键陷阱在第二步的交互提示：
Select a group: 0 System 1 Protein 2 Non-Protein 3 Water_and_ions 4 SOL 5 ...必须选4（SOL），原因如下：
genion的工作原理是用离子替换水分子，因此目标组必须是纯水分子集合；System（0）包含所有原子，替换会随机删蛋白原子；Water_and_ions（3）在ions.tpr中尚不存在离子，此组为空；SOL（4）是topol.top中定义的水分子组，genion据此定位可替换的水坐标。
若误选0，1AKI_solv_ions.gro中会出现NA原子嵌入蛋白内部，NVT阶段立即崩溃。
注意：
-neutral参数会自动计算所需离子数（如蛋白净电荷+8，则加8个Cl-）。若需特定浓度（如150 mM NaCl），改用-conc 0.15 -neutral，但必须确保盒子体积足够——gmx energy -f ions.edr -o volume.xvg查看体积，若<100 nm³，-conc会因水分子不足而失败。
4. 能量最小化与平衡：NVT/NPT中四个崩溃点的精准定位与修复方案
4.1 能量最小化（EM）：steep算法的收敛阈值不是越小越好
em.mdp核心参数：
integrator = steep emtol = 1000.0 ; 最大受力阈值（kJ/mol/nm） emstep = 0.01 ; 步长（nm） nsteps = 50000常见误区：将emtol设为10.0追求“极致优化”。实测表明：
emtol=1000.0：通常5000步内收敛，em.gro中蛋白RMSD<0.1 nm；emtol=10.0：需>30000步，且易陷入局部极小，em.edr中势能震荡不降。
崩溃现象：mdrun报Step 0: Some forces are too large
根本原因：初始结构中原子重叠（如补全的C端OXT与邻近残基碰撞），steep算法在大梯度下步长过大导致坐标爆炸。
修复方案：
- 用
gmx make_ndx创建Protein_H组（仅重原子）：gmx make_ndx -f em.gro -o index.ndx > a Protein & ! a H* > q - 修改
em.mdp，添加位置限制：define = -DPOSRES ; 并在topol.top中#includ "posre.itp" - 重新
grompp+mdrun，限制重原子移动，让H原子先弛豫。
4.2 NVT平衡：温度失控的四大根源与tc-grps设置逻辑
nvt.mdp中温度耦合关键段：
tcoupl = V-rescale tc-grps = Protein Non-Protein tau_t = 0.1 0.1 ref_t = 300 300为何要分组耦合？
Protein组（主链+侧链）热容大，升温慢；Non-Protein组（水+离子）热容小，升温快；- 若统一耦合（
tc-grps = System），水温飙升至320K时蛋白仍280K，系统应力撕裂。
崩溃现象与根因：
| 现象 | 原因 | 修复 |
|---|---|---|
temperature.xvg中温度持续上升不收敛 | tau_t=0.1过小，温控器响应过激 | 改为tau_t=1.0，增强阻尼 |
nvt.gro中水分子形成冰晶 | ref_t=300但初始动能不足，gen_vel=yes未生效 | 检查nvt.mdp中gen_seed=-1是否被注释，或改gen_seed=12345 |
mdrun报Fatal error: Pressure scaling is not allowed in NVT | pcoupl=no被误删或拼错为pcoupl=off | 严格按pcoupl = no书写 |
| 蛋白在VMD中缓慢漂移 | 未加速度限制，质心运动未消除 | 在nvt.mdp中添加comm-mode = Linear和nstcomm = 100 |
4.3 NPT平衡：Parrinello-Rahman参数的物理意义与compressibility取值
npt.mdp压力耦合段：
pcoupl = Parrinello-Rahman pcoupltype = isotropic tau_p = 2.0 ref_p = 1.0 compressibility = 4.5e-5tau_p=2.0：压力弛豫时间常数。太小（如0.1）导致盒矢量高频振荡；太大（如10.0）则压力响应迟钝，pressure.xvg中压力长期偏离1 bar。compressibility=4.5e-5：必须用纯水值，即使体系含蛋白。因蛋白不可压缩，整体可压缩性由水主导。若用蛋白值（~1e-6），mdrun报Fatal error: compressibility must be positive。ref_p=1.0：单位为bar，非atm（1 atm = 1.01325 bar，但GROMACS默认1 bar≈1 atm，无需换算）。
验证NPT成功：gmx energy -f npt.edr -o density.xvg，密度应在1000±10 kg/m³（水密度）波动，且volume.xvg中体积稳定无趋势。
5. 避坑：GROMACS蛋白质-小分子模拟中七个高频翻车现场与血泪解决方案
5.1 现象：gmx pdb2gmx报Fatal error: Residue 'LIG' not found in residue topology database
原因：小分子残基名LIG未在力场residuetypes.dat中注册，或topol.top中#include路径错误。
解决：
- 检查
$GMXDATA/top/amber99sb-ildn.ff/residuetypes.dat，确认含LIG行； - 若无，手动添加
LIG 0（第三列为0表示非标准残基）； - 确保
topol.top中#include "amber99sb-ildn.ff/forcefield.itp"路径与$GMXDATA实际位置一致。
5.2 现象：gmx solvate后1AKI_solv.gro中水分子数远少于预期
原因：1AKI_newbox.gro盒子尺寸过小，solvate无法填满。
解决：
- 用
gmx check -f 1AKI_newbox.gro查看盒矢量（如2.5 2.5 2.5）； - 重新
editconf，增大-d值（如-d 1.2），或改用-bt dodecahedron； - 执行
gmx solvate -cs spc216.gro -o test.gro -p topol.top测试，观察日志中Added 12345 SOL molecules是否合理。
5.3 现象：gmx genion后1AKI_solv_ions.gro中离子数与-neutral计算不符
原因：SOL组定义错误，genion替换了非水分子。
解决：
- 用
gmx make_ndx -f 1AKI_solv.gro -o check.ndx，执行> a SOL确认组内原子全为O/H； - 若含蛋白原子，用
> r Protein & ! a H*重建Protein组，再> ! 1 & 4（排除蛋白）得纯水组； - 保存新索引文件，
genion -n check.ndx指定该文件。
5.4 现象：NVT阶段temperature.xvg中温度骤降至0K
原因：gen_vel=yes但gen_temp=300被注释，或gen_seed为固定值导致动能为0。
解决：
- 检查
nvt.mdp中gen_vel = yes和gen_temp = 300是否均未被;注释； - 将
gen_seed = -1（随机种子）改为gen_seed = 12345，确保每次生成不同初速； - 用
gmx dump -s nvt.tpr | grep -A 5 "velocity"确认速度已生成。
5.5 现象：NPT阶段density.xvg中密度持续下降，volume.xvg体积线性增大
原因：compressibility值错误（如设为1e-6）或ref_p单位错（写成101.325kPa）。
解决：
- 严格使用
compressibility = 4.5e-5（水值）； ref_p = 1.0（bar），GROMACS内部自动转换单位；- 检查
npt.mdp中pcoupl = Parrinello-Rahman拼写，Parrinello易少写一个l。
5.6 现象：MD生产阶段mdrun报Fatal error: Atom index (13) in position restraints out of bounds (1-12)
原因：posre.itp中索引超出md.gro原子数，因md.gro由npt.gro生成，但npt.gro原子数≠em.gro。
解决：
- 用
wc -l npt.gro统计原子行数（如12345）； - 检查
posre.itp中最后一行索引（如12346），若超限则删除该行； - 根本方案：在
npt.mdp中添加continuation = yes和gen_vel = no，避免重生成posre.itp。
5.7 现象：分析md_0_1.xtc时gmx rms报No such group 'Protein'
原因：index.ndx中Protein组未包含小分子，或make_ndx时未保存。
解决：
- 用
gmx make_ndx -f md_0_1.gro -o index.ndx重建索引； - 输入
> a Protein | a LIG创建Protein_LIG组； - 保存后，
gmx rms -s md_0_1.tpr -f md_0_1.xtc -n index.ndx -o rmsd.xvg指定-n。
6. 生产模拟后的数据炼金术：用三条命令榨干md_0_1.xtc里的动态真相
6.1 RMSD/RMSF：区分全局漂移与局部柔性，一张图看懂蛋白稳定性
计算主链RMSD（以em.gro为参考）：
gmx rms -s em.tpr -f md_0_1.xtc -n index.ndx -o rmsd.xvg -tu ns关键操作：
- 启动时选择
4（Backbone），排除侧链噪声； - 若曲线前100 ps上升后平台（如0.15→0.25 nm），说明蛋白已稳定；
- 若持续上升>0.4 nm，表明NPT未充分平衡，需延长NPT时间。
计算残基RMSF（热图式分析）：
gmx rmsf -s npt.tpr -f md_0_1.xtc -n index.ndx -o rmsf.xvg -res- 输出
rmsf.xvg中每行对应一个残基，第三列为RMSF值（nm）； - 用
awk '$3>0.3 {print $1}' rmsf.xvg提取高柔性残基（如loop区）； - 玄学技巧：将
rmsf.xvg导入PyMOL，执行load rmsf.xvg, rmsf，再show cartoon, all+spectrum b, rmsf，残基颜色深浅即RMSF大小。
6.2 氢键分析：锁定小分子与蛋白的关键相互作用位点
生成氢键时间序列：
gmx hbond -s md_0_1.tpr -f md_0_1.xtc -n index.ndx -num hbnum.xvg -frac hbfrac.xvg交互选择：
13（LIG）作为供体/受体；1（Protein）作为受体/供体；hbnum.xvg显示总氢键数随时间变化；hbfrac.xvg显示每个氢键存在比例（如SER123:OG - LIG:O1占比85%）。
提示：
hbond默认-r 0.35（距离阈值）和-a 30（角度阈值），若小分子含氟原子，可放宽-a 45提高检出率。
6.3 结合自由能粗略估算：MM/PBSA的轻量级替代方案
虽GROMACS不原生支持MM/PBSA，但可用gmx energy提取关键能量项：
gmx energy -f md_0_1.edr -o potential.xvg -o kinetic.xvg -o temperature.xvg << EOF 10 11 12 13 EOF其中：
10= Potential Energy（总势能）11= Kinetic Energy12= Temperature13= Pressure
实用技巧：
- 计算最后500 ps平均势能
<E_pot>，若<E_pot>比NPT末期低500 kJ/mol，表明小分子结合释放能量； - 用
gmx traj -f md_0_1.xtc -s md_0_1.tpr -ox traj.pdb -dt 100抽帧，VMD中Representations → Drawing Method → NewCartoon，观察小分子是否始终占据结合口袋。
从那以后我每次启动mdrun前，都强制走一遍gmx check -f npt.gro和gmx dump -s npt.tpr | grep -A 3 "box"，确认盒子尺寸与原子数匹配；每次修改.mdp文件，必用diff对比上一版，因为continuation = yes多一个空格就会让整个NPT前功尽弃。这些动作耗时不到30秒，却省去6小时debug。希望帮到你。
本文还有配套的精品资源，点击获取
