GROMACS Colvars距离与谐约束

希望跟踪两个原子或原子组的距离，或者给这个距离施加谐约束，可以使用GROMACS内置的Colvars接口。本篇从“只监测”开始，再比较两个约束强度。选用水中两颗氧原子作清楚、短小的操作示例，不把它解释为蛋白质反应坐标或自由能结果。

下载完整输入、真实轨迹与绘图资料包，约147 MiB。资料包含32个短算例及6次同构型电场力对照。解压进入资料包目录，先查看已有结果：

bash view.sh

四张轨迹数据图位于 analysis/figures，电场力图位于 field-analysis/figures；每个条件的CSV和数值汇总位于 analysis/条件名。脚本调用集群已有MDAnalysis2.9.0、NumPy、Matplotlib环境，画已有数据时不提交计算。

本例使用348个SPC/E水分子，共1044个原子，初始坐标来自已有水盒的NVT平衡末帧。新计算重新生成速度，300 K、固定盒子NVT、步长1 fs、长度20 ps，每0.1 ps保存一帧，共201帧完整位置/速度TRR。两种速度种子13579、24680分别在P100和V100 16G计算。不同平台重复相同种子，不当作四个独立初始种子。

GROMACS2026.4 thread-MPI，GNU13.3、CUDA12.8、FFTW3.3.10。短程非键和PME使用GPU，键相互作用及坐标更新使用CPU；Colvars集体变量的计算不因此变成GPU内核。

1 先确认模块，再准备和提交#

module load gromacs/2026.4-p100-threadmpi
gmx --version
python3 prepare.py
cd cases/p100-colvars-K100-s13579
cat colvars.in
sbatch job.slurm
squeue -u "$USER"

版本输出应有 Colvars support: enabled。本集群版本显示Colvars2025-10-13；V100入口为 gromacs/2026.4-v100-threadmpi，进入对应的 v100-... 条件目录。附件另有 colvars-monitor、colvars-K400 两类配置。

本例流程：复制已有水盒坐标 → 选择距离与约束强度 → grompp把Colvars配置写入TPR → 20 ps NVT → 读取CV、偏置力及能量 → 画图对照。

模板每项申请1张对应GPU和1个CPU核心，使用核时分区；套餐用户按账户规则修改分区。不要在登录节点直接运行完整MD。已有 cases 或结果目录时，重新解压到新目录准备，避免覆盖。

任务实际执行的主命令如下，供已分配资源中的手动操作参考：

gmx grompp -f md.mdp -c initial.gro -p topol.top -o md.tpr   -po md-processed.mdp -pp md-processed.top
gmx mdrun -deffnm md -ntmpi 1 -ntomp 1 -nb gpu -pme gpu   -bonded cpu -update cpu

拓扑完整分子表如下；原子电荷、刚性水几何由所引用的SPC/E模型给出：

#include "oplsaa.ff/forcefield.itp"
#include "oplsaa.ff/spce.itp"

[ system ]
SPC/E liquid transport tutorial in water

[ molecules ]
SOL               348

2 Colvars输入怎么写#

colvars.in完整内容：

colvarsTrajFrequency 100
colvar {
 name d
 width 1.0
 outputAppliedForce on
 distance {
  group1 { atomNumbers 1 }
  group2 { atomNumbers 4 }
 }
}
harmonic {
 name restrain
 colvars d
 centers .5
 forceConstant 100
 outputEnergy on
}

atomNumbers 1 和 4 是从1开始编号的两颗水氧。单原子组的质心就是它本身；换成多原子组时，要核对组的定义和质量。距离按周期边界的最小镜像处理。

centers .5 为0.5 nm，colvarsTrajFrequency 100 每100步输出一次，1 fs步长下间隔0.1 ps。outputAppliedForce输出该集体变量上的偏置力，outputEnergy输出谐偏置能。

这里 width 1.0 固定为1 nm，forceConstant 100 对应有效刚度100 kJ/(mol·nm²)。改变width时，不能继续把同一个forceConstant数值当成相同有效刚度，需按Colvars的width归一约定重新换算。

只监测时，保留 colvar 块，删除整个 harmonic 块；附件的monitor条件已经这样准备。K400条件只把forceConstant改成400。这里没有metadynamics、ABF或PMF重建。

3 完整MDP#

integrator = md
dt = .001
nsteps = 20000
cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
rcoulomb = .9
rvdw = .9
DispCorr = EnerPres
fourierspacing = .12
pme-order = 4
pbc = xyz
constraints = none
nstenergy = 100
nstcalcenergy = 100
nstlog = 1000
nstxout = 100
nstvout = 100
nstfout = 0
nstxout-compressed = 0
tcoupl = v-rescale
tc-grps = System
tau-t = .5
ref-t = 300
pcoupl = no
comm-mode = Linear
nstcomm = 100
continuation = no
gen-vel = yes
gen-temp = 300
gen-seed = 13579
ld-seed = 13680
colvars-active = yes
colvars-configfile = colvars.in

最后两个选项启用Colvars并指定配置文件。grompp时要能找到该文件，内容会随TPR保存。水的SETTLE写在拓扑内，constraints = none 不会关闭拓扑中已有的刚性水定义。

4 任务结束后提取结果#

在当前算例目录执行：

module load gromacs/2026.4-p100-threadmpi
python3 ../../export_native.py
ls native-analysis

V100更换为对应模块。这个步骤读取已有轨迹和能量，输出到新建的 native-analysis，不会覆盖同名分析目录。全部参考结果及绘图脚本可在资料包根目录运行 bash view.sh 查看。

md.colvars.traj表头说明各列：step、d、fa_d；有谐约束时还包括E_restrain。step是步数，本例乘0.001得到ps。偏置力和能量应满足：

F_bias = -K * (d - 0.5)
V_bias = 0.5 * K * (d - 0.5)^2

这是集体变量上的偏置力，不是某颗原子的总受力。监测模式的偏置力为0，距离仍会因真实分子运动改变。

在当前条件目录画自己的新结果：

module load ase/3.29.0
python - <<'PY'
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
d=np.loadtxt("md.colvars.traj",comments="#",ndmin=2)
plt.plot(d[:,0]*0.001,d[:,1])
plt.xlabel("Time (ps)");plt.ylabel("O1-O4 distance (nm)")
plt.tight_layout();plt.savefig("new-colvars-distance.png",dpi=180)
PY

这段程序直接读取刚生成的原生数据；PNG保存在当前目录。XVG以 # 和 @ 开头的行是说明/列名，不能当作数值读入。

5 实际数据与图#

K400把距离限制在0.5 nm附近；K100的涨落更宽。监测模式没有把距离拉向0.5 nm的偏置。表中是有限窗口平均值，不能把它直接当成无偏距离分布或自由能。

每个条件的完整轨迹和输入都随附件保留。更复杂的变量、偏置和收敛分析，应按Colvars GROMACS手册选择，并核对GROMACS接口说明。

6 Colvars、PLUMED和QM/MM怎样选#

PLUMED的实际输入和提交见GROMACS PLUMED偏置与集体变量，更多偏置方法见PLUMED专题。普通thread-MPI模块未启用PLUMED接口；加载独立PLUMED本身不会自动改变GROMACS的编译能力。独立CP2K可用也不代表GROMACS已链接CP2K。

配置文件能被 grompp 读取，不等于运行阶段已支持所选QM/MM接口。使用专用联用前，先核对实际模块的编译支持与对应接口文档。