We present a method based on kinetic molecular theory that identifies reactions of various molecularities in molecular dynamics (MD) simulations of bulk gases. The method allows characterization of the thermodynamic conditions at which higher than bimolecular reactions are a factor in the mechanisms of complex gas-phase chemistry. Starting with Bodenstein’s definition of termolecular collisions we derive analytical expressions for the frequency of higher molecularity collisions. We have developed a relationship for the ratio of the frequencies of termolecular to bimolecular collisions in terms of the temperature, density, and collision times. To demonstrate the method we used ReaxFF in LAMMPS to carry out MD simulations for NVT ensembles of mixtures of H2:O2 over the density range 120.2 kg m-3 to 332.7 kg m-3 and temperature range 3,000 K to 5,000 K. The simulations yield ReaxFF-based predictions of the relative importance of termolecular collisions O2∙∙∙H2∙∙∙O2 and bimolecular collisions O2∙∙∙H2 in the early chemistry of hydrogen combustion.