LAMMPS ReaxFF物种与键级统计

ReaxFF不要求事先为每条O–H写入永久拓扑键。若要知道轨迹中有哪些分子，需要读取连续键级并选择识别阈值。本例在四水分子短轨迹上同时输出原生物种、键级邻接和完整精度原子电荷，演示如何相互检查。

下载三篇共用的完整输入、反应模板、原生科学输出和绘图包。按共同运行步骤提交一次即可。对应文件位于 inputs/reaxff-species/；更详细的ReaxFF/QEq和COMB3差别见既有电荷教程。

1. 输入和原生统计命令#

12原子、30 Å周期盒，4个水分子；势参数是配套CHO文件，类型1为O、类型2为H。初速度100 K、NVE、每步0.05 fs、20步共1 fs。这个长度只够展示文件和统计流程，不用于推断水的平衡性质或反应速率。

cat inputs/reaxff-species/input.lammps
cat inputs/reaxff-species/water.data

完整输入：

units real
atom_style charge
boundary p p p
read_data water.data
pair_style reaxff NULL
pair_coeff * * ffield.reax.cho O H
fix charges all qeq/reaxff 1 0.0 10.0 1e-10 reaxff maxiter 500
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
fix species all reaxff/species 1 1 1 species.out element O H cutoff * * 0.3
fix orders all reaxff/bonds 1 bondorders.out
compute qa all property/atom q
compute totalq all reduce sum c_qa
velocity all create 100 46631 mom yes rot no dist gaussian
fix integrator all nve
timestep 0.05
thermo 1
thermo_style custom step temp pe ke etotal c_totalq
thermo_modify norm no format float %.17g
dump atoms all custom 1 atoms.dump id type q x y z fx fy fz
dump_modify atoms sort id format float %.17g
run 20
write_data final.data nocoeff

qeq/reaxff每步解电荷；reaxff/species 1 1 1每步采样一次、每次平均1个样本、每步写出物种；reaxff/bonds 1每步保存键级邻接。cutoff * * 0.3指定物种分组阈值，element O H与类型顺序对应。生产任务可以使用多帧平均减小阈值附近抖动，同时要遵守采样/平均/输出步数约束，并注意原生统计对邻居重建频率的要求。

2. 先直接查看输出#

head -n 6 inputs/reaxff-species/species.out
head -n 20 inputs/reaxff-species/bondorders.out

本次原生物种输出开头：

#  Timestep    No_Moles    No_Specs         H2O
          0           4           1           4
#  Timestep    No_Moles    No_Specs         H2O
          1           4           1           4
#  Timestep    No_Moles    No_Specs         H2O
          2           4           1           4

No_Moles在这里表示识别出的分子个数，不是物质的量mol；No_Specs是不同物种数。全21帧均为4个H₂O、1种物种，没有在本段轨迹中发生解离。元素排序由输出约定决定，分析时不要只凭固定列位置读取所有任务。

3. 用键级图独立检查物种#

bondorders.out每行包含原子ID、类型、邻接数、邻居ID、对应键级及电荷等。脚本逐帧将大于0.3的无向连接建图，求连通分量，再按每个分量的真实元素数生成化学式，与原生 species.out比较；本次21帧全部一致。

注意该键级文件只列出高于粗阈值0.3的邻居，且键级是有限小数的文本。本例物种阈值也为0.3，O–H键级远离阈值，能清楚识别；不能从这个文件凭空重建更低阈值下未输出的弱键，也不能把被舍入的边界数值当成精确的键开关。

4. 总电荷为什么不能用三位小数相加#

原生键级文件的电荷仅打印约三位小数，显示的每个O约−0.525、每个H约+0.263，直接相加会产生舍入残差。判断系统是否守恒，要使用 atoms.dump 的17位电荷或thermo中 c_totalq；本次完整精度净电荷最大绝对值约 8.33e-16 e。

module load sisl/0.16.4
python analyze.py
python plot.py

这只重新读取结果并画图，不重跑MD。脚本检查全部原子ID、21帧时间、双向键级、物种、完整精度电荷和有限数值。部分电荷会随几何变化，分子总电荷并不由全体系中性自动保证；本例整体中性也不表示一般反应中每个物种都中性。

需要研究实际反应物种演化时，应延长轨迹、验证势参数和电荷求解，检查不同识别阈值/时间平均下的稳定性，并保留短寿命中间体的处理规则。本篇未提供真实化学反应速率、平衡组成或电荷转移的定量结论。官方说明：物种统计、键级输出。