LAMMPS自由能微扰与BAR实操

本组在鸢算科技集群实际运行了5个108粒子LJ流体窗口，耦合参数λ依次为0.8、0.9、1.0、1.1、1.2；每个窗口先平衡5000步，再取样10000步。同时用双原子原生单点做独立数值积分基准。

使用LJ约化单位：温度1.5、数密度0.7、dt=0.002，截断半径2.5并平移势能，统一不加长程尾修正。本文自由能为NVT条件下的Helmholtz自由能差ΔA，总体系量；数值不是eV、kcal/mol或每原子结果，也不直接对应某种化学反应。

下载完整输入、实际原始结果、结构和分析工具：

mkdir free-energy-demo
cd free-energy-demo
curl -fL https://yuansuan.top/assets/downloads/full-function/lammps-free-energy-wave86-67e57d88b57e.zip -o examples.zip
unzip examples.zip
sbatch submit.slurm

脚本按顺序计算五个窗口和双原子基准，再生成CSV、XYZ和四张图。模块为 lammps/30Sep2026-p100-gpu，构建使用GCC13.3、OpenMPI5.0.10、CUDA12.8；数据工具使用 ase/3.29.0。这套TI/FEP脚本走CPU计算路径，提交模板申请一个CPU核与一个P100配额。以下直接调用lmp的命令在已分配计算资源中执行。

计算完成后，在解压目录运行：

module purge
module load ase/3.29.0
python analyze.py work
python export_windows.py work
python plot_results.py work

参考数据在reference/，后处理参数work可换成reference；再次计算请解压到新目录，保留既有work。

1 原生FEP输出什么#

每个固定λ窗口都计算ε的±0.1扰动，原子坐标保持原帧位置。本例ε=λ，因此这个数值也对应相邻λ窗口。

variable delta equal .1
variable reverse equal -.1
compute forward all fep 1.5 pair lj/cut epsilon 1 1 v_delta tail no volume no
compute reverse all fep 1.5 pair lj/cut epsilon 1 1 v_reverse tail no volume no

c_forward[1] 是U新−U旧，c_forward[2] 是 exp(−ΔU/T)；第三分量为盒体积。本例NVT使用volume no，温度为LJ约化1.5。tail no与各窗口的势能约定一致。

完整模型、平衡和取样输入见 TI教程；资料包的每个窗口已同时输出前向与反向数据。

2 单帧扰动与自由能不同#

自由能用整个原态系综的指数平均：

ΔA_A→B = −T ln <exp[−(U_B−U_A)/T]>_A

平均ΔU与该公式不同；也不能先平均ΔU再取指数。配套工具使用稳定的log-mean-exp，保留不同帧的权重。

先查看原生列：

head -8 work/cases/window-0.8/samples.dat
head -8 work/cases/window-0.9/samples.dat

用第一段0.8→0.9直接计算：

python - <<'PYDATA'
import numpy as np
x=np.loadtxt('work/cases/window-0.8/samples.dat',comments='#')
x=x[x[:,0]>0]
logw=-x[:,2]/1.5
m=logw.max()
forward=-1.5*(m+np.log(np.exp(logw-m).mean()))
print('Forward Delta A =',forward)
PYDATA

本次前向结果约−39.884761。反向应从0.9窗口的负扰动取样，再把B→A结果变号，才能与A→B同向比较。

3 BAR如何使用两边的数据#

analyze.py 同时读取相邻两个窗口的原生ΔU，用Bennett接受率方程求根。此处两个方向各500个样本，采用等样本数的两态公式；实现随包提供，不另要求安装pymbar。

python analyze.py work
cat work/analysis/fep-segments.csv

ESS衡量指数权重是否集中于少数帧。这里按重加权权重计算，未扣除轨迹的时间相关，不能将它直接当作独立样本数。

4 四种总自由能估计#

四段结果求和，本次得到：

图左以BAR估计为画图参照，仅便于看差异；BAR也需要充分取样和相邻状态重叠。四种结果接近，说明这段数据适合展示操作流程，不能代替时长、窗口和统计误差收敛研究。

python - <<'PYDATA'
import json
from pathlib import Path
a=json.loads(Path('work/analysis/FREE-ENERGY.json').read_text())['fluid']
for k in ['TI_Simpson','FEP_forward','FEP_reverse','BAR']:
    print(k,a[k])
PYDATA

5 换自己的体系时改哪里#

先定义两个目标Hamiltonian和同一单位/能量约定，再选择耦合路径、窗口和每窗平衡/取样。保留两边原生ΔU，观察权重集中度与方向差异；必要时增加中间窗口和取样。替换成电荷或其他势参数前，检查该pair style的参数提取支持及体系电中性、长程算法等条件。

双原子方法基准与结构工具提供可独立核对的数值例子。参数及支持项见 LAMMPS compute fep。

位置优化算法与微扰精度#

五种位置优化算法给出完整输入、实际迭代点和停止条件；GPU/Kokkos最小化说明GPU力与算法后缀、受约束力与物理力。当前30Sep2026已移除fire/old，请使用fire。

TI/FEP的GPU精度用CPU、混合精度GPU和双精度Kokkos对照实际参数变化，微扰步长选择比较ε和σ的正反能量差。固定坐标的能量导数与平衡集合的自由能导数分别说明；120个单点不当作平衡取样。