LAMMPS刚体平移与转动实操

把一组原子作为刚体时，内部相对几何保持不变，体系仍可以整体平移与转动。本文用4粒子的自由方形转子，实际对照 rigid、rigid/nve、rigid/small、rigid/nve/small，以及两个Kokkos版本，学习如何定义刚体、避免重复积分并读取轨迹。

1. 下载与提交#

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

mkdir rigid-demo
cd rigid-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。

2. 准备一个有解析参照的刚体#

4粒子质量均为1，初始位置组成边长1的正方形；质心初速度 (0.1,0.05,0)，绕z轴角速度0.2。相互作用ε设为零，没有外力或外力矩，因此质心应匀速运动，内部几何与角动量保持不变。它是算法控制模型，不是某种材料的真实受力或GPU性能测量。

inputs/rigid-cpu/model.data 是完整原生数据文件：

LAMMPS teaching model; LJ reduced units

4 atoms
0 bonds
0 angles

1 atom types

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

Masses

1 1

Atoms # molecular

1 1 1 -0.5 -0.5 0
2 1 1 0.5 -0.5 0
3 1 1 0.5 0.5 0
4 1 1 -0.5 0.5 0

Velocities

1 0.2 -0.05 0
2 0.2 0.15 0
3 0 0.15 0
4 0 -0.05 0

每个原子的molecule-ID都为1；后续 /small 输入按molecule识别同一个刚体。真实多分子体系应给各个刚体正确且不同的molecule-ID，不能把所有分子误并成一个整体。

3. 完整输入与独立运行#

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.0 1.0
neighbor .3 bin
neigh_modify every 1 delay 0 check yes
fix integrator all rigid single
timestep 0.001
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 3000
print "COMPLETE" file complete.dat screen no

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

fix ... rigid single 本身负责该组原子的积分，不要再对相同原子加 fix nve，否则会重复积分。刚体以整体平移、转动的自由度运动，控温时也应采用与刚体自由度匹配的设置；这里不加恒温器，观察能量守恒。

rigid/nve 使用另一种刚体NVE积分方案；rigid/small 和 rigid/nve/small 使用适合大量小刚体的数据分配方式，不表示把刚体拆成可变形原子。本文只有一个小刚体，不能据此给这些算法做速度排名。

4. Kokkos版本怎么调用#

包内 inputs/rigid-small-kk/in.lammps 使用 atom_style molecular/kk、pair_style lj/cut/kk、run_style verlet/kk，以及 fix integrator all rigid/small/kk molecule。半邻居表与Newton on配置和其他Kokkos例子一致：

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

第二个Kokkos输入为 rigid/nve/small/kk。这里只对照这两个命名版本，不把其他刚体后缀的可用性混为同一个结论。

5. 读取实际守恒量与轨迹#

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'] == 'rigid':
            print(r['style'], r['max_geometry_error'], r['max_energy_drift'])
PYCODE

初始总动能0.065：质心平动贡献0.025、转动贡献0.04，角动量Lz=0.4。6个条件跑3000步、dt=0.001，六条内部距离的最大偏差约1.19×10⁻¹³，最大能量漂移约3.38×10⁻¹⁴。数据工具另用匀速平移和解析转动公式核对坐标。

这些很小的误差来自无外力、对称转子的控制条件，不保证任意复杂刚体、外部相互作用、温压控制或更大时间步都同样准确。需要约束特定键角时看SHAKE/RATTLE，多时间步看RESPA。

原生选项：fix rigid。