PLUMED专题

CP2K 联用 PLUMED：从输入到偏置力

使用两个 Ar 原子的经典模型，准备 CP2K 与 PLUMED 输入、提交六组短轨迹，提取距离、偏置能与真实受力。

更新于 2026-11-07约 7 分钟阅读按任务与环境选择适用配置

阅读本篇目录

本例通过一个很小的 Ar 二聚体说明：先准备 CP2K 结构和动力学输入，再让 PLUMED 跟踪原子间距并施加谐偏置，最后查看 COLVAR、轨迹和力。它是理解接口、单位和文件流程的例子，不用于判断实际材料的热力学性质。

下载完整输入、六组真实参考轨迹与绘图程序。使用 CP2K 2026.2 CPU FIST，oneAPI 2024.2 / Intel MPI 2021.13，以及其配套 PLUMED 2.10.1；没有电子结构自洽计算。

六组实际短轨迹的距离变化和谐偏置能

1 准备工作目录#

unzip plumed-four-engines.zip
cd plumed-four-engines-wave278
ls inputs/cp2k

三个初始距离分别为 3.40、3.82、4.50 Å，每个距离各有无偏置对照和谐偏置两组。zero 组仍启用 PLUMED，但 KAPPA 为 0，以便用相同接口记录距离。

设置	本例取值
原子、质量	2 个 Ar，39.948 u
晶胞	30 × 30 × 30 Å，三维周期
势函数	Lennard-Jones，ε/kB = 119.8 K，σ = 3.405 Å
截断	7.5 Å，截断处平移能量
动力学	NVE，20 步，步长 0.1 fs，总长 2 fs
初始速度	两原子都显式设为零
谐偏置	中心 4 Å，KAPPA = 2.5 eV/Å²；对照组为 0

30 Å 周期盒与 7.5 Å 截断使本例不与周期镜像在截断范围内相互作用。这里显式指定零速度；输入中的温度不是一次随机抽样，也没有用恒温器把体系维持在该温度。

2 完整 CP2K 输入#

下面是 inputs/cp2k/ar3p4-bias/input.inp。其余五组只改变原子初始间距、项目名和 PLUMED 的 KAPPA。

&GLOBAL
 PROJECT ar3p4-bias
 RUN_TYPE MD
 PRINT_LEVEL MEDIUM
&END GLOBAL
&MOTION
 &MD
  ENSEMBLE NVE
  STEPS 20
  TIMESTEP [fs] 0.1
  TEMPERATURE 1
  TEMP_TOL 0
 &END MD
 &FREE_ENERGY
  &METADYN
   USE_PLUMED .TRUE.
   PLUMED_INPUT_FILE plumed.dat
  &END METADYN
 &END FREE_ENERGY
 &PRINT
  &TRAJECTORY
   FORMAT XYZ
   &EACH
    MD 1
   &END EACH
  &END TRAJECTORY
  &FORCES
   FORMAT XYZ
   &EACH
    MD 1
   &END EACH
  &END FORCES
 &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
  &CELL
   ABC 30 30 30
   PERIODIC XYZ
  &END CELL
  &COORD
   Ar 10 10 10
   Ar 13.400000000000 10 10
  &END COORD
 &VELOCITY
  0 0 0
  0 0 0
 &END VELOCITY
 &KIND Ar
  ELEMENT Ar
  MASS 39.948
 &END KIND
  &TOPOLOGY
   CONN_FILE_FORMAT OFF
  &END TOPOLOGY
 &END SUBSYS
 &PRINT
  &FORCES ON
  &END FORCES
 &END PRINT
&END FORCE_EVAL

METHOD FIST 选择经典力场；DO_ELECTROSTATICS .FALSE. 和零电荷关闭本例不需要的静电项。PERIODIC XYZ 与晶胞共同确定边界。

MOTION/FREE_ENERGY/METADYN 中的两行负责接入 PLUMED；具体施加哪种偏置由 plumed.dat 决定。

3 PLUMED 输入与力的含义#

UNITS LENGTH=A ENERGY=eV TIME=fs
d: DISTANCE ATOMS=1,2
b: RESTRAINT ARG=d AT=4 KAPPA=2.5
PRINT ARG=d,b.bias STRIDE=1 FILE=COLVAR FMT=%24.16e
DUMPATOMS ATOMS=1,2 STRIDE=1 FILE=plumed-coordinates.xyz PRECISION=16 UNITS=A
DUMPFORCES ARG=d STRIDE=1 FILE=cv-force.dat FMT=%24.16e
DUMPDERIVATIVES ARG=d STRIDE=1 FILE=cv-derivatives.dat FMT=%24.16e

DISTANCE 读取第 1、2 个原子，编号从 1 开始。
RESTRAINT 使用 V = ½K(d-4)²，不是 K(d-4)²。
DUMPFORCES ARG=d 输出作用在距离这个 CV 上的广义力 -K(d-4)，不是直接输出每个原子的三维力。
DUMPDERIVATIVES 输出 CV 对坐标及晶胞的导数。取前六个原子笛卡尔坐标分量，与 CV 力相乘，得到两个原子的偏置力；后面的九个晶胞导数不能误当成另三个原子。

设单位向量从原子 1 指向 2，则原子 1 的偏置力是 K(d-4) r̂，原子 2 为相反方向。距离小于中心时拉开，距离大于中心时拉近。

4 提交新计算#

资料包提供单核 CPU 脚本，依次运行六个小算例：

sbatch cp2k-job.sh
squeue -u "$USER"

CPU 分区暂时不能调度时，可以先执行下一节的参考数据分析，待资源可用后再提交。任务脚本调用的实际计算命令是：

module load cp2k/2026.2-8173m
mpirun -np 1 cp2k.psmp -i input.inp -o cp2k.out

cp2k.sh 会把 inputs/cp2k 复制成新的 work-cp2k，保留参考文件。若 work-cp2k 已存在，脚本会停止并提醒换新目录，防止覆盖已有结果。

ls work-cp2k/ar3p4-bias
head work-cp2k/ar3p4-bias/COLVAR
tail work-cp2k/ar3p4-bias/cp2k.out

5 从真实结果看偏置作用#

初始距离 / Å	KAPPA / eV·Å⁻²	2 fs 距离 / Å	首次记录偏置 / eV
3.40	0	3.40007228	0
3.40	2.5	3.40144835	0.44999973
3.82	0	3.82000010	0
3.82	2.5	3.82041295	0.04050000
4.50	0	4.49999376	0
4.50	2.5	4.49884695	0.31249998

3.40 和 3.82 Å 的偏置组都向 4 Å 增大，4.50 Å 的偏置组向 4 Å 减小。2 fs 很短，图只显示响应方向和输入含义；不能据此认为体系已经平衡或已经充分采样。

ar3p4-bias/COLVAR 开头是：

#! FIELDS time d b.bias
0.100000  3.4000001807103786  0.4499997289344730
0.200000  3.4000079686562401  0.4499880470950142
0.300000  3.4000233637283372  0.4499649550898239

本例 CP2K 原生坐标含 0 fs 初始帧，共 21 帧；PLUMED 首次接收坐标在 0.1 fs，共 20 行。因此要把 COLVAR 与 CP2K 轨迹的第 1 步开始对齐。初始半步只使用初始化时的 LJ 力，不能把初始 0 fs 当成已经施加并输出了 PLUMED 偏置的帧。

谐偏置力示意及真实轨迹反推出的力与解析力对照

图中使用连续三个位置和速度 Verlet 的位置关系反推出原子力，再与同一时刻的 LJ 力加 PLUMED 偏置力比较。第一处有偏置的实际力及全部 19 个可重建时刻均一致，最大偏差约 6.0×10⁻⁵ eV/Å，主要受 CP2K 轨迹坐标打印精度影响。

不要把宿主输出中的某一列势能直接认作“已经含 PLUMED 的总势能”。本例将原生 LJ 部分与 COLVAR 中的 b.bias 分开读取，偏置是否作用通过实际轨迹的受力响应判断。

6 重画与分析自己的结果#

bash run.sh
head analysis/cp2k/cp2k-interface-summary.csv
ls figures

这一步不启动新动力学。三张 PNG 分别显示四种接口的原生数据、六组 Ar 轨迹、以及偏置力示意和数值对照。CSV 保存距离、偏置能、偏置力和数值误差，方便修改配色或整理自己的图。

新任务结束后可单独处理：

module load ase/3.29.0
python analyze_cp2k.py work-cp2k analysis/my-result

输出目录自己取一个新名字。换成自己的体系时，也要修改原子选择、势函数、参考中心、KAPPA、时间步长和采样长度；这个二原子 LJ 检查程序的解析力公式不能直接拿去判断 DFT、水或复杂材料的总力。

继续看四种动力学软件的接口区别、PLUMED 谐偏置到原子力和伞形采样。

CP2K 官方接口参数；PLUMED 官方 DUMPFORCES。

需要根据你的体系选择资源？联系讨论计算需求。实际可用软件与权限，以当前账户和运行环境为准。