LAMMPS多时间步与能量守恒实操

RESPA把不同力放在不同时间尺度计算，能否使用更大的外层步长仍需要数值核对。本文用同一个“强谐振键＋较弱LJ相互作用”双粒子模型，对照Verlet和两层RESPA，并用步长减半与独立高精度常微分方程积分检查结果。

1. 下载、提交与输入模型#

下载完整输入、原生轨迹和绘图工具：

mkdir respa-demo
cd respa-demo
curl -fL https://yuansuan.top/assets/downloads/full-function/lammps-constraints-respa-wave89-a2947ce45eb6.zip -o examples.zip
unzip examples.zip
sbatch submit.slurm

提交脚本申请1个P100和1个CPU核，顺序完成10个CPU条件与3个Kokkos对照，再生成表格、XYZ与4张图。每次复算使用全新解压目录。直接执行lmp的命令在已有计算资源中运行；登录节点可做后处理。

采用 lammps/30Sep2026-p100-gpu、lammps/30Sep2026-p100-kokkos，构建为GCC13.3、OpenMPI5.0.10、CUDA12.8；数据工具使用 ase/3.29.0。所有数值采用LJ约化单位，不能直接换标签当成Å、eV或ps。

模型两粒子质量均为1，初始距离1.1，初速度为零。谐振键的平衡距离1、LAMMPS力常数K=100；非键LJ的ε=0.01、σ=1、截断2.5且平移能量。special_bonds lj 1 1 1 保留相连粒子的LJ作用，不能忽略这行而仍比较同一个模型。

两种积分器均演化到LJ时间4；外层dt=0.002时2000步，dt=0.001时4000步。比较的是相同物理时长。

完整 inputs/respa-dt0.002/model.data：

LAMMPS teaching model; LJ reduced units

2 atoms
1 bonds
0 angles

1 atom types
1 bond types

-6 6 xlo xhi
-6 6 ylo yhi
-6 6 zlo zhi

Masses

1 1

Atoms # molecular

1 1 1 -0.55 0 0
2 1 1 0.55 0 0

Velocities

1 0 0 0
2 0 0 0

Bonds

1 1 1 2

2. 两层RESPA完整输入#

units lj
atom_style molecular
atom_modify map array sort 0 0.0
boundary p p p
read_data model.data
pair_style lj/cut 2.5
pair_coeff * * 0.01 1.0
neighbor .3 bin
neigh_modify every 1 delay 0 check yes
pair_modify shift yes
special_bonds lj 1.0 1.0 1.0
bond_style harmonic
bond_coeff 1 100.0 1.0
run_style respa 2 2 bond 1 pair 2
fix integrator all nve
timestep 0.002
thermo 100
thermo_style custom step time pe ke etotal temp
thermo_modify norm no format float %.16g
dump trajectory all custom 10 trajectory.dump id mol type xu yu zu vx vy vz fx fy fz
dump_modify trajectory sort id format float %.16g
run 2000
print "COMPLETE" file complete.dat screen no

run_style respa 2 2 bond 1 pair 2 表示两层积分：键力位于第1层、LJ pair力位于第2层；外层每步执行2个内层子步。timestep 0.002 是外层步长，内层键力步长为0.001。本例原生日志相应打印 Time steps : 1:0.001 2:0.002。

Verlet对照输入只使用 run_style verlet，所有力在同一时间步更新。两层RESPA仍需要 fix nve 完成积分，不要再额外加第二个NVE fix。

在已分配资源中单独运行：

module purge
module load lammps/30Sep2026-p100-gpu
cd inputs/respa-dt0.002
lmp -in in.lammps
cd ../..

3. 逐步提取实际轨迹与能量#

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

本次原生输出在 work/cases/，结果表和XYZ轨迹在 work/analysis/，高清图在 work/figures/。包内 reference/ 保留本文实算结果；未重新计算时，可以把以上 work 改成 reference 重新分析。XYZ中的X是约化粒子占位符，lj_velocity 保存原生速度，不假定真实元素或物理单位。

python - <<'PYCODE'
import csv
with open('work/analysis/summary.csv') as f:
    for r in csv.DictReader(f):
        if r['kind'] == 'respa':
            print(r['name'], r['max_energy_drift'], r['max_reference_trajectory_error'])
PYCODE

同一种方法将dt减半，能量误差约降为四分之一。RESPA在这个双力级模型中因为内层键力步长更小，能量误差较低；这不等于它在任意力场、分层方式或大步长下都更稳定，也不是时间性能结论。

4. 结果是否真的对应同一个力场#

LAMMPS谐振键能量为 K(r−r0)^2，这里不是另加1/2的写法。分析脚本按这个约定重建势能与动能，并逐点核对原生热力学表；check_results.py 从标量总势能做坐标中心差分，独立检查原始输出的力。

相对距离的独立参照用DOP853积分、相对容差1e-12。两粒子同质量且只有径向内力，相对坐标方程为 r'' = −2 dU/dr。这是数值参照，不是用连续解覆盖LAMMPS轨迹；原生半步误差保留在CSV中。

迁移到真实分子时，先确定哪些项适合较小内层步长，核对势函数、长程项、约束与所选RESPA配置的兼容性，再做代表性NVE能量和时间步测试。不要把此处2:1的层级直接当作所有体系的推荐设置。

原生参数：run_style。