LAMMPS + ReaxFF

When to Use

• User needs reactive molecular dynamics (bonds can break and form)
• User wants to study combustion, pyrolysis, oxidation, or corrosion
• User needs to model chemical reactions at the atomic scale without QM
• User has a ReaxFF parameter file for the relevant elements

Prerequisites

1. LAMMPS compiled with REAXFF package (lmp -h | grep REAXFF)
2. ReaxFF parameter file (ffield.reax.xxx) for the element combination
3. Initial structure — can be built with Packmol (see data/packmol/SKILL.md)

Workflow Steps

1. Create workflow

catgo_workflow_engine(action="create", params={"name": "ReaxFF combustion - methane/O2"})

2. Add ReaxFF MD task

catgo_workflow_engine(action="add_task", params={
  "workflow_id": "wf_xxx",
  "task_type": "shell",
  "name": "reaxff_md",
  "command": "lmp -in lammps.in > lammps.log 2>&1",
  "input_files": {
    "lammps.in": "<input script>",
    "ffield.reax.CHO": "<reaxff params>"
  },
  "system_name": "CH4_combustion"
})

3. Analyze with ReacNetGenerator (optional)

After MD completes, use analysis/reacnetgen/SKILL.md to extract reaction networks.

LAMMPS Input Template — ReaxFF NVT

units           real
boundary        p p p
atom_style      charge

read_data       structure.lmp

pair_style      reaxff NULL
pair_coeff      * * ffield.reax.CHO C H O

neighbor        2.0 bin
neigh_modify    every 10 delay 0 check yes

# Charge equilibration
fix             qeq all qeq/reaxff 1 0.0 10.0 1.0e-6 reaxff

# Velocities
velocity        all create 2000.0 12345 dist gaussian

# NVT
fix             1 all nvt temp 2000.0 2000.0 100.0

# Timestep (fs in real units)
timestep        0.25

# Output
thermo          100
thermo_style    custom step temp pe ke etotal press vol density

dump            1 all custom 100 traj.lammpstrj id type q x y z
dump_modify     1 sort id element C H O

# Bond order output (for reaction analysis)
fix             bonds all reaxff/bonds 1000 bonds.reax

run             500000

Important ReaxFF-Specific Settings

Charge Equilibration (MANDATORY)

ReaxFF requires dynamic charge equilibration at every step:

fix  qeq all qeq/reaxff 1 0.0 10.0 1.0e-6 reaxff

Never omit this fix. Without it, charges are not updated and results are physically meaningless.

Small Timestep

ReaxFF requires small timesteps due to reactive dynamics:

• 0.25 fs — safe for most systems (recommended)
• 0.5 fs — acceptable for heavy-atom systems without H
• 1.0 fs — only for non-reactive equilibration phases

Bond Order Output

For reaction analysis, write bond information:

fix  bonds all reaxff/bonds 1000 bonds.reax

Feed bonds.reax to ReacNetGenerator for automated reaction network extraction.

Parameter File Sources

File	Elements	Reference
ffield.reax.CHO	C, H, O	Chenoweth et al. (2008) — combustion
ffield.reax.Fe_O_C_H	Fe, O, C, H	Aryanpour et al. (2010) — iron oxidation
ffield.reax.AB	A, B (various)	Check doi.org/10.1021/acs.jctc.* for specific systems

Always cite the original ReaxFF parameterization paper.

Parameter Guidance

Parameter	Typical value	Notes
timestep	0.25 fs	Real units (fs); must be small for reactive dynamics
NVT temp damp	100 fs	100-200 fs typical in real units
qeq tolerance	1.0e-6	Charge equilibration convergence
Temperature	1500-3000 K	Reactive events need high T to overcome barriers
Run length	500K-5M steps	125 ps - 1.25 ns at 0.25 fs timestep
Dump freq	100-500	For trajectory and reaction analysis

Common Pitfalls

1. Missing qeq/reaxff fix — charges are not static in ReaxFF. Omitting charge equilibration gives garbage results.
2. Wrong units — ReaxFF uses units real (kcal/mol, Ang, fs). Never use units metal.
3. Timestep too large — 0.25 fs is the safe maximum. Larger timesteps cause energy conservation failure.
4. Wrong element order — pair_coeff * * ffield.reax.CHO C H O must match atom types in data file.
5. Parameter file incompatibility — ReaxFF parameters are NOT transferable between element combinations. Use only validated parameter sets.
6. System too small — ReaxFF MD of reactions needs at least 500-1000 atoms for meaningful statistics.
7. Temperature too low — reactions have barriers. Below 1000 K, you may see no reactive events in affordable simulation time.