GROMACS 积分器与速度定义实操

多数普通 MD 任务从 integrator=md 开始即可。选择积分器时，除了看能不能运行，还要理解输出速度在哪个时刻、动能怎样计算，以及能与哪些温压耦合方式组合。

1. 从同一套输入建立独立分支#

本例使用同一已平衡 SPC/E 水盒：边长 2.5 nm，510 个水分子、1530 个原子。生产阶段用 NVE，每组 20 ps；比较 md、md-vv、md-vv-avek，各用 0.5、1、2 fs 三个步长。

在分配好的任务里，分别执行：

bash run.sh md-dt.001
bash run.sh md-vv-dt.001
bash run.sh md-vv-avek-dt.001

每条命令创建自己的 work/名称，不覆盖其他分支。程序先用该分支 MDP、原拓扑和原检查点生成 TPR，再积分。核心步骤为：

gmx grompp -f md.mdp -c start.gro -t start.cpt \
  -p topol.top -o md.tpr
gmx mdrun -deffnm md -ntmpi 1 -ntomp 2 \
  -nb gpu -pme gpu -bonded cpu -update cpu -pin off -notunepme

这些不同积分器读取相同的数值坐标和速度，但速度对应的时间点可能不同，不能称为完全相同的严格动力学起点；这里比较各自实际输出和能量变化，不要求轨迹逐帧重叠。

integrator = md
dt = .001
nsteps = 20000
continuation = yes
tcoupl = no
pcoupl = no
gen-vel = no
comm-mode = Linear
nstcomm = 100
nstenergy = 200
nstcalcenergy = 200
nstlog = 1000
nstxout = 200
nstvout = 200
nstfout = 0
nstxout-compressed = 200
compressed-x-precision = 1000
cutoff-scheme = Verlet
nstlist = 20
verlet-buffer-tolerance = -1
rlist = 1.0
coulombtype = PME
rcoulomb = .8
vdwtype = Cut-off
vdw-modifier = Potential-shift
rvdw = .8
DispCorr = no
fourierspacing = .12
pme-order = 4
pbc = xyz
constraints = h-bonds
constraint-algorithm = lincs
lincs-order = 4
lincs-iter = 1

dt 单位是 ps。1 fs 为 0.001 ps，20 ps 对应 20000 步。gen-vel=no 不重新随机生成速度，continuation=yes 读取已经满足约束的检查点状态。TRR 每 0.2 ps 保存坐标和速度。

2. 三种积分器怎样选#

本版 md-vv-avek 即使做 NVE，也要求输入中明确 nsttcouple=1、nstpcouple=1；下载包已写好。这不是给 NVE 打开恒温或恒压，实际仍为 tcoupl=no、pcoupl=no。

表中的能量是整个水盒的 kJ/mol 数值，不是单个水分子的值。图按原子数归一化便于阅读。这里短轨迹中的能量变化没有随步长严格单调减小，不能据此拟合一个通用误差阶数；输入起点、不同速度时刻、固定截断/缓冲、约束和混合精度都参与了结果。正式 NVE 仍需用自己的体系检查能量漂移和步长。

3. 为什么 TRR 速度算出的动能可能不同#

printf 'Potential\nKinetic-En.\nTotal-Energy\nTemperature\nPressure\n0\n' | \
  gmx energy -f work/md-vv-dt.001/md.edr -o vv-energy.xvg
bash analyze.sh work

脚本从原生能量文件取动能，并从 TRR 速度另算 0.5×Σmv²。md-vv 的整步速度结果在本例中仅差约 0.002 kJ/mol 量级；md 的半步速度和 md-vv-avek 的半步平均定义不能这样直接替换，差异可达数 kJ/mol。不要把这种定义差异误判成“能量文件错误”。

4. SD 与 BD 不是给普通 MD 换个名字#

bash run.sh water-sd
bash run.sh water-bd

sd 是带惯性的随机动力学，tau-t=2 ps 控制摩擦相关时间，ld-seed 固定随机流。bd 是 Brownian/位置 Langevin 动力学，本例显式给 bd-fric=1000；它的速度、摩擦和动力学含义与普通 NVE 不同。不要要求两者逐帧相同，也不要把 BD 当作保守动力学来检查总能量守恒。

integrator = sd
dt = .001
nsteps = 20000
continuation = yes
tcoupl = no
pcoupl = no
gen-vel = no
comm-mode = Linear
nstcomm = 100
nstenergy = 200
nstcalcenergy = 200
nstlog = 1000
nstxout = 200
nstvout = 200
nstfout = 0
nstxout-compressed = 200
compressed-x-precision = 1000
cutoff-scheme = Verlet
nstlist = 20
verlet-buffer-tolerance = -1
rlist = 1.0
coulombtype = PME
rcoulomb = .8
vdwtype = Cut-off
vdw-modifier = Potential-shift
rvdw = .8
DispCorr = no
fourierspacing = .12
pme-order = 4
pbc = xyz
constraints = h-bonds
constraint-algorithm = lincs
lincs-order = 4
lincs-iter = 1
tc-grps = System
tau-t = 2
ref-t = 300
ld-seed = 2612201

integrator = bd
dt = .001
nsteps = 20000
continuation = yes
tcoupl = no
pcoupl = no
gen-vel = no
comm-mode = Linear
nstcomm = 100
nstenergy = 200
nstcalcenergy = 200
nstlog = 1000
nstxout = 200
nstvout = 200
nstfout = 0
nstxout-compressed = 200
compressed-x-precision = 1000
cutoff-scheme = Verlet
nstlist = 20
verlet-buffer-tolerance = -1
rlist = 1.0
coulombtype = PME
rcoulomb = .8
vdwtype = Cut-off
vdw-modifier = Potential-shift
rvdw = .8
DispCorr = no
fourierspacing = .12
pme-order = 4
pbc = xyz
constraints = h-bonds
constraint-algorithm = lincs
lincs-order = 6
lincs-iter = 1
tc-grps = System
tau-t = 2
ref-t = 300
bd-fric = 1000
ld-seed = 2612202

本例人为选择了摩擦参数，只演示输入与输出；不能用这 20 ps 报告真实水的扩散常数。

官方入口：积分器与 MDP 选项。

下载与重新运行#

下载完整输入、实际结果、MDP、CSV和五张图的绘图脚本。本包直接提供已经平衡的 510 分子 SPC/E 水盒和 1UBQ 蛋白水溶液输入；柔性水分支另有 10 ps 准备阶段。先查看现成数据：

bash analyze.sh

结果在 reference/analysis，图片在 reference/figures。在新解压目录提交整套新计算：

bash submit.sh
squeue -u "$USER"

作业使用实际分区 gpu-8p100-hs，每组申请一张 P100 和两个 CPU 核；双精度对照使用该任务的 CPU 核计算，不使用 GPU 力计算。22 个任务脚本包括一个柔性水准备任务、20 个生产条件和一个后处理任务，依赖关系已写好。需要先看脚本而不提交时，在另一个新解压目录用 bash submit.sh --dry-run。jobs、work 已存在时不重复提交整套任务。

正文中的单条 bash run.sh ... 只在已经分配到对应资源的任务中执行，不在登录/管理节点运行 MD。输入、拓扑、检查点和力场要配套；不要拿重新建模后的拓扑继续读取旧检查点。

混合精度 GPU 路径使用 gromacs/2026.4-p100-threadmpi：GROMACS 2026.4、GNU13.3、CUDA12.8、FFTW3.3.10、thread-MPI/OpenMP。双精度路径使用 gromacs/2026.4-cpu-double，命令 gmx_d，固定 v2026.4 源码构建，程序标记 2026.4-dev，CPU-only、GNU13.3、AVX2、FFTPACK。不是 Intel Xeon 8173M 的性能数据。

提取使用 gmx_d，绘图用 NumPy/Matplotlib；结构与轨迹也通过 MDAnalysis 核对。脚本保留双精度 TRR 的实际精度，避免把双精度坐标先缩成 float32 再判断极紧约束误差。参考数据均为真实计算，但 20 ps 水与 5 ps 蛋白仅用于学会操作，不作为收敛热力学、扩散或性能结论。