LAMMPS温控方法与速度重标定实操

选择恒温器时，要分清谁负责时间积分、谁只改变力或速度，以及目标温度是否应在每一时刻都精确相等。本文在鸢算科技集群实际运行了同一108粒子LJ模型的8种温度/积分配置，并配套GPU和Kokkos数值对照。

1. 下载、提交与统一模型#

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

mkdir thermostat-demo
cd thermostat-demo
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节点性能数据。

初始结构为3×3×3 FCC晶胞、108粒子，质量1、数密度0.7，初始速度温度1.8；LJ ε=σ=1，截断2.5，势能平移到截断处为零。所有量是LJ约化单位；dt=0.002，每个条件10000步，总时间20。

除NVE外，本篇控温目标为1.0；用时间10—20的数据作操作性均值对照，不作为已收敛的统计物性。初始FCC和末态如下，不从短轨迹判定材料相图。

2. CSVR完整输入：积分器与控温fix分别设置#

inputs/csvr-cpu/in.lammps：

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 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

fix nve 完成积分，fix temp/csvr 随机重标定速度以控温，二者分工不同。Tdamp=0.2 是LJ时间，不是0.2个时间步；随机种子为20261102。更换为真实体系时按实际units解释阻尼时间。

最后两行 unfix 后的 run 0 只在同一坐标上重算纯势力，没有继续推进MD，便于检查物理力；主MD数据取前面的10000步段。

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

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

3. 按目的切换方法#

NVT fix本身已积分，不要再给同一组叠加NVE。其他只提供温控的fix才另配积分器。周期重标定和Berendsen可以帮助理解输入行为、温度松弛；不要把它们当作严格正则系综涨落采样的通用推荐。

4. 从实际轨迹看温度#

NVE末态温度不等于初始1.8，本例中部分动能转为势能；这不等于总能量不守恒。CPU、GPU、Kokkos的NVE最大相对总能量变化约1.50×10⁻⁴、1.76×10⁻⁴、1.15×10⁻⁴；仍需通过自己的步长和更长轨迹判断研究所需精度。

表中均值和标准差都来自短窗口。rescale每10步输出，其他方法每100步输出，时间相关和取样间隔不同，不能用表中的标准差给恒温器作“涨落优劣”排序。

5. 周期重标定为什么不能只看重设时刻#

temp/rescale 100 ... 每100步重设温度；如果也每100步打印thermo，看到的恰好都是重设后的T=1，会漏掉期间变化。本例最终输入改为每10步打印，保留两次重设之间的波动。不能据恒定的抽样点声称整段轨迹没有温度变化。

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'] == 'cpu':
            print(r['method'], r['mean_temperature_last_half'], r['temperature_sd_last_half'])
PYCODE

更换阶段先 unfix 旧积分器/恒温器，再设置新的方法。需要允许晶胞变化时见压控与晶胞，GPU/Kokkos与力精度见后端实操。本例不宣称任意恒温器在任意体系下都能充分遍历正则分布。

原生参数：Nosé–Hoover、CSVR/CSLD、Berendsen、周期rescale。