简介：《分子动力学步骤说明.pdf》面向刚接触分子动力学（MD）模拟的初学者，以及需要撰写论文方法部分的科研人员，以线性叙述梳理从体系构建到数据分析的完整流程，帮助读者先建立起可操作的模拟整体框架。内容涵盖加溶剂与加离子、能量最小化、NVT与NPT预平衡、生产项模拟及轨迹分析，并说明orca与GROMACS在流程中的定位；对RMSD、回旋半径Rg、SASA等常规分析指标，也给出使用提醒与常见误区。压缩包共1个PDF文件，约6.87MB，篇幅紧凑，适合作为上机前的流程导读，也可在写作方法部分时作为对照清单。已有444人学习，读者可据此理解各阶段的目的与判据，避免跳过预平衡、误读RMSD趋于平稳即收敛等问题，并理清论文书写顺序与实操流程之间的对应关系。
1. 拿到一份分子动力学步骤说明，先别急着敲命令
很多人第一次接触分子动力学模拟，是从一份“步骤说明”开始的——里面写着体系构建、能量最小化、NVT、NPT、生产模拟、数据分析六段，每段都像天书。手头这份 PDF 的价值不在于教你某个软件按钮怎么点，而在于把“论文 METHOD 怎么写”和“命令行里实际怎么跑”这两条线并到了一起。它适合两类人：一是刚进组、被导师扔了一个蛋白-配体复合物、需要在几周内跑出可发表轨迹的新手；二是跑过几轮但总在预平衡崩、RMSD 图不会读、MMPBSA 能量项符号看不懂的老手。下面按这份说明的六个阶段，把每一步的命令、参数和翻车点拆开讲。
2. 体系构建与能量最小化：从 pdb2gmx 到 emtol 的四个关键参数
2.1 力场、水模型与质子化状态：先定这三样再动手
分子动力学模拟的能量和力都来自力场，力场选错后面全白搭。蛋白用 amber14sb 是这几年的稳妥选择，小分子配体走 GAFF2，水模型配 TIP3P，三者是配套的。原文给的组合就是 amber14sb + GAFF2 + TIP3P，这个组合在蛋白-配体复合物里兼容性好，文献里也常见。换力场不是不行，但换之前要想清楚：MMPBSA 的熵项和溶剂化项参数是按力场给的，中途混搭容易在能量分解那一步对不上号。
质子化状态是新手最容易忽略的一步。pdb2gmx 会按力场默认的 pKa 猜，但 HIS 残基经常猜错——它可能被质子化成 HID、HIE 或 HIP 三种形式，不同形式影响氢键网络和口袋电荷。我的习惯是先看蛋白口袋附近有没有 HIS，有的话用 PROPKA 或 H++ 先跑一遍 pKa 预测，再手动在 pdb2gmx 交互界面里指定。
小分子这边，流程比蛋白多两步：先用量化软件优化几何构型、算 RESP 电荷，再用 antechamber 生成 GAFF2 参数。原文用的是 Orca 在 B3LYP/6-31G* 下优化，RESP2.0 算电荷。这套流程在 gmx_MMPBSA 计算里对电荷一致性要求高，别偷懒直接用 BCC 电荷凑合——BCC 快，但和 RESP 的静电势拟合精度差一截，结合自由能的绝对值会有偏移。
# 蛋白拓扑生成：-ff 指定力场，-water 指定水模型，-ignh 让工具按力场重新加氢 gmx pdb2gmx -f protein.pdb -o protein_processed.gro -water tip3p -ff amber14sb -ignh # 小分子配体参数：先假定已完成 Orca 优化并拿到 RESP 电荷文件 # -c resp 表示读取 RESP 电荷，-fi/-fo 是输入输出格式 antechamber -i lig_resp.mol2 -fi mol2 -o lig_gaff.mol2 -fo mol2 -c resp -s 2 parmchk2 -i lig_gaff.mol2 -f mol2 -o lig.frcmod参数说明：-ignh会丢掉 PDB 里的氢，按力场模板重新加，避免实验结构里的氢位置和力场冲突；-ff换成 amber99sb-ildn 也能跑，但小分子参数库匹配度要自己验；antechamber的-s 2是输出中间文件方便排查电荷是否合理。代码跑完先看protein_processed.gro末尾的盒子尺寸和topol.top里的分子数，确认蛋白链数没被拆错。
2.2 加水、加离子与中和电荷：盒子尺寸和盐浓度怎么定
水盒子的大小直接决定周期边界条件会不会引入自相互作用。经验值是溶质到盒子边缘至少留 1.0 nm，做蛋白-配体复合物建议 1.2 nm 起步。盒子类型选十二面体比立方体省水，同样距离下体积小一截，算起来快。加离子分两步：先中和体系净电荷，再按生理条件补到 0.15 M。原文写的是“添加钠离子平衡体系电荷”，实际操作里 genion 会根据体系电荷自动选 NA 或 CL。
有一个容易被忽略的点：genion 是随机替换水分子放离子的，它不管附近有没有原子冲突。所以加完离子不能直接跑动力学，必须过一遍能量最小化。这也是原文把能量最小化单独列一步的原因之一。
# 定义盒子：-c 居中，-d 溶质到盒边距离，-bt 盒子类型 gmx editconf -f protein_processed.gro -o box.gro -c -d 1.2 -bt dodecahedron # 填充水分子：-cs 指定水模型坐标文件 gmx solvate -cp box.gro -cs spc216.gro -o solv.gro -p topol.top # 先做一次极简 grompp 生成 tpr，genion 需要读 tpr gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr # 加离子：-neutral 只中和，-conc 0.15 补盐到生理浓度 gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral -conc 0.15参数说明：-d 1.2是溶质到盒边的距离，单位 nm，不是盒子总边长；-bt dodecahedron十二面体盒子在可视化时看起来是斜的，不影响计算；-conc 0.15是额外加的盐浓度，-neutral会先中和再加盐；genion 运行时会问替换哪一组溶剂，选SOL。跑完检查topol.top里[ molecules ]段，SOL 数量应该减少，NA/CL 出现对应行。
2.3 能量最小化：为什么不能直接开跑，emtol 怎么看
前面三步留下的坑都在这儿集中爆发：pdb2gmx 按模板加的氢可能和重原子距离过近，solvate 填的水分子取向未必是最佳氢键结构，genion 随机放的离子可能和蛋白表面原子重叠。这些重叠在动力学里就是巨大的力，一步积分就能把体系炸飞，日志里常见Atom X has NaN coordinates或者Segmentation fault。
能量最小化的目的不是找到全局能量最低点——那在有限步数里做不到，也没必要——而是把体系拉到“没有原子重叠、几何构型合理”的起始点。最陡下降法适合前期消除大冲突，共轭梯度法适合接近极小值后精细收敛。原文用的是最陡下降法 50000 步，emtol 1000 kJ/mol/nm，这个阈值在实际项目里够用，跑完看em.log最后一行 Fmax 是否低于 emtol 即可。
# em.mdp 关键参数（写进文件，不是命令行） # integrator = steep ; 最陡下降 # nsteps = 50000 ; 最大步数 # emtol = 1000.0 ; 收敛阈值 kJ/mol/nm # emstep = 0.01 ; 初始步长 # coulombtype = PME ; 静电用 PME # rcoulomb = 1.4 ; 库伦截断 nm # rvdw = 1.4 ; 范德华截断 nm gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr gmx mdrun -v -deffnm em参数说明：emtol不是越小越好，1000 已经能保证没有明显重叠，设到 100 以下会多跑很多步但收益有限；emstep是最陡下降的初始步长，体系冲突严重时可以降到 0.005；rcoulomb和rvdw统一设 1.4 nm 是 amber 力场的常见搭配，和后续 NVT/NPT 保持一致，别在不同阶段改截断距离。
3. NVT 与 NPT 预平衡：为什么不能跳过恒温直接上恒压
3.1 两步预平衡的物理理由：先稳温度，再稳密度
能量最小化只保证了“没有原子重叠”，但溶剂分子和离子的取向、溶质周围的氢键网络还没弛豫。直接开无限制动力学，体系会在最初几百皮秒里剧烈调整，轻则轨迹前段不可用，重则盒子变形、蛋白解折叠。原文说得很直白：不经过预平衡就直接跑非限制动力学，体系可能崩溃。
NVT 阶段的任务是把体系拉到目标温度，让溶质在溶剂里找到合理取向。这一步保留位置限制（-DPOSRES），蛋白重原子被约束在初始位置附近，只有水和离子自由运动。等温度稳定了，再进 NPT 阶段稳住压力。原文特别提醒：在恒压条件下同时对可能远离平衡的体系施加速度产生和坐标缩放，会导致不稳定。这就是为什么不能一步到位直接上 NPT——恒压器在密度还差很远的时候会疯狂缩放盒子，和温控器打架。
NPT 阶段用 Parrinello-Rahman 恒压器比 Berendsen 更严格，但 Berendsen 在预平衡阶段更稳。我的习惯是 NVT 和 NPT 预平衡都用 Berendsen，生产阶段再切 Parrinello-Rahman。这份说明里没写具体恒压器型号，但原理上“预平衡用温和的、生产用严格的”是通用做法。
3.2 mdp 文件里真正要改的五行
预平衡的 mdp 不用从头写，拿官方模板改几行就行。真正影响结果的是温控器、恒压器、时间常数和位置限制开关。tc-grps把蛋白和非蛋白分成两组分别控温，是因为两者的热容和弛豫时间不一样，混在一起控温会让蛋白温度被水拖偏。
# nvt.mdp 关键段 # define = -DPOSRES ; 保留位置限制 # tcoupl = V-rescale ; 温控器 # tc-grps = Protein Non-Protein # tau-t = 0.1 0.1 ; 温度耦合时间常数 ps # ref-t = 300 300 ; 目标温度 K # pcoupl = no ; NVT 阶段不控压 # npt.mdp 关键段 # define = -DPOSRES ; 仍然保留位置限制 # pcoupl = Berendsen ; 预平衡阶段用温和恒压器 # pcoupltype = isotropic # tau-p = 1.0 ; 压力耦合时间常数 ps # ref-p = 1.0 ; 目标压力 bar # compressibility = 4.5e-5 ; 水的等温压缩系数参数说明：tau-t 0.1是 V-rescale 的常见值，太大温度响应慢，太小会震荡；tau-p 1.0在预平衡阶段比生产阶段的 2.0 更激进，目的是快速把密度拉到位；compressibility 4.5e-5是水的值，换其他溶剂要查手册。两个阶段都要保留-DPOSRES，生产阶段才去掉。跑之前确认posre.itp在拓扑里被#include了，否则位置限制不生效，预平衡等于白做。
3.3 预平衡成功的判断：温度、密度、压力三条线
判断预平衡是否成功，不看步数跑完没有，看三条曲线是否走平。温度应该在目标值上下 5 K 内波动，密度应该接近 1000 kg/m³（水体系），压力在 1 bar 附近有波动但均值不漂。提取这些量用gmx energy，它能从 edr 文件里读任意能量项。
有一个常见误判：只看最后一步的值。正确做法是看后 20% 的时间段均值，前段还在爬升不算数。如果密度一直往下掉，说明恒压器在压缩盒子但体系没稳住，回头检查能量最小化是不是没收敛，或者初始盒子给得太小导致水分子被挤得太密。
# 提取 NPT 阶段的温度、密度、压力，用于判断是否平衡 printf "Temperature\nDensity\nPressure\n" | gmx energy -f npt.edr -o npt_check.xvg # 快速看均值：-b 指定从哪个时间点开始统计，避开前段爬升 printf "Density\n" | gmx energy -f npt.edr -o density.xvg -b 200参数说明：-b 200表示从 200 ps 开始统计，具体数值按你的预平衡长度定，一般取总时长的一半之后；npt_check.xvg可以直接用文本编辑器打开看两列数据，也可以丢进绘图工具。压力波动大是正常的，瞬时压力在几百 bar 上下跳，看的是均值。如果温度或密度在最后 20% 还在单调漂移，延长预平衡或者回头查参数，别硬着头皮进生产阶段。
4. 生产模拟与轨迹分析：RMSD、RMSF、Rg、SASA 的读图方法
4.1 去 PBC：trjconv 的三步走与“破碎分子”的由来
生产模拟跑完，轨迹里的分子经常看起来“碎了”——蛋白一半在盒子这边、一半在那边，或者整个蛋白跑到盒子外面。这不是模拟出错，是周期边界条件的可视化表现。mdrun 为了计算效率，让分子跨边界时保持“破碎”状态输出，后处理再拼回来。原文说得很清楚：与其在运行中花时间重映像，不如让 mdrun 输出破碎结构，后面再连起来。
去 PBC 的标准流程是三步：先消除整体平动和转动，再把蛋白放回盒子中心，最后按需要输出特定组。第一步的-pbc mol是按分子整体处理，不是按原子；-center把选定组居中；-ur compact用紧凑表示。第二步的-fit rot+trans做叠合，消除转动和平动，做 RMSD 之前必须做这一步，否则 RMSD 里混进了整体漂移。
# 第一步：按分子整体去 PBC，并居中 gmx trjconv -s md.tpr -f md.xtc -o md_center.xtc -center -pbc mol -ur compact # 第二步：叠合到参考结构，消除转动和平动 gmx trjconv -s md.tpr -f md_center.xtc -o md_fit.xtc -fit rot+trans # 第三步：只输出蛋白组，减小文件体积（可选） gmx trjconv -s md.tpr -f md_fit.xtc -o protein_only.xtc -n index.ndx参数说明：第一步会问选哪个组做居中，通常选Protein或System；第二步会问叠合用哪个组，选Backbone更稳；-ur compact只影响输出表示，不影响物理。跑完用可视化工具打开md_fit.xtc确认蛋白是完整的，如果还有碎片，检查-pbc mol有没有加对。
4.2 RMSD 与 RMSF：一个看整体稳定性，一个看局部柔性
RMSD 是把每一帧和参考结构做最小二叠合后算的均方根偏差，反映整体结构随时间偏离多少。RMSF 是每个残基相对于平均位置的波动，反映局部柔性。两者互补：RMSD 平台不代表体系“收敛”，原文专门提醒过这一点——结构可以在 RMSD 平台期继续发生大幅构象变化，只是这些变化在叠合后相互抵消了。所以 RMSD 要结合 Rg、SASA 一起看，不能单指标下结论。
RMSF 的读法要结合结构。原文给的阈值是 0.3 nm，超过这个值的残基算高波动。高波动区域可能是 N/C 端、loop 区，也可能是口袋附近的柔性残基。如果高波动正好落在结合口袋周围，要小心——可能是配体没结合稳，也可能是口袋本来就柔性大。我的习惯是把 RMSF 高值残基标在结构图上，和口袋残基列表对照，看是“功能相关柔性”还是“模拟没稳住”。
# RMSD：先选 Backbone 做叠合，再选 Backbone 算偏差 gmx rms -s md.tpr -f md_fit.xtc -o rmsd.xvg -tu ns # RMSF：-res 按残基平均输出 gmx rmsf -s md.tpr -f md_fit.xtc -o rmsf.xvg -res参数说明：-tu ns把时间轴单位换成 ns，图更好读；gmx rms运行时会问两次组选择，第一次选叠合组（Backbone），第二次选计算组（Backbone 或 Protein）；gmx rmsf -res输出的是每个残基一个值，不加-res是每个原子一个值。RMSF 图里横坐标是残基号，纵坐标是 nm，0.3 nm 的参考线可以自己画上去。
4.3 Rg、SASA、氢键与 DSSP：四个指标联动读构象变化
Rg 衡量蛋白的密实程度，折叠良好的蛋白 Rg 应该基本不变，去折叠时 Rg 会增大。SASA 是溶剂可及表面积，SASA 下降通常意味着蛋白紧缩或者口袋被配体盖住。氢键数在蛋白-配体界面稳定时应该平稳，大幅下降可能提示结合模式变了。DSSP 把二级结构分成螺旋、折叠、卷曲等类别，看各类比例随时间的变化。
这四个指标要联动读。原文给了一个很好的例子：Rg 在 7 ns 和 24 ns 各有一个波动，正好对应 RMSD 的两次波折，说明这两个时间点体系在经历构象转变。SASA 在 0-50 ns 波动下降，提示蛋白逐步紧缩。DSSP 显示 coil 减少、β-sheet 和 α-helix 增加，和 SASA 下降互相印证。单看任何一个指标都容易过度解读，四个放一起才能说“体系在往哪个方向变”。
# 回旋半径 gmx gyrate -s md.tpr -f md_fit.xtc -o rg.xvg # 溶剂可及表面积，-surface 指定计算组 gmx sasa -s md.tpr -f md_fit.xtc -o sasa.xvg -surface 'Protein' # 氢键数量 gmx hbond -s md.tpr -f md_fit.xtc -num hbnum.xvg # DSSP 二级结构，需要单独安装 dssp 可执行文件 gmx do_dssp -s md.tpr -f md_fit.xtc -o dssp.dat -sc dssp_sc.xvg参数说明：gmx gyrate默认对 Protein 组计算，输出的是整体 Rg；gmx sasa -surface 'Protein'指定计算蛋白表面，不指定会算整个体系；gmx hbond的-num输出每帧氢键数，默认按蛋白-配体或蛋白-蛋白内部算，具体组选择在交互里定；gmx do_dssp需要系统里有 dssp 程序，版本不匹配会报错，常见做法是用 conda 单独装一个。
5. 避坑与排查：五条从崩溃日志里抠出来的经验
5.1 能量最小化第一步就报 “Atom X has NaN coordinates”
现象：grompp 能过，mdrun 跑几步就崩，日志里出现 NaN 坐标或者 “There are 1 atoms with infinite coordinates”。
原因：几乎都是初始结构里有原子重叠。pdb2gmx 加的氢和重原子距离过近、genion 随机放的离子和蛋白表面原子重叠、或者输入 PDB 本身有 altloc 没处理干净。
解决：先跑gmx editconf -f xxx.gro -o check.pdb用可视化工具看有没有明显重叠；把 emstep 降到 0.005 再跑一次；如果还崩，用gmx trjconv -pbc mol检查盒子有没有问题。实在找不到，把 genion 那一步重做，换一个随机种子（-seed）再试。血泪经验：NaN 出现之前通常有一步 Fmax 异常大，看 em.log 里前几步的 Fmax，能提前定位是哪个区域出的问题。
5.2 NVT 跑着跑着温度飙到几百 K
现象：NVT 预平衡前几十 ps 温度正常，之后突然往上飙，最终稳定在 400 K 以上下不来。
原因：最常见的是位置限制没生效。posre.itp没被 include 进拓扑，或者define = -DPOSRES拼写错了，蛋白重原子自由运动，和溶剂碰撞产生大量热。其次是时间步长太大，2 fs 配氢键约束（constraints = h-bonds）才安全，没加约束还用 2 fs 会积累误差。
解决：检查topol.top里有没有#include "posre.itp"，检查 mdp 里define那行有没有拼错；确认constraints = h-bonds和dt = 0.002配套。如果温度还是高，把tau-t从 0.1 调到 0.5 让耦合更温和，先把温度拉下来再逐步收紧。
5.3 NPT 密度一路往下掉，盒子越缩越小
现象：NPT 预平衡开始后密度从 1000 一路掉到 800 以下，盒子体积持续缩小，压力均值偏离 1 bar 很远。
原因：恒压器参数过激或者初始盒子太小。tau-p设得太小（比如 0.1）会让恒压器反应过度，盒子剧烈震荡；初始盒子留边不足 1.0 nm，水分子被挤得太密，恒压器一放开就拼命扩。也有可能是能量最小化没收敛，残留的原子重叠在 NPT 里被放大。
解决：先把tau-p调回 1.0 或 2.0，compressibility确认是 4.5e-5；回头检查 editconf 的-d是不是至少 1.0；重跑能量最小化，确认 em.log 里 Fmax 低于 emtol 再进 NPT。如果密度已经掉下去了，别在同一个 tpr 上接着跑，从头做一遍预平衡比硬救省时间。
5.4 RMSD 图前面一段是断的或者跳变
现象：gmx rms出来的 xvg 文件里，前几 ns 的 RMSD 值异常大或者直接是断点，后面才正常。
原因：轨迹没做叠合就去算 RMSD，整体平动和转动混进了偏差里；或者-pbc mol没做，蛋白跨盒子边界导致坐标不连续；还有一种可能是参考结构（-s的 tpr）和轨迹的原子数不一致，比如生产模拟用了不同的拓扑。
解决：按 4.1 的三步走先做-pbc mol和-fit rot+trans，再算 RMSD；确认-s用的 tpr 和生产模拟是同一个；如果轨迹是分段跑的（比如续跑），用gmx trjcat拼起来再做后处理。RMSD 前段跳变不一定是模拟问题，先排除后处理顺序。
5.5 MMPBSA 算出来 TOTAL 是正的
现象：gmx_MMPBSA 跑完，FINAL_RESULTS_MMPBSA.dat里 TOTAL 是正值，按预期蛋白-配体结合应该是负的。
原因：最常见的是轨迹没做去 PBC 和叠合，结合界面的原子坐标是错的；其次是拓扑里配体电荷和 tpr 里不一致，MMPBSA 读到的电荷和实际模拟用的对不上；还有可能是选的帧窗口包含了未平衡的前段，构象还没稳定。
解决：确认输入给 MMPBSA 的轨迹是md_fit.xtc而不是原始md.xtc；检查配体参数文件里的电荷和topol.top里一致；用-startframe跳过前 20% 的帧。如果 TOTAL 还是正，单独看 VDWAALS 和 EEl 是不是负的，如果这两项也是正的，基本可以确定是坐标或电荷出了问题，不是方法本身的问题。
6. 进阶：MMPBSA 结合自由能分解与写作骨架
MMPBSA 的读法有优先级：先看 TOTAL 的符号和量级，再看 VDWAALS 和 EEl 的贡献，最后看溶剂化项。原文给的例子 TOTAL 是 -19.28 kcal/mol，VDWAALS 和 EEl 都是负的，说明疏水作用和静电作用都有利；GSOLV 是正的，其中 EGB 大于 0 不利于结合，ESURF 略有利。这个符号模式是大多数蛋白-配体复合物的常态，如果 EGB 特别大把 TOTAL 拉正了，要回头看配体是不是太极性、口袋是不是太 solvent-exposed。
| 能量项 | 含义 | 有利结合的符号 |
|---|---|---|
| VDWAALS | 范德华能 | 负 |
| EEl | 静电能 | 负 |
| EGB | 极性溶剂化能 | 正（不利于） |
| ESURF | 非极性溶剂化能 | 负（略有利） |
| GGAS | 气相自由能（VDWAALS+EEl） | 负 |
| GSOLV | 溶剂化自由能（EGB+ESURF） | 视体系 |
| TOTAL | 结合自由能（GGAS+GSOLV） | 负 |
残基分解的读法和稳态分析不一样，原文专门提醒过：能量分解用的是全轨迹，稳态分析用的是稳定帧，两者侧重不同。分解里正贡献的残基是“干扰作用”，负贡献的是“有利作用”，但别只看绝对值排序，结合结构看哪些残基在界面上、哪些只是长程静电贡献。
# gmx_MMPBSA 典型输入：单轨迹法，跳过前 20% 帧 gmx_MMPBSA -O -i mmpbsa.in -cs md.tpr -ci index.ndx -cg 1 13 \ -ct md_fit.xtc -cp topol.top -o FINAL_RESULTS_MMPBSA.dat \ -eo FINAL_RESULTS_MMPBSA.csv -nogui参数说明：-cg 1 13是受体组和配体组的索引，具体编号看index.ndx；-ct必须是做过 PBC 和叠合的轨迹；-nogui跳过图形界面直接出结果。mmpbsa.in 里设startframe=1, endframe=500, interval=5控制采样密度，帧数不是越多越好，采太密相邻帧高度相关，统计误差反而估不准。
写作骨架按原文的节奏走：方法段写清力场、水模型、软件、系综、截断、时长；结果段第一段宏观描述 RMSD/RMSF/Rg/SASA/氢键的联动，第二段 DSSP 二级结构变化，第三段结合模式（氢键供体受体具体残基），第四段能量分解。图的黑白 600 dpi、彩色 300 dpi 是投稿常规要求，引用文献按原文给的清单对号入座，力场、水模型、PME、MMPBSA 各有对应出处。
从那以后我每次跑完生产模拟，都强制先走一遍trjconv三步再碰任何分析命令，宁可多花十分钟也不让 PBC 的坑混进 RMSD 图里。希望帮到你。
本文还有配套的精品资源，点击获取
