VASP Blue Moon：从约束力积分自由能

本系列用鸢算科技集群实际完成的H₂短轨迹，演示固定键长的约束分子动力学、Blue Moon平均力积分与保留完整状态续算。五个距离，每个距离两条不同随机种子的轨迹；各200步、步长0.5 fs，另有五个独立静态力参照和一条100步续算。

下载输入、原始输出与绘图脚本，保存为example.zip：

python3 -m zipfile -e example.zip .
cd vasp-blue-moon-radial-example
bash results.sh
bash read.sh
bash plot.sh

这三条bash命令读取已经完成的实际结果，不提交新计算。生成的图片在figures/，结果表在analysis-read/。集群环境为VASP 6.6.1标准GPU版，NVHPC 25.3、CUDA 12.8、HPC-X 2.22.1、MKL 2023.2；PBE PAW64 H，ENCUT=520 eV、EDIFF=1E−6，Γ点。P100和V100 16G各承担一套随机种子的窗口，一张卡一个任务。

这是经典原子核的教学模型，使用300 K固定体积、Nosé–Hoover恒温器；100 fs窗口用于展示操作，不能据此宣称任意体系的自由能已经收敛。分子零点能、量子核运动和解离反应势垒不在本例内。

1. 从REPORT取出四列#

本例b_m>之后四列依次为约束乘子λ、质量度量权重|Z|⁻¹ᐟ²、几何项GkT、加权梯度。对于本例的一个距离约束，平均自由能梯度按：

dA/dr = average(column4) / average(column2)

计算。不要只取乘子的平均值，也不要对加权的单帧比值直接做平均来代替一般的加权比值。多维集体坐标和其他类型约束要按对应理论处理。

grep 'b_m>' reference-results/NH-r0.80-s1/REPORT | head
bash read.sh
head analysis-read/BlueMoon-window-gradients.csv

read.py逐个读取原生REPORT和vasprun.xml，舍弃前50步，再计算平均梯度、五段块统计和实际键长偏差。块标准误仅帮助观察波动，不证明这些短轨迹已经独立或采样收敛。

2. 按距离顺序积分#

五个点必须按0.70→0.90 Å排序。脚本以0.70 Å为零点，用相邻窗口平均梯度做梯形积分：

A(r[j+1])-A(r[j])
  = (g[j]+g[j+1])/2 * (r[j+1]-r[j])

梯度为eV/Å、距离为Å，所以积分结果为eV。两套随机种子分别积分，避免只展示一条看上去平滑的曲线。

cat analysis-read/BlueMoon-radial-freeenergy.csv
bash plot.sh

3. 独立参照为什么不是单纯的静态能#

对于本例的孤立双原子经典径向分布，球坐标的体积因子是r²，所以参照为：

A(r) = U(r) - 2*kB*T*ln(r) + constant
dA/dr = dU/dr - 2*kB*T/r

U(r)来自五个独立静态VASP点，dU/dr由沿键方向的原子力得到。只把静态电子能曲线与径向自由能直接相减，会漏掉几何项。这一解析参照仅适合这里的双原子距离模型，不能把-2kBT ln(r)随意加到所有复杂反应坐标上。

本次十个窗口的最大键长偏差约1.37×10⁻⁷ Å；平均梯度与独立径向参照的最大差约0.0322 eV/Å，积分曲线的最大差约0.00729 eV。它们是本次短模型的实际比较值，不是对任意体系的误差保证。

直接下载窗口表 · 下载自由能结果

4. 分析自己新算的结果#

先完成十个窗口和五个静态参照，按完整状态续算完成额外100步，再运行配套完整读取入口：

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

这里read.sh同时检查续算文件；没有完成续算时会指出缺少的文件。参考数据与自己的结果分开保存。

本次先前还尝试过短Andersen轨迹，与径向参照的偏差超过预设比较范围，因此没有纳入本页主结果。这不表示所有Andersen应用都不适用；具体体系仍需比较热化、时间步、恒温参数和采样误差。本篇主例采用上文明确的Nosé–Hoover输入。

更多定义见VASP Blue Moon ensemble。