LAMMPS自由能基准与数据工具

本组在鸢算科技集群实际运行了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 为什么另做一个双原子参照#

108粒子流体需要有限轨迹估计系综平均。为了先核对能量导数、指数平均和BAR实现，本包另用双原子径向模型：相对距离积分区间为0.9–2.4，温度1.5，ε从0.8到1.2。

参考配分函数的相对坐标积分是：

Z(ε) ∝ ∫[0.9,2.4] r² exp[−Uε(r)/T] dr
ΔA = −T ln[Z(1.2)/Z(0.8)]

r²来自三维相对坐标的径向Jacobian，角向与不变的动能因子在比值中抵消。该壳层积分定义了参照模型的相空间；每个积分节点实际调用原生LAMMPS求能量、力、TI导数和FEP扰动量。它与108粒子流体是两个不同模型，参照数值不用于替代流体的ΔA。

2 原生单点如何给出基准数据#

完整输入在 inputs/radial-reference/in.lammps，5个ε各64个径向节点，共320个原生单点。下面给出初始化和第一个实际节点，后续节点由包内输入继续执行：

units lj
atom_style atomic
atom_modify map array sort 0 0.0
boundary f f f
region box block -4 4 -4 4 -4 4 units box
create_box 1 box
create_atoms 1 single 0 0 0 units box
create_atoms 1 single 1 0 0 units box
mass 1 1
pair_style lj/cut 2.5
pair_coeff * * 1 1
pair_modify shift yes
variable lam equal .8
fix alchemy all adapt/fep 0 pair lj/cut epsilon 1 1 v_lam scale no reset no
thermo 1000000
variable one equal 1.0
variable delta equal .1
variable reverse equal -.1
compute derivative all ti lj/cut 1 v_lam v_one
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
thermo_modify norm no format float %.16g
variable lam equal 0.8
set atom 2 x 0.9005212186981708
reset_timestep 1
run 0
print "0.8 0.9005212186981708 $(pe:%.16g) $(c_derivative:%.16g) $(c_forward[1]:%.16g) $(c_forward[2]:%.16g) $(c_reverse[1]:%.16g) $(c_reverse[2]:%.16g) $(fx[2]:%.16g)" file points.dat screen no

每点使用独立Step编号，让原生compute在新坐标上重新求值。输入文件中有完整节点，QUADRATURE.json保存其数值权重；这里的节点序列是积分数据。

在已分配的计算资源内单独运行：

module purge
module load lammps/30Sep2026-p100-gpu
mkdir manual-reference
cp inputs/radial-reference/in.lammps manual-reference/
cp inputs/radial-reference/QUADRATURE.json manual-reference/
cd manual-reference
lmp -in in.lammps
cd ..

整包输出：

head -6 work/cases/radial-reference/points.dat
head -6 work/analysis/radial-native-points.csv

列依次为λ、r、原生能量、TI导数、前向ΔU、前向指数项、反向ΔU、反向指数项以及第二个粒子的径向力。

3 本次独立数值比较#

TI保留有限λ积分间隔带来的约3.55×10⁻⁸差异，没有把所有方法强行写成相同数值。原生能量和径向力还与同一截断平移LJ表达式核对，力也用坐标能量差分交叉比较。

python - <<'PYDATA'
import json
from pathlib import Path
a=json.loads(Path('work/analysis/FREE-ENERGY.json').read_text())['radial_benchmark']
for k in ['partition256','native_FEP_forward','native_FEP_reverse','native_BAR','TI_five_state_Simpson']:
    print(k,a[k])
PYDATA

4 实际108粒子结构导出#

python export_windows.py work
head -8 work/analysis/window-1.0-final.xyz
python - <<'PYDATA'
from ase.io import read
p=read('work/analysis/window-1.0-final.xyz')
print('Particles:',len(p))
print('First position:',p.positions[0])
print('Coordinate units:',p.info['coordinate_units'])
PYDATA

每个窗口导出实际终态的108个粒子、周期盒、particle_id和lj_force属性。X为占位物种，坐标与力保持约化单位；没有转换成真实元素的Å或eV/Å。

5 工具与输入之间的关系#

• analyze.py：读原生样本与单点，生成TI窗口表、FEP区间表、BAR和径向基准。
• export_windows.py：从实际终态dump导出可打开的XYZ结构，保留粒子ID和单位标记。
• plot_results.py：从实际CSV/JSON和结构出四张PNG。

python analyze.py work
python export_windows.py work
python plot_results.py work
ls work/figures

这些工具不改变原生轨迹和输入参数。进一步操作见 TI窗口与完整输入 和 FEP/BAR指数平均。