简介：这份PDF围绕分子动力学模拟的完整流程展开，面向刚接触MD、需要搭建模拟体系并撰写论文方法部分的研究生与科研人员。内容从体系构建、加溶剂与带电离子、能量最小化讲起，依次梳理NVT与NPT预平衡、生产项模拟和数据分析，并提及ORCA、GROMACS等工具的使用场景；分析部分不仅覆盖重包装、RMSD、回旋半径、SASA与二级结构等指标，还提示在模拟前规划分析项目、避免仅凭单一指标判断收敛。资源包仅含1个PDF文件，大小约6.87MB，便于集中阅读和按流程查阅。已有444人学习，适合作为入门导读与方法写作的衔接材料，帮助读者把线性流程、专业术语和论文Method书写对应起来，形成可执行的学习路线。

1. 从一帧“炸开”的结构说起：MD 流程里最容易翻车的从来不是“跑”

分子动力学模拟最反直觉的地方在于：一次跑崩的模拟，很少是“跑”的那一步出的错。跑 50 ns 生产项轨迹，只是把前面每一步的误差放大——加氢时选错的质子化状态、溶剂化时没留够的盒子边距、能量最小化时没降下去的那个力，都会在几百皮秒之后以 LINCS 警告、NaN 坐标或者一条莫名断裂的轨迹的形式找上门。

真正决定成败的是流程里那些看不见的环节：每个原子上的力怎么算、动能怎么换算成温度、周期性边界怎么把无限体系塞进一个有限盒子。搞懂这条链路上的机制，才看得懂 RMSD 曲线为什么在 7 ns 和 24 ns 处拐弯，也才写得出方法部分那段看似平淡、实则每个参数都有出处的文字。

2. 体系构建与能量最小化：为什么“结构合理”比“能量最低”更要紧

把一管蛋白丢进水里，第一步永远不是让它动起来，而是让它“待得下去”。体系构建和能量最小化这两个阶段，本质都在回答同一个问题：任何一个原子，在当前位置上受到的其他原子给的力，是不是它能承受的。如果某个水分子被随手放在了离蛋白侧链 0.08 nm 的地方，这个力会大到足以在第一个积分步就把结构撕开。后面所有漂亮的 RMSD 曲线，前提都是这一关过了。

2.1 质子化状态：一个氢原子的取舍如何决定整条轨迹

力场要的是一套完整的、电荷自洽的原子拓扑，不是晶体结构里那份“缺氢”的坐标。pdb2gmx 这一步做的事情，是按照你指定的力场模板，把每个残基缺失的氢补齐，并顺便分配每类原子的电荷和原子类型。用-ignh是让工具丢掉从晶体结构里带过来的氢重新排布，因为 X 射线解析出的氢位置本来就不可靠，命名规范也未必和力场模板对得上。

真正的决策点在质子化状态。组氨酸有三种常见处理：HID（δ 位带氢，中性）、HIE（ε 位带氢，中性）、HIP（双质子化，带 +1 电荷）。同一个残基，选哪一种，侧链的电荷分布、能形成的氢键网络、甚至整个口袋的静电势都会变。天冬氨酸和谷氨酸默认去质子化（−1），赖氨酸和精氨酸默认质子化（+1），半胱氨酸要看是否形成二硫键。这些判断没有“标准答案”，只有当你在某个 pH 下做实验、或者有明确的催化机理证据时才有依据。

血泪经验是：这个小选择往往被当成事后才想起来确认的细节。等到轨迹跑到 20 ns，突然发现那个本该和配体成盐桥的残基带着相反的电荷，整段模拟的结论就得推倒重来，而这中间根本没有什么“后悔药”。

# 生成拓扑：指定 Amber 蛋白力场与 TIP3P 水模型 # -ignh 丢弃输入结构里的氢，由力场模板重新生成 # -his 可以交互指定每个组氨酸的质子化形式，避免默认值选错 gmx pdb2gmx -f complex.pdb -o complex.gro -ff amber14sb -water tip3p -ignh -his

这段命令的输出不只是坐标文件，还有topol.top和若干个#include的残基拓扑——后面所有力、能量、约束的计算，都是从这份拓扑里长出来的。-water tip3p必须和后面溶剂化用的水模型一致，混用两个水模型是新手最常埋下的第一颗雷。

2.2 溶剂化与离子中和：电荷为零不等于局部不打架

盒子怎么定，不是随手画个立方体。editconf -bt dodecahedron生成的十二面体盒，在同样的溶质-盒壁距离下体积只有立方体的约 0.71 倍，能省下近三成水分子，也就省下三成的计算量。参数-d 1.2指的是溶质原子到盒子最近壁面的距离，这个数字必须和后面实空间截断半径对得上：周期性边界下的最小镜像约定要求盒子任意方向的高度都大于两倍截断半径，否则一个原子会通过镜像和自己发生非物理的相互作用。

原文里静电和范德华的截断都取 1.4 nm，这意味着溶质到盒壁的距离至少也要到 1.4 nm 才安全。习惯上写 1.2 nm 也行，因为十二面体的最小像距离是按面间距算的，斜对角方向更宽裕，但取值低于截断半径就一定要警惕。

填水用的是预平衡的格点水盒子，工具会把整个溶质包起来再裁掉溢出的部分。接着是中和电荷：genion会根据拓扑里的总电荷数，自动算出需要多少个 Na⁺ 或 Cl⁻，然后从体系里随机挑水分子，把它们原地替换成离子。注意“随机”两个字——它只看电荷总量，完全不管这个离子落点会不会正好压在蛋白的酸性残基上、会不会和某个水分子位置重叠。

提示：genion的替换位置可以用-random控制随机种子。做同一体系的多条重复轨迹时，换一个种子能避免几条轨迹共享同一个初始离子构型而带来的系统性偏差。

# 定义十二面体盒子，溶质到壁面保留 1.4 nm，与截断半径一致 gmx editconf -f complex.gro -o box.gro -bt dodecahedron -d 1.4 -c # 用预平衡的格点水填充盒子，水模型由 -cs 指定，需与力场匹配 gmx solvate -cp box.gro -cs spc216.gro -o solv.gro -p topol.top # 中和体系总电荷，Na+ 替换水分子；-neutral 只补到电中性，不额外加盐 gmx genion -s ions.tpr -o ion.gro -p topol.top -pname NA -nname CL -neutral

小分子那一侧的准备逻辑是平行的，但更麻烦：它不是标准残基，力场模板里没有它的原子类型和电荷。通用的做法是先用一款量子化学软件在 B3LYP/6-31G* 水平上做几何优化，修正键长、键角、二面角，再用 RESP 方法拟合出一套固定电荷，最后交给通用力场（比如 GAFF2）去匹配原子类型。基组这件事属于领域共识，就像勾股定理，写进论文时不必单独引一篇文献来证明它存在。

2.3 能量最小化：最陡下降法与共轭梯度法分别在做什么

到这一步，体系里一定还藏着若干“不合理的接触”：pdb2gmx 按模板摆的氢原子、solvate 按格点填的水、genion 随机塞进去的离子，三者都只保证局部几何在自己的逻辑里成立，没考虑彼此。能量最小化的任务是把这些接触的力降到一个积分器能承受的量级。

原理上，每个原子受的力就是势能对坐标的负梯度：F = −∇U。最陡下降法沿这个方向挪一步，步长由当前最大的力分量归一化决定。它的好处是绝对稳健，地形再差也不会往上走；坏处是接近极小点时会在谷底来回震荡，收敛极慢。共轭梯度法把上一步的搜索方向按 Fletcher-Reeves 之类的公式揉进当前方向，避免重复走同一个方向，收敛快得多，但如果初始结构离极小点太远，它反而可能跑飞。

所以常规套路是两段式：先跑几百到几千步最陡下降，把最大的那几个力压下去，再切到共轭梯度收尾。原文里写的“最陡下降法，最大步数 50000 步”就是第一段的配置，收敛判据通常用emtol，单位是 kJ/(mol·nm)，1000 是个常见起点，想把结构压得更紧可以降到 100 甚至更低。

; 能量最小化：先最陡下降，把最大力压下去 integrator = steep emtol = 1000 ; 最大力收敛阈值，单位 kJ/(mol·nm) nsteps = 50000 ; 步数上限，不是必须跑满 nstlist = 10 cutoff-scheme = Verlet coulombtype = PME rcoulomb = 1.4 rvdw = 1.4

跑完之后最该看的不是“有没有跑满 50000 步”，而是最后一步的Fmax。它降到了emtol以下，说明结构已经站在一个局部极小点附近，动力学可以开始了；它卡在几千不动，说明体系里有个持续存在的坏接触，多半是某个残基侧链和离子或者水分子贴死了，得回去改构型，而不是加大步数硬扛。注意，这里找到的从来不是全局能量最低点，也不需要是——它的唯一职责，是给积分器一个不会立刻发散的起点。

3. 两段预平衡：NVT 稳温度、NPT 稳压力，顺序为什么不能颠倒

最小化结束时，体系在一个零温的势能面上待着，溶剂分子的取向、离子的分布，全是按各自模板摆出来的“自认为合理”的样子。预平衡要做的事，是让这些东西在有限温度下真正动起来并互相适应。原文把这段拆成 NVT 和 NPT 两步，并且强调顺序，这背后有明确的物理理由，不是流程上的仪式感。

3.1 温度不是设出来的，是速度算出来的

热力学温度在分子模拟里没有独立定义，它来自动能。对 N 个自由度的体系，T = 2·E_kin / (N_df·k_B)，其中 E_kin 是所有原子 ½mv² 的和。所以“把体系加热到 300 K”在代码里实际做的事，是按 Maxwell-Boltzmann 分布给每个原子抽一个初速度，让速度平方的统计平均刚好对应 300 K。

这也解释了一个常被忽略的细节：温度耦合的组怎么分。蛋白和溶剂（水加离子）通常分成两组，各自用tau-t控制耦合强度。tau-t是恒温器把温度拉回参考值的时间常数，取 0.1 ps 意味着几十到上百步内完成一次有效修正；取太大温度飘，取太小会人为压制体系的自然涨落。

恒温器的选择是有取舍的。Berendsen 恒温器实现简单、压制温度偏差极快，但它产生的不是标准正则系综，涨落被强行压制，用于预平衡无妨，用于需要严格系综平均的生产项就不合适。V-rescale 在 Berendsen 的基础上加了一个随机项补偿，能够正确采样正则系综，代价是要固定随机种子才能复现。Nose-Hoover 系综正确，但温度会以某个周期振荡，短时间窗口里看着不收敛，容易误判。

3.2 同时上恒温器和恒压器为什么会把体系拽散

刚完成最小化的体系，密度大概率是偏的——十二面体盒是按你给的-d算的，水位未必就是 1000 kg/m³ 该有的位置，内部还留着若干空洞或者过密的区域。恒压器的工作方式是按压缩系数缩放整个盒子的三个方向（或者各向同性缩放），把测到的瞬时压力推向参考压力 1 bar。

问题在于，压力本身是个涨落极大的量。一个 5 万原子的体系，瞬时压力在几百 bar 上下跳是常态，恒压器必须靠时间常数（tau-p，常取 2 ps）把它平滑掉。如果此时恒温器也刚启动，两组缩放同时在动坐标和速度，一个在改盒子尺寸、一个在改速度大小，两者的修正会互相叠加，初期的扰动被放大，体系就可能出现某个水分子被挤到蛋白里、进而触发 LINCS 报错。

正确的顺序是先只开恒温器。NVT 阶段体系体积固定，恒温器只需要处理一件事：把速度分布调整到目标温度，让溶剂分子在溶质周围重排、让离子扩散到能量更合适的位置。原文提到可以先加位置限制（-DPOSRES），把蛋白重原子用谐振子势拴在初始坐标上，这样预平衡期间溶质不会因为溶剂初期的冲击而大幅变形，等于把“溶剂弛豫”和“溶质构象变化”两件事分开处理。

注意：位置限制的力常数不能太大。1000 kJ/(mol·nm²) 是常见值，太大等于把蛋白冻成刚体，平衡出来的溶剂结构对放开限制后的构象并不适用；太小又起不到保护作用。

3.3 NPT 阶段与生产项：什么算“平衡好了”

NVT 把温度稳住之后，再叠上恒压器。这一阶段要盯的不是压力曲线本身，而是密度曲线。压力会一直抖，密度不会——密度从初始值逐步收敛到一个平台，上下波动在 1% 以内，才说明盒子尺寸和内部填充达到了一致。恒压器同样有取舍：Berendsen 恒压器收敛快，但系综不对；Parrinello-Rahman 系综正确，初期振荡大，得配合更长的tau-p；C-rescale 是折中方案，收敛和系综都说得过去。

; NPT 预平衡：温度继续耦合，同时打开压力耦合 define = -DPOSRES ; 保留位置限制，仅放开溶剂 tcoupl = V-rescale tc-grps = Protein Non-Protein tau-t = 0.1 0.1 ref-t = 300 300 pcoupl = C-rescale tau-p = 2.0 ref-p = 1.0 compressibility = 4.5e-5 ; 水的等温压缩系数，单位 bar^-1 gen-vel = no ; 速度从 NVT 的续跑结果继承 dt = 0.002 ; 2 fs 步长，依赖 LINCS 约束含氢键 constraints = h-bonds

生产项要做的第一件事就是去掉-DPOSRES，让蛋白彻底自由。这一步之后收集的轨迹才是可用于回答科学问题的数据。判断“能不能开跑”的几个硬指标：温度在目标值附近涨落且没有长期漂移、密度已经到平台、势能在预平衡最后一段的移动平均基本水平、RMSD 在最初几纳秒后进入平缓区。这四条里有任何一条不满足就开跑，后面的分析基本是白做。

至于 2 fs 的时间步长，它不是拍脑袋定的。含氢键的振动周期在 10 fs 量级，积分步长要到这个量级的十分之一以下才稳定，也就是 1 fs 左右；用 LINCS 把含氢的键长约束住，把这些最快的振动模式从积分里摘出去，步长才能放宽到 2 fs。如果不加约束硬上 2 fs，很快就能看到能量漂移。

4. 轨迹分析：RMSD、RMSF、Rg、SASA 各自在测什么

到这里模拟跑完了，硬盘上躺着一个.xtc或.trr。分析阶段的黑匣子感最强，因为每个指标都有现成命令，但每个指标也都能被读错。原文特意强调“不应该简单地执行一个模拟，并希望随着模拟进行会自动出现有趣的现象”，说的就是这件事：分析项目应该在开跑之前就设计好，而不是跑完了挨个指标试一遍看图说话。

4.1 PBC 重包装：为什么轨迹里的分子是碎的

先处理一个纯表示层的问题。mdrun 在计算力的时候用的是最小镜像约定，它知道每个原子的“真实”邻居是谁，但输出坐标时不做还原，直接写原始坐标。于是分子在盒子里扩散、跨越边界时，在可视化里就会看到蛋白被切成两半，一半在这头、一半在那头，或者整个分子“跑”到盒子外面去了。

这不是物理错误，无限体系里本来没有“外面”这个概念。但如果直接拿这种轨迹去做 RMSD 叠合，叠合算法会试图把两半分别塞回原位，结果完全失真。所以分析的第一步永远是重包装：以某个组（通常是蛋白）为中心，把分子拆散的部分按周期性边界重新接回一起。

# 先按分子完整性重包装，再以蛋白为中心居中 # -pbc mol 保证每个分子不被 PBC 切开，-center 把参考组移到盒子中央 gmx trjconv -s md.tpr -f md.xtc -o whole.xtc -pbc mol -center -ur compact # 可选：把体系压回盒子原点附近，避免长轨迹里坐标数值越滚越大 gmx trjconv -s md.tpr -f whole.xtc -o final.xtc -pbc whole

-pbc mol和-pbc whole的区别值得记一下：前者保证每个分子内部的原子连在一起，后者把分子当作整体移回主盒。做蛋白-配体复合物时，如果配体被单独切断，还要再加一个-pbc cluster或按残基索引单独做一次。这一步做没做，后面的 RMSD 曲线会长得完全不一样。

4.2 RMSD 与最小二乘叠合：一个被过度解读的数字

RMSD 的计算分两步：先把当前帧和参考结构做最小二乘叠合（求解一个最优旋转和平移，使对应原子的距离平方和最小，本质是 Kabsch 算法），再算叠合后对应原子距离的均方根。它衡量的是“构象差异”，不是“结构好坏”。

原文里那句提醒很关键：RMSD 是一个简化的度量，单一数值可以有多种解读，因此它本身并不足以说明问题。曲线进入平台，只能说明结构相对于起始构型不再发生大的整体偏移；它完全可能是一个持续缓慢漂移的结构，也可能是一个在两种构象之间反复切换、平均下来看起来稳定的结构。仅凭 RMSD 平台就宣布“模拟收敛”并据此终止，是这类工作里最常见的过度解读。

读法上，RMSD 更应该被当作一个“事件探测器”。原文提到的 7 ns 和 24 ns 两处波动，幅度不超过 0.2 nm，这种量级的转折才是值得回去看轨迹的地方——它对应当地发生了侧链重排、loop 区摆动还是配体轻微滑移，需要配合 Rg 和 RMSF 一起判断。

4.3 RMSF、Rg、SASA：三个互相印证的指标

RMSF 测的是每个原子（或残基）相对于时间平均位置的涨落幅度，它的时间平均对象是整条轨迹，而不是某一帧。所以它天然是“柔性图谱”：数值高的残基是 N 端、C 端、loop、以及没有配体压住的表面侧链。原文用 0.3 nm 作为截止值把高波动区段列出来，分别对应 A 链和 B 链上的若干区段，这个做法的意义在于把“哪些区域在动”变成一个可枚举的清单，再去和口袋附近的残基位置做对照。

Rg 是回旋半径，衡量的是原子相对质心的质量加权均方距离的平方根，物理含义就是“有多紧凑”。一个折叠良好的球蛋白，Rg 会在一个稳定值附近小幅涨落；Rg 持续上升意味着结构在伸展甚至去折叠。原文提到 Rg 在 7 ns 和 24 ns 各有一个波动，正好和 RMSD 的转折对得上，这种跨指标的对应关系才是可信的信号——单独一条曲线上的小起伏，很可能是噪声。

SASA 是溶剂可及表面积，用探针球在分子表面滚一遍算出来的暴露面积。原文观察到 SASA 在 50 ns 内波动下降，据此推断蛋白逐步紧缩。这个推断成立，但要注意它的限度：SASA 下降也可能只是配体盖住了某个区域，或者表面 loop 塌陷，不一定是整体构象变紧。所以 SASA 最好和 Rg 绑在一起看——两者同向变化，结论才稳。

4.4 DSSP 与氢键统计：二级结构是在数氢键模式

DSSP 不是“看形状”判断二级结构，它是按主链 C=O 与 N-H 之间形成的氢键模式做规则匹配：连续的 i→i+4 氢键定义 α 螺旋，i→i+3 定义 3₁₀ 螺旋，β 折叠靠桥接模式识别，剩下的归入 turn、bend、coil。所以一套 DSSP 时间序列，本质是一张“每个残基每一帧属于哪类二级结构”的表格，再做统计就能得到各类结构的占比随时间的变化。

原文对比初始帧和末尾帧，发现总的结构数量变少、coil 减少、β 折叠和 bend 增多、α 螺旋增多。这类描述在论文里很常见，但要多问一句：变化发生在哪些残基上？如果只是末端几个残基从 coil 变成了 bend，那属于边缘噪声；如果核心区的螺旋含量真的上升，那是构象重排。氢键数量（hbnum）统计同理，原文观察到在 1 到 2 之间波动，这个数字看着小，是因为它统计的是某个特定界面上的氢键，不是全体系——写的时候必须把统计范围交代清楚，否则读者会误以为是全蛋白氢键数。

4.5 MMPBSA 的能量分解为什么和稳态帧分析对不上

结合自由能用的是单轨迹法：从复合物轨迹里把受体和配体分别抽出来，当作两个游离态，用同一批构象去算三个状态的能量再做差。分子内项（键、角、二面角）在三个状态里因为坐标相同而基本抵消，所以实际算的是四块：范德华项 VDWAALS、静电项 EEl、极性溶剂化项 EGB、非极性溶剂化项 ESURF。气相自由能 GGAS = VDWAALS + EEl，溶剂化自由能 GSOLV = EGB + ESURF，总结合自由能 TOTAL = GGAS + GSOLV。

读这张表有个固定顺序。先看 TOTAL 的符号和量级，负值说明结合在能量上有利；原文给出的平均值约 −19.28 kcal/mol，这是一个偏强的数值，值得和实验亲和力做一次交叉核对。再看分解：VDWAALS 为负且贡献可观，说明疏水作用在起主要作用；EEl 为负但被 EGB 的正值大幅抵消，这是典型的“静电在真空中很强、在水中被去溶剂化惩罚吃掉大半”的图景；ESURF 通常是小负值，代表疏水表面积的贡献。

一个必然会被问到的问题：残基分解出来的关键残基，为什么和稳态帧里看到的那两个成氢键的残基对不上？原因很简单——残基分解是对整条轨迹的所有帧取平均，每个残基的贡献是“在所有构象里平均下来它帮了多少忙”；而稳态帧分析看的是某一个代表性构象里的具体相互作用。前者侧重过程与统计，后者侧重结果与结构细节。两者不矛盾，反而互为补充：稳态帧告诉你“靠什么稳住”，残基分解告诉你“整体上谁贡献最大”。

5. 常见问题与排查：五类翻车现场的现象、原因与解法

模拟跑不起来的时候，日志里的报错信息往往只指向最后一步，真正的原因在前面好几步。下面这几条是从构型到能量一路上最常撞见的，按现象、原因、解法来拆。

5.1 构建与最小化阶段：结构还没开始动就已经错了

第一类：pdb2gmx报残基未知或者原子名不匹配。现象是一条Residue 'XXX' not found in residue topology database，或者跑完了但拓扑里的原子数和坐标文件对不上。原因是输入的 PDB 里有非标准残基、缺失原子、或者原子命名不符合力场库里的大小写规范。解法是先用工具把缺失原子补全、把非标准残基的处理方式确认清楚（是当配体单独建模，还是当作修饰残基手工建模板），别指望加个-ignh就能糊过去。

第二类：能量最小化跑满了 50000 步，Fmax卡在几千下不来。现象是势能曲线在最后几千步几乎是平的，但最大力始终高于emtol。原因通常是体系里存在一个持续存在的坏接触，比如 genion 把钠离子塞进了某个酸性残基的氧原子旁边、或者晶体结构里带的结晶水和蛋白侧链重叠。解法是把最小化结束的结构导出来，按力的大小给原子排序，找出那几个“重灾区”残基，用可视化工具看一眼，手动挪开或者直接删掉冲突的水分子重跑。单纯把emtol调大是掩盖问题，不是解决问题。

第三类：最小化之后紧接着跑动力学，第一步就 LINCS 报错。现象是Step X, time Y: LINCS WARNING relative constraint deviation。原因是最小化没真正收敛，或者时间步长给得太大。解法是回头检查最小化的最终Fmax，确认它低于阈值；确认constraints = h-bonds已经打开；如果结构里有配体，还要确认配体的键参数在力场里是完整的，一个参数缺失的键会在积分时表现成无限硬的弹簧。

5.2 跑动阶段：能量曲线是最好的报警器

第四类：温度和压力曲线看着“不收敛”。现象是 NVT 阶段温度在 300 K 上下大幅跳，或者 NPT 阶段密度一直往下掉。前者的原因多半是tau-t取太小，恒温器强行压制自然涨落，或者体系的初始速度分布和实际自由度不匹配（比如把水分子当成刚体处理时自由度算错了）；后者的原因通常是位置限制没打开、溶质在压力耦合下被压变形，或者盒子尺寸本身就不够。解法是先看势能曲线的移动平均有没有水平，再看密度是否收敛到一个平台，别只盯着瞬时曲线做判断。

第五类：轨迹跑到一半突然出现 NaN 坐标或者体系“炸开”。现象是能量曲线某个点竖直上冲，之后所有分析都作废。原因可能是某个原子跑到了另一个原子的 LJ 势的极小点以内（近距离的排斥力按 1/r¹² 增长，一步就能把速度推到极大），也可能是恒压器在初期把盒子缩放得过头。解法是回到 NPT 的最后一段结构重新起跑，同时确认verlet-buffer-tolerance和邻域搜索列表的更新频率设置合理；如果问题反复出现，就要怀疑力场参数本身，尤其是配体的电荷和原子类型。

5.3 分析阶段：读错图比算错数更常见

第六类：RMSD 曲线出现莫名其妙的台阶。现象是曲线在某个时间点突然跳了 0.5 nm 以上再落回来。原因几乎一定是没做 PBC 重包装，叠合算法在某一帧把分子接到了错误的位置。解法是回到重包装那一步，确认-pbc mol加-center都做了，参考组选的是蛋白而不是整个体系。

第七类：把 RMSD 平台当成“收敛”的证明。现象是论文里写“模拟在 30 ns 后达到平衡”，审稿人问依据。原因就是前面说的，RMSD 是单一简化指标，能给出平台的结构有很多种。解法是补上至少两个独立证据：Rg 在同一时间段是否稳定、SASA 是否不再单调变化、关键残基之间的距离或二面角是否收敛。多条曲线指向同一个时间点，那句话才站得住。

6. 同一份轨迹两种结论：稳态帧与全轨迹平均的取舍

同一批模拟数据，用不同的取样方式能得到两套看起来互相矛盾的结论，这不是错误，而是方法本身的取舍。残基分解用的是整条轨迹的系综平均，回答的是“在整个过程中，哪个残基平均贡献了最多的结合能”；稳态帧分析挑的是能量最稳的那一段或者代表性构象，回答的是“在这个具体结构里，谁和谁通过什么方式咬在一起”。前者对构象采样敏感、数值平滑但丢失细节，后者细节清晰但依赖你挑的那一帧是否具有代表性。写结果的时候把两套结果并列呈现、分别标注取样范围，比强行让它们一致要诚实得多，也更经得起追问。

第二个取舍在自由能计算的采样策略上。单轨迹法只用复合物一条轨迹，受体和配体的游离态构象直接从复合物里剥离，好处是分子内项大部分抵消、统计误差小、算得快，代价是它假设结合不引起显著的构象重排。如果这个体系的结合明显伴随口袋的开合，单轨迹法会系统性低估构象熵的代价。多轨迹法分别对三个状态各跑一条轨迹，能捕捉构象变化，但引入的统计噪声大得多，收敛所需的采样量成倍上升。

第三个取舍在写作顺序和操作顺序之间。实际操作里，能量最小化之后通常是先 NVT 稳温度、再 NPT 稳压力，然后才是生产项；但方法部分为了叙述连贯，常常把两个系综的平衡合并成一句话带过。这不影响结果，前提是你自己清楚每一步的真实执行顺序和它解决了什么问题，否则在被问到“为什么先 NVT”时会答不上来。参数也一样：截断半径、步长、恒温器类型这些，写进论文时每一个都要能对应到一次具体的选择，而不是从模板里抄来的默认值。

从那以后我每次开跑新体系，都会强制走一遍固定流程：先确认质子化状态和力场一致，再看最小化的最终 Fmax，再盯 NPT 的密度平台，最后分析前一定先做 PBC 重包装。这四步里任何一步没打勾，后面所有的曲线都不看。希望帮到你。

本文还有配套的精品资源，点击获取