LAMMPS温压控制的GPU与Kokkos操作

同样的温控/压控输入切到GPU或Kokkos时，需要配套修改势、积分fix和启动参数；混合精度与双精度也会影响能量、力及后续长轨迹。本文使用12个方法、3条后端路径的真实数据，演示如何操作并解释差别。

1. 下载与复算#

下载全部原生输入、实际轨迹、CPU单点对照与绘图工具：

mkdir thermostat-backends
cd thermostat-backends
curl -fL https://yuansuan.top/assets/downloads/full-function/lammps-thermostats-wave91-df7f1ec315d1.zip -o examples.zip
unzip examples.zip
sbatch submit.slurm

脚本申请1个P100和1个CPU核，顺序完成12个CPU计算条件、12个GPU力对照和12个Kokkos对照；GPU末态另做相同坐标CPU单点。每次复算请使用新的解压目录。直接调用lmp的命令在已分配的计算资源中执行，登录节点可做后处理。

模块为 lammps/30Sep2026-p100-gpu、lammps/30Sep2026-p100-kokkos，编译使用GCC13.3、OpenMPI5.0.10、CUDA12.8；Python数据工具使用 ase/3.29.0。本组CPU执行路径用于数值对照，不是CPU节点性能数据。

总计36个MD条件，每个108粒子、10000步。CPU的12个条件作为主要方法示例，GPU与Kokkos为后端对照；相同坐标的12个CPU单点是力检查步骤，不另当作12个新MD体系。

2. GPU包：CPU侧控温配合GPU势力#

以 inputs/csvr-gpu/in.lammps 为例：

package gpu 1
newton off
units lj
atom_style atomic
boundary p p p
lattice fcc .7
region box block 0 3 0 3 0 3
create_box 1 box
create_atoms 1 box
mass 1 1
pair_style lj/cut/gpu 2.5
pair_coeff * * 1.0 1.0
pair_modify shift yes
neighbor .3 bin
neigh_modify every 1 delay 0 check yes
velocity all create 1.8 20261029 mom yes rot no dist gaussian
fix integrator all nve
fix bath all temp/csvr 1.0 1.0 .2 20261102
timestep .002
thermo 100
thermo_style custom step time temp pe ke etotal press vol density ecouple econserve
thermo_modify norm no format float %.16g
dump trajectory all custom 1000 trajectory.dump id type xu yu zu vx vy vz fx fy fz
dump_modify trajectory sort id format float %.16g
run 10000
unfix integrator
unfix bath
run 0
write_dump all custom physical.dump id type x y z fx fy fz modify sort id format float %.16g
print "PHYSICAL $(pe:%.16g) $(vol:%.16g)" file physical-energy.dat screen no
print "COMPLETE" file complete.dat screen no

GPU包的势为 lj/cut/gpu，控温仍用 temp/csvr 配合 nve；不假设存在每个fix的 /gpu 版本。这里用 package gpu 1 与 newton off。

module purge
module load lammps/30Sep2026-p100-gpu
cd inputs/csvr-gpu
lmp -in in.lammps
cd ../..

3. Kokkos：显式后缀与启动命令#

inputs/csvr-kk/in.lammps：

package kokkos neigh half
newton on
units lj
atom_style atomic/kk
boundary p p p
lattice fcc .7
region box block 0 3 0 3 0 3
create_box 1 box
create_atoms 1 box
mass 1 1
pair_style lj/cut/kk 2.5
pair_coeff * * 1.0 1.0
pair_modify shift yes
neighbor .3 bin
neigh_modify every 1 delay 0 check yes
velocity all create 1.8 20261029 mom yes rot no dist gaussian
run_style verlet/kk
fix integrator all nve/kk
fix bath all temp/csvr/kk 1.0 1.0 .2 20261102
timestep .002
thermo 100
thermo_style custom step time temp pe ke etotal press vol density ecouple econserve
thermo_modify norm no format float %.16g
dump trajectory all custom 1000 trajectory.dump id type xu yu zu vx vy vz fx fy fz
dump_modify trajectory sort id format float %.16g
run 10000
unfix integrator
unfix bath
run 0
write_dump all custom physical.dump id type x y z fx fy fz modify sort id format float %.16g
print "PHYSICAL $(pe:%.16g) $(vol:%.16g)" file physical-energy.dat screen no
print "COMPLETE" file complete.dat screen no

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

输入配套使用 atomic/kk、lj/cut/kk、nve/kk、temp/csvr/kk 和 verlet/kk，采用半邻居表/Newton on。其他方法的实际后缀见包内各份完整输入，不只在命令行添加一个开关。

4. 为什么同一种方法最后一帧会不同#

微小舍入差异在非线性MD中会逐渐放大，随机控温路径也未必产生逐位相同的随机序列。不能仅因最后一帧坐标或温度不同，就判定某条路径物理上错误；也不能仅因均值相近就宣称严格系综、长时间统计或材料物性已经收敛。

这里分别核对原生输出的动能/温度自由度、体积/密度、固定或变化的晶胞，以及实际末态的势能和力；没有要求36条轨迹逐帧重合。这些数据也不是速度基准。

5. 用相同坐标检查混合精度#

各自末态坐标与解析周期平移LJ模型对照，CPU最大能量/力绝对偏差约1.20×10⁻¹²、1.09×10⁻¹²；Kokkos约2.28×10⁻¹³、1.10×10⁻¹²。混合精度GPU的最大能量偏差约2.15×10⁻⁵，最大力分量偏差约5.46×10⁻⁴，单位都是本例LJ约化单位。

这些偏差不能直接换算成eV/Å，也不意味着任意势、原子数或力阈值都适用。如果研究需要更严格的力或小能量差，优先做相同坐标CPU或双精度Kokkos核对，结合实际势函数测试精度。

每个GPU结果另有 cpu-control.in，从原生末态dump读取相同盒子和坐标，用CPU双精度重新算一次势能和力。在已分配资源内可独立执行：

module purge
module load lammps/30Sep2026-p100-gpu
cd work/cases/csvr-gpu
lmp -in cpu-control.in -log cpu-control-manual.log
cd ../../..

这一步没有MD推进，只是固定坐标单点；脚本先移除温控/压控fix再取物理力，避免把随机热浴力当作势能梯度。

6. 提取实际数据与图片#

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

原生输出在 work/cases/，逐时CSV、汇总和含完整周期晶胞的XYZ在 work/analysis/，6张图在 work/figures/。reference/ 是本文实际数据；未重算时，可以把以上 work 改为 reference 重现图表。X为LJ粒子占位符，XYZ坐标、晶胞和 lj_velocity 仍是LJ约化单位，不能直接当成Å、eV或真实元素。

python - <<'PYCODE'
import csv
with open('work/analysis/summary.csv') as f:
    for r in csv.DictReader(f):
        if r['backend'] == 'gpu':
            print(r['method'], r['final_energy_formula_error'], r['final_force_formula_max_error'])
PYCODE

选择具体控温方法看温控操作，设置压控方向和倾斜自由度看压控与晶胞。

原生后端说明：GPU包、Kokkos。