LAMMPS 伞形采样：窗口到WHAM自由能

1. 下载并看参考结果#

下载完整输入、22个窗口的实际采样、WHAM脚本和图片，保存为 example.zip：

python3 -m zipfile -e example.zip .
cd lammps-umbrella-examples
bash results.sh
bash plot.sh

本例使用LAMMPS30Sep2026 Intel CPU MPI构建，模块 lammps/30Sep2026-8173m；后处理用集群 ase/3.29.0 中的NumPy、SciPy、Matplotlib。实际采样采用兼容CPU上的同一程序，不把小模型耗时当作CPU节点性能数据。

2. 先确定反应坐标和原势#

一个粒子只沿x运动，约化温度T=0.2，原势为 U(x)=16*((x-1.5)^2-0.25)^2。两势阱x=1和2、势垒x=1.5，解析势垒为1。它是一维教学模型，不是材料力场或真实离子对；没有球坐标r²测度，也就没有额外径向Jacobian修正。

窗口中心从0.9到2.1、间隔0.12，共11个；每个窗口添加 0.5*K*(x-x0)^2，K=40。这里的K与某些命令直接写 K*(x-x0)^2的约定不同，不能混用。

variable V atom 16*((x-1.5)^2-0.25)^2
variable Fx atom -64*(x-1.5)*((x-1.5)^2-0.25)
fix external all addforce v_Fx 0 0 energy v_V
fix_modify external energy yes
fix umbrella all spring tether 40.0 0.90000000 NULL NULL 0.0
fix_modify umbrella energy yes
fix planar all setforce NULL 0.0 0.0
fix integrate all nve
fix bath all langevin .2 .2 .2 74511 zero no
timestep .002
run 10000
reset_timestep 0
variable X equal x[1]
fix samples all ave/time 10 1 10 v_X f_umbrella &
    file samples.dat format " %.16g"
run 200000

完整输入还创建粒子、定义一维温度和输出，直接查看：

cat inputs/w00-s74511/in.lmp

每窗先平衡10000步，再采样200000步，每10步保存一次。原生文件还写入time=0这一行；分析明确去掉它，每窗使用20000个生产样本，两套种子共440000个样本。

3. 提交各窗口#

cat job.slurm
sbatch job.slurm
squeue -u "$USER"

数组0–21对应11个窗口的两套种子，最多两个同时运行。CPU分区按现行56物理核规则申请，小模型只启动1个MPI进程；资源繁忙时按Slurm排队，不要直接在登录终端运行整组任务。各窗口写到独立 work/wXX-sYYYYY，不会混用文件。

4. 先看重叠，再拼接#

参考数据已随包给出，马上可以重算后处理：

module purge
module load ase/3.29.0
python3 collect.py
head -n 8 reference-results/analysis/umbrella-pmf.csv

自己22个窗口都完成后：

python3 collect.py --data-root work
head -n 8 work/analysis/umbrella-pmf.csv

WHAM使用0.8–2.2范围、0.01间隔的直方图，同时考虑各窗口偏置和样本数，迭代求出无偏概率P(x)，再算 F(x)=-T*ln(P(x))+常数。不能把各窗口的 -ln(histogram)直接首尾相接。

左图只显示一套种子的11个直方图；右图的主曲线使用两套共22窗。自由能零点通过一个共同常数对齐解析曲线，只移动纵轴，不能改变势垒或曲线形状。

5. 这次得到什么#

解析参考势垒为1.000；RMS在x=0.98–2.02的已采样区间计算。两套种子势垒差约0.0351，说明只看一条平滑曲线会隐藏采样差异。

这个短例子可用于学会全流程。真实材料还需要有意义的反应坐标、足够平衡、窗口重叠、相关时间和误差估计；距离坐标、角度或约束坐标应使用对应的测度约定。自由能势垒本身也不等于反应速率。偏置定义可查LAMMPS spring。