Also check that your force field satisfies the following physical constraints:
- ATM r_vdw ≥ r_s ≥ r_p ≥ r_pp
- ATM r_vdw > rcore2
- ATM ecore2 > epsilon
- ATM eta > 7.2 gamma (QEQ), eta > 8.13 gamma (ACKS2)
- ATM bcut_acks2 ≥ max(r_s)
- electronegativity hierarchy ie. chi(O) > chi(N) > chi(C) > chi(H)
- BND De_s ≥ De_p ≥ De_pp
- BND p_be2 ≥ p_bo2, p_bo4, p_bo6
- BND p_bo6 ≥ p_bo4 ≥ p_bo2
- OFD r_pp ≥ r_p ≥ r_s
- OFD.X.Y.alpha ≥ max(ATM.X.alpha, ATM.Y.alpha)
- ANG p_pen1 ≥ |p_val1|
- ANG p_pen1 ≥ |p_val2|
- TOR V1 ≥ |V2|
- TOR V1 ≥ |V3|
- HBD r0_hb > r_s (ensure H-bond length exceeds σ-bond radius)
#4 is the most important, its called the “polarization catastrophe” [Troubleshooting and warnings — ReaxFF 2025.1 documentation]. the closer you get to that bound the more ill-conditioned the QEQ matrix gets and the more likely you are to have some numerical issues with your simulation.
a lot of force fields in the past have cheated on those physical constraints to overfit to a very specific system under study in a given paper. if the system you are simulating doesnt match almost exactly the training data of the original paper then that force field is likely not transferable to your new application.
one more word of caution, check that pi or pi-pi coeffs are properly trained in your force field if they are relevant in your application. in my personal experience, there are many potentials with C/H/O/N atoms but not all have the pi-bond parameters trained so a benzene molecule might behave in a completely unphysical manner, ie. hydrolyze in a few femtoseconds when benzene has a half-life in water measured in days or even weeks.