LAMMPS自由能微扰的GPU精度对照

计算自由能微扰前，先确认参数变化确实生效，并检查GPU数值精度对小能量差的影响。本文在鸢算科技集群完成CPU、混合精度GPU与双精度Kokkos的同坐标TI/FEP对照：ε扫描108个单点、σ扫描12个单点，共120个单点，归为一套完整数值对照流程。

这些点的坐标是指定的，不是热平衡取样，不能把它们直接平均成自由能。真正的NVT窗口、TI积分和FEP/BAR指数平均见TI工作流和FEP/BAR工作流。

1. 准备、提交与结果位置#

下载完整输入、实际结果和绘图工具。在自己的工作目录解压并提交：

mkdir fep-precision-demo
cd fep-precision-demo
curl -fL https://yuansuan.top/assets/downloads/full-function/lammps-fep-precision-wave88-f278e5aaeb2c.zip -o examples.zip
unzip examples.zip
sbatch submit.slurm

提交脚本申请1个P100和1个CPU核，顺序执行CPU、GPU和Kokkos对照，再生成CSV和图片。每次重新计算请使用新的解压目录，避免覆盖已有结果。直接调用 lmp 的命令在已分配的计算资源中执行。

使用 lammps/30Sep2026-p100-gpu 与 lammps/30Sep2026-p100-kokkos，GCC13.3、OpenMPI5.0.10、CUDA12.8；后处理使用 ase/3.29.0。这里比较数值结果，双原子小模型不适合评价GPU加速比。

参数使用LJ约化单位，截断2.5，势能平移为零，无尾修正；FEP的温度参数T=1.5。ε对照包含距离0.95、1.2、2.1，ε=0.8、1.0、1.2与4个微扰幅度；σ对照固定距离1.2、ε=1、σ=1。

2. CPU原生输入：每个微扰值都重新求值#

下面是完整的σ扫描输入 inputs/r1.20-cpu-sigma/in.lammps：

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.2 0 0 units box
mass 1 1
pair_style lj/cut 2.5
pair_coeff * * 1.0 1.0
pair_modify shift yes
variable lam equal .8
variable one equal 1
variable delta equal .1
variable reverse equal -v_delta
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
compute forward all fep 1.5 pair lj/cut sigma 1 1 v_delta tail no volume no
compute reverse all fep 1.5 pair lj/cut sigma 1 1 v_reverse tail no volume no
thermo 1000000
thermo_modify norm no format float %.16g
variable lam equal 1.0
variable delta equal 0.01
reset_timestep 1
run 0
print "1.0 1.2 0.01 $(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
variable lam equal 1.0
variable delta equal 0.001
reset_timestep 2
run 0
print "1.0 1.2 0.001 $(pe:%.16g) $(c_derivative:%.16g) $(c_forward[1]:%.16g) $(c_forward[2]:%.16g) $(c_reverse[1]:%.16g) $(c_reverse[2]:%.16g) $(fx[2]:%.16g)" append points.dat screen no
variable lam equal 1.0
variable delta equal 1e-05
reset_timestep 3
run 0
print "1.0 1.2 1e-05 $(pe:%.16g) $(c_derivative:%.16g) $(c_forward[1]:%.16g) $(c_forward[2]:%.16g) $(c_reverse[1]:%.16g) $(c_reverse[2]:%.16g) $(fx[2]:%.16g)" append points.dat screen no
variable lam equal 1.0
variable delta equal 1e-07
reset_timestep 4
run 0
print "1.0 1.2 1e-07 $(pe:%.16g) $(c_derivative:%.16g) $(c_forward[1]:%.16g) $(c_forward[2]:%.16g) $(c_reverse[1]:%.16g) $(c_reverse[2]:%.16g) $(fx[2]:%.16g)" append points.dat screen no
print "COMPLETE" file complete.dat screen no

fix adapt/fep ... epsilon ... v_lam scale no reset no 把ε赋为绝对值。本例循环改变ε时，不能把这行改成 scale yes 后仍解释为绝对ε，否则多次求值可能改变耦合路径。

compute ... ti lj/cut 1 v_lam v_one 计算的是对ε的能量导数。即使当前文件扫描σ，它仍是ε导数，不能拿它冒充σ导数。σ前向与反向能量变化分别来自 compute fep ... sigma ... v_delta 和 v_reverse。

每个点使用不同的 reset_timestep，再执行 run 0，确保本次参数和计算结果重新求值。没有MD积分，也没有温度平衡过程。

3. 切换GPU与Kokkos#

实际GPU输入使用 package gpu 1、newton off、pair_style lj/cut/gpu 2.5；TI/FEP/adapt里的势名称也改为 lj/cut/gpu。Kokkos版本配套使用 atomic/kk、lj/cut/kk、半邻居表、Newton on和 verlet/kk，并用 -k on g 1 启动。

module purge
module load lammps/30Sep2026-p100-gpu
cd inputs/r1.20-gpu-sigma
lmp -in in.lammps
cd ../..
module purge
module load lammps/30Sep2026-p100-kokkos
cd inputs/r1.20-kk-sigma
lmp -k on g 1 -in in.lammps
cd ../..

这些显式输入包含在下载包中；可以直接查看，避免只加启动后缀却漏改 compute 或 fix 中的势名称：

cat inputs/r1.20-gpu-sigma/in.lammps
cat inputs/r1.20-kk-sigma/in.lammps

当前 lj/cut/gpu 实现支持这些系数的重新初始化；本例用每个点的解析能量、力和实际±δ能量差逐项对照，确认参数修改后结果生效。其他势能否用于同类FEP，需要核对其参数与重新初始化支持。

4. 实际结果与能量差精度#

大尺度能量曲线基本重合，不意味着任意小的能量差都可靠。ε微扰δ=1e-7时，混合精度GPU部分前向能量差变成零，前向导数最大相对误差为1；CPU/Kokkos的中心导数最大相对误差约6.81×10⁻¹⁰。GPU的ΔU绝对误差最大约1.98×10⁻⁷，除以1e-7后就会被明显放大。

这描述的是本例GPU包混合精度、势参数、坐标和微扰尺度的组合，不意味着所有GPU自由能计算都不可靠。选择更合适的δ、检查精度路径，并在代表性构型上做CPU或双精度对照。

5. 终端提取、表格导出与画图#

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

work/cases/ 保存本次原始输出，work/analysis/ 保存表格，work/figures/ 保存高清图。包内 reference/ 是本文实算结果；未重算时，可以把上述命令中的 work 改为 reference 来重新分析和绘图。

head -n 4 work/cases/r1.20-gpu-sigma/points.dat
python - <<'PYCODE'
import csv
with open('work/analysis/derivative-errors.csv') as f:
    for r in csv.DictReader(f):
        if r['parameter'] == 'epsilon' and float(r['delta']) == 1e-7:
            print(r['backend'], r['max_forward_relative_derivative_error'], r['max_central_relative_derivative_error'])
PYCODE

points.dat 的10列依次为ε、距离、δ、基态能量、TI的ε导数、前向ΔU、前向Boltzmann因子、反向ΔU、反向Boltzmann因子、粒子2的x方向力。

export_fep_tables.py 从总表无损拆出 epsilon-native-controls.csv 的108行和 sigma-native-controls.csv 的12行。不要对这些不同坐标点做均匀平均来报告自由能；后续应在真实平衡轨迹上进行适当的集合平均、重叠和取样收敛检查。微扰步长如何选见差分步长实操。

原生说明：compute fep、compute ti、fix adapt/fep。