简介一份面向分子动力学初学者的GROMACS学习笔记系统梳理蛋白质-小分子复合物从结构准备到动力学模拟的完整流程。笔记覆盖AMBER、OPLS-AA、GROMOS、CHARMM、MARTINI等常用力场以及SPC、SPC/E、TIP3P、TIP4P、TIP5P水模型的选择逻辑结合具体命令讲解了pdb2gmx生成拓扑、editconf定义盒子、solvate添加溶剂、后续NPT模拟等关键环节。同时包含PDB预处理、电荷状态分配、缺失原子补全、氢原子忽略等实操细节并介绍了RMSD、RMSF、氢键分析等常用数据分析方法帮助研究者从模拟轨迹中提取结构变化与相互作用信息。压缩包含单个docx文档共2.39MB内容以分阶段笔记形式呈现适合终端操作时对照查阅。已有12046人学习对于希望快速建立GROMACS整体框架、减少踩坑经历的研究者很有参考价值。1. 这套 GROMACS 学习笔记把蛋白质-小分子模拟的全流程串成了能落地的操作手册分子动力学模拟这行有个很现实的问题资料满天飞但你搜到的教程要么只讲某一步要么默认你已经会了结果是把 VMD 打开结构、力场文件怎么改、为什么 NVT 老崩这些最磨人的环节全跳过去了。这套关于 GROMACS 中蛋白质-小分子复合物动力学模拟的笔记正好补上这个缺口——它不是讲原理的大部头而是从下载结构、补全缺失原子开始一路把 pdb2gmx、建盒子、加溶剂加离子、能量最小化、NVT/NPT 平衡到正式模拟的每一行命令都记录下来的实战流程。如果你是刚接触 GROMACS 的从业者或研究生照着这份笔记的步骤能把你之前零零散散拼出来的模拟流程理顺如果你已经跑过几个体系笔记里关于 .mdp 参数怎么调、体系崩溃的四个原因、ligand.itp 和力场文件冲突这几块也能帮你省下不少排查时间。2. 预处理与拓扑生成从 PDB 原始结构到 gro/top/posre 三件套2.1 输入文件的四个角色gro、top、itp、ndx 谁管什么GROMACS 的模拟体系里文件类型常让新手晕头转向但拆开看其实就四类活。结构文件负责提供原子坐标常见的有 pdb、gro、g96它们记录的是每个原子叫什么、属于哪个残基、坐落在什么位置拓扑文件 top 定义的是整个系统的化学连接性包括原子类型、电荷、键、角度和二面角也就是分子力场描述的全部内容itp 是拓扑的包含文件最典型的就是 pdb2gmx 生成的 posre.itp位置限制文件以及小分子配体的参数文件ndx 索引文件则是在分析或耦合时指定原子分组的工具比如把蛋白质和小分子各归为一组。初次跑流程的人最容易踩的坑是搞混 gro 和 top 的关系——gro 只是坐标top 才是化学身份。你换了结构文件但没同步更新拓扑后面 grompp 必然报坐标数量不匹配。这套笔记在开头就点明了运行需要的文件组合初始构型gro/pdb、参数控制文件mdp、力场文件top/itp、索引文件ndx。我的建议是每个模拟项目单独建目录把结构文件、拓扑、mdp、力场子目录分开放这样排查问题的时候不用在一堆 1AKI_processed.gro、1AKI_newbox.gro、solv_ions.gro 里大海捞针。2.2 pdb2gmx 的参数细节水模型、-ignh 与 NMR 结构处理拿到蛋白质结构文件之后第一步是清理。从蛋白质数据库下载的原始 pdb 通常带着结晶水还可能缺原子。笔记里提到的做法我比较认同用分子可视化软件直接选水分子删除然后用缺失检测找出断裂的残基。检查缺失时找结构里面标 MISSING 的位置补原子可以用专业建模软件处理。补完之后记得查看 C 端必要时把缺失的主链原子加上这一布偷懒的话后面 pdb2gmx 会报原子缺失错误。清完结构就轮到核心命令 pdb2gmx它干三件事补氢、根据所选力场生成拓扑、输出 GROMACS 格式的结构文件。基本形式是gmx pdb2gmx -f 1AKI_clean.pdb -o 1AKI_processed.gro -water spce -p topol.top -r posre.itp参数含义-f 指定输入结构-o 指定输出结构-water 指定水模型-p 指定拓扑输出文件名-r 指定位置限制文件。如果不加 -water程序会在终端里逐条问你想用哪个模型interactive 模式对新手反而更友好能让你一步步理解每个选择在干什么。注意 pdb2gmx 本身不改盒子所以生成的 gro 文件最后三行盒向量可能显示为零这个等 editconf 阶段再处理。关于水模型的选择GROMACS 自带五种常用模型SPC、SPC/E、TIP3P、TIP4P、TIP5P。SPC 是最早的简单点电荷模型运算最快但介电性质一般SPC/E 就是它的修正版把极化效应近似进了有效电荷里体相水性质更接近实验TIP3P 是三点模型的另一个代表很多力场都拿它做参数化基准TIP4P 增加了第四个虚位点把电荷放在 M 点上计算量上升但取向性更好TIP5P 是五点模型精度最高。做蛋白质-小分子复合物模拟时我们实验室的常规做法是用 OPLS-AA 力场配 SPC/E 或 TIP3P用 AMBER 力场就配 TIP3P。三点水模型的盒子文件可以直接用 spc216.gro 替代GROMACS 官网给的教程也是这么干的。2.3 力场选型OPLS-AA、AMBER、CHARMM、MARTINI 的边界在哪笔记里列出了 GROMACS 支持的多种力场AMBER、OPLS-AA、GROMOS、CHARMM、MARTINI这个选型直接影响后面所有模拟结果值得单独说。OPLS-AA全原子优化势液相模拟对有机液体和蛋白质溶液的表面积累很丰富参数覆盖面广做蛋白质-小分子复合物时小分子参数相对好找AMBER 系列在生物大分子尤其是核酸和蛋白质方面有大量验证配套的 GAFF 小分子力场参数也很通用很多发表的工作都用的 AMBER 力场配合 TIP3P 水结果可比性最强CHARMM 36 是凝聚相模拟的另一个主流选择和 CHARMM-GUI 的在线工具配合度高但参数文件和 GROMACS 的匹配要额外处理GROMOS 是联合原子力场氢原子是隐式的计算量小不少适合长时间尺度的膜蛋白模拟MARTINI 则是粗粒化力场一个珠子代表四到五个重原子适合做脂质自组装、膜蛋白聚集这类大尺度现象代价是牺牲化学细节。对做蛋白质-小分子复合物的从业者我的选型标准很简单如果配体里有常见的药物分子骨架优先用 OPLS-AA 或 AMBER因为小分子参数能从自动化工具或已有文献里直接拿到。如果主要研究蛋白质本身构象变化选 CHARMM36 和 AMBER 更稳妥。如果上万纳秒的聚集行为才是重点才考虑 MARTINI。选力场还要注意一件事pdb2gmx 选择的力场会同时决定水模型和离子参数SPC/E 和 TIP3P 不能混用不同力场的参数化逻辑改来改去最容易让结果变得不伦不类。pdb2gmx 里有两个选项要特别留意。-ignh 表示忽略 pdb 里的氢原子这对 NMR 解析的结构尤为重要因为 NMR 结构里氢原子的命名约定很可能和 GROMACS 力场不一致。晶体结构通常不含氢不需要这个选项。-ter 是交互式分配 N 端和 C 端电荷状态-inter 则进一步分配 Glu、Asp、Lys、Arg、His 这些残基的质子化状态以及是否形成二硫键这一步对 pH 环境和氧化还原状态敏感的体系很关键。处理完生成 topol.top 和 posre.itp 之后建议先把 top 文件打开看一眼确认分子类型数量、原子总数和后面 solvate/genion 更新后一致。3. 体系构建盒子尺寸、溶剂化与离子中和的参数选择3.1 editconf 建盒子为什么 -d 1.0 和 -c 不能省蛋白质有了干净的拓扑接下来要把它放进一个「房间」里。这个房间就是模拟盒子需要在里面装满水而且要让蛋白质和它自己的周期镜像永远隔开足够距离否则边界那边的镜像分子会对当前的蛋白质施加虚假的力。这就是最小镜像约定的含义。实现这一步的命令是 editconfgmx editconf -f 1AKI_processed.gro -o 1AKI_newbox.gro -c -d 1.0 -bt cubic参数说明-c 把蛋白质放到盒子中心-d 1.0 表示蛋白质任意原子到盒子边缘的最小距离是 1.0 nm-bt cubic 指定立方体盒子。为什么要强调 -c如果不加蛋白质可能靠近盒子一角溶剂化时一侧水太多一侧太少最麻烦的是如果蛋白质穿出盒子边界后面运行会直接报错或产生完全无意义的结果。关于 -d 1.0 这个值它保证蛋白质任意两个周期镜像之间的距离至少 2.0 nm这个尺度对 GROMACS 常用的 1.0 nm 截断半径已经足够。如果你用更长的截断或 PME 的格点间距设得更大这个值还得往上加。盒子类型的选择直接影响计算量。立方盒最直观但空间利用率低菱形十二面体的体积大约是同周期距离立方盒的 71%意味着溶剂分子数量少计算量明显下降截断八面体也类似。笔记里作者在某个项目里用了 cubic但更建议用 dodecahedron特别是体系稍大或者打算跑上百纳秒的时候。体积缩小 29% 带来的加速非常可观。3.2 solvate 加水spc216.gro 为什么能通用于三点模型盒子定好之后往里灌水这一步由 solvate 完成gmx solvate -cp 1AKI_newbox.gro -cs spc216.gro -o 1AKI_solv.gro -p topol.top参数说明-cp 是溶质的坐标文件protein-cs 是溶剂坐标文件-o 是输出结构-p 指定要更新的拓扑文件。solvate 会自动把水分子数量和类型追加到 top 的 molecules 段里所以不用手动改拓扑。把 spc216.gro 作为三种三点水模型SPC、SPC/E、TIP3P的通用溶剂盒子是官方教程的常规做法因为三点水坐标形式相同力场参数在后续 grompp 时会被 .mdp 和它p文件中设定的模型替换掉。加完溶剂后强烈建议可视化检查一次这一步看起来多余但省掉它风险很大盒子太小导致部分水分子被挤到蛋白质内部或者蛋白质本身跑到了盒子边界。笔记给的方法是把 gro 转成 pdb 用分子可视化软件看或者直接用 VMD 打开 gro 文件。密度检查也有用如果体系密度明显偏离 0.9~1.1 g/cm³ 的范围盒子尺寸大概率出了问题。注意这个检查的价值后面正式模拟发现的异常结构百分之九十在前期已经埋下了隐患。3.3 genion 加离子-neutral 背后到底做了什么蛋白质表面带有净电荷直接跑模拟不带离子会导致体系整体带电不仅在物理上不真实在 Ewald 求和下还会造成严重伪影。GROMACS 的加离子工具是 genion它的工作原理是通过 .tpr 文件读取拓扑把指定组里的一些水分子替换成离子。命令如下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执行 gmx genion 的时候终端会列出系统里的所有分组——Protein、Non-Protein、Water、SOL 等必须输入盐溶液那组的编号通常是 SOL 的编号否则离子可能替换掉蛋白质原子后面的模拟过程会直接崩溃。这里还有个小问题grompp 要用一个 mdp 文件生成 tpr这个 mdp 选什么样其实不太重要重要的是能让系统完成原子描述。笔记里用了 steepest descent 的配置先做能量最小化再插离子也是可行的不过更常规的做法是直接在原来的体系上跑 grompp 生成 tpr不关心 tpr 对应的动力学是否合理。-neutral 这个参数的作用是计算系统当前净电荷然后自动添加足够数量的反离子——正离子和负离子成对加入并不会影响总量。比如系统净电荷 8genion 就会添加 8 个 Cl⁻。不需要电中性的话也可以手动指定 -nn 或 -np 来控制离子数量。关于离子浓度默认加的只是中和电荷的量生理条件下还需要额外加 0.15 mol/L 的 NaCl可以结合手动指定来实现。3.4 自定义力场文件ligand.itp 的[ atomtypes ]替换与模块裁剪做蛋白质-小分子复合物模拟最痛点永远是小分子的参数。GROMACS 没有现成的小分子力场参数需要从其他工具生成然后合并到力场文件里。笔记里记录了一个经典的合并方式把力场文件所在目录复制到当前工作目录然后修改 forcefield.itp把 ligand.itp 里的 [ atomtypes ] 部分替换掉原力场对应部分同时删掉 ligand.itp 里的 [ defaults ]、[ system ]、[ molecules ] 三个段落。为什么要删这三个段落因为 [ defaults ] 定义的非键相互作用规则和 [ atomtypes ] 只应出现一次[ system ] 和 [ molecules ] 是顶层拓扑的段落重复出现会造成多个 system 定义和分子数量重复grompp 和 mdrun 直接报语法错误。复杂的是ligand.itp 里的 [ atomtypes ] 如果和 forcefield.itp 里已有的类型雷同替换时会把所有引用该类型的原子团参数一起覆盖。最稳妥的做法其实是用工具生成小分子参数时选一个和主力场兼容的原子类型命名比如 OPLS-AA 体系的配体参数生成后再核对键、角、二面角类型是否被主力场支持。我在实操中遇到过不少次配体参数里出现主力场没有的 LJ 参数类型grompp 报错的方式还特别隐蔽——只在 mdrun 阶段以 LINCS 警告形式出现。4. 平衡模拟避坑EM/NVT/NPT 的 .mdp 参数与四个崩溃原因4.1 能量最小化steep 算法与 emtol1000 的收敛判据体系加入溶剂和离子后初始结构可能因为原子重叠或水分子取向不合理而存在很高的能量如果直接开始动力学模拟这些热点会让积分器失效。所以第一步先做能量最小化EM让体系从局部高位点滑到最近的势能谷底。典型 mdp 文件长这样integrator steep emtol 1000.0 emstep 0.01 nsteps 50000 cutoff-scheme Verlet coulombtype PME rcoulomb 1.0 rvdw 1.0 pbc xyz递推逻辑integrator 用 steep最陡下降法是 EM 的默认选项它沿负梯度方向迭代虽然收敛慢但非常稳健emtol 是停止条件当体系最大受力小于 1000.0 kJ/mol/nm 时就认为已收敛这个标准是官方教程的预设值emstep 是最小化步长0.01 nm 很保守主要防止步长跨越能量尖峰nsteps 给了 50000 步上限防止极端体系停不下来。跑 EM 的命令是gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr gmx mdrun -v -deffnm em跑完检查 em.log 里最后几行的收敛信息。判断 EM 是否真正成功有两条一是 log 文件没有报错二是 em.gro 原子坐标没有异常重叠。如果 nsteps 用完还没收敛通常是初始结构有硬伤比如溶剂分子被挤进了蛋白质内部或者盒子太小导致原子跟周期镜像重叠。4.2 NVT 平衡V-rescale 温控与位置限制的配合EM 之后体系能量已经在势阱里但温度还没有意义因为粒子的速度分布是人为指定的。NVT 平衡要做的事情是让温度达到目标值并稳定下来同时对蛋白质重原子施加位置限制POSRES免得它因为溶剂环境突变而剧烈构象变化。NVT 平衡的 mdp 关键参数define -DPOSRES integrator md nsteps 50000 dt 0.002 constraints h-bonds constraint_algorithm lincs tcoupl V-rescale tc-grps Protein Non-Protein tau_t 0.1 0.1 ref_t 300 300 pcoupl no gen_vel yes gen_temp 300值得解释的参数-DPOSRES 这个 define 会让 pdb2gmx 生成的 posre.itp 生效把蛋白质重原子拴在初始位置配体如果没有单独的 posre_ligand.itp需要先跑一遍 genrestr 再把它 include 进 topV-rescale 是扩展的 Berendsen 温控器它比原版 Berendsen 多一个随机项能正确还原动能涨落tc-grps 分成 Protein 和 Non-Protein 两组让蛋白质和溶剂分别耦合到 300 K比整个体系一组更接近真实行为tau_t 是温度耦合时间常数0.1 ps 是常规值gen_vel yes 从 Maxwell 分布随机生成初始速度gen_seed 设成 -1 的话每次跑都会得到不同初始速度想精确复现某个模拟就固定成一个正整数。执行命令gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr gmx mdrun -v -deffnm nvt-c 指定输入坐标-r 指定位置限制参考结构通常就和 -c 用同一个文件-t 在续跑时需要第一步平衡不用加。NVT 跑完后检查温度曲线是否爬升后稳定在 300 K 附近上下浮动几个 K 都正常。4.3 NPT 平衡Parrinello-Rahman 恒压器与 continuationyesNVT 稳定了温度接下来要稳定压力。NPT 系综里盒子体积可变密度和压力关系密切这一步的 mdp 参数变化集中在压力耦合块和安全重启设置continuation yes gen_vel no pcoupl Parrinello-Rahman pcoupltype isotropic tau_p 2.0 ref_p 1.0 compressibility 4.5e-5continuationyes 表示从 NVT 阶段继续此时不能再从 Maxwell 分布生成速度要靠 nvt.cpt 检查点文件把此前的位置和速度读进来所以要加 -t 参数。Parrinello-Rahman 是一个能正确产生 NPT 系综的恒压器但在模拟初期体系离平衡太远时容易振荡所以正规流程是在 NPT 先用 Berendsen 恒压器跑一段再切换到 PR不过 100 ps 直接从 PR 开始一般也能扛住。tcoupl 沿用 V-rescale 不变温度组别和 NVT 保持一致。gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr gmx mdrun -v -deffnm npt跑完后用 gmx energy 看压力曲线目标是在 1 bar 附近稳定波动。注意 NPT 检查盒子的尺寸变化来自 log 文件里的 box 记录如果体积变化超过百分之几说明之前的盒子尺寸设置不合理。4.4 四个崩溃原因排查初始态、收敛、盒子、耦合参数笔记里明确记录了 NVT 阶段体系崩溃的四个原因这也是整个平衡流程最容易翻车的集中区。按实际踩坑频率排个序第一个坑是初始态不合理。最典型的表现是体系初始化后原子重叠或者水分子插进蛋白质内部导致 LINCS 警告刷屏、能量爆到 10 的 6 次方以上。解决手段是回到 EM 之前的 gro 文件重新可视化检查溶剂和溶质是否有穿透检查 -d 1.0 是否加够了如果配体结合腔里有水分子遗留也要手动删除或重排。第二个坑是 EM 没有真正收敛就强行进了 NVT。判断方法是看 em.log 末尾的 Steepest Descents 是否停在 Converged to machine precision但机器精度严格收敛对复杂体系本来就是奢求临床上更常用的是看最大受力是否降到 emtol 以下。如果 nsteps 用完了但受力还是非常大优先检查是否有原子类型参数缺失比如配体原子类型没在主力场里定义全。第三个坑是盒子大小。周期性边界条件下盒子边长必须大于两倍截断半径加溶质直径如果 -d 设太小或者蛋白质本身很大周期镜像之间会发生原子重叠。检查方法跑 NVT 前用 editconf 输出的 gro 量一次蛋白质最长对角线确保盒子边长减掉蛋白质直径后仍大于 2.0 nm。第四个坑是温度耦合参数设置不合理tau_t 取值过大导致控温跟不上或者 ref_t 设定和 gen_temp 不一致导致初始冲击太大。步长也要注意常规 2 fs 时必须开启 h-bonds 约束如果没有 constraints h-bonds2 fs 的步长对水分子键振动来说太大了体系很快发散。碰到这种情况先改成 1 fs 跑一小段看能不能稳住能稳住再考虑逐步加回 2 fs。5. 正式 MD 与结果验证从 mdrun 到 RMSD/RMSF/氢键分析5.1 生产模拟的 .mdp 调整关 TRR、压缩坐标输出频率NPT 结束后体系在目标温度和压力下已经平衡接下来释放位置限制跑生产模拟收集数据。生产模拟的 mdp 和 NPT 很接近几个关键差异需要注意。第一define 里不再有 -DPOSRES如果想对配体或者某些 loop 区继续限制也可以局部保留。第二nsteps 按目标时长换算比如 2 fs 步长跑 100 ns 就是 50000000 步但通常分多个段跑更稳妥。第三输出控制的差异很关键nstxout 0 nstvout 0 nstfout 0 nstxout-compressed 5000 nstlog 5000 nstenergy 5000trr 文件是完整精度的坐标轨迹每步都写体积膨胀极快生产模拟里一般关掉只保留压缩轨迹文件xtcnstxout-compressed 5000 表示每 10 ps 存一帧100 ns 就是 10000 帧对 RMSD、RMSF、氢键分析完全够用。体系的能量和 log 记录频率是 5000 步一次对应的也是 10 ps 间隔。生产模拟的命令是gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md_0_1.tpr gmx mdrun -v -deffnm md_0_1注意这里多了 -n index.ndx 参数。正式模拟阶段索引文件比之前更重要因为分析和耦合分组都要用它来定义组。如果你在 genion 之后手动修改过 ndx 里的离子组归属这里更要用更新过的索引。5.2 gmx energy 验证温度、压力、势能是否在平台期模拟结束后第一件事不是急着看蛋白质动了多少而是检查轨迹的热力学性质是否正常。这个环节需要的能量文件是 .edr命令如下gmx energy -f md_0_1.edr -o temperature.xvg程序会列出文件里所有能量项输入编号或名称选择 Temperature输出 xvg 后用绘图工具画出温度随时间的变化。如果温度稳定在 300 K 附近做小幅热涨落说明体系没有异常发热或冷却。同样方法再检查 PressureNPT 段、Potential势能。势能曲线如果出现明显的台阶式跳跃通常意味着体系发生了不希望的构象跃迁或原子穿透。如果温度曲线偏离目标 10 K 以上检查是不是恒温器耦合组没设对或者水盒子在 NPT 阶段发生了过度膨胀/收缩导致密度异常。5.3 RMSD/RMSF/氢键分析判断复合物到底稳不稳确认热力学正常后接下来看结构稳定性。RMSD 衡量的是蛋白质主链随时间相对初始结构的偏离程度gmx rmsd -s md_0_1.tpr -f md_0_1.xtc -o rmsd.xvg程序会让你选择计算 RMSD 的原子组和参考组通常选 Backbone 作为分析组和参考组配体可以单独选出来算配体 RMSD。如果蛋白质 RMSD 在模拟后段趋于平台比如 0.2~0.3 nm 附近波动说明构象已经稳定如果持续上升没有收敛迹象可能模拟时长不够或者体系本身存在较大构象漂移。RMSF 则是按残基逐一计算平均涨落幅度gmx rmsf -s md_0_1.tpr -f md_0_1.xtc -o rmsf.xvgRMSF 数值高的区域通常是柔性 loop 或末端结合位点周围的残基如果波动很大说明配体结合不稳定。氢键分析则直接关注蛋白质和小分子间的相互作用持久性gmx hbond -s md_0_1.tpr -f md_0_1.xtc -num hbond.xvg选择蛋白质组和配体组后程序会统计随时间变化的氢键数量。挑几个能量最低的氢键用 gmx hbond 结合距离和角度过滤条件再看它们在模拟时间内的占有率就能判断相互作用到底是稳定存在还是偶发接触。这套笔记对生产模拟后的处理写得比较精简我按自己的习惯补一个标准动作每一段模拟结束后强制跑一遍「温度-势能-压力-RMSD」四件套全部正常才继续跑下一段任何一个指标异常就停下排查。早期跑一个上百纳秒的膜蛋白模拟时因为偷懒没看势能曲线跑了十几纳秒才发现体系里有一个水分子被挤进了蛋白疏水核心整个轨迹作废重跑。从那以后我每次拿到一段新轨迹都强制走一遍这四步验证流程确认没有热力学异常和构象漂移才进入分析阶段。希望这套预处理到验证的完整流程能帮你省下一些弯路。本文还有配套的精品资源点击获取