Transferable Reactive Force Fields: Extensions of ReaxFF-lg to

Subscriber access provided by UNIVERSITY OF SOUTH CAROLINA LIBRARIES

Article

Transferable Reactive Force Fields: Extensions of ReaxFF-lg to Nitromethane James P Larentzos, and Betsy M. Rice J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.6b11761 • Publication Date (Web): 08 Feb 2017 Downloaded from http://pubs.acs.org on February 8, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Transferable Reactive Force Fields: Extensions of ReaxFF-lg to Nitromethane James P. Larentzos* and Betsy M. Rice U.S. Army Research Laboratory Aberdeen Proving Ground, Maryland 21005, USA Abstract Transferable ReaxFF-lg models of nitromethane that predict a variety of material properties over a wide range of thermodynamic states are obtained by screening a library of ~6600 potentials that were previously optimized through the Multiple Objective Evolutionary Strategies (MOES) approach using a training set that included information for other energetic materials composed of carbon, hydrogen, nitrogen and oxygen. Models that best match experimental nitromethane lattice constants at 4.2 K, 1 atm are evaluated for transferability to high pressure states at room temperature and are shown to better predict various liquid and solid-phase structural, thermodynamic and transport properties as compared to the existing ReaxFF and ReaxFF-lg nitromethane parameterizations. Although demonstrated for an energetic material, the library of ReaxFF-lg models are supplied to the scientific community to enable new research explorations of complex reactive phenomena in a variety of materials research applications. 1. Introduction: The emergence of force fields that describe bond breaking and bond formation for use in reactive molecular dynamics (RMD) simulations 1-7 has made possible atomic-scale investigations of complex physical and chemical phenomena that cannot be easily obtained through experimental means. Probably the most widely-used reactive force field is ReaxFF, developed in the early 2000s and implemented in the massively-parallel molecular dynamics software package LAMMPS 8. Since the inception of ReaxFF, a modification to the original force field (denoted ReaxFF-lg 9) and different parameterizations have been obtained for materials and processes for which the original parameterization was not transferable (see the comprehensive review by Senftle et al 10 for the complete history of the ReaxFF formalism and how it has evolved to its current state). We, along with others 11-13, have developed algorithms that would allow for reparameterization of this force field to different fitting sets 14. Our approach, the Multiple Objective Evolutionary Strategies (MOES), was previously applied to parameterize ReaxFF-lg reactive potentials suitable for molecular dynamics (MD) simulations of energetic molecular crystals 14-15. In our previous study, a subset of the 600+ ReaxFF-lg 9 adjustable parameters that affect the non-bonded van der Waals, electrostatic, hydrogen bond, and dispersion interactions were systematically varied through an evolutionary strategies algorithm while simultaneously optimizing multiple objective functions defined within the training set. The training set consisted of relevant ab initio structural and reaction pathway energetic data of various chemical systems, including high energy explosive materials such as 1,3,5-trinitro-1,3,5-triazine (RDX), 1,1-diamino-2,2-dinitroethene (FOX-7), octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine (HMX), nitromethane (NM), pentaerythritol tetranitrate (PETN), 1,3,5-triamino-2,4,6-trinitrobenzene (TATB), and triacetonetriperoxide (TATP). The training set was subdivided into physically relevant objective functions that consisted of structural, energetic or chemical properties of the materials. A complete description of the training set, adjustable parameters, and objective functions was described in our original publication 14. Multiple objective optimizations do not yield a single solution (i.e., parameterization), but rather are used to populate the solutions that lie on an N-dimensional Pareto surface in objective space, where N is the number of objective functions. A Pareto-optimal solution is defined by the condition that no objective function can be improved without making another objective function worse. Thus, all Pareto-optimal solutions are considered to be equally good although there is a trade-off in better satisfying one objective over another. The computationally intensive search through parameter space to find the Pareto-optimal

*Corresponding Author. Email: [email protected]. Tel: (410) 306-0809. ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

solutions is only performed once for a given set of adjustable parameters and training set. Our earlier study produced ~6600 ReaxFF-lg Pareto-optimal solutions (i.e., candidate force fields). The ~6600 parameterizations derive from seven independent Pareto surfaces, six of which correspond to the ReaxFFlg models presented in Table 2 of Ref. 15 (Models 3-8), and the seventh corresponds to a ReaxFF-lg model that was subsequently obtained through a 6-objective optimization that varied the 46 van der Waals and electrostatic parameters referenced. The term “model” used in this Table denotes the ReaxFF or ReaxFF-lg functions that were optimized using N={5,6} objective functions, while varying different numbers of parameters (i.e., M={8, 46, 54, 62} of the long-range ReaxFF-lg van der Waals, electrostatic and/or dispersion parameters). The objective functions were composed of charges, geometries, heats of formation, reaction pathway energies and dimer energies. We demonstrated the most promising of the force fields at predicting the ambient state crystallographic parameters of the energetic molecular crystal RDX in isostress-isothermal simulations, and identified three ReaxFF-lg parameterizations that performed as well as or better than the original ReaxFF-lg parameterization in describing RDX 14. We used these newly-developed force fields to explore transferability, limited only to ambient crystal structures, to other energetic molecular crystals HMX, 2,4,6,8,10,12-hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane (CL-20), PETN, FOX-7 and TATB 15. Unfortunately, none of these models could describe FOX-7 and TATB, demonstrating that the parameterizations had very limited transferability to a few systems at a single state point (ambient). While disappointing, this is not surprising, since the selected force fields were based on their ability to predict RDX at a single state point. Had the selection criteria been expanded to include other materials, different force fields might have been chosen that demonstrated better transferability. Nonetheless, our earlier efforts generated thousands of candidate force fields that might be transferable across chemical systems and/or thermodynamic states that were not included in the fitting set. Of particular interest to us are force fields that will predict dynamic response of energetic materials subjected to thermal or shock initiation when used in molecular simulation. In these simulations, the energetic material samples a wide range of temperatures and pressures as it undergoes rapid, exothermic chemical reactions. Therefore, obtaining adequate force fields to describe the variety of material states in such simulations is a challenge. To explore the possibility that our collection of candidate force fields might contain transferable force fields for use in such explorations, we sampled the same library of optimized potentials produced from our previous work and identified appropriate parameterizations for the prototypical energetic molecular system NM over a wide range of thermodynamic states. Our choice of NM is predicated on our current efforts to investigate shear-induced chemical reaction of an energetic material under pressure 16. The advantages of choosing NM are its small size, chemical simplicity, solid state phase stability over a large pressure range, and availability of a wealth of experimental information over a wide span of temperatures and pressures (including shock). In this work, we describe a two-stage selection process in which ~6600 force fields are screened, and a subset are then subjected to a second down-selection process in which the models are assessed in ability to predict pressure-dependent lattice parameters. After this, we show a rigorous assessment of transferability of the models by calculating properties not used in the down selection process, including equilibrium properties in the liquid and solid phases, solid-liquid phase coexistence (i.e., melting point), shock Hugoniot equation-of-state (EOS), transport properties and reaction chemistry. Identification and validation of ReaxFF-lg models that well describe NM over a wide range of thermodynamic conditions is not the singular goal of this paper. A second goal of equal importance is to provide the scientific community with the library of optimized ReaxFF-lg parameterizations for use in other applications that depend on user needs. The library of parameterizations that lie on the Pareto surfaces have been trained to reproduce structural, energetic and chemical properties of hydrocarbons, nitramines as well as various energetic materials; thus, these parameterizations could be applied to other C-H-N-O chemical systems. While distributing these force fields will not result in the complete elimination of the need to reparameterize reactive force fields, it is hoped that the need will be reduced by providing sufficiently accurate force fields for other applications. The work here provides a

ACS Paragon Plus Environment

Page 2 of 32

Page 3 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


straightforward approach to sample and assess the utility of of ReaxFF-lg force fields within an existing library. The paper is structured as follows. We first describe the model selection process in Section 2, followed by a detailed examination of the molecular structural parameters, liquid-phase densities, melting points, shock Hugoniot EOS, and reaction chemistry predicted by the best models in the remaining sections. 2. Computational Details: At ambient conditions, NM exists in the liquid-phase and crystallizes at a temperature of ~244 K 17. Trevino et al 18 solved the crystal structure of NM at 4.2 K and 1 atm through single crystal x-ray diffraction and neutron powder diffraction. The crystal structure is orthorhombic with space group P212121, lattice constants a = 5.1832 Å, b = 6.2357 Å, c= 8.5181 Å and Z = 4 molecules in the unit cell; a representation is shown in Figure 1. Additional experimental studies have since been performed to solve the crystal structures at 298 K and pressures ranging from 1.1 GPa to 15.3 GPa 19-21 through x-ray diffraction techniques.

Figure 1: Representation of the NM crystal unit cell determined by Trevino et al. 18 at 4.2 K, 1 atm. Atom labels are consistent with those used to describe H-C-N-O dihedral angles hereafter. Numerous studies have examined NM using methods that allow for simulation of reactivity, i.e. ab initio molecular dynamics 22-31 and reactive molecular dynamics (RMD) simulations 9, 32-40. Various parameterizations of ReaxFF 40-44 and the ReaxFF-lg dispersion-corrected extension 9 have been applied to RMD simulations; however, the quality of these parameterizations with regard to describing NM has not been extensively assessed. In this work, our aim is to search our library of Pareto-optimal parameterizations to find a suitable reactive model of NM that is transferable over the thermodynamics states measured experimentally and compare them to the existing ReaxFF 43 and ReaxFF-lg 9 parameterizations. We limit our search to the ReaxFF-lg models, mainly because our previous studies 15 have indicated that we achieved better overall structural models when the lg dispersion term is included. The complete library of Pareto-optimal ReaxFF-lg parameterizations is summarized and provided in the Supporting Information (Table S1).



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2.1 First stage screening: A MD screening-phase simulation is conducted for each of the ~6600 ReaxFF-lg parameterizations to compute the equilibrium crystal structure of NM at 4.2 K and 1 atm. The model NM structure used in the simulation is a supercell composed of 6x5x4 unit cells. The velocities of the atoms are initialized according to a Gaussian distribution at a temperature of 4.2 K. The screening-phase consists of a 2.5 ps MD simulation in the NVE ensemble with velocity rescaling at every step to a temperature of 4.2 K, followed by a 2.5 ps NVT simulation and subsequent 10.0 ps anisotropic-NPT simulation (1 time step = 0.5 fs). Damping parameters are 1.0 ps for the barostat and 0.1 ps for the thermostat. The time at which the system reaches an equilibrated state is determined by block-averaging the anisotropic-NPT portion of the equilibration phase through the Flyvbjerg-Petersen method 45-46. The lattice constants are timeaveraged over the period from which the system reaches the equilibrated state until the end of the 15 ps equilibration phase, then compared with experiment at 4.2 K to compute an absolute average lattice constant error. In total, there were ~340 parameterizations that provided smaller average lattice constant errors than both the ReaxFF and ReaxFF-lg models. This simple, yet efficient screening-phase consumed a total of ~160,000 CPU-hours using multi-core parallel processing on a Cray XC30 supercomputer, confirming the feasibility to exhaustively screen the library of optimized potentials through a brute-force procedure with modern supercomputing power. It is to be noted that it is the user’s decision as to what properties should be computed in the initial screening phase and what metrics should define the “best” parameterizations. The criteria can be as simple or as complex as the user chooses. For this work, we chose to examine structural properties and used the absolute average lattice constant error as the metric to evaluate the parameterizations. Choosing a slightly different metric, such as a root mean square error, would no doubt provide a different subset and/or ranking of models to evaluate. Alternatively, one could choose to screen with a different property (or multiple properties), which should be chosen to align with the goals of the application. The focus in this work is not to determine the optimal screening properties or metrics that should be used, but rather demonstrate the ease in which reasonable, transferable parameterizations can be obtained for use in predictive simulations. 2.2 Second stage down-selection: High pressure experimental x-ray diffraction data19-21 for solid NM at 298 K and pressures ranging from 1.1 GPa to 15.3 GPa 19-21 are used in the second stage down-selection process. For each of the seven Pareto surfaces, the three models that best predicted the lattice constants after the first screening stage (a total of 21 models) have been selectively chosen for further evaluation and comparison with experimental structural data at higher pressure states (see Supporting Information, Table S2 for more information). In principle, any of the ~6600 Pareto-optimal solutions or ~340 solutions identified in the first stage screening could be selected for testing during this second stage, and may lead to even better predictive models than are identified in this manuscript. Rather than attempt to do an exhaustive search over all optimized solutions (being mindful of time and computational resource constraints), our focus is to demonstrate a reasonable pathway to obtain a model to reliably predict the properties of interest for a given application. The simulated structures were obtained by restarting the 15.0 ps screening-phase calculations and slowly ramping the temperature and pressure in succession to the experimentally-measured state points at 298 K, 1 GPa, 2 GPa, 3.5 GPa, 6 GPa, 7.8 GPa, 9.5 GPa, 11.6 GPa, 13.4 GPa and 15.3 GPa.. At each state point of interest, equilibration and production phase anisotropic-NPT simulations (1 time step = 0.25 fs) were conducted for at least 50 ps prior to ramping to the next higher pressure state. Upon completion of the 50 ps NPT calculations at each state point, an additional, independent MD calculation was conducted for 50 ps in the isothermal-isostress (NsT) ensemble to allow volumetric and cell vector fluctuations of the


Page 4 of 32

Page 5 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


simulation cell. The predicted lattice constants and cell angles were obtained by averaging over the final 37.5 ps of the NsT simulation. Each parameterization results in an orthogonal crystal structure, where on average, the cell angles (α, β, γ) are 90.0 ± 0.01°. A residual sum of squares (RSS) of the a, b can c lattice constants as compared to the experimental lattice constants (, and ) over the entire pressure range is defined as:

= 〈〉 − + 〈〉 − + 〈〉 − ̅

(1)

where the angle brackets denote the ensemble average of the respective lattice constants. The 21 MOES-optimized parameterizations were scored based on the RSS at each state point and compared to the ReaxFF and ReaxFF-lg parameterizations. Parameterizations 6-46-lg-0498, 6-46-lg0524 and 6-54-lg-0077 (hereafter referred to as Models 0498, 0524 and 0077) have the smallest RSS errors and are further evaluated for the remainder of this work; lattice constants over the entire pressure range predicted by these models are compared to experiment as well as the ReaxFF and ReaxFF-lg models and are shown in Figure 2 (left frame). Percent errors of the respective lattice constants for all of the models are also shown in Figure 2 (right frame). These figures clearly show that Models 0498, 0524 and 0077 have an overall significant improvement over the ReaxFF and ReaxFF-lg parameterizations for each lattice constant over the range of pressures, especially for lattice constant b. Between the two original parameterizations, ReaxFF has a consistently better performance in predicting lattice constants a and c than ReaxFF-lg, whereas ReaxFF-lg better predicts lattice constant b for pressures less than 11 GPa. Models 0498, 0524 and 0077 predicted lattice constants to within 3.8% of experiment over the entire pressure range. Conversely, 19 out of 27 predictions made by ReaxFF-lg have percent deviations from experiment greater than 3.8%. ReaxFF performs slightly better than ReaxFF-lg, with 12 of the 27 predictions being greater than 3.8%. The dependence of the lattice constants and percent errors on pressure and temperature is provided in Tables S3-S10 and Figure S1 of the Supporting Information.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 32

Figure 2: Lattice constant dependence on pressure (left frame) and % error (right frame). All parameterizations are fitted to the Murnaghan equation 21, 47: = − 1

(2)

where P(V) is the pressure at volume V, V0 is the volume at P = 0 GPa, B0 is the bulk modulus at zero pressure, and is dB0/dP. Because NM exists as a liquid at 298 K, 0 GPa, a 3rd-order polynomial is fitted to the P-V data collected at 298 K, 1.0-15.3 GPa and used to extrapolate to zero pressure. Table 1 compares the extrapolated lattice constants and volume at 298 K, zero pressure (V0) to the fitted volume and bulk modulus obtained through a non-linear least squares analysis. Table 1: The computed lattice constants, cell volume, bulk moduli and derivative values at 298 K, 0 GPa Parameterization a0 (Å) b0 (Å) c0 (Å) V0 (Å3) V0-fit (Å3) B0 (GPa) Model 0498 5.192 6.223 8.702 281.17 285.4 11.1 4.8 Model 0524 5.177 6.200 8.667 278.03 287.5 8.7 5.2 Model 0077 5.172 6.509 8.712 293.27 309.2 6.3 5.2 ReaxFF-lg 4.992 6.591 8.195 269.64 292.7 3.5 6.9 ReaxFF 5.038 7.013 8.345 294.86 342.6 2.7 5.5 18 Trevino et al 5.270 6.375 8.832 296.72 Cromer et al 21 5.227 6.330 8.758 289.77 292.7 7.0 5.7 Yarger and Olinger 20 5.197 6.292 8.747 286.02 291.9 9.1 6.0 Model 0524 is the only parameterization that falls within the range of bulk modulus values reported experimentally. Model 0077 moderately underestimates the lower experimental value, while ReaxFF


Page 7 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and ReaxFF-lg have a larger underestimation. Model 0498 overestimates the upper experimental measurement, but agrees very well with the predicted lattice constants and cell volume at zero pressure. 3. Results: The remainder of the paper will describe the ability of the three MOES-optimized models (i.e., models 0498, 0524, 0077), as well as ReaxFF-lg and ReaxFF to predict other properties that were not used in either the fitting or selection. 3.1 Molecular crystallographic and structural assessment: The time- and spatially-averaged molecular crystallographic and molecular structural details for each molecule in the unit cells were predicted using models 0498, 0524 and 0077, ReaxFF-lg and ReaxFF over a range of temperatures and pressures. These results were obtained from system configurations recorded at 0.25 ps intervals during each production trajectory and compared with experimental information using the procedures described in Section 2.5 of our previous publication 15. The information, given in Supporting Information Tables S4-S8 and S11-S16, includes the differences in atomic positions of each predicted molecule from its experimental counterpart due to molecular deformation and to relative orientation within the unit cell and fractional positions of the mass centers (sx, sy, and sz) for each thermodynamic condition studied here. Unlike Ref. 48, however, we performed these structural analyses for the heavy atoms (i.e., C, N, O) only, as we observed that the methyl groups for the molecules within the supercell undergo rotations (sometimes freely or slightly hindered) about the C-N bonds at some of the thermodynamic conditions we studied. In such cases, our method of analyses would produce incorrect and misleading time- and spatially- averaged positions for the hydrogen atoms. Instead, we have compiled cumulative distributions of the Hi-C1-N2-O3 (i = 5-7) dihedral angles calculated from the molecular configurations in each snapshot of the production trajectories, to provide a comparison between the predicted equilibrium orientation of the methyl group with the experimental structure. The pressure dependence of the rms deviation of the atomic positions of the predicted structures for all models from those of the experimental structures due to molecular structural differences only is shown in Figure 3 (top-left frame). Also, the range of values for the CN and NO bonds and the ONO angle predicted by the models at all pressures and temperatures studied is given for comparison with experiment in Supporting Information Table S17 and S18. The ReaxFF model had a very large rms deviation of molecular structure in the unit cell, with the most significant molecular structural difference being in the N-O bonds and O-N-O angles. The experimental values for the O-N bonds ranged from 1.17 to 1.25 Å, while ReaxFF predicted values ranging from 0.873 to 1.056 Å. Also, experimental values for the O-N-O angles ranged from 121.4 to 126.6º, whereas the ReaxFF values ranged from 103.7 to 112.3°. The molecular structural parameters for all other models at all pressures were in fairly good agreement with experiment, with Model 0498 performing slightly worse, as reflected in rms deviations ranging from 0.03 to 0.06 Å. Interestingly, the ReaxFF-lg and ReaxFF models were in closer agreement with experiment for the CN bond distances; however, overall, Model 0077 appears to have the smallest deviation from experiment due to molecular deformation.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3: Pressure and temperature dependence of the rms deviations of predicted atomic positions from experiment due to molecular deformation (top frames); the rms deviations of predicted atomic positions from experiment due to relative orientation in the unit cell (middle frames); the absolute error of the positions of the center-of-mass fractional positions, averaged over molecules and dimensions (bottom frames). Data for ReaxFF at 228K are not included (see text). The pressure dependence of the orientational rms deviation of the atomic positions of the predicted structures from experiment is shown in Figure 3 (middle-left frame). This figure indicates that the structures predicted by ReaxFF and ReaxFF-lg have a substantially larger reorientation of the molecules within the unit cell than that of the other models. Similarly, there is a larger deviation of the locations of the mass centers (in fractional units) in the unit cells for the ReaxFF and ReaxFF-lg models compared to those predicted by the other models, as shown in Figure 3 (bottom-left frame). The small absolute average deviations of the fractional positions for these models, especially Model 0498, indicate that the locations of the molecular mass centers within the unit cells are fairly close to the experimental positions at all pressures. The temperature-dependent rms deviation of predicted atomic positions from experiment due to structural deformation and orientation is shown in Figure 3 (right frames) with all models having similar error in molecular structure across the temperature range, except for ReaxFF at 228K. In simulations at 228 K using the ReaxFF model, the NO2 groups of the molecules in the simulation cell are freely rotating about the C-N bond, resulting in predicting incorrect time- and spatially-averaged atomic positions for the oxygens. Therefore, this structure is not included in the analyses. Differences in molecular orientation of the molecules in the cells from experiment are similar in magnitude to the pressure-dependent differences for Models 0524 and 0498, and ReaxFF-lg and ReaxFF. Model 0077, however has a larger deviation


Page 8 of 32

Page 9 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


across the temperature range, and is intermediate between the low-error 0524 and 0498 models and the original ReaxFF and ReaxFF-lg parameterizations. A similar overall behavior of all of the models is evidenced in the average error in position of molecular mass centers within the unit cells; the ReaxFF and ReaxFF-lg models are in larger disagreement from experiment than the remaining models, but Models 0524 and 0498 have the smallest error, and Model 0077 is intermediate. The temperature and pressure-dependent behavior of the methyl groups can be determined from the cumulative distributions of the Hi-C1-N2-O3 (i=5,6,7) dihedral angles compiled during the trajectory integrations. Plots of the distributions as functions of temperature (upper frames) and pressure (lower frames) for Model 0077, ReaxFF-lg, and ReaxFF are provided in Figure 4. Focusing first on the temperature dependence of the methyl group orientation, at 4.2 K, 1 atm, the distributions for all models have three peaks that are separated from one another by ~120º, which is consistent with “three nuclei related to one another by a threefold rotation about the C-N bond” [quote from Trevino and Rymes 49]. The locations of the three peaks in the distributions for Models 0524 (-89º, 33º, 150º) and 0498 (-88º, 32º, 150º) are in almost exact agreement with the experimental values (-89º, 32º, 151º), while those of Model 0077 (-96º, 25º, 143º), ReaxFF-lg (-98º, 24º, 141º) and ReaxFF (-97º, 24º, 142º) are within 10º of the experimental values. Also, consistent with experiment observations 18, the distributions broaden with temperature, indicating an increase in librational motion about the equilibrium position with increasing temperature. By 228 K, the distribution for ReaxFF indicates that the methyl group are freely rotating, while the distributions for Models 0077, 0524 and 0498 indicate a very slightly-hindered methyl group rotation at this temperature. The distribution for ReaxFF-lg at 228K shows only infrequent methyl group reorientation.

Figure 4: The temperature-dependent (top frames) and pressure-dependent (bottom frames) H-C-N-O dihedral angle distributions. Models 0498 and 0526 are omitted for clarity, but are similar to Model 0077. The dashed, vertical lines are provided to illustrate peak shifts with pressure. With regard to the pressure-dependent orientation of the methyl groups, Cromer et al. 21 suggest that upon conversion of NM to a solid at 0.3 GPa, the methyl groups freely rotate about the C-N bond at pressures less than 0.6 GPa. They also suggest that for pressures ranging from 0.6 to 3.5 GPa, the methyl group rotation is slightly hindered. At 3.5 GPa, the methyl group becomes fixed, and is rotated about the C-N bond ~45º from its orientation at 4.2 K, 1 atm. Citroni et al 31 describes this reoriented structure as “eclipsed” relative to the low-temperature, low-pressure “staggered” configuration, and is singularly observed at 11.0 and 15.0 GPa. The results of Citroni et al 31 at 1.1 and 3.2 GPa were best described by a system of molecules having methyl groups distributed in both the staggered and eclipsed orientations.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 32

Although the 7.6 GPa model could also be described by such a system having “2-fold disordered methyl orientation”, it was equally well described by a model wherein the methyl groups of all molecules were oriented in the eclipsed position. While all of the models show slightly hindered methyl group rotations at 1 and 2 GPa, only Models 0077 and ReaxFF-lg show that at pressures of 3.5 GPa and higher, the methyl groups are essentially fixed in place. The methyl groups in Models 0498 and 0524 are fixed at pressures of 6 GPa and higher, but the methyl groups in ReaxFF undergo hindered rotation at all pressures. All of the models except for ReaxFF have only three peaks in their distributions at all pressures, with the locations of the peaks at 1 GPa close to the values at 4.2 K 1 atm, and, with the exception of Model 0077, shift with increasing pressure. Estimates of the location of the peaks for ReaxFF-lg, Models 0524 and 0498 at 15 GPa show shifts of only ~3-11° from the low temperature values. None of these models show a transition to the eclipsed configuration observed in experiment. The cumulative distributions for ReaxFF are similar to the other models for pressures up to 3.5 GPa, except the peaks are much broader, and there is a greater degree of methyl group libration/rotation about the C-N bond. The distributions for this model are distinctly different from those of the other models at pressures above 3.5 GPa; they have six peaks rather than three, with the smaller peaks corresponding to molecules in an eclipsed position and the larger peaks corresponding to molecules in the low-temperature the staggered position. 3.2 Liquid phase property assessment: The liquid-phase equilibrium and transport properties are evaluated for all parameterizations. Isothermalisobaric (NPT) ensemble simulations are conducted over a temperature range from 255 K to 400 K at 0 GPa. The simulation cell volume is permitted to fluctuate isotropically and is controlled through a NoséHoover thermostat and barostat with damping parameters of 0.1 ps and 1.0 ps, respectively. Initially, the structure was equilibrated for 50 ps at the targeted state point, followed by a production-phase simulation over which the thermodynamic properties were sampled every 0.1 ps for 100 ps. The equilibrium densities are computed for all parameterizations over the 100 ps production-phase simulations. Additional 100 ps NVE simulations are conducted at the equilibrium temperatures and densities to compute the self-diffusion coefficients. The liquid-phase density dependence on temperature is shown in Figure 5 and compared with experimental measurements 50. The experiments were performed at temperatures ranging between 293 K and 353 K. The ReaxFF-lg parameterization overestimates the liquid-phase densities over the entire temperature range by >25%, on average. The densities obtained from Model 0498 and Model 0524 are also larger than experiment by on average ~7.8% and ~9.0%, respectively, over the entire temperature range. Model 0077 is in good agreement with experiment, underestimating the experimental liquid-phase density experiments by ~1.7% over the entire temperature range. ReaxFF also underestimates the liquid phase density by on average ~8.1%. Overall, all parameterizations qualitatively capture the trend of linearly decreasing density with temperature, as shown in Figure 5.


Page 11 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5: Liquid-phase density dependence of NM on temperature of the various models compared with experiment. The self-diffusion coefficient was computed as function of temperature and is shown in Figure 6 for each parameterization and compared to the 17O NMR diffusion experiments of Price et al 51 and 1H and 2H NMR studies of Holz et al 52. The self-diffusivity was computed through the Einstein relation by computing the mean square displacement of the liquid phase trajectories at temperatures ranging from 255 K to 400 K. Overall, the predicted values from ReaxFF, ReaxFF-lg, Model 0498 and Model 0524 simulations are smaller than the experimental measurements, whereas Model 0077 values are slightly larger than the experimental measurements. The temperature dependence can be described through an Arrhenius relationship. The activation energies are determined by fitting linear trend lines to the data presented in Figure 6 and measuring the slopes (see Table 2 for the comparison with experiment). Models 0077, 0498, 0524 and ReaxFF best match the self-diffusion of NM, while ReaxFF-lg deviates significantly from experiment. Table 2: Comparison of the activation energies for the self-diffusion coefficients Parameterization Ea (kJ/mol) Model 0498 11.50 Model 0524 9.22 Model 0077 8.36 ReaxFF-lg 17.66 ReaxFF 11.69 Experiment - Price et al 51 9.85 Experiment - Holz et al 52 10.17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 32

Figure 6: Semi-log plot of liquid-phase NM self-diffusivity as a function of temperature compared with experiment 51-52. 3.3 Melting point assessment: The best-performing parameterizations were further tested by comparing the melting point with the experimental measurement of 244.73 K at 1 atm of pressure 17. In this work, the direct melting point method of Morris et al 53-54 is applied to efficiently screen the quality of the candidate ReaxFF-lg parameterizations of NM. In this method, the location of the solid-liquid interface is monitored directly at various temperatures over the timespan of the simulation. Temperatures above the melting point result in continued melting as the interface moves into the solid phase, whereas temperatures below the melting point result in continued crystallization that diminishes the liquid phase. At the melting point, the solid and liquid phases coexist in equal proportion and exhibit a constant temperature profile. A 3D-periodic model is constructed by creating a 6 x 5 x 8 rectangular supercell with dimensions 31.098 x 31.18 x 68.144 (960 NM molecules). The solid-liquid interface is created by freezing the molecules on one-half of the simulation cell to their lattice positions, while permitting the molecules in the other half to move freely without constraint. The sample is heated to 298 K, 1 atm for 125.0 ps in the NPT ensemble (dt = 0.5 fs), thus creating a solid-liquid interface at the center of the simulation cell and a second interface at the edge of the simulation cell as a result of periodic boundary conditions. The constraints on the frozen molecules are then removed and the system is further equilibrated at 1 atm of pressure and the target temperature (ranging from 200K to 290 K in 5 K increments) over a 75 ps time period (dt = 0.25 fs) through the following simulation protocol: 12.5 ps NVE ensemble with temperature rescaling followed by 12.5 ps in the NVT ensemble and 50 ps in the NPT ensemble. After equilibration, production calculations are performed in the NVE ensemble for a total of 200 ps. The temperature of the system is monitored and the melting point is estimated by identifying the system that exhibits a constant temperature profile over the timespan of the simulation. The melting point ranges for each parameterization are provided in Table 3. While none of the parameterizations provided a precise prediction of the experimental measurement of 244 K, we note that


Page 13 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the uncertainty in the predictions from the direct melting method can be rather large, and thus we cannot determine unambiguously whether the imprecision is due to the model, the method, or a combination of both. More rigorous (and resource intensive) free energy based approaches are able to precisely compute the thermodynamic melting points55-56 and are more appropriate for evaluation of the models. In this work, estimates of the melting point are obtained from a rather simple, efficient and straightforward approach and are shown to be within ~45 K of the experimental measurement for all models. Table 3: Estimated melting points for various NM parameterizations. Parameterization Melting Point Model 0498 270-275K Model 0524 210-215K Model 0077 195-200K ReaxFF-lg 260-265 K ReaxFF 210-215 K Experiment 17 244 K 3.4 Shock Hugoniot assessment: The Hugoniot curve is the ensemble of accessible states that a system can attain in sustaining a shocked state. The initial unshocked and shocked states are thermodynamically linked by the mass, momentum and energy conservation relations: 1 ! − ! = − $ − $ 2

(3)

%&' − &( = % &'

(4)

− = % &' &(

(5)

where ! the specific internal energy (sum of the kinetic and potential energy), ) is the pressure, % is the specific density, $ = * is the specific volume, and the subscript 0 refers to the initial, unshocked state. &' is the shock velocity in the material and &( is the particle velocity. The Hugoniot curve consists of the set of (PVT) points for which the Hugoniot expression: 1 (6) +, = ! − ! − + − 2 is zero. The procedure developed by Erpenbeck 57 was used to determine the Hugoniot curve for NM, which involves defining the initial, unshocked state of unreacted liquid NM and performing several separate canonical ensemble (NVT) simulations at appropriately chosen temperatures and densities. The initial, unshocked state of liquid NM was determined by equilibrating at - = 298 K and = 0 GPa. To compute the Hugoniot curve, a series of separate canonical ensemble (NVT) simulations are conducted over a range of temperatures and densities. Prior to performing the canonical simulations, the initial equilibrium densities at various states are determined by compressing a liquid sample at 400 K to pressures of 0.0, 0.1, 0.25, 0.5, 1.0, 2.5, 5.0, 7.5 and 10.0 GPa. After obtaining these equilibrated configurations, the densities are constrained and a series of NVT simulations are conducted for 50 ps at



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 32

temperatures of 300 K, 400 K, 500 K, 600 K, 700 K and 800 K. The equilibrium pressures are determined and the Hugoniot relation in Equation (6 is evaluated. The Hugoniot value, Hg, is plotted as a function of temperature, then extrapolated to zero to determine the corresponding Hugoniot temperature (-./, ). The Hugoniot temperature can then be used to determine the Hugoniot pressure (./, ) at the Hugoniot density. Figure 7 shows the pressure-volume and pressure-temperature Hugoniot curves and compares the various parameterizations to the available experimental data 58-60. Overall, the ReaxFF and MOES-optimized parameterizations are shown to match the P-V Hugoniot experimental data. For the P-T Hugoniot curve, the MOES-optimized models generally predict higher pressures than experiment. The ReaxFF model best captures the P-T Hugoniot curve over the temperature range. The ReaxFF-lg model is unable to capture either the P-V or the P-T Hugoniot behavior. Overall, ReaxFF and Model 0077 best match the experimental data for the P-V and P-T Hugoniot curves.

Figure 7: Pressure-volume (left) and pressure-temperature (right) shock Hugoniot equation of state for the various models. The experimental pressure-volume data of Marsh 58 and pressure-temperature data of Delpuech and Menil 59are denoted by black dots. The pressure-temperature Hugoniot equation-of-state of Lynse and Hardesty 60 is denoted by the black line. 3.5 Reaction chemistry assessment: The literature on the NM chemical kinetics experiments is rather extensive, providing information on the reaction pathways, reaction rate constants and decomposition products of the gas, liquid and solid phases (see Ref. 61 for a detailed review of the previous experiments). However, many questions remain regarding the details of the complex condensed phase reactions of NM that occur under high pressures and temperatures, conditions for which a number of theoretical studies have examined the initial thermal decomposition events in the liquid 33, 35 and solid phases 26, 36. Unfortunately, experimental evidence confirming the theoretical observations at these conditions is lacking, and experimental rate data are available only at lower temperatures. Conventional MD cannot approach the requisite time scales to simulate reactions at these temperatures, although development of promising theoretical methods to explore the reaction chemistry at longer time scales is ongoing 62. Due to the lack of definitive experimental evidence of reaction mechanisms at the time scales for which we can utilize atomistic simulations, we are unable to validate any predictions of reaction pathways and intermediate species reaction rates. Instead, we limit our discussion to NM decomposition events at high temperatures for the purpose of providing qualitative comparisons among the models discussed in this paper.


Page 15 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


A series of independent RMD simulations were conducted to investigate the NM decomposition rate at high temperatures predicted by the ReaxFF, ReaxFF-lg and MOES-optimized NM models. For each of the simulations performed at target temperatures ranging from 2500 K to 4000 K, we followed the approach of Rom et al. 35 In this approach, a crystalline sample of NM at a density of 1.40 g/cm3 (480 molecules) was equilibrated at 78 K, 0 GPa for 200 ps. The equilibrium system was then rapidly heated to the target temperature over a 100-fs time scale through NVT MD simulations with a 0.1 fs time step. Once the system reached the target temperature, an additional 500 ps NVT simulation was conducted. The ReaxFF bond-order and chemical species analysis tool provided within the LAMMPS distribution is used to identify whether bonds have been created or broken, i.e., whether chemical reactions have occurred. Pre-defined cutoff values are specified for each pair of atoms to compute the bond order. In this work, we used the predefined cutoff values defined in the supporting information of Rom et al. 35 An average bond order was computed every 100 time steps to determine the chemical species present in the system. Bond orders were sampled every 10 time steps, providing 10 samples every 100 time steps that were averaged. While the simulations simulate system heating at a consistent rate to the target temperatures, the inherent differences in the models prevent a rigorous, quantitative comparison of the resulting reaction chemistries. Because the intermolecular interactions differ among the models, the system pressures evolve at different rates as heating occurs, and converge to different pressures. The resulting pressures influence the decomposition of NM and subsequent formation of the product gases (see Table 4, for example, where the pressures during the RMD simulations at 3500 K are reported). Noting that these differences are present, we provide comparisons of the relative rates of decomposition, the initial decomposition reaction mechanisms, and the final products that are formed over long simulation times. Table 4: Decomposition and equilibrium pressures during the RMD simulation at 3500 K. Parameterization Pressure at 5% NM Decomposition (GPa)a Equilibrium Pressure (GPa)b Model 0498 5.47 ± 0.93 7.98 ± 0.02 Model 0524 6.68 ± 0.84 8.21 ± 0.02 Model 0077 7.17 ± 1.32 10.23 ± 0.04 ReaxFF-lg 5.61 ± 0.45 9.35 ± 0.01 ReaxFF 7.51 ± 0.70 9.35 ± 0.01 a The average pressure during NM decomposition is computed over a 100 fs time period, centered about the time at which 5% NM has decomposed. b The equilibrium pressure and standard deviation is computed though block averaging over 100 ps intervals over the simulation period of 200 ps to 500 ps. The amount of reacted NM is monitored at various target temperatures ranging from 2500 K to 4000 K (see Figure 8). Over the entire temperature range, Models 0498 and 0524 are observed to have the fastest decomposition rate, followed by ReaxFF and ReaxFF-lg with intermediate decomposition rates and Model 0077 with the slowest decomposition rate.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 32

Figure 8: NM decomposition at 2500 K, 3000 K, 3500 K and 4000 K, 0 GPa. The plots are normalized by the initial number of NM molecules and are offset from one another by 1.0 unit to illustrate the dependence on temperature. The reaction intermediates and products were monitored during the aforementioned simulations, and those resulting from the 3500 K simulation are shown in Figure 9. For all models, the formation of CH3 and NO2 occurs within a few picoseconds of the NM decomposition. The CH3 and NO2 species form in a 1:1 ratio, indicating that a scission of the C-N bond is the initial reaction step. Shortly thereafter, NO, CH2O, CH3NO, and CH4 species begin to form, followed by HNO formation. The models agree qualitatively in forming similar reaction intermediates, with differences in the rates and relative species concentrations. At longer simulation times, the reaction products converge to constant concentrations, as shown in Figure 9 (right frames). In general, all models result in the production of N2, H2, H2O, CO2 and NH3. The ReaxFF and ReaxFF-lg show consistent behavior, where the reaction rates and the relative concentrations are nearly identical. Water is observed to be the most abundant product species, followed by H2, N2, CO2 and NH3, respectively. Model 0524 also produces a large amount of H2O, but the relative proportion of H2 is less than that observed with the ReaxFF and ReaxFF-lg models. Model 0498 follows a similar trend, but the final products are produced at a slower rate than the ReaxFF and ReaxFF-lg models. Model 0077 differs from all other models, where large amounts of H2 and CO2 are produced, yet only small quantities of H2O are formed. While the temperature and pressure dependence of the NM decomposition mechanism may differ for each of the models, the early-time NM decomposition mechanisms appear to be qualitatively similar across the models, suggesting the pertinent chemical features are retained upon re-parameterization. The discrepancies, however, clearly highlight the need for


Page 17 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


information to refine the model to predict the correct overall chemistry, as well as to provide model validation.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 9: Reaction intermediates (left frames) and reaction products (right frames) after NM decomposition at 3500 K.


Page 18 of 32

Page 19 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4. Conclusions: In this work, we have provided the scientific community with the library of thousands of optimized ReaxFF-lg parameterizations for use in reactive molecular simulations obtained in an earlier study 14, and demonstrated a protocol for screening these for particular applications. In our case, we required a model of NM that describes material properties over different phases and a large range of temperatures and pressures, for the purpose of investigating the effect of shear on energetic materials under pressure 16. While we did not exhaustively search the entire library for the most transferable NM model, we demonstrated that a relatively simple screening approach could be utilized to obtain multiple models that perform equally as well as or better than existing ReaxFF and ReaxFF-lg parameterizations commonly used in the literature for certain applications. The models were compared to existing experimental structural, transport and thermodynamic properties of NM and were shown to better predict many material properties across a wide variety thermodynamic states than the existing models. We further compared the initial thermal decomposition events of NM against the original ReaxFF and ReaxFF-lg models to examine any effects that the parameter optimization process may have had on the reaction chemistry. In general, all models show an initial C-N bond scission and produce similar reaction intermediates and products, albeit with varying reaction rates and relative chemical compositions, indicating that the MOES-optimized models retain the pertinent reaction chemistry that is qualitatively consistent with the existing models. While quantitative differences between the models were observed when examining the temperature-dependence on the reaction rates and the relative species concentrations as the reactions progress, experimental results are not available to determine which, if any, model represents the actual system. This clearly highlights the need for advanced experimental and theoretical treatments to refine and validate reactive models. Ideally, there will always be the hope to develop models that are transferable to all chemical species with C-H-N-O chemistry as well as to all thermodynamic states. In practice, the models will generally perform well for the chemical systems, properties and conditions that lie within a region of phase space encompassing the training set data. ReaxFF and ReaxFF-lg have done a remarkable job of capturing chemistry of various chemical systems at ambient conditions, but may not always be appropriate for all chemical systems or for systems under extreme conditions. The modeler must always be cognizant of the strengths and limitations of a given model. Care should always be taken to faithfully assess whether the models are valid for the chemical systems or states to which they will be applied and whether reparameterization is required for use in applications for which the model has not been trained. In cases where the confidence in the model’s predictive ability is diminished, the approach outlined in this work addresses the need to tune the model parameters towards one’s specific needs and/or applications using an existing library of potentials that are suitable for C-H-N-O chemistry 14-15. If a suitable candidate is not obtained, the user can then reparameterize the force field with an augmented training set using published methods 14, 17. We do not make the claim that there is one model that will fit for all applications. Instead, we clearly demonstrate that these models are more suitable and transferable to high pressure states than the existing ReaxFF and ReaxFF-lg models. The rigorous comparison of the models against existing experimental measurements was conducted to highlight the strengths and limitations of each model, including commonly used ReaxFF and ReaxFF-lg models, in capturing material properties of different phases under various conditions. While no model is found to be remarkably superior for all properties tested (as expected), the assessment does provide insight into which properties or thermodynamic states the models are expected to be valid, and those which they are not. Ultimately, the choice of the best model will depend on the user’s application and the conditions under which the models will be applied.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 32

Acknowledgments: We wish to thank the High Performance Computing Modernization Program for providing the supercomputing resources used in this work. All the calculations presented were run at the DoD Supercomputing Resource Centers (DSRCs) at the U.S. Army, Air Force and Naval Research Laboratories. Supporting Information: The library of ReaxFF-lg potentials are supplied as a TXT file. Figure S1 and Tables S1-S18 described within the text are supplied in the PDF file. This material is available free of charge via the Internet at http://pubs.acs.org. References: 1. Farah, K.; Müller-Plathe, F.; Böhm, M. C., Classical Reactive Molecular Dynamics Implementations: State of the Art. ChemPhysChem 2012, 13 (5), 1127-1151. 2. Abell, G. C., Empirical Chemical Pseudopotential Theory of Molecular and Metallic Bonding. Phys. Rev. B 1985, 31 (10), 6184-6196. 3. Brenner, D. W., Empirical Potential for Hydrocarbons for Use in Simulating the Chemical Vapor Deposition of Diamond Films. Phys. Rev. B 1990, 42 (15), 9458-9471. 4. Brenner, D. W.; Shenderova, O. A.; Harrison, J. A.; Stuart, S. J.; Ni, B.; Sinnott, S. B., A SecondGeneration Reactive Empirical Bond Order (REBO) Potential Energy Expression for Hydrocarbons. J. Phys.: Condens. Matter 2002, 14 (4), 783. 5. Stuart, S. J.; Tutein, A. B.; Harrison, J. A., A Reactive Potential for Hydrocarbons with Intermolecular Interactions. J. Chem. Phys. 2000, 112 (14), 6472-6486. 6. Tersoff, J., Empirical Interatomic Potential for Carbon, with Applications to Amorphous Carbon. Phys. Rev. Lett. 1988, 61 (25), 2879-2882. 7. van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A., ReaxFF: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 105 (41), 9396-9409. 8. Plimpton, S., Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117 (1), 1-19. 9. Liu, L.; Liu, Y.; Zybin, S. V.; Sun, H.; Goddard, W. A., ReaxFF-lg: Correction of the ReaxFF Reactive Force Field for London Dispersion, with Applications to the Equations of State for Energetic Materials. J. Phys. Chem. A 2011, 115 (40), 11016-11022. 10. Senftle, T. P.; Hong, S.; Islam, M. M.; Kylasa, S. B.; Zheng, Y.; Shin, Y. K.; Junkermeier, C.; Engel-Herbert, R.; Janik, M. J.; Aktulga, H. M.; Verstraelen, T., et al., The ReaxFF Reactive Force-Field: Development, Applications and Future Directions. npj Comput. Mater. 2016, 2, 15011. 11. Jaramillo-Botero, A.; Naserifar, S.; Goddard, W. A., General Multiobjective Force Field Optimization Framework, with Application to Reactive Force Fields for Silicon Carbide. J. Chem. Theory Comput. 2014, 10 (4), 1426-1439. 12. Larsson, H. R.; van Duin, A. C. T.; Hartke, B., Global Optimization of Parameters in the Reactive Force Field ReaxFF for SiOH. J. Comput. Chem. 2013, 34 (25), 2178-2189. 13. Pahari, P.; Chaturvedi, S., Determination of Best-Fit Potential Parameters for a Reactive Force Field Using a Genetic Algorithm. J. Mol. Model. 2012, 18 (3), 1049-1061. 14. Larentzos, J. P.; Rice, B. M.; Byrd, E. F. C.; Weingarten, N. S.; Lill, J. V., Parameterizing Complex Reactive Force Fields Using Multiple Objective Evolutionary Strategies (MOES). Part 1: ReaxFF Models for Cyclotrimethylene Trinitramine (RDX) and 1,1-Diamino-2,2-Dinitroethene (FOX-7). J. Chem. Theory Comput. 2015, 11 (2), 381-391. 15. Rice, B. M.; Larentzos, J. P.; Byrd, E. F. C.; Weingarten, N. S., Parameterizing Complex Reactive Force Fields Using Multiple Objective Evolutionary Strategies (MOES): Part 2: Transferability of ReaxFF Models to C–H–N–O Energetic Materials. J. Chem. Theory Comput. 2015, 11 (2), 392-405.


Page 21 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


16. Steele, B. A.; Larentzos, J. P. Molecular Simulations of Shear-Induced Dynamics in Nitromethane; ARL Technical Note ARL-TN-0784; U.S. Army Research Laboratory: September, 2016. 17. Jones, W. M.; Giauque, W. F., The Entropy of Nitromethane. Heat Capacity of Solid and Liquid. Vapor Pressure, Heats of Fusion and Vaporization. J. Am. Chem. Soc. 1947, 69 (5), 983-987. 18. Trevino, S. F.; Prince, E.; Hubbard, C. R., Refinement of the Structure of Solid Nitromethane. J. Chem. Phys. 1980, 73 (6), 2996-3000. 19. Citroni, M.; Datchi, F.; Bini, R.; Di Vaira, M.; Pruzan, P.; Canny, B.; Schettino, V., Crystal Structure of Nitromethane up to the Reaction Threshold Pressure. J. Phys. Chem. B 2008, 112 (4), 10951103. 20. Yarger, F. L.; Olinger, B., Compression of Solid Nitromethane to 15 GPa at 298 K. J. Chem. Phys. 1986, 85 (3), 1534-1538. 21. Cromer, D. T.; Ryan, R. R.; Schiferl, D., The Structure of Nitromethane at Pressures of 0.3 to 6.0 GPa. J. Phys. Chem. 1985, 89 (11), 2315-2318. 22. Seminario, J. M.; Concha, M. C.; Politzer, P., A Density Functional/Molecular Dynamics Study of the Structure of Liquid Nitromethane. J. Chem. Phys. 1995, 102 (20), 8281-8282. 23. Tuckerman, M. E.; Klein, M. L., ab initio Molecular Dynamics Study of Solid Nitromethane. Chem. Phys. Lett. 1998, 283 (3–4), 147-151. 24. Reed, E. J.; Joannopoulos, J. D.; Fried, L. E., Electronic Excitations in Shocked Nitromethane. Phys. Rev. B 2000, 62 (24), 16500-16509. 25. Reed, E. J.; Riad Manaa, M.; Fried, L. E.; Glaesemann, K. R.; Joannopoulos, J. D., A Transient Semimetallic Layer in Detonating Nitromethane. Nat. Phys. 2008, 4 (1), 72-76. 26. Manaa, R. M.; Reed, E. J.; Fried, L. E.; Galli, G.; Gygi, F., Early Chemistry in Hot and Dense Nitromethane: Molecular Dynamics Simulations. J. Chem. Phys. 2004, 120 (21), 10146-10153. 27. Liu, H.; Zhao, J.; Wei, D.; Gong, Z., Structural and Vibrational Properties of Solid Nitromethane under High Pressure by Density Functional Theory. J. Chem. Phys. 2006, 124 (12), 124501. 28. Conroy, M. W.; Oleynik, I. I.; Zybin, S. V.; White, C. T., Density Functional Theory Calculations of Solid Nitromethane under Hydrostatic and Uniaxial Compressions with Empirical van der Waals Correction. J. Phys. Chem. A 2009, 113 (15), 3610-3614. 29. Xu, K.; Wei, D.-Q.; Chen, X.-R.; Ji, G.-F., Thermal Decomposition of Solid Phase Nitromethane under Various Heating Rates and Target Temperatures Based on ab initio Molecular Dynamics Simulations. J. Mol. Model. 2014, 20 (10), 1-10. 30. Zhang, J.-D.; Kang, L.-H.; Cheng, X.-L., Theoretical Study of the Reaction Mechanism of CH3NO2 with NO2, NO and CO: The Bimolecular Reactions That Cannot Be Ignored. J. Mol. Model. 2015, 21 (1), 1-6. 31. Citroni, M.; Bini, R.; Pagliai, M.; Cardini, G.; Schettino, V., Nitromethane Decomposition under High Static Pressure. J. Phys. Chem. B 2010, 114 (29), 9420-9428. 32. van Duin, A. C. T.; Zybin, S. V.; Chenoweth, K.; Zhang, L.; Han, S. P.; Strachan, A.; Goddard, W. A., Reactive Force Fields Based on Quantum Mechanics for Applications to Materials at Extreme Conditions. AIP Conf. Proc. 2006, 845 (1), 581-584. 33. Hervouët, A.; Desbiens, N.; Bourasseau, E.; Maillet, J.-B., Microscopic Approaches to Liquid Nitromethane Detonation Properties. J. Phys. Chem. B 2008, 112 (16), 5070-5078. 34. Guo, F.; Zhang, H.; Cheng, X., Molecular Dynamic Simulations of Solid Nitromethane under High Pressures. J. Theor. Comput. Chem. 2010, 09 (01), 315-325. 35. Rom, N.; Zybin, S. V.; van Duin, A. C. T.; Goddard, W. A.; Zeiri, Y.; Katz, G.; Kosloff, R., Density-Dependent Liquid Nitromethane Decomposition: Molecular Dynamics Simulations Based on ReaxFF. J. Phys. Chem. A 2011, 115 (36), 10181-10202. 36. Han, S.-p.; van Duin, A. C. T.; Goddard, W. A.; Strachan, A., Thermal Decomposition of Condensed-Phase Nitromethane from Molecular Dynamics from ReaxFF Reactive Dynamics. J. Phys. Chem. B 2011, 115 (20), 6534-6540.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

37. Guo, F.; Cheng, X.-l.; Zhang, H., Reactive Molecular Dynamics Simulation of Solid Nitromethane Impact on (010) Surfaces Induced and Nonimpact Thermal Decomposition. J. Phys. Chem. A 2012, 116 (14), 3514-3520. 38. Fileti, E. E.; Chaban, V. V.; Prezhdo, O. V., Exploding Nitromethane in Silico, in Real Time. J. Phys. Chem. Lett. 2014, 5 (19), 3415-3420. 39. Guo, F.; Zhang, H.; Hu, H.-Q.; Cheng, X.-L.; Zhang, L.-Y., Hugoniot Curve Calculation of Nitromethane Decomposition Mixtures: A Reactive Force Field Molecular Dynamics Approach. Chin. Phys. B 2015, 24 (11), 118201. 40. Wood, M. A.; Strachan, A., Nonequilibrium Reaction Kinetics in Molecular Solids. J. Phys. Chem. C 2016, 120 (1), 542-552. 41. Nielson, K. D.; van Duin, A. C. T.; Oxgaard, J.; Deng, W.-Q.; Goddard, W. A., Development of the ReaxFF Reactive Force Field for Describing Transition Metal Catalyzed Reactions, with Application to the Initial Stages of the Catalytic Formation of Carbon Nanotubes. J. Phys. Chem. A 2005, 109 (3), 493-499. 42. Strachan, A.; van Duin, A. C. T.; Chakraborty, D.; Dasgupta, S.; Goddard, W. A., Shock Waves in High-Energy Materials: The Initial Chemical Events in Nitramine RDX. Phys. Rev. Lett. 2003, 91 (9), 098301. 43. Strachan, A.; Kober, E. M.; van Duin, A. C. T.; Oxgaard, J.; Goddard, W. A., Thermal Decomposition of RDX from Reactive Molecular Dynamics. J. Chem. Phys. 2005, 122 (5), 054502. 44. Zhang, L.; Zybin, S. V.; van Duin, A. C. T.; Dasgupta, S.; Goddard, W. A.; Kober, E. M., Carbon Cluster Formation During Thermal Decomposition of Octahydro-1,3,5,7-Tetranitro-1,3,5,7-Tetrazocine and 1,3,5-Triamino-2,4,6-Trinitrobenzene High Explosives from ReaxFF Reactive Molecular Dynamics Simulations. J. Phys. Chem. A 2009, 113 (40), 10619-10640. 45. Flyvbjerg, H.; Petersen, H. G., Error Estimates on Averages of Correlated Data. J. Chem. Phys. 1989, 91 (1), 461-466. 46. Frenkel, D.; Smit, B., Understanding Molecular Simulation : From Algorithms to Applications. 2nd ed.; Academic Press: San Diego, Calif, 2002. 47. Murnaghan, F. D., Finite Deformation of an Elastic Solid. Wiley: New York, 1951. 48. Agrawal, P. M.; Rice, B. M.; Thompson, D. L., Molecular Dynamics Study of the Melting of Nitromethane. J. Chem. Phys. 2003, 119 (18), 9617-9627. 49. Trevino, S. F.; Rymes, W. H., A Study of Methyl Reorientation in Solid Nitromethane by Neutron Scattering. J. Chem. Phys. 1980, 73 (6), 3001-3006. 50. Rabinovich, I. B., Influence of Isotopy on the Physicochemical Properties of Liquids. American ed.; Consultants Bureau: New York,, 1970. 51. Price, W. S.; Ide, H.; Arata, Y., Translational and Rotational Motion of Isolated Water Molecules in Nitromethane Studied Using 17O NMR. J. Chem. Phys. 2000, 113 (9), 3686-3689. 52. Holz, M.; Mao, X. a.; Seiferling, D.; Sacco, A., Experimental Study of Dynamic Isotope Effects in Molecular Liquids: Detection of Translation‐Rotation Coupling. J. Chem. Phys. 1996, 104 (2), 669679. 53. Morris, J. R.; Song, X., The Melting Lines of Model Systems Calculated from Coexistence Simulations. J. Chem. Phys. 2002, 116 (21), 9352-9358. 54. Morris, J. R.; Wang, C. Z.; Ho, K. M.; Chan, C. T., Melting Line of Aluminum from Simulations of Coexisting Phases. Phys. Rev. B 1994, 49 (5), 3109-3115. 55. Eike, D. M.; Brennecke, J. F.; Maginn, E. J., Toward a Robust and General Molecular Simulation Method for Computing Solid-Liquid Coexistence. J. Chem. Phys. 2005, 122 (1), 014115. 56. Eike, D. M.; Maginn, E. J., Atomistic Simulation of Solid-Liquid Coexistence for Molecular Systems: Application to Triazole and Benzene. J. Chem. Phys. 2006, 124 (16), 164503. 57. Erpenbeck, J. J., Molecular Dynamics of Detonation. I. Equation of State and Hugoniot Curve for a Simple Reactive Fluid. Phys. Rev. A 1992, 46 (10), 6406-6416. 58. Marsh, S. P., LASL Shock Hugoniot Data. Univ of California Press: 1980; Vol. 5.


Page 22 of 32

Page 23 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


59. Delpuech, A.; Mentil, A.; Pouligny, B., Raman Scattering Temperature Measurement Behind a Shock Wave a Study of Solid Explosives. In Shock Waves in Condensed Matter; Springer US: Boston, MA, 1986; pp 877-882. 60. Lysne, P. C.; Hardesty, D. R., Fundamental Equation of State of Liquid Nitromethane to 100 kbar. J. Chem. Phys. 1973, 59 (12), 6512-6523. 61. Blais, N. C.; Engelke, R.; Sheffield, S. A., Mass Spectroscopic Study of the Chemical Reaction Zone in Detonating Liquid Nitromethane. J. Phys. Chem. A 1997, 101 (44), 8285-8295. 62. Cheng, T.; Jaramillo-Botero, A.; Goddard, W. A.; Sun, H., Adaptive Accelerated ReaxFF Reactive Dynamics with Validation from Simulating Hydrogen Combustion. J. Am. Chem. Soc. 2014, 136 (26), 9434-9442. TOC Graphic:



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1: Representation of the NM crystal unit cell determined by Trevino et al. 18 at 4.2 K, 1 atm. Atom labels are consistent with those used to describe H-C-N-O dihedral angles hereafter. Figure 1 82x50mm (300 x 300 DPI)


Page 24 of 32

Page 25 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2: Lattice constant dependence on pressure (left frame) and % error (right frame). Figure 2 177x114mm (300 x 300 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3: Pressure and temperature dependence of the rms deviations of predicted atomic positions from experiment due to molecular deformation (top frames); the rms deviations of predicted atomic positions from experiment due to relative orientation in the unit cell (middle frames); the absolute error of the positions of the center-of-mass fractional positions, averaged over molecules and dimensions (bottom frames). Data for ReaxFF at 228K are not included (see text). Figure 3 177x128mm (300 x 300 DPI)


Page 26 of 32

Page 27 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4: The temperature-dependent (top frames) and pressure-dependent (bottom frames) H-C-N-O dihedral angle distributions. Models 0498 and 0526 are omitted for clarity, but are similar to Model 0077. The dashed, vertical lines are provided to illustrate peak shifts with pressure. Figure 4 177x78mm (300 x 300 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 5: Liquid-phase density dependence of NM on temperature of the various models compared with experiment. Figure 5 82x64mm (300 x 300 DPI)


Page 28 of 32

Page 29 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6: Semi-log plot of liquid-phase NM self-diffusivity as a function of temperature compared with experiment 51-52. 82x64mm (300 x 300 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 7: Pressure-volume (left) and pressure-temperature (right) shock Hugoniot equation of state for the various models. The experimental pressure-volume data of Marsh 56 and pressure-temperature data of Delpuech and Menil 57 are denoted by black dots. The pressure-temperature Hugoniot equation-of-state of Lynse and Hardesty 58 is denoted by the black line. Figure 7 177x71mm (300 x 300 DPI)


Page 30 of 32

Page 31 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 8: NM decomposition at 2500 K, 3000 K, 3500 K and 4000 K, 0 GPa. The plots are normalized by the initial number of NM molecules and are offset from one another by 1.0 unit to illustrate the dependence on temperature. Figure 8 82x123mm (300 x 300 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 9: Reaction intermediates (left frames) and reaction products (right frames) after NM decomposition at 3500 K. Figure 9 139x241mm (300 x 300 DPI)


Page 32 of 32

Transferable Reactive Force Fields: Extensions of ReaxFF-lg to

Recommend Documents