LAMMPS键长键角约束实操

几何约束用于保持选定键长、键角；把弹簧力常数增大只是另一种带柔性势的模型，不能替代对约束算法的正确设置。本文在鸢算科技集群实际运行SHAKE、RATTLE与Kokkos SHAKE，给出拓扑、输入、fix顺序和位置/速度残差的读取方法。

1. 下载并运行#

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

mkdir constraints-demo
cd constraints-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. 三粒子教学模型与拓扑#

中心粒子质量16，两端质量1，两个键长均为1，夹角104.52°；输入包含两条bond和一条angle。虽然几何类似三点水，这里没有电荷和水的真实非键力场，使用LJ约化单位并令非键ε=0，不能用于预测水的结构、介电或扩散性质。

初速度对应质心平移和整体转动。键/角势提供约束的平衡几何，SHAKE/RATTLE选择这些拓扑类型后保持相应长度与角度。

inputs/rattle-cpu/model.data：

LAMMPS teaching model; LJ reduced units

3 atoms
2 bonds
1 angles

2 atom types
1 bond types
1 angle types

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

Masses

1 16
2 1

Atoms # molecular

1 1 1 0 0 0
2 1 2 1 0 0
3 1 2 -0.2507179360729566 0.9680601822879179 0

Velocities

1 0.1107562242476435 0.04167464373414397 0
2 0.1107562242476435 0.241674643734144 0
3 -0.08285581220994007 -0.008468943480447355 0

Bonds

1 1 1 2
2 1 1 3

Angles

1 1 2 1 3

3. 完整RATTLE输入：注意fix顺序#

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
bond_style harmonic
bond_coeff 1 100.0 1.0
angle_style harmonic
angle_coeff 1 100.0 104.52
fix integrator all nve
fix constraints all rattle 1e-10 100 0 b 1 a 1
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

fix integrator all nve 完成速度Verlet积分，fix constraints all rattle ... 放在时间积分fix之后。RATTLE要在合适的积分阶段修正位置和速度，不能随意交换两行。

1e-10 是约束求解容差，100 是单步最大迭代次数，0 表示不周期打印约束统计；b 1 a 1 选择bond type 1与angle type 1。它不是把第1个原子固定，也不是将所有原子速度归零。

在已有计算资源中单独运行：

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

4. SHAKE与Kokkos输入#

CPU SHAKE的完整输入在 inputs/shake-cpu/in.lammps，约束行为 fix constraints all shake 1e-10 100 0 b 1 a 1。Kokkos版本使用 shake/kk、nve/kk、molecular/kk、harmonic/kk 和 lj/cut/kk：

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

SHAKE主要满足位置几何约束；RATTLE额外显式满足键方向相对速度约束 r_ij · v_ij = 0。不能仅凭坐标看起来相同，就认为速度处理也完全相同。

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'] == 'constraint':
            print(r['name'], r['max_geometry_error'], r['max_angle_error_deg'], r['max_bond_velocity_dot'])
PYCODE

三个条件均为3000步、dt=0.001，最大键长偏差约6.67×10⁻¹⁶、最大角度偏差约4.27×10⁻¹⁴度。RATTLE的最大 |r_ij·v_ij| 约9.72×10⁻¹⁷，SHAKE在本例约2.14×10⁻¹³，Kokkos SHAKE约2.31×10⁻¹³。这是特定自由转动控制模型的实际数值，不能认为所有SHAKE轨迹的速度残差都会如此小。

analyze.py 从未包裹坐标和原生速度重算几何、动能、总动量与角动量，并与LAMMPS原生热力学表对照。查看 work/analysis/rattle-cpu-timeseries.csv 可逐时读取；XYZ仅作为结构与速度的可视化副本。

迁移到真实分子时，先检查bond/angle类型、质量、力场平衡几何和适用的约束簇，再测试时间步；不要对同一自由度重复添加刚体与SHAKE。本文只演示MD几何约束，最小化中同名fix可能采用强谐振近似，不能照搬为完全相同的算法。

原生说明：SHAKE与RATTLE。