Classical MD models as implemented in highly popular MD software such as Amber , LAMMPS , GROMACS and NAMD are based on the assumption of static chemical bonds and, in general, static charges. Therefore, they are not applicable to study phenomena where chemical reactions and charge polarization effects play a significant role. To address this gap, reactive force fields (e.g., ReaxFF , REBO , Tersoff ) have been developed. These bond order potentials allow bonds to form and break throughout the simulation and they can dynamically assign partial charges to atoms using suitable charge models such as the electronegativity equalization method (EEM) . The functional forms for reactive potentials are significantly more complex than their classical counterparts due to the presence of dynamic bonds and charges. For instance, ReaxFF has a formulation that contains more than 100 parameters for a simulation with 3 elements, and is about two orders of magnitude more expensive than a typical Lennard-Jones potential. Consequently, training reactive force fields is an even more difficult task due to the need to capture complex phenomena such as charge distributions and reactions, and due to the large number of parameters involved.