CP2K硬距离约束与SHAKE

严格固定键长或原子间距离时，可以使用CP2K原生集体变量约束。它在动力学积分中投影位置与速度；谐和伞形窗口则允许距离围绕中心波动。本篇用三个Ar原子的实际轨迹，演示输入、运行和结果检查。

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.9 Å的近等边三角形，固定原子1和2的距离，原子3仍能运动：

Ar 10.000000000000 10.000000000000 10.000000000000
Ar 13.900000000000 10.000000000000 10.000000000000
Ar 11.950000000000 13.377499074759 10.000000000000

这里选择三原子模型，使固定一个距离后仍有内部运动自由度。原子编号始终对应输入结构顺序。迁移到分子或界面时，应先确定要约束的物理过程与自由度，再选择原子和CV。

3 完整输入与关键参数#

下载包中实际使用的inputs/hard-Ar3/input.inp如下：

&GLOBAL
 PROJECT hard-distance-Ar3
 RUN_TYPE MD
 PRINT_LEVEL MEDIUM
 SEED 20261111
&END GLOBAL
&MOTION
 &MD
  ENSEMBLE LANGEVIN
  STEPS 1000
  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
  SHAKE_TOLERANCE 1e-8
  CONSTRAINT_INIT TRUE
  &COLLECTIVE
   INTERMOLECULAR TRUE
   COLVAR 1
   TARGET 7.369931918254007
  &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
  &VELOCITIES
   FORMAT XYZ
   UNIT angstrom/fs
   FILENAME =velocities.xyz
   &EACH
    MD 5
   &END EACH
  &END VELOCITIES
  &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
   Ar 11.950000000000 13.377499074759 10
  &END COORD
  &TOPOLOGY
   CONN_FILE_FORMAT OFF
  &END TOPOLOGY
 &END SUBSYS
 &PRINT
  &FORCES ON
  &END FORCES
 &END PRINT
&END FORCE_EVAL

COLVAR 1引用SUBSYS/COLVAR中的距离。TARGET采用该版本默认的Bohr单位，7.369931918254007 Bohr对应3.9 Å。这个硬约束块没有RESTRAINT子块；添加谐和项会改变所执行的物理操作。

CONSTRAINT_INIT TRUE在开始时处理约束初态，SHAKE_TOLERANCE 1e-8设置迭代精度。本例使用Langevin、10 K、0.5 fs步长、1000步；在此版本输入中使用SHAKE_TOLERANCE。位置、速度、力每5步输出，速度显式采用Å/fs，力显式采用Ha/Bohr。

4 查看实际距离与速度#

下面命令直接读取已计算的字段，输出固定距离误差与径向相对速度的最大绝对值：

source env.sh
"$YS_PY" - <<'PYCODE'
import numpy as np
x = np.load('results/hard-Ar3/fields.npz')
print('r12 max error / Angstrom:', abs(x['distance12_A']-3.9).max())
print('radial relative velocity / Angstrom fs^-1:',
      abs(x['radial_relative_velocity_A_fs']).max())
PYCODE

本次计算固定距离的最大残差约9.60×10⁻¹¹ Å，径向相对速度最大约2.16×10⁻⁷ Å/fs。位置约束和速度约束都需要检查；只看某一帧距离接近目标不能说明整段轨迹都符合设置。

图中r₁₂保持3.9 Å，r₁₃和r₂₃仍有变化，体现了约束选定距离后其他自由度继续运动。结构图使用实际轨迹的x-y投影，并保持两个坐标方向相同的长度比例。

5 迁移到自己的体系#

先替换结构、原子种类和势能设置；再修改DISTANCE/ATOMS与TARGET，检查约束是否和初态结构一致。保留速度与位置输出，先跑短任务确认无异常约束漂移，再延长生产段。使用DFT时，还需要独立设置基组、赝势、电子收敛和截断等参数，本篇的经典Ar势场不适用于材料DFT。

硬约束得到的是受约束轨迹。若要进一步求平均约束力或自由能导数，还涉及约束力输出、统计采样和相应几何修正，不能仅由这里的固定距离直接得到势垒。窗口联合去偏置可参考伞形采样与WHAM教程。

参考：CP2K CONSTRAINT、COLLECTIVE。