Run reactive molecular dynamics simulations in LAMMPS with the ReaxFF potential, including preparing input scripts (pair_style reaxff + fix qeq/reaxff), mapping LAMMPS atom types to elements via pair_coeff, choosing ensembles (NVE/NVT/NPT), and adding common ReaxFF diagnostics such as species analysis. Use when the user wants LAMMPS+ReaxFF workflows or needs a working, annotated `input.lammps` template.
Use this skill when the user wants to run molecular dynamics in LAMMPS with a ReaxFF force field, prepare or explain an input.lammps file, and set up charge equilibration (QEq) correctly.
ffield.reax.*). Do not guess which file is appropriate.
potentials/ffield.reax.* at https://github.com/lammps/lammps/tree/develop/potentials).data.system) and the atom type → element mapping needed by pair_coeff.atom_style charge or atom_style full, and ensure charges are initialized either from the data file (with a charge column compatible with the chosen ) or via explicit commands (e.g. or equal-style variables). Do rely on as a substitute for a real charge field used by ReaxFF/QEq.atom_stylesetfix property/atom qfix qeq/reaxff, unless the user explicitly requests otherwise.lmp -h output before execution.Ask only for what is missing:
ffield.reax...)pair_coeff * * ffield ...)uv is available)Use:
uvx --from 'lammps[mpi]' lmp -in input.lammps
Notes:
error while loading shared libraries: libmpi.so..., you likely installed an MPI-linked lmp without MPI runtime libraries. Prefer uvx --from 'lammps[mpi]' ... (bundles MPI runtime), or load/install MPICH/OpenMPI via system packages/conda/HPC module.Do not invent the executable. Ask which command should be used, e.g.:
lmp -in input.lammpsmpirun -np 32 lmp_mpi -in input.lammpssrun lmp -in input.lammpsSee also assets/input.reaxff.nvt.lammps.
# --------- user knobs ---------
variable NSTEPS equal 200000
variable THERMO equal 200
variable DUMP equal 1000
variable TEMP equal 300.0
variable TAU_T equal 100.0
# Timestep (fs for units real). For high-T / reactive runs, 0.1 fs is often safer.
variable DT equal 0.25
# QEq parameters
variable QEQ_EVERY equal 1
variable QEQ_TOL equal 1.0e-6
variable QEQ_CUTLO equal 0.0
variable QEQ_CUTHI equal 10.0
units real
boundary p p p
atom_style charge
read_data data.system
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
# ReaxFF potential
pair_style reaxff NULL
pair_coeff * * ffield.reax C H O
# Charge equilibration (required for most ReaxFF parameterizations)
fix fqeq all qeq/reaxff ${QEQ_EVERY} ${QEQ_CUTLO} ${QEQ_CUTHI} ${QEQ_TOL} reaxff
# (`reaxff` here means QEq parameters are extracted from the ReaxFF force field file.)
# Thermo and trajectory
thermo_style custom step temp pe ke etotal press vol density
thermo ${THERMO}
dump 1 all custom ${DUMP} traj.lammpstrj id type q x y z
# Dynamics
velocity all create ${TEMP} 12345 mom yes rot yes dist gaussian
fix fnvt all nvt temp ${TEMP} ${TEMP} ${TAU_T}
timestep ${DT}
run ${NSTEPS}
units real is a common choice for ReaxFF (time in fs). Many published ReaxFF workflows use real, but the correct choice depends on the parameterization and your conventions.0.25 fs may be fine for moderate temperatures, but for high-temperature ReaxFF (especially with H present) it is common to reduce to 0.1 fs (or even 0.05 fs if needed). A quick sanity check is a short NVE segment to verify total-energy drift before running long NVT/NPT.atom_style charge is used because ReaxFF and QEq require per-atom charges.pair_style reaxff NULL uses default ReaxFF control settings. If you have a ReaxFF control file, replace NULL with its filename.pair_coeff * * ffield.reax C H O:
fix qeq/reaxff ... reaxff uses QEq parameters extracted from the ReaxFF force field file.etotal drift is reasonable (and that the run does not blow up).Example (units real):
reset_timestep 0
unfix fnvt
fix fnve all nve
# high-T ReaxFF often needs a smaller timestep
# (common choices: 0.1 fs; if needed 0.05 fs)
timestep 0.1
run 2000
maxiter (see LAMMPS fix qeq/reaxff).If the user wants reaction product tracking, add fix reaxff/species (see references/reaxff-workflow.md). This writes time series counts of detected molecular species using bond-order cutoffs.
After a run, report at least:
log.lammps)