ReaxFF vs Classical Force Fields: A Comprehensive Guide for Combustion Chemistry Simulations

Abstract

This article provides a detailed comparative analysis of the Reactive Force Field (ReaxFF) and classical force fields for modeling combustion chemistry. Tailored for researchers and scientists, it explores the foundational principles of both methods, with ReaxFF's bond-order formalism enabling the simulation of reactive events. It covers key methodological approaches and applications in pyrolysis and combustion of aerospace fuels and energetic materials. The content further addresses common challenges, performance benchmarks, and strategies for force field selection and optimization. By synthesizing insights from validation and comparative studies, this guide aims to equip professionals with the knowledge to effectively apply and troubleshoot these computational tools in energy and combustion research.

Understanding the Core Principles: From Rigid Bonds to Reactive Chemistry

In combustion chemistry research, molecular dynamics (MD) simulations are indispensable for probing reaction mechanisms, intermediate species, and energy release processes at the atomic scale. The accuracy of these simulations hinges entirely on the underlying molecular forcefield—a mathematical model that describes the potential energy of a system as a function of its atomic coordinates [1]. Force fields calculate the forces acting on atoms, enabling the simulation of their motion according to Newton's laws [1]. Among these models, classical harmonic force fields have been widely successful for simulating biomolecules, polymers, and molecular crystals. However, they possess a fundamental limitation: they require predefined atomic connectivity, preventing the simulation of chemical reactions where bonds break and form [2] [3].

This constraint becomes particularly significant in combustion chemistry, where complex reaction networks involve continuous bond dissociation and formation. Classical force fields treat molecules as having fixed bond topologies, with energy calculated as a sum of bonded terms (bonds, angles, dihedrals) and non-bonded terms (electrostatics, van der Waals) [4] [1]. This formulation is inadequate for modeling reactive events, creating a persistent gap in our ability to simulate combustion processes across relevant time and length scales using traditional molecular mechanics approaches [2] [3].

The Classical Force Field Paradigm and Its Functional Forms

Energy Formulation of Additive Force Fields

Classical force fields, often termed "additive" or "non-reactive" force fields, employ a potential energy function that separates contributions into bonded and non-bonded terms. A typical expression is shown below [4]:

U(r) = ∑bondskb(b - b0)2 + ∑angleskθ(θ - θ0)2 + ∑dihedralskχ(1 + cos(nχ - δ)) + ∑vdW,i≠jεij[(Rmin,ij/rij)12 - 2(Rmin,ij/rij)6] + ∑elec,i≠j(qiqj)/(4πε0rij)

where the symbols represent:

• b₀ and θ₀: Equilibrium bond length and valence angle
• kb, kθ, k_χ: Force constants for bonds, angles, and dihedrals
• n: Dihedral multiplicity; δ: Dihedral angle phase
• εij: van der Waals well depth; Rmin,ij: Distance at minimum energy
• qi, qj: Partial atomic charges; r_ij: Interatomic distance

Table 1: Key Components of Classical Force Field Energy Functions

Energy Component	Functional Form	Atoms Involved	Role in Combustion Chemistry Limitations
Bond Stretching	Harmonic: kb(b - b0)2	2 (bonded)	Prevents bond dissociation at high temperatures
Angle Bending	Harmonic: kθ(θ - θ0)2	3 (bonded)	Restricts molecular flexibility during pyrolysis
Dihedral Torsion	Cosine series: kχ[1 + cos(nχ - δ)]	4 (bonded)	Limits conformational sampling in reactive intermediates
van der Waals	Lennard-Jones 6-12	All non-bonded	Models dispersion/repulsion but without reactivity
Electrostatics	Coulomb's law: (qiqj)/(4πε0rij)	All non-bonded	Fixed charge distribution ignores electronic changes

The Fixed Topology Requirement and Its Consequences

The fundamental limitation of classical force fields stems from their requirement of predefined atomic connectivity. The chemical structure—which atoms are bonded to which others—must be specified before simulation begins and remains fixed throughout [3] [5]. This fixed topology approach has two critical consequences for combustion chemistry research:

• Inability to Simulate Reaction Chemistry: Combustion processes involve complex networks of bond-breaking and bond-forming events. The fixed connectivity of classical force fields precludes simulations of these essential chemical transformations [2] [3].

• Environmental Transferability Issues: Classical force fields use atom types that are more specific than element types alone (e.g., distinguishing sp² vs. sp³ carbon) to achieve transferability across molecular environments [1]. However, this approach cannot accommodate changes in hybridization or bonding that occur during reactions.

The following diagram illustrates the fundamental architecture of classical force fields and where the connectivity constraint arises:

Reactive Force Fields: Overcoming the Connectivity Limitation

ReaxFF: A Bond-Order Approach to Reactivity

The Reactive Force Field (ReaxFF) methodology was developed specifically to bridge the gap between quantum mechanical accuracy and classical force field efficiency for simulating reactive systems [2]. ReaxFF eliminates the fixed connectivity constraint through a bond-order formalism that dynamically describes chemical bonding based on interatomic distances [2]. This approach implicitly captures electronic effects governing bond formation and dissociation without expensive quantum calculations.

In ReaxFF, the bond order between atoms i and j is empirically calculated from interatomic distance using the formula [2]:

BOij = BOijσ + BOijπ + BOijππ = exp[pbo1(rij/roσ)^pbo2] + exp[pbo3(rij/roπ)^pbo4] + exp[pbo5(rij/roππ)^pbo6]

where:

• BO_ij: Total bond order between atoms i and j
• BOij^σ, BOij^π, BO_ij^ππ: Sigma, pi, and double pi bond order components
• r_ij: Interatomic distance
• ro^σ, ro^π, r_o^ππ: Equilibrium bond lengths for different bond types
• p_bo*: Empirical parameters

This bond-order calculation is continuous and contains no discontinuities through transitions between different bond types, creating a differentiable potential energy surface required for calculating interatomic forces [2].

Comprehensive Energy Formulation in ReaxFF

The ReaxFF potential energy includes multiple contributions that depend on the dynamic bond order [2]:

Esystem = Ebond + Eover + Eangle + Etors + EvdWaals + ECoulomb + ESpecific

where:

• E_bond: Energy associated with bond formation based on bond order
• E_over: Penalty energy preventing atomic over-coordination
• E_angle: Three-body valence angle strain
• E_tors: Four-body torsional angle strain
• E_vdWaals: van der Waals dispersive interactions
• E_Coulomb: Electrostatic interactions
• E_Specific: System-specific terms (lone-pair, conjugation, etc.)

Table 2: Comparison of Classical Force Fields vs. ReaxFF for Combustion Chemistry

Characteristic	Classical Force Fields	ReaxFF
Atomic Connectivity	Fixed and predefined	Dynamic, based on interatomic distances
Bond Representation	Harmonic springs	Bond-order formalism
Reaction Capability	None without manual intervention	Full bond breaking/formation
Parameterization	Fitted to molecular properties	Trained against QM data including reaction pathways
Computational Cost	Lower (baseline)	~30x higher than classical FFs [3]
Combustion Applications	Limited to non-reactive molecular dynamics	Hydrocarbon oxidation, pyrolysis, catalyst interactions
Electrostatics	Fixed partial charges	Charge equilibration at each step
Element Transferability	Limited by atom type definitions	Single parameter set for elements across phases

Methodologies: Parameterization and Validation Approaches

Force Field Parameterization Strategies

The development of robust force fields for combustion research requires careful parameterization against reliable reference data:

Classical Force Field Parameterization: Classical force fields are typically parameterized using a combination of experimental data (structural, vibrational, thermodynamic) and quantum mechanical calculations [4] [6]. The process involves:

• Target Data Selection: Crystal structures, gas-phase geometries, vibrational frequencies, sublimation enthalpies, and liquid properties
• Parameter Optimization: Iterative adjustment of parameters to minimize differences between calculated and target properties
• Validation: Testing against experimental data not included in the parameterization

ReaxFF Parameterization: ReaxFF employs a more comprehensive training approach against quantum mechanical data [2]:

• Training Set Construction: Includes QM data for bond dissociation curves, angle bending, reaction barriers, and energies of reaction
• Parameter Optimization: Uses the Metropolis Monte Carlo approach to optimize parameters to reproduce QM training data
• Validation: Testing on properties not included in training set to ensure transferability

Experimental Protocol for Combustion Chemistry Validation

To validate force field performance for combustion systems, researchers employ several computational experimental protocols:

Reactive Molecular Dynamics Simulations:

• System Preparation: Construct simulation cells containing fuel and oxidizer molecules at relevant densities
• Equilibration: Run initial MD simulation to equilibrate system at target temperature and pressure
• Production MD: Perform extended MD simulations while monitoring reaction events
• Analysis: Track reaction products, intermediate species, and temporal evolution of chemical composition

Quantitative Validation Metrics:

• Ignition Delay Times: Compare simulated and experimental ignition delays
• Species Profiles: Validate intermediate and product concentrations against experimental measurements
• Reaction Pathways: Analyze dominant reaction mechanisms and compare to QM calculations and literature
• Thermodynamic Properties: Validate energies of reaction and formation

The following diagram illustrates the contrasting computational workflows for classical and reactive molecular dynamics:

Table 3: Research Reagent Solutions for Reactive Molecular Dynamics

Tool/Resource	Function	Relevance to Combustion Chemistry
ReaxFF Parameter Sets	Element-specific parameters for reactive simulations	Enable simulation of hydrocarbon oxidation mechanisms
QM Reference Data	High-level quantum calculations for parameterization	Provide training data for bond dissociation and reaction barriers
LAMMPS	Open-source MD package with ReaxFF implementation	Primary simulation engine for large-scale reactive MD
PuReMD	Purdue Reactive Molecular Dynamics code	Optimized for efficient ReaxFF simulations
CHARMM, AMBER, GAFF	Classical force fields for non-reactive components	Model reactant molecules before ignition
VMD/OVITO	Visualization and analysis tools	Track reaction progress and species evolution
AIMD	Ab Initio Molecular Dynamics	Generate reference data for validation of ReaxFF results

The predefined atomic connectivity required by classical force fields represents a fundamental limitation for combustion chemistry research, where bond-breaking and bond-forming events are central to the processes of interest. ReaxFF addresses this limitation through its bond-order formalism, enabling dynamic connectivity that evolves during simulations. While computationally more demanding than classical force fields, ReaxFF provides a crucial bridge between quantum accuracy and classical efficiency for simulating complex reactive systems.

As combustion research increasingly focuses on optimizing efficiency and reducing emissions, the ability to simulate detailed reaction mechanisms at realistic time and length scales becomes ever more critical. Reactive force fields like ReaxFF, especially when combined with emerging machine learning approaches [7], offer a powerful toolkit for unraveling the complex chemistry of combustion processes and designing next-generation combustion systems.

In the field of combustion chemistry research, computational scientists face a persistent challenge: the trade-off between quantum mechanical (QM) accuracy and classical molecular dynamics scale. Quantum mechanical methods, while providing valuable electronic-level guidance, are often too computationally intensive for simulating the full dynamic evolution of a system, particularly at the scales relevant for combustion processes [2]. Alternatively, traditional empirical interatomic potentials require significantly fewer computational resources, enabling larger-scale and longer-timescale simulations, but they typically require predefined connectivity between atoms, which prevents the simulation of reactive events where bonds break and form [2] [8]. This fundamental limitation has historically restricted the application of molecular dynamics to equilibrium states rather than reactive processes central to combustion.

The Reactive Force Field (ReaxFF) methodology, developed by Adri van Duin, William A. Goddard, III, and co-workers, represents a transformative approach that bridges this gap [9]. By casting the empirical interatomic potential within a bond-order formalism, ReaxFF implicitly describes chemical bonding without expensive QM calculations, enabling the simulation of reactive events at a fraction of the computational cost [2]. This breakthrough has opened unprecedented opportunities for studying complex combustion phenomena at the atomic scale, providing insights into reaction mechanisms, intermediate species, and dynamic processes that are often difficult to observe experimentally. For combustion researchers investigating hydrocarbon oxidation, pollutant formation, or energetic material decomposition, ReaxFF offers a powerful tool that balances computational efficiency with chemical accuracy.

The Core Breakthrough: Bond-Order Formalism in ReaxFF

Fundamental Principles

The revolutionary aspect of ReaxFF lies in its departure from traditional fixed-bond approaches. Whereas classical force fields rely on explicitly defined bonds that remain constant throughout a simulation, ReaxFF eschews explicit bonds in favor of continuously calculated bond orders derived from interatomic distances [9]. This fundamental shift in perspective allows the force field to naturally accommodate bond formation and breaking as chemical environments change during dynamic processes such as combustion reactions.

The bond order between atoms i and j is calculated empirically from interatomic distance using the formula [2]: BOij = BOij^σ + BOij^π + BOij^ππ = exp[pbo1(rij/r0^σ)^pbo2] + exp[pbo3(rij/r0^π)^pbo4] + exp[pbo5(rij/r0^ππ)^pbo6] where BO is the bond order between atoms i and j, r_ij is interatomic distance, r0 terms are equilibrium bond lengths, and pbo terms are empirical parameters [2]. This continuous function contains no discontinuities through transitions between σ, π, and ππ bond character, yielding a differentiable potential energy surface essential for calculating interatomic forces [2].

Comprehensive Energy Description

The ReaxFF potential incorporates multiple energy contributions that collectively describe the complex energy landscape of reactive systems [2] [10]: Esystem = Ebond + Eover + Eangle + Etors + EvdWaals + ECoulomb + ESpecific

Table: Energy Components in the ReaxFF Potential

Energy Term	Physical Significance	Role in Reactive Simulations
Ebond	Energy associated with bond formation/breaking	Core reactive term based on bond order
Eover	Over-coordination penalty	Prevents unrealistic atomic valences
Eangle	Three-body angle strain	Maintains molecular geometry
Etors	Four-body torsional strain	Governs conformational flexibility
EvdWaals	van der Waals interactions	Describes dispersion forces
ECoulomb	Electrostatic interactions	Handles charge-charge interactions
ESpecific	System-specific corrections	Addresses unique chemical environments

Each term contributes to accurately modeling both reactive and non-reactive interactions. For instance, the over-coordination penalty (Eover) applies a stiff energy penalty if an atom forms more bonds than allowed by its atomic valence rules (e.g., a carbon atom forming more than four bonds) [2]. The system-specific term (ESpecific) represents specialized corrections that are included only when necessary to capture properties particular to certain systems, such as lone-pair interactions, conjugation effects, hydrogen binding, and C2 corrections [2].

Figure 1: The ReaxFF Computational Workflow illustrating how interatomic distances are continuously converted into bond orders, which then determine energy and forces in molecular dynamics simulations

Comparative Analysis: ReaxFF vs. Classical Force Fields in Combustion Research

Fundamental Methodological Differences

The distinction between ReaxFF and classical force fields becomes particularly significant in combustion chemistry applications, where bond rearrangement is fundamental to the processes being studied. Classical force fields maintain fixed bonding connectivity throughout simulations, making them unsuitable for modeling chemical reactions [9] [8]. In contrast, ReaxFF's bond-order formalism allows it to automatically handle the bond formation and breaking events that characterize combustion processes, from initial fuel decomposition to final oxidation products.

Table: Comparative Analysis of Force Field Approaches for Combustion Chemistry

Feature	Classical Force Fields	ReaxFF
Bond Treatment	Fixed, explicitly defined	Dynamic, based on interatomic distances
Reaction Capability	Cannot simulate reactions	Models bond formation/breaking
Charge Treatment	Fixed or geometry-based	Dynamic (charge equilibration)
Computational Cost	Low	Moderate (10-100x classical)
System Preparation	Requires predefined connectivity	No connectivity specification needed
Transferability	System-specific	Broad across phases
Combustion Application	Limited to non-reactive processes	Full reaction mechanism elucidation

Advantages for Combustion Systems

ReaxFF's design provides several distinct advantages for combustion research. Its phase transferability means that an oxygen atom is treated with the same mathematical formalism whether it is in the gas phase as O₂, in the liquid phase within an H₂O molecule, or incorporated in a solid oxide [2]. This is particularly valuable in combustion systems where interfaces between phases play crucial roles in reaction pathways. Additionally, ReaxFF's parameterization against QM data ensures that it captures the essential chemistry of reactive processes while remaining computationally efficient enough to simulate the time and length scales relevant to combustion phenomena [9] [2].

The force field's ability to simulate reactive events without predefined reaction pathways makes it exceptionally valuable for discovering previously unknown reaction mechanisms [10]. In complex combustion processes involving multiple intermediate species and competing pathways, this unbiased approach can reveal unexpected chemistry that might be overlooked in traditional kinetic modeling based on predetermined reaction sets.

Figure 2: Comparative capabilities of Classical Force Fields versus ReaxFF for combustion research, highlighting the fundamental differences in bond treatment and resultant applications

Parameterization and Training of ReaxFF

Training Methodology

The development of an accurate ReaxFF force field requires extensive parameterization against high-quality quantum mechanical data. The process begins with constructing a comprehensive training set that covers the relevant chemical phase space, including bond and angle stretches, activation and reaction energies, equation of state, surface energies, and other essential properties [9]. Density functional theory (DFT) calculations typically serve as the source for this training data, providing a pragmatic balance between accuracy and computational feasibility [9].

Due to the complexity of the force field with its many parameters, global optimization techniques offer the best chance to obtain a parameter set that closely describes the training data [9]. Recent advances have employed Monte Carlo and evolutionary algorithms to systematically optimize ReaxFF parameters [9]. The parameterization process must carefully balance the description of various chemical environments to ensure transferability across different molecular contexts and phases.

Specialized Force Fields for Combustion Chemistry

The ReaxFF framework has been parameterized for numerous systems relevant to combustion chemistry:

• Hydrocarbon Oxidation (CHO.ff): Parameterized for simulating hydrocarbon oxidation, this force field was trained against DFT calculations on dissociation energies for various bonds containing carbon, oxygen, and hydrogen, along with angle and torsion potentials for C/H/O interactions [11].
• High Energy Materials (HE.ff, HE2.ff): Developed for simulating energetic materials like RDX and HMX, these force fields were trained against thousands of QM calculations including bond breaking/compression curves, angle and torsion bending data, and crystal data [11].
• Recent Combustion-Specific Parameterizations: The ReaxFFCHO-S22 force field represents a recent advancement specifically re-parameterized against high-level quantum mechanical calculations to capture detailed combustion behaviors of hydrocarbon fuels [10].

Table: Key ReaxFF Parameterizations for Combustion-Related Chemistry

Force Field	Elements	Combustion Applications	Reference
CHO.ff	C, H, O	Hydrocarbon oxidation	Chenoweth et al., 2008 [11]
HE.ff	C, H, O, N	Energetic materials decomposition	Zhang et al., 2009 [11]
HCONSB.ff	H, C, O, N, S, B	Coal combustion, soot formation	Castro-Marcano et al., 2012 [11]
ReaxFFCHO-S22	C, H, O	Hydrocarbon combustion mechanisms	Xiao et al., 2025 [10]

Experimental Protocols: ReaxFF Molecular Dynamics in Practice

Simulation Setup and Workflow

Implementing ReaxFF molecular dynamics simulations for combustion research follows a systematic protocol:

• System Construction: Build the initial molecular system containing reactants (e.g., fuel and oxidizer molecules) in appropriate stoichiometry and configuration. For the cyclohexane combustion study by Xiao et al. [10], researchers designed a periodic box measuring 33.34 Å on each side, populated with 12 cyclohexane and 108 O₂ molecules, representing an equivalence ratio (Φ) of 1.0 for complete combustion.

• Force Field Selection: Choose an appropriate parameter set for the chemical system under investigation. The ReaxFFCHO-S22 force field was specifically developed for hydrocarbon combustion simulations and provides accurate description of C/H/O interactions [10].

• Equilibration: Perform initial equilibration under the target conditions (temperature, density) using non-reactive MD or preliminary ReaxFF MD with restricted bonding.

• Production Simulation: Conduct the main ReaxFF MD simulation with appropriate parameters. Typical ReaxFF simulations use a time step of 0.25 fs, though smaller time steps may be necessary for high-temperature studies (>1500 K) [8]. The covalent interaction range is typically set to 5 Å, sufficient for most elements to capture even the weakest covalent interactions [2].

• Trajectory Analysis: Analyze the resulting trajectory to identify reaction events, intermediates, products, and mechanisms using specialized tools such as the Reaction Analysis and Visualization tool [10].

Computational Implementation

ReaxFF is implemented in several computational packages, each offering different capabilities:

• LAMMPS: The open-source Large-scale Atomic/Molecular Massively Parallel Simulator includes ReaxFF implementation suitable for large-scale combustion simulations [2].
• PuReMD: The Purdue Reactive Molecular Dynamics code is optimized specifically for ReaxFF simulations [2] [12].
• AMS: The Amsterdam Modeling Suite provides a user-friendly interface for ReaxFF simulations [11].
• ReaxFF/AMBER: A hybrid implementation that enables ReaxFF simulations within the AMBER package, particularly useful for biomolecular systems or when combining reactive and non-reactive regions [8].

Figure 3: Standard workflow for ReaxFF molecular dynamics simulations in combustion research

Software and Computational Tools

Table: Essential Software Tools for ReaxFF Combustion Simulations

Tool Name	Function	Application in Combustion Research
LAMMPS	Molecular dynamics simulator	Large-scale parallel ReaxFF simulations of combustion systems
PuReMD	Reactive molecular dynamics	Optimized ReaxFF implementation for efficient reaction sampling
AMS/GUI	Graphical interface for ReaxFF	System setup, visualization, and preliminary simulations
OVITO	Visualization and analysis	Trajectory analysis and reaction mechanism identification
OMEAN	Reaction analysis	Automated identification of reaction events from trajectories

Specialized Force Fields for Combustion Systems

Combustion researchers can leverage several specifically parameterized force fields:

• CHO.ff: The foundational hydrocarbon oxidation force field parameterized against DFT calculations on dissociation energies for bonds containing carbon, oxygen, and hydrogen, along with angle and torsion potentials [11].
• ReaxFFCHO-S22: A recently re-parameterized force field for hydrocarbon combustion that accurately captures primary intermediates and reaction pathways observed in experimental studies [10].
• HCONSB.ff: Extends capabilities to nitrogen- and sulfur-containing species important in fossil fuel combustion and pollutant formation [11].

Application Case Study: Cyclohexane Combustion Mechanisms

Research Context and Methodology

A recent application of ReaxFF to cyclohexane combustion demonstrates the methodology's power for elucidating complex reaction networks [10]. Cyclohexane serves as a typical surrogate for hydrocarbon fuels in engine emissions research, and understanding its combustion is crucial for improving efficiency and reducing pollutants. The study employed the specialized ReaxFFCHO-S22 force field in MD simulations across various temperatures, densities, and equivalence ratios while also considering the influence of water [10].

The simulation protocol involved constructing a system with 12 cyclohexane and 108 O₂ molecules in a periodic box at a density of 0.2 g/cm³, achieving a target temperature of 2800 K—a condition commonly adopted in reactive MD simulations to accelerate the reaction timeline within computationally feasible simulation times [10]. The production simulations ran for 1 ns, during which researchers tracked the temporal evolution of molecular species and reaction events [10].

Key Findings and Mechanistic Insights

The ReaxFF simulations revealed that cyclohexane combustion predominantly begins with hexagon ring opening via homolytic C-C bond cleavage, followed by subsequent oxidation of smaller hydrocarbon radicals [10]. This primary mechanism aligns well with experimental observations, validating the ReaxFF approach. The simulations identified significant intermediates and products including C₂H₄, CH₂O, CO, CO₂, and H₂O, providing atomic-level insight into the formation pathways of both desirable complete combustion products and potentially harmful intermediates [10].

Notably, the study found that elevated temperatures and higher densities facilitate the oxidation process of cyclohexane, while the incorporation of small amounts of H₂O promoted the conversion of CO to CO₂—a crucial finding for understanding and optimizing clean combustion mechanisms [10]. The ReaxFF methodology enabled tracking of the complete reaction network without predefined pathways, demonstrating its value for discovering non-intuitive reaction mechanisms in complex combustion environments.

Future Directions and Emerging Alternatives

Advancements in ReaxFF Methodology

The ReaxFF method continues to evolve with several promising directions enhancing its applicability to combustion research. Recent developments include hybrid QM/ReaxFF approaches that combine the accuracy of quantum mechanics for the reactive center with the efficiency of ReaxFF for the environment [8]. Additionally, improved parameter optimization techniques using machine learning and global optimization algorithms are producing more accurate and transferable force fields [9]. Extensions to new elements and chemical environments continue to expand the range of combustion systems accessible to ReaxFF simulation.

Neural Network Potentials as Emerging Alternatives

While ReaxFF represents a significant advancement over classical force fields, recent years have seen the emergence of neural network potentials (NNPs) that offer potentially higher accuracy while maintaining computational efficiency. The EMFF-2025 model represents a general NNP for C, H, N, O-based high-energy materials that achieves DFT-level accuracy in predicting structures, mechanical properties, and decomposition characteristics [13].

Such machine learning potentials leverage transfer learning with minimal data from DFT calculations, offering a promising alternative that may address some limitations of ReaxFF in describing reaction potential energy surfaces [13]. However, NNPs currently face challenges in automated exploration of reactive chemical space during sampling, as they require simultaneous exploration of molecular species changes and structural variations associated with non-equilibrium thermodynamic processes [13]. For the foreseeable future, ReaxFF remains a robust and widely validated approach for combustion mechanism discovery, particularly for systems where extensive training data for NNPs is not yet available.

ReaxFF's bond-order formalism represents a fundamental breakthrough in molecular simulation that has dramatically expanded the scope of atomic-scale investigations in combustion chemistry. By enabling dynamic bond formation and breaking within a computationally efficient framework, ReaxFF bridges the critical gap between quantum mechanical accuracy and classical molecular dynamics scale. The methodology has proven invaluable for elucidating complex reaction networks in hydrocarbon oxidation, pollutant formation, and energetic material decomposition—systems where traditional force fields fall short and quantum methods remain computationally prohibitive.

While emerging neural network potentials show promise for future applications, ReaxFF remains a mature, validated, and widely accessible tool for combustion researchers. Its continued development through improved parameterizations and hybrid approaches ensures that it will remain relevant for addressing fundamental challenges in understanding and optimizing combustion processes across energy, propulsion, and environmental applications.

Atomistic-scale computational techniques provide a powerful means for exploring and optimizing properties of novel materials. In combustion chemistry research, where processes involve complex reactions across gas, liquid, and solid phases, accurately simulating bond formation and breaking is paramount. Methods based on quantum mechanics (QM) offer valuable theoretical guidance at the electronic level but are often too computationally intense for simulations considering the full dynamic evolution of a system. In contrast, empirical interatomic potentials based on classical principles require significantly fewer computational resources, enabling simulations of dynamic processes over longer timeframes and larger scales. However, such classical methods typically require predefined connectivity between atoms, precluding simulations involving reactive events. The Reactive Force Field (ReaxFF) method was developed to help bridge this gap, approaching from the classical side by casting the empirical interatomic potential within a bond-order formalism, thus implicitly describing chemical bonding without expensive QM calculations [2].

For combustion research, this classical treatment of reactive chemistry has opened the door for numerous studies of phenomena occurring on scales previously inaccessible to computational methods. ReaxFF enables simulations involving reactive events at interfaces between solid, liquid, and gas phases because its description of each element is transferable across phases. For example, an oxygen atom is treated with the same mathematical formalism whether in the gas phase as O₂, in the liquid phase within an H₂O molecule, or incorporated in a solid oxide [2]. This transferability, coupled with lower computational expense allowing for longer simulation timescales, enables ReaxFF to consider phenomena dependent not only on species reactivity but also on dynamic factors like diffusivity and solubility that affect how species migrate through complex combustion systems [2] [14].

Theoretical Foundation: Core Energy Contributions in ReaxFF

The ReaxFF potential employs a sophisticated combination of energy terms to describe both reactive and non-reactive interactions between atoms. This allows ReaxFF to accurately model both covalent and electrostatic interactions for diverse materials encountered in combustion environments, from hydrocarbon fuels to their oxidation products and soot precursors [2] [15]. The total system energy in ReaxFF is described by the equation:

[E{\text{system}} = E{\text{bond}} + E{\text{over}} + E{\text{angle}} + E{\text{tors}} + E{\text{vdWaals}} + E{\text{Coulomb}} + E{\text{Specific}}]

where each term addresses specific physical interactions, with the bond-order, overcoordination, and polarizable charge contributions being particularly crucial for modeling chemical reactivity [2] [16].

Table: Core Energy Contributions in ReaxFF

Energy Term	Symbol	Physical Significance	Role in Combustion Chemistry
Bond Energy	(E_{\text{bond}})	Energy from bond formation/breaking	Models fuel decomposition and product formation
Overcoordination Penalty	(E_{\text{over}})	Prevents atomic overcoordination	Maintains proper valence during radical reactions
Valence Angle Strain	(E_{\text{angle}})	Three-body angle strain energy	Controls molecular geometry evolution
Torsional Strain	(E_{\text{tors}})	Four-body torsional strain energy	Affects conformational changes in hydrocarbons
van der Waals	(E_{\text{vdWaals}})	Dispersive interactions	Important for molecular association and soot precursor formation
Coulomb	(E_{\text{Coulomb}})	Electrostatic interactions	Critical for charge transfer in oxidation reactions
Specific Corrections	(E_{\text{Specific}})	System-specific terms (lone-pair, conjugation, etc.)	Enhances accuracy for specific functional groups

Bond-Order Formalism: The Cornerstone of Reactivity

The bond-order formalism represents the foundational innovation that enables ReaxFF to simulate chemical reactions without predefined connectivity. In ReaxFF, bond order is empirically calculated from interatomic distances using a continuous function:

[ BO{ij} = BO{ij}^\sigma + BO{ij}^\pi + BO{ij}^{\pi\pi} = \exp\left[p{bo1}\left(\frac{r{ij}}{r0^\sigma}\right)^{p{bo2}}\right] + \exp\left[p{bo3}\left(\frac{r{ij}}{r0^\pi}\right)^{p{bo4}}\right] + \exp\left[p{bo5}\left(\frac{r{ij}}{r0^{\pi\pi}}\right)^{p{bo6}}\right] ]

where (BO) is the bond order between atoms i and j, (r{ij}) is interatomic distance, (r0) terms are equilibrium bond lengths, and (p_{bo}) terms are empirical parameters [2]. This equation is continuous, containing no discontinuities through transitions between σ, π, and ππ bond character, which yields a differentiable potential energy surface required for calculating interatomic forces [2].

This bond-order formula accommodates long-distance covalent interactions characteristic in transition state structures, allowing the force field to accurately predict reaction barriers. This covalent range is typically taken to be 5 Ångstrom—sufficient for most elements to capture even the weakest covalent interactions—but can be extended for elements with very large covalent radii [2]. This innovative approach replaces the fixed bond topologies of classical force fields with dynamic bond orders that continuously evolve during simulations, enabling the simulation of bond formation and breaking that is fundamental to combustion chemistry [9] [16].

Overcoordination and Undercoordination Treatments

In the context of ReaxFF, overcoordination and undercoordination refer to situations where atoms are either bonded with more or fewer neighbors than expected based on bond formation rules [16]. Overcoordination occurs when an atom forms more bonds than typical for its valence, which can happen when bond orders are artificially inflated, resulting in unrealistic geometries. Undercoordination happens when an atom has fewer bonds or neighbors than expected based on its chemical environment, which can make structures energetically unfavorable [16].

ReaxFF addresses these challenges through specific energy terms. The overcoordination penalty ((E_{over})) applies a stiff energy penalty if an atom forms more bonds than allowed by atomic valence rules (e.g., a carbon atom forming more than four bonds) [2]. This is complemented by undercoordination terms that correct for insufficient bonding. These corrections are essential for maintaining chemical realism during reactive events in combustion, such as when hydrocarbon radicals undergo decomposition or recombination reactions where traditional valence rules might be temporarily violated [2] [16].

The parameters governing these coordination-dependent terms are carefully optimized during force field development to prevent unrealistic atomic configurations while allowing legitimate reactive processes to occur. This ensures that intermediate species formed during combustion reactions, such as the transition states in fuel decomposition pathways, are properly stabilized without allowing unphysical structures to persist [2].

Polarizable Charge Descriptions

ReaxFF employs charge equilibration schemes to represent the dynamics of electron density during simulations [16]. Unlike fixed partial charges in classical force fields, ReaxFF supports several charge equilibration methods:

• Charge Equilibration (QEq): Adjusts atomic charges to equalize electronegativity across the system
• Atom-Condensed Kohn-Sham DFT to second order (ACKS2): Provides a more sophisticated description of charge distribution
• Charge Transfer and Polarization in Equilibrium (QTPIE): Handles charge transfer and polarization effects [16]

These polarizable charge descriptions are critical for combustion simulations where charge distribution evolves continuously as bonds form and break. The Coulombic energy ((E_{Coulomb})) is calculated between all atoms, regardless of connectivity and bond-order, using these dynamic charges [2]. This approach ensures that electrostatic interactions, which drive many chemical processes in combustion including proton transfers, ion-molecule reactions, and oxidation processes, are properly represented throughout the simulation.

Comparative Analysis: ReaxFF vs. Classical Force Fields in Combustion Research

The fundamental differences between ReaxFF and classical force fields have significant implications for their application in combustion chemistry research. Classical force fields rely on fixed bonding topologies, making them unsuitable for simulating chemical reactions where bonds form and break. In contrast, ReaxFF's bond-order formalism and dynamic charge equilibration enable the simulation of complex reactive processes central to combustion [2] [9].

Table: ReaxFF vs. Classical Force Fields for Combustion Applications

Feature	Classical Force Fields	ReaxFF
Bond Treatment	Fixed connectivity, predefined bonds	Dynamic bond orders based on interatomic distances
Chemical Reactions	Cannot simulate bond formation/breaking	Explicitly models reaction chemistry
Charge Treatment	Fixed partial charges	Polarizable charge equilibration (QEq, ACKS2, QTPIE)
Phase Transferability	Limited, often parameterized for specific phases	High, consistent description across solid, liquid, gas phases
Combustion Applications	Limited to physical processes (mixing, adsorption)	Full reaction mechanisms, soot formation, catalytic combustion
Computational Cost	Lower	Higher than classical but significantly lower than QM methods

For combustion research, this capability translates to the ability to simulate complete reaction networks from initial fuel decomposition through intermediate oxidation to final products. Recent studies have demonstrated this power in applications including:

• Fuel Pyrolysis and Oxidation: Revealing initial reaction mechanisms of toluene, ethanol, and n-heptane mixed fuels [17]
• Soot Formation: Providing atomistic insights into polycyclic aromatic hydrocarbon (PAH) formation in novel biofuels and JP-10 jet fuel [15]
• Nanoparticle Synthesis: Simulating flame synthesis of nanomaterials including boron nitride nanostructures [14] [15]
• Supercritical Combustion: Investigating n-dodecane and toluene combustion under supercritical conditions where traditional kinetic mechanisms face limitations [15]

The transferability across phases inherent in ReaxFF is particularly valuable for combustion systems, which often involve interfaces between gaseous fuels, liquid droplets, and solid surfaces or soot particles [2] [14]. This enables researchers to study complex multiphase combustion phenomena that were previously inaccessible to atomistic simulation.

Practical Implementation: Parameters and Force Field Development

The ReaxFF method requires careful parameterization to ensure accurate description of chemical reactivity. The force field utilizes numerous parameters organized in specific sections within force field files: General parameters (41), Atoms (32 per atom type), Bonds (16), Off-diagonal (6), Angles (7), Torsions (7), and Hydrogen bonds (4) [18].

Key Parameters Governing Core Energy Contributions

Table: Selected ReaxFF Parameters for Core Energy Contributions

Parameter Category	Key Parameters	Function	Equation Reference
Bond Order	(r0^\sigma), (r0^\pi), (r0^{\pi\pi}), (p{bo1})-(p_{bo6})	Define bond order calculation from interatomic distances	Eq. 2 [2]
Overcoordination	(p{boc1}), (p{boc2})	Control energy penalty for overcoordinated atoms	Eq. 4c, 4d [18]
Valence	(Vali), (Vali^e)	Define atomic valency and number of valence electrons	Eq. 3a, 4b, 5, 9a [18]
van der Waals	(r{vdW}), (D{ij}), (alpha_{ij})	Govern dispersive interactions	Eq. 23a, 23b [18]
Charge Equilibration	(gamma_i), (chiEEM), (etaEEM)	Control charge distribution and electrostatic interactions	Eq. 24 [18]

Two parameters of particular importance for practical simulations are the upper taper radius (parameter #13), which describes the non-bonded cutoff radius, and the bond order cutoff (parameter #30), which describes the bond order threshold above which atoms are considered connected. Both parameters significantly impact ReaxFF calculation speed; decreasing the taper radius or increasing the bond order cutoff can make ReaxFF run considerably faster. However, these parameters substantially influence force description and should not be changed without re-parameterization [18].

Force Field Development and Validation

Developing ReaxFF parameters requires extensive training sets covering relevant chemical phase space, including bond and angle stretches, activation and reaction energies, equation of state, surface energies, and more [9]. Training data is typically generated with electronic structure methods, often DFT calculations using accurate functionals [9]. Due to the complexity of the ReaxFF functional form, global optimization techniques offer the best chance to obtain parameter sets that closely describe training data [9].

FitSNAP-ReaxFF provides a workflow for optimizing ReaxFF parameters using the Covariance Matrix Adaptation Evolution Strategy (CMAES) algorithm, which minimizes the sum of squared errors between DFT reference data and predicted energies/forces given current parameter values [16]. This approach enables researchers to develop specialized force fields for specific combustion applications, such as the reparameterization for gas-phase B-N chemistry to model formation of boron nitride nanostructures [15].

Experimental Protocols and Methodologies

Typical ReaxFF-MD Simulation Workflow for Combustion Studies

Implementing ReaxFF molecular dynamics (ReaxFF-MD) simulations for combustion chemistry follows a systematic protocol:

• System Initialization:

  • Construct initial simulation cell with fuel and oxidizer molecules at appropriate densities
  • For mixed fuel studies like toluene/ethanol/n-heptane, use weight percentages relevant to applications (e.g., 10%/10%/80%) [17]

• Equilibration:

  • Apply energy minimization to remove bad contacts
  • Equilibrate system at target temperature using thermostats (e.g., Berendsen thermostat)
  • Use temperature ranges appropriate for combustion phenomena (typically 2000-4500K for high-temperature pyrolysis and oxidation) [17]

• Production Simulation:

  • Conduct reactive MD simulations with time steps typically around 0.1-0.5 fs
  • Run multiple trajectories (e.g., 20 trajectories per condition) to improve statistics [19]
  • Simulate for durations sufficient to observe phenomena of interest (e.g., several hundred picoseconds for fullerene formation) [19]

• Analysis:

  • Track molecular populations over time to determine reaction rates
  • Analyze chemical species evolution and intermediate formation
  • Characterize reaction networks and pathways
  • Calculate activation energies from temperature-dependent kinetics [17]

ReaxFF-MD Simulation Workflow

Specific Combustion Chemistry Application: Mixed Fuel Oxidation

A detailed ReaxFF-MD methodology for studying mixed fuel combustion illustrates the application of these core energy contributions [17]:

System Setup:

• Create simulation cell containing toluene, ethanol, and n-heptane molecules in appropriate proportions
• Add oxygen molecules to represent oxidizer environment
• Implement fuel-rich conditions with specific oxygen-carbon ratios

Simulation Parameters:

• Temperature range: 2000-4500 K with 500 K intervals
• Time step: 0.1 fs for numerical stability with reactive events
• Simulation duration: Sufficient to observe major decomposition and intermediate formation (typically 100-500 ps)
• Multiple trajectories (5-10) at each temperature for statistical reliability

Analysis Methods:

• Track reactant consumption rates to determine initiation mechanisms
• Identify key intermediates (C₂H₄, CH₃, OH radicals)
• Monitor product formation (CO, CO₂)
• Calculate activation energies from Arrhenius plots of temperature-dependent rate constants
• Analyze effects of ethanol/toluene addition on reaction pathways, particularly their influence on OH radical concentrations and C₂H₄ consumption [17]

Table: Computational Tools for ReaxFF Combustion Research

Tool Category	Specific Software/Resources	Combustion Application
MD Engines	LAMMPS (KOKKOS package), PuReMD, ADF, Materials Studio	High-performance ReaxFF-MD simulations of combustion processes
Force Field Databases	ReaxFF parameter sets for C/H/O/N, hydrocarbons, biofuel components	Provides validated parameters for specific fuel types
Analysis Tools	Custom scripts for species tracking, reaction analysis, visualization packages	Post-processing of simulation trajectories to extract chemical insights
Quantum Chemistry	DFT codes (VASP, Quantum ESPRESSO) for training data generation	Force field parameterization and validation for novel fuel molecules
Optimization Frameworks	FitSNAP-ReaxFF with CMAES algorithm	Development of new force fields for specialized combustion applications

Specialized Computational Approaches for Combustion

Combustion simulations often require specialized computational approaches to drive reactions and analyze results:

Hydrogen Removal Strategies: In fullerene formation studies, the hydrogen to carbon (H/C) ratio is gradually reduced during MD simulations to drive the combustion process, mimicking actual combustion conditions where hydrogen is preferentially oxidized [19].

Temperature Programming: Studies employ multiple temperature points (2000K, 2500K, 3000K, 3500K, 4000K, 4500K) to establish Arrhenius behavior and extract activation energies comparable to experimental values (e.g., 43.36-45.19 kcal mol⁻¹ for mixed fuel decomposition) [17].

Reactive Environment Modeling: Systems include relevant oxidizers (O₂) and potential catalysts to reproduce realistic combustion environments, enabling study of complex processes like soot formation and catalytic combustion [14] [15].

ReaxFF has established itself as a powerful computational tool for combustion research, providing atomic-level insights into complex reactive processes that bridge the gap between quantum chemistry and macroscopic simulations. The core energy contributions—bond-order formalism, overcoordination treatments, and polarizable charges—enable physically realistic simulation of reaction chemistry across multiple phases relevant to combustion systems.

Future developments in ReaxFF will likely focus on improving accuracy for specific reaction classes important in combustion, extending to broader elements and materials, and enhancing computational efficiency through machine learning approaches and advanced computing architectures [2] [14]. As computational power grows and force fields become more refined, ReaxFF is poised to play an increasingly important role in predicting combustion behavior, optimizing fuel formulations, and designing cleaner combustion systems through atomic-level insights unavailable from experimental methods alone.

The integration of ReaxFF with multiscale modeling frameworks represents a particularly promising direction, where atomic-scale insights from ReaxFF can inform mesoscale and continuum models used in engineering design of combustion devices [14]. This multiscale approach will further enhance the value of ReaxFF simulations in translating atomic-level understanding to practical improvements in combustion efficiency and emissions reduction.

In the computational study of combustion chemistry, the selection of a molecular dynamics (MD) force field is a foundational decision that directly determines the phenomena a researcher can observe. This choice often narrows to a dichotomy between classical (harmonic) force fields, which excel at simulating systems in or near equilibrium, and reactive force fields, which are engineered to model the bond-breaking and bond-forming events characteristic of non-equilibrium processes like combustion. Framed within the specific context of ReaxFF versus classical force fields, this guide provides a structured philosophy for method selection. We dissect the underlying theories, provide direct performance comparisons, and outline detailed protocols for applying these methods to combustion-related research, empowering scientists to align their computational tools with their research objectives effectively.

Theoretical Foundations and Force Field Philosophies

The core distinction between classical and reactive force fields lies in their treatment of chemical bonds and their resulting applicability to different thermodynamic states.

Classical Force Fields: The Equilibrium Specialist

Classical force fields, such as PCFF, OPLS-AA, AMBER, and CHARMM, are built on a philosophy of computational efficiency and accuracy for stable molecular systems [3] [20]. They use a fixed-bond topology, meaning the chemical connectivity between atoms is defined at the beginning of a simulation and does not change. The energy calculation is typically a sum of harmonic potentials for bonded interactions (bonds, angles, dihedrals) and pair-wise potentials for non-bonded interactions (van der Waals, Coulombics) [20] [21].

• Functional Form Philosophy: The energy for a bond between atoms i and j is typically described by a harmonic potential: ( E{bond} = k{ij}(r{ij} - r{0,ij})^2 ), where ( k{ij} ) is the force constant and ( r{0,ij} ) is the equilibrium bond length [3]. This parabola is an excellent approximation near the energy minimum (equilibrium) but becomes increasingly unphysical as the bond is stretched, never allowing it to break.
• Applicability: This formulation makes classical FFs ideal for calculating equilibrium properties such as density, volumetric behavior, diffusion coefficients, and structural dynamics in polymers, biomolecules, and liquids where no chemical reactions occur [22] [21]. Their computational cost is low, enabling simulations of large systems (>100,000 atoms) for long timescales (nanoseconds to microseconds) [20].

Reactive Force Fields (ReaxFF): The Non-Equilibrium Generalist

ReaxFF was developed to bridge the gap between quantum mechanical accuracy and the scale of classical MD for reactive systems [2]. Its philosophy centers on dynamic bonding. Instead of pre-defined connections, it uses a bond-order formalism derived from interatomic distances to determine bond strength and existence on-the-fly [2] [23].

• Functional Form Philosophy: The bond order ( BO{ij} ) is a continuous function of the interatomic distance ( r{ij} ): ( BO{ij} = \exp[p{bo1}(\frac{r{ij}}{r{0,\sigma}})^{p{bo2}}] + \exp[p{bo3}(\frac{r{ij}}{r{0,\pi}})^{p{bo4}}] + \exp[p{bo5}(\frac{r{ij}}{r{0,\pi\pi}})^{p_{bo6}}] ) This empirical formula captures σ, π, and ππ bond contributions [2]. The bond order then feeds into the calculation of bond energy, angle strain, and torsion energy, creating a potential energy surface that is sensitive to the chemical environment and allows for bond dissociation.
• Applicability: This makes ReaxFF powerful for non-equilibrium processes like combustion, pyrolysis, catalysis, and fracture, where reaction pathways and products are not known a priori [23] [14]. It can simulate complex reaction networks involving radical species, which are ubiquitous in combustion, as shown in its application to bituminous coal combustion [23].

Emerging Hybrid and Alternative Approaches

The landscape of force fields is evolving. The Reactive INTERFACE Force Field (IFF-R) has been recently proposed as a hybrid approach [3]. It replaces the harmonic bond potential in classical FFs with a Morse potential, ( E{bond} = D{ij}[\exp(-\alpha{ij}(r{ij}-r{0,ij})) - 1]^2 ), which includes a finite bond dissociation energy ( D{ij} ) [3]. This method aims to maintain the simplicity and interpretability of classical FFs while adding bond-breaking capabilities, reportedly offering a ~30x speedup over ReaxFF [3]. Furthermore, Machine Learning Force Fields (ML-FFs) are emerging as a powerful tool, using data from quantum mechanics to achieve near-DFT accuracy at a fraction of the cost, though they currently require significant data and computational resources for inference [20] [21].

Table 1: Philosophical and Formal Comparison of Force Field Types

Feature	Classical Force Fields	Reactive Force Fields (ReaxFF)	Machine Learning Force Fields
Bond Treatment	Fixed, harmonic potentials	Dynamic, bond-order formalism	Learned from QM data
Reactivity	Non-reactive	Fully reactive	Potentially reactive
Computational Cost	Low	High (~30x Classical)	Very High (~100x Classical)
Primary Domain	Equilibrium properties, structure, dynamics	Non-equilibrium processes, reaction mechanisms	High-accuracy property prediction
Key Strength	Speed, stability for large systems	Captures complex chemistry without pre-definition	Near-QM accuracy for complex PES
Key Weakness	Cannot model chemical reactions	Computationally expensive; complex parameterization	Data hunger; high computational cost

Quantitative Comparison and Decision Framework

Selecting the appropriate force field requires a balance of accuracy, computational cost, and the specific scientific question. The data in Table 2 provides a quantitative basis for this decision.

Table 2: Performance Benchmarking for Combustion-Related Systems

System / Property	Classical FF (e.g., PCFF, OPLS)	ReaxFF	Reference Method	Notes & Context
Polyethylene Energy/Forces	OPLS & PCFF show closest agreement with DFT	Shows deviation from DFT	DFT Calculations [22]	Under extreme conditions away from equilibrium
Polyethylene Melting Behavior	PCFF agrees with experiment	PCFF, TraPPE-UA, ReaxFF perform equally well vs. experiment	Experimental Phase Diagram [22]	For phase behavior near equilibrium
Coal Combustion Pathways	Not applicable	Tracks oxidation & cracking of ethers, carbonyls, etc.	DFT Validation [23]	Identifies key radicals (·OH, ·HO₂) and intermediates
Computational Speed (vs. Classical)	1x (Baseline)	~0.03x (30x slower)	N/A	General scaling estimate [3] [20]
H₂ Combustion PES Sampling	Not applicable	Suitable for MD	ωB97X-V/cc-pVTZ [24]	ReaxFF trained on such QM data for reactivity

The data leads to a clear decision framework, visualized below.

Diagram 1: Force Field Selection Workflow

Experimental Protocols for Combustion Chemistry

To ensure reproducible and meaningful results, following a rigorous protocol is essential. Below are detailed methodologies for setting up and validating simulations for combustion chemistry.

Protocol 1: ReaxFF MD for Fuel Combustion Pathway Analysis

This protocol is adapted from studies investigating the microscopic reaction behaviors of bituminous coal combustion [23].

• System Construction:

  • Obtain or construct a validated molecular model of the fuel (e.g., the Wiser bituminous coal model).
  • Place a single fuel molecule or a small cluster in a large simulation box (e.g., > 50 Å side length).
  • Add a sufficient number of O₂ molecules to achieve the desired oxygen-deficient or stoichiometric conditions. Fill the remainder of the box with N₂ or Ar to represent an inert atmosphere and control density.

• Force Field and Parameterization:

  • Select a ReaxFF parameter set from the "combustion branch" (e.g., CHO.ff or HCONSB.ff) that is trained for hydrocarbon oxidation [11] [2].
  • Energy minimization: Perform a steepest descent or conjugate gradient minimization to remove bad contacts and high-energy configurations in the initial structure.

• Equilibration and Production Run:

  • Equilibration: Run an NVT (canonical ensemble) simulation for 10-50 ps at the target high temperature (e.g., 2000-3000 K) using a thermostat like Berendsen or Nosé-Hoover. This allows the system to reach thermal equilibrium.
  • Production: Switch to an NVE (microcanonical ensemble) simulation for 100-500 ps to model the reactive dynamics without artificial temperature control influencing the reaction energetics. Use a time step of 0.1-0.25 fs.

• Trajectory Analysis:

  • Use bond-order analysis (built into tools like LAMMPS or AMS) to track the formation and dissociation of specific functional groups (e.g., ether bonds, carbonyl groups, hydroxyls).
  • Identify and count the population of key radical species (·OH, ·HO₂) and stable products (CO, CO₂, H₂O) over time.
  • Perform a reaction pathway analysis by tracing the chemical evolution of specific atoms or molecular fragments through the simulation.

Protocol 2: Classical MD for Thermodynamic and Transport Properties

This protocol is suitable for calculating properties like density, volumetric changes, and diffusion in fuels or polymers under pre-combustion conditions [22].

• System Construction:

  • Build an amorphous cell containing multiple polymer chains (e.g., polyethylene) or fuel molecules, targeting an experimentally known density.
  • Use the PACKMOL or Amorphous Cell module in materials studio to generate a statistically representative initial configuration.

• Force Field Selection:

  • Choose a Class II force field like PCFF or COMPASS, which includes cross-coupling terms and is known for good accuracy in predicting condensed-phase properties of polymers and organics [22] [20].

• Equilibration and Production:

  • Energy Minimization: Thoroughly minimize the energy of the initial structure.
  • Density Equilibration: Run an NPT (isothermal-isobaric ensemble) simulation for 1-5 ns at the target temperature and pressure (e.g., 300 K, 1 atm) using a barostat (e.g., Berendsen, Parrinello-Rahman) to relax the density.
  • Production Run: Perform a longer NPT or NVT simulation (5-50 ns) with a 1 fs time step. Save the trajectory every 1-10 ps for analysis.

• Property Calculation:

  • Density: Calculate as the average mass/volume over the production run.
  • Glass Transition Temperature (Tg): Run a series of NPT simulations at decreasing temperatures and plot specific volume vs. temperature; the intersection of the linear fits identifies Tg.
  • Diffusion Coefficient: Calculate the Mean Squared Displacement (MSD) of molecules from the trajectory and use the Einstein relation: ( D = \frac{1}{6N{a}}\lim{t \to \infty} \frac{d}{dt} \sum{i=1}^{N{a}} \langle | \vec{r}i(t) - \vec{r}i(0) |^2 \rangle ), where ( N{a} ) is the number of atoms, and ( \vec{r}i ) is their position [21].

Validation and Benchmarking Strategy

Regardless of the method chosen, validation against experimental or high-level theoretical data is crucial.

• For ReaxFF: Validate by comparing the bond dissociation energies (BDEs) and valence angle bend energies of small, relevant molecules (e.g., H₂O, CH₄, O₂) against DFT calculations [23]. As shown in Table 2, this ensures the force field describes the fundamental chemistry correctly before applying it to a complex system.
• For Classical FFs: Validate by comparing predicted mass densities, enthalpies of vaporization, and elastic moduli against available experimental data [22] [3]. For polyethylene, PCFF has been shown to exhibit close agreement with experimental melting behavior [22].

This section details the essential "reagents" and computational tools required for research in this field.

Table 3: Essential Research Reagents and Resources

Item Name	Function / Description	Example Use Case
ReaxFF Combustion Branch	A set of parameters (e.g., CHO.ff) optimized for gas-phase hydrocarbon oxidation reactions.	Simulating the combustion pathways of hydrocarbon fuels [11] [2].
ReaxFF Aqueous Branch	A set of parameters with a focus on aqueous-phase chemistry and solvation effects.	Studying corrosion or electrochemical processes in combustion-related environments [11].
PCFF Force Field	A Class II classical force field with cross-coupling terms for accurate condensed-phase properties.	Predicting density, mechanical properties, and phase behavior of polymers and organic fuels [22] [20].
Wiser Bituminous Coal Model	A widely accepted 2D molecular structure model for bituminous coal.	Serving as a representative, complex fuel structure for ReaxFF MD combustion studies [23].
LAMMPS MD Package	A highly versatile and open-source MD simulation package that supports both classical and ReaxFF.	Performing large-scale, high-performance production MD simulations [2] [14].
AMS/ReaxFF Package	A commercial software (SCM) with a dedicated implementation of ReaxFF and a graphical interface.	Setting up, running, and analyzing ReaxFF simulations with a user-friendly GUI [11].
QM Benchmark Dataset	High-quality QM data (e.g., H₂ combustion PES with ~290,000 energies) for training/validation.	Validating ReaxFF results or training new ML force fields [24].

Implementing ReaxFF in Combustion Research: From Aerospace Fuels to Energetic Materials

Combustion research requires a detailed understanding of complex chemical reactions, including pyrolysis and oxidation pathways across diverse hydrocarbon fuels. While comprehensive kinetic models exist for smaller hydrocarbons (typically C3 or lower), developing and validating reaction mechanisms for larger hydrocarbons presents a daunting challenge due to the complexity of their reaction networks [25]. Molecular dynamics (MD) simulations have emerged as powerful computational tools for fundamental research across science and engineering disciplines, with particular relevance for combustion and energy systems [14]. Classical MD methods, while computationally efficient, typically require predefined atomic connectivity and cannot simulate chemical reactions where bonds break and form [2]. This limitation has driven the development of reactive force fields (ReaxFF), which bridge the gap between quantum mechanical methods that are accurate but computationally expensive, and classical force fields that are efficient but non-reactive [2].

The ReaxFF method, first introduced in 2001, employs a bond-order formalism that enables the simulation of chemical reactions within a molecular dynamics framework [2]. By calculating bond orders empirically from interatomic distances, ReaxFF implicitly describes chemical bonding without expensive quantum mechanical calculations, making it possible to simulate reactive events for larger systems and longer timescales than would be feasible with quantum methods [2]. This article examines the development and application of specific ReaxFF parameter sets for combustion chemistry, focusing on the CHO-2008 and CHO-2016 parameterizations, their comparative strengths and limitations, and their role in advancing combustion research.

ReaxFF Fundamentals: Bridging Quantum and Classical Approaches

Theoretical Framework of Reactive Force Fields

The ReaxFF method utilizes an empirical interatomic potential based on a bond-order formalism, allowing it to describe reactive events without the computational expense of quantum mechanical calculations [2]. Unlike classical force fields that maintain fixed atomic connectivity, ReaxFF employs a continuous function to calculate bond orders from interatomic distances, enabling seamless transitions between bonded and non-bonded states [2]. The total system energy in ReaxFF is described by multiple contributions:

Esystem = Ebond + Eover + Eangle + Etors + EvdWaals + ECoulomb + ESpecific [2]

Where Ebond represents bond energy, Eover is an over-coordination penalty, Eangle and Etors account for angle and torsion strain, EvdWaals and ECoulomb describe non-bonded interactions, and ESpecific includes system-specific terms [2]. This comprehensive energy description allows ReaxFF to accurately model both covalent and electrostatic interactions across diverse materials and chemical environments.

A key advantage of ReaxFF is its transferability across phases - the same mathematical formalism describes an oxygen atom whether it appears in gaseous O2, liquid H2O, or solid oxides [2]. This transferability, combined with computational efficiency that enables longer simulation timescales, allows ReaxFF to investigate phenomena dependent on both reactivity and dynamic factors such as diffusivity and solubility [2].

Comparison with Classical and Alternative Reactive Methods

Classical force fields for molecular materials typically employ harmonic bond potentials with fixed connectivity, making them unsuitable for simulating chemical reactions [6]. While these methods provide computational efficiency for studying systems at equilibrium, they cannot model bond dissociation or formation events essential for combustion chemistry [6].

Recent alternatives to ReaxFF have emerged, such as the Reactive INTERFACE Force Field (IFF-R), which replaces harmonic bond potentials with reactive Morse potentials [3]. This approach maintains compatibility with existing force fields like CHARMM and AMBER while enabling bond dissociation through interpretable Morse parameters [3]. Although IFF-R offers approximately 30-fold faster computation compared to ReaxFF [3], ReaxFF remains widely validated for complex combustion systems and provides a more comprehensive description of chemical environments through its bond-order formalism.

Hybrid approaches have also been developed, such as ReaxFF/AMBER, which integrates ReaxFF's reactive capabilities within the established AMBER molecular dynamics package [8]. This framework enables study of local reactive events in large systems at a fraction of the computational cost of QM/MM models [8].

The CHO Parameter Sets: Evolution from CHO-2008 to CHO-2016

The CHO-2008 Combustion Force Field

The initial ReaxFF combustion force field, developed by Chenoweth et al. (commonly referred to as CHO-2008), represented a significant advancement for atomistic-level investigation of combustion chemistry [25]. This parameter set achieved significant popularity for studying large hydrocarbon combustion and opened new avenues for researchers to investigate combustion chemistry from the atomistic level [26]. The CHO-2008 parameterization demonstrated particular utility for studying pyrolysis and oxidation of complex hydrocarbon fuels, providing reaction kinetics for complex fuel and fuel mixtures with accuracy approaching ab-initio-based methods but at significantly lower computational expense [25].

Despite its successes, the CHO-2008 description exhibited two significant limitations. First, it failed to accurately describe the chemistry of small hydrocarbon oxidation, especially the conversion of CO to CO2, which is highly relevant to syngas combustion [25]. Second, the parameterization obtained faster than expected hydrogen abstraction by O2 from hydrocarbons, thus underestimating the oxidation initiation temperature [25]. These limitations restricted its applicability across the full spectrum of combustion-relevant chemistries.

The CHO-2016 Advancements and Improvements

To address the limitations of CHO-2008, Ashraf and van Duin systematically improved the description, resulting in the CHO-2016 parameter set [25]. The retraining process expanded the original DFT-based training set by including reactions and transition state structures specifically relevant to syngas combustion and oxidation initiation pathways [25]. The aim was to retain the accuracy of the 2008 description for larger hydrocarbons while improving performance for small molecule chemistry and oxidation initiation.

Validation studies demonstrated that the CHO-2016 parameter set significantly improved the description of C1 chemistry and resolved the low-temperature oxidation initiation problem [25]. When applied to syngas and methane oxidation simulations, the new parameters correctly described CO to CO2 conversion kinetics [25]. Additionally, Arrhenius parameters for JP-10 decomposition and initiation mechanism pathways of n-butylbenzene pyrolysis maintained good agreement with both experimental data and CHO-2008 simulation results, confirming the transferability of the CHO-2016 description across a wide range of hydrocarbon chemistry [25].

Table 1: Comparison of CHO-2008 and CHO-2016 ReaxFF Parameter Sets

Feature	CHO-2008	CHO-2016
Development Basis	Original parameterization by Chenoweth et al.	Retrained with expanded training set including syngas pathways
Small Molecule Oxidation	Poor description of CO to CO2 conversion	Significantly improved C1 chemistry
Oxidation Initiation	Underestimated initiation temperature	Accurate initiation temperature prediction
Large Hydrocarbon Performance	Accurate for complex fuels	Maintained accuracy for fuels like JP-10 and n-butylbenzene
Syngas Combustion	Limited applicability	Greatly improved performance
Training Set	Original DFT-based training	Expanded with syngas and oxidation initiation transition states

Practical Implementation: Methodologies and Protocols

Simulation Setup and Workflow

Implementing ReaxFF combustion simulations requires careful attention to parameter selection, system setup, and simulation protocols. The fundamental workflow begins with selecting the appropriate parameter set based on the specific chemical system under investigation, with CHO-2016 generally preferred for systems involving small molecule oxidation or syngas chemistry.

Table 2: Research Reagent Solutions for ReaxFF Combustion Simulations

Component	Function	Implementation Considerations
ReaxFF Code	Provides core reactive MD capability	Implement via LAMMPS, PuReMD, or AMBER integration
Parameter Set	Defines force field parameters	CHO-2016 for general combustion; CHO-2008 for specific large hydrocarbon systems
System Builder	Creates initial molecular configurations	In-house scripts or molecular builder tools
Thermostat	Controls simulation temperature	NVT ensemble commonly used for combustion simulations
Analysis Tools	Processes trajectory and reaction data	Custom scripts for species tracking; visualization software

The simulation workflow involves multiple stages, from initial system construction through to reaction analysis, as illustrated in the following diagram:

Diagram 1: ReaxFF Combustion Simulation Workflow

Validation Protocols and Performance Assessment

Rigorous validation is essential for ensuring the accuracy and reliability of combustion simulations. The CHO-2016 parameter set was validated through high-temperature NVT-MD simulations studying oxidation and pyrolysis of four different hydrocarbon fuels: syngas, methane, JP-10, and n-butylbenzene [25]. This multi-fuel approach ensured transferability across different hydrocarbon chemistries.

For syngas and methane oxidation simulations, key validation metrics included accurate prediction of CO to CO2 conversion kinetics and oxidation initiation temperatures [25]. For larger hydrocarbons like JP-10 and n-butylbenzene, validation focused on matching Arrhenius parameters for decomposition and initiation mechanism pathways against both experimental data and CHO-2008 simulation results [25]. This comprehensive validation strategy confirmed that CHO-2016 retained the strengths of CHO-2008 for larger hydrocarbons while significantly improving small molecule chemistry description.

Simulation technical parameters typically employ a time step of 0.25 femtoseconds for most simulations, with smaller time steps required for higher temperature studies (>1500K) [8]. Bond orders are evaluated at each time step to determine atomic connectivity within a predefined distance cutoff, typically 5 Ångstroms [8].

Comparative Analysis: ReaxFF vs. Classical Force Fields in Combustion Research

Capability Assessment Across Simulation Domains

The choice between reactive and classical force fields depends significantly on the specific research objectives and chemical processes under investigation. The following table summarizes key comparative aspects:

Table 3: Functional Comparison of Force Field Approaches for Combustion Chemistry

Simulation Aspect	Classical Force Fields	ReaxFF
Bond Breaking/Formation	Not available without predefined reactions	Fully dynamic and automated
Computational Cost	Lower computational requirements	Higher than classical but lower than QM
System Size Scaling	Suitable for larger systems	Practical for medium to large systems
Chemical Transferability	Limited to parameterized chemistries	High transferability across phases
Reaction Mechanism Discovery	Cannot discover new pathways	Capable of revealing unexpected mechanisms
Charge Distribution	Fixed or limited polarization	Dynamic charge equilibration at each step
Parameterization Complexity	Simplified for specific chemistries	Complex training against QM and experimental data

Application-Specific Considerations

Classical force fields remain valuable for studying non-reactive aspects of combustion systems, such as fuel transport, phase behavior, and thermal properties, where their computational efficiency provides significant advantages [6]. However, ReaxFF enables investigation of reactive processes fundamental to combustion, including ignition chemistry, reaction pathways, and pollutant formation [14].

The ReaxFF method has been successfully deployed to gain fundamental insights into pyrolysis and oxidation of gas/liquid/solid fuels, revealing detailed energy changes and chemical pathways [14]. Furthermore, complex physico-chemical dynamic processes in catalytic reactions, soot formation, and flame synthesis of nanoparticles can be examined from an atomistic perspective [14]. These capabilities make ReaxFF particularly valuable for investigating complex, coupled processes where chemistry interacts with transport phenomena.

Future Directions and Emerging Methodologies

The continued evolution of reactive force field methodologies includes several promising directions. Machine learning approaches are being integrated to enhance parameterization and improve accuracy [14]. Multiscale modeling frameworks that combine ReaxFF with both quantum methods and coarse-grained approaches are extending simulation capabilities across broader temporal and spatial scales [14].

Emerging hybrid approaches like ReaxFF/AMBER demonstrate potential for studying local reactive events in large, complex systems, such as enzymatic processes or material interfaces [8]. The implementation of ReaxFF within widely distributed MD packages like AMBER also facilitates enhanced sampling methods, including umbrella sampling for mapping reaction profiles [8].

Recent developments in alternative reactive methodologies, such as IFF-R with its Morse potential-based approach, offer promising routes to combining the simplicity of harmonic potentials with reactive capabilities [3]. While these methods currently lack the extensive validation history of ReaxFF for combustion systems, they represent active areas of development that may offer advantages in specific applications.

The development of specialized parameter sets like CHO-2008 and CHO-2016 has significantly advanced the capabilities of reactive force fields for combustion research. The CHO-2016 parameter set in particular addresses critical limitations in small molecule oxidation and initiation chemistry while maintaining the strengths of its predecessor for larger hydrocarbons. As computational resources continue to expand and methodologies refine, ReaxFF simulations will play an increasingly valuable role in elucidating complex combustion phenomena, complementing both experimental investigations and higher-level quantum calculations. The appropriate selection of parameter sets based on specific chemical systems, coupled with rigorous validation protocols, enables researchers to leverage these powerful tools for advancing fundamental understanding and practical applications in combustion science and energy technology.

Probing Pyrolysis and Combustion Mechanisms of Aviation and Aerospace Fuels

The development of high-performance aviation and aerospace fuels requires a deep, mechanistic understanding of their pyrolysis and combustion behavior at the molecular level. These processes involve complex reaction networks with numerous parallel and sequential steps, including bond breaking, radical formation, and oxidation pathways, which occur at unprecedented speeds under extreme temperature and pressure conditions. Experimental techniques alone struggle to capture the fleeting intermediates and radical species that govern these reactions. Consequently, computational modeling has become an indispensable tool for probing these mechanisms. The central methodological choice for atomistic simulations lies in the selection of an appropriate force field—an empirical model that describes the potential energy of a system as a function of atomic coordinates. This technical guide examines the capabilities, applications, and limitations of two predominant approaches: classical non-reactive force fields and the reactive force field (ReaxFF), framing this comparison within the broader thesis of their respective roles in advancing combustion chemistry research.

Classical force fields, rooted in pre-defined molecular connectivity, utilize harmonic approximations for bond stretching and angle bending, combined with van der Waals and Coulombic potentials for non-bonded interactions [27]. While computationally efficient, their inherent inability to simulate bond breaking and formation limits their direct application to combustion chemistry. In contrast, ReaxFF, first developed in 2001, introduces a bond-order formalism derived from interatomic distances, allowing for dynamic bond formation and dissociation during simulations without the prohibitive computational cost of quantum mechanical (QM) methods [2]. By bridging the gap between accuracy and scale, ReaxFF has emerged as a powerful computational tool for exploring, developing, and optimizing material properties in reactive systems, opening the door to atomistic simulations of complex combustion phenomena over meaningful timescales [14] [2].

Theoretical Foundations: Force Fields from Classical to Reactive

Classical Non-Reactive Force Fields

Classical molecular dynamics (MD) relies on force fields constructed with a fixed bonding topology. The total potential energy ((E_{system})) is typically described as a sum of several terms:

(E{system} = E{bond} + E{angle} + E{torsion} + E{vdW} + E{Coulomb})

Here, (E{bond}) and (E{angle}) represent harmonic potentials for bond and angle strain, respectively; (E{torsion}) accounts for dihedral angle energetics; while (E{vdW}) (van der Waals) and (E_{Coulomb}) (electrostatic) describe non-bonded interactions [27]. The primary limitation for combustion studies is the absence of a term describing chemical reactivity. The fixed connectivity prevents simulations of pyrolysis (thermal decomposition) or oxidation, rendering classical MD suitable primarily for studying non-reactive physical properties in combustion, such as fuel droplet evaporation, transport phenomena, and phase behavior [14].

The Reactive Force Field (ReaxFF) Methodology

ReaxFF was developed to transcend the limitations of classical force fields by enabling reactive chemistry within a MD framework. Its core innovation is the use of a bond-order formalism, which implicitly describes chemical bonding without expensive QM calculations [2]. In ReaxFF, the bond order ((BO{ij})) between atoms (i) and (j) is calculated directly from their interatomic distance ((r{ij})) using an empirical formula that considers (\sigma), (\pi), and (\pi\pi) bonds:

(BO{ij} = BO{ij}^{\sigma} + BO{ij}^{\pi} + BO{ij}^{\pi\pi} = \exp\left[p{bo1}\left(\frac{r{ij}}{r0^{\sigma}}\right)^{p{bo2}}\right] + \exp\left[p{bo3}\left(\frac{r{ij}}{r0^{\pi}}\right)^{p{bo4}}\right] + \exp\left[p{bo5}\left(\frac{r{ij}}{r0^{\pi\pi}}\right)^{p{bo6}}\right])

This bond order is continuous and differentiable, allowing for smooth transitions during bond formation and dissociation [2]. The system's potential energy in ReaxFF incorporates bond-order-dependent terms (e.g., bond energy, angle strain, torsion) and bond-order-independent terms (e.g., van der Waals, Coulombic interactions). A critical component is an over-coordination penalty ((E_{over})), which prevents atoms from exceeding their chemical valence, and a charge equilibration scheme (QEq) that dynamically calculates atomic charges at each step [2]. This framework allows ReaxFF to accurately model the making and breaking of chemical bonds, fundamental to simulating pyrolysis and combustion.

Table 1: Fundamental Comparison Between Classical and Reactive Force Fields.

Feature	Classical Force Fields	Reactive Force Field (ReaxFF)
Bond Treatment	Pre-defined, fixed connectivity	Dynamic, based on continuous bond-order
Reactivity	Cannot simulate chemical reactions	Explicitly simulates bond breaking/formation
Energy Calculation	Harmonic potentials for bonds/angles	Bond-order dependent energy terms
Atomic Charges	Often fixed or based on simple models	Dynamically updated via charge equilibration
Computational Cost	Relatively low	Higher than classical, but far lower than QM
Primary Combustion Application	Physical properties, phase changes	Pyrolysis pathways, oxidation kinetics, ignition

Research Applications in Aviation and Aerospace Fuels

ReaxFF MD has been extensively applied to unravel the complex reaction mechanisms governing the behavior of advanced fuels. Its ability to provide atomistic insights into reaction pathways, ignition delay, and catalytic enhancement has made it a cornerstone of modern combustion research.

Pyrolysis and Oxidation of High-Density Fuels

High-density aerospace fuels, such as exo-tetrahydrodicyclopentadiene (exo-THDCPD, the primary component of JP-10), are prized for their high volumetric energy density, which is crucial for hypersonic propulsion. However, their stable, caged molecular structures lead to long ignition delays. ReaxFF MD simulations have been pivotal in understanding their decomposition. For instance, one study simulated the combustion of exo-THDCPD dispersed with Pt-graphene nanocatalysts [28]. The simulations, conducted at temperatures between 2600 K and 3000 K, tracked fuel and oxygen consumption rates, radical formation, and intermediate species. The results demonstrated that the Pt-graphene catalyst significantly reduced ignition delay by accelerating radical formation (e.g., H, OH) and enhanced combustion efficiency by promoting the oxidation of CO to CO₂, thereby suppressing the accumulation of intermediate hydrocarbons [28].

Another study on the high-strain cage hydrocarbon norbornane employed a complementary approach, using a neural network-assisted molecular dynamics (Deep Potential) to achieve high-precision simulation of its pyrolysis [29]. The analysis revealed that at lower temperatures (1273 K), hydrogen abstraction reactions dominated fuel consumption (69%), whereas at higher temperatures (1373 K), the mechanism shifted overwhelmingly (79%) towards C-C bond scission and ring-opening reactions [29]. These findings highlight how ReaxFF and next-generation MLIPs can identify temperature-dependent mechanistic shifts that are critical for optimizing fuel formulation for specific operating conditions.

Modeling Sustainable Aviation Fuels (SAFs)

The transition to sustainable aviation fuels (SAFs), particularly hydrotreated esters and fatty acids (HEFA) from waste cooking oil, is a key pathway for decarbonizing aviation. Understanding their combustion chemistry is essential for their integration into existing engines. A study on two CAAC-certified SAFs utilized a Hybrid Chemistry (HyChem) approach to develop their reaction mechanisms [30]. The HyChem model is a physics-based framework that separates the global combustion process into lumped reaction steps for initial fuel pyrolysis followed by detailed oxidation of the resulting fragments (e.g., C₂H₄, CH₄, C₆H₆) using a foundational mechanism like USC Mech II [31]. While not a direct ReaxFF application, the HyChem philosophy aligns with the multi-scale spirit of ReaxFF by simplifying the complex parent fuel chemistry. The study used a multi-objective genetic algorithm (NSGA-II) to optimize the lumped pyrolysis steps and found that reactions with H and OH radicals played the most important roles in system reactivity, a conclusion that can be validated and explored at the atomistic level using ReaxFF [30].

Energetic Nanocatalysts and Additives

The dispersion of nanocatalysts in fuels is a promising strategy for tuning combustion performance. ReaxFF MD is uniquely suited to probe the atomistic interactions between fuel molecules and catalyst surfaces. The study on Pt-graphene in JP-10 is a prime example [28]. Reaction pathway analysis revealed that the nanocatalyst shifted the fuel breakdown mechanism towards oxidation-driven pathways from the very initial stages, rather than pure thermal pyrolysis. This catalytic effect provides a molecular-level explanation for the observed reduction in ignition delay and more complete combustion, offering insights for designing next-generation catalytic additives for high-speed propulsion [28]. A recent review on ReaxFF MD applications in aviation/aerospace fuels and energetic additives further underscores its utility in screening and understanding the effects of various additives on fuel decomposition kinetics and reaction networks [32].

Experimental Protocols and Computational Methodologies

To ensure robustness and reliability, research in this field combines advanced simulations with rigorous experimental validation. Below are detailed methodologies for key techniques.

ReaxFF Molecular Dynamics Simulation Protocol

A typical ReaxFF MD simulation for probing combustion mechanisms follows a structured workflow, from system preparation to analysis [14] [28].

• System Initialization: An initial simulation cell is constructed containing fuel and oxidizer molecules (e.g., O₂) at a density corresponding to the target pressure. For catalytic studies, a catalyst nanoparticle or surface (e.g., Pt-graphene) is included.
• Force Field Parameterization: A validated set of ReaxFF parameters is selected. These parameters are trained against a large set of QM calculations and/or experimental data to ensure accuracy for the specific elements and reactions under study (e.g., C/H/O for hydrocarbon combustion) [2] [33].
• Equilibration: The system is first equilibrated in the NVT (constant Number, Volume, Temperature) or NPT (constant Number, Pressure, Temperature) ensemble at a lower temperature to achieve a stable starting configuration.
• Production Run: The equilibrated system is then subjected to high-temperature NVT simulations (e.g., 2000-3000 K) to mimic combustion conditions. The simulation is run for a sufficient time (nanoseconds) to observe decomposition and ignition.
• Trajectory Analysis: The atomic trajectory is analyzed to extract:
  • Fuel Consumption: The number of intact fuel molecules over time.
  • Species Progression: The formation and consumption profiles of key radicals (H, O, OH), intermediates (CO, C₂H₂), and final products (CO₂, H₂O).
  • Reaction Pathways: By tracking the evolution of specific atomic bonds, dominant reaction pathways can be identified and visualized.

Table 2: Key Research Reagents and Computational Tools in Combustion Chemistry.

Item / Software	Function / Description	Example in Use
ReaxFF Force Field	Reactive potential for simulating bond breaking/formation.	Modeling JP-10 pyrolysis pathways [28].
Classical Force Fields	Non-reactive potential for structure and property analysis.	Studying fuel droplet evaporation [14].
LAMMPS/PuReMD	Molecular dynamics software packages.	Running large-scale ReaxFF MD simulations [2].
HyChem Model	Hybrid (lumped + detailed) chemistry model for real fuels.	Developing reaction mechanisms for SAFs [30] [31].
Pt-Graphene Nanocatalyst	Catalytic additive to enhance ignition and combustion.	Accelerating THDCPD oxidation in ReaxFF MD [28].
Flow Reactor	Experimental apparatus for obtaining species time-history.	Validating pyrolysis mechanisms of SAFs [30].

Integration with Experiment and Machine Learning

The predictive power of any simulation is validated through comparison with experimental data. Key validation experiments include:

• Flow Reactor Studies: These provide speciation data (concentrations of reactants, intermediates, and products as a function of reaction time) during pyrolysis and oxidation at atmospheric pressure and temperatures around 1000 K, which serve as a direct benchmark for ReaxFF mechanisms [30].
• Ignition Delay Time (IDT) Measurements: Data from shock tube experiments on the time to auto-ignition under high-temperature, high-pressure conditions are a critical global parameter for validating the kinetic model's predictive capability [30].

Furthermore, the field is rapidly advancing through integration with machine learning (ML). ML techniques are now being used to optimize ReaxFF parameters more efficiently, improving the accuracy of energy and force predictions, particularly near transition states [33]. Another innovative approach involves pre-training Machine Learned Interaction Potentials (MLIPs) with massive amounts of cheap, classically generated non-reactive data before fine-tuning with a small set of high-quality QM data. This "teacher-student" strategy enhances the robustness and stability of MD simulations when exploring new regions of the potential energy surface, a common challenge in reactive system modeling [34].

Comparative Analysis: ReaxFF vs. Classical Force Fields in a Research Thesis

The choice between ReaxFF and classical force fields is not a matter of superiority but of application context. Each occupies a distinct niche in the multiscale modeling ecosystem for combustion research, as summarized in the diagram below.

Diagram 1: A multiscale modeling framework for combustion, showing the complementary roles of different computational methods.

Synergistic Roles in a Multiscale Framework

As illustrated in Diagram 1, ReaxFF and classical MD serve as connecting bridges between different scales of modeling:

• Classical MD is ideal for calculating non-reactive physical properties critical to engine performance, such as transport coefficients (diffusivity, viscosity), thermal conductivity, and liquid-vapor phase behavior of fuel droplets [14]. These properties are essential inputs for larger-scale Computational Fluid Dynamics (CFD) simulations of entire combustors. Its high computational efficiency allows it to probe these phenomena over longer timescales and larger system sizes than ReaxFF.

• ReaxFF MD operates at the crucial mesoscale, directly simulating the complex reaction networks that constitute pyrolysis and oxidation. It provides a fundamental understanding that is often inaccessible by experiment alone. The mechanistic insights and kinetic data generated by ReaxFF, such as the discovery of new reaction pathways or the effect of a nanocatalyst, can be used to inform and refine higher-level models like HyChem and detailed kinetic models [28] [32]. These models, in turn, provide the chemical reaction mechanisms for efficient CFD simulations.

Limitations and Future Directions

Despite its power, ReaxFF has limitations. Its accuracy is inherently tied to the quality and breadth of the training data used to parameterize it [2] [33]. The computational cost, while far lower than QM, is still significant, typically limiting simulations to nanosecond timescales, which requires elevated simulation temperatures to observe ignition events [28]. The future of force fields in combustion research lies in overcoming these challenges through:

• Machine-Learning Accelerated Force Fields: MLIPs are showing promise in achieving QM-level accuracy at a computational cost approaching that of classical MD [34]. Their development and training, potentially using ReaxFF data as a baseline, will further expand the accessible time and length scales for accurate reactive simulations.
• Automated Parameter Optimization: The use of ML-driven and genetic algorithm-based approaches (e.g., NSGA-II, JAX-ReaxFF) will streamline and improve the force field parameterization process, making it more robust and transferable [30] [33].
• Tighter Experiment-Simulation Integration: The synergy between flow reactor/shock tube experiments and ReaxFF simulations will continue to be vital for model validation and development, ensuring that theoretical predictions accurately reflect real-world chemistry [30].

The quest to probe the pyrolysis and combustion mechanisms of advanced aviation and aerospace fuels relies heavily on sophisticated computational tools. Within this landscape, ReaxFF has established itself as a transformative methodology, uniquely capable of providing atomistic-level insight into reactive processes that are central to fuel performance, ignition, and emissions. While classical force fields remain indispensable for simulating non-reactive physical phenomena, it is the reactive capability of ReaxFF that has enabled a fundamental shift from phenomenological observation to mechanistic prediction in combustion chemistry. The ongoing integration of machine learning and high-performance computing promises to further enhance the accuracy and scope of these simulations, solidifying the role of reactive molecular dynamics as a cornerstone for the rational design of next-generation fuels and propulsion systems.

Studying the Decomposition and Reactivity of Energetic Molecular Crystals

The study of energetic molecular crystals (EMCs), such as RDX, HMX, and TATB, is crucial for advancing applications in propulsion, explosives, and pyrotechnics. Understanding their decomposition and reactivity at the atomic level is fundamental to designing safer, high-performance materials. Computational modeling plays a pivotal role in elucidating complex reaction mechanisms that are difficult to probe experimentally. This technical guide focuses on the core methodologies of classical and reactive force fields (ReaxFF), framing their application within combustion chemistry research. It provides a comparative analysis, detailed protocols, and practical resources for researchers investigating the intricate chemistry of EMCs.

Force Field Fundamentals: A Comparative Framework

Molecular dynamics (MD) simulations compute atomic trajectories by solving Newton's equations of motion and are a cornerstone for investigating material properties at the atomistic level. The accuracy of these simulations is fundamentally governed by the force field (FF) employed, which describes the potential energy surface of the system. For the study of EMCs, FFs can be broadly categorized into two paradigms: classical non-reactive FFs and reactive FFs [14] [6].

Classical Force Fields utilize harmonic potentials for bonded interactions (bond stretching, angle bending) in conjunction with non-bonded potentials (e.g., van der Waals, Coulombic interactions). They are computationally efficient and excellent for predicting structural, mechanical, and thermal properties of EMCs, such as lattice parameters, density, elastic constants, and thermal expansion [6]. However, their fixed bonding connectivity prevents the simulation of chemical reactions, such as bond breaking and formation, which are central to decomposition and combustion processes [14] [2].

The Reactive Force Field (ReaxFF) bridges this gap. Developed by van Duin and colleagues, ReaxFF employs a bond-order formalism, where the bond order is empirically calculated from interatomic distances and is continuously updated during the simulation [2]. This allows bonds to break and form dynamically. ReaxFF incorporates polarizable charge descriptions and includes energy terms for over-coordination penalties, angle strain, and torsional strain, enabling it to simulate complex chemical reactions within large-scale MD frameworks with an accuracy close to quantum mechanics (QM) but at a fraction of the computational cost [14] [2].

Table 1: Core Comparison of Classical vs. Reactive Force Fields for Energetic Materials.

Feature	Classical Force Fields	Reactive Force Field (ReaxFF)
Bonding Representation	Fixed, harmonic bonds [6]	Dynamic, bond-order based [2]
Chemical Reactivity	Not possible [14]	Models bond breaking/formation [2]
Primary Applications	Polymorph prediction, crystal morphology, density, mechanical properties, thermal expansion [6]	Thermal decomposition pathways, reaction kinetics, combustion chemistry, catalytic effects [35] [23]
Computational Cost	Relatively low [6]	Higher than classical FFs, but lower than QM [2]
Key Training Data	Crystal structures, vibrational frequencies, elastic constants [6]	QM-derived reaction energies, barriers, bond dissociation energies [35]
Example EMC FF	SRT, SBFF, Boyd's potential [6]	ReaxFF-CHON, ReaxFF-lg [2]

Quantitative Performance Comparison

The selection of a force field is often dictated by the specific property under investigation. Quantitative benchmarks against experimental data and high-level QM calculations are essential for validating a force field's predictive capability.

Table 2: Accuracy Benchmarks for Force Field Predictions on Energetic Materials.

Property	Classical FF Performance	ReaxFF Performance	Notes
Crystal Density	High accuracy (deviations <0.5% for refined FFs) [6]	Good agreement with experiment	Critical for detonation performance
Mechanical Properties (Bulk/Shear Modulus)	Accurately reproduced for RDX, HMX, CL-20 [6]	Not a primary application	SRT and SBFF potentials are well-tested
Thermal Decomposition Products	Not applicable	Accurately identifies major products (e.g., N₂, H₂O, NO₂, N₂O) [35]	Provides time-resolved product evolution
Bond Dissociation Energy	Not applicable	High agreement with DFT (e.g., ~2.3% error for MgC bond) [35]	Validates the training of the force field parameters
Effect of Nano-Catalysts	Not applicable	Elucidates metal-fuel interactions (e.g., MgO bond formation with RDX) [35]	Reveals atomistic mechanisms of combustion enhancement

Experimental and Simulation Protocols

ReaxFF Molecular Dynamics for Thermal Decomposition

This protocol outlines the procedure for simulating the thermal decomposition of an EMC, such as RDX, in the presence of a metal nanoparticle catalyst using ReaxFF [35].

Objective: To investigate the initial decomposition pathways, reaction kinetics, and catalytic effects of metal nanoparticles on the thermal decomposition of RDX.

Methodology Details:

• System Construction:
  • Build a crystalline unit cell of the EMC (e.g., RDX).
  • Insert a metal nanoparticle (e.g., Mg nanoparticle) into the simulation box.
  • Use periodic boundary conditions to simulate a bulk-like environment.
• Force Field Parameterization:
  • Employ a trained ReaxFF parameter set encompassing all relevant elements (e.g., Mg/C/H/O/N for RDX-Mg systems) [35].
  • The parameter set must be trained and validated against Density Functional Theory (DFT) data, including bond dissociation energies and crystal equations of state, to ensure reliability [35].
• Simulation Setup:
  • Energy minimization: Relax the initial structure to remove high-energy contacts.
  • Equilibration: Run an NVT or NPT simulation at the target temperature (e.g., 1500-3500 K) to bring the system to a thermally equilibrated state. High temperatures are often used to accelerate reactions within computationally feasible timescales.
• Production Run:
  • Perform ReaxFF MD in the NVT ensemble for a duration sufficient to observe decomposition (tens to hundreds of picoseconds).
  • Save trajectory files at frequent intervals for post-processing.
• Trajectory Analysis:
  • Species Analysis: Use bond order analysis to track the consumption of reactants (e.g., RDX molecules) and the formation of intermediate and final products (e.g., NO₂, N₂, H₂O, MgO, Mg₃N₂).
  • Reaction Pathways: Identify key initial decomposition steps, such as N–NO₂ bond homolysis in RDX, and subsequent radical reactions.
  • Catalyst Evolution: Monitor structural changes in the metal nanoparticle, such as phase separation and the formation of oxides/nitrides [35].

The workflow for this protocol is summarized in the diagram below.

Classical MD for Mechanical and Thermal Properties

This protocol describes the use of classical MD to determine key physicochemical properties of an EMC.

Objective: To predict the mechanical properties and thermal expansion behavior of a pure EMC crystal, such as HMX or RDX.

Methodology Details:

• System Construction:
  • Build a sufficiently large supercell of the EMC crystal structure.
• Force Field Selection:
  • Choose a refined classical potential validated for the target EMC, such as the Smith-Bharadwaj (SB) potential for HMX and RDX or a potential specifically developed for TATB [6].
• Simulation Setup:
  • Energy minimization of the supercell.
  • Equilibration in the NPT ensemble at the target temperature and pressure (e.g., 300 K, 1 atm) to obtain the correct equilibrium density.
• Property Calculation:
  • Elastic Constants: Apply small strain deformations to the equilibrated cell and calculate the stress tensor. The elastic stiffness matrix (Cᵢⱼ) is derived from the linear relationship between stress and strain.
  • Thermal Expansion: Perform a series of NPT simulations across a temperature range. The volumetric thermal expansion coefficient (αᵥ) is calculated from the slope of the volume vs. temperature plot.
  • Mechanical Moduli: The bulk (K) and shear (G) moduli are computed from the calculated elastic constants using Voigt-Reuss-Hill averaging schemes [6].

Successful simulation studies require both robust software tools and well-defined molecular models. The following table lists key "research reagents" in the computational chemist's toolkit.

Table 3: Essential Resources for Simulating Energetic Molecular Crystals.

Item Name	Function/Brief Explanation	Relevance to EMC Research
ReaxFF Parameter Set (e.g., Mg/C/H/O/N)	Defines reactive interactions; trained on QM data [35].	Enables simulation of decomposition and combustion of EMCs with additives.
Wiser Bituminous Coal Model	A widely accepted molecular model for complex carbonaceous materials [23].	Serves as a reference for developing models for coal-based energetic systems or soot formation.
Classical Force Fields (e.g., SBFF, SRT-AMBER)	Pre-parameterized potentials for specific EMCs like RDX and HMX [6].	Essential for calculating crystal packing, mechanical properties, and thermal expansion.
MD Simulation Software (e.g., LAMMPS)	Open-source code that implements both classical and ReaxFF MD [14].	The computational engine for performing the dynamics simulations.
Visualization & Analysis Tools (e.g., OVITO, VMD)	Software for visualizing trajectories and analyzing results.	Critical for interpreting simulation outcomes, such as identifying reaction events and measuring structures.

The choice between classical and reactive force fields is not one of superiority but of application. Classical force fields offer high efficiency and accuracy for predicting the physicochemical properties of intact EMC crystals. In contrast, ReaxFF provides a powerful, QM-informed methodology for unraveling the complex chemical reaction networks that underpin decomposition and combustion, particularly in multi-component systems involving catalysts and oxidizers. A robust research strategy often leverages both paradigms: using classical FFs to establish initial structures and properties, and employing ReaxFF MD to probe the ensuing reactive chemistry. This combined approach provides a comprehensive atomic-scale understanding of energetic materials, from their solid-state characteristics to their gas-phase reactivity, ultimately guiding the design of next-generation materials with tailored performance and safety profiles.

In combustion chemistry research, a fundamental challenge exists: the need to model complex chemical reactions with quantum mechanical accuracy while simultaneously simulating the behavior of large molecular systems over experimentally relevant timescales. Traditional quantum mechanics (QM) methods, though valuable for electronic-level insights, suffer from severe computational cost limitations that restrict their application to small systems and short timescales, hampering our theoretical understanding of the dynamic evolution of combustion processes [2]. Classical molecular mechanics (MM) force fields, while computationally efficient, typically require predefined atomic connectivity and are fundamentally unable to simulate reactive events where bonds break and form [8]. The Reactive Force Field (ReaxFF) methodology, developed to bridge this gap, approaches reactive simulation from the classical side by casting the empirical interatomic potential within a bond-order formalism, thus implicitly describing chemical bonding without expensive QM calculations [2].

Hybrid ReaxFF/classical frameworks represent a sophisticated evolution in this field, enabling researchers to embed highly reactive regions modeled with ReaxFF within larger systems described by non-reactive classical potentials. This approach is particularly valuable for combustion systems where reactive events are often localized—such as at fuel-oxidizer interfaces, within catalyst active sites, or during initial bond cleavage in pyrolysis—while the majority of the system (solvent environments, material matrices, or bulk phases) can be efficiently treated with classical potentials. By strategically allocating computational resources, these hybrid implementations make it feasible to simulate reactive chemistry in complex, multiphase combustion environments that were previously inaccessible to computational methods [14] [8].

Theoretical Foundation: The ReaxFF Methodology

Bond-Order Formalism and Energy Composition

The ReaxFF method employs a bond-order formalism in conjunction with polarizable charge descriptions to describe both reactive and non-reactive interactions between atoms. This approach allows ReaxFF to accurately model both covalent and electrostatic interactions for a diverse range of materials relevant to combustion chemistry [2]. The fundamental concept revolves around calculating bond order (BO)—a measure of bond strength—directly from interatomic distances using an empirical formula:

[BO{ij} = BO{ij}^\sigma + BO{ij}^\pi + BO{ij}^{\pi\pi} = \exp\left[p{bo1}\left(\frac{r{ij}}{r0^\sigma}\right)^{p{bo2}}\right] + \exp\left[p{bo3}\left(\frac{r{ij}}{r0^\pi}\right)^{p{bo4}}\right] + \exp\left[p{bo5}\left(\frac{r{ij}}{r0^{\pi\pi}}\right)^{p{bo6}}\right]]

where (BO{ij}) is the bond order between atoms i and j, (r{ij}) is interatomic distance, (r0) terms are equilibrium bond lengths, and (p{bo}) terms are empirical parameters [2]. This continuous, differentiable function contains no discontinuities through transitions between σ, π, and ππ bond character, yielding a smooth potential energy surface required for calculating interatomic forces.

The total system energy in ReaxFF is composed of multiple contributions summarized by:

[E{\text{system}} = E{\text{bond}} + E{\text{over}} + E{\text{angle}} + E{\text{tors}} + E{\text{vdWaals}} + E{\text{Coulomb}} + E{\text{Specific}}]

where (E{\text{bond}}) is the energy associated with forming bonds between atoms, (E{\text{angle}}) and (E{\text{tors}}) are the energies associated with three-body valence angle strain and four-body torsional angle strain, and (E{\text{over}}) is an energy penalty preventing over-coordination of atoms based on atomic valence rules [2]. (E{\text{Coulomb}}) and (E{\text{vdWaals}}) represent electrostatic and dispersive contributions calculated between all atoms, while (E_{\text{Specific}}) encompasses system-specific terms such as lone-pair, conjugation, hydrogen binding, and C₂ corrections that may be required for specific combustion systems [2].

Charge Equilibration and Environmental Transfers

A critical feature of ReaxFF is its dynamic charge distribution scheme, which calculates partial atomic charges at each molecular dynamics step using charge equilibration methods (such as QEq or EEM) [8]. This approach, while essential for accurately modeling charge transfer during reactions, introduces significant computational complexity and becomes particularly challenging in hybrid implementations where charge consistency must be maintained across different force field regimes [36]. The method demonstrates remarkable transferability across phases, treating elements with the same mathematical formalism regardless of whether they are in gas, liquid, or solid phases—a crucial capability for combustion systems involving multiphase interactions [2].

Implementation Architectures for Hybrid Frameworks

The ReaxFF/AMBER Framework

The ReaxFF/AMBER hybrid implementation represents a significant advancement for studying local reactive events in large biomolecular or condensed-phase systems at a fraction of the computational costs of QM/MM models [8]. This framework introduces ReaxFF capabilities to capture bond breaking and formation within the established AMBER molecular dynamics package, leveraging AMBER's well-validated non-reactive force fields for the majority of the system while reserving ReaxFF for chemically active regions. The implementation allows researchers to exploit AMBER's extensive suite of analysis tools, including its potential of mean force (PMF) capabilities for mapping reaction profiles of organic transformations in solution [8].

In practice, the ReaxFF/AMBER framework partitions the system into reactive and non-reactive regions, with careful attention to managing the boundary between these regions. The implementation requires specialized parameters for handling interactions across different force field domains, particularly for the charge equilibration scheme which must be constrained to prevent unphysical charge transfer between regions [8] [36]. Validation studies, such as simulations of a benzene molecule solvated in water, have demonstrated the framework's ability to maintain energy conservation and proper thermodynamic behavior across the hybrid interface [8].

Integration with LAMMPS and PuReMD

The Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) provides another platform for hybrid ReaxFF implementations, offering extensive capabilities for simulating complex material systems with various boundary conditions and coupling schemes [2]. Similarly, the Purdue Reactive Molecular Dynamics (PuReMD) code has been optimized specifically for ReaxFF simulations, providing efficient algorithms for bond order calculations and charge equilibration [2]. However, practical implementations reveal significant challenges, particularly regarding the integration of ReaxFF with other potentials in hybrid pair styles [36].

As noted in discussions of LAMMPS implementations, "The ReaxFF model by its very nature is not a good choice for use in a hybrid pair style" due primarily to inconsistencies in charge equilibration when atoms modeled by different potentials interact [36]. These challenges manifest as unphysical energy changes when non-ReaxFF atoms are introduced into the system, even when those atoms carry zero charge, suggesting fundamental incompatibilities in how electrostatic interactions are handled across different force field regimes [36].

Table 1: Comparison of Hybrid ReaxFF Implementation Platforms

Platform	Advantages	Limitations	Typical Applications
ReaxFF/AMBER	Access to AMBER's MM force fields and PMF capabilities; Simplified setup for biological systems	Limited to AMBER-compatible force fields; Boundary management challenges	Solvated organic reactions; Enzymatic catalysis; Biomaterial interfaces
LAMMPS	High parallel efficiency; Extensive boundary condition options; Multiple thermostat/barostat choices	Significant charge equilibration issues in hybrid styles [36]; Complex parameter synchronization	Materials interfaces; Nanocomposites; Solid-fluid reactions
PuReMD	Optimized specifically for ReaxFF; Efficient neighbor listing; Fast charge equilibration	Limited MM force field options; Fewer hybrid implementation examples	Combustion chemistry; Reactive nanoparticle synthesis; Catalyst screening

Performance Benchmarks and Comparative Analysis

Computational Efficiency and System Scaling

The computational cost of ReaxFF simulations, while significantly lower than QM methods, remains substantially higher than classical non-reactive force fields due to the complex bond order calculations and charge equilibration required at each time step [8]. This computational differential becomes particularly important in hybrid implementations, where the ReaxFF region typically dictates the overall simulation time step and limits parallel efficiency due to load imbalance between processors handling reactive and non-reactive regions.

Recent developments in reactive simulation methods have introduced alternative approaches that aim to address these computational challenges. The Reactive INTERFACE Force Field (IFF-R), which replaces harmonic bond potentials with reactive, energy-conserving Morse potentials, claims to offer approximately 30 times faster performance than ReaxFF while maintaining compatibility with existing force fields for organic and inorganic compounds [3]. This performance advantage potentially positions IFF-R as a competitive alternative for certain classes of combustion simulations where complex bond formation pathways are less critical.

Table 2: Performance Comparison of Reactive Simulation Methods

Method	Computational Cost Relative to MM	Key Strengths	Key Limitations
QM/MM	100-10,000x	High accuracy for reaction barriers; Electronic structure detail	Limited to small QM regions; Picosecond timescales
ReaxFF	10-100x	Full reactivity; Bond formation/breaking; Charge transfer	High parameterization burden; Complex functional form
IFF-R	3-30x	Simple parameterization; Backward compatibility; Rapid simulations	Limited bond formation capabilities; Template-based reactions [3]
Classical MM	1x (reference)	Extremely fast; Large spatiotemporal scales	No reactivity; Fixed bonding topology

Accuracy and Transferability in Combustion Systems

The performance of ReaxFF potentials has been extensively evaluated for various material classes relevant to combustion chemistry. For sp² carbon systems—including graphene, carbon nanotubes, and fullerenes that appear in soot formation pathways—multiple ReaxFF parameterizations have demonstrated reasonable agreement with DFT calculations and experimental data for structural properties, energetics, and mechanical behavior under strain [5]. However, systematic evaluations have identified specific limitations, including underestimation of Young's modulus, overestimation of Poisson's ratio, and anomalous stress-strain curve behavior in certain deformation regimes [5].

For energetic materials relevant to combustion and propulsion applications, ReaxFF has shown particular utility in modeling complex decomposition pathways and sensitivity properties. The method has been successfully applied to nitramine compounds including RDX and HMX, providing insights into initiation mechanisms and reaction kinetics under various temperature and pressure conditions [2] [6]. The transferability of ReaxFF across phases enables simulations of solid-state decomposition, melt-phase reactions, and gas-phase oxidation processes within a consistent theoretical framework [2].

Practical Implementation Guide

System Setup and Region Partitioning

Successful implementation of hybrid ReaxFF/classical simulations requires careful system setup and region partitioning. The following workflow outlines the critical steps for preparing a hybrid simulation:

Figure 1: Workflow for hybrid ReaxFF simulation setup.

• System Preparation: Begin with a fully atomistic model of the complete system, ensuring proper solvation, ionization, and initial geometry. For combustion systems, this may include fuel clusters, oxidizer molecules, and catalyst surfaces in appropriate periodic boundary conditions.

• Region Partitioning: Identify and tag atoms belonging to reactive (ReaxFF) and non-reactive (MM) regions. The reactive region should include all atoms directly involved in bond-breaking/forming events plus a sufficient buffer zone (typically 5-10 Å) to ensure proper electrostatic and van der Waals interactions across the boundary.

• Parameter Assignment: Assign ReaxFF parameters to reactive regions and classical MM parameters to non-reactive regions. Special attention must be paid to the boundary region, where specialized coupling parameters may be required to prevent unphysical forces or energy drift.

• Equilibration: Perform careful equilibration of both regions, typically beginning with energy minimization followed by gradual heating and density equilibration. Constraint algorithms may be necessary to maintain proper geometry at the interface during this phase.

• Production Simulation: Conduct the production simulation with appropriate thermodynamic ensemble (NVT, NPT) settings, ensuring that data collection includes detailed bonding information and reaction analysis specific to the ReaxFF region.

Research Reagent Solutions: Computational Tools

Table 3: Essential Software Tools for Hybrid ReaxFF Simulations

Tool Name	Function	Implementation Notes
AMBER	Molecular dynamics package with ReaxFF/MM hybrid support	Provides fix qeq/reaxff for charge equilibration; PMF capabilities for reaction profiling [8]
LAMMPS	Open-source MD simulator with ReaxFF implementation	Supports hybrid pair styles but with noted limitations for ReaxFF combinations [36]
PuReMD	Optimized reactive molecular dynamics code	Efficient for pure ReaxFF simulations; Limited hybrid capabilities [2]
ADF	Quantum chemistry package with ReaxFF module	Parametrization tools for developing new ReaxFF parameter sets [37]
Materials Studio	Commercial modeling environment with ReaxFF	User-friendly interface for simulation setup and analysis [2]

Boundary Management and Charge Consistency

The most significant technical challenge in hybrid ReaxFF implementations involves managing the boundary between reactive and non-reactive regions, particularly regarding charge equilibration. The following diagram illustrates the charge equilibration challenge in hybrid systems:

Figure 2: Charge equilibration challenges and solutions in hybrid ReaxFF.

Practical solutions to these challenges include:

• Charge Constraint Algorithms: Application of harmonic constraints or Lagrange multipliers to fix atomic charges in the MM region while allowing dynamic charge equilibration in the ReaxFF region.

• Boundary Buffering: Implementation of a sufficiently large buffer region between reactive and non-reactive zones to minimize unphysical boundary effects on the electronic structure.

• Neutral Group Partitioning: Division of the system into chemically intuitive neutral groups that span the boundary region, ensuring that charge transfer occurs primarily within recognizable chemical motifs.

• Extended Lagrangian Approaches: Implementation of extended Lagrangian methods for charge equilibration that improve energy conservation and numerical stability across force field boundaries.

Applications in Combustion Chemistry

Pyrolysis and Oxidation Pathways

ReaxFF has been successfully deployed to gain fundamental insights into pyrolysis and oxidation of various gas, liquid, and solid fuels, revealing detailed energy changes and chemical pathways at the atomistic level [14]. Hybrid implementations enable researchers to embed reactive fuel molecules within larger non-reactive environments, such as fuel droplets in air or solid fuel particles in oxidizer streams. These simulations have provided unprecedented details of initial bond cleavage events, radical formation, and subsequent oxidation pathways that initiate and sustain combustion processes.

For hydrocarbon combustion systems, the CHO-2008 parameterization has demonstrated particular utility in modeling complex reaction networks involving hundreds of simultaneous chemical transformations [5]. The hybrid approach allows focusing computational resources on the reactive core while treating the bulk environment with efficient classical potentials, enabling simulation timescales relevant to actual combustion processes.

Soot Formation and Nanoparticle Synthesis

The complex physicochemical dynamic processes involved in soot formation and flame synthesis of nanoparticles have been made plainly visible from an atomistic perspective through ReaxFF simulations [14]. Hybrid frameworks are particularly valuable for these applications, as they allow modeling of the gas-phase reaction chemistry that generates molecular precursors while simultaneously capturing the surface growth and aggregation processes that lead to nanoparticle formation.

Studies of sp² carbon systems, including graphene, fullerenes, and carbon nanotubes, have benefited from specialized ReaxFF parameterizations that accurately reproduce the structural, energetic, and mechanical properties of these materials [5]. The hybrid approach enables researchers to simulate the interface between growing carbon nanostructures and their reactive environments at unprecedented scales, providing insights into growth mechanisms and structural defects that influence material properties.

Catalytic Combustion and Surface Reactions

Heterogeneous combustion systems involving catalytic surfaces represent an ideal application for hybrid ReaxFF/classical frameworks. In these systems, the reactive events are highly localized at the catalyst surface and the immediate interface with the gas phase, while the bulk catalyst material and the majority of the gas phase can be treated with classical potentials. This partitioning enables atomistic modeling of catalytic combustion processes with realistic surface geometries and gas-phase compositions.

ReaxFF has been applied to various catalyst systems relevant to combustion, including metal oxides, supported metal nanoparticles, and zeolitic materials [2] [37]. The method's ability to model charge transfer and bond rearrangement at catalyst surfaces provides unique insights into adsorption energies, reaction barriers, and product distributions that control combustion efficiency and emissions.

Future Perspectives and Development Directions

The field of hybrid reactive molecular dynamics continues to evolve rapidly, with several promising directions for future development. Machine learning approaches are being integrated with traditional force fields to improve accuracy while maintaining computational efficiency, potentially offering solutions to current transferability limitations [14]. Multiscale modeling frameworks that couple ReaxFF with coarse-grained methods and continuum approaches are extending the accessible length and timescales for combustion simulation [14].

Recent developments in explicit electron treatment within ReaxFF (eReaxFF) enable more accurate modeling of redox processes and electron transfer reactions that are fundamental to electrochemical combustion and battery applications [37]. Additionally, advanced sampling methods and accelerated dynamics techniques are being integrated with hybrid frameworks to overcome the timescale limitations inherent in molecular dynamics, particularly for rare events such as ignition processes or catalyst deactivation [8] [37].

As computational power continues to grow and algorithms become more sophisticated, hybrid ReaxFF/classical frameworks are poised to become increasingly central to combustion chemistry research, providing a unique bridge between electronic-scale accuracy and device-scale simulation that ultimately enables predictive design of combustion systems with optimized efficiency and reduced environmental impact.

Navigating Challenges and Best Practices in Reactive Molecular Dynamics

In the realm of computational combustion research, the selection between classical force fields (FFs) and reactive force fields (ReaxFF) represents a critical methodological crossroads, with significant implications for predictive accuracy and simulation scalability. Classical FFs, which rely on fixed bonding relationships and simplified functional forms, offer computational efficiency but fall short in simulating reactive processes fundamental to combustion, such as bond breaking and formation [2]. In contrast, ReaxFF employs a bond-order formalism that enables dynamic connectivity, allowing it to model complex reaction pathways and intermediate species characteristic of hydrocarbon oxidation and other combustion phenomena [14] [2]. However, this enhanced capability comes with increased complexity, as ReaxFF typically contains hundreds of parameters compared to the 10-100 parameters found in classical FFs, creating substantial challenges in parametrization and transferability [38] [39]. These challenges form the core focus of this technical guide, which examines the prevalent pitfalls in force field application and outlines systematic approaches to overcome them within the specific context of combustion chemistry research.

Fundamental Differences: Classical FFs vs. ReaxFF

The divergence between classical and reactive force fields stems from their underlying mathematical frameworks and treatment of chemical bonding. Understanding these foundational differences is essential for recognizing the specific limitations and failure modes associated with each approach.

Table 1: Core Functional Components of Classical and Reactive Force Fields

Energy Component	Classical Force Fields	Reactive Force Fields (ReaxFF)
Bonding	Predefined harmonic/quadratic springs	Bond-order derived from interatomic distances
Angles	Harmonic angle bending	Bond-order dependent angle strain
Torsions	Cosine-based dihedral potentials	Bond-order dependent torsion energy
Non-bonded	Fixed partial charges (Coulomb), Lennard-Jones	Charge equilibration (QEq), bond-order corrected vdWaals
Reactivity	Not available	Over-coordination penalty, conjugation, lone pairs
Transferability	Limited to trained molecular environments	Higher, but requires element-specific parametrization

Classical force fields calculate system energy using simplified interatomic potential functions well-suited for nonreactive interactions, such as bond stretching represented by harmonic functions and dispersion represented by van der Waals potentials [2]. This formulation depends on predefined connectivity, rendering them incapable of simulating chemical reactions where bonds break and form [2]. The energy expression for classical FFs typically follows a simple functional relationship that maps system energy to atomic positions and charges, with parameters often possessing clear physical meanings and high interpretability [38].

ReaxFF addresses this limitation through a bond-order formalism where bond order is empirically calculated from interatomic distances [2]. The total energy in ReaxFF incorporates multiple contributions:

[E{\text{system}} = E{\text{bond}} + E{\text{over}} + E{\text{angle}} + E{\text{tors}} + E{\text{vdWaals}} + E{\text{Coulomb}} + E{\text{Specific}}]

where bond order ((BO_{ij})) is calculated as:

[BO{ij} = BO{ij}^\sigma + BO{ij}^\pi + BO{ij}^{\pi\pi} = \exp\left[p{bo1}\left(\frac{r{ij}}{ro^\sigma}\right)^{p{bo2}}\right] + \exp\left[p{bo3}\left(\frac{r{ij}}{ro^\pi}\right)^{p{bo4}}\right] + \exp\left[p{bo5}\left(\frac{r{ij}}{ro^{\pi\pi}}\right)^{p{bo6}}\right]]

This approach enables ReaxFF to dynamically describe reactivity while maintaining a differentiable potential energy surface necessary for force calculations [2]. The bond-order formalism allows ReaxFF to accommodate long-distance covalent interactions characteristic of transition state structures, enabling reasonable prediction of reaction barriers [2].

Diagram 1: The ReaxFF parametrization workflow showing the iterative refinement process based on quantum mechanical data and experimental validation.

Critical Challenges in Force Field Transferability

Transferability—the ability of a force field to accurately describe systems and conditions beyond those explicitly included in its training set—represents a fundamental challenge in computational combustion science. The specialized parametrization of force fields for specific molecular environments often limits their applicability to broader chemical contexts.

The Branching Problem: Combustion vs. Aqueous Descriptions

A prominent manifestation of transferability limitations in ReaxFF is the existence of separate parameter branches optimized for different chemical environments. The combustion branch focuses on accurately describing water as a gas-phase molecule and hydrocarbon oxidation processes, while the aqueous branch targets liquid-phase chemistry and solvation effects [11]. This division creates significant pitfalls when applying force fields outside their intended domain. For instance, the O/H parameters differ substantially between these branches, leading to potentially unrealistic results when used interchangeably [11].

Elemental and Phase Transferability Gaps

Transferability failures frequently occur when force fields encounter elements or phases not adequately represented in their training data. For example, the ReaxFF description for H/C/O compounds developed by Chenoweth et al. (2008) performs excellently in combustion chemistry but demonstrates notable inaccuracies in describing the mechanical properties of graphite and diamond [39]. Similarly, force fields parameterized for gas-phase interactions may fail to capture condensed-phase behavior, limiting their utility for simulating droplet combustion or soot formation processes [14].

Functional Form Limitations

The mathematical structure of force fields inherently constrains their transferability. Classical FFs with fixed bonding terms cannot adapt to changes in hybridization or bond order that occur during combustion reactions [6]. While ReaxFF's bond-order formalism provides greater flexibility, it still relies on parameterized functions that may not adequately capture all relevant potential energy surfaces, particularly for transition metals or exotic chemistries encountered in catalytic combustion [2] [39].

Table 2: Common Transferability Pitfalls and Their Consequences in Combustion Simulations

Pitfall Category	Specific Examples	Impact on Combustion Research
Elemental Coverage	Missing parameters for new elements (e.g., metals in catalysts)	Inaccurate modeling of catalytic combustion, nanoparticle synthesis
Phase Transferability	Aqueous branch parameters applied to gas-phase combustion	Incorrect reaction barriers, unrealistic intermediate stability
Reactivity Limits	Inadequate description of transition states for novel reactions	Wrong prediction of ignition timing, combustion efficiency
Property-specific Errors	Accurate energies but poor mechanical property prediction	Unreliable modeling of material stress in combustion environments

Parametrization represents the most nuanced aspect of force field development, where methodological choices can introduce systematic errors that propagate through subsequent simulations. The training data selection directly determines a force field's strengths and limitations, with insufficient coverage of chemical space leading to parametric insufficiencies.

Training Set Deficiencies

A fundamental source of parametrization error stems from incomplete or unrepresentative training data. Traditional parametrization of intermolecular interactions has primarily relied on pure component properties such as liquid density (ρl) and enthalpy of vaporization (ΔHvap) [40]. This approach, while computationally efficient, fails to adequately capture mixture interactions essential for predicting combustion behavior in complex fuel blends [40]. For example, force fields accurately reproducing pure component properties may still incorrectly predict azeotropic behavior in binary mixtures, leading to errors in vapor-liquid equilibrium calculations relevant to fuel evaporation [40].

The limited availability of high-quality experimental data further compounds these challenges. The NIST ThermoML Archive contains approximately 500 ΔHvap data points compared to over 60,000 measurements of pure densities, creating an inherent imbalance in parametrization constraints [40]. This data scarcity problem is particularly acute for reactive intermediates and transition states central to combustion chemistry, which are difficult to characterize experimentally yet crucial for accurate reaction kinetics prediction.

Optimization Method Limitations

The optimization algorithms used in parameter fitting introduce their own limitations. Traditional sequential one-parameter parabolic interpolation (SOPPI) methods train each parameter sequentially, resulting in excessive time consumption and susceptibility to local minima [39]. While genetic algorithms (GA) can avoid local optima, they involve complex selection and crossover operators and suffer from premature convergence issues [39]. Recent approaches combining simulated annealing (SA) and particle swarm optimization (PSO) with a concentrated attention mechanism (CAM) show improved efficiency but still require careful implementation to avoid parametric compensation effects [39].

The high dimensionality of the optimization problem presents additional challenges. With ReaxFF containing hundreds of parameters, the risk of overfitting to the training set increases significantly, potentially reducing transferability to unseen systems [39]. This problem is exacerbated when training data is sparse or lacks diversity in chemical environments, a common scenario in combustion chemistry where certain reaction classes may be underrepresented.

Methodological Solutions for Enhanced Transferability

Advanced Parametrization Strategies

Addressing parametrization errors requires methodological innovations that expand beyond traditional approaches. Incorporating binary mixture properties such as densities (ρl(x)) and enthalpies of mixing (ΔHmix(x)) provides enhanced constraints on intermolecular parameters, particularly for capturing solute-solvent interactions critical to spray combustion and fuel droplet evaporation [40]. This approach directly addresses transferability gaps by ensuring force fields accurately represent the cross-interactions between different molecular species present in practical combustion systems.

Diagram 2: Evolution of parametrization methodologies showing the progression from classical pre-training to machine learning potentials with enhanced validation.

Multi-stage training protocols that separate physical preconditioning from chemical refinement offer promising avenues for improving parametrization robustness. Recent research demonstrates that pre-training machine-learned potentials using classical force field data followed by fine-tuning with high-quality ab initio data enhances stability and transferability [34]. This approach leverages the comprehensive phase space coverage of classical FFs while maintaining the quantum mechanical accuracy needed for reactive simulations, effectively addressing the "PES holes" problem where ML potentials encounter unphysical configurations during dynamics [34].

Systematic Validation Frameworks

Comprehensive validation against diverse data types represents a critical safeguard against parametrization errors. Beyond traditional benchmarks against energies and structures, validation should include mechanical properties, vibrational spectra, and transport coefficients to ensure broad transferability [6]. For combustion applications, particular emphasis should be placed on validating reaction barriers and rates against experimental measurements or high-level quantum calculations, as these directly impact predictive accuracy for ignition timing and combustion efficiency.

The development of branch-consistent parameters that maintain accuracy across different chemical environments (e.g., gas, liquid, interface) represents an ongoing challenge. Regular cross-validation between combustion and aqueous branches can identify systematic deviations and guide parametrization improvements [11]. Additionally, establishing community-wide standards for force field validation specific to combustion systems would facilitate more systematic assessment of transferability limitations.

Emerging Approaches: Machine Learning Potentials

Machine learning potentials (MLPs) represent a paradigm shift in force field development, offering the potential to overcome fundamental limitations of both classical and reactive FFs. MLPs trained on large-scale quantum mechanical datasets can achieve density functional theory (DFT) level accuracy at a fraction of the computational cost, enabling high-fidelity simulations of complex combustion processes [13].

Recent developments such as the EMFF-2025 model for C/H/N/O-based energetic materials demonstrate MLPs' capability to predict mechanical properties and decomposition characteristics with quantum accuracy [13]. These models leverage transfer learning strategies, where a pre-trained base model is fine-tuned with minimal system-specific data, significantly reducing the computational expense associated with traditional parametrization [13]. For combustion research, this approach enables rapid adaptation to novel fuel molecules while maintaining accuracy across different reaction classes.

Foundation models like MACE demonstrate remarkable performance in modeling complex proton transport phenomena, achieving quantitative agreement with AIMD simulations after minimal fine-tuning with as little as 1 picosecond of reference data [41]. This capability is particularly valuable for combustion systems involving proton-coupled electron transfer or acid-catalyzed reactions, where traditional FFs often struggle with accurate description.

Table 3: Comparison of Force Field Types for Combustion Applications

Characteristic	Classical FFs	ReaxFF	Machine Learning Potentials
Computational Cost	Low	Medium	Medium-High (training), Low (inference)
Reactivity Description	None	Bond-order formalism	DFT-level accuracy
Parametrization Effort	Low (10-100 parameters)	High (100+ parameters)	Very High (data-driven)
Transferability	Limited to trained systems	Moderate, with branches	High with sufficient training data
Best Applications	Non-reactive molecular dynamics, transport properties	Reaction network exploration, pyrolysis	High-accuracy kinetics, complex interfaces

The Scientist's Toolkit: Essential Methodologies

Research Reagent Solutions

Table 4: Essential Computational Tools for Force Field Development and Validation

Tool Category	Specific Examples	Primary Function	Application in Combustion Research
Quantum Chemistry Codes	VASP, CP2K, Gaussian, Q-Chem	Generate training data	Reference energies, forces, reaction barriers
ReaxFF Parametrization	Standalone ReaxFF, SOPPI, SA+PSO+CAM	Force field optimization	Develop system-specific parameters
MD Simulation Engines	LAMMPS, PuReMD, AMS	Dynamics simulations	Pyrolysis, oxidation, catalyst studies
Validation Tools	Custom scripts, ThermoML	Property calculation	Compare with experimental measurements
MLP Frameworks	DEEPMD, MACE, EMFF-2025	Machine learning potentials	High-accuracy reactive simulations

Experimental Protocols for Force Field Validation

Protocol 1: Combustion Reaction Kinetics Validation

• System Setup: Create simulation cells containing fuel-oxidizer mixtures at relevant combustion conditions (temperatures: 300-2500K, pressures: 0.1-100 atm).
• Simulation Procedure: Perform NVT molecular dynamics simulations using candidate force field, monitoring reaction progression.
• Data Collection: Track formation rates of key intermediates (e.g., OH, HO2, CH2O) and final products (CO2, H2O).
• Benchmarking: Compare temporal concentration profiles against experimental measurements from shock tube studies or flow reactor data.
• Validation Metrics: Calculate root-mean-square deviations for ignition delay times and major species concentrations.

Protocol 2: Transferability Assessment Across Phases

• Multi-phase Systems: Construct simulation cells containing gas-liquid interfaces (e.g., fuel droplets in air).
• Property Evaluation: Calculate surface tension, evaporation rates, and interfacial mass transfer coefficients.
• Comparative Analysis: Benchmark against experimental interfacial data and high-level quantum mechanics/molecular mechanics (QM/MM) calculations.
• Error Quantification: Establish acceptable deviation thresholds for phase-equilibrium properties (<5% for vapor-liquid equilibrium).

The challenges of force field transferability and parametrization represent significant but addressable barriers to accurate computational combustion research. Navigating these pitfalls requires careful consideration of force field limitations, systematic validation across multiple properties, and adoption of emerging methodologies that enhance parametric robustness. The integration of machine learning approaches with traditional force field development holds particular promise for overcoming current limitations, potentially enabling a new generation of combustion models with unprecedented accuracy and transferability. As computational methods continue to evolve, the rigorous assessment of parametric uncertainties and systematic addressing of transferability gaps will remain essential for reliable simulation of complex combustion phenomena across the diverse range of conditions encountered in practical energy systems.

The accurate simulation of combustion chemistry is a cornerstone in the development of advanced propulsion systems and cleaner energy technologies. Molecular dynamics (MD) simulations provide a powerful tool for probing these complex reactive processes at the atomic scale. For combustion researchers, the choice of force field—the mathematical model describing atomic interactions—represents a critical decision point balancing computational cost against physical accuracy. This technical guide examines the performance landscape of reactive force fields, with a specific focus on documented failures of the Reactive Force Field (ReaxFF) methodology against the emerging context of classical and machine learning alternatives. Framed within a broader thesis evaluating ReaxFF against classical force fields for combustion chemistry, this analysis synthesizes recent benchmarking studies to provide researchers with a clear assessment of capabilities and limitations.

The ReaxFF method has been extensively applied to simulate combustion processes, offering a computationally efficient alternative to quantum mechanical (QM) methods by employing bond-order formalism to enable bond breaking and formation. [28] However, its parametrized nature raises fundamental questions about transferability and accuracy across different chemical systems. Concurrently, classical force fields augmented with reactive capabilities and machine learning potentials have emerged as promising alternatives, claiming significant improvements in both accuracy and computational efficiency. [3] [13] This whitepaper provides a comprehensive technical examination of documented ReaxFF failures, quantitative benchmarking data, and methodological frameworks for evaluating force field performance in combustion systems.

Force Field Architectures in Combustion Research

Methodological Spectrum and Key Differentiators

Table 1: Classification of Force Field Methods for Combustion Chemistry

Force Field Type	Reactive Capability	Computational Cost	Parametrization Complexity	Key Characteristics
Classical Harmonic	Non-reactive	Low	Low (10-100 parameters)	Harmonic bond potentials; unsuitable for bond breaking/formation [3] [38]
Reactive (ReaxFF)	Bond-order dependent	Moderate-High	High (100+ parameters)	Complex functional forms; parametrized potential [38] [42]
Morse Potential (IFF-R)	Enabled via Morse potentials	Moderate (~30x faster than ReaxFF)	Moderate (3 parameters per bond)	Clean harmonic replacement; interpretable parameters [3]
Machine Learning (ML)	Inherent via ML architecture	Variable (DFT-level accuracy)	Very High (training data intensive)	Universal MLIPs (uMLIPs) trained on diverse datasets [13] [43]

The ReaxFF Formalism and Combustion Applications

ReaxFF employs a bond order formalism to enable bond dissociation and formation during molecular dynamics simulations. The force field utilizes bond-order-dependent polarization charges to model both reactive and non-reactive atomic interactions, with parameter sets typically derived from quantum mechanical calculations for specific chemical systems. [28] This approach has been applied to various combustion scenarios, including hydrocarbon fuel decomposition and nanocatalyst-enhanced combustion processes. [28] For instance, Pt-graphene catalyzed combustion of high-energy-density fuels has been investigated using ReaxFF to analyze ignition delay and combustion efficiency improvements. [28]

However, the method's fundamental limitation lies in its empirical parametrization, which may not fully capture the complexity of potential energy surfaces (PES) across diverse reaction pathways. As noted in benchmarking studies, "ReaxFF still struggles to achieve the accuracy of density functional theory (DFT) in describing reaction potential energy surfaces (PES), particularly when applied to new molecular systems." [13] This shortcoming becomes particularly problematic in combustion environments where complex reaction networks and radical interactions dominate the chemical dynamics.

Documented ReaxFF Failures: Quantitative and Qualitative Evidence

Systematic Benchmarking in Hydrogen Combustion

A rigorous benchmark study evaluating ReaxFF performance for hydrogen combustion systems revealed both quantitative and qualitative failures when compared to higher-level quantum mechanical methods. [42] The investigation analyzed several common parametrizations of ReaxFF, demonstrating instances where these potentials failed to accurately describe reactive events fundamental to hydrogen combustion chemistry. The study employed quantum mechanical calculations as reference data, comparing reaction energetics, kinetics, and dynamics predicted by ReaxFF against these high-accuracy benchmarks.

Table 2: Documented ReaxFF Failure Modes in Combustion Systems

Failure Type	Specific Manifestation	System Impact	Experimental Evidence
Quantitative	Inaccurate reaction energy barriers	Erroneous reaction rates and kinetics	Deviation from QM references in H₂ combustion [42]
Quantitative	Incorrect force predictions	flawed atomic trajectories and dynamics	MAE exceeding acceptable thresholds [33]
Qualitative	Misrepresentation of reaction pathways	Fundamentally incorrect mechanism prediction	Altered product distribution in hydrocarbon oxidation [13]
Qualitative	Failure to describe transition states	Inaccurate reaction coordinate mapping	Discrepancies in ignition delay predictions [42]

Performance Limitations in Complex Combustion Environments

Beyond fundamental benchmarking, ReaxFF demonstrates limitations in applied combustion scenarios. While the force field has shown utility in modeling nanocatalyst-enhanced combustion—such as Pt-graphene hybrids in exo-tetrahydrodicyclopentadiene (exo-THDCPD) combustion [28]—its generalizability across diverse fuel types and conditions remains constrained. The parametrization-specific nature of ReaxFF means that parameters optimized for one class of compounds (e.g., hydrocarbons) may perform poorly for others (e.g., nitrogen-containing energetic materials), necessitating system-specific reparameterization. [13]

This transferability challenge is particularly acute in combustion systems where complex mixtures and evolving intermediate species create a constantly changing chemical environment. The inability of fixed parametrizations to adapt to these dynamic conditions represents a fundamental limitation of the ReaxFF approach compared to more flexible machine learning alternatives that can potentially learn and adapt to diverse chemical environments. [13] [43]

Experimental Benchmarking Methodologies

Workflow for Force Field Validation

Reference Data Generation and Property Calculations

The foundation of reliable force field benchmarking lies in the generation of high-quality quantum mechanical reference data. Methodologies should employ high-level ab initio methods such as CCSD(T), MP2, or carefully validated density functionals for calculating reference reaction energies, barriers, and forces. [3] [33] For combustion systems, particular attention should be paid to radical reaction pathways, transition state geometries, and bond dissociation curves, as these represent critical failure points for empirical potentials.

Key properties for benchmarking include:

• Bond Dissociation Curves: Comparing energy versus bond length profiles against QM references across the entire reaction coordinate. [3]
• Reaction Energy Barriers: Calculating activation energies for elementary combustion steps (e.g., H-abstraction, radical recombination). [33]
• Atomic Forces: Evaluating force prediction accuracy under various molecular configurations. [13]
• Reaction Pathways: Mapping complete potential energy surfaces for representative combustion reactions. [38]

The JAX-ReaxFF framework offers a modern approach for systematic ReaxFF optimization using gradient-based algorithms, though fundamental limitations in the functional form may persist even after parametrization improvements. [33]

Emerging Alternatives and Comparative Performance

Morse Potential Approaches

The Reactive INTERFACE Force Field (IFF-R) represents a significant methodological alternative to bond-order potentials by replacing classical harmonic bond potentials with reactive, energy-conserving Morse potentials. [3] This approach maintains compatibility with existing biomolecular force fields (CHARMM, AMBER, OPLS-AA) while enabling bond dissociation through three interpretable parameters per bond type. The Morse potential formulation provides a more physically justified description of bond dissociation compared to harmonic approximations, with demonstrated applications in bond breaking for polymers, carbon nanostructures, and composite materials. [3]

Critically, IFF-R claims a 30-fold speed increase compared to ReaxFF while maintaining accuracy relative to non-reactive force fields. [3] This performance advantage stems from the simpler functional form that avoids the complex bond-order calculations inherent to ReaxFF. The method also supports bond formation through template-based approaches such as the REACTER toolkit, providing comprehensive reactivity within a more computationally efficient framework. [3]

Machine Learning Force Fields

Machine learning potentials represent a paradigm shift in force field development, offering DFT-level accuracy with significantly reduced computational cost compared to direct quantum mechanical calculations. [13] Frameworks such as Deep Potential (DP) and universal machine learning interatomic potentials (uMLIPs) have demonstrated remarkable capabilities in modeling complex reactive processes while maintaining first-principles accuracy. [13] [43]

The EMFF-2025 model, for instance, provides a specialized neural network potential for C, H, N, O-based high-energy materials, predicting structures, mechanical properties, and decomposition characteristics with DFT-level accuracy. [13] Similarly, universal MLIPs like M3GNet, MACE, and ORB models have shown excellent performance across diverse chemical systems, though accuracy can degrade for lower-dimensional structures. [43] These models achieve energy errors below 10 meV/atom and errors in atomic positions of 0.01–0.02 Å across diverse dimensionalities, demonstrating sufficient accuracy for most combustion applications. [43]

Performance Comparison Across Methodologies

Table 3: Comparative Performance Metrics for Combustion Force Fields

Methodology	Energy Error Range	Force Error Range	Computational Speed	Combustion Application Readiness
ReaxFF	Variable; system-dependent	Variable; system-dependent	Moderate	Mature but with documented failures [42]
IFF-R (Morse)	Maintains base FF accuracy	Maintains base FF accuracy	~30x faster than ReaxFF	Emerging validation [3]
ML Potentials	MAE < 0.1 eV/atom [13]	MAE < 2 eV/Å [13]	DFT-level accuracy at lower cost	Rapidly advancing [13] [43]

Research Reagent Solutions: Computational Tools for Combustion Chemistry

Table 4: Essential Computational Tools for Combustion Force Field Research

Tool/Resource	Function	Application Context
ReaxFF Parametrizations	System-specific reactive MD	Combustion pathway exploration [28]
IFF-R Parameters	Morse potential implementation	Bond breaking in polymers, composites [3]
JAX-ReaxFF Framework	ReaxFF optimization platform	Force field parameter refinement [33]
Deep Potential Generator	ML potential training	Developing system-specific MLFFs [13]
Universal MLIPs (MACE, ORB)	Pre-trained ML potentials	Transferable reactive simulations [43]
QM Reference Data	High-level calculation outputs	Force field benchmarking and training [42]

This technical assessment demonstrates that while ReaxFF has provided valuable insights into combustion mechanisms, its documented quantitative and qualitative failures necessitate careful validation for any new application system. The emerging landscape of force field methodologies offers promising alternatives, with Morse potential approaches (IFF-R) providing interpretability and speed advantages, while machine learning potentials deliver unprecedented accuracy with improving computational efficiency.

For combustion researchers, the selection of an appropriate force field must balance system specificity, computational resources, and accuracy requirements. Benchmarking against high-level quantum mechanical references remains essential, particularly for reactions central to combustion mechanisms where ReaxFF has demonstrated fundamental limitations. The continued development of transferable, accurate, and computationally efficient potentials will increasingly enable reliable simulation of complex combustion environments across multiple scales and chemical systems.

In the pursuit of understanding complex combustion processes, molecular dynamics (MD) simulations provide an essential bridge between atomic-scale interactions and macroscopic observable phenomena. Within this domain, the Reactive Force Field (ReaxFF) has emerged as a powerful methodology that transcends the limitations of classical non-reactive force fields by enabling the simulation of bond formation and dissociation during chemical reactions. Unlike classical harmonic force fields such as CHARMM, AMBER, or PCFF, which maintain fixed bonding connectivity and are thus unsuitable for modeling chemical reactions, ReaxFF employs a bond-order formalism that dynamically describes covalent bonding interactions [10] [3]. This foundation allows ReaxFF to simulate complex reaction networks in hydrocarbon oxidation, pyrolysis, and detonation processes without predefined reaction pathways, making it particularly valuable for combustion chemistry research where reaction mechanisms are often too complex to enumerate completely [10].

The core strength of ReaxFF lies in its comprehensive energy function, which partitions the total system energy into multiple components including bond energy, lone-pair energy, over-coordination energy, valence angle strain, torsion energy, and non-bonded van der Waals and Coulomb interactions [10] [39]. This sophisticated framework, while more computationally expensive than classical force fields, provides the chemical fidelity necessary to model reactive processes with accuracy approaching quantum mechanical methods for large-scale systems containing thousands of atoms over meaningful timescales [39]. However, achieving accurate and reliable results with ReaxFF requires careful attention to three critical optimization parameters: time step integration, temperature control, and training data selection for force field parameterization. These choices collectively determine the numerical stability, thermodynamic accuracy, and predictive transferability of ReaxFF simulations in combustion applications.

Critical Optimization Parameters and Protocols

Time Step Selection and Numerical Integration

The choice of time step (Δt) in ReaxFF MD simulations represents a critical balance between numerical stability and computational efficiency. For modeling combustion reactions involving rapid bond breaking and formation, particularly with light hydrogen atoms exhibiting high-frequency vibrations, a time step of 0.1-0.25 femtoseconds (fs) is generally recommended [44]. This conservative time step ensures proper resolution of the fastest molecular vibrations and prevents spurious energy accumulation that can lead to simulation instability. The ReaxFF methodology employs specific algorithms to handle the dynamic changes in bonding connectivity that occur during reactions. The bond order calculation, updated every iteration based on interatomic distances, enables the force field to naturally describe bond dissociation and formation without predefined reaction templates [10]. For systems containing hydrogen bonds or involving high-temperature combustion conditions, the shorter end of this time step range (0.1 fs) is often necessary to maintain energy conservation and numerical accuracy throughout the simulation.

Table 1: Time Step Guidelines for ReaxFF Combustion Simulations

System Type	Recommended Time Step	Rationale	Typical Simulation Length
Hydrocarbon oxidation (with H atoms)	0.1 fs	Resolves high-frequency C-H/O-H vibrations	10-100 ps
Nanofluid fuel combustion (MNP/C₂H₅OH)	0.25 fs	Balance of stability and efficiency for metal-organic interfaces	50-200 ps
Energetic materials decomposition	0.1-0.25 fs	Captures rapid bond cleavage in NO₂ groups	10-50 ps
Cyclohexane combustion mechanism	0.25 fs	Suitable for ring-opening and oxidation pathways	50-100 ps

The computational expense of ReaxFF simulations, which can be approximately 30 times slower than the recently introduced IFF-R method that uses Morse potentials, necessitates careful consideration of the time step selection to maximize productive sampling of reactive events [3]. For complex combustion systems requiring longer simulation times, a hierarchical approach can be implemented where shorter time steps are used during phases of rapid reaction propagation, with slightly longer time steps employed during equilibrium or slow-combustion phases, though this requires careful validation to ensure artifact-free dynamics.

Temperature Control and Thermodynamic Ensembles

Proper temperature control is paramount in ReaxFF combustion simulations, as temperature directly governs reaction rates, ignition timing, and product distributions. The Nosé-Hoover thermostat has demonstrated excellent performance for combustion applications, providing improved thermal sampling compared to simpler thermostat algorithms [44]. For the Nosé-Hoover thermostat, a damping constant of 100 fs represents a reasonable default value, though this parameter should be adjusted based on system size to match the period of internal oscillations within the simulated system [44]. Combustion simulations typically employ elevated temperatures (2000-3500 K) to accelerate reaction kinetics within computationally feasible timescales, though this represents a compromise between observation of rare events and physiological relevance.

Table 2: Temperature Parameters in Representative ReaxFF Combustion Studies

Combustion System	Simulation Temperature	Thermostat Type	Key Findings
Methane/Oxygen mixture [44]	3500 K	Nosé-Hoover (100 fs damping)	Accelerated combustion pathway observation
Cyclohexane combustion [10]	2800 K	Not specified	Facilitated ring-opening and oxidation kinetics
Magnesium/Ethanol nanofluid [45]	Gradient: 1000-3500 K	Not specified	Ignition delay reduction with nanoparticle addition

The elevated temperatures commonly used in ReaxFF combustion studies (typically 2800-3500 K) serve to accelerate reaction kinetics beyond experimentally relevant conditions to observe complete combustion pathways within computationally feasible simulation timescales (picoseconds to nanoseconds) [44] [10]. This approach represents a recognized compromise where qualitative mechanistic insights are prioritized over quantitative rate predictions. For specific investigations of low-temperature combustion chemistry or ignition phenomena, temperature ramping protocols can be implemented, gradually increasing from ambient to high temperatures to observe spontaneous reaction initiation [45].

Training Data Selection and Force Field Parameterization

The accuracy and transferability of ReaxFF simulations are fundamentally determined by the quality and comprehensiveness of the training data used during parameter optimization. The parameterization process typically employs a multi-objective optimization approach that minimizes the discrepancy between ReaxFF predictions and high-level quantum mechanical (QM) reference data across diverse chemical configurations [39]. An effective training set must encompass bond dissociation curves for all relevant bond types, angle bending potentials, torsion profiles, reaction energies and barriers, and crystal structure data for condensed-phase systems [11]. For combustion chemistry specifically, the training must emphasize accurate representation of transition states and reaction enthalpies for key elementary reactions including hydrogen abstraction, radical recombination, and oxidation pathways.

Modern parameter optimization strategies combine global search algorithms like Simulated Annealing (SA) with local refinement methods such as Particle Swarm Optimization (PSO) to efficiently navigate the high-dimensional parameter space while avoiding local minima [39]. The introduction of a Concentrated Attention Mechanism (CAM) that weights chemically significant configurations (e.g., transition states, optimal structures) more heavily during optimization has been shown to further improve parameter accuracy [39]. It is crucial to recognize that specialized ReaxFF parameter sets exist for different chemical domains, primarily divided between the "combustion branch" and "aqueous branch," with the former optimized for gas-phase oxidation chemistry and the latter for condensed-phase aqueous environments [11]. Selecting the appropriate force field branch for the target application is essential for obtaining physically meaningful results.

Experimental Protocols and Validation Methodologies

Protocol: Methane Combustion Simulation

The following step-by-step protocol outlines a typical ReaxFF simulation for methane combustion, based on established computational methodologies [44]:

• System Construction: Create an initial simulation cell containing methane (CH₄) and oxygen (O₂) molecules at a stoichiometric ratio of 1:2.5 for complete combustion (CH₄ + 2O₂ → CO₂ + 2H₂O). Utilize the "Generate Molecules" function in molecular builder software to pack approximately 500 atoms at a density of 1.163 g/cm³ with random non-overlapping positions and orientations.

• Force Field Selection: Choose the appropriate ReaxFF parameter set for hydrocarbon oxidation (e.g., CHO.ff from the combustion branch) [11]. This force field has been specifically parameterized against DFT calculations for C/H/O bond dissociation energies, angle distortions, and reaction barriers relevant to combustion processes.

• Simulation Parameters: Configure the molecular dynamics parameters as follows:

  • Integration time step: 0.25 fs
  • Total simulation steps: 30,000 (7.5 ps)
  • Thermostat: Nosé-Hoover with damping constant of 100 fs
  • Temperature: 3500 K (elevated to accelerate reactions)
  • Ensemble: NVT (constant number, volume, and temperature)

• Simulation Execution and Monitoring: Run the simulation while monitoring progress through:

  • Real-time visualization of molecular structures
  • Temperature and conserved energy trajectories
  • Molecular population analysis tracking reactants (CH₄, O₂) and products (H₂O, CO₂)

• Trajectory Analysis: Upon completion, analyze the simulation trajectory to identify:

  • Reaction initiation time and ignition delay
  • Intermediate species formation (radicals, formaldehyde, carbon monoxide)
  • Final product distribution and combustion completeness
  • Dominant reaction pathways through sequence analysis

Validation and Analysis Techniques

Validating ReaxFF simulation outcomes against experimental and high-level computational reference data is essential for establishing credibility. Several analytical approaches facilitate this verification process:

• Chemical Trajectory Analysis: Tools like AMSmovie or custom scripts enable tracking of molecular populations throughout the simulation, providing kinetic profiles of reactant consumption and product formation [44]. For cyclohexane combustion, this approach has confirmed the production of expected intermediates (C₂H₄, CH₂O) and final products (CO, CO₂, H₂O) consistent with experimental observations [10].

• Reaction Pathway Mapping: By analyzing bond formation and cleavage sequences, researchers can reconstruct detailed reaction mechanisms. In magnesium/ethanol nanofuel simulations, this technique revealed that magnesium nanoparticles reduce ignition delay by 62.0% and enhance hydrogen production by 3.4 times through facilitated dehydrogenation pathways [45].

• Quantitative Benchmarking: Comparing computed properties such as activation energies, product distributions, and crystal parameters with experimental measurements or DFT calculations provides quantitative validation. For instance, the EMFF-2025 neural network potential demonstrated mean absolute errors within ±0.1 eV/atom for energies and ±2 eV/Å for forces when validated against DFT calculations for high-energy materials [13].

Table 3: Essential Resources for ReaxFF Combustion Simulations

Resource Category	Specific Examples	Function/Purpose
ReaxFF Force Fields	CHO.ff (hydrocarbon oxidation) [11], HE.ff (energetic materials) [11], ReaxFFCHO-S22 (cycloalkane combustion) [10]	Provide parameterized potentials for specific chemical systems and elements
Software Packages	LAMMPS [45], AMS [44], ReaxFF module in SCM [11]	MD engines with ReaxFF implementation for simulation execution
Analysis Tools	AMSmovie [44], Reaction Analysis tools [10], Python/MATLAB scripts	Visualization, trajectory analysis, and reaction pathway identification
Reference Data	DFT calculations [39], Experimental combustion data [10], Quantum mechanical benchmarks [13]	Training and validation targets for force field parameterization
Computational Resources	High-performance computing clusters, Parallel computing architectures [10]	Enable large-scale (1000+ atoms) simulations over relevant timescales

The effective application of ReaxFF molecular dynamics to combustion chemistry requires meticulous attention to three interdependent optimization domains: temporal discretization, thermal regulation, and parameter training. The selection of an appropriate time step (0.1-0.25 fs) ensures numerical stability while capturing the rapid bond rearrangements characteristic of combustion reactions. Implementation of proper temperature control through advanced thermostats like Nosé-Hoover maintains realistic thermodynamic conditions while potentially accelerating reactions to observable timescales. Most fundamentally, the choice of specialized force fields from the combustion branch, parameterized against comprehensive quantum mechanical training data specific to hydrocarbon oxidation chemistry, establishes the physical accuracy and predictive capability of the simulations.

When strategically implemented, ReaxFF provides a powerful computational platform that surpasses the limitations of classical non-reactive force fields for combustion research, enabling atomic-scale insights into complex reaction networks, ignition phenomena, and nanoscale additive effects. The continuing development of optimized parameter sets, enhanced sampling algorithms, and more efficient optimization frameworks promises to further expand the utility of ReaxFF for probing the fundamental chemical mechanisms underlying combustion processes across increasingly complex and technologically relevant fuel systems.

In computational chemistry, the choice of simulation method is fundamentally a trade-off between the electronic-level accuracy of quantum mechanics (QM) and the computational feasibility of simulating large systems over meaningful timescales. [2] For researchers studying complex processes in combustion chemistry, such as hydrocarbon oxidation or coal fire dynamics, this balance is critical. Quantum mechanical methods, while offering valuable theoretical guidance, are often too computationally intense for simulations that consider the full dynamic evolution of a system, particularly when dealing with the large molecular structures and extended reaction pathways typical in combustion research. [2] [38]

Classical force fields provide a computationally efficient alternative but typically require predefined connectivity between atoms, precluding simulations of reactive events where bonds break and form. [2] The Reactive Force Field (ReaxFF) methodology was developed specifically to bridge this gap, approaching it from the classical side by casting the empirical interatomic potential within a bond-order formalism. [2] This allows ReaxFF to describe chemical bonding implicitly without expensive QM calculations, enabling simulations of reactive chemistry at scales previously inaccessible to computational methods. [2]

Theoretical Foundations: Methodological Approaches Compared

The ReaxFF Formalism

ReaxFF employs a bond-order formalism in conjunction with polarisable charge descriptions to describe both reactive and non-reactive interactions between atoms. [2] The total system energy in ReaxFF is described by the equation:

Esystem = Ebond + Eover + Eangle + Etors + EvdWaals + ECoulomb + ESpecific [2]

Where Ebond represents bond energy, Eangle and Etors account for angle and torsion strain, Eover prevents atom over-coordination, and EvdWaals and ECoulomb describe non-bonded interactions. [2] The key innovation lies in the bond-order calculation, which is empirically determined from interatomic distances:

BOij = BOijσ + BOijπ + BOijππ = exp[pbo1(rijroσ)pbo2] + exp[pbo3(rijroπ)pbo4] + exp[pbo5(rijroππ)pbo6] [2]

This continuous, differentiable function contains no discontinuities through transitions between σ, π, and ππ bond character, enabling the simulation of bond breaking and formation. [2]

Classical Force Fields and Emerging Alternatives

Traditional classical force fields (e.g., CHARMM, AMBER, OPLS-AA) utilize harmonic bond potentials and static atomic charges, making them computationally efficient but incapable of modeling chemical reactions. [3] These force fields typically contain between 10-100 parameters with clear physical interpretations. [38]

More recently, machine learning potentials (MLPs) have emerged as powerful alternatives, with demonstrated capability to achieve density functional theory (DFT)-level accuracy at significantly reduced computational cost. [13] For instance, the EMFF-2025 neural network potential developed for energetic materials with C, H, N, and O elements can predict structures, mechanical properties, and decomposition characteristics with DFT-level accuracy while being more efficient than traditional quantum chemical calculations. [13]

The Reactive INTERFACE Force Field (IFF-R) represents another approach, replacing harmonic bond potentials with reactive Morse potentials while maintaining compatibility with existing force fields like CHARMM and PCFF. [3] This method reportedly offers about 30 times faster computational speed compared to prior reactive simulation methods like ReaxFF. [3]

Quantitative Comparison: Accuracy, Cost, and Applicability

Table 1: Comparison of Computational Methods for Molecular Dynamics Simulations

Method	Computational Cost	System Size Limit	Time Scale	Reactivity	Parameter Interpretability
Quantum Mechanics (DFT)	Very High	~100-1,000 atoms	Picoseconds to nanoseconds	Intrinsic	High (electronic structure)
ReaxFF	Moderate	>1,000,000 atoms	Nanoseconds to microseconds	Bond-order based	Moderate (physically motivated terms)
Classical FF	Low	>10,000,000 atoms	Nanoseconds to milliseconds	Non-reactive	High (clear physical meaning)
Machine Learning FF	Variable (training high, inference moderate)	~100,000 atoms	Nanoseconds to microseconds	Depends on training data	Low (black-box nature)
IFF-R	Low to Moderate	Similar to classical FF	Microseconds	Morse potential-based	High (interpretable parameters)

Table 2: Performance Metrics for ReaxFF in Combustion Applications

Application	System Size (atoms)	Simulation Time	Accuracy vs DFT	Key Validation
Hydrocarbon Oxidation [2]	Thousands	Nanoseconds	Close agreement	Reaction energies, barriers
Bituminous Coal Combustion [23]	50,789	Not specified	Validated for bond dissociation	Bond dissociation energies, valence angles
Lignite Combustion [46]	>1,000	Hundreds of picoseconds	Quantitative product analysis	Product distribution, reaction pathways
Mg/H₂O Reaction [47]	Not specified	Not specified	Trained against QM data	Reaction energies, equation of state

Practical Implementation: Experimental Protocols and Workflows

ReaxFF Molecular Dynamics Simulation Protocol

A typical ReaxFF MD simulation for combustion chemistry follows these methodological steps:

• System Initialization: Construct the molecular system using established molecular models. For coal combustion studies, this often involves well-characterized models like the Wiser bituminous coal model or the Illinois No. 6 coal model containing up to 50,789 atoms. [23]

• Force Field Selection: Choose appropriate ReaxFF parameters. For combustion applications, the "combustion branch" parameters (e.g., CHO.ff, HCONSB.ff) are typically employed, as they are optimized for gas-phase reactions as opposed to the "aqueous branch" used for solution chemistry. [11] [2]

• Energy Minimization: Perform initial energy minimization to relieve bad contacts and high-energy configurations in the starting structure.

• Equilibration Phase: Conduct equilibration MD in the NVT (constant Number, Volume, Temperature) or NPT (constant Number, Pressure, Temperature) ensemble to reach thermodynamic equilibrium. Typical combustion simulations employ temperatures from 1000-3000K to accelerate reactive events. [23] [46]

• Production Simulation: Run extended MD simulations in the NVE (constant Number, Volume, Energy) or NVT ensembles to collect trajectory data for analysis. Simulation times typically range from hundreds of picoseconds to nanoseconds, depending on system size and available computational resources. [46]

• Trajectory Analysis: Analyze the simulation trajectory using specialized tools to identify reaction products, intermediates, and pathways through bond order analysis, species tracking, and reactive event detection.

Parameterization and Validation

The development of ReaxFF parameters follows a rigorous protocol:

• Training Set Construction: A quantum mechanical training set is developed containing data on bond dissociations, angle distortions, reaction energies and barriers, and crystal properties relevant to the target application. [11]

• Parameter Optimization: Force field parameters are optimized to reproduce the QM training data using various algorithms. Recent advances include methods combining simulated annealing and particle swarm optimization with concentrated attention mechanisms (SA+PSO+CAM), which have demonstrated improved efficiency and accuracy compared to traditional approaches. [39]

• Validation: The parameterized force field is validated against additional QM and experimental data not included in the training set. For combustion applications, this typically includes validation of bond dissociation energies, reaction barriers, and product distributions. [23] [47]

Diagram 1: ReaxFF Methodology Workflow. This flowchart illustrates the comprehensive workflow for conducting ReaxFF studies, from system setup through parameter development, simulation execution, and final validation.

Case Studies in Combustion Chemistry

Coal Combustion Mechanisms

ReaxFF MD simulations have provided unprecedented atomic-level insights into coal combustion processes. In studies of bituminous coal combustion under oxygen-deficient conditions, ReaxFF simulations revealed that combustion proceeds through a coupling of oxidation and cracking reactions. [23] The simulations identified that ether bonds are the most active functional groups, with alkyl groups adjacent to ether bonds departing early in the simulation, followed by remaining oxygen atoms being attacked by free radicals to form carbon oxides or water. [23] In contrast, carbonyl groups were found to be the most stable functional groups. [23]

These simulations also quantified temperature and oxygen concentration effects, demonstrating that higher temperatures (2000-3000K) promote coal molecule cracking, while oxygen concentration strongly affects free radical production and consequently the formation pathways of carbon oxides and water. [23] The identification of ·OH, ·HO2, and C2O2 as important radicals in coal combustion provided mechanistic insights that would be difficult to obtain experimentally. [23]

Electric Field Effects on Combustion

Recent ReaxFF investigations have extended to exotic conditions such as combustion under external electric fields (EEF). Simulations of lignite combustion under EEF revealed that electric field intensity significantly influences combustion efficiency and product distribution. [46] At an equivalence ratio of 1, increasing EEF intensity from 0 to 1 V/Å reduced C-containing products by 41% while increasing semi-coke proportion threefold, indicating an inhibitory effect on the combustion process. [46]

The simulations provided mechanistic understanding of these effects, showing that EEF delays reaction times of functional groups - extending the time for –OH to generate H2O by 34.8 ps and delaying –COOH conversion to CO2 by 42 ps under 0.5 V/Å field. [46] This demonstrates ReaxFF's capability to capture subtle field effects on reaction kinetics that would be challenging to probe experimentally.

Table 3: Research Reagent Solutions for Combustion Chemistry Simulations

Resource	Type	Function	Application Examples
CHO.ff Force Field [11]	Parameter Set	Describes hydrocarbon oxidation	Combustion of hydrocarbons, soot formation
HCONSB.ff Force Field [11]	Parameter Set	Extended parameters for C/H/O/N/S/B systems	Coal char combustion, soot incandescence
JAX-ReaxFF Framework [47]	Optimization Tool	Gradient-based force field parameter optimization	Mg/O/H parameter development for nanoparticle combustion
SA+PSO+CAM Algorithm [39]	Optimization Method	Multi-objective parameter optimization	Automated parameter development for new element systems
AMS ReaxFF Module [37]	Simulation Software	Parallelized ReaxFF implementation with advanced analysis	Polymer combustion, surface chemistry, battery materials

Emerging Methodologies

The field of reactive molecular dynamics continues to evolve with several promising directions:

Machine Learning Force Fields: ML-based potentials like EMFF-2025 are achieving DFT-level accuracy for energetic materials while offering significant computational advantages. [13] These models can predict mechanical properties and decomposition behavior across a wide chemical space, though they currently face limitations in interpretability and require extensive training data. [13]

Hybrid Approaches: Methods like IFF-R that replace specific terms in traditional force fields with reactive potentials offer an attractive middle ground, maintaining the interpretability and efficiency of classical force fields while enabling bond dissociation. [3]

Advanced Optimization Algorithms: New parameter optimization approaches like the SA+PSO+CAM method are addressing the critical bottleneck of force field development, potentially expanding ReaxFF's applicability to novel material systems. [39]

For combustion chemistry researchers, the choice between ReaxFF and classical force fields depends fundamentally on the research question and scale requirements. ReaxFF provides the unique capability to simulate reactive events in systems of thousands to millions of atoms, offering atomic-level insights into complex combustion processes that bridge the gap between quantum accuracy and classical efficiency. The methodology has demonstrated particular value in studying coal combustion mechanisms, pollutant formation pathways, and complex reaction networks under various conditions.

As computational power increases and methodologies continue to advance, the balance between computational cost and accuracy will inevitably shift toward enabling more accurate simulations of larger systems for longer timescales. However, the fundamental tradeoff will remain, requiring researchers to strategically select computational methods aligned with their specific scientific objectives in combustion chemistry.

Benchmarking Performance Against Quantum Mechanics and Experiments

ReaxFF, standing for "reactive force field," represents a fundamental shift from traditional force fields in molecular dynamics simulations. Whereas traditional force fields are unable to model chemical reactions due to their dependence on pre-defined bonds, ReaxFF eschews explicit bonds in favor of bond orders, which allows for continuous bond formation and breaking during simulations [9]. This capability makes it particularly valuable for studying complex chemical processes in fields such as combustion chemistry, where reaction pathways are complex and dynamic. The foundation of ReaxFF's accuracy lies in its parameterization against high-fidelity quantum mechanical (QM) data, enabling it to achieve a balance between computational efficiency and quantum-level accuracy for large-scale systems that would be prohibitively expensive for pure QM methods [47] [14].

The development and validation of ReaxFF force fields constitute a critical methodology for researchers investigating reactive systems. This technical guide details the paradigms through which ReaxFF parameters are trained, optimized, and validated against quantum mechanical reference data, with specific consideration for applications in combustion research. The process ensures that the force field accurately describes both energetic relationships and dynamic reaction processes, providing researchers with a reliable tool for exploring chemical phenomena at the atomic scale.

Theoretical Foundation of ReaxFF

Bond Order Formalism and Energy Contributions

At the core of ReaxFF is a bond order-based formalism that dynamically describes covalent interactions between atoms. Unlike traditional molecular mechanics force fields that use fixed connectivity, ReaxFF calculates bond orders (BO) directly from interatomic distances, allowing them to evolve continuously during a simulation [47]. This bond order is subsequently used to calculate the energy associated with covalent bonds. The total system energy in ReaxFF is partitioned into several contributions, reflecting a comprehensive physical model:

E_total = E_bond + E_ov + E_under + E_lp + E_val + E_pen + E_tors + E_conj + E_vdWaals + E_Coulomb

These terms account for a wide range of atomic interactions, including valence angle strain (Eval), torsional angle forces (Etors), van der Waals interactions (EvdWaals), and Coulombic forces (ECoulomb) with dynamic, geometry-dependent charges [47] [14]. The sophisticated parameterization of these energy terms against QM data enables ReaxFF to accurately model complex chemical environments, including those found in combustion systems such as hydrocarbon oxidation and high-energy material decomposition [11] [14].

Charge Equilibration Methods

For describing electrostatic interactions, ReaxFF utilizes charge equilibration methods such as Electronegativity Equalization Method (EEM) or the more advanced Atom-Condensed Kohn-Sham DFT approximated to second order (ACKS2) [48]. These methods calculate dynamic, geometry-dependent atomic charges that respond to the chemical environment during simulations. The parameters for these charge equilibration schemes are critically derived from QM calculations to ensure accurate representation of molecular polarization and charge transfer in reactive processes, which are essential for modeling combustion intermediates and transition states [48].

Training ReaxFF Force Fields

Quantum Mechanical Training Set Construction

The development of a robust ReaxFF force field begins with the construction of an extensive and representative QM training set. This dataset must encompass the relevant chemical phase space that the force field is intended to describe [9]. For combustion chemistry applications, this typically includes systems containing C, H, O, and N elements [11]. The training set comprises various QM calculations designed to probe different aspects of the potential energy surface (PES):

• Bond dissociation curves for all possible bonds to characterize bond breaking and formation
• Angle and torsion bending scans to parameterize valence and dihedral interactions
• Reaction energy pathways for key chemical reactions, including transition states
• Crystal structure data and equations of state for condensed phases
• Interaction energies between molecules and surfaces for catalytic systems

For example, in developing a force field for hydrocarbon oxidation (CHO.ff), density functional theory (DFT) calculations were performed on various systems including "dissociation energies for various bonds containing carbon, oxygen, and hydrogen," with "dissociation curves calculated by constraining only the bond length of interest and re-optimization of the remaining internal coordinates" [11]. Similarly, for a high-energy materials force field (HE.ff), the training set contained "bond breaking and compression curves for all possible bonds, angle and torsion bending data for all possible cases, as well as crystal data" [11].

Parameter Optimization Process

The parameter optimization process involves systematically adjusting the hundreds of parameters in the ReaxFF energy expression until the force field's description of the training set properties matches the QM data as closely as possible. This is typically framed as a minimization problem where the objective is to reduce the difference between ReaxFF-predicted and QM reference values [47].

Modern optimization approaches utilize advanced algorithms to navigate the high-dimensional parameter space. For instance, the JAX-ReaxFF framework, based on a gradient descent algorithm, has been recently employed for efficient parameter optimization [47]. This framework allows for systematic optimization of interaction parameters against a QM training set, as demonstrated in the development of a Mg/O/H force field for studying magnesium nanoparticles in water systems [47].

Table: Representative Quantum Mechanical Data in ReaxFF Training Sets

Data Category	Specific Calculations	Purpose in Force Field Development	Example Systems
Bond Dissociation	Single, double, triple bond dissociation curves	Parameterize bond order relationships	H-H, C-C, C-O, O-O [11]
Angle/Torsion	Angle bending, torsion scans	Fit valence angle and dihedral parameters	H-O-H, H-C-H, C-O-H [11]
Reaction Pathways	Reaction energies and barriers	Capture chemical reactivity	CH₄ decomposition, H₂O dissociation [47] [49]
Crystal Properties	Equation of state, elastic constants	Reproduce condensed phase behavior	MgH₂, Mg(OH)₂ crystals [47]
Surface Interactions	Adsorption energies, migration barriers	Model heterogeneous catalysis	Hydrocarbons on Ni surfaces [49]

The parameter optimization is an iterative process, often requiring multiple rounds of refinement and validation. The complexity of this process is underscored by the fact that ReaxFF is "a fairly complex force field with many parameters," necessitating "an extensive training set covering the relevant chemical phase space" [9].

Validation Protocols for ReaxFF Force Fields

Energetic and Structural Validation

Once a ReaxFF force field has been trained, it must undergo rigorous validation against both the training data and additional QM calculations not included in the training set. Energetic validation involves comparing ReaxFF-predicted energies and forces with QM reference data across various chemical configurations. Effective validation demonstrates that the force field maintains chemical accuracy across a diverse range of structures and reactions.

A key metric for validation is the mean absolute error (MAE) for energies and forces. For instance, in the development of the EMFF-2025 neural network potential (a related but more advanced approach), the MAE for energy was predominantly within ±0.1 eV/atom, and the MAE for force was mainly within ±2 eV/Å [13]. Similar accuracy standards are applied to ReaxFF force fields during validation.

Specific validation tests include:

• Formation energies of molecules and solids compared to QM calculations
• Reaction energy profiles for key chemical processes
• Geometric parameters (bond lengths, angles) of stable structures
• Vibrational frequencies for small molecules

For example, in the development of a ReaxFF force field for Fe/Al/Ni alloys, the validation process confirmed that "ReaxFF accurately reproduces activation energies, demonstrating values of 18.4 kcal/mol, in good agreement with quantum results of 18.9 kcal/mol and experimental data of 17.7 kcal/mol" [49].

Reactive Dynamics Validation

Beyond static properties, ReaxFF must be validated for its performance in molecular dynamics simulations of reactive processes. This involves comparing reaction rates, mechanisms, and products against both QM data and experimental observations where available.

For combustion applications, common validation scenarios include:

• Pyrolysis pathways of hydrocarbons and energetic materials
• Oxidation mechanisms and intermediate species formation
• Catalytic processes on metal surfaces
• High-temperature decomposition behavior

In a study of hydrocarbon reactions on nickel surfaces, ReaxFF was validated by comparing its prediction of methane dissociation on different nickel surfaces against QM calculations. The force field correctly reproduced the observed reactivity trend: "the (111) surface is the least reactive, the (100) surface has the fastest reaction rates, and the stepped (111) surface has an intermediate reaction rate" [49]. This validation confirmed ReaxFF's ability to capture the importance of surface defects in accelerating reaction rates.

Table: Validation Metrics for ReaxFF in Combustion-Related Systems

Validation Category	Key Metrics	Acceptance Criteria	Example from Literature
Energetic Accuracy	Formation energies, reaction barriers	MAE < 0.1 eV/atom for energy	Mg/O/H force field validation [47]
Structural Fidelity	Bond lengths, angles, crystal parameters	Deviation < 1-2% from QM	CHO.ff for hydrocarbon oxidation [11]
Reaction Dynamics	Product distribution, reaction rates	Qualitative/quantitative match with experiment	PS hydrothermal gasification [50]
Surface Chemistry	Adsorption energies, reaction pathways	Activation energies within 1-2 kcal/mol of QM	Hydrocarbon dissociation on Ni [49]
Transferability	Performance on unseen molecules/systems	Maintains accuracy outside training set	HE.ff application to novel explosives [11]

Case Study: Combustion Chemistry Application

Hydrocarbon Oxidation with CHO Force Field

The development and validation of the CHO.ff force field for hydrocarbon oxidation exemplifies the ReaxFF training paradigm. This force field was specifically parameterized for combustion applications, with training data derived from DFT calculations on relevant systems [11]. The training process involved calculating "dissociation energies for various bonds containing carbon, oxygen, and hydrogen," with "dissociation curves calculated by constraining only the bond length of interest and re-optimization of the remaining internal coordinates" [11].

Validation of the CHO force field included its application to methane combustion simulations, where it successfully reproduced expected reaction pathways and product distributions. The force field has since become a foundation for many combustion-related ReaxFF studies and has been extended to include additional elements important in combustion systems, such as nitrogen and sulfur in the HCONSB.ff force field [11].

Magnesium Combustion in Water Atmosphere

A recent example demonstrates the development of a Mg/O/H ReaxFF force field to study the combustion of magnesium nanoparticles in water atmosphere for hydrogen production [47]. The force field parameters were optimized using the JAX-ReaxFF framework based on a QM training set that included "the interactions between Mg/O/H, as well as the equation of state of MgH₂ and Mg(OH)₂ crystals" [47].

The validation process confirmed that the optimized force field accurately reproduced the QM-derived potential energy surfaces for various molecular interactions. Subsequently, molecular dynamics simulations revealed the detailed mechanism of magnesium-water combustion: "H₂O dissociates on the surface of Mg nanoparticles to form Mg-H, Mg-OH, and Mg-O bonds," and "structural evolutions of Mg nanoparticles depend on the temperature and the density of H₂O" [47]. This case study illustrates how a properly validated ReaxFF force field can provide atomic-level insights into combustion processes that are difficult to observe experimentally.

Research Reagent Solutions: Essential Tools for ReaxFF Studies

Table: Key Computational Tools for ReaxFF Development and Application

Tool/Solution	Function	Application in ReaxFF Workflow
Quantum Chemistry Software	Generates reference training data	Calculation of bond dissociation energies, reaction pathways, and property surfaces [47]
ReaxFF Parameterization Tools	Optimizes force field parameters	Fitting ReaxFF parameters to QM data using algorithms like JAX-ReaxFF [47]
Molecular Dynamics Engines	Performs ReaxFF simulations	Production MD simulations of reactive systems [44] [50]
Visualization Software	Analyzes simulation trajectories	Identification of chemical species and reaction mechanisms [44]
Validation Scripts	Compares ReaxFF and QM results	Quantitative assessment of force field accuracy [48]

Workflow Diagram: ReaxFF Development and Validation

The following diagram illustrates the comprehensive workflow for developing and validating a ReaxFF force field, from initial QM data generation to final application in molecular dynamics simulations:

ReaxFF Force Field Development and Validation Workflow

The training and validation of ReaxFF force fields against quantum mechanical data represent a rigorous methodology for developing accurate computational tools for combustion chemistry research. The process involves careful construction of QM training sets that encompass the relevant chemical space, systematic parameter optimization using advanced algorithms, and comprehensive validation against both QM data and experimental observations. When properly developed and validated, ReaxFF provides a powerful platform for investigating complex reactive processes in combustion systems at the atomic level, offering insights that bridge the gap between quantum accuracy and molecular dynamics scale. As computational resources advance and training methodologies refine, ReaxFF continues to evolve as an indispensable tool in the computational chemist's toolkit for understanding and predicting combustion phenomena.

The accurate prediction of reaction barriers and kinetics is a cornerstone of computational chemistry, directly impacting the design of novel materials and therapeutic compounds. In the specialized field of combustion chemistry, where reactions involve complex networks of bond breaking and formation, the choice of computational method is critical. While classical force fields offer computational efficiency for non-reactive molecular dynamics, their fixed bonding representations are fundamentally unsuited for modeling chemical reactions [14]. In contrast, reactive force fields (ReaxFF) employ a bond-order formalism to enable dynamic bond breaking and formation, bridging the gap between quantum mechanical accuracy and the scale of molecular dynamics simulations [2]. This technical analysis provides a direct performance comparison of these methodologies, evaluating their accuracy in predicting reaction barriers and kinetics, with specific applications in combustion chemistry.

Fundamental Methodological Differences

The core distinction between classical and reactive force fields lies in their treatment of chemical bonds, which directly dictates their capability to model reaction pathways.

Classical Force Fields

Classical force fields calculate system energy using simplified interatomic potential functions, typically comprising:

• Harmonic potentials for bond stretching and angle bending
• Lennard-Jones potentials for dispersion forces
• Coulombic interactions for electrostatic forces [38]

These force fields contain 10-100 parameters with high interpretability, as each term corresponds to a clear physical quantity. However, they maintain predefined connectivity between atoms, precluding simulations of bond breaking and formation. This limitation restricts their application to non-reactive interactions, such as adsorption, diffusion, and conformational changes, where chemical identity remains constant [38] [14].

Reactive Force Fields (ReaxFF)

ReaxFF introduces a bond-order formalism derived from interatomic distances, creating a continuous potential energy surface that allows bonds to break and form dynamically. The total system energy in ReaxFF includes multiple contributions:

[E{\text{system}} = E{\text{bond}} + E{\text{over}} + E{\text{angle}} + E{\text{tors}} + E{\text{vdWaals}} + E{\text{Coulomb}} + E{\text{Specific}}]

where the bond order ( BO_{ij} ) is empirically calculated as:

[BO{ij} = BO{ij}^{\sigma} + BO{ij}^{\pi} + BO{ij}^{\pi\pi} = \exp\left[p{bo1}\left(\frac{r{ij}}{ro^{\sigma}}\right)^{p{bo2}}\right] + \exp\left[p{bo3}\left(\frac{r{ij}}{ro^{\pi}}\right)^{p{bo4}}\right] + \exp\left[p{bo5}\left(\frac{r{ij}}{ro^{\pi\pi}}\right)^{p{bo6}}\right]] [2]

This approach requires 100-1,000 parameters and incorporates polarizable charge descriptions to model both covalent and electrostatic interactions for diverse materials [38]. The method's parameterization is trained against quantum mechanical data and optionally experimental data, enabling it to reproduce reaction energies and barriers with reasonable fidelity to quantum mechanics while accessing larger scales [47] [2].

Emerging Alternatives

Recent advancements have introduced hybrid and machine-learning approaches:

• IFF-R: Replaces harmonic bond potentials with Morse potentials to enable bond dissociation while maintaining compatibility with classical force fields, offering ~30x speedup over ReaxFF [3].
• Machine Learning Force Fields (MLFFs): Use active learning frameworks to autonomously construct training sets from quantum mechanical data, achieving quantum accuracy with significantly improved computational efficiency [51] [52].

Table 1: Fundamental Characteristics of Force Field Types

Feature	Classical Force Fields	Reactive Force Fields (ReaxFF)	Machine Learning Force Fields
Bond Treatment	Fixed connectivity	Dynamic via bond order	Learned from QM data
Parameter Count	10-100	100-1,000	1,000-100,000+
Reactive Capability	No	Yes	Yes
Computational Cost	Low	Moderate to High	Variable (often High for training)
Primary Applications	Conformational changes, diffusion	Chemical reactions, combustion, catalysis	Complex PES, specific reactive systems

Quantitative Performance Comparison

Accuracy in Reaction Barrier Prediction

The accuracy of reaction barrier prediction is crucial for modeling chemical kinetics. ReaxFF has demonstrated capability in reproducing quantum mechanical reaction barriers for various systems, though with variable precision:

In combustion applications, ReaxFF simulations of hydrogen-oxygen systems have successfully identified key chain reactions such as:

• H₂ + O₂ → HO₂ + H
• H + O₂ → HO₂
• HO₂ + H → H₂O₂

These simulations revealed inhibition effects of CF₃CHCl₂ (HCFC-123), showing how F and Cl radicals compete with chain propagation by consuming H and OH radicals through reactions like H₂ + Cl → HCl + H and H + F → HF [53]. The calculated activation energy for the one-step H₂-O₂ oxidation reaction increased from 33.77 kcal/mol to 37.68 kcal/mol with CF₃CHCl₂ addition, demonstrating ReaxFF's sensitivity to chemical inhibitors [53].

For magnesium-water systems, a recently developed Mg/O/H ReaxFF force field accurately simulated hydrogen production mechanisms, revealing that H₂ production trails H₂O consumption and occurs through inward diffusion of O atoms and outward diffusion of Mg atoms from magnesium hydride [47].

Machine learning force fields have demonstrated exceptional accuracy in catalytic systems. For CO₂ hydrogenation on indium oxide, MLFFs achieved energy barriers within 0.05 eV of Density Functional Theory (DFT) calculations, enabling discovery of an alternative reaction path with a 40% reduction in activation energy [51].

Computational Efficiency and System Size

Classical force fields maintain significant advantages in computational efficiency, enabling microsecond simulations of systems containing 100,000+ atoms [38]. Their simplified functional forms allow rapid calculation of forces, making them ideal for studying large-scale molecular motion without reactive events.

ReaxFF simulations are ~30 times slower than classical MD but substantially faster than quantum mechanical methods [14] [3]. Typical ReaxFF applications handle systems of 1,000-10,000 atoms for time scales up to nanoseconds, sufficient to observe complex reaction phenomena in heterogeneous systems [2].

The recently developed IFF-R method addresses this efficiency gap, offering ~30x speedup over ReaxFF while maintaining reactive capabilities through Morse potentials instead of bond-order formalism [3].

Table 2: Direct Performance Comparison for Combustion Chemistry

Performance Metric	Classical Force Fields	ReaxFF	Reference
Reaction Barrier Accuracy	Not applicable	± few kcal/mol for trained systems	[53]
Activation Energy Sensitivity	Not applicable	Detects ~4 kcal/mol changes from inhibitors	[53]
Typical System Size	100,000+ atoms	1,000-10,000 atoms	[38] [2]
Simulation Timescale	Microseconds	Nanoseconds	[38] [2]
H₂-O₂ Combustion Mechanism	Cannot simulate	Identifies key radicals (HO₂, H₂O₂) and inhibition pathways	[53]

Experimental Protocols for Validation

ReaxFF Training and Validation Workflow

The development of a reliable ReaxFF parameter set follows a rigorous protocol:

Training Set Construction: The process begins with creating a quantum mechanical (QM) training set containing:

• Bond dissociation curves for all relevant bonds
• Angle and torsion bending data
• Reaction energies and barriers for key reaction steps
• Crystal data for condensed phases [11]

For example, the development of the Ammonia Borane (AB.ff) force field included QM data describing single, double, and triple bond dissociations for all B/N/O/H combinations, angular distortions in small AB-related molecules, and reaction barriers for key steps like H₂ release from AB [11].

Parameter Optimization: Using the JAX-ReaxFF framework with gradient descent algorithms, parameters are optimized to minimize the difference between ReaxFF predictions and QM training data [47]. This includes optimizing interaction parameters between element pairs (e.g., Mg-H, Mg-O, H-O for magnesium-water systems) and validating against equations of state for relevant crystals (e.g., MgH₂ and Mg(OH)₂) [47].

Validation: The force field is validated against QM and experimental data not included in the training set, ensuring transferability across different molecular environments and phases [2].

Reactive Molecular Dynamics Simulation Protocol

Once validated, ReaxFF molecular dynamics simulations follow this standardized protocol:

System Construction: Build the initial configuration containing reactants in appropriate stoichiometry. For combustion studies, this typically involves fuel and oxidizer molecules in a periodic simulation box [53].

Equilibration: Perform energy minimization followed by thermal equilibration at the target temperature using thermostats (e.g., Berendsen or Nosé-Hoover).

Production Simulation: Conduct the reactive MD simulation while monitoring:

• Temperature and pressure evolution
• Chemical species concentrations over time
• Bond formation and breaking events
• Radical formation and consumption

Analysis: Identify reaction products, calculate ignition delay times, and determine activation energies through Arrhenius analysis [53]. For example, in H₂-O₂ combustion studies with inhibitors, analysis includes tracking the formation of HO₂ and H₂O₂ radicals and measuring the system's ignition delay time [53].

Application in Combustion Chemistry

Hydrogen-Oxygen Combustion Mechanisms

ReaxFF has provided atomic-level insights into H₂-O₂ combustion, particularly above the secondary explosion limit where HO₂ and H₂O₂ chemistry dominates. Simulations have revealed:

• The initial chain reactions (H₂ + O₂ → HO₂ + H, H + O₂ → HO₂, and HO₂ + H → H₂O₂) are critical in the early combustion stages
• At high pressures, OH generation occurs mainly through H₂O₂ → OH + OH and H + HO₂ → OH + OH [53]
• The ignition delay time increases from 150 ps to 550 ps with 10% CF₃CHCl₂ inhibitor addition
• The activation energy for one-step oxidation increases from 33.77 kcal/mol to 37.68 kcal/mol with inhibitor addition [53]

Inhibition Mechanisms

ReaxFF-MD simulations have elucidated the inhibition mechanism of CF₃CHCl₂ (HCFC-123) in H₂-O₂ combustion through:

• Physical effects: Dilution and endothermic decomposition reducing system reactivity
• Chemical interference: F and Cl radicals competing for H and OH radicals through reactions:
  • H₂ + Cl → HCl + H
  • H + F → HF
  • Cl + OH → HCl + O [53]

These reactions reduce the concentration of free radicals (H, OH, O), thereby inhibiting the chain-branching reactions essential for flame propagation.

Magnesium-Water Reaction for Hydrogen Production

In energy applications, ReaxFF simulations have revealed the structural evolution of magnesium nanoparticles in water atmospheres:

• H₂O dissociation on Mg surfaces forms Mg-H, Mg-OH, and Mg-O bonds
• The inward diffusion rate of H atoms surpasses that of O atoms, forming a Mg hydride core and Mg oxide shell
• H₂ production correlates with H₂O consumption but trails behind it
• The inward diffusion of O atoms and outward diffusion of Mg atoms facilitate H₂ release from magnesium hydride [47]

Table 3: Essential Resources for Reactive Molecular Dynamics Studies

Resource Category	Specific Examples	Function/Application	Relevance to Combustion Chemistry
ReaxFF Force Fields	CHO.ff (Hydrocarbon oxidation), HCONSB.ff (Soot formation), HE.ff (High explosives)	Element-specific parameter sets for different chemical environments	CHO.ff for fuel oxidation; HCONSB.ff for soot formation pathways [11]
ReaxFF Software Implementations	LAMMPS, ADF, PuReMD, Materials Studio	Molecular dynamics engines with ReaxFF integration	LAMMPS widely used for combustion simulations [2]
Quantum Chemistry Codes	VASP, CP2K, Gaussian, Q-Chem	Generate training data for force field parameterization	Provide reference data for reaction barriers and energies [38]
System Preparation Tools	AMSinput, REACTER toolkit	Model building and simulation setup	AMSinput enables force field selection and system construction [11]
Analysis Tools	TRAVIS, VMD, OVITO	Trajectory analysis and visualization	Identify reaction products and intermediates from MD trajectories [53]

The direct performance comparison between classical and reactive force fields reveals a fundamental trade-off between computational efficiency and chemical accuracy. Classical force fields remain indispensable for large-scale non-reactive simulations but cannot describe bond breaking/formation, making them unsuitable for predicting reaction barriers and kinetics. ReaxFF provides a practical compromise, enabling nanosecond-scale simulations of reactive processes with reasonable accuracy for trained systems, as demonstrated in combustion studies of H₂-O₂ systems and Mg-water reactions. Emerging methodologies like IFF-R and machine learning force fields offer promising directions, with MLFFs achieving near-DFT accuracy for specific catalytic systems. For combustion chemistry research, ReaxFF currently represents the most balanced approach, providing atomistic insights into complex reaction networks while accessing time and length scales relevant to experimental phenomena.

In the computational modeling of molecular systems, the choice of force field dictates the boundaries of what phenomena can be simulated. For combustion chemistry research—a field defined by complex reaction networks, bond cleavage, and species transformation—this choice is critical. Classical molecular mechanics (MM) force fields and the reactive force field (ReaxFF) represent two philosophically distinct approaches to calculating system energy and atomic interactions. Classical force fields, such as those defined in CHARMM, AMBER, and OPLS-AA, utilize a harmonic potential framework that excels at simulating systems near equilibrium but is fundamentally incapable of modeling chemical reactions due to its predefined, static connectivity [54] [3]. In stark contrast, ReaxFF employs a bond-order formalism that dynamically describes covalent interactions, thereby enabling the simulation of bond formation and breaking during molecular dynamics (MD) trajectories [2] [8]. This review provides a technical dissection of these methodologies, evaluating their respective strengths and limitations within the specific context of combustion research, where the inability to model reactivity renders classical force fields of limited utility for probing the core chemical phenomena.

Fundamental Principles and Functional Forms

The Harmonic Foundation of Classical Force Fields

Classical force fields are built upon a Class I additive potential energy function, which partitions the total system energy into bonded and nonbonded terms. The mathematical expression is given by:

[ E{system} = E{bonded} + E_{nonbonded} ]

Where the components are explicitly defined as follows [54]:

• Bonded Terms: ( E{bonded} = \sum{bonds} Kb(b-b0)^2 + \sum{angles} K\theta(\theta-\theta0)^2 + \sum{dihedrals} \sum{n} K{\phi,n}(1+\cos(n\phi-\deltan)) + \sum{improper} K{\varphi}(\varphi-\varphi0)^2 )
• Nonbonded Terms: ( E{nonbonded} = \sum{ij} \frac{qi qj}{4\pi D r{ij}} + \sum{ij} \varepsilon{ij} \left[ \left( \frac{R{min,ij}}{r{ij}} \right)^{12} - 2 \left( \frac{R{min,ij}}{r_{ij}} \right)^6 \right] )

The bonded terms rely on harmonic potentials for bonds and angles, with fixed equilibrium values ((b0), (\theta0)) and force constants ((Kb), (K\theta)). This formulation is computationally efficient but requires predefinition of atomic connectivity, creating an insurmountable barrier to simulating chemical reactions where bonds break and form [54] [3]. The nonbonded terms use fixed partial charges ((qi), (qj)) and Lennard-Jones potentials, which lack the polarization effects crucial for accurately modeling charge transfer during reactions in combustion environments [55].

The Bond-Order Formalism of Reactive Force Fields

ReaxFF bridges the gap between quantum mechanical accuracy and classical efficiency by introducing a bond-order formalism that dynamically describes covalent interactions. The overall system energy in ReaxFF contains significantly more terms than classical force fields [8]:

[ E{system} = E{bond} + E{lp} + E{over} + E{under} + E{val} + E{pen} + E{coa} + E{C2} + E{tors} + E{conj} + E{Hbond} + E{vdWaals} + E{Coulomb} ]

The most critical innovation is that bond order ((BO{ij})), which determines the strength of the connection between atoms (i) and (j), is calculated directly from the interatomic distance (r{ij}) [2]:

[ BO{ij} = BO{ij}^\sigma + BO{ij}^\pi + BO{ij}^{\pi\pi} = \exp\left[p{bo1}\left(\frac{r{ij}}{r0^\sigma}\right)^{p{bo2}}\right] + \exp\left[p{bo3}\left(\frac{r{ij}}{r0^\pi}\right)^{p{bo4}}\right] + \exp\left[p{bo5}\left(\frac{r{ij}}{r0^{\pi\pi}}\right)^{p{bo6}}\right] ]

This distance-dependent bond order allows bonds to form, break, and transition between single, double, and triple character seamlessly during simulations. Additionally, ReaxFF employs a charge equilibration method that dynamically calculates atomic partial charges at each time step, capturing charge transfer effects that are critical in combustion reactions [2] [8].

Comparative Visualization of Methodologies

The diagram below illustrates the fundamental differences in how classical and reactive force fields calculate system energy and determine atomic interactions.

Quantitative Comparison of Force Field Capabilities

Performance Metrics for Combustion Chemistry

Table 1: Technical comparison of classical and reactive force fields for combustion-relevant properties

Property	Classical Force Fields	ReaxFF	Implications for Combustion Research
Bond Dissociation	Not possible; fixed connectivity	Continuous description via bond order [2]	Enables simulation of initiation, propagation, termination steps
Charge Transfer	Fixed partial charges [54]	Dynamic at each step via QEq/EEM [8]	Critical for radical chemistry and oxidation reactions
Parameterization	Multiple atom types per element [3]	Single atom type per element [8]	ReaxFF more transferable but requires element-specific parameterization
Reaction Barriers	Not applicable	Approximated through training [2]	Reasonable kinetics without QM computational cost
Computational Cost	1x (reference) [3]	~30x classical FFs [3]	Limits system size and simulation timescale
Combustion Branch	Not applicable	Specifically parameterized [56]	Accurate for gas-phase hydrocarbon oxidation [56]

Specialized ReaxFF Parameterizations for Combustion

The ReaxFF framework has been specifically adapted for combustion chemistry through dedicated parameterizations. The combustion branch of ReaxFF, distinct from the aqueous branch, focuses on accurately describing gas-phase hydrocarbon oxidation processes [56]. Key parameterizations include:

• CHO.ff: Specifically designed for hydrocarbon oxidation, parameterized against dissociation energies, angular distortions, and torsional profiles for C/H/O systems [56].
• ReaxFFCHO-S22: A recently developed force field re-parameterized against high-level quantum mechanical calculations that successfully simulated cyclohexane combustion mechanisms, including ring-opening via C-C bond cleavage and subsequent oxidation pathways [10].
• HCONSB.ff: Extends capabilities to nitrogen and sulfur-containing compounds, enabling simulation of fuel impurities and nitrogen oxide formation [56].

These specialized parameterizations allow ReaxFF to capture complex combustion phenomena such as polycyclic aromatic hydrocarbon (PAH) formation, soot precursor chemistry, and pollutant formation mechanisms that are completely inaccessible to classical force fields [14] [10].

Methodologies for Combustion Chemistry Simulation

Protocol for ReaxFF Combustion Simulations

Implementing ReaxFF molecular dynamics for combustion research requires specific methodological considerations:

• System Preparation:

  • Construct simulation cell with fuel and oxidizer molecules at appropriate equivalence ratio
  • For cyclohexane combustion studies, a typical setup includes 12 cyclohexane and 108 O₂ molecules in a periodic box (33.34 Å side length) at density of 0.2 g/cm³ [10]

• Force Field Selection:

  • Select appropriate combustion-branch parameterization (e.g., CHO.ff, ReaxFFCHO-S22)
  • Verify element compatibility (C/H/O for basic hydrocarbons; N/S for advanced fuels)

• Simulation Parameters:

  • Use smaller timesteps (0.25 fs) compared to classical MD (1-2 fs) due to bond order calculations [8]
  • Set temperature typically between 2500-3000K to accelerate reaction kinetics within feasible simulation time [10]
  • Implement NVT or NVE ensembles depending on research question

• Simulation Execution:

  • Run for sufficient duration (typically 1 ns for combustion studies) to observe complete reaction pathways [10]
  • Monitor bond orders, species formation, and reaction intermediates

• Trajectory Analysis:

  • Use specialized tools (e.g., R.A.V.E.) to identify chemical species and reactions from bond order data [10]
  • Track evolution of key intermediates (C₂H₄, CH₂O, CO, CO₂, H₂O) and radical species

Table 2: Essential computational tools for reactive molecular dynamics in combustion research

Tool/Resource	Function	Availability
LAMMPS	Open-source MD engine with ReaxFF implementation [2]	Publicly available
PuReMD	Purdue Reactive MD, optimized for ReaxFF simulations [2]	Publicly available
ReaxFF Parameter Sets	Combustion-specific parameterizations (CHO.ff, S22) [56] [10]	SCM, literature
R.A.V.E.	Reaction Analysis and Visualization tool [10]	Various implementations
AMSinput	GUI for setting up ReaxFF calculations [56]	Commercial (SCM)
CHARMM, AMBER	Classical force field software for non-reactive properties [54]	Academic/commercial

Emerging Methodologies and Future Directions

While ReaxFF has significantly advanced reactive simulation capabilities, several emerging methodologies show promise for addressing its limitations:

• Neural Network Potentials (NNPs): Models like EMFF-2025 offer DFT-level accuracy for C/H/N/O systems with improved computational efficiency, showing particular promise for energetic materials and complex decomposition pathways [13].
• Hybrid ReaxFF/MM Approaches: The ReaxFF/AMBER framework enables reactive simulations of localized chemical events in large biomolecular systems, potentially applicable to biofuel combustion chemistry [8].
• Simplified Reactive Methods: The IFF-R approach replaces harmonic bonds with Morse potentials, offering a 30-fold speed advantage over ReaxFF while maintaining reasonable accuracy for bond dissociation [3].

These methodologies represent the ongoing effort to balance computational cost with chemical accuracy in reactive simulations, potentially offering new approaches for multiscale combustion modeling that bridges from electronic structure to macroscopic phenomena.

For combustion chemistry research, the verdict on classical force fields is unequivocal: their fundamental architectural limitations render them unsuitable for investigating chemical reactions central to combustion processes. The harmonic potential framework, with its fixed connectivity and static charge distributions, cannot describe the bond breaking, formation, and complex rearrangement events that define fuel oxidation and pollutant formation pathways. ReaxFF, with its bond-order formalism, dynamic charge equilibration, and specialized combustion parameterizations, provides the necessary theoretical framework to simulate these complex reactive phenomena. While the computational cost of ReaxFF remains significantly higher than classical force fields, and parameterization for new elements requires substantial effort, its capacity to simulate chemical reactions without predefined pathways makes it an indispensable tool for combustion researchers. The ongoing development of neural network potentials and hybrid methodologies promises to further expand the accessible timescales and system sizes for reactive simulation, potentially enabling fully predictive computational combustion science from molecular to device scales.

Combustion science, fundamental to energy conversion, aerospace propulsion, and environmental engineering, requires a deep understanding of complex chemical processes spanning multiple temporal and spatial scales. The accurate simulation of combustion dynamics has been a longstanding challenge due to the intricate interplay of chemical reactions, heat transfer, mass transfer, and fluid flow. Computational methods, particularly molecular dynamics (MD) simulations, have emerged as crucial tools for probing these processes at the atomic level. However, traditional approaches to modeling the potential energy surface (PES)—which determines how atoms interact and react—have faced a persistent trade-off between computational efficiency and physical accuracy. Classical force fields, while computationally efficient, cannot simulate bond breaking and formation essential for combustion chemistry. Reactive force fields like ReaxFF address this limitation but require complex parameterizations and remain computationally demanding. The emerging paradigm of machine learning force fields (MLFFs) is now transforming this landscape by offering a path to quantum-mechanical accuracy with near-classical computational efficiency, positioning them as powerful tools for unraveling combustion mechanisms and accelerating the development of cleaner combustion technologies. [57] [38]

This whitepaper provides an in-depth technical examination of MLFFs within the specific context of combustion science. It contrasts this emerging approach with established reactive methodologies, details the construction and application of specialized MLFFs, and presents a practical research framework for their implementation, supported by quantitative comparisons and experimental protocols.

The Force Field Landscape in Combustion Research

The computational study of combustion chemistry relies on various force fields, each with distinct capabilities and limitations. The following table summarizes the core characteristics of the three main classes.

Table 1: Comparison of Force Field Types for Combustion Applications

Feature	Classical Force Fields	Reactive Force Fields (ReaxFF)	Machine Learning Force Fields (MLFFs)
Reactivity	Non-reactive (fixed bonding)	Reactive (dynamic bonding via bond-order)	Reactive (inherently reactive PES)
Functional Form	Simple, physics-based potentials (e.g., harmonic)	Complex, with many energy terms (>20) [3]	Data-driven, complex non-linear functions (e.g., neural networks)
Accuracy	Low for reactions; good for structures	Medium; system-dependent [13]	High; can achieve DFT-level accuracy [13]
Computational Cost	Lowest	High (~30x classical) [3]	Variable (can be lower than ReaxFF) [13]
Parameterization	10-100 interpretable parameters [38]	100+ parameters, complex fitting [38]	"Black-box" models; require large QM datasets
Development Focus	Reproducing equilibrium properties	Capturing reaction pathways	Interpolating and generalizing from QM data
Key Combustion Limitation	Cannot simulate chemical reactions	Transferability; accuracy of reaction barriers [13]	Data hunger; extrapolation reliability

Classical force fields use simplified analytical functions to describe bonded (e.g., bonds, angles) and non-bonded (e.g., van der Waals, electrostatic) interactions. While computationally efficient, their fixed bonding topology makes them unsuitable for simulating the bond-breaking and formation events at the heart of combustion. [38]

Reactive force fields, most notably ReaxFF, overcome this by employing a bond-order formalism, where the energy terms depend on the instantaneous distances between atoms, allowing bonds to form and break dynamically. While this enables MD simulations of complex reactions without predefined pathways, it comes at a significant cost. ReaxFF's functional form contains many more energy terms than classical force fields, making parameterization complex and computationally expensive, with simulations running approximately 30 times slower than non-reactive classical simulations. [3] [10] Furthermore, its accuracy in describing reaction potential energy surfaces, especially for new systems, can struggle to match that of quantum mechanics. [13]

MLFFs represent a paradigm shift. They bypass predefined physical models and instead use machine learning to directly learn the relationship between a system's atomic configuration and its potential energy and forces from high-quality quantum mechanical (QM) data. This approach allows them to potentially achieve density functional theory (DFT)-level accuracy while being orders of magnitude faster, making high-fidelity, reactive molecular dynamics simulations of combustion systems more accessible than ever before. [13] [38]

Machine Learning Force Fields: Methodology and Workflow

The development of a robust MLFF follows a structured pipeline designed to capture the underlying physics of combustion systems accurately. The core components and workflow are detailed below.

Core Architectures and Physical Constraints

MLFFs rely on an atomic descriptor that converts the positions and chemical species of an atom and its neighbors into a fixed-length, rotationally and translationally invariant vector. This descriptor serves as the input for a machine learning model, typically a neural network, which predicts the atomic energy. The total energy of the system is the sum of all atomic energies, and forces are obtained via automatic differentiation. [58]

Key architectures include:

• Deep Potential (DP): A widely used framework that employs a local descriptor and a deep neural network to represent the atomic potential, demonstrating exceptional capabilities in modeling complex reactive processes. [13]
• Message Passing Neural Networks (MPNNs): Graph-based models that iteratively pass information between bonded (or nearby) atoms, effectively capturing complex many-body interactions. [58]
• Atomic Cluster Expansion (ACE): Provides a complete basis in the space of atomic environments, forming a rigorous foundation for constructing accurate and transferable potentials. [58]

A critical advancement is the integration of physical priors through Physics-Informed Machine Learning (PIML). Physical knowledge, such as underlying symmetries and conservation laws, can be embedded into the ML model through its architecture (e.g., built-in equivariance) or as soft constraints in the loss function. This enhances the model's physical consistency and reliability, especially in extrapolative scenarios common in combustion. [57]

The MLFF Development Pipeline

The process of creating a specialized MLFF for combustion research involves several key stages, visualized in the workflow below.

Title: MLFF Development and Deployment Workflow

• Data Generation: The foundation of any MLFF is a high-quality, diverse dataset of atomic configurations with corresponding energies and forces computed from a high-level QM method (e.g., DFT, CCSD(T)). For combustion, this set must include reactants, products, transition states, intermediates, and various distorted molecular geometries sampled from QM-based molecular dynamics (QM-MD) or path-based methods like Nudged Elastic Band (NEB). [13]

• Model Training: A machine learning model (e.g., a neural network) is trained to minimize the difference between its predictions and the QM reference data. The loss function (L) is typically a weighted sum of energy and force errors: L = w_E|E_pred - E_QM|² + w_FΣ|F_pred - F_QM|² Ensuring energy conservation is a critical aspect, sometimes requiring specific parameterizations or distillation techniques from non-conservative teachers. [59]

• Active Learning: This iterative process is crucial for robustness. The initial MLFF is used to run MD simulations. Configurations where the model is uncertain (e.g., with high prediction variance) are selected, their QM properties are computed, and they are added to the training set. This loop continues until the model performs reliably across all relevant states. Tools like DP-GEN automate this process. [13]

• Validation and Deployment: The final model is validated by comparing MLFF-MD results against key experimental observables (e.g., reaction product distributions, activation energies) and QM data not seen during training. Successful validation allows deployment for large-scale production MD simulations. [13]

MLFFs in Practice: Applications and Protocols for Combustion

Case Study: EMFF-2025 for High-Energy Materials

The EMFF-2025 model exemplifies a specialized MLFF for C/H/N/O-based high-energy materials (HEMs), which are directly relevant to propellants and explosives in combustion. This model was developed using a transfer learning strategy, building upon a pre-trained model (DP-CHNO-2024) with minimal additional DFT data. [13]

Key Quantitative Results:

• Accuracy: Achieved a mean absolute error (MAE) for energy predominantly within ± 0.1 eV/atom and for force within ± 2 eV/Å, demonstrating chemical accuracy. [13]
• Performance: Successfully predicted crystal structures, mechanical properties, and thermal decomposition behaviors of 20 different HEMs. [13]
• Scientific Insight: Uncovered that most HEMs follow similar high-temperature decomposition mechanisms, challenging the conventional view of material-specific behavior. This demonstrates the power of MLFFs to reveal new chemical insights. [13]

Table 2: Key Research Reagents and Computational Tools for MLFF Development

Item / Resource	Type	Primary Function in MLFF Development
VASP / CP2K	Software	First-principles DFT code for generating quantum mechanical training data (energies, forces).
DP-GEN	Software	Automated active learning platform for generating robust and general-purpose MLFFs.
JAX-ReaxFF	Software	Framework for efficient, gradient-based optimization of ReaxFF parameters using ML. [33]
LAMMPS	Software	High-performance MD simulator that integrates with various MLFF formats for production runs.
Pyiron	Software	Integrated development environment (IDE) that manages the complete MLFF workflow. [58]
ReaxFFCHO-S22	Force Field	A specialized ReaxFF parameter set for hydrocarbon combustion, used for benchmarking. [10]
Cyclohexane / n-Dodecane	Chemical System	Common surrogate fuel molecules for validating combustion mechanisms. [10]

Detailed Protocol: Thermal Decomposition Simulation

The following workflow diagram and protocol outline the process of using an MLFF to study the thermal decomposition of a fuel molecule, a core process in combustion.

Title: MLFF Simulation and Analysis Protocol

Objective: To simulate the thermal decomposition and initial oxidation reactions of a hydrocarbon fuel (e.g., cyclohexane) using a trained MLFF and analyze the reaction mechanism.

Materials:

• Software: LAMMPS (with MLFF plugin) or equivalent MD engine; visualization/analysis tools (e.g., OVITO, MDAnalysis).
• Force Field: A validated MLFF for C/H/O systems (e.g., EMFF-2025 or a custom model).
• Initial Configuration: A periodic simulation box containing fuel molecules (e.g., 12 cyclohexane) and oxidizer (e.g., 108 O₂ molecules) at a target density (e.g., 0.2 g/cm³). [10]

Methodology:

• System Preparation and Equilibration:
  • Energy-minimize the initial configuration to remove bad contacts.
  • Perform a short (e.g., 10-50 ps) equilibration in the NVT ensemble to bring the system to the target simulation temperature.

• Production Reactive MD Simulation:

  • Switch to the NVE (microcanonical) ensemble to model an adiabatic reaction process.
  • Run the simulation for a duration sufficient to observe decomposition and reaction (e.g., 1 ns at 2800 K to accelerate kinetics). [10]
  • Save the atomic trajectory at frequent intervals.

• Data Analysis:

  • Species Population Analysis: Use bond analysis algorithms (e.g., based on cutoff distances or common neighbor analysis) to identify molecules and track the population of key species (e.g., C₂H₄, CH₂O, CO, CO₂, H₂O) over time. [10]
  • Reaction Pathway Identification: Manually inspect the trajectory or use automated reaction analysis tools (e.g., RxnAnalyzer) to identify sequences of bond-breaking and formation events. The initial step for cyclohexane is often homolytic C-C bond cleavage for ring-opening. [10]
  • Kinetic Extraction: From the species time profiles, extract apparent rate constants for the consumption of fuel or formation of major products.

Future Outlook and Strategic Implementation

The field of MLFFs is rapidly evolving. Key future directions include the development of foundation models for broad regions of chemical space, which can then be efficiently specialized for specific combustion systems via knowledge distillation, creating smaller, faster "simulation engines" that are up to 20 times faster than the original model. [59] Furthermore, the integration of multi-fidelity learning and uncertainty quantification will be critical for building reliable and trustworthy models for engineering design. [58]

For research teams aiming to integrate MLFFs into their combustion chemistry workflow, the following strategic steps are recommended:

• Start with a Pre-trained Model: Begin efforts with existing, well-validated foundation models (e.g., for C/H/N/O) to minimize initial data generation costs.
• Implement a Robust Active Learning Pipeline: Prioritize the establishment of an automated active learning workflow (e.g., using DP-GEN) to ensure the model's robustness and transferability.
• Validate Extensively Against Key Experiments: Always benchmark MLFF predictions against critical experimental data, such as ignition delay times and product species distributions, to establish credibility.
• Focus on Knowledge Distillation: For production-level MD simulations, invest in distilling large foundation models into faster, specialized potentials to achieve the required throughput for meaningful statistical analysis. [59]

Machine learning force fields represent a transformative advancement in computational combustion science. By effectively bridging the critical gap between the high accuracy of quantum mechanics and the practical scalability of classical force fields, they empower researchers to simulate complex reactive processes with unprecedented fidelity. While reactive force fields like ReaxFF have paved the way for studying bond-breaking events in MD, MLFFs are poised to supersede them for many applications by offering superior accuracy and, increasingly, competitive computational speed. As foundation models improve and specialization techniques mature, MLFFs will become an indispensable tool in the scientist's toolkit, driving forward the understanding and optimization of next-generation combustion systems.

Conclusion

The choice between ReaxFF and classical force fields is not a matter of superiority but of application. Classical force fields remain efficient for simulating non-reactive physical properties and equilibria, but ReaxFF is indispensable for modeling the complex, dynamic chemistry of combustion processes, such as bond breaking and radical reactions. While ReaxFF bridges a critical gap between quantum accuracy and classical efficiency, users must be aware of its limitations, including potential parametrization failures and high computational cost. Future directions point toward more accurate and transferable parametrizations, increased use of hybrid QM/ReaxFF/MM methods for multi-scale problems, and the exploration of machine learning force fields to further expand the scope and accuracy of combustion simulations. This progression will ultimately enable more predictive and reliable computational design in energy conversion and propulsion technologies.

References

• 1. molecular dynamics - What does "force field" mean? [https://mattermodeling.stackexchange.com/questions/7003/what-does-force-field-mean]
• 2. The ReaxFF reactive force-field: development, ... - Nature [https://www.nature.com/articles/npjcompumats201511]
• 3. Implementing reactivity in molecular dynamics simulations ... [https://www.nature.com/articles/s41467-024-50793-0]
• 4. Force fields for small molecules - PMC [https://pmc.ncbi.nlm.nih.gov/articles/PMC6733265/]
• 5. Evaluating the performance of ReaxFF potentials for sp2 ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC9465816/]
• 6. Molecular Forcefield Methods for Describing Energetic ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC8912029/]
• 7. Hybrid classical/machine-learning force fields for the accurate ... [https://pubs.rsc.org/en/content/articlehtml/2023/sc/d3sc04317g]
• 8. ReaxFF/AMBER—A Framework for Hybrid Reactive/ ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC8145783/]
• 9. ReaxFF [https://en.wikipedia.org/wiki/ReaxFF]
• 10. S22 Reactive Molecular Dynamics Simulations [https://www.sciencedirect.com/science/article/abs/pii/S174396712500371X]
• 11. Included force fields — ReaxFF 2025.1 documentation - SCM [https://www.scm.com/doc/ReaxFF/Included_Forcefields.html]
• 12. The ReaxFF reactive force-field: Development, applications ... [https://www.osti.gov/pages/biblio/1261375]
• 13. EMFF-2025: a general neural network potential for ... [https://www.nature.com/articles/s41524-025-01809-w]
• 14. Classical and reactive molecular dynamics: Principles and ... [https://www.sciencedirect.com/science/article/pii/S036012852300014X]
• 15. APPLICATIONS OF REAXFF FOR COMBUSTION ... [https://etda.libraries.psu.edu/catalog/23467adl233]
• 16. 6. ReaxFF Models — FitSNAP documentation [https://alphataubio.com/fitsnap/reaxff.html]
• 17. ReaxFF molecular dynamics simulation on the combustion ... [https://www.sciencedirect.com/science/article/abs/pii/S030101042400017X]
• 18. Force field format specification — ReaxFF 2025.1 ... - SCM [https://www.scm.com/doc/ReaxFF/ffield_descrp.html]
• 19. Reactive Molecular Dynamics Simulation of Fullerene Combustion ... [https://pubmed.ncbi.nlm.nih.gov/26606475/]
• 20. Development of force fields used in molecular dynamics ... [https://www.jsol-cae.com/en/tech-blog/20230518.html]
• 21. Modeling the potential energy surface by force for... fields [https://pubs.rsc.org/en/content/articlehtml/2025/sc/d5sc02715b]
• 22. Accuracy of classical for polyethylene structures away... force fields [https://link.springer.com/article/10.1557/s43579-023-00503-6]
• 23. Reactive force field (ReaxFF) molecular dynamics ... [https://www.sciencedirect.com/science/article/abs/pii/S0016236122004483]
• 24. A benchmark dataset for Hydrogen Combustion - PMC [https://pmc.ncbi.nlm.nih.gov/articles/PMC9114378/]
• 25. Extension of the ReaxFF Combustion Force Field toward ... [https://pubmed.ncbi.nlm.nih.gov/28072539/]
• 26. Extension of the ReaxFF Combustion Force Field toward ... [https://pure.psu.edu/en/publications/extension-of-the-reaxff-combustion-force-field-toward-syngas-comb/]
• 27. Force fields for classical molecular dynamics [https://pubmed.ncbi.nlm.nih.gov/23034750/]
• 28. ReaxFF molecular dynamics simulations of high-energy ... [https://www.sciencedirect.com/science/article/abs/pii/S0010218025005115]
• 29. 清华大学能动系和燃烧能源中心杨斌课题组肖航等人在空天 ... [https://www.cce.tsinghua.edu.cn/info/1012/4069.htm]
• 30. A study of pyrolysis reaction mechanisms for two different ... [https://www.sciencedirect.com/science/article/abs/pii/S016523702500422X]
• 31. HyChem Model [https://web.stanford.edu/group/haiwanglab/HyChem/]
• 32. Recent ReaxFF MD studies on pyrolysis and combustion ... [https://xlink.rsc.org/?doi=D2YA00285J&newsite=1]
• 33. Machine Learning Driven ReaxFF Optimization for ... [https://cjcp.ustc.edu.cn/hxwlxb/article/doi/10.1063/1674-0068/cjcp2503023]
• 34. Pre-training machine learned interaction potentials with ... [https://arxiv.org/html/2509.14205v1]
• 35. ReaxFF force field development for the thermal ... [https://www.sciencedirect.com/science/article/abs/pii/S000926142500510X]
• 36. ReaxFF and Hybrid Pair style - LAMMPS General Discussion [https://matsci.org/t/reaxff-and-hybrid-pair-style/43254]
• 37. ReaxFF - SCM [https://www.scm.com/amsterdam-modeling-suite/reaxff/]
• 38. Modeling the potential energy surface by force fields for ... [https://pubs.rsc.org/en/content/articlehtml/2025/sc/d5sc02715b?page=search]
• 39. An improved reactive force field parameter optimization ... [https://www.sciencedirect.com/science/article/abs/pii/S0927025625001193]
• 40. Improving force field accuracy by training against condensed ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC9254460/]
• 41. Modelling complex proton transport phenomena [https://arxiv.org/html/2501.04876v1]
• 42. Benchmarking the Performance of the ReaxFF Reactive ... [https://impact.ornl.gov/en/publications/benchmarking-the-performance-of-the-reaxff-reactive-force-field-o/]
• 43. Universal Machine Learning Potential for Systems with ... [https://arxiv.org/html/2508.15614v1]
• 44. Combustion simulation (ReaxFF) - SCM [https://www.scm.com/doc/Tutorials/MolecularDynamicsAndMonteCarlo/BurningMethane.html]
• 45. ReaxFF-MD study on the combustion characteristics and ... [https://www.sciencedirect.com/science/article/abs/pii/S0016236125025098]
• 46. Reactive force field molecular dynamics (ReaxFF-MD ... [https://www.sciencedirect.com/science/article/abs/pii/S0016236123027989]
• 47. Development of a ReaxFF reactive force-field modeling for ... [https://www.sciencedirect.com/science/article/abs/pii/S0169433225009213]
• 48. GitHub - tovrstra/ validation _eem_acks2_ reaxff : A reference... [https://github.com/tovrstra/validation_eem_acks2_reaxff]
• 49. (PDF) Development and Validation of a ReaxFF Reactive Force Field... [https://www.academia.edu/11494639/Development_and_Validation_of_a_ReaxFF_Reactive_Force_Field_for_Fe_Al_Ni_Alloys_Molecular_Dynamics_Study_of_Elastic_Constants_Diffusion_and_Segregation]
• 50. ReaxFF molecular dynamics studies on the impact of ... [https://www.sciencedirect.com/science/article/pii/S2666821125000134]
• 51. Accurate energy barriers for catalytic reaction pathways [https://www.nature.com/articles/s41524-023-01124-2]
• 52. Active learning of reactive Bayesian force fields applied to ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC9440250/]
• 53. Inhibition mechanism of CF 3 CHCl 2 (HCFC-123) on H 2 [https://www.sciencedirect.com/science/article/abs/pii/S0360319925023523]
• 54. Molecular Mechanics - PMC - PubMed Central - NIH [https://pmc.ncbi.nlm.nih.gov/articles/PMC4026342/]
• 55. Addressing pathological deficiencies in the simulation of ... [https://www.nature.com/articles/s41598-025-04482-7]
• 56. Included Forcefields — OldReaxFF 2025.1 documentation [https://www.scm.com/doc/OldReaxFF/Included_Forcefields.html]
• 57. Physics-informed machine learning for combustion: A review [https://arxiv.org/html/2509.03347v1]
• 58. Machine Learning Force Fields • IMSI [https://www.imsi.institute/activities/data-driven-materials-informatics/machine-learning-force-fields/]
• 59. Towards Fast, Specialized Machine Learning Force Fields [https://www2.eecs.berkeley.edu/Pubs/TechRpts/2025/EECS-2025-41.html]