LAMMPS微扰与差分步长选择
ε与σ两类参数的原生正反微扰、中心差分、截断误差和GPU有限精度。
阅读本篇目录
微扰和数值差分不是步长越小越好。本文使用同坐标的真实LAMMPS输出，分别改变线性参数ε与非线性参数σ，把截断误差、有限精度和自由能采样的含义分开说明。
1. 取得数据与重新运行#
下载完整输入、实际结果和绘图工具。在自己的工作目录解压并提交：
mkdir fep-step-demo
cd fep-step-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加速比。
完整CPU/GPU/Kokkos输入、10列原生输出的含义见GPU精度对照。本例是LJ约化模型，不对应真实物质的自由能；ε扫描9个固定构型/参数状态，σ扫描1个固定状态，每条路径使用4个δ。
2. 从原生能量变化提取导数#
令原生前向输出为 ΔU+ = U(p+δ)−U(p)，反向输出为 ΔU− = U(p−δ)−U(p)。前向差分是 ΔU+/δ，中心差分是 (ΔU+−ΔU−)/(2δ)；反向ΔU已经带有实际符号，不要再次误取负号。
关键原生配置为：
variable delta equal 0.001
variable reverse equal -v_delta
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
包中 analyze.py 用实际输入的平移LJ势解析导数做参照，并保留±δ两侧原生能量变化，没有用解析能量代替LAMMPS数据。
先生成主表和图：
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 来重新分析和绘图。
python - <<'PYCODE'
import csv
with open('work/analysis/fep-precision.csv') as f:
for r in csv.DictReader(f):
if r['parameter'] == 'sigma' and r['backend'] == 'cpu':
d = float(r['delta'])
up = float(r['native_forward_deltaU'])
down = float(r['native_reverse_deltaU'])
print(d, up/d, (up-down)/(2*d))
PYCODE
3. 线性ε：主要观察小差值的有限精度#
平移LJ势在固定r、σ和截断下对ε线性。因此本例ε的中心与前向差分在精确算术中都等于解析导数；看到小δ误差增长，主要是在观察能量差的数值消减，而不是复杂的积分截断误差。
| 路径 | 微扰量δ | 前向导数最大相对误差 | 中心导数最大相对误差 |
|---|---|---|---|
| cpu | 0.1 | 4.82879e-15 | 3.59352e-15 |
| cpu | 0.001 | 2.8021e-13 | 1.63698e-13 |
| cpu | 1e-05 | 2.04986e-11 | 9.20389e-12 |
| cpu | 1e-07 | 1.94989e-09 | 6.80554e-10 |
| gpu | 0.1 | 8.51807e-07 | 8.51807e-07 |
| gpu | 0.001 | 0.000147543 | 7.93961e-05 |
| gpu | 1e-05 | 0.0241966 | 0.0116863 |
| gpu | 1e-07 | 1 | 0.397108 |
| kk | 0.1 | 4.82879e-15 | 3.59352e-15 |
| kk | 0.001 | 2.8021e-13 | 1.63698e-13 |
| kk | 1e-05 | 2.04986e-11 | 9.20389e-12 |
| kk | 1e-07 | 1.94989e-09 | 6.80554e-10 |
表中每个误差是9个状态中的最大值。GPU包混合精度在较大的δ下仍能得到接近的导数，小到1e-7时前向差分可能完全丢失变化量。CPU与Kokkos在这幅图上几乎重合。
4. 非线性σ：截断误差与舍入误差竞争#
σ改变的是相互作用长度尺度，势能包含σ的6次和12次方。此时 U/ε 仍是ε导数，不是σ导数。对σ需使用原生±δ的能量变化，或对应的正确解析导数。
| 路径 | 微扰量δ | 前向导数最大相对误差 | 中心导数最大相对误差 |
|---|---|---|---|
| cpu | 0.01 | 0.0410696 | 0.00282624 |
| cpu | 0.001 | 0.00384475 | 2.82489e-05 |
| cpu | 1e-05 | 3.8167e-05 | 2.81947e-09 |
| cpu | 1e-07 | 3.79895e-07 | 5.24405e-10 |
| gpu | 0.01 | 0.0410713 | 0.00282717 |
| gpu | 0.001 | 0.00381222 | 1.19294e-05 |
| gpu | 1e-05 | 0.000197902 | 0.00136363 |
| gpu | 1e-07 | 0.0674145 | 0.0674145 |
| kk | 0.01 | 0.0410696 | 0.00282624 |
| kk | 0.001 | 0.00384475 | 2.82489e-05 |
| kk | 1e-05 | 3.8167e-05 | 2.81947e-09 |
| kk | 1e-07 | 3.79895e-07 | 5.24405e-10 |
在双精度路径上，中心差分能明显减少单边差分的截断误差；继续减小δ最终也受有限精度限制。GPU混合精度在δ=1e-7时中心导数相对误差约0.0674。不要从这个双原子模型推出适合所有体系的统一δ值。
5. 迁移到自己的自由能计算#
先选若干代表性构型，固定坐标扫描多组δ，并与CPU或双精度能量/力交叉核对。优先寻找结果相对稳定的步长区间，再检查真实采样中的相邻窗口重叠、相关时间和窗口数；中心差分不能替代统计收敛。
这里得到的是固定坐标的能量参数导数。NVT热力学积分还需要每个参数窗口内的平衡集合平均，FEP需要对Boltzmann权重做集合平均；不能把120个单点当成一条平衡轨迹。下一步可按TI窗口和FEP/BAR的完整操作继续。