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 明确λ如何进入势函数#

本例直接设置ε=λ，σ=1保持不变。因此 Uλ=λU0，每帧的能量导数是 dUλ/dλ=Uλ/λ。TI对各窗口的平均导数积分，得到0.8→1.2的ΔA。

左图来自λ=1.0实际终点，右图是同一结构的最小镜像粒子对距离计数；右图没有作RDF归一化。所有坐标均为LJ约化长度，结构中的X为占位标签。

2 λ=1.0窗口的完整输入#

units lj
atom_style atomic
atom_modify map array sort 0 0.0
boundary p p p
lattice fcc .7
region box block 0 3 0 3 0 3 units lattice
create_box 1 box
create_atoms 1 box
mass 1 1
pair_style lj/cut 2.5
pair_coeff * * 1 1
pair_modify shift yes
neighbor .3 bin
neigh_modify every 1 delay 0 check yes
variable lam equal 1
fix alchemy all adapt/fep 0 pair lj/cut epsilon 1 1 v_lam scale no reset no
velocity all create 1.5 5295 mom yes rot no dist gaussian
fix bath all nvt temp 1.5 1.5 .2
timestep .002
thermo 1000
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
run 5000
reset_timestep 0
fix sample all ave/time 20 1 20 c_derivative c_forward[1] c_forward[2] c_reverse[1] c_reverse[2] c_thermo_pe c_thermo_temp file samples.dat format " %.16g"
dump trajectory all custom 200 trajectory.dump id type x y z vx vy vz
dump_modify trajectory sort id format float %.16g
run 10000
unfix sample
unfix bath
run 0
write_dump all custom final.dump id type x y z fx fy fz modify sort id format float %.16g
print "FINAL $(pe:%.16g) $(c_derivative:%.16g)" file final-energy.dat screen no

其他窗口的完整输入在 inputs/window-0.8/ 至 inputs/window-1.2/。同一窗口内λ固定，先平衡，再把Step重置为0开始记录；Step0样本单独保留，统计使用后面的500个样本。

3 绝对赋值与多段run#

variable lam equal 1.0
fix alchemy all adapt/fep 0 pair lj/cut epsilon 1 1 v_lam scale no reset no
compute derivative all ti lj/cut 1 v_lam v_one

scale no 把变量值直接作为新ε，因此平衡、生产和最后run0都保持ε=λ。scale yes 使用初始化记录参数的相对倍数；跨多段run复用时，需要认真核对该基准。reset no 表示运行结束后保留当前参数。本例选择绝对赋值，避免重复缩放改变窗口含义。

adapt/fep 的N=0表示在运行初始化时赋值；compute ti 中 v_lam 为能量缩放函数、v_one=1 为其λ导数。这里的U/λ关系适用于该线性ε路径；改σ、改变多项参数或换势函数时，应重新定义耦合路径和导数。

4 单独运行一个窗口#

module purge
module load lammps/30Sep2026-p100-gpu
mkdir manual-window
cp inputs/window-1.0/in.lammps manual-window/
cd manual-window
lmp -in in.lammps
cd ..

在整包运行后查看数据：

head -8 work/cases/window-1.0/samples.dat
cat work/analysis/ti-windows.csv

每20步记录一次。关键列为Step、dU/dλ、前向ΔU、前向指数权重、反向ΔU、反向指数权重、势能和温度。

5 本次五窗口的实际均值#

每窗500个生产样本分成10块，图中的误差线用于观察这段数据的块间变化；它没有证明时间相关已经消除或自由能已经收敛。

这里λ步长0.1，五点复合Simpson公式为：

ΔA ≈ (0.1/3) × (m0 + 4m1 + 2m2 + 4m3 + m4)

其中mi为各窗口平均导数。本次得到 ΔA≈−165.595522，是108粒子体系的约化总自由能差。

python - <<'PYDATA'
import numpy as np
x=np.genfromtxt('work/analysis/ti-windows.csv',delimiter=',',names=True)
value=(0.1/3)*np.dot([1,4,2,4,1],x['mean_derivative'])
print('TI Delta A =',value)
PYDATA

使用相同截断、平移和尾修正约定，才能比较各窗口。后续应检查平衡时长、生产时长、λ间隔和不同随机初态；FEP/BAR教程给出另一组同数据估计，双原子基准与工具说明独立参照。

原生参数见 compute ti 和 fix adapt/fep。

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

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

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