简介：围绕分子动力学模拟全流程整理的PDF文档，面向化学、物理、生物等方向初入MD模拟的研究生与科研人员，也适合需要把模拟流程规范写入论文方法部分的学习者。内容从体系构建、加溶剂与离子、能量最小化，到NVT/NPT预平衡、生产项模拟和轨迹数据分析逐步展开，并补充RMSD、回旋半径、SASA、二级结构等常用分析指标及GROMACS、ORCA工具说明。压缩包共1个PDF文件，约6.87MB，便于集中阅读和随查随用。已有444人学习，说明其在入门导读与流程梳理方面具有一定参考价值。文档还结合规范化写作案例，帮助读者在理解模拟步骤的同时兼顾论文方法部分的表达，对初次接触专业术语、需要建立线性流程认知的人较为友好。
1. 为什么崩掉的总是能量最小化那一步
如果你问做过分子动力学模拟的人，一次 50 ns 的跑动里最容易翻车的是哪一段，多数答案会指向最开始的那几皮秒——不是生产项，而是能量最小化。这份分子动力学模拟步骤说明把整个流程拆成体系构建、能量最小化、NVT 预平衡、NPT 预平衡、生产项、数据分析六段，顺序不是排版习惯，而是一条逐层排除失稳源的因果链。搞懂这条链，你才明白每个 mdp 参数在替哪一步擦屁股，也才知道 3 ps 就报 LINCS warning 时该回头改哪一段，而不是盲目把 dt 从 0.002 调到 0.001 硬扛。
2. 力场与积分器：一个构型的能量和受力是怎么算出来的
分子动力学的主循环朴素得近乎无聊：给定一组原子坐标 r，算出体系能量 U(r) 和每个原子受到的力 F = -∂U/∂r，再在离散时间步长 δt 上把这些力积分成新的坐标，然后重复。所有复杂度都压在第一件事上——U(r) 到底怎么算。经验力场给出的答案是一个把能量拆成若干可加项的近似式，而积分器负责的只是"拿到力之后怎么走下一步"。这两件事分开看都简单，合在一起才决定了你的体系是跑满 50 ns，还是 3 ps 就炸。
2.1 从原子坐标到能量：力场的分项拆解
一个典型的生物分子力场把总能量分成成键项和非键项两大块。成键项描述拓扑上直接相连的原子之间的相互作用，包括键长伸缩、键角弯折、正常二面角扭转和异常二面角（improper）——前两者用谐振子近似，二面角用周期函数，异常二面角主要用来维持手性和平面性。非键项则覆盖所有不通过一两个键直接相连的原子对，由 Lennard-Jones 12-6 势和库仑势两部分组成：前者管范德华吸引与近距离排斥，后者管静电。
这套拆法的关键含义是：能量对坐标的导数可以逐项解析求出来，所以力的计算和能量的计算是同一次遍历完成的。你在 log 里看到的 Epot 和 Fmax 不是两个独立指标，而是同一个函数值和它梯度范数的最大值。理解这一点，后面看能量最小化的收敛判据就不会犯迷糊。
参数体系必须配套。蛋白走 amber14sb，小分子走 GAFF2，水走 TIP3P，这套组合在生物分子体系里是稳妥的默认。配套不只是"名字好听"，还涉及 1-4 相互作用的缩放系数：amber 系列把 1-4 静电和范德华分别乘以 1/1.2 和 1/2.0，这个值写在拓扑的[ defaults ]里，一旦蛋白和小分子的 itp 用了不同的 fudgeQQ，整条能量曲线就会系统性偏移，跑起来看着没事，分析出来的结合自由能却对不上。
小分子的电荷不属于力场自带内容，必须由使用者提供。原文给的路线是先用 ORCA 在 B3LYP/6-31G* 水平做量子化学几何优化，修正键长、键角、二面角，再用 RESP2.0 拟合固定电荷。这里有两个容易被略过的设计动机：一是 6-31G* 这个基组在 RESP 拟合里是历史惯例，它对气相静电势的估计偏高，恰好补偿了从气相到凝聚相时电荷分布被环境极化掉的那部分；二是 RESP 采用两阶段拟合，第一阶段用弱约束把所有原子拉向一个合理初值，第二阶段只对非等价原子放开约束，等价原子（比如甲基上的三个氢）强制同电荷，避免拟合出物理上不合理的电荷分布。
从量级上感受一下：一个键伸缩项偏离平衡位置 0.01 nm，能量大概是几个到几十个 kJ/mol；一对非键原子的 LJ 项在正常距离下也只有几个 kJ/mol。但原子数动辄十万量级，任何一对原子如果坐标重叠到 0.05 nm 以内，LJ 的 1/r¹² 项就能把总能量推到 10⁷ 以上。这就是为什么"体系崩了"在数值上的表现几乎总是能量爆表和力发散，而不是某个键断了。
2.2 2 fs 时间步长：约束算法替你挡掉了什么
Verlet 类积分器的稳定性有个经验约束：时间步长要小于体系最高频振动周期的十分之一左右。生物体系里频率最高的自由度是含氢键的伸缩振动，周期大约 10 fs，所以不加任何处理时，dt 只能取 1 fs 甚至更小。想让 dt 翻倍到 2 fs，办法是不再让这些键真的振动，而是用约束算法把键长固定住——自由度被删掉，最高频率随之消失，积分器就能迈大步。
; 生产项 mdp 关键片段 integrator = md dt = 0.002 ; 2 fs，前提是下面的约束真正生效 constraints = h-bonds ; 只约束含氢键，重原子键仍按力场自由振动 constraint-algorithm = lincs lincs-iter = 1 ; 迭代次数，1 次通常足够 lincs-order = 4 ; 展开阶数，阶数越高单步越贵但越稳 cutoff-scheme = Verlet nstlist = 10 ; Verlet 方案下由程序自动调整，这里只是初值constraints = h-bonds只冻结含氢的键，重原子之间的键依然按谐振子振动，这样既拿到了 2 fs 的步长，又保留了大部分慢自由度的真实动力学。lincs-order = 4指的是约束求解时对逆质量矩阵做四阶展开，阶数太低在大步长下解不干净，阶数太高纯属浪费——对 2 fs 来说四阶是久经检验的甜点。lincs-iter = 1配合四阶展开，在大多数生物体系里一次迭代就能把约束残差压到容差以内。
沿着这条思路再往前一步，有人会把constraints设成all-bonds并把 dt 提到 0.004，这需要额外的氢质量重分配（把氢的质量部分转移给相连的重原子），否则积分误差会累积到不可接受。这不是不能做，但它改变了体系的惯性分布，对需要精确振动谱或者动力学细节敏感的体系并不合适。
LINCS warning这条报错值得单独说一句：它几乎从来不是约束算法本身的问题，而是坐标里存在真实冲突——两个原子靠得太近，约束求解在给定容差内根本解不出来。把 dt 调小往往能让它"跑过去"，但那只是把冲突藏起来，后面的轨迹依旧是脏的。
2.3 PME 与 1.4 nm：长程静电为什么不能简单截断
库仑相互作用按 1/r 衰减，衰减得极慢。如果你只算截断半径以内的静电，等于人为规定"超过这个距离的电荷彼此看不见"，其后果是体系内部出现虚假的电荷分层、界面处的偶极取向失真，以及和截断半径强相关的能量漂移。对带净电荷的蛋白-配体复合物，这种误差足以让结合自由能的数值失去意义。
粒子网格 Ewald（PME）的处理方式是把这个慢衰减的 1/r 拆成两部分：短程部分用误差函数补余（erfc）的形式，衰减快，可以在实空间用截断半径直接算；长程部分光滑且周期性强，转到倒空间用 FFT 在网格上求和。两部分之和精确等于原始的库仑和，代价是多了一次网格化和一次傅里叶变换。
| 参数 | 取值 | 作用与调整建议 |
|---|---|---|
coulombtype | PME | 生物体系标配，别用 cut-off |
rcoulomb | 1.4 nm | 实空间静电截断，与 rvdw 保持一致 |
rvdw | 1.4 nm | 范德华截断，1.4 比常见的 1.0–1.2 更贵更稳 |
fourierspacing | 0.12 nm | 网格间距，越小越准越慢，0.12–0.16 常用 |
pme-order | 4 | 插值阶数，4 阶足够 |
ewald-rtol | 1e-5 | 倒空间相对误差容限 |
原文把库仑和范德华的截断都定在 1.4 nm，这个选择偏保守。截断从 1.2 提到 1.4，实空间邻接表变大，单步耗时明显上升，但换来的是界面上水分子取向更自然、蛋白表面带电残基的静电环境更接近真实。如果你的体系是小型球状蛋白、预算又紧，1.2 nm 也能跑；如果是膜体系或者高度带电的核酸复合物，1.4 nm 基本是底线。
3. 六个阶段的顺序账：为什么必须 EM → NVT → NPT → PMD
六阶段流程之所以是这个顺序，是因为每一步都在消掉一个特定的失稳源，而且这些失稳源之间有依赖关系：几何冲突不消掉，温度就控制不住；温度不先稳住，压力耦合就会把盒子推向震荡。跳过任何一步，后面都要用更大的代价补回来。
3.1 体系构建：pdb2gmx 补的氢、solvate 填的水、genion 撒的离子
体系构建这一步做了三件事，每一件都只保证"化学上说得通"，不保证"几何上不打架"。pdb2gmx按力场预设的几何构型给蛋白补氢，同时根据你选的质子化状态决定 His、Asp、Glu、Cys 上的氢怎么加——这一步的输出里如果有 warning 说某个残基的质子化状态需要人工确认，千万别略过，它直接影响后续的静电环境。
solvate是在格点上铺水分子，格点间距决定了初始水的密度和取向；铺出来的水在局部未必是最优氢键网络，有些分子会正好卡在蛋白表面的凹槽里，距离某几个重原子近到不合理。genion更粗暴，它在体系里随机挑位置把水替换成钠离子或氯离子来中和电荷，完全不考虑这些离子会不会正好贴着蛋白的带电残基。
# 补氢并按力场生成拓扑；-ignh 表示忽略输入 pdb 里已有的氢 gmx pdb2gmx -f complex.pdb -o complex.gro -water tip3p -ff amber14sb -ignh # 定义盒子：溶质距盒边至少 1.2 nm，十二面体比立方体省约 30% 的水 gmx editconf -f complex.gro -o box.gro -c -d 1.2 -bt dodecahedron # 铺水 gmx solvate -cp box.gro -cs spc216.gro -o solv.gro -p topol.top # 先 grompp 出一个 tpr，genion 需要它来读体系信息 gmx grompp -f em.mdp -c solv.gro -p topol.top -o ions.tpr # 中和电荷并补到 0.15 mol/L 的盐浓度 gmx genion -s ions.tpr -o ionized.gro -p topol.top -pname NA -nname CL -neutral -conc 0.15-d 1.2是溶质到盒子边界的最近距离，这个值影响两件事：太小会导致周期镜像之间直接发生非键相互作用，太大则水分子数暴涨、算力白白烧掉。-bt dodecahedron用截角八面体，同样的最小镜像距离下体积只有立方体的约 71%。-neutral只加刚好中和体系净电荷的离子，-conc 0.15在此之上再补到生理盐浓度，两个参数一起用才是"电中性 + 生理离子强度"。
注意genion执行后会修改topol.top里的分子计数，跑完要回头确认 SOL 的数量减少了对应的离子数。这一步没对上，后面 grompp 会报分子数不匹配，但报错信息往往指向别处，排查起来很绕。
3.2 能量最小化：目标不是全局极小，而是"能起跑"
前面三种操作留下的几何冲突，全靠能量最小化消掉。这里最容易误解的一点是：最小化不是要找到势能面的全局最低点，那个点在有限步数内根本到不了，而且也没有物理意义——真实体系在有限温度下本来就不在极小点上。目标只是把那些会让积分器发散的超近距离接触消掉，得到一个"力的大小在可接受范围内"的起始构型。
integrator = steep ; 最陡下降，沿负梯度方向走，稳但慢 emtol = 1000.0 ; 收敛判据：最大力小于 1000 kJ/mol/nm emstep = 0.01 ; 初始步长 nm nsteps = 50000 ; 步数上限，跑不完也没关系emtol是最大力判据，单位是 kJ/mol/nm，不是能量差——很多人第一次看这个参数会以为是能量收敛阈值。1000 这个值听起来很松，但它是"能安全起跑"的合理门槛；把emtol压到 10 以下，往往会看到 Fmax 卡在几百上下反复横跳，因为体系已经掉进某个局部极小，再优化也只是在谷底磨蹭，对后续毫无收益。
steep（最陡下降）的优点是每一步都保证能量不升，缺点是接近极小时步长越来越小、收敛极慢。常见的做法是先用steep跑几百到几千步把最恶劣的冲突碾平，再切到cg（共轭梯度）用历史梯度信息加速收尾。原文写的是"分别通过 Verlet 和 cg 算法"完成最小化，实际操作里就是这两段接力：steep负责活着走出危险区，cg负责把残余力压到判据以下。
判断最小化是否成功的标准，看 log 里的 Epot 曲线和 EMTEST 输出：Epot 应该单调下降然后趋平，最终 Epot 在 10⁵ 到 10⁶ kJ/mol 量级是正常的；如果 Epot 停在 -10⁷ 这种量级，说明还有原子重叠没解决，硬着头皮往下跑必炸。
3.3 NVT 先稳温度、NPT 再稳压力：顺序反了会怎样
最小化只优化了势能，没有动动能。要让体系进入设定温度，得给原子分配初速度，然后靠恒温器把温度维持在目标值附近。这一步放在 NVT 系综下做，盒子固定，只让温度到位、溶剂和离子围绕溶质弛豫。
为什么不直接上 NPT：恒压器的工作原理是根据瞬时压力的偏差去缩放盒子，而瞬时压力对坐标变化极其敏感。如果此时体系密度还偏离平衡值、溶剂还没弛豫、恒温器正在猛推速度，两个控制回路会互相激发——盒子缩放导致压力突变，压力突变又触发更大的缩放。原文对此的表述很克制："在恒压条件下，对一个可能远离平衡的体系，同时应用速度产生和坐标缩放可能导致体系不稳定。"说白了就是两个控制器打架，体系在几十皮秒内就能崩给你看。
; NVT 预平衡：盒子固定，只控温 integrator = md dt = 0.002 nsteps = 250000 ; 500 ps tcoupl = V-rescale ; 能产生正确正则系综的随机恒温器 tc-grps = Protein_LIG Water_and_ions tau-t = 0.1 0.1 ; 耦合时间常数 ps ref-t = 310 310 pcoupl = no ; 关键：这一步不控压 constraints = h-bonds define = -DPOSRES ; 对溶质重原子加位置限制define = -DPOSRES会去读拓扑里的posre.itp，对溶质重原子施加谐振位置限制，典型力常数 1000 kJ/mol/nm²。这个限制的作用是让溶剂和离子先围绕一个"冻住"的溶质弛豫——如果此时溶质可以自由移动，尚未平衡的溶剂会推着它乱跑，等溶剂平衡了，溶质已经到了一个扭曲的构象上。位置限制是预平衡阶段的后悔药，生产项开始前必须去掉。
恒温器的选择有讲究。Berendsen 恒温器温和、抗震荡，但它不产生正确的正则系综，只会把温度"拉"向目标值，动力学细节被平滑掉了，只适合早期预平衡。V-rescale 在 Berendsen 的基础上加了一个随机项来补偿能量，能严格产生正则系综，预平衡和生产项都能用。Nose-Hoover 更严格，但需要调耦合频率，参数没调好会出现温度振荡。
NPT 阶段开恒压器，pcoupl = Parrinello-Rahman或C-rescale，tau-p = 2.0ps，ref-p = 1.0bar。tau-p不能设得太小，小于 0.5 时盒子响应过快，会在目标体积附近来回震荡；也不能太大，否则压力平衡要跑很久。C-rescale 相比 Berendsen 的优势是能产生正确的等温等压系综，同时保持了 Berendsen 的数值稳定性，现在基本是生产项的首选。
这一阶段的验收标准很具体：温度曲线在 310 K 上下小幅波动（±5 K 以内），密度曲线在 1000 kg/m³ 附近收敛，压力在 1 bar 上下波动但不漂移。三个条件都满足，才谈得上进入生产项。
3.4 生产项：去掉位置限制，轨迹才开始有科学意义
预平衡结束后，溶剂弛豫好了，温度和压力都稳了，这时去掉POSRES，以无偏方式跑 MD 收集数据。所谓"无限制 MD"就是这个意思——溶质不再受任何人为约束，轨迹反映的是体系在给定系综下的真实涨落。
生产项的时长由科学问题决定，原文给的是 50 ns。这个长度对球状蛋白的构象稳定性研究够用，但对涉及大尺度构象变化或者慢解离过程的问题可能远远不够。轨迹输出频率要权衡：nstxout-compressed = 5000配合 2 fs 步长，就是每 10 ps 存一帧，50 ns 得到 5000 帧，对 RMSD、RMSF 这类分析绰绰有余。把输出频率设成每步都写，只会在 I/O 上烧掉大量时间，分析时还得再降采样。
能量和日志的输出频率同理，nstenergy = 5000、nstlog = 5000足够画出像样的收敛曲线。真出问题时你会回头翻 log，但不需要每步一行。
4. 数据分析：RMSD 平了就等于体系收敛了吗
跑完 50 ns 只是拿到了原料。数据分析是整条流程里最需要判断力的部分，也是最容易做成"跑完一套脚本、贴几张图"的部分。原文有一句话说得很直接：不应该简单地执行一个模拟，然后指望随着模拟进行会自动出现有趣的现象。哪些数据重要，必须在开跑之前就想清楚。
4.1 PBC 重包装：看起来碎掉的分子不是 bug
周期性边界条件让 mdrun 只计算每个原子与最近镜像的相互作用，但输出的坐标是原始坐标，不会自动帮你把跨边界的分子拼回去。于是你在可视化软件里会看到一个蛋白质被切成两半，一部分在盒子左边，一部分在右边；或者一个配体飘到盒子外面去了。
这在物理上完全没问题——对无限周期体系而言不存在"外面"——但在分析上必须处理，否则任何基于坐标的计算都会被这个人为的断裂污染。
# 第一步：把分子整体平移回盒子内，并让指定组居中 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-pbc mol按分子整体做平移，保证一个分子的所有原子一起移动，不会出现分子内部被拆散的情况；-center把选定的组（通常选 Protein）移到盒子中心；-ur compact把坐标压到紧凑表示，方便后续处理和可视化。这两条命令的顺序不能反：先做 pbc 和居中，再做-fit rot+trans。如果先叠合，跨边界的原子会把最小二乘拟合带偏，得到的旋转矩阵本身就是错的。
-pbc mol对完整分子有效，但如果体系里有被周期性边界切断的长期存在的东西（比如一条跨越盒子的膜），需要用-pbc whole或者专门的膜处理选项。选错选项的典型症状是轨迹里某些分子在相邻帧之间瞬移。
4.2 RMSD、RMSF、Rg：三个指标各自的盲区
RMSD 是把当前结构与参考结构做最小二乘叠合后算的坐标偏差，通常以时间为横轴。它的价值在于发现"跳变"，而不是确认"平稳"。原文给的例子很典型：蛋白和配体在经历初期短暂扰动后，大约在 8 ns、24 ns、200 ns、278 ns 处出现波动，但幅度很小，上下不超过 0.2 nm。这种量级的波动更像是侧链重排或溶剂重组的正常涨落，不是构象转变。
原文特意强调了两点误解：RMSD 是简化的度量，单个数值可以有多种解读；不能因为 RMSD 趋于平稳就宣布模拟收敛，因为结构相对于起始结构仍可能发生大幅变化。这句话背后的机制是：RMSD 是全局平均量，一个结构域的相对旋转可以被另一个结构域的反向运动抵消掉，曲线看起来平，实际构象空间还在漂。
import numpy as np def read_xvg(path): """读取 GROMACS 的 .xvg，跳过注释行与元数据行""" rows = [] with open(path) as f: for line in f: if line.startswith(("#", "@")) or not line.strip(): continue rows.append([float(x) for x in line.split()[:2]]) return np.array(rows) rmsd = read_xvg("rmsd.xvg") # 第 0 列是时间 ps，第 1 列是 RMSD nm time = rmsd[:, 0] / 1000.0 # 换成 ns 便于和论文图对齐 values = rmsd[:, 1] tail = values[len(values) // 3 * 2:] # 取后 1/3 轨迹作为平台期样本 print(f"plateau mean = {tail.mean():.3f} nm, std = {tail.std():.3f} nm")这段代码只做一件事：用后 1/3 轨迹的均值和标准差来量化"平坦"的程度，而不是用眼睛看末帧。均值和标准差一起看才有意义——均值稳定但标准差大，说明体系仍在大幅涨落；两者都小，才算真正进入平台期。这个判据比"曲线看起来平了"可靠，但它依然不能替代多指标交叉验证。
RMSF 衡量每个残基在模拟过程中的涨落幅度，必须基于已经叠合过的轨迹计算，且要按残基做平均而不是按原子。原文以 0.3 nm 为截止值标记高波动区域，得到 A 链的 2-11、36-94、105-107、244-257、482-508 区间，B 链的 2-14、26-35、44-101、245-263、464-530 区间。这些数字的意义不在绝对值，而在分布：高波动区通常对应溶剂暴露的 loop 区或末端，低波动区对应核心和二级结构。如果某个本应稳定的β折叠区突然出现在高波动列表里，那才是需要回头查轨迹的信号。
Rg（回转半径）衡量密实度，对折叠良好的球状蛋白基本是常数，因此它对"折叠是否稳定"这个问题并不敏感——原文明确说了这一点。它的价值在别处：当体系发生去折叠、结构域间相对运动或整体紧缩时，Rg 会给出 RMSD 看不出来的信息。原文的例子是 7 ns 和 24 ns 两处 Rg 有明显变化，正好对应 RMSD 的两次波动，这种多指标在同一时间点上的咬合，才是可信的构象变化证据，而不是某张图的孤立特征。
4.3 SASA、氢键数与 DSSP：用交叉验证替代单一指标
SASA 是溶剂可及表面积，用探针球在分子表面滚一遍算出来的。它常被用来判断蛋白是紧缩还是松散。原文的数据是 0–50 ns 内 SASA 波动下降，解释为蛋白逐步紧缩、突变有利。这个解释合理，但要加两个前提：SASA 对探针半径敏感（默认 0.14 nm，接近一个水分子半径），也对输出组敏感——只算蛋白和算蛋白加配体，数值差别不小，论文里必须写清楚算的是哪个组。
氢键数量用gmx hbond统计，原文的结果是蛋白-配体氢键数在 1–2 之间波动，没有重大变化。这个工具默认的几何判据是供体-受体距离 0.35 nm、角度 30°，判据一改，计数就变。所以这个数字的意义不在于"有几个氢键"，而在于"在同一判据下随时间是否稳定"。想拿它写进论文，判据参数最好一并列出。
DSSP 把每一帧的每个残基映射到二级结构类别，再统计各类别随时间的占比。原文对比了初始帧和末帧：coil 卷曲变少、β-sheet 变多、bend 变多、α-helix 变多，总体结构元件在向更有序的方向迁移。这种"占比迁移"比单帧截图有说服力得多，因为它排除了"某一帧恰好看起来规整"的偶然性。新版 GROMACS 里gmx dssp已经内置，老版本叫do_dssp且依赖外部可执行文件，跑之前先确认版本，否则会浪费半天在路径配置上。
4.4 MMPBSA：TOTAL 为负说明什么，残基分解为什么和稳态分析对不上
MMPBSA 用单轨迹法估计结合自由能：从生产项轨迹里抽一批帧，对每一帧分别计算复合物、受体、配体在溶剂化状态下的分子力学能量和溶剂化能，三者相减得到该帧的结合自由能，再对帧取平均。拆分出来的能量项是理解结果的关键：
| 能量项 | 含义 | 论文里的解读 |
|---|---|---|
| VDWAALS | 范德华能 | 负值表示疏水接触有利 |
| EEL | 静电能 | 负值表示静电吸引有利 |
| EGB | 极性溶剂化能 | 正值表示去溶剂化要付出代价 |
| ESURF | 非极性溶剂化能 | 通常略负，空腔项有利 |
| GGAS | 气相自由能 | 等于 VDWAALS + EEL |
| GSOLV | 总溶剂化能 | 等于 EGB + ESURF |
| TOTAL | 结合自由能 | 等于 GSOLV + GGAS |
原文的数据是 TOTAL 平均 -19.28 kcal/mol，负值表示结合有利。往下拆：VDWAALS 和 EEL 都是负的，说明疏水作用和静电作用共同驱动结合；GSOLV 是正的，主要来自 EGB 的去溶剂化惩罚——把配体从水里捞出来塞进口袋，必须破坏它原来的水化壳，这笔能量账要还；ESURF 略有利，因为它反映的是埋藏疏水表面减少空腔的贡献。
残基分解把总能量摊到每个残基上，看谁在帮忙、谁在拖后腿。原文的例子里，A 链 ALA-1 和 ARG-21 的贡献为正，是干扰项；PRO-23、MET-233、PRO-234 以及 B 链的 LYS-329、SER-350、GLU-353、GLN-354 贡献为负，有利于相互作用。
这里有一个非常值得记住的细节：残基分解的关键残基，和稳态姿态分析里靠氢键锁定小分子的残基（原文是 ARG-129 作为氢供体、LEU-352 作为氢受体）对不上，是完全正常的。原因不是算错了，而是两者用的数据范围不同——分解用的是整段轨迹的平均，稳态姿态用的是某一帧的代表性构象。一个讲全过程的热力学贡献，一个讲某一刻的结构特征，结论本来就不该一一对应。想通这一点，才不会在写结果的时候硬把两套数据凑成一致。
# 单轨迹法：-cg 依次给复合物组和受体组，配体由差集自动得到 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.dat -eo FINAL_RESULTS.csv-cs和-ct必须来自同一次运行，且轨迹必须已经做过 pbc 处理和叠合；-cg的两个组号要和index.ndx里的组号严格对应，组选错了会直接报找不到原子。帧间隔也要控制：50 ns 的轨迹如果每 100 ps 取一帧，就是 500 帧，这个量级既能体现构象分布的宽度，又不至于让计算排到天亮。
5. 常见问题排查：那些让体系在第一步就崩掉的细节
下面这些是跑分子动力学时反复出现的翻车现场，按阶段排。每一条都按"现象 → 原因 → 解决"来写，方便你直接在日志里对号入座。
5.1 体系构建与能量最小化阶段
现象一：grompp 或 mdrun 报 LINCS warning，log 里能量出现 -10⁷ 这种量级。原因通常是原子重叠——genion加的离子正好贴在某几个原子上，或者输入的 pdb 里有缺失原子被补到了重合位置，或者水分子被铺进了蛋白内部的空腔。解决方式是从源头堵：先确认pdb2gmx输出的 warning 都看过，再用较大的-d重铺盒子，离子浓度改用-conc 0手动放置，最后在最小化阶段确认 Epot 曲线单调下降到 10⁵–10⁶ 量级再往下走。
现象二：grompp报某个原子名找不到，或者拓扑里分子数不匹配。原因是小分子的 itp 原子命名和 gro 文件里的名字不一致，或者genion改过拓扑后没重新检查 SOL 计数。解决方式是把小分子生成流程里的命名规范固定下来，生成 itp 后手动 diff 一遍原子名；genion之后打开topol.top确认离子数和 SOL 减量对得上。
现象三：能量最小化跑满 50000 步还没到emtol，Fmax 停在几千。原因可能是emtol设得太严，体系已经进了局部极小；也可能是拓扑里存在无法通过位移消解的冲突，比如两个本不该成键的原子被错误地连上了。解决方式是先看 Epot 是否还在单调下降——如果已经趋平，把emtol放宽到 1000 就够；如果曲线在震荡，回去查拓扑和连接关系，别继续加步数硬扛。
5.2 平衡与生产阶段
现象一：NVT 阶段温度冲到几百 K 之后体系炸掉。原因常见于恒温器耦合时间常数设得太小（tau-t = 0.01），恒温器每步都在猛推速度；或者预平衡时没加POSRES，溶剂尚未弛豫就把溶质带跑。解决方式是把tau-t设到 0.1–1.0，先跑一段带位置限制的预平衡，让溶剂和离子安顿下来再放开。
现象二：NPT 阶段盒子持续收缩或者反复膨胀，密度不收敛。原因是初始密度偏离平衡值太多，或者tau-p设得太小导致恒压器响应过快。解决方式是在 NVT 阶段多跑一段把温度彻底稳住，NPT 用tau-p = 2.0，然后盯住密度曲线——它应该在 1000 kg/m³ 附近收敛，而不是单调地往上或往下爬。
现象三：生产项跑得异常慢，单步耗时远超预期。原因可能是截断半径设得过大、PME 网格过密，或者nstlist被手动固定成了不合适的值。解决方式是让cutoff-scheme = Verlet自动管理邻接表更新，把fourierspacing
本文还有配套的精品资源，点击获取
