这些年带过不少学生和合作者从实验转来做计算我发现大家拿到一个蛋白-配体复合物最常问的一句话就是“我们做个分子动力学模拟验证一下结合稳定性”。但真正打开Gromacs面对一堆mdp文件、力场、拓扑和不断报错的命令行很多人第一周基本是懵的。我自己当初从纯实验背景转到计算模拟完整跑通第一轮蛋白-配体动力学模拟加上结果分析前前后后折腾了将近两周其中大部分时间花在配体拓扑生成和一堆莫名其妙的报错上。这篇教程不打算讲太多理论公式而是把我现在做项目时一定会走的完整流程按实际操作顺序写出来从结构文件准备、配体拓扑生成、体系搭建、能量最小化、NVT/NPT平衡到成品模拟和轨迹分析。每一步都附上可复制的命令和参数同时解释“为什么要这么做”以及哪些地方容易翻车。适合刚接触Gromacs、准备拿自己的蛋白-配体体系练手的人也适合已经跑过标准流程但想在结果分析里多挖点信息的人。1. 模拟前先想清楚这套体系要跑多久、用什么力场1.1 计算成本估算2 fs时间步长和100 ns模拟意味着什么很多人第一步就直接下载结构、写mdp、开跑跑到一半发现算得太慢或者资源不够再来调整体系大小和模拟时长非常被动。我先养成的习惯是拿到体系先数原子估算一下需要多少核时再决定跑多长。分子动力学的时间步长一般设成2 fs0.002 ps。为什么是2 fs而不是更大因为生物分子里最快的运动是C-H、O-H这类含氢键的伸缩振动周期在10 fs量级。时间步长必须远小于最快振动周期否则积分会发散。Gromacs里用LINCS算法把所有含氢键的伸缩约束住才能把步长从1 fs提高到2 fs。再往上加哪怕用更大的步长体系能量也会迅速漂移轨迹基本不可用。所以模拟时长的计算很简单100 ns 100000 ps 50000000步。也就是说如果你定义nsteps 50000000dt 0.002跑完就是100 ns。至于需要多少天取决于体系原子数和硬件。一套大约4万个原子蛋白300个残基加配体加水和离子的体系用单张主流GPU跑大约每天能跑100到200 ns如果只有CPU多核同样体系每天大概10到30 ns。这个差别直接决定了你的模拟方案是否可行。我一般建议先跑一个20到30 ns的预实验看RMSD曲线是否趋于平稳温度和势能是否稳定再决定要不要跑更长。不要一上来就提交200 ns的任务否则很可能跑了三天发现平衡阶段就有问题全部浪费。1.2 力场和水模型选型蛋白和配体必须“一家亲”力场是分子动力学里最核心的“游戏规则”。Gromacs本身不自带力场参数它只是按你指定的力场去读取参数文件。蛋白常用的有AMBER99SB-ILDN、CHARMM36m、OPLS-AA等配体则通常用GAFF配合AMBER力场或CGenFF配合CHARMM力场。这里最容易犯的错误是蛋白用AMBER力场配体却拿CGenFF参数强行混在一起。Gromacs语法上不一定报错因为原子类型名称不同但两类参数的数学形式和加和规则存在差异能量算出来就是四不像。我的原则是蛋白和配体必须来自同一个“家族”。具体来说两条推荐路线蛋白用AMBER99SB-ILDN配体用GAFF通过ACPYPE或antechamber生成拓扑蛋白用CHARMM36m配体用CGenFF通过CGenFF服务器生成拓扑。水模型也要和力场匹配。AMBER系列通常配TIP3P水CHARMM36m同样常配TIP3P你可以在pdb2gmx时指定-water tip3p。如果选了SPC/E或OPC水模型原则上也可以但需要确认配套的离子参数新手没必要在这个环节给自己挖坑直接选默认搭配最省心。另外如果你的蛋白有金属离子、二硫键、特殊修饰残基选力场时要特别注意CHARMM36m内置了较多翻译后修饰残基AMBER这边往往需要手动处理。所以带金属蛋白或修饰蛋白我默认优先考虑CHARMM36m路线。2. 结构文件准备蛋白、配体、复合物的三个入口2.1 蛋白结构从PDB下来后要处理什么从RCSB PDB数据库下载结构后不要直接丢进pdb2gmx。PDB文件里通常包含很多干扰信息晶体水分子、缓冲液分子、金属离子、多个构象、缺失残基等。pdb2gmx遇到不认识的残基名或重复原子轻则警告重则直接报错退出。我的标准处理步骤是在PyMOL或VMD里打开PDB文件看看蛋白有几条链、有没有缺失的loop区域、配体名称是什么。删除晶体水分子和所有非必须的异源分子。注意如果某个离子比如锌离子是蛋白功能必需的要保留因为去掉后蛋白结构可能在模拟中不稳定。如果结构有缺失残基在PDB文件里会显示为空白或REMARK 465记录。残基缺失比较长时建议用Modeller或AlphaFold补全缺失只有几个残基且位于表面loop有时可以直接保留缺失状态继续模拟但你需要明白这可能影响局部构象。用文本编辑器或grep看一下蛋白链末尾有没有TER标记多链蛋白链间是否清楚。pdb2gmx命令会做加氢、分配质子化状态、生成拓扑的工作。加氢这里有个坑如果PDB里已经带了氢有些高分辨率结构有建议加-ignh参数忽略原有氢原子让pdb2gmx按力场规则重新加氢这是最稳妥的做法。2.2 配体结构3D构象、质子化状态和坐标来源配体结构通常从PubChem下载但PubChem的2D结构的SDF里没有3D坐标。推荐的做法是从PubChem拿到SMILES或SDF后用OpenBabel生成3D构象并加氢再用Gaussian或xtb做一次半经验/DFT几何优化获得低能构象。这一步不是必须的但配体的初始构象如果是一个高能折叠状态模拟前几百ps可能会发生明显的构象重排影响你对结合稳定性的判断。质子化状态非常关键。配体有个可电离的羧基生理pH 7.4下大概率是去质子化的如果是氨基则大概率是质子化的。用OpenBabel默认状态经常不准。我会在配体准备阶段就明确pH条件必要时用ChemAxon Marvin或者RDKit的pKa工具判断。另外如果你是用分子对接得到复合物坐标那么配体的3D构象已经被对接软件处理过一般不需要重新生成3D结构直接沿用对接结果里的配体坐标和质子化状态即可。但要注意两点对接软件产生的配体构象可能带有不合理的键长键角最好也用半经验方法快速优化一下以及对接结果文件的原子命名和拓扑生成时的命名要保持一致。2.3 复合物坐标对接结果和晶体结构的区别晶体结构如果有共晶配体那是最好的起点因为实验坐标直接告诉你配体在蛋白里的结合姿态。但很多情况下你的配体是全新的分子没有共晶结构只能依赖对接。对接结果需要注意“打分最高不代表结合模式正确”。我一般会看对接pose在已知关键残基附近有没有形成氢键或疏水接触比如某个已知的催化残基或突变实验证实影响活性的残基。如果配体和这些残基没有接触即使对接分数很高我也会多保留几个pose做模拟分别看稳定性。从晶体结构出发时复合物坐标通常直接用晶体里的配体坐标。但有些晶体结构里的配体与蛋白之间有较多空间冲突或极近接触这和晶体分辨率、对称分子有关。可以先在PyMOL里拉一下距离手动调整冲突原子避免模拟第一步就因斥力爆炸。3. 配体拓扑生成整个流程最容易翻车的地方3.1 为什么pdb2gmx搞不定配体pdb2gmx的本质是从力场的残基库中寻找你提供的氨基酸或核酸序列为它们加氢、生成键连关系、输出拓扑。它对“非标准残基”——配体小分子——根本无能为力通常直接跳过或报错。原因是力场参数是一套完整的势函数参数集每个原子有原子类型、电荷每对原子之间有键伸缩参数、键角参数、二面角参数还有非键相互作用参数。氨基酸的标准参数被预先存在力场文件里pdb2gmx只需要查表。而配体是任意分子没有现成条目必须用另一套工具为它“量身定做”全套参数。这就是配体拓扑生成。配体拓扑生成工具有很多常用的三条路线CHARMM-GUI的Ligand Reader Modeler输出CGenFF参数支持Gromacs格式CGenFF在线服务器加cgenff_charmm2gmx.py脚本输出Gromacs可直接使用的itp文件ACPYPE配合Antechamber基于GAFF力场生成Gromacs拓扑。下面展开两条最常用的路线。3.2 CGenFF在线服务器与cgenff_charmm2gmx.py的完整流程CGenFF路线和CHARMM36m蛋白力场天然匹配是我做蛋白-配体模拟的默认选择。流程如下准备一个配体结构的PDB或MOL2文件可以带氢也可以不带。确保原子命名规范。打开CGenFF服务器上传结构填好分子名称提交。服务器返回结果包含一个penalty值。这是参数可靠性评分penalty越低说明原子类型和参数越可信低于10通常认为很可靠10到50需要检查一下特殊二面角项高于50的话要非常谨慎最好手动检查甚至更换力场路线。下载结果中的.str文件和带氢的PDB文件。在本地用cgenff_charmm2gmx.py脚本转换python cgenff_charmm2gmx.py LIG LIG.str lig.pdb运行后会生成LIG.itp和一个包含配体坐标的LIG.pdb名字可能有版本差异。要注意cgenff_charmm2gmx.py这个脚本可以从GitHub上找不同版本对Python2和Python3的兼容性不同。我遇到过好几个朋友在Windows上跑这个脚本因为环境缺少numpy或者路径问题卡了半天。建议直接在Linux环境里用conda建一个Python环境来跑。成功生成itp后把LIG.itp里[ atoms ]段的原子顺序和配体PDB里的原子顺序一一对照确认没有错位。CGenFF脚本在原子排列上发生错乱的情况不常见但也确实出现过尤其是存在同分异构或对称基团时。3.3 ACPYPE/GAFF路线与两条路线的取舍如果你蛋白选的是AMBER力场配体那边对应的是GAFF力场通常用ACPYPE生成拓扑。流程大致是用antechamber把配体结构转成带GAFF原子类型的mol2文件并计算AM1-BCC电荷antechamber -i lig.pdb -fi pdb -o lig.mol2 -fo mol2 -c bcc -nc 0 -at gaff参数-nc配体的净电荷带-1电荷羧基就是-nc -1。生成frcmod文件parmchk2 -i lig.mol2 -f mol2 -o lig.frcmod。用tleap生成AMBER拓扑再转成Gromacs格式更省事的方式是直接装好acpype用acpype -i lig.mol2 -o gmx一步到位。ACPYPE的优势是自动生成Gromacs的itp和posre文件操作简单。虽然GAFF路线也能用但要注意AMBER力场现在很多新版蛋白力场使用ff19SB等配体使用GAFF2更匹配不同版本的组合需要注意电荷模型和原子类型的一致性。如果你只是学生交作业或者做初步筛选AMBER99SB-ILDN加GAFF的组合完全够用如果目标是发文章而且要求参数高度一致建议统一走CGenFF加CHARMM36m。3.4 生成后必做的三项检查原子名、电荷、残基名生成配体拓扑之后我强烈建议不要直接往下跑先做三个检查检查原子名是否和坐标文件完全对应。这是最常见的报错来源。Gromacs在gmx grompp的时候会比对拓扑里的原子和坐标文件里的原子名称、数量任何一个对不上都会报“Atom names in ... do not match”之类的错误。配体itp里的原子是带残基名的比如LIG而坐标文件里配体残基名也要改成LIG否则残基名不一致同样报错。检查电荷是否合理。CGenFF会给出部分电荷比如-0.834这种这是正常的Gromacs完全支持部分电荷。但如果某个碳原子带着-2.5的大电荷那基本说明参数分配出了问题需要回到CGenFF结果检查。AM1-BCC电荷同样可能出现非整数需要确认总电荷接近整数比如-1.006这种微小偏差是正常的。检查是否需要加#include和分子条目。在topol.top文件里需要加一行#include LIG.itp然后在[ molecules ]段里加一行LIG 1。很多新手忘了这两步导致配体电荷和原子坐标对不上。记得在加完拓扑后再跑一次gmx grompp确认没有警告和错误。4. 体系搭建从两组坐标到可以跑模拟的盒子4.1 pdb2gmx生成蛋白拓扑把配体坐标合回去配体拓扑单独生成后体系搭建的第一步是给蛋白生成拓扑gmx pdb2gmx -f protein.pdb -o protein_processed.gro -ff charmm36m -water tip3p -ignh执行后会提示选择二硫键配对如果有的话比如有两对Cys之间形成二硫键会列出编号让你选。这里把成对的Cys序号选上拓扑中会生成相应的二硫键约束。pdb2gmx会生成posre.itp、topol.top和去掉了配体的蛋白坐标。因为pdb2gmx不认识配体复合物文件直接丢进去会把配体丢掉。所以常规操作顺序是先单独处理蛋白再把配体坐标合回处理后的蛋白坐标文件。合回坐标有两种方式如果原始复合物是从晶体或对接得到先用PyMOL把蛋白和配体分别导出为pdb蛋白走完pdb2gmx后用PyMOL再把处理后的蛋白gro和配体pdb合并——注意PyMOL导出gro格式的兼容性一般不如直接把配体的坐标行追加到protein_processed.gro末尾同时把原子总数改成“蛋白原子数配体原子数”。更稳妥的方式是用Gromacs自带的gmx insert-molecules但这个方法适合随机放置配体如果你希望配体保持在结合位点应该手动合并坐标。我自己最常用的方式是写一个简单脚本把pdb2gmx输出gro文件最后一行原子总数修改并把配体部分坐标从配体PDB中提取的原子行追加进去然后手动编辑gro末尾的原子数。gro格式的原子坐标精度到小数点后三位有时候配体坐标精度不够合回后会有少量原子重叠但之后能量最小化会处理不必太过担心。4.2 盒子形状与加水菱形十二面体为何比立方体香体系盒子大小的常规要求是蛋白表面到盒子边界至少留1.0 nm这个距离对应的是非键截断半径1.2 nm加上缓冲避免蛋白镜像之间的周期性相互作用影响到真实部分。gmx editconf -f complex.gro -o complex_box.gro -bt dodecahedron -d 1.0为什么用菱形十二面体而不是立方体因为相同周期距离下菱形十二面体的体积比立方体少约30%相当于能省约30%的水分子直接减少溶剂自由度显著降低计算量。可溶蛋白体系我都习惯用dodecahedron只有膜蛋白或者特定各向异性体系才需要别的盒子类型。加水的命令gmx solvate -cp complex_box.gro -cs spc216.gro -o complex_solv.gro -p topol.topspc216.gro是Gromacs自带的一个已平衡水盒子同时适用于TIP3P、SPC等多种水模型里面没有氢原子位置冲突问题可以直接用。加完水后topol.top的[ molecules ]部分会自动追加大量SOL条目你不需要手动改水分子数目。4.3 添加反离子到生理浓度纯水体系里蛋白表面带电荷体系净电荷往往不为零。分子动力学模拟中人为引入净电荷会造成程度不等的静电问题而且PME算法在周期性盒子中处理净电荷也需要补偿背景电荷。所以必须加反离子中和体系同时如果需要模拟生理条件还要加NaCl到0.15 M。先把带离子的体系做一次grompp生成一个临时tpr文件再用genion把溶剂分子替换成离子gmx grompp -f ions.mdp -c complex_solv.gro -p topol.top -o ions.tpr gmx genion -s ions.tpr -o complex_ions.gro -p topol.top -pname NA -nname CL -neutral -conc 0.15执行genion时它会让你选择要被替换的原子组通常会列出SOL、Protein、LIG等。此时输入SOL表示从水分子中随机挑一部分替换成离子。系统会额外输出“Number of NA ions ... Number of CL- ions ...”你可以根据盒子体积大致估算浓度是否合理。净电荷不为0时-neutral会补足中和所需的反离子-conc 0.15则是在中和基础上额外加入生理浓度盐。ions.mdp这个文件不需要特殊内容实际上grompp只是需要一个合法的mdp文件来构建tpr可以用一个极简的空mdp或最小化mdp来跑。5. 能量最小化先让体系“坐下来”再谈模拟5.1 为什么体系在模拟第0步就会炸当你把蛋白、配体、水、离子拼成一个盒子后这个体系的初始结构几乎必然存在局部空间冲突。原因很简单PDB结构或对接结构来源于实验或算法原子的坐标本身就有不确定性而你在加水时水分子是根据网格随机放置的很可能有一个水分子正好落在蛋白侧链上。如果用这样的结构直接开始MD最大问题是力场里原子间排斥力随距离急剧上升两原子距离比平衡键长小0.1 nm时排斥能量可能高达几千甚至几万kJ/mol这个巨大的力会让原子在一步之内获得极高速度体系温度瞬间爆表模拟直接崩溃。能量最小化的目的就是在模拟开始前用最速下降等算法找到势能面的局部极小点把不合理的接触“推开”让最大受力降到可接受范围。5.2 em.mdp参数与收敛判断能量最小化mdp参数一般长这样integrator steep emtol 1000 emstep 0.01 nsteps 50000 cutoff-scheme Verlet nstlist 10 rlist 1.2 rcoulomb 1.2 rvdw 1.2 coulombtype PMEemtol 1000的单位是kJ/mol/nm意思是当体系中最大受力小于1000时最小化收敛。这个阈值对应能量梯度的判断等价于一个比较“宽松”但足够安全的收敛标准。emstep是最陡下降的初始步长取了0.01 nm如果算法发现能量上升会自动减半步长。执行gmx grompp -f em.mdp -c complex_ions.gro -r complex_ions.gro -p topol.top -o em.tpr gmx mdrun -deffnm em -v-r用于生成位置约束文件虽然最小化这里用不到但后面平衡需要posre文件所以建议加上。如果用了-r complex_ions.groposre.itp会在grompp阶段自动读进来。最小化完成后看em.log末尾的势能正常情况下势能应该是负的且收敛输出提示达到Fmax收敛。也可以用gmx energy -f em.edr -o potential.xvg选择Potential项画图。一个合理的能量最小化结果是从一个较高的正值或较小的负值开始快速下降到稳定的负值最后几乎不再变化。如果在能量轨迹里看到一千多甚至上万的势能且剧烈波动说明需要检查结构。5.3 最小化不收敛时的两条出路能量最小化不收敛通常是两种原因一是结构里有严重的原子重叠二是力场参数出了问题尤其是配体拓扑。第一种情况可以在最小化之前先用PyMOL或者在gro文件里检查配体和周围残基的距离。如果发现有原子距离小于0.1 nm先手动调整配体位置或删除个别重叠水分子。注意pdb2gmx和solvate添加的水分子编号可能在gro文件里难找你也可以直接用gmx select选距离蛋白太近的水分子输出编号手动删掉。第二种情况配体itp里有不合理的键长、键角或原子类型最小化时会有巨大内应力。可以尝试先把配体从体系中拿出来单独对配体做气相最小化看它能否正常收敛。如果配体单独也发不了最小化问题出在配体拓扑需要回到参数生成环节。还有一种临时方案是给配体加位置约束比如在em.mdp里define -DFLEXIBLE配合配体posre文件先把蛋白和水调整好再释放配体做第二轮最小化。这个方法虽然有效但也会掩盖一部分配体拓扑问题所以只建议作为排查手段。6. NVT/NPT平衡把体系带到300K和1bar6.1 NVT升温控温器选择与位置约束能量最小化后体系处于势能局部极小的低温状态需要升温到目标温度同时让蛋白骨架不能乱动给水分子和侧链充分松弛的时间。这个过程就是NVT平衡即保持原子数、体积、温度不变。NVT平衡的mdp核心参数integrator md dt 0.002 nsteps 25000 tcoupl V-rescale tc-grps System tau_t 0.1 ref_t 300 pcoupl no constraints h-bonds constraint_algorithm lincs continuation nonsteps 25000乘以dt 0.002就是50 ps。作为升温阶段这个长度足够让体系从初始温度到达300K并稳定下来。NVT平衡里最重要的是位置约束蛋白重原子被固定在初始位置只允许水、离子和蛋白侧链小范围移动。这么做是因为直接从最小化后的结构开始放开全部自由度蛋白可能会在温度上升的瞬间产生明显构象漂移尤其在溶剂还没完全平衡时这种漂移很容易把蛋白推向不合理的结构。位置约束通过define -DPOSRES启用。在topol.top里protein和配体通常会各有一个#ifdef POSRES的posre引入段。-DPOSRES会让蛋白重原子受到一个力常数为1000 kJ/mol/nm²的谐振子约束。NVT阶段加上这个定义跑50 ps后体系温度通常会稳定在300K附近动能和势能也会趋于平缓。6.2 NPT加压为什么要看密度是否收尾NVT平衡完成后体系温度已经正确但压力还不正确——因为在固定体积下溶剂密度未必对应1 bar压力下的密度。NPT平衡在保持温度的同时通过压力耦合调节盒子体积让体系密度和压力达到目标值。对于等温等压条件这更接近一个真实的生理环境。NPT平衡mdp参数integrator md dt 0.002 nsteps 500000 tcoupl V-rescale tc-grps System tau_t 1.0 ref_t 300 pcoupl Parrinello-Rahman pcoupltype isotropic tau_p 2.0 ref_p 1.0 compressibility 4.5e-5 constraints h-bonds continuation yesNPT平衡时长建议至少500 ps有的体系甚至需要1 ns才能让密度稳定。这里为什么tau_p 2.0而不是0.1压力耦合的响应时间比温度慢得多时间常数太短会导致盒子体积剧烈震荡密度曲线噪声很大。NPT平衡结束后检查gmx energy里的Density项水的理论密度约1000 kg/m³蛋白和配体的存在会让体系密度略高通常在1000到1100 kg/m³之间。如果密度在这个范围内且曲线在平衡后期基本水平说明体系已经达到正确的压力条件。如果密度还在缓慢上升或下降说明平衡不够要延长NPT时间。6.3 平衡成功的三个检查每次平衡结束我都会例行检查三样东西缺一不可温度曲线是否稳定在300K附近。如果温度持续偏离目标好几K说明控温参数或初始速度有问题先别急着跑生产。压力曲线是否在平均值附近波动。Parrinello-Rahman的输出压力噪声较大是正常的关键是看平均值是否在1 bar附近。如果压力一直在几十bar甚至上百bar波动盒子尺寸可能不合适要检查周期距离是否足够。密度是否收敛到一个稳定值。密度的收敛意味着溶剂已经被“压实”到正确的密度这是整个NPT阶段成功的核心指标。另外可以看一下势能曲线平衡后期势能不应该有明显上升趋势。7. 成品模拟参数文件、GPU加速与续跑7.1 production mdp参数逐行平衡结束后可以进行成品模拟。生产模拟的mdp和NPT平衡有很多相似之处但有一个关键点继续使用Parrinello-Rahman压力耦合时如果你是从NPT平衡的cpt文件续跑可以无缝衔接如果生产模拟单独从头开始则需要先跑一小段Berendsen压力耦合的预平衡。最稳妥的方式是直接续跑NPT平衡的检查点文件。一个常用的生产mdpintegrator md dt 0.002 nsteps 50000000 nstxout-compressed 5000 cutoff-scheme Verlet rlist 1.2 rcoulomb 1.2 rvdw 1.2 coulombtype PME fourierspacing 0.12 tcoupl V-rescale tc-grps System tau_t 1.0 ref_t 300 pcoupl Parrinello-Rahman pcoupltype isotropic tau_p 2.0 ref_p 1.0 compressibility 4.5e-5 constraints h-bonds continuation yes关键行解释nsteps 50000000对应100 ns模拟。nstxout-compressed 5000即每5000步也就是10 ps存一帧压缩轨迹。100 ns会产生1万帧文件体积和分辨率足够做绝大多数分析。如果只想看大趋势可以设为10000甚至20000但做氢键和距离分析时帧数太少会影响统计。continuation yes告诉Gromacs这是一段续跑不要重新生成初始速度。PME的fourierspacing 0.12表示倒空间格点间距约0.12 nm精度足够没必要更小否则会显著增加计算量。运行命令gmx grompp -f md.mdp -c npt.gro -t npt.cpt -r npt.gro -p topol.top -o md.tpr gmx mdrun -deffnm md -v7.2 GPU加速和速度评估Gromacs的性能优化做得非常出色支持GPU加速后可溶蛋白体系的模拟速度比纯CPU快一个数量级以上。在mdrun时指定gmx mdrun -deffnm md -v -nb gpu -pme gpu -bonded gpu -update gpu-nb是非键计算-pme是长程静电的PME计算-bonded是键合相互作用-update是坐标更新。如果你的GPU支持一般NVIDIA的Pascal架构以后都支持把四个都设为gpu可以让非键和键合计算都跑到GPU上。Gromacs的日志和终端会输出“GPU”相关性能数据例如Performance: 185 ns/day之类的指标。这里提醒一个容易忽略的问题-update gpu会要求GPU计算和CPU通信有较高带宽某些老显卡或者集群节点配置不当时反而更慢。可以先试跑几千步比较不同设置的性能再决定要不要用-update gpu。另外盒子里如果有大量水分子非键计算是主要瓶颈GPU提升非常明显。但如果体系很小比如只有几千个原子GPU的空闲等待可能吃掉提升这时纯CPU反而更稳。具体看性能测试结果。7.3 模拟中断后的续跑跑100 ns甚至更长的模拟中途集群重启、作业超时、硬件故障几乎是必然事件。不要慌Gromacs的检查点机制设计得很好。模拟中断后目录里有md.cpt文件续跑命令gmx mdrun -deffnm md -cpi md.cpt -s md.tpr -vGromacs会从检查点状态继续而不是从头开始。注意续跑时不要让两个mdrun同时操作同一个tpr和cpt文件否则会写坏检查点。集群环境下尤其要注意作业管理系统是否允许你在上次作业结束后再次提交。还有一个实用技巧运行期间定期把md.cpt备份到另一个目录因为如果磁盘满了cpt写不进去模拟会中止备份可以避免丢失太长的轨迹。模拟稳定跑完的标志是目录下出现md.gro、md.xtc和完整的md.log日志末尾会显示“Finished mdrun on node”之类的信息。8. 轨迹分析从PBC处理到RMSD/RMSF/氢键/自由能8.1 第一件事永远是去周期性与居中模拟是在周期性盒子中进行的分子在模拟中可能穿越盒子边界直接分析原始md.xtc会看到蛋白链“断裂”成两半配体位置在盒子两端跳跃。所有轨迹分析的第一步都是把周期性效应处理干净。推荐的处理方式gmx trjconv -s md.tpr -f md.xtc -o md_whole.xtc -pbc whole gmx trjconv -s md.tpr -f md_whole.xtc -o md_center.xtc -pbc mol -center第一条命令-pbc whole把每个分子写完整不跨越边界。第二条命令-pbc mol -center先把体系居中再把分子置于盒子中心。执行第二条时Gromacs会询问选择哪个组做居中一般选Protein因为以蛋白为中心来分析配体周围环境最合理。处理完成后所有后续分析都基于md_center.xtc不要再回头用原始轨迹。经验不足的人最容易在这里翻车RMSD曲线看起来大幅度跳变其实只是分子在盒子边界来回穿行造成的假象。8.2 RMSD和RMSF稳定性与残基柔性RMSD是最基础也是最重要的稳定性指标。它计算每个帧的蛋白骨架原子相对于参考结构的均方根位移。参考结构的选择有讲究如果想知道模拟过程中蛋白是否发生了大幅构象变化用能量最小化后的结构或PDB晶体结构作参考如果想知道模拟是否达到稳定平衡也可以把平衡之后的轨迹平均结构作为参考。常用命令gmx rms -s md.tpr -f md_center.xtc -o rmsd.xvg -tu ns交互选择时拟合组选Backbone计算组也选Backbone。-tu ns让x轴时间以ns为单位方便作图。RMSD曲线通常在初始阶段快速上升然后在某个平台值附近波动。平台值因蛋白大小而异小蛋白可能在0.1到0.2 nm大蛋白或柔性结构在0.3到0.5 nm都很常见。如果RMSD在模拟后期仍然持续上升且没有平台趋势说明体系没有真正稳定可能需要重新评估模拟条件和时长。RMSF反映每个残基在模拟中的波动幅度gmx rmsf -s md.tpr -f md_center.xtc -o rmsf.xvg -res-res会把多个原子的RMSF按残基平均得到一个“每残基一个值”的序列方便绘制曲线。看RMSF时重点关注结合位点残基如果配体结合口袋附近残基的RMSF明显高于蛋白整体平均水平可能说明配体结合得不够稳定结合pose在模拟中有明显晃动如果口袋残基RMSF较低说明结合模式比较刚性是一个稳定结合的信号。8.3 配体结合相关的分析氢键、距离、接触RMSD和RMSF不够直接回答“配体结合得稳不稳”还需要看分子层面的相互作用。氢键分析gmx hbond -s md.tpr -f md_center.xtc -num hbond_num.xvg -life hbond_life.xvg交互选择两个组一个是Protein另一个是配体组LIG。-num输出的是每个时间点的氢键数目随时间的变化曲线-life输出氢键存在时间的分布。平均氢键数量越多、寿命越长说明配体和蛋白之间的极性相互作用越稳定。注意氢键的定义默认是供体-受体距离小于0.35 nm且角度偏离直线小于30度这个标准可以在mdp参数里调整但对新手先用默认即可。关键距离监控gmx distance -s md.tpr -f md_center.xtc -select atomname ...-select的语法需要指定具体的原子。比如配体的某个氧原子和蛋白某个残基侧链氮原子之间的距离。这类“看特定距离是否稳定地维持在一定范围”的分析比看一大堆氢键更有针对性尤其是已知关键残基时。选定原子编号可以从gmx select交互式命令行查询。接触面积分析可以用gmx mindist或第三方工具如PyMOL的contact命令统计配体原子在0.3到0.4 nm范围内接触到的蛋白残基。这个指标能直观说明配体周围有没有形成稳定的疏水口袋。8.4 结合自由能MM/PBSA的跨界提醒如果想更进一步估算结合自由能主流的做法是用MM/PBSA方法。Gromacs本身不带现成的MM/PBSA实现现在常用第三方工具gmx_MMPBSA。基本流程用复合物tpr文件和轨迹文件分别计算复合物、受体和配体在气相中的分子力学能量然后用隐式溶剂模型估算溶剂化自由能二者相减得到结合自由能。命令大致是gmx_MMPBSA -O -i mmpbsa.in -cs complex.tpr -ct complex.xtc -cg 1 13 -cp complex.top ...-cg指定受体和配体分别对应哪两个index组-cp提供复合物的拓扑文件。实际操作中坑非常多轨迹里必须包含蛋白和配体的所有原子而topol.top需要同时包括蛋白和配体参数受体和配体的坐标要从原始轨迹中提取并分别生成tpr。建议初学者先把前面的基础分析做好再考虑MM/PBSA因为算出来的绝对值对参数和方法选择非常敏感。一个重要提醒MM/PBSA给出的结合自由能是“能量估计”不是实验结合能。它的价值在于比较一系列相似配体或不同突变体之间的相对趋势或者作为结合模式排序的参考不要过度解读绝对值。误差范围内你只能看出“强结合”“弱结合”这样的定性结论不要拿它和小数点后两位的实验IC50去对表。9. 跑完这个流程之后我踩过的坑和一些习惯第一批模拟能顺利跑完结果分析也画出了图很多人就觉得大功告成了。但作为已经跑过数不清个蛋白-配体体系的人我始终有几件事每次都会做一遍并且建议你也养成习惯。第一跑完立即回头看日志里的警告。md.log里会有大量注释有些是无关紧要的但如果出现“turned around, warning, LINCS warning”或者“Step ... too long”的字样意味着模拟中有局部原子受力过大或者约束算法出问题。哪怕模拟最终正常结束这些警告也会影响轨迹质量一定要回到对应时间段检查。第二把轨迹按时间段拆开看RMSD。不要只看整条轨迹的平均值要分早期、中期、后期分别看。一个常见的假象是整段RMSD曲线平均看起来挺平但前20 ns一直在爬升后20 ns趋于稳定。这时候如果你把前60 ns也算进“平衡后”区间配体结合分析会包含大量构象搜索阶段的数据结果会偏乐观或偏模糊。我的习惯是去掉前20%作为“预平衡”丢弃只分析后面的稳定段。第三选原子组时谨慎使用System组。很多分析默认可以选择System但你真正关心的可能是蛋白骨架或者是配体周围一定范围内的残基。用gmx make_ndx建好常用组比如Protein_Backbone、LIG、Protein_LIG、within 0.5 of LIG这类索引组能让后续分析高效很多也能避免选错组导致的数据无效。第四轨迹可视化一定要做。再多的数值分析都无法代替亲眼看一下配体在结合口袋里到底怎么动。用VMD或PyMOL加载md_center.xtc和md.tpr慢速播放最后20 ns的轨迹。我曾经遇到过RMSD、氢键、接触面积全部正常的体系最后盯着轨迹看才发现配体在模拟后期旋转了180度原来结合模式早已改变。数值指标只是帮你定位问题视觉确认才能让你真正理解体系行为。第五保存好所有输入文件和参数版本。Gromacs版本更新很快同一套mdp在不同大版本下的默认值可能有差异。我在每个项目目录下都会保留一份version.txt记录Gromacs版本、mdp文件、力场版本和PDB编号。这样几个月后回来看结果还能复现当初的完整流程。新手尤其容易忽略这一点等你论文返修要求补跑模拟时发现自己忘了当初用的参数才是最痛苦的。最后说一句实操体会蛋白-配体动力学模拟这个流程本质上是一项需要反复调试和验证的技术活第一次跑通花了快两周第二次可能只要半天第五次的时候我已经能在一个小时内完成从原始结构到成品轨迹的全部准备。关键在于每跑完一个体系都把当时的命令、参数和踩坑记录整理成自己的checklist。等你的checklist越积越多就会发现Gromacs并没有想象中那么神秘——它只是要求你在每个环节都足够细心和诚实。 
长河君 · 2026/10/3 4:17:59 · 阅读 1.2k Gromacs蛋白-配体模拟实操教程：从拓扑生成到轨迹分析
这些年带过不少学生和合作者从实验转来做计算我发现大家拿到一个蛋白-配体复合物最常问的一句话就是“我们做个分子动力学模拟验证一下结合稳定性”。但真正打开Gromacs面对一堆mdp文件、力场、拓扑和不断报错的命令行很多人第一周基本是懵的。我自己当初从纯实验背景转到计算模拟完整跑通第一轮蛋白-配体动力学模拟加上结果分析前前后后折腾了将近两周其中大部分时间花在配体拓扑生成和一堆莫名其妙的报错上。这篇教程不打算讲太多理论公式而是把我现在做项目时一定会走的完整流程按实际操作顺序写出来从结构文件准备、配体拓扑生成、体系搭建、能量最小化、NVT/NPT平衡到成品模拟和轨迹分析。每一步都附上可复制的命令和参数同时解释“为什么要这么做”以及哪些地方容易翻车。适合刚接触Gromacs、准备拿自己的蛋白-配体体系练手的人也适合已经跑过标准流程但想在结果分析里多挖点信息的人。1. 模拟前先想清楚这套体系要跑多久、用什么力场1.1 计算成本估算2 fs时间步长和100 ns模拟意味着什么很多人第一步就直接下载结构、写mdp、开跑跑到一半发现算得太慢或者资源不够再来调整体系大小和模拟时长非常被动。我先养成的习惯是拿到体系先数原子估算一下需要多少核时再决定跑多长。分子动力学的时间步长一般设成2 fs0.002 ps。为什么是2 fs而不是更大因为生物分子里最快的运动是C-H、O-H这类含氢键的伸缩振动周期在10 fs量级。时间步长必须远小于最快振动周期否则积分会发散。Gromacs里用LINCS算法把所有含氢键的伸缩约束住才能把步长从1 fs提高到2 fs。再往上加哪怕用更大的步长体系能量也会迅速漂移轨迹基本不可用。所以模拟时长的计算很简单100 ns 100000 ps 50000000步。也就是说如果你定义nsteps 50000000dt 0.002跑完就是100 ns。至于需要多少天取决于体系原子数和硬件。一套大约4万个原子蛋白300个残基加配体加水和离子的体系用单张主流GPU跑大约每天能跑100到200 ns如果只有CPU多核同样体系每天大概10到30 ns。这个差别直接决定了你的模拟方案是否可行。我一般建议先跑一个20到30 ns的预实验看RMSD曲线是否趋于平稳温度和势能是否稳定再决定要不要跑更长。不要一上来就提交200 ns的任务否则很可能跑了三天发现平衡阶段就有问题全部浪费。1.2 力场和水模型选型蛋白和配体必须“一家亲”力场是分子动力学里最核心的“游戏规则”。Gromacs本身不自带力场参数它只是按你指定的力场去读取参数文件。蛋白常用的有AMBER99SB-ILDN、CHARMM36m、OPLS-AA等配体则通常用GAFF配合AMBER力场或CGenFF配合CHARMM力场。这里最容易犯的错误是蛋白用AMBER力场配体却拿CGenFF参数强行混在一起。Gromacs语法上不一定报错因为原子类型名称不同但两类参数的数学形式和加和规则存在差异能量算出来就是四不像。我的原则是蛋白和配体必须来自同一个“家族”。具体来说两条推荐路线蛋白用AMBER99SB-ILDN配体用GAFF通过ACPYPE或antechamber生成拓扑蛋白用CHARMM36m配体用CGenFF通过CGenFF服务器生成拓扑。水模型也要和力场匹配。AMBER系列通常配TIP3P水CHARMM36m同样常配TIP3P你可以在pdb2gmx时指定-water tip3p。如果选了SPC/E或OPC水模型原则上也可以但需要确认配套的离子参数新手没必要在这个环节给自己挖坑直接选默认搭配最省心。另外如果你的蛋白有金属离子、二硫键、特殊修饰残基选力场时要特别注意CHARMM36m内置了较多翻译后修饰残基AMBER这边往往需要手动处理。所以带金属蛋白或修饰蛋白我默认优先考虑CHARMM36m路线。2. 结构文件准备蛋白、配体、复合物的三个入口2.1 蛋白结构从PDB下来后要处理什么从RCSB PDB数据库下载结构后不要直接丢进pdb2gmx。PDB文件里通常包含很多干扰信息晶体水分子、缓冲液分子、金属离子、多个构象、缺失残基等。pdb2gmx遇到不认识的残基名或重复原子轻则警告重则直接报错退出。我的标准处理步骤是在PyMOL或VMD里打开PDB文件看看蛋白有几条链、有没有缺失的loop区域、配体名称是什么。删除晶体水分子和所有非必须的异源分子。注意如果某个离子比如锌离子是蛋白功能必需的要保留因为去掉后蛋白结构可能在模拟中不稳定。如果结构有缺失残基在PDB文件里会显示为空白或REMARK 465记录。残基缺失比较长时建议用Modeller或AlphaFold补全缺失只有几个残基且位于表面loop有时可以直接保留缺失状态继续模拟但你需要明白这可能影响局部构象。用文本编辑器或grep看一下蛋白链末尾有没有TER标记多链蛋白链间是否清楚。pdb2gmx命令会做加氢、分配质子化状态、生成拓扑的工作。加氢这里有个坑如果PDB里已经带了氢有些高分辨率结构有建议加-ignh参数忽略原有氢原子让pdb2gmx按力场规则重新加氢这是最稳妥的做法。2.2 配体结构3D构象、质子化状态和坐标来源配体结构通常从PubChem下载但PubChem的2D结构的SDF里没有3D坐标。推荐的做法是从PubChem拿到SMILES或SDF后用OpenBabel生成3D构象并加氢再用Gaussian或xtb做一次半经验/DFT几何优化获得低能构象。这一步不是必须的但配体的初始构象如果是一个高能折叠状态模拟前几百ps可能会发生明显的构象重排影响你对结合稳定性的判断。质子化状态非常关键。配体有个可电离的羧基生理pH 7.4下大概率是去质子化的如果是氨基则大概率是质子化的。用OpenBabel默认状态经常不准。我会在配体准备阶段就明确pH条件必要时用ChemAxon Marvin或者RDKit的pKa工具判断。另外如果你是用分子对接得到复合物坐标那么配体的3D构象已经被对接软件处理过一般不需要重新生成3D结构直接沿用对接结果里的配体坐标和质子化状态即可。但要注意两点对接软件产生的配体构象可能带有不合理的键长键角最好也用半经验方法快速优化一下以及对接结果文件的原子命名和拓扑生成时的命名要保持一致。2.3 复合物坐标对接结果和晶体结构的区别晶体结构如果有共晶配体那是最好的起点因为实验坐标直接告诉你配体在蛋白里的结合姿态。但很多情况下你的配体是全新的分子没有共晶结构只能依赖对接。对接结果需要注意“打分最高不代表结合模式正确”。我一般会看对接pose在已知关键残基附近有没有形成氢键或疏水接触比如某个已知的催化残基或突变实验证实影响活性的残基。如果配体和这些残基没有接触即使对接分数很高我也会多保留几个pose做模拟分别看稳定性。从晶体结构出发时复合物坐标通常直接用晶体里的配体坐标。但有些晶体结构里的配体与蛋白之间有较多空间冲突或极近接触这和晶体分辨率、对称分子有关。可以先在PyMOL里拉一下距离手动调整冲突原子避免模拟第一步就因斥力爆炸。3. 配体拓扑生成整个流程最容易翻车的地方3.1 为什么pdb2gmx搞不定配体pdb2gmx的本质是从力场的残基库中寻找你提供的氨基酸或核酸序列为它们加氢、生成键连关系、输出拓扑。它对“非标准残基”——配体小分子——根本无能为力通常直接跳过或报错。原因是力场参数是一套完整的势函数参数集每个原子有原子类型、电荷每对原子之间有键伸缩参数、键角参数、二面角参数还有非键相互作用参数。氨基酸的标准参数被预先存在力场文件里pdb2gmx只需要查表。而配体是任意分子没有现成条目必须用另一套工具为它“量身定做”全套参数。这就是配体拓扑生成。配体拓扑生成工具有很多常用的三条路线CHARMM-GUI的Ligand Reader Modeler输出CGenFF参数支持Gromacs格式CGenFF在线服务器加cgenff_charmm2gmx.py脚本输出Gromacs可直接使用的itp文件ACPYPE配合Antechamber基于GAFF力场生成Gromacs拓扑。下面展开两条最常用的路线。3.2 CGenFF在线服务器与cgenff_charmm2gmx.py的完整流程CGenFF路线和CHARMM36m蛋白力场天然匹配是我做蛋白-配体模拟的默认选择。流程如下准备一个配体结构的PDB或MOL2文件可以带氢也可以不带。确保原子命名规范。打开CGenFF服务器上传结构填好分子名称提交。服务器返回结果包含一个penalty值。这是参数可靠性评分penalty越低说明原子类型和参数越可信低于10通常认为很可靠10到50需要检查一下特殊二面角项高于50的话要非常谨慎最好手动检查甚至更换力场路线。下载结果中的.str文件和带氢的PDB文件。在本地用cgenff_charmm2gmx.py脚本转换python cgenff_charmm2gmx.py LIG LIG.str lig.pdb运行后会生成LIG.itp和一个包含配体坐标的LIG.pdb名字可能有版本差异。要注意cgenff_charmm2gmx.py这个脚本可以从GitHub上找不同版本对Python2和Python3的兼容性不同。我遇到过好几个朋友在Windows上跑这个脚本因为环境缺少numpy或者路径问题卡了半天。建议直接在Linux环境里用conda建一个Python环境来跑。成功生成itp后把LIG.itp里[ atoms ]段的原子顺序和配体PDB里的原子顺序一一对照确认没有错位。CGenFF脚本在原子排列上发生错乱的情况不常见但也确实出现过尤其是存在同分异构或对称基团时。3.3 ACPYPE/GAFF路线与两条路线的取舍如果你蛋白选的是AMBER力场配体那边对应的是GAFF力场通常用ACPYPE生成拓扑。流程大致是用antechamber把配体结构转成带GAFF原子类型的mol2文件并计算AM1-BCC电荷antechamber -i lig.pdb -fi pdb -o lig.mol2 -fo mol2 -c bcc -nc 0 -at gaff参数-nc配体的净电荷带-1电荷羧基就是-nc -1。生成frcmod文件parmchk2 -i lig.mol2 -f mol2 -o lig.frcmod。用tleap生成AMBER拓扑再转成Gromacs格式更省事的方式是直接装好acpype用acpype -i lig.mol2 -o gmx一步到位。ACPYPE的优势是自动生成Gromacs的itp和posre文件操作简单。虽然GAFF路线也能用但要注意AMBER力场现在很多新版蛋白力场使用ff19SB等配体使用GAFF2更匹配不同版本的组合需要注意电荷模型和原子类型的一致性。如果你只是学生交作业或者做初步筛选AMBER99SB-ILDN加GAFF的组合完全够用如果目标是发文章而且要求参数高度一致建议统一走CGenFF加CHARMM36m。3.4 生成后必做的三项检查原子名、电荷、残基名生成配体拓扑之后我强烈建议不要直接往下跑先做三个检查检查原子名是否和坐标文件完全对应。这是最常见的报错来源。Gromacs在gmx grompp的时候会比对拓扑里的原子和坐标文件里的原子名称、数量任何一个对不上都会报“Atom names in ... do not match”之类的错误。配体itp里的原子是带残基名的比如LIG而坐标文件里配体残基名也要改成LIG否则残基名不一致同样报错。检查电荷是否合理。CGenFF会给出部分电荷比如-0.834这种这是正常的Gromacs完全支持部分电荷。但如果某个碳原子带着-2.5的大电荷那基本说明参数分配出了问题需要回到CGenFF结果检查。AM1-BCC电荷同样可能出现非整数需要确认总电荷接近整数比如-1.006这种微小偏差是正常的。检查是否需要加#include和分子条目。在topol.top文件里需要加一行#include LIG.itp然后在[ molecules ]段里加一行LIG 1。很多新手忘了这两步导致配体电荷和原子坐标对不上。记得在加完拓扑后再跑一次gmx grompp确认没有警告和错误。4. 体系搭建从两组坐标到可以跑模拟的盒子4.1 pdb2gmx生成蛋白拓扑把配体坐标合回去配体拓扑单独生成后体系搭建的第一步是给蛋白生成拓扑gmx pdb2gmx -f protein.pdb -o protein_processed.gro -ff charmm36m -water tip3p -ignh执行后会提示选择二硫键配对如果有的话比如有两对Cys之间形成二硫键会列出编号让你选。这里把成对的Cys序号选上拓扑中会生成相应的二硫键约束。pdb2gmx会生成posre.itp、topol.top和去掉了配体的蛋白坐标。因为pdb2gmx不认识配体复合物文件直接丢进去会把配体丢掉。所以常规操作顺序是先单独处理蛋白再把配体坐标合回处理后的蛋白坐标文件。合回坐标有两种方式如果原始复合物是从晶体或对接得到先用PyMOL把蛋白和配体分别导出为pdb蛋白走完pdb2gmx后用PyMOL再把处理后的蛋白gro和配体pdb合并——注意PyMOL导出gro格式的兼容性一般不如直接把配体的坐标行追加到protein_processed.gro末尾同时把原子总数改成“蛋白原子数配体原子数”。更稳妥的方式是用Gromacs自带的gmx insert-molecules但这个方法适合随机放置配体如果你希望配体保持在结合位点应该手动合并坐标。我自己最常用的方式是写一个简单脚本把pdb2gmx输出gro文件最后一行原子总数修改并把配体部分坐标从配体PDB中提取的原子行追加进去然后手动编辑gro末尾的原子数。gro格式的原子坐标精度到小数点后三位有时候配体坐标精度不够合回后会有少量原子重叠但之后能量最小化会处理不必太过担心。4.2 盒子形状与加水菱形十二面体为何比立方体香体系盒子大小的常规要求是蛋白表面到盒子边界至少留1.0 nm这个距离对应的是非键截断半径1.2 nm加上缓冲避免蛋白镜像之间的周期性相互作用影响到真实部分。gmx editconf -f complex.gro -o complex_box.gro -bt dodecahedron -d 1.0为什么用菱形十二面体而不是立方体因为相同周期距离下菱形十二面体的体积比立方体少约30%相当于能省约30%的水分子直接减少溶剂自由度显著降低计算量。可溶蛋白体系我都习惯用dodecahedron只有膜蛋白或者特定各向异性体系才需要别的盒子类型。加水的命令gmx solvate -cp complex_box.gro -cs spc216.gro -o complex_solv.gro -p topol.topspc216.gro是Gromacs自带的一个已平衡水盒子同时适用于TIP3P、SPC等多种水模型里面没有氢原子位置冲突问题可以直接用。加完水后topol.top的[ molecules ]部分会自动追加大量SOL条目你不需要手动改水分子数目。4.3 添加反离子到生理浓度纯水体系里蛋白表面带电荷体系净电荷往往不为零。分子动力学模拟中人为引入净电荷会造成程度不等的静电问题而且PME算法在周期性盒子中处理净电荷也需要补偿背景电荷。所以必须加反离子中和体系同时如果需要模拟生理条件还要加NaCl到0.15 M。先把带离子的体系做一次grompp生成一个临时tpr文件再用genion把溶剂分子替换成离子gmx grompp -f ions.mdp -c complex_solv.gro -p topol.top -o ions.tpr gmx genion -s ions.tpr -o complex_ions.gro -p topol.top -pname NA -nname CL -neutral -conc 0.15执行genion时它会让你选择要被替换的原子组通常会列出SOL、Protein、LIG等。此时输入SOL表示从水分子中随机挑一部分替换成离子。系统会额外输出“Number of NA ions ... Number of CL- ions ...”你可以根据盒子体积大致估算浓度是否合理。净电荷不为0时-neutral会补足中和所需的反离子-conc 0.15则是在中和基础上额外加入生理浓度盐。ions.mdp这个文件不需要特殊内容实际上grompp只是需要一个合法的mdp文件来构建tpr可以用一个极简的空mdp或最小化mdp来跑。5. 能量最小化先让体系“坐下来”再谈模拟5.1 为什么体系在模拟第0步就会炸当你把蛋白、配体、水、离子拼成一个盒子后这个体系的初始结构几乎必然存在局部空间冲突。原因很简单PDB结构或对接结构来源于实验或算法原子的坐标本身就有不确定性而你在加水时水分子是根据网格随机放置的很可能有一个水分子正好落在蛋白侧链上。如果用这样的结构直接开始MD最大问题是力场里原子间排斥力随距离急剧上升两原子距离比平衡键长小0.1 nm时排斥能量可能高达几千甚至几万kJ/mol这个巨大的力会让原子在一步之内获得极高速度体系温度瞬间爆表模拟直接崩溃。能量最小化的目的就是在模拟开始前用最速下降等算法找到势能面的局部极小点把不合理的接触“推开”让最大受力降到可接受范围。5.2 em.mdp参数与收敛判断能量最小化mdp参数一般长这样integrator steep emtol 1000 emstep 0.01 nsteps 50000 cutoff-scheme Verlet nstlist 10 rlist 1.2 rcoulomb 1.2 rvdw 1.2 coulombtype PMEemtol 1000的单位是kJ/mol/nm意思是当体系中最大受力小于1000时最小化收敛。这个阈值对应能量梯度的判断等价于一个比较“宽松”但足够安全的收敛标准。emstep是最陡下降的初始步长取了0.01 nm如果算法发现能量上升会自动减半步长。执行gmx grompp -f em.mdp -c complex_ions.gro -r complex_ions.gro -p topol.top -o em.tpr gmx mdrun -deffnm em -v-r用于生成位置约束文件虽然最小化这里用不到但后面平衡需要posre文件所以建议加上。如果用了-r complex_ions.groposre.itp会在grompp阶段自动读进来。最小化完成后看em.log末尾的势能正常情况下势能应该是负的且收敛输出提示达到Fmax收敛。也可以用gmx energy -f em.edr -o potential.xvg选择Potential项画图。一个合理的能量最小化结果是从一个较高的正值或较小的负值开始快速下降到稳定的负值最后几乎不再变化。如果在能量轨迹里看到一千多甚至上万的势能且剧烈波动说明需要检查结构。5.3 最小化不收敛时的两条出路能量最小化不收敛通常是两种原因一是结构里有严重的原子重叠二是力场参数出了问题尤其是配体拓扑。第一种情况可以在最小化之前先用PyMOL或者在gro文件里检查配体和周围残基的距离。如果发现有原子距离小于0.1 nm先手动调整配体位置或删除个别重叠水分子。注意pdb2gmx和solvate添加的水分子编号可能在gro文件里难找你也可以直接用gmx select选距离蛋白太近的水分子输出编号手动删掉。第二种情况配体itp里有不合理的键长、键角或原子类型最小化时会有巨大内应力。可以尝试先把配体从体系中拿出来单独对配体做气相最小化看它能否正常收敛。如果配体单独也发不了最小化问题出在配体拓扑需要回到参数生成环节。还有一种临时方案是给配体加位置约束比如在em.mdp里define -DFLEXIBLE配合配体posre文件先把蛋白和水调整好再释放配体做第二轮最小化。这个方法虽然有效但也会掩盖一部分配体拓扑问题所以只建议作为排查手段。6. NVT/NPT平衡把体系带到300K和1bar6.1 NVT升温控温器选择与位置约束能量最小化后体系处于势能局部极小的低温状态需要升温到目标温度同时让蛋白骨架不能乱动给水分子和侧链充分松弛的时间。这个过程就是NVT平衡即保持原子数、体积、温度不变。NVT平衡的mdp核心参数integrator md dt 0.002 nsteps 25000 tcoupl V-rescale tc-grps System tau_t 0.1 ref_t 300 pcoupl no constraints h-bonds constraint_algorithm lincs continuation nonsteps 25000乘以dt 0.002就是50 ps。作为升温阶段这个长度足够让体系从初始温度到达300K并稳定下来。NVT平衡里最重要的是位置约束蛋白重原子被固定在初始位置只允许水、离子和蛋白侧链小范围移动。这么做是因为直接从最小化后的结构开始放开全部自由度蛋白可能会在温度上升的瞬间产生明显构象漂移尤其在溶剂还没完全平衡时这种漂移很容易把蛋白推向不合理的结构。位置约束通过define -DPOSRES启用。在topol.top里protein和配体通常会各有一个#ifdef POSRES的posre引入段。-DPOSRES会让蛋白重原子受到一个力常数为1000 kJ/mol/nm²的谐振子约束。NVT阶段加上这个定义跑50 ps后体系温度通常会稳定在300K附近动能和势能也会趋于平缓。6.2 NPT加压为什么要看密度是否收尾NVT平衡完成后体系温度已经正确但压力还不正确——因为在固定体积下溶剂密度未必对应1 bar压力下的密度。NPT平衡在保持温度的同时通过压力耦合调节盒子体积让体系密度和压力达到目标值。对于等温等压条件这更接近一个真实的生理环境。NPT平衡mdp参数integrator md dt 0.002 nsteps 500000 tcoupl V-rescale tc-grps System tau_t 1.0 ref_t 300 pcoupl Parrinello-Rahman pcoupltype isotropic tau_p 2.0 ref_p 1.0 compressibility 4.5e-5 constraints h-bonds continuation yesNPT平衡时长建议至少500 ps有的体系甚至需要1 ns才能让密度稳定。这里为什么tau_p 2.0而不是0.1压力耦合的响应时间比温度慢得多时间常数太短会导致盒子体积剧烈震荡密度曲线噪声很大。NPT平衡结束后检查gmx energy里的Density项水的理论密度约1000 kg/m³蛋白和配体的存在会让体系密度略高通常在1000到1100 kg/m³之间。如果密度在这个范围内且曲线在平衡后期基本水平说明体系已经达到正确的压力条件。如果密度还在缓慢上升或下降说明平衡不够要延长NPT时间。6.3 平衡成功的三个检查每次平衡结束我都会例行检查三样东西缺一不可温度曲线是否稳定在300K附近。如果温度持续偏离目标好几K说明控温参数或初始速度有问题先别急着跑生产。压力曲线是否在平均值附近波动。Parrinello-Rahman的输出压力噪声较大是正常的关键是看平均值是否在1 bar附近。如果压力一直在几十bar甚至上百bar波动盒子尺寸可能不合适要检查周期距离是否足够。密度是否收敛到一个稳定值。密度的收敛意味着溶剂已经被“压实”到正确的密度这是整个NPT阶段成功的核心指标。另外可以看一下势能曲线平衡后期势能不应该有明显上升趋势。7. 成品模拟参数文件、GPU加速与续跑7.1 production mdp参数逐行平衡结束后可以进行成品模拟。生产模拟的mdp和NPT平衡有很多相似之处但有一个关键点继续使用Parrinello-Rahman压力耦合时如果你是从NPT平衡的cpt文件续跑可以无缝衔接如果生产模拟单独从头开始则需要先跑一小段Berendsen压力耦合的预平衡。最稳妥的方式是直接续跑NPT平衡的检查点文件。一个常用的生产mdpintegrator md dt 0.002 nsteps 50000000 nstxout-compressed 5000 cutoff-scheme Verlet rlist 1.2 rcoulomb 1.2 rvdw 1.2 coulombtype PME fourierspacing 0.12 tcoupl V-rescale tc-grps System tau_t 1.0 ref_t 300 pcoupl Parrinello-Rahman pcoupltype isotropic tau_p 2.0 ref_p 1.0 compressibility 4.5e-5 constraints h-bonds continuation yes关键行解释nsteps 50000000对应100 ns模拟。nstxout-compressed 5000即每5000步也就是10 ps存一帧压缩轨迹。100 ns会产生1万帧文件体积和分辨率足够做绝大多数分析。如果只想看大趋势可以设为10000甚至20000但做氢键和距离分析时帧数太少会影响统计。continuation yes告诉Gromacs这是一段续跑不要重新生成初始速度。PME的fourierspacing 0.12表示倒空间格点间距约0.12 nm精度足够没必要更小否则会显著增加计算量。运行命令gmx grompp -f md.mdp -c npt.gro -t npt.cpt -r npt.gro -p topol.top -o md.tpr gmx mdrun -deffnm md -v7.2 GPU加速和速度评估Gromacs的性能优化做得非常出色支持GPU加速后可溶蛋白体系的模拟速度比纯CPU快一个数量级以上。在mdrun时指定gmx mdrun -deffnm md -v -nb gpu -pme gpu -bonded gpu -update gpu-nb是非键计算-pme是长程静电的PME计算-bonded是键合相互作用-update是坐标更新。如果你的GPU支持一般NVIDIA的Pascal架构以后都支持把四个都设为gpu可以让非键和键合计算都跑到GPU上。Gromacs的日志和终端会输出“GPU”相关性能数据例如Performance: 185 ns/day之类的指标。这里提醒一个容易忽略的问题-update gpu会要求GPU计算和CPU通信有较高带宽某些老显卡或者集群节点配置不当时反而更慢。可以先试跑几千步比较不同设置的性能再决定要不要用-update gpu。另外盒子里如果有大量水分子非键计算是主要瓶颈GPU提升非常明显。但如果体系很小比如只有几千个原子GPU的空闲等待可能吃掉提升这时纯CPU反而更稳。具体看性能测试结果。7.3 模拟中断后的续跑跑100 ns甚至更长的模拟中途集群重启、作业超时、硬件故障几乎是必然事件。不要慌Gromacs的检查点机制设计得很好。模拟中断后目录里有md.cpt文件续跑命令gmx mdrun -deffnm md -cpi md.cpt -s md.tpr -vGromacs会从检查点状态继续而不是从头开始。注意续跑时不要让两个mdrun同时操作同一个tpr和cpt文件否则会写坏检查点。集群环境下尤其要注意作业管理系统是否允许你在上次作业结束后再次提交。还有一个实用技巧运行期间定期把md.cpt备份到另一个目录因为如果磁盘满了cpt写不进去模拟会中止备份可以避免丢失太长的轨迹。模拟稳定跑完的标志是目录下出现md.gro、md.xtc和完整的md.log日志末尾会显示“Finished mdrun on node”之类的信息。8. 轨迹分析从PBC处理到RMSD/RMSF/氢键/自由能8.1 第一件事永远是去周期性与居中模拟是在周期性盒子中进行的分子在模拟中可能穿越盒子边界直接分析原始md.xtc会看到蛋白链“断裂”成两半配体位置在盒子两端跳跃。所有轨迹分析的第一步都是把周期性效应处理干净。推荐的处理方式gmx trjconv -s md.tpr -f md.xtc -o md_whole.xtc -pbc whole gmx trjconv -s md.tpr -f md_whole.xtc -o md_center.xtc -pbc mol -center第一条命令-pbc whole把每个分子写完整不跨越边界。第二条命令-pbc mol -center先把体系居中再把分子置于盒子中心。执行第二条时Gromacs会询问选择哪个组做居中一般选Protein因为以蛋白为中心来分析配体周围环境最合理。处理完成后所有后续分析都基于md_center.xtc不要再回头用原始轨迹。经验不足的人最容易在这里翻车RMSD曲线看起来大幅度跳变其实只是分子在盒子边界来回穿行造成的假象。8.2 RMSD和RMSF稳定性与残基柔性RMSD是最基础也是最重要的稳定性指标。它计算每个帧的蛋白骨架原子相对于参考结构的均方根位移。参考结构的选择有讲究如果想知道模拟过程中蛋白是否发生了大幅构象变化用能量最小化后的结构或PDB晶体结构作参考如果想知道模拟是否达到稳定平衡也可以把平衡之后的轨迹平均结构作为参考。常用命令gmx rms -s md.tpr -f md_center.xtc -o rmsd.xvg -tu ns交互选择时拟合组选Backbone计算组也选Backbone。-tu ns让x轴时间以ns为单位方便作图。RMSD曲线通常在初始阶段快速上升然后在某个平台值附近波动。平台值因蛋白大小而异小蛋白可能在0.1到0.2 nm大蛋白或柔性结构在0.3到0.5 nm都很常见。如果RMSD在模拟后期仍然持续上升且没有平台趋势说明体系没有真正稳定可能需要重新评估模拟条件和时长。RMSF反映每个残基在模拟中的波动幅度gmx rmsf -s md.tpr -f md_center.xtc -o rmsf.xvg -res-res会把多个原子的RMSF按残基平均得到一个“每残基一个值”的序列方便绘制曲线。看RMSF时重点关注结合位点残基如果配体结合口袋附近残基的RMSF明显高于蛋白整体平均水平可能说明配体结合得不够稳定结合pose在模拟中有明显晃动如果口袋残基RMSF较低说明结合模式比较刚性是一个稳定结合的信号。8.3 配体结合相关的分析氢键、距离、接触RMSD和RMSF不够直接回答“配体结合得稳不稳”还需要看分子层面的相互作用。氢键分析gmx hbond -s md.tpr -f md_center.xtc -num hbond_num.xvg -life hbond_life.xvg交互选择两个组一个是Protein另一个是配体组LIG。-num输出的是每个时间点的氢键数目随时间的变化曲线-life输出氢键存在时间的分布。平均氢键数量越多、寿命越长说明配体和蛋白之间的极性相互作用越稳定。注意氢键的定义默认是供体-受体距离小于0.35 nm且角度偏离直线小于30度这个标准可以在mdp参数里调整但对新手先用默认即可。关键距离监控gmx distance -s md.tpr -f md_center.xtc -select atomname ...-select的语法需要指定具体的原子。比如配体的某个氧原子和蛋白某个残基侧链氮原子之间的距离。这类“看特定距离是否稳定地维持在一定范围”的分析比看一大堆氢键更有针对性尤其是已知关键残基时。选定原子编号可以从gmx select交互式命令行查询。接触面积分析可以用gmx mindist或第三方工具如PyMOL的contact命令统计配体原子在0.3到0.4 nm范围内接触到的蛋白残基。这个指标能直观说明配体周围有没有形成稳定的疏水口袋。8.4 结合自由能MM/PBSA的跨界提醒如果想更进一步估算结合自由能主流的做法是用MM/PBSA方法。Gromacs本身不带现成的MM/PBSA实现现在常用第三方工具gmx_MMPBSA。基本流程用复合物tpr文件和轨迹文件分别计算复合物、受体和配体在气相中的分子力学能量然后用隐式溶剂模型估算溶剂化自由能二者相减得到结合自由能。命令大致是gmx_MMPBSA -O -i mmpbsa.in -cs complex.tpr -ct complex.xtc -cg 1 13 -cp complex.top ...-cg指定受体和配体分别对应哪两个index组-cp提供复合物的拓扑文件。实际操作中坑非常多轨迹里必须包含蛋白和配体的所有原子而topol.top需要同时包括蛋白和配体参数受体和配体的坐标要从原始轨迹中提取并分别生成tpr。建议初学者先把前面的基础分析做好再考虑MM/PBSA因为算出来的绝对值对参数和方法选择非常敏感。一个重要提醒MM/PBSA给出的结合自由能是“能量估计”不是实验结合能。它的价值在于比较一系列相似配体或不同突变体之间的相对趋势或者作为结合模式排序的参考不要过度解读绝对值。误差范围内你只能看出“强结合”“弱结合”这样的定性结论不要拿它和小数点后两位的实验IC50去对表。9. 跑完这个流程之后我踩过的坑和一些习惯第一批模拟能顺利跑完结果分析也画出了图很多人就觉得大功告成了。但作为已经跑过数不清个蛋白-配体体系的人我始终有几件事每次都会做一遍并且建议你也养成习惯。第一跑完立即回头看日志里的警告。md.log里会有大量注释有些是无关紧要的但如果出现“turned around, warning, LINCS warning”或者“Step ... too long”的字样意味着模拟中有局部原子受力过大或者约束算法出问题。哪怕模拟最终正常结束这些警告也会影响轨迹质量一定要回到对应时间段检查。第二把轨迹按时间段拆开看RMSD。不要只看整条轨迹的平均值要分早期、中期、后期分别看。一个常见的假象是整段RMSD曲线平均看起来挺平但前20 ns一直在爬升后20 ns趋于稳定。这时候如果你把前60 ns也算进“平衡后”区间配体结合分析会包含大量构象搜索阶段的数据结果会偏乐观或偏模糊。我的习惯是去掉前20%作为“预平衡”丢弃只分析后面的稳定段。第三选原子组时谨慎使用System组。很多分析默认可以选择System但你真正关心的可能是蛋白骨架或者是配体周围一定范围内的残基。用gmx make_ndx建好常用组比如Protein_Backbone、LIG、Protein_LIG、within 0.5 of LIG这类索引组能让后续分析高效很多也能避免选错组导致的数据无效。第四轨迹可视化一定要做。再多的数值分析都无法代替亲眼看一下配体在结合口袋里到底怎么动。用VMD或PyMOL加载md_center.xtc和md.tpr慢速播放最后20 ns的轨迹。我曾经遇到过RMSD、氢键、接触面积全部正常的体系最后盯着轨迹看才发现配体在模拟后期旋转了180度原来结合模式早已改变。数值指标只是帮你定位问题视觉确认才能让你真正理解体系行为。第五保存好所有输入文件和参数版本。Gromacs版本更新很快同一套mdp在不同大版本下的默认值可能有差异。我在每个项目目录下都会保留一份version.txt记录Gromacs版本、mdp文件、力场版本和PDB编号。这样几个月后回来看结果还能复现当初的完整流程。新手尤其容易忽略这一点等你论文返修要求补跑模拟时发现自己忘了当初用的参数才是最痛苦的。最后说一句实操体会蛋白-配体动力学模拟这个流程本质上是一项需要反复调试和验证的技术活第一次跑通花了快两周第二次可能只要半天第五次的时候我已经能在一个小时内完成从原始结构到成品轨迹的全部准备。关键在于每跑完一个体系都把当时的命令、参数和踩坑记录整理成自己的checklist。等你的checklist越积越多就会发现Gromacs并没有想象中那么神秘——它只是要求你在每个环节都足够细心和诚实。