CP2K静态重加权与WHAM

伞形采样给体系施加了偏置，轨迹里的距离分布因此改变。得到无偏分布时，需要记录每个窗口的CV、偏置和温度，再正确组合。本篇接着四个Ar₂窗口，展示从轨迹到PLUMED重加权、原生WHAM和图像的完整流程。

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模块。

如果前两篇已在同一目录完成，直接查看这里的结果即可。四个PLUMED窗口每个3 ps，去除前0.5 ps后各保留1000帧，总计4000帧；原生CP2K的四条对应轨迹用于接口对照，不重复加入WHAM。

2 用已有轨迹提取CV、偏置与静态权重#

以3.9 Å窗口为例，离线控制文件包含：

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

原生工具的调用形式是plumed driver --plumed driver.dat --ixyz trajectory.xyz --length-units A --timestep 0.0005 --trajectory-stride 5，在对应的分析目录中执行。driver的--timestep采用ps，本例0.0005 ps=0.5 fs；--trajectory-stride 5表示XYZ每5个MD步一帧。输入中UNITS TIME=fs控制打印时间，最终相邻输出应间隔2.5 fs。操作包已逐帧核对实际时间，避免把输出步号当成物理时间。

REWEIGHT_BIAS输出的是未归一化权重的对数Vbias/kBT，后续要归一化再用于统计。高偏置区域权重会较大，并不意味着采样充分。这里是固定谐和偏置；随时间变化的MetaD需要相应的时间依赖偏置修正，不能直接照抄这套静态公式。

3 联合四个窗口运行原生WHAM#

在results/driver/wham/full/中可以看到脚本生成的完整控制文件：

UNITS LENGTH=A ENERGY=eV TIME=fs
x: READ FILE=all.COLVAR VALUES=x IGNORE_TIME IGNORE_FORCES
b: RESTRAINT ARG=x AT=@replicas:{3.7,3.8,3.9,4.0} KAPPA=0.05
h: WHAM_HISTOGRAM ARG=x BIAS=b.bias TEMP=10 GRID_MIN=3.3 GRID_MAX=4.8 GRID_BIN=150
f: CONVERT_TO_FES GRID=h TEMP=10
DUMPGRID GRID=h FILE=histogram.dat FMT=%.16e
DUMPGRID GRID=f FILE=fes.dat FMT=%.16e
WHAM_WEIGHTS BIAS=b.bias TEMP=10 FILE=weights.dat FMT=%.16e

对应原生调用为mpirun -np 4 plumed driver --noatoms --plumed wham.dat --multi 4 --timestep 1。这里的4个副本按不同中心评价同一批组合样本的偏置；它们是后处理进程，与前面四条MD轨迹分阶段运行。all.COLVAR第一列是后处理索引，真实来源窗口和物理时间保留在samples.csv中，不能把拼接后的索引称为连续物理轨迹。

该版本WHAM收集阶段会跳过分析步0，脚本因此在文件开头保留一行初始化占位，再写4000个真实样本。占位行不计入4000帧；请连同控制文件一起保留脚本生成的数据，不要随意删第一行。这里只去掉伞形偏置，体系的LJ相互作用仍保留。

4 读结果、画图与检查统计#

source env.sh
"$YS_PY" - <<'PYCODE'
import json
from pathlib import Path
root = Path('results/driver/wham')
a = json.loads((root/'WHAM-VALIDATION.json').read_text())
for row in a['results']:
    print(row['block'], 'frames=', row['samples'],
          'neighbor overlap=', row['neighbor_overlap'])
print('FES:', root/'full/fes.0.dat')
print('weights:', root/'full/weights.0.dat')
PYCODE

fes.0.dat第一列为距离Å，第二列为自由能eV；自由能可整体平移零点。未访问网格保持未知/非有限值，画图时不能补成零自由能。weights.0.dat保存逐帧WHAM权重；不能用单窗口exp(Vbias/kBT)直接代替多窗口求出的联合权重。

本次4000帧样本相邻窗口重叠指标约0.833、0.833、0.822，访问到57个网格。这个重叠指标只描述直方图覆盖。图中的前后两半各2000帧，曲线仍有差别，表明短轨迹的统计起伏不能忽略。

对自由的三维二原子距离，径向分布含r²几何因子，因此已知模型的参考径向自由能为U(r)−2kBT ln r加任意常数；图中的虚线采用这一表达式，而不是只画LJ势能。短时间数据与参考形状有偏离，不能据此报告高精度自由能势垒。

5 换成自己的生产计算#

先选具有物理意义的CV和窗口范围，再用短预跑调整窗口间距、弹簧强度和采样时长，使邻窗有效重叠。正式计算应延长平衡与生产段、比较独立种子、检查自相关并增加时间块；必要时检查窗口内是否存在其他慢自由度。不同温度、不同CV定义或不同偏置参数的文件不能未经处理直接混合。

相关内容：伞形与PLUMED接口、固定距离约束、MetaD续算与偏置历史、PLUMED伞形采样、PLUMED WHAM。

参考：PLUMED driver、REWEIGHT_BIAS、WHAM_HISTOGRAM。