PLUMED谐偏置：从RESTRAINT到原子力

本篇用水体系氧1到氧4的距离，说明RESTRAINT产生什么偏置，以及怎样检查力的方向和单位。资料包同时给出实际轨迹与独立几何扰动；后者用于求导检查，不是新的MD采样。

下载输入、偏置能、原子力与绘图程序。

1 完整复跑#

将下载文件保存为plumed-example.zip，在新的CPU工作目录执行：

unzip plumed-example.zip -d plumed-force
cd plumed-force
bash run.sh

默认入口不提交MD任务。偏置力结果位于：

head results/force-fixture/COLVAR
head results/force-fixture/forces.xyz
head results/force-fixture/reference-vs-FD.csv

force-fixture从实际MD起点取前7个原子，依次沿每个原子坐标做±1e-5 nm扰动，连同起点共43帧，用于独立有限差分。

2 先写清楚势能和单位#

UNITS LENGTH=nm TIME=ps ENERGY=kj/mol
d: DISTANCE ATOMS=1,4
r: RESTRAINT ARG=d AT=0.5 KAPPA=100
PRINT ARG=d,r.bias STRIDE=1 FILE=COLVAR FMT=%24.16e

本例的谐势为：

V(d) = 1/2 × k × (d - d0)^2
k = 100 kJ·mol^-1·nm^-2
d0 = 0.5 nm

起始距离约0.997882 nm，偏置能约12.394330 kJ/mol。当距离大于目标时，偏置力倾向于拉近这两个氧；小于目标时，方向相反。

这是有限强度的谐偏置，距离仍会随温度和分子力场波动。RESTRAINT不是严格固定几何距离，也不能仅凭一条短偏置轨迹得到收敛的自由能曲线。

3 用driver输出原子偏置力#

资料包入口实际执行的核心命令为：

plumed driver --plumed plumed.dat --ixyz frames.xyz --length-units nm --box 3.0067100524902344,3.0067100524902344,3.0067100524902344 --dump-forces forces.xyz

这里使用原生起点的周期晶胞。完整脚本还设置了输出有效数字。forces.xyz给出这个几何输入在PLUMED偏置下的原子力；driver没有执行TIP3P力场，因此它不是GROMACS总原子力。

不要把DUMPFORCES ARG=d输出的CV广义力直接当成某个原子的x、y、z力。CV广义力需要通过CV对原子坐标的导数传递到原子；DUMPFORCES官方说明专门区分了两者。

4 两种独立检查#

解析检查对最小镜像距离求梯度：

r = R1 - R4（取最小镜像）
d = |r|
F1 = -k × (d - d0) × r/d
F4 = -F1

其他原子没有直接受到这个距离偏置，所以本例的对应偏置力为0；这不表示它们在真实MD中没有力。两个受力原子的合力也应为0。

第二种检查只使用相邻扰动帧的偏置能：

Fi,x ≈ -[V(Ri,x+h) - V(Ri,x-h)] / (2h)
h = 1e-5 nm

独立差分不知道PLUMED内部导数，因而可以帮助发现符号、单位或原子编号的错误。

本次43帧中，PLUMED力与解析式的最大绝对差约1.28e-12 kJ·mol^-1·nm^-1；起点21个原子坐标分量的能量差分与PLUMED力最大差约3.99e-8 kJ·mol^-1·nm^-1。7个原子的结果按原子1x、1y、1z、原子2x…排列。右图显示的是实际差分误差，没有把它画成总物理力。

5 怎样用于自己的体系#

先确认原子编号和周期边界，再选择有明确物理含义的距离、角度或其他CV。目标与力常数应结合体系和单位确定；过强的偏置可能要求更小时间步，过弱则不能有效限制采样范围。

需要umbrella sampling时，通常还要建立多个中心、分别平衡和采样、检查相邻窗口重叠，再进行WHAM等无偏重建。本篇只提供单中心偏置和导数操作参考，不把100步示例当作这些统计步骤已完成。多窗口实操可以继续参考GROMACS伞形采样与WHAM。

回到PLUMED入口查看完整输入与提交流程，或按driver教程提取现有轨迹。