GPUMD 联用 PLUMED：距离与谐偏置

本例先选两个 Si 原子之间的距离作为集体变量，再对这个距离施加谐偏置。按准备结构、提交任务、查看距离、检查偏置的顺序操作，后续接MetaD 与 WTMetaD及重加权。

下载完整输入、实际轨迹、COLVAR、HILLS 与绘图程序。资料包含 P100 和 V100 16G 的实际结果；每组动态任务 8000 步、步长 0.25 fs，共 2 ps，用于演示操作。

1 先看结构和编号#

晶格常数为 5.43 Å，金刚石常规胞沿三个方向各复制两次，得到 64 原子、边长 10.86 Å 的周期晶胞。原始模型文件以如下内容开头，完整 64 行坐标在资料包中：

64
Lattice="10.86 0 0 0 10.86 0 0 0 10.86" Properties=species:S:1:pos:R:3 pbc="T T T"
Si 0 0 0
Si 0 2.715 2.715
Si 2.715 0 2.715
Si 2.715 2.715 0
Si 1.3575 1.3575 1.3575

PLUMED 的 ATOMS=1,5 从 1 开始编号，对应第 1、5 行原子坐标。初始距离为 2.351259 Å。DISTANCE 使用周期边界下的距离，不能只取坐标直接相减而忽略跨边界情况。

势函数是 GPUMD 官方公开的 Si_2022_NEP4_4body.txt，在算例中命名为 nep.txt。它用于本例 Si 体系；换元素、相态或研究范围时，需要重新选择并检查适用的势函数。

2 先复现配套结果，再提交计算#

把下载文件保存为 gpumd-plumed.zip，进入新的工作目录：

unzip gpumd-plumed.zip -d gpumd-plumed
cd gpumd-plumed
bash view.sh
ls analysis/figures

view.sh 调用 ASE、PLUMED 的 driver 和 sum_hills，处理已提供的真实结果，不提交新的 MD。新计算按下面的步骤进行：

python3 prepare.py
cd cases/p100
sbatch job.slurm
squeue -u "$USER"

V100 16G 使用 cases/v100/job.slurm。这两个脚本各申请 1 张 GPU、1 个 CPU，先跑两组固定构型对照，再跑无偏置、谐偏置、普通 MetaD、WTMetaD 四种短采样。脚本中的分区应与自己的资源权限一致。

当前 GPUMD 5.9.1 构建使用 GCC 13.3、CUDA 12.9；联用的是 GPUMD 私有的 PLUMED 2.10.1。先加载相应 GPUMD 模块，后处理工具使用：

source /etc/profile
module load cuda/12.9 gpumd/5.9.1-p100
PLUMED="$GPUMD_PLUMED_ROOT/bin/plumed"
"$PLUMED" info --version

不要在这个环境中另外加载公共默认 PLUMED 2.9.4，以免工具、内核和输入语法混用。

3 距离和谐偏置分别怎么写#

只记录距离时，plumed.dat 为：

UNITS LENGTH=A TIME=fs ENERGY=eV
d: DISTANCE ATOMS=1,5
PRINT ARG=d FILE=COLVAR STRIDE=10 FMT=%24.16e

增加谐偏置后，完整文件为：

UNITS LENGTH=A TIME=fs ENERGY=eV
d: DISTANCE ATOMS=1,5
r: RESTRAINT ARG=d AT=2.40 KAPPA=20
PRINT ARG=d,r.bias FILE=COLVAR STRIDE=10 FMT=%24.16e

d 是距离变量，r.bias 是谐偏置能。这里使用 Å、fs、eV，力常数 KAPPA=20 的单位为 eV/Ų。

V(d) = 1/2 × 20 × (d - 2.40)^2

距离偏离 2.40 Å 越多，偏置越强。RESTRAINT 是有限强度的谐势，不会把距离严格固定在 2.40 Å；热运动与原有 Si 势函数仍共同决定实际分布。

对应的 run.in 为：

potential nep.txt
velocity 300 seed 13579
ensemble nvt_bdp 300 300 100
time_step 0.25
plumed plumed.dat 1 0
dump_thermo 10
dump_xyz 10 trajectory.xyz velocity force precision double
run 8000

plumed plumed.dat 1 0 表示读取该文件、每一步调用、从新计算开始。施加偏置时使用每步调用的 1，让偏置力进入积分过程。最后的重启参数涉及 PLUMED 文件的续读，不能据此认为坐标、速度和恒温器的全部状态都已恢复。

4 想修改自己的偏置时，单独准备干净目录#

回到资料包根目录，只复制输入，避免沿用已有输出：

mkdir my-umbrella
cp reference-results/p100/umbrella/{model.xyz,nep.txt,run.in,plumed.dat} my-umbrella/
cd my-umbrella
nano plumed.dat

例如修改 AT 为自己的目标距离，再新建 job.slurm：

#!/bin/bash
#SBATCH --job-name=Si-umbrella
#SBATCH --partition=gpu-8p100-hs
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=1
#SBATCH --gres=gpu:p100-16g:1
#SBATCH --time=00:20:00
#SBATCH --output=md-%j.out
set -eo pipefail
source /etc/profile
module purge
module load cuda/12.9 gpumd/5.9.1-p100
export OMP_NUM_THREADS=1
cd "$SLURM_SUBMIT_DIR"
gpumd

sbatch job.slurm

使用 V100 16G 时，分区、GPU 类型和 GPUMD 模块一起替换为对应入口。资料包的 run_native.py 会按教程参数重新生成四种输入；单独修改偏置时，用本节独立目录更直接。

5 任务结束后先看哪些文件#

head COLVAR
head thermo.out
head trajectory.xyz

COLVAR 的表头为 time d r.bias，第一行是初始构型，此后每 10 步输出一次，即 2.5 fs 一行；轨迹从第一次输出的 2.5 fs 开始。比较两者时应按时间匹配，不能直接把第一行与第一帧对应。

当前 GPUMD 接口的热力学势能和 XYZ energy 是物理势能，PLUMED 偏置能另见 r.bias。需要总势能时，要先匹配同一时刻，再考虑相加；不要把偏置前后的 thermo.out 势能差直接称为偏置能。

本例初始距离为 2.351259 Å，偏置能为 0.0237569 eV。固定构型的两次单点使用同一组坐标：无偏置力与有偏置力相减，再与谐势解析求导比较，最大分量差约 2.36×10⁻⁷ eV/Å。这些固定构型检查使用 time_step 0，不代表额外的 MD 采样。

新任务结束后，在资料包根目录处理自己的四种轨迹：

source /etc/profile
module load ase/3.29.0 gpumd/5.9.1-p100
python3 analyze.py --cases cases --case p100 --output analysis-new

分析脚本针对本资料包的固定输入和输出组织；单独修改过长度、输出频率或目录后，应相应调整脚本的读取范围。