CP2K伞形采样与PLUMED接口

想让原子间距离围绕某个值波动，可以施加谐和偏置并设置多个窗口。这里用非周期Ar₂小模型演示：先定义距离，再分别运行原生约束和PLUMED版本，最后检查窗口分布。它适合熟悉操作流程，迁移到真实反应或吸附体系时需要重新选择CV、势能方法、窗口和采样时长。

1 取得输入并运行#

下载约束、伞形与重加权操作包，保存为cp2k-umbrella.zip。解压目录提供job.slurm，普通作业可用sbatch job.slurm提交。下面短命令用于已经分配至少4个物理CPU核心的计算环境；四个独立窗口各用1个MPI进程、1线程，后处理的四副本WHAM在窗口全部结束后运行。

python3 -m zipfile -e cp2k-umbrella.zip .
cd cp2k-umbrella
bash setup.sh

bash run.sh
bash verify.sh
bash plot.sh

run.sh依次完成八段伞形轨迹和一段三原子约束轨迹，再从已有轨迹提取高精度CV、偏置、权重并调用WHAM。完整结果放在results/，三张图片在plots/。已完成且输入相同的目录会复用结果；修改物理参数后请另建工作目录。

环境为CP2K2026.2 CPU模块cp2k/2026.2-8173m，Intel ifx2024.2.1、IntelMPI2021.13、MKL2024.2。MD接入的是该CP2K编译配套的PLUMED2.10.1；离线分析在独立进程加载plumed/2.9.4-oneapi2024.2。脚本已分开这两个环境，运行MD时请直接使用CP2K模块。

2 先看窗口与模型#

四个窗口中心为3.7、3.8、3.9、4.0 Å，偏置写为Vbias = ½ κ (r−r₀)²，κ=0.05 eV/Ų。每个窗口6000步、步长0.5 fs，共3 ps；Langevin温度10 K、摩擦系数0.01 fs⁻¹。每5步输出一次位置和力，包含初态共1201帧。前0.5 ps不参加后面的WHAM统计。

相互作用为CP2K FIST/Lennard–Jones：ε=119.8 K、σ=3.405 Å、截断7.5 Å、截断处平移能量，30 Å非周期盒子。四个窗口的起始距离各等于窗口中心。每对原生/PLUMED输入使用相同随机种子，便于比较接口；这两条轨迹不当作两份独立统计样本。

3 原生谐和窗口的完整输入#

下面是3.9 Å窗口inputs/native-r3.9/input.inp。COLVAR按原子编号1、2定义距离，COLVAR 1引用第一项CV。

&GLOBAL
 PROJECT native-r3.9
 RUN_TYPE MD
 PRINT_LEVEL MEDIUM
 SEED 20261111
&END GLOBAL
&MOTION
 &MD
  ENSEMBLE LANGEVIN
  STEPS 6000
  TIMESTEP [fs] 0.5
  TEMPERATURE 10
  &LANGEVIN
   GAMMA [fs^-1] 0.01
  &END LANGEVIN
  &PRINT
   &ENERGY
    FILENAME =energy.dat
    &EACH
     MD 1
    &END EACH
   &END ENERGY
  &END PRINT
 &END MD
 &CONSTRAINT
  &COLLECTIVE
   INTERMOLECULAR TRUE
   COLVAR 1
   TARGET 7.369931918254007
   &RESTRAINT
    K 0.0002572714783335379
   &END RESTRAINT
  &END COLLECTIVE
 &END CONSTRAINT
 &PRINT
  &TRAJECTORY
   FORMAT XYZ
   FILENAME =trajectory.xyz
   &EACH
    MD 5
   &END EACH
  &END TRAJECTORY
  &FORCES
   FORMAT XYZ
   UNIT hartree/bohr
   FILENAME =forces.xyz
   &EACH
    MD 5
   &END EACH
  &END FORCES
  &RESTART
   BACKUP_COPIES 0
   &EACH
    MD 1000
   &END EACH
  &END RESTART
 &END PRINT
&END MOTION
&FORCE_EVAL
 METHOD FIST
 &MM
  &FORCEFIELD
   &SPLINE
    EPS_SPLINE 1E-10
   &END SPLINE
   DO_ELECTROSTATICS .FALSE.
   SHIFT_CUTOFF .TRUE.
   &CHARGE
    ATOM Ar
    CHARGE 0.0
   &END CHARGE
   &NONBONDED
    &LENNARD-JONES
     ATOMS Ar Ar
     EPSILON [K_e] 119.8
     SIGMA [angstrom] 3.405
     RCUT [angstrom] 7.5
    &END LENNARD-JONES
   &END NONBONDED
  &END FORCEFIELD
  &POISSON
   &EWALD
    EWALD_TYPE NONE
   &END EWALD
  &END POISSON
 &END MM
 &SUBSYS
 &COLVAR
  &DISTANCE
   ATOMS 1 2
  &END DISTANCE
 &END COLVAR
  &CELL
   ABC 30 30 30
   PERIODIC NONE
  &END CELL
  &COORD
   Ar 10 10 10
   Ar 13.900000000000 10 10
  &END COORD
  &TOPOLOGY
   CONN_FILE_FORMAT OFF
  &END TOPOLOGY
 &END SUBSYS
 &PRINT
  &FORCES ON
  &END FORCES
 &END PRINT
&END FORCE_EVAL

原生TARGET的默认距离单位为Bohr，3.9 Å在本版本中为7.369931918254007 Bohr。原生RESTRAINT K采用K(s−s₀)²，没有前面的½。因此这里K=0.0002572714783335379 Ha/Bohr²，对应PLUMED的KAPPA=0.05 eV/Ų；不要只把数值0.05原样抄过去。

4 改用PLUMED接口#

配套的inputs/plumed-r3.9/input.inp保留相同势场、温度、步数和起始结构，将MOTION/CONSTRAINT中的谐和项替换为：

&FREE_ENERGY
 METHOD METADYN
 &METADYN
  USE_PLUMED TRUE
  PLUMED_INPUT_FILE plumed.dat
 &END METADYN
&END FREE_ENERGY

同一目录的plumed.dat为：

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

UNITS明确采用Å、eV、fs，RESTRAINT按½κΔr²施加偏置。FMT保留足够有效位，便于读取较小偏置；STRIDE=5与位置输出频率一致。更换CV时，需要同时修改距离定义、偏置参数和后处理脚本。

5 看实际轨迹与窗口重叠#

去掉前0.5 ps后，原生窗口平均距离依次为3.796311、3.841601、3.891010、3.947582 Å，标准差约0.0772–0.1010 Å。平均值不必等于窗口中心，因为体系本身的相互作用也施加力。图中相邻窗口分布有重叠，后续才能联合去偏置。

用终端查看四个窗口的实际输出：

source env.sh
for d in results/cases/plumed-*; do
  echo "$d"
  head -n 5 "$d/COLVAR.plumed"
done

每个目录还包含energy.dat、trajectory.xyz、forces.xyz、output.log和检查点。**本次CP2K版本的外部PLUMED接口在.ener中记录体系本身的势能，偏置能单独见COLVAR.plumed的r.bias；原生谐和约束的势能列则包含约束偏置。**比较总势能时，应统一成相同口径，并先换算Hartree和eV。

6 接下来怎么做#

先阅读静态重加权与WHAM组合窗口。若希望严格固定某个距离，阅读硬距离约束；若希望随时间添加高斯偏置，阅读原生MetaD。真实体系还需增加平衡和生产时间、查看自相关、多次独立采样和分块稳定性；本例3 ps主要展示命令和文件之间的衔接。

参考：CP2K约束输入、PLUMED RESTRAINT、本版本PLUMED接口实现。