PLUMED偏置力链：从集体变量到原子力

集体变量数值正确后，还应确认偏置能量对原子坐标的导数。对自定义函数、虚拟原子或多层组合，这一步能帮助发现变量、单位及周期性设置的问题。本篇接着几何与函数教程，使用同一个资料包和同一次driver结果。

1. 查看已经计算的偏置与力#

下载四原子几何与偏置资料包，保存为example.zip：

python3 -m zipfile -e example.zip .
cd plumed-geometry-bias-chain
bash results.sh
bash read.sh
bash plot.sh

先用保存的真实结果查看图表，不需要重新计算。figures/是图片，data/是CSV，reference-results/是原生结果。

head reference-results/COLVAR
head -12 reference-results/forces.xyz
head reference-results/derivative-check.dat

COLVAR的b.bias是总谐偏置能量，单位kJ/mol；forces.xyz保存每个真实原子的三个偏置力分量，单位kJ/(mol·nm)。这里没有材料势或DFT，输出力只来自本输入的偏置。

forces.xyz每帧仍有XYZ样式的两行头部；第二行含有维里相关值，不能当成总能量读取。能量请从COLVAR的b.bias列取。

2. 谐偏置如何回传到原子#

输入末尾：

b: RESTRAINT ARG=r,a,t,rg,n,cd,fd,custom,combo,less AT=0.12,1.4,1.2,0.15,1.5,0.1,0.15,0.02,0.1,0.7 KAPPA=3,2,1,4,1.5,2,1,2,1,1

每个参数按同一位置一一对应。偏置为V=Σ K_i(s_i−s_i,0)²/2，再沿集体变量链求导，得到F=−∂V/∂R。二面角差值要用其周期性处理，不能简单跨越−π/π相减。

KAPPA的单位随变量变化：距离项使用kJ/(mol·nm²)，角度项使用kJ/(mol·rad²)；对于自己定义的组合变量，应按该变量的定义确定常数单位。这里的常数用于导数演示，不是可直接用于任意材料的推荐参数。

本例包含相对固定原点的距离fd，因此外部锚点可以对体系施加净力。不能据此要求总偏置力恒为零。

3. 原生导数检查#

run.sh中的：

"$plumed_prefix/bin/plumed" driver   --plumed plumed.dat --ixyz scan.xyz   --length-units nm --box 1,1,1   --dump-forces forces.xyz --dump-forces-fmt %.12e   --debug-forces derivative-check.dat

会同时保存原子力和PLUMED原生数值导数对照。按上一教程使用sbatch run.slurm运行即可。--debug-forces会额外扰动坐标，适合小体系排查；通常不需要加在正式长轨迹的每次后处理中。

4. 独立能量差分#

bash read.sh
head data/PLUMED-atom-force-check.csv
cat analysis-read/DATA-VALIDATION.json

读取程序使用独立的距离、角度、二面角、几何中心、开关和谐偏置公式，以h=10⁻⁶ nm对所有原子坐标做中心差分：F≈−[V(R+h)−V(R−h)]/(2h)，不直接复用PLUMED输出力作为参考。

两种差分过程不同，误差不必完全一致。这里说明这套具体变量链在给定坐标和参数下的数值一致性；不代表材料势精度、正式MD时间步稳定性或自由能已收敛。

自己的结果可以用：

bash read.sh work analysis-new
bash plot.sh analysis-new figures-new work

若误差明显增大，先查看坐标单位、原子编号、周期性、开关参数及角度跨周期情况，再选择少量帧做更细差分。derivative-check.dat还可能含盒子相关分量，本例只把前12个笛卡尔坐标分量逐帧与原子力比较。