PLUMED 二肽集体变量实操

这一篇从一个小分子模型开始，把 PLUMED 的集体变量计算接到实际动力学中：准备结构与力场，写 plumed.dat，提交任务，得到 COLVAR，最后画出距离、角度和分子构象图。

下载完整结构、拓扑、输入、四组真实轨迹与分析程序。模型是常用的丙氨酸二肽 ACE–ALA–NME，共22个原子。结构与拓扑取自 GROMACS 官方测试数据，使用 AMBER99SB-ILDN；本例在真空模型中运行，未加入水或离子。P100、V100 16G 各计算两组不同初速度种子的短轨迹，每组5 ps，用来演示操作，不据此给出生物分子的平衡构象占比或自由能面。

1 先复现参考结果#

把下载文件保存为 plumed-alanine.zip，解压后运行：

unzip plumed-alanine.zip
cd plumed-alanine-wave280
bash run.sh
ls analysis
ls analysis/figures

run.sh 读取附件已有的实际数据，不启动新动力学。它会生成 CSV 和三张 PNG：集体变量随时间的变化、周期扭转角、以及带编号的真实分子结构。

head analysis/summary.csv
head analysis/p100-seed41-CV.csv

表中距离单位是 nm，角度原始单位是 rad，时间单位是 ps。图中的键角和扭转角有些转换成了度，读图时看清坐标轴。

2 找准原子编号#

PLUMED 的原子编号从1开始。本例沿用 initial.gro 的22个原子顺序，不另外重排：

head -n 26 inputs/initial.gro
head -n 42 inputs/topol.top

topol.top 的 [ atoms ] 列表给出编号、原子名、电荷和质量。换成自己的分子后，先从结构与拓扑确认编号，再修改这些选择，不能只把体系文件换掉而保留旧编号。

3 写 PLUMED 输入#

附件中的完整 inputs/plumed.dat：

UNITS LENGTH=nm ENERGY=kj/mol TIME=ps
WHOLEMOLECULES ENTITY0=1-22
d: DISTANCE ATOMS=1,19
a: ANGLE ATOMS=5,7,9
phi: TORSION ATOMS=5,7,9,15
psi: TORSION ATOMS=7,9,15,17
c: COORDINATION GROUPA=6,16 GROUPB=7,17 R_0=0.35 NN=6 MM=12
PRINT ARG=d,a,phi,psi,c STRIDE=10 FILE=COLVAR FMT=%24.16e
DUMPATOMS ATOMS=1-22 STRIDE=10 FILE=plumed-coordinates.xyz PRECISION=16 UNITS=nm

WHOLEMOLECULES 按编号把分子恢复成连续构型，便于处理跨周期边界的分子；距离、键角和扭转角仍按当前边界处理。PRINT 每10个动力学步保存一行，同时保存对应时刻的高精度坐标，方便把数据与结构对应起来。

这里只观察集体变量，没有施加偏置。要在跑通本例后加入谐偏置或元动力学，再接着看谐偏置与原子力、伞形采样及WTMetaD重加权。

4 准备自己的新任务#

bash prepare.sh

按提示输入 1 选择 P100，或输入 2 选择 V100 16G，直接回车默认 P100。程序创建新目录并给出下一步命令，例如：

cd work/p100-seed41
sbatch run.sh
squeue -u "$USER"

准备程序不会覆盖已有的同名工作目录。套餐用户按自己的分区权限调整提交脚本，保留本例的单GPU、单MPI和配套环境。CPU或GPU环境的区别见四种引擎接口。

本例采用 gromacs/2026.4-p100-mpi-plumed 或对应 V100 模块，配套 PLUMED 2.9.4 的 GNU/OpenMPI 5.0.10 构建。加载 GROMACS 模块后使用它带入的 PLUMED 内核。

任务脚本先做200步陡降预处理，然后执行5 ps NVT。预处理用于整理起始构型，本例不将其称为已经按力阈值收敛的结构优化。正式研究时应按照自己的体系检查结构准备、平衡过程和采样长度。

5 动力学输入和实际命令#

md.mdp 主要内容如下；附件保存完整文件：

integrator = md
dt = 0.001
nsteps = 5000
cutoff-scheme = Verlet
nstlist = 10
rlist = 1.4
coulombtype = Cut-off
rcoulomb = 1.4
vdwtype = Cut-off
rvdw = 1.4
pbc = xyz
constraints = h-bonds
constraint-algorithm = lincs
lincs-order = 4
continuation = no
gen-vel = yes
gen-temp = 300
gen-seed = 41
tcoupl = v-rescale
tc-grps = System
tau-t = 0.1
ref-t = 300
pcoupl = no
comm-mode = Linear
nstcomm = 100
nstcalcenergy = 10
nstenergy = 10
nstlog = 100
nstxout = 10
nstvout = 0
nstfout = 0
nstxout-compressed = 0

3 nm周期盒、1.4 nm截断和真空分子是本例的模型设定，不是水溶液蛋白体系的通用参数。温度设为300 K，22原子小体系的瞬时温度会有明显涨落。

提交脚本内实际执行：

gmx_mpi grompp -f em.mdp -c initial.gro -p topol.top -o em.tpr
mpirun -np 1 gmx_mpi mdrun -deffnm em -ntomp 1 -nb gpu -pme cpu -update cpu
gmx_mpi grompp -f md.mdp -c em.gro -p topol.top -o md.tpr
mpirun -np 1 gmx_mpi mdrun -deffnm md -ntomp 1 -nb gpu -pme cpu -update cpu -plumed plumed.dat

最后一行的 -plumed plumed.dat 把集体变量文件接入动力学。上面每一步都使用同一套结构、拓扑和模块；grompp 报出的结构或拓扑错误应先解决，再进入下一步。

6 看 COLVAR 与实际变化#

head -n 10 COLVAR
tail COLVAR
ls -lh plumed-coordinates.xyz

P100种子41的开头数据（略去额外周期范围注释，数字做适当显示舍入）：

#! FIELDS time d a phi psi c
0.000 0.72506033 2.18917271 -2.81571730 2.88927334 2.77265824
0.010 0.72625384 2.22915929 -2.78645845 2.83315187 2.75189979
0.020 0.72787454 2.20850536 -2.79672828 2.85061883 2.71222692

每条轨迹含0 ps初始帧，共501行。4组实际结果中，端点距离变化范围如下：

c 不是整数，也不是化学键价或溶剂配位数。这里有四对 O–N 接触，每一对按距离给出一个平滑权重，再把权重相加。R_0=0.35 nm 控制接触距离尺度；改变它会改变计数含义。配套脚本按当前PLUMED版本的实际有理切换函数处理自动截断与归一化。

7 正确处理周期扭转角#

φ和ψ原始输出在 [-π, π] 内。画时间曲线时可以转换为度，也可以展开角度；画φ–ψ分布时仍应保留周期变量的含义。

例如，在自己的 Python 分析中：

import numpy as np

values = np.loadtxt("COLVAR")
phi_rad = values[:, 3]
psi_rad = values[:, 4]
phi_deg = np.rad2deg(phi_rad)
psi_deg = np.rad2deg(psi_rad)
phi_trace_deg = np.rad2deg(np.unwrap(phi_rad))

比较两个角度时，用圆周差 atan2(sin(a-b), cos(a-b))，而不是直接用 a-b。这样不会把靠近 -π 和 π 的两个相近角度当作相差整整一圈。

本例用实际坐标重新计算距离、键角、扭转角和平滑接触，与原生输出逐帧对应。结构图、曲线和CSV均来自这些实际保存帧。

8 分析自己的新轨迹#

回到资料包根目录，给新结果指定一个新的输出目录：

module load ase/3.29.0
python analyze.py work analysis-my-run

work 下保留 prepare.sh 创建的一个或两个算例目录即可；每个目录需要 COLVAR、plumed-coordinates.xyz 和 energy.xvg。任务脚本已经在结束时从EDR提取温度和势能。

head analysis-my-run/summary.csv
ls analysis-my-run/figures

脚本基于本例22原子的编号和四对接触编写。换分子、编号、输出列或盒形后，先相应修改选择与几何重建部分，再分析自己的结果。

结构与拓扑出处：GROMACS官方测试数据。附件保留原始拓扑及其LGPL 2.1许可说明。更多输入组合见 PLUMED专题。