Large scale molecular dynamics (MD) simulations are performed to study the oxidation of highly oriented pyrolytic graphite (HOPG) by hyperthermal atomic oxygen beam (5 eV). Simulations are performed using the ReaxFF classical reactive force field. We present here additional evidence that this method accurately reproduces ab initio derived energies relevant to HOPG oxidation. HOPG is modeled as multilayer graphene and etch-pit formation and evolution is directly simulated through a large number of sequential atomic oxygen collisions. The simulations predict that an oxygen coverage is first established that acts as a precursor to carbon-removal reactions, which ultimately etch wide but shallow pits, as observed in experiments. In quantitative agreement with experiment, the simulations predict the most abundant product species to be O2 (via recombination reactions), followed by CO2, with CO as the least abundant product species. Although recombination occurs all over the graphene sheet, the carbon-removal reactions occur only about the edges of the etch pit. Through isolated defect analysis on small graphene models as well as trajectory analysis performed directly on the predicted etch pit, the activation energies for the dominant reaction mechanisms leading to O 2, CO2, and CO product species are determined to be 0.3, 0.52, and 0.67 eV, respectively. Overall, the qualitative and quantitative agreement between MD simulation and experiment is very promising. Thus, the MD simulation approach and C/H/O ReaxFF parametrization may be useful for simulating high-temperature gas interactions with graphitic materials where the microstructure is more complex than HOPG.