有人发来一张温度曲线，最后一个点在 295 K 左右，目标温度是 300 K。只看末帧，似乎体系已经足够接近目标。我们用 GROMACS 跑了一个仅有二百二十一枚水分子的盒子，最后一帧正好符合这个描述，但第一皮秒的平均温度只有约 250 K，整个短窗口还经历了很大的压力变化。这个小例子可以练习三个常被混在一起的判断。最小化是否达到预先设定的力阈值，控温是否经过明显的初始暂态，以及观测长度是否足以支持准备回答的物理问题。
从坐标数回到拓扑中的分子数
准备体系时，坐标文件给出原子的位置，拓扑描述分子由哪些原子组成、怎样相互作用，参数文件决定积分、边界条件和输出频率。运行预处理时，三者要放在一起核对。用自带水盒填充一个边长二纳米的立方盒，并让溶剂化工具同步更新拓扑的分子段。如果手工复制了溶剂化前的拓扑，就可能看到水分子已经挤满盒子，预处理却报告坐标与拓扑的原子总数不符。处理时先查分子表和水模型的原子数，再检查是否引用了旧文件。通过提高允许警告的数量来继续，不会纠正原子排列或相互作用的错配。GROMACS 预处理器说明
这次用随软件分发的 SPC 刚性水拓扑和水盒模板。填充后数得六百六十三个原子，恰好等于二百二十一枚三位点水分子。三条盒边都是二纳米，短程截断设为零点八纳米，小于最短盒边的一半，长程静电用粒子网格 Ewald 求和。将来替换成蛋白、配体或膜时，还要核对溶质与水模型的参数兼容性、溶质与其周期镜像的距离、离子数及体系电荷。仅凭输出文件中写着“水盒”，还无法知道是否使用了预期的组合。
| 文件或量 | 本次实际读到什么 | 可以回答什么 |
|---|---|---|
| 最小化日志 | 32 步达到最大力阈值 | 几何松弛是否符合本次停止条件 |
| 能量文件中的温度 | 101 个导出时点 | 不同时段是否经历温度暂态 |
| 恒容阶段的压力 | 从约 −5870 到 142 bar | 固定体积下这段短轨迹的波动 |
| 导出的末帧 | 约 295.36 K | 该时刻的状态，仍需前面的时间序列 |
实跑经验｜同一轨迹，换一个时间窗就换了一个结论
我们在 Ubuntu 安装 GROMACS 2023.3，用一个进程、一个线程完成溶剂化、能量最小化和恒粒子数、恒体积、恒温的短动力学运行。最小化采用最速下降法，将最大力阈值设为每纳米每摩尔二百千焦。日志记录第三十二步达到停止条件，最终最大力约为 159.8，落在编号三百零四的原子上。导出的首末势能约为每摩尔负一千九百五十三和负九千六百五十九千焦。这里应一并报告停止原因和最大力。单看势能明显下降，无法区分达到阈值与用尽允许步数两种情形。
随后以两飞秒为步长积分两千步，全程四皮秒，温度耦合使用速度重标定，目标三百开尔文，关闭压力耦合。每二十步写一次能量，包含零时刻共得到一百零一个温度记录。从这些导出记录作等权算术平均，整段约为 281.29 K。若把前四皮秒依次按左闭右开的一皮秒窗口分成四段，每段二十五个时点，均值依次约为 249.80、295.21、294.82 和 284.76 K。四皮秒处的最终记录单独为 295.36 K，未计入上述四个分段。说明窗口端点和采样间隔，别人才能从同一份文件得到同一组数字。
# md.edr 是动力学运行的能量文件。按名称选择，避免照搬另一体系的菜单编号。
printf 'Temperature\n0\n' | gmx energy -f md.edr -o temperature.xvg -xvg none
printf 'Pressure\n0\n' | gmx energy -f md.edr -o pressure.xvg -xvg none
# 计算第一个左闭右开窗口的导出样本均值，xvg 第一列为 ps
awk '$1 >= 0 && $1 < 1 {sum += $2; n++} END {print n, sum/n}' temperature.xvg
命令行的能量工具还能读取储存在能量文件里的累积统计量，打印的统计值未必等于屏幕上有限个导出点的手算平均。本文的数字明确来自导出点，方便读者直接复查。如果自己的均值对不上，先核对使用的是程序打印的平均、导出记录的平均还是选定区间的数据，同时检查有没有把零时刻和最后一帧都算了两次。能量工具对统计量和可导出帧的说明
还可以用这份水盒练习审核交接记录。最小化的输出坐标要成为动力学的输入，初次赋予速度的参数应与随后从检查点续跑的参数分开保存。延长一条已有轨迹时，如果又生成了一套初速度，之前记录的轨迹和新段就不能当成同一条连续演化来解释。把最小化输入、动力学输入和相应检查点分别命名，再比较最后一个时间戳与下一段第一个时间戳，比看到两个同名文件后凭修改日期拼接可靠。
压力均值为什么不能拿来证明已经平衡
本例恒容阶段逐点压力最低约为负五千八百七十、最高约为一百四十二巴，导出点的算术平均约为负一千五百八十八巴。这个模拟没有启用压力耦合，也没有用足够长的平衡轨迹来估计体相性质。把负值直接解释成某种生物机制，或者据它反推真实液态水在室温下的平衡压力，都会高估这个算例能提供的证据。排查真实体系可以先问，压力是受控量还是仅记录的量，盒体积是否允许变化，温度轨迹是否正经历暂态。随后再选合适的平衡方案并监测与研究问题相关的量。GROMACS 动力学与温压耦合算法
轨迹长度是否够，要回到要估计的量。只想确认拓扑、单位和输出频率是否设对，短水盒很合适。要报告蛋白构象转换、配体结合稳定性或密度，就应检查相应观测量随时间的变化、相邻记录的相关性和独立重复之间的差异。把同一段轨迹保存得更密，只会提高取样频率，无法自动增加独立样本。一个可交接的做法，是同时留存参数、预处理后的运行输入、最小化日志和能量文件，并为每个统计量写清时间窗和单位。若下一位同事只拿到经裁剪的曲线图片，他很难判断曲线之前发生过什么。
扩大到多种构象和多个独立重复前，先拿一个代表体系确认输入、输出频率与平衡判据。需要集中运行这些轨迹、保存中间结果时，可了解云生信的生物信息分析服务器。
