PLUMED实操：集体变量、偏置与后处理
从水体系的真实GPU短任务开始，准备plumed.dat、提交GROMACS任务、查看COLVAR，再做离线分析。
阅读本篇目录
PLUMED用来定义和分析集体变量，也可以给分子动力学施加偏置。本专题先用一个可以完整复跑的小例子讲清楚操作顺序：准备输入 → 提交原生MD → 查看COLVAR → driver离线重算 → 检查偏置能和力。更复杂的自由能计算需要另外设置采样与统计。
1 本例算什么#
体系为884个TIP3P水分子，共2652个原子。使用GROMACS 2026.4外部MPI+PLUMED构建，分别在P100和V100 16G运行100步，时间步长0.002 ps，总时长0.2 ps。软件采用GNU GCC/G++ 13.3.0、Open MPI 5.0.10、CUDA 12.8和PLUMED 2.9.4。
两个条件使用不同随机种子，演示同一套操作流程；不把两条轨迹作硬件逐帧相等或性能评测。100步也不足以给出平衡分布、结合自由能或收敛PMF。
图中仅画出原生起点的三个水分子，原子编号仍使用完整体系的编号。蓝线是氧1和氧4的距离，橙线标明二面角3–1–4–5。距离初始约0.998 nm，目标设为0.5 nm；这是偏置操作示范，不是氢键水对结合自由能模型。
2 先建立新算例#
将下载文件保存为plumed-example.zip，解压后选择P100示例：
unzip plumed-example.zip -d plumed-example
cd plumed-example
python new-case.py p100
cd work/p100
新目录包含md.gro、topol.top、run.mdp、plumed.dat和run.sh。如果已有同名工作目录，程序会退出；保留原结果，在新的解压目录重建即可。V100示例把p100改成v100。
cat plumed.dat
cat run.mdp
sbatch run.sh
提交需要对应GPU资源权限。run.sh通过Slurm申请一张GPU和两个CPU线程，不要在登录节点直接执行完整MD。
3 plumed.dat逐行解释#
UNITS LENGTH=nm TIME=ps ENERGY=kj/mol
d: DISTANCE ATOMS=1,4
a: ANGLE ATOMS=2,1,3
t: TORSION ATOMS=3,1,4,5
c: COORDINATION GROUPA=1 GROUPB=4,7 R_0=0.35 NN=6 MM=12
r: RESTRAINT ARG=d AT=0.5 KAPPA=100
PRINT ARG=d,a,t,c,r.bias STRIDE=1 FILE=COLVAR FMT=%24.16e
DUMPATOMS ATOMS=1-7 FILE=native-positions.xyz STRIDE=1 PRECISION=16 UNITS=nm
| 设置 | 本例含义 |
|---|---|
UNITS |
长度nm、时间ps、能量kJ/mol；kj/mol按原拼写使用 |
DISTANCE |
氧1到氧4的最小镜像距离 |
ANGLE |
原子2–1–3夹角，中间的1是顶点；输出弧度 |
TORSION |
3–1–4–5的有符号二面角，输出弧度 |
COORDINATION |
氧1与氧4、7之间的平滑接触数；R_0是开关函数尺度，不是硬切断半径 |
RESTRAINT |
只对距离d加谐势，目标0.5 nm，力常数100 kJ·mol⁻¹·nm⁻² |
PRINT |
每步打印这几个值和偏置能，保留较多有效数字 |
DUMPATOMS |
每步保存前7个原子的高精度nm坐标，供离线对照 |
角度和配位数在这里用于观察，没有额外加偏置。二面角原子顺序会影响符号；所选两平面还应避免接近共线，以免坐标舍入引起明显角度差。
4 原生程序如何启动#
提交脚本已写好以下步骤，普通用户通常只需sbatch run.sh：
module load gromacs/2026.4-p100-mpi-plumed
gmx_mpi grompp -f run.mdp -c md.gro -p topol.top -o run.tpr
mpirun -np 1 gmx_mpi mdrun -deffnm run -ntomp 2 -nb gpu -pme gpu -update cpu -plumed plumed.dat
这是外部MPI构建，使用gmx_mpi和mpirun，不要再加thread-MPI的-ntmpi。V100脚本加载对应的2026.4-v100-mpi-plumed模块。非键和PME走GPU，坐标更新在CPU，与当前PLUMED联用路径匹配。
本例run.mdp的关键设置为：
integrator = md
dt = 0.002
nsteps = 100
constraints = h-bonds
tcoupl = v-rescale
tc-grps = System
tau-t = 0.5
ref-t = 300
pcoupl = no
nstxout = 1
nstvout = 0
nstfout = 0
完整文件还包含截断、PME、随机种子和输出设置，直接使用下载包里的完整输入。topol.top引用已安装GROMACS的amber99sb-ildn.ff目录中的TIP3P定义，不需要用户重新编写水力场。
5 最后查看哪些文件#
head COLVAR
tail COLVAR
ls run.trr native-positions.xyz
COLVAR字段为time d a t c r.bias，本次保存0至0.2 ps共101帧。run.trr用于读取完整体系轨迹，native-positions.xyz保存这几个CV所需原子的坐标。
二面角绘图做了连续展开，便于观察跨越±π的变化；原始COLVAR仍保留PLUMED的有符号角度，不把展开后的数值当成另一种角度定义。距离没有在100步内固定到0.5 nm，这符合有限强度谐偏置的含义，不能把RESTRAINT当作严格几何约束。
下一步按离线driver教程重算同一批帧；需要理解偏置力时看谐势与原子力。已有GROMACS联用、LAMMPS联用和GPUMD联用也保留各自输入与环境要求。