LAMMPS ReaxFF与COMB3电荷实操

下载完整输入、势文件、实际轨迹和绘图脚本，解压进入目录，先查看已有结果：

bash run.sh

输出 analysis/ 和 figures/。重新计算时先核对脚本里的分区适合自己的账户，再提交：

sbatch submit.slurm

这套脚本申请1个P100和1个CPU，依次运行普通计算路径、适用的GPU/Kokkos路径，再分析自己的结果。已有 work/analysis/figures 时请另解压一份，不要覆盖。直接执行 bash run.sh compute 或文中的 lmp 命令时，必须已经处于分配好的计算资源中。

实际软件为 LAMMPS 30Sep2026，GCC 13.3、OpenMPI 5.0.10、CUDA 12.8；GPU包采用Mixed Precision，Kokkos 5.2.1采用双精度。程序链接现有PLUMED 2.9.4，但本文不使用偏置。数据后处理使用ASE 3.29环境中的NumPy与Matplotlib。普通路径的数值对照在同一P100计算资源的CPU上完成，不是Intel Xeon 8173M性能测试。

1. 先区别两个模型与单位#

本例使用30 Å周期盒中的四个水分子，共12原子。改变第一个分子的O–H1输入距离，其他输入坐标保持相同，并在每个条件下重新平衡电荷。全胞初始与计算后总电荷为零；不要求每个分子的部分电荷分别为零。

• ReaxFF采用原配套 ffield.reax.cho，units real：能量kcal/mol、时间fs。
• COMB3采用原配套 ffield.comb3，units metal：能量eV、时间ps。

下面图表把能量统一换成eV，1 kcal/mol约为0.0433641 eV；原始日志继续保留各自的原生单位。两套势是独立参数化，绝对能量零点不能直接互比。这些曲线是有限周期盒中的全胞能量变化，不是严格孤立水分子的解离能或DFT反应路径。

2. ReaxFF与QEq完整输入#

下面是实际Kokkos输入 inputs/reaxff-force-0-kk/in.lammps：

package kokkos neigh half
newton on
units real
atom_style charge/kk
boundary p p p
read_data water.data
pair_style reaxff/kk NULL
pair_coeff * * potentials/ffield.reax.cho O H
fix charges all qeq/reaxff/kk 1 0.0 10.0 1e-10 reaxff maxiter 500
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
compute pea all pe/atom
compute pesum all reduce sum c_pea
compute qatom all property/atom q
compute qsum all reduce sum c_qatom
timestep 0.1
thermo 100
thermo_style custom step time temp pe ke etotal press vol c_pesum c_qsum
thermo_modify norm no format float %.16g
run_style verlet/kk
run 0
write_dump all custom final.dump id type x y z q vx vy vz fx fy fz c_pea modify sort id format float %.16g
print "$(pe:%.16g) $(vol:%.16g) $(c_pesum:%.16g) $(c_qsum:%.16g)" file energy.dat screen no

配套 water.data 原子段格式为 id type q x y z，type 1为O，type 2为H；完整12原子输入在下载包中。pair_coeff 的 O H 必须与类型顺序一致。这里用 atom_style charge 对应变体，不需要先为每条O–H定义永久拓扑键；反应势的相互作用来自自己的键级形式。

fix qeq/reaxff 每步更新电荷，给出距离范围、求解阈值和参数来源。它与COMB的 fix qeq/comb 不是同一个命令。Kokkos输入配套 qeq/reaxff/kk 与 run_style verlet/kk，原生启动方式：

module purge
module load lammps/30Sep2026-p100-kokkos
cd inputs/reaxff-force-0-kk
cp -r ../../potentials .
lmp -k on g 1 -in in.lammps
cat energy.dat

3. COMB3与初始化#

inputs/comb3-force-0-cpu/in.lammps：

newton on
units metal
atom_style charge
boundary p p p
read_data water.data
pair_style comb3 polar_off
pair_coeff * * potentials/ffield.comb3 O H
fix charges all qeq/comb 1 1e-8
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
compute pea all pe/atom
compute pesum all reduce sum c_pea
compute qatom all property/atom q
compute qsum all reduce sum c_qatom
timestep 5e-05
thermo 100
thermo_style custom step time temp pe ke etotal press vol c_pesum c_qsum
thermo_modify norm no format float %.16g
run 0
run 0
print "MEASUREMENT_BEGIN"
run 0
write_dump all custom final.dump id type x y z q vx vy vz fx fy fz c_pea modify sort id format float %.16g
print "$(pe:%.16g) $(vol:%.16g) $(c_pesum:%.16g) $(c_qsum:%.16g)" file energy.dat screen no

COMB3使用自己的势文件和 polar_off 选择；不要与第二代COMB的 ffield.comb 互换。这里使用普通计算路径，没有COMB3的GPU加速声明。

本例在正式记录前执行两次 run 0。实际初始首个评价的势能约-33.6354486 eV，后续固定坐标评价稳定在约-38.2341155 eV；不能把这种初始化偏移当成分子运动释放的物理能量。第三次评价或短MD开始后再记录数据。两次预评价是本例的操作方式，其他体系仍需检查自己的电荷与能量稳定，并完成适当的结构优化和平衡。

4. 改变一条O–H距离后得到什么#

实际全胞能量相对各自0.95 Å输入的变化，以及第一分子的部分电荷：

部分电荷是该经典势及电荷求解方法下的量，不是形式氧化态、Bader电荷或任意电子结构分析的同义词。伸长一条键也可能改变与其他原子及周期镜像的相互作用，不应把这张有限模型曲线直接报告成实验键解离能。

5. 用能量差分核对力，先说明电荷路径#

本例ReaxFF正常QEq参考的H1-x方向原生力为约0.60863234 eV/Å。每个偏移点都重新平衡电荷时，能量差分给出约0.60199 eV/Å；它们在这个实际设置中并不相同。

为看清差别，脚本从本次计算自己的参考电荷派生两类额外输入：固定参考电荷时作位移差分，以及在参考几何上只沿实际电荷变化方向取能量差。结果为：

固定电荷差分与原生力一致；后一个电荷变化项能解释两条路径的差异。这里不把重新平衡后的能量曲线直接当作固定电荷力的高精度差分检验，也不宣称这个具体势/QEq设置给出严格等价的变分力。

诊断输入用 checkqeq no 和自己刚得到的电荷，仅用于这项比较；正常反应MD仍保留对应QEq，不要通过关闭电荷平衡来掩盖差异。generate_reax_diagnostics.py 自动读本次的 final.dump，不会把其他计算预先保存的电荷套进自己的体系。

6. 查看输出与继续正式计算#

bash run.sh
cat analysis/forces.csv
cat analysis/backend-errors.csv

analysis/conditions.csv列出全胞能量、电荷和几何，RESULTS.json还有两条电荷路径与短MD的完整数值。ReaxFF另以0.1 fs、COMB3以0.00005 ps做1000步短NVE，用于检查实际轨迹和能量变化。正式反应、输运或长时间统计需要继续收敛时间步、模型、边界和电荷设置，不从这些短轨迹给出速率或热导结论。

官方说明：ReaxFF、COMB/COMB3、COMB电荷平衡。其他经典势操作见Si共价势与AIREBO碳片。