The Adaptive Intermolecular Reactive Empirical Bond Order potential (AIREBO) for hydrocarbons has been widely used to study dynamic bonding processes under ambient conditions. However, its intermolecular interactions are modeled by a Lennard-Jones (LJ) potential whose unphysically divergent power-law repulsion causes AIREBO to fail when applied to systems at high pressure. We present a modified potential, AIREBO-M, where we have replaced the singular Lennard-Jones potential with a Morse potential. We optimize the new functional form to improve intermolecular steric repulsions, while preserving the ambient thermodynamics of the original potentials as much as possible. The potential is fit to experimental measurements of the layer spacing of graphite up to 14 GPa and first principles calculations of steric interactions between small alkanes. To validate AIREBO-M’s accuracy and transferability, we apply it to a graphite bilayer and orthorhombic polyethylene. AIREBO-M gives bilayer compression consistent with quantum calculations, and it accurately reproduces the quasistatic and shock compression of orthorhombic polyethlyene up to at least 40 GPa.

## I. INTRODUCTION

Extreme pressures can activate many mechanical and chemical processes in solids that are usually dormant at ambi ent conditions. A wide array of complex phenomena is common at 10 to 100 GPa including: defect nucleation in metals, structural phase transitions in molecular crystals, and bond rehybridization and wear in covalent solids.^{1–3} Such dynamic complexity makes high pressure material behavior a rich av enue for research. However, the extreme conditions necessary are difficult to maintain experimentally, making robust numerical tools desirable.

Quantum mechanics (QM) based numerical methods, such as Density Functional Theory (DFT), have been successfully applied to explore the dynamic chemistry and energetics of many materials at high pressure and temperatures.^{4–6} However, the computational expense and purely covalent form of such methods prevent them from describing the structural and mechanical properties of many extended molecular materials, such as graphitic and hydrocarbon solids, which are a rapidly growing focus of industrial applications. Dispersion interactions can be included in DFT calculations,^{7–10} for example, by adding semi-empirical potentials,^{10–13} but this increases the computational expense. More efficient methods are needed to study dynamic problems in large systems.

An attractive alternative to QM calculations is the family of reactive bond order potentials (BOPs). Traditional interaction potentials define and maintain a fixed topology of chemical bonds. The development of BOPs, such as the Adaptive Intermolecular Reactive Empirical Bond Order potential (AIREBO)^{14} and the Reactive Force Field (REAXFF),^{15} has allowed modelers to explore dynamic bonding processes. The complexity of the many-particle interactions in BOPs makes parameterization difficult. As a result, current BOPs quickly lose fidelity outside a limited range of ambient pressures and temperatures, corresponding to near equilibrium states of their molecular and chemical bonds. In order to accurately describe far from equilibrium pressures and dynamics, existing BOPs must be modified to improve their description of highly compressed molecular configurations. For AIREBO and REAXFF, these limitations are well illustrated by the recent work of Mattsson *et al.* on polyethylene (PE),^{4} which concluded that both models greatly over-predict the stiffness of PE compared to experimental data at pressures above a few GPa. For AIREBO, the excess stiffness results from the use of Lennard-Jones (LJ) interactions between molecules that diverge rapidly at small separations.

In this work, we replace the Lennard-Jones interactions in AIREBO by a nondivergent Morse potential. The Morse potential has an additional parameter that allows the repulsive region to be softened with little change to the minimum and attractive regions. This allows us to extend the validity of the potential to pressures of about 40 GPa without affecting the behavior at ambient pressure that has been successfully captured by the original AIREBO.^{14} The Morse potential for carbon-carbon interactions is fit to x-ray experiments on graphite and tested against theoretical calculations for bilayers. Similar x-ray data for hydrocarbons do not separately constrain the form of C–H and H–H interactions. This is accomplished by fitting to post-Hartree Fock quantum calculations for a range of configurations. The resulting parameters are then validated against x-ray and shock Hugoniot data for crystalline PE.

In what follows, we discuss the underlying theory and the parameterization and validation of the new model. Section II discusses the theoretical underpinning and algorithmic details of the simulations performed. Section III explains the rationale for our parameterization protocol and reports the resulting Morse parameters. Section IV presents results from the new model for the graphene bilayer and the isothermal compression and shock properties of orthorhombic PE.

## II. POTENTIALS AND SIMULATION METHODS

### A. Potentials

The Second Generation Reactive Empirical Bond Order Potential (REBO2) of Brenner *et al.*^{16} is a reactive potential that dynamically describes covalent hydrocarbon chemistry with classical many-body potential energy functions. Brenner *et al.*’s potential, originally intended to study network carbon solids like diamond, does not contain any energy terms for intermolecular interactions. Thus, REBO2 cannot describe the structure of liquids or molecular solids with any accuracy.

The AIREBO potential of Stuart *et al.*^{14} adds nonbonded/intermolecular interactions for hydrocarbons to the purely covalent REBO2 via a set of Lennard-Jones potentials for the three interaction types: C–C, C–H, and H–H. The literature contains many applications of AIREBO to a wide variety of carbon and hydrocarbon allotropes—ranging from gas phase hydrocarbon combustion to phase separation in organic fluids to exotic carbon-nanotubes.^{17–19} The AIREBO potential is openly available in Sandia’s LAMMPS software, which is used for all calculations presented in this paper. A thorough and concise description of the model can be found in the appendix of Stuart’s original article.^{14}

AIREBO’s nonbonded interactions are modeled by a smoothly truncated 12-6 LJ potential

where the *ij* indices indicate the chemical species (C or H) of the two interacting atoms, *ϵ* defines a depth for the interaction potential well, and 2^{1/6}*σ* defines the location of the minimum energy. The potential is smoothly set to zero by a third-order spline at a long-range cutoff value of 3*σ*. For separations less than 2^{1/6}*σ*, similar splines reduce the strength of intermolecular repulsions if it is favorable for a pair of atoms to bond chemically.^{14} When covalent bonding is unfavorable, the two atoms interact with the full LJ.

The use of a repulsive *r*^{−12} power law in the LJ potential was motivated more by computational convenience than physical reasoning. It causes the potential to rise rapidly as molecular separation decreases past the minimum.^{4} The anomalously high repulsive forces cause AIREBO to fail to describe the details of many high pressure systems, even while accurately describing the same systems in ambient conditions. To improve the model’s accuracy at high pressures, we have replaced its Lennard-Jones potentials with a Morse potential

where *ϵ* and *r ^{eq}* define the depth and location of the minimum energy, and the new parameter

*α*modifies the curvature of the potential energy at its minimum separation. The Morse potential is truncated with a third-order spline at the same cutoffs used in AIREBO, and the splines for adaptive repulsion are implemented by the same technique as AIREBO. Specifically, the outer and inner cut-offs for these splines are set, respectively, to the location of the energy minimum and the location where the energy crosses zero.

^{14}

The Morse potential’s parameters can be adjusted to match the equilibrium properties of hydrocarbons that are captured by the LJ potentials in AIREBO. One approach is to fit the location, depth, and second derivative of the potential energy minimum for each interaction type. This gives a mapping from (*ϵ*, *σ*) to (*ϵ*, *r ^{eq}*,

*α*)

The resulting unoptimized parameters are listed in Table I. However, the presence of a third parameter in the Morse potential allows more freedom in fitting a wide range of separations than in the LJ potential. We use Eq. (3) and AIREBO’s LJ parameters to obtain initial values of Morse parameters and then adjust *α* to fit repulsive interactions without affecting the equilibrium spacing and energy. In this way, we improve the high pressure response while retaining agreement with the region of the AIREBO potential that determines equilibrium thermodynamic properties at ambient pressures.

### B. Simulation methods

In the model parameterization and validation that follows, we perform two types of simulation protocol. When making comparisons to experimental measurements for pressure-volume (PV) curves of graphite and PE or the shock Hugoniot equation of state for PE, we utilize NPT and NVT ensemble molecular dynamics (MD) with periodic systems of the corresponding bulk material. When making comparison to first principle calculations of the microscopic interactions between ethane and pentane dimers or a graphene bilayer, we simply perform static energy calculations on pregenerated rigid configurations of the molecules. We next briefly describe the generation of these structures and the numerical details of the simulation protocols. Snapshots of simulation cells are shown in the figures presenting the associated properties.

#### 1. Structure generation

The bulk graphite configuration consists of ten graphene layers, each layer containing 416 carbon atoms, shown in Fig. 1. The layers are approximately square and arranged in an alternating AB stacking, with a period *c* equal to twice the experimental interlayer spacing of 3.35 angstroms. We create the rigid graphene bilayer by taking two layers of the bulk graphite system and making the simulation box periodic for the in-plane directions only. The trajectory is generated by varying the layer separation.

The polyethylene configuration consists of 5 by 7 by 12 unit cells of the orthorhombic (Pnam) PE crystal.^{20} Thus, the simulation cell has 70 PE chains with 12 monomers along the chain per period of the simulation box. Configurations undergo geometry optimization with the Fast Inertial Relaxation Engine (FIRE) minimizer,^{21} before undergoing further MD treatment.

We use ethane and pentane molecules to generate dimer configurations exploring coordination registries relevant to bulk PE, while still being small enough molecules for use with accurate quantum mechanical techniques. We create six dimer trajectories, four for pentane-pentane dimers and two for ethane-ethane, illustrated in Fig. 2. Each trajectory is generated by rigidly rotating and translating a second ethane or pentane molecule with respect to the first, varying their separation along a fixed axis. We generate initial geometries for the alkane molecules with geometry optimization at the DFT level (M06/aug-cc-pvTZ).^{22}

#### 2. Molecular dynamics methods

Aside from the specific atomic configurations used, the protocols for computing the equation of state for both graphite and PE systems are identical. All MD simulations use the LAMMPS simulation package.^{23} The velocity-Verlet integra tor with a time-step of 1 fs is used to solve the equations of motion within a periodic simulation box. A Nosé-Hoover thermostat with time constant 0.25 ps is used for constant temperature simulations. To maintain a specified hydrostatic pressure, the three periods of the simulation cell are adjusted independently by a Nose-Hoover barostat with time constant 1.0 ps.^{24}

For pressure-volume curve generation, systems equilibrate to their initial isotherm temperature (300 or 553 K) and pressure (0 GPa) during 20-30 ps of simulation in the NPT ensemble. Next, the system pressure is increased by 2 GPa over 2 ps, followed by equilibration at fixed pressure for 4 ps. Such 6 ps pressure steps continue up to a pressure of 40 GPa, with PV points generated from the average of the last 2 ps of equilibration. Tests with longer equilibration times gave equivalent results. The Hugoniot equation of state for PE is produced from such MD generated isotherms via the quasi-static method of Erpenbeck. The method is described in detail in Ref. 25 and has been shown to be insensitive to the system size and simulation parameters.^{26,27}

#### 3. Quantum calculation methods

Dispersion corrected DFT methods are the most commonly used quantum mechanical techniques for studying molecular crystals and large clusters.^{28} However, despite their successes in predicting equilibrium structures for models of crystalline PE, the same methods predict bulk moduli that are significantly too large.^{28,29} DFT methods without dispersion corrections predict crystal structures more accurately as pressure increases,^{20,29} but at pressures of 2-3 GPa, the error in unit cell vectors may be as large as 0.2 Å.^{20} These findings suggest that more accurate *ab initio* methods may be necessary to study the steric interactions of hydrocarbon solids at extreme pressures.

*Ab initio* calculations on small models of *n*-alkane dimers^{30} have revealed that the second-order Møller-Plesset (MP2) values of equilibrium interaction energies obtained by using the correlation consistent basis set of triple zeta quality (aug-cc-pvTZ)^{31} are within 0.01 eV, or better, of the MP2 energies at the basis set limit. Such MP2 energies are also very close to results obtained from the more accurate coupled-cluster with singles, doubles, and approximated triples (CCSD(T)) calculation method.^{30,32} The CCSD(T) method is prohibitively expensive for some of our non-symmetrical pentane dimmers and therefore, in this work we use the MP2/aug-cc-pvTZ approach. Interaction energies are further corrected for basis set superposition errors (BSSE)^{33} by using the counterpoise (CP) method.^{34}

Using the MP2/aug-cc-pvTZ/CP model, we compute dimer interaction energy curves from each frame of the ethane and pentane dimer trajectories, which range in separation from the limit of molecular dissociation to large steric repulsions at short distances. As expected, we find the BSSE corrections to be critically important and they account for up to 50% of the total binding energy at equilibrium distances. While accurate dispersion forces play a diminishing role at high pressure, where steric interactions dominate, they are essential for accurately defining the equilibrium separations, which distinguish between repulsion and attraction dominated regions. Our fitting strategy is based on the fact that the unmodified AIREBO captures the ambient thermodynamic and structural energetics of many alkanes. Thus, we use our MP2 calculations to improve the steric rise in energy with small separation, while deferring to the unmodified AIREBO for the details of interaction energy wells.

## III. PARAMETERIZATION

In the original AIREBO paper, Lennard-Jones potentials were fit in two stages. C–C interactions were first fit to the properties of graphite, then C–H and H–H interactions were fit simultaneously to the x-ray structure and thermodynamics of small alkane liquids. To fit our Morse parameters, we follow a similar process. First, Morse C–C interactions are tuned to accurate high pressure experiments on graphite by Hanfland *et al.*^{35} This approach is difficult to apply for hydrocarbons because structural measurements do not isolate the separate contributions of C–H and H–H intermolecular bonds. Thus, we use MP2 calculations for a range of ethane and pentane dimer geometries to determine the C–H and H–H Morse potentials.

### A. Fitting C–C interactions

For typical hydrocarbon systems, carbon atoms rarely come within steric interaction distances of one another due to the steric repulsion of surrounding hydrogens. This makes it difficult to use hydrocarbon data to determine C–C interactions. Using pure carbon allotropes, such as graphene, carbon nanotubes, or graphite, isolates the effect of C–C interactions. These interactions directly determine the interlayer spacing of graphite, which has been measured with high pressure x-ray scattering by Hanfland *et al.*^{35} Over the studied range of pressures, the in-plane lattice constant is nearly constant. Since it is mostly determined by covalent interactions, AIREBO and AIREBO-M give consistent in-plane results.

Hanfland *et al.*^{35} used x-ray scattering to track the equilibrium layer separation *c*/2 of a graphite sample at room temperature from 0 to 14 GPa. To reproduce this experiment, we perform finite temperature MD in the NPT ensemble at pressures across this range, for a periodic bulk graphite system. We measure the layer spacing of the system as a function of pressure and modify the C–C Morse parameters until our system reproduces the equilibrium spacing and compression profile seen by Hanfland *et al.*

We found that the parameters in Table I did not give the measured lattice constants of graphite at ambient pressure. Thus, in addition to optimizing the new degree of freedom, *α*, we also slightly reduced the length scale of the interaction, *r ^{eq}*. We believe this is necessary because AIREBO’s LJ parameters for C–C were originally fit to graphite using an analytic expression that excluded the effects of AIREBO’s splines for adaptive repulsion. Plotting the resulting curves in Fig. 1, we see the original AIREBO overestimates both graphite’s equilibrium layer spacing at 0 GPa as well as its mechanical response under compression for T = 300 K. Note that this difference will have little effect on simulations of hydrocarbons.

The new Morse potential parameters for C–C may be found in the first row of Table II. The model with the new parameters, which we call AIREBO-M, produces equilibrium and response properties for graphite layer compression in excellent agreement with the experiment of Hanfland *et al.* (Fig. 1).

### B. Fitting C–H and H–H interactions

Compared to pure carbon, we find hydrogen’s dispersion forces are substantially weaker, and steric repulsion dominates its intermolecular interactions. Such steric repulsions are due to Pauli exclusion effects, which quantum mechanical methods tend to capture very well. Additionally, C–H and H–H interactions tend to be sensitive to the geometry and relative orientations of participating hydrocarbon molecules. Such details are beyond the scope of most experimental techniques, which typically only provide average lattice constants. In contrast, first principle calculations provide independent information about the energies of different geometries. To this end, we have used MP2,^{31} a post Hartree-Fock quantum chemistry method, to construct a training set of interaction energies for rigidly translated ethane and pentane dimers, illustrated in Fig. 2.

The MP2 calculations compute the total energy as a function of separation for each dimer, and since the individual molecules are held fixed, changes in the energy correspond to intermolecular interactions only. Analogous energy vs. separation curves are calculated with AIREBO-M for different parameters to determine the best fit. Since we are primarily concerned with improving the interaction’s high pressure response, and we wish to preserve the accuracy of AIREBO in ambient conditions, we limit the range of separations used for calculating the fit to those where steric repulsions dominate. We define the steric region of each trajectory as any separation where the MP2 energy is greater than 0.1 eV (with 0 eV at infinite separation). Fig. 2 shows data to 0.5 eV, but fits extended up to energies where covalent bonding changed, typically 1.5 eV.

To minimize changes from AIREBO, we optimize by scaling the *α*’s by the same factor for C–H and H–H—i.e., $ \alpha ij =\kappa \alpha ij 0 $ where the $ \alpha ij 0 $ are defined in Table I. The least mean squared fit gives *κ* = 0.877. The results are in excellent agreement with MP2 in four of the panels in Fig. 2 but give a higher binding energy and shorter distance in panels (b) and (d). Fits with other *α* and *r ^{eq}* do not improve the agreement in these panels without producing comparable deviations in the other panels. We expect that the difficulty in simultaneously fitting all configurations may reflect directional bonding that is not captured in the LJ or Morse potential. Given that varying

*r*and

^{eq}*ϵ*does not improve the global fit and that the AIREBO values have been validated in ambient conditions, we leave them unchanged. The Morse potential parameters for C–H and H–H may be found in the second and third rows of Table II.

Figure 3 illustrates the effect of rescaling *α* in the manner just described. It compares the LJ potential to the Morse potential with *α* chosen to match the second derivative at the minimum, via Eq. (3), as well as with *α* reduced by the factor *κ* = 0.877 used in Fig. 2. Note that reducing *κ* from 1.0 to 0.877 produces a much better fit of the entire attractive portion of the LJ potential. Although *κ* was chosen to match the repulsive region where the LJ potential rises too rapidly, it produces a Morse potential that is very close to providing the best fit to the attractive region of the LJ potential from *r*/*σ* = 2^{1/6} to 3. Over this range, the rms deviation from the LJ potential is only 0.0004*ϵ*. Thus, this choice of *κ* satisfies both of our objectives for the revised potential.

## IV. VALIDATION

In the parameterization protocol, we applied a mixed training set of experimental data and theoretical quantum chemistry calculations chosen to distinguish most cleanly and accurately the interactions for different elements. To validate the fits, we use information from complementary techniques. Our fits of C–C interactions to finite temperature x-ray experiments are compared to theoretical DFT calculations of the energy of graphene bilayers. Similarly, to validate our fits of C–H and H–H interactions to quantum chemistry calculations of the energy of small alkanes, we compare to finite temperature experiments on cyrstalline PE at high pressure and during shocks.

### A. Compressing a graphene bilayer

Understanding the planar interactions of two graphene sheets is important for the development of exfoliation-based processing and nanoscale fabrication in general. The literature contains many high quality studies on the layer interactions of graphene, including a recent dispersion-corrected DFT based analysis by Reguzzoni *et al.*^{36} Their work calculates the interaction energy as a function of the separation of two flat graphene sheets, using the DFT-D method of Grimme.^{10} Grimme’s method for dispersion corrected DFT is widely used^{10–13} and we consider the DFT-D calculations of Reguzzoni *et al.* to be a well-characterized theoretical reference for testing the compression mode our C–C interaction potential is fit to.

In the analysis of Reguzzoni *et al.*, they compute an energy versus separation curve for two rigid and periodic graphene sheets. To compare to their data, we generate an analogous configuration and repeat their calculation with our potential. Fig. 4 compares the interaction energy of the bilayer planes for AIREBO-M (red circles), AIREBO (blue squares), and DFT-D (black diamonds). We see that AIREBO-M produces a rise in steric energy under compression in very good agreement with the quantum calculation. AIREBO-M also reproduces the location and curvature of the DFT-D minimum energy, even though we did not intentionally seek to reproduce these features. While this agreement may be due to similar training systems for AIREBO-M and the DFT-D dispersion corrections, we consider AIREBO-M’s excellent reproduction of the bilayer mechanics calculated with the more theoretically sophisticated and expensive DFT-D to be a strong validation of our experimental fits for C–C interactions. Note that the repulsive interactions at the smallest separations correspond to significantly larger pressures (∼18 GPa) than achieved in the experiments of Hanfland *et al.*

The three methods show different behaviors at large distances and thus have different binding energies. Both theoretical and experimental values of the cohesion of a graphene bilayer show variations that are comparable to the difference between the three values.^{37} For this reason, we are satisfied that our model produces a cohesive energy comfortably within the range of current experimental values.

### B. Isothermal and Shock compression of crystalline polyethylene

Highly aligned polyethylene has become a focus of research for wide-scale industrial application due to its high tensile strength and relatively low mass density.^{38,39} Ultra-High Molecular Weight Polyethylene (UHMWPE), in the form of fibers and composite fabrics, is already used to create cables, helmets, body armor, and other products, where a high strength to weight ratio is critical. Building an understanding of PE’s structural behavior at high pressures is critical to improve its application; however, there are relatively few detailed experiments available in the literature.

One such experiment is the thorough x-ray analysis by Fontana *et al.* of PE brought to high pressures in a diamond anvil cell.^{20} Their work explores the anisotropic deformation of the orthorhombic PE crystal from 2 to 40 GPa at 553 K, as well as its eventual structural transition to higher-pressure monoclinic phases around 16 GPa. It is an ideal experiment to validate the AIREBO-M’s theoretically fit C–H and H–H interactions, which dominate the mechanical response of the molecular crystal. To compare to the experiments of Fontana *et al.*, we again perform finite temperature MD in the NPT ensemble, analogous to that used for bulk graphite, to generate a series of states along the T = 553 K PV isotherm. From these states, we determine the deformation of the PE unit cell for all three lattice directions.

In Fig. 5(a), we plot the change in the lattice vectors of orthorhombic PE with increasing pressure. The experimental data show that deformation primarily occurs in the lateral (**a** and **b**) lattice directions, where the crystal is stabilized by weaker intermolecular interactions. The axial (**c**) direction, dominated by covalent carbon-carbon bonds, responds negligibly to the pressures explored. Based on their data, Fontana *et al.* predict that the orthorhombic PE phase becomes metastable at pressures higher than 16 GPa, in favor of a high-pressure monoclinic phase. Beyond 25 GPa, the x-ray experiments were no longer able to distinguish the orthorhombic signature. After losing the orthorhombic signature, Fontana *et al.* continued tracking higher pressure crystal phases, and used them to generate PV data up to 40 GPa, reproduced in Fig. 5(b).

AIREBO-M (red circles) shows excellent agreement with the Fontana experiment (black diamonds) for deformation in all three lattice directions. It accurately distinguishes between the structural differences in the lateral lattice directions and reproduces the shape of their mechanical response. AIREBO-M also predicts a structural phase transition to a new orthorhombic phase, seen as a jump in the lattice parameters, between 18 and 20 GPa. This new phase corresponds to a flattening of the “herringbone” pattern of the orthorhombic crystal as the constituent chains come in closer contact. AIREBO-M predicts this structural transition soon after the 16 GPa point, where Fontana *et al.* predict the low pressure orthorhombic phase becomes metastable and begins to coexist with a high pressure monoclinic crystal phase. Whether this new phase is compatible with the x-ray data of Fontana *et al.*, or may be another metastable, transitionary phase to the reported monoclinic phase is a subject of current research. In contrast, AIREBO (blue squares) over-predicts the stiffness for deformations in both lateral directions. This is caused by the diverging repulsive interactions of the LJ potentials. Consequently, it washes out much of the underlying molecular structure and does not show any structural transitions.

We have performed similar calculations with the REAXFF parameterization of Chenoweth *et al.*^{15} (green triangles). We use this parameterization because of its recent use in studies of similar systems by Mattsson *et al.*,^{4} and because it is the version that is distributed in current and previous versions of LAMMPS. Our REAXFF results agree with the experimental data at the lowest pressures (2 GPa), but quickly lose fidelity as pressure increases. REAXFF also shows a structural transition to a new orthorhombic phase, similar to that of AIREBO-M, but it occurs near 5 GPa, far below the experimentally proposed metastability point at 16 GPa. Additionally, the REAXFF shows extension in the axial, **c**, direction with increasing compressive load. This counterintuitive expansion under compression is not seen in experiments and calls into question the applicability of this parameterization for hydrocarbon solids.

In Fig. 5(b), we see that AIREBO-M reproduces the entire experimental PV curve of Fontana *et al.* very well. Note that both curves evolve smoothly through the transition between crystalline phases. As before, the AIREBO potential is too stiff and dramatically under-predicts the compression. The REAXFF results agree well with experiment at higher pressure, despite the unrealistic expansion along the axial direction with increasing pressure. They gradually rise above experiments as P decreases below 15 GPa.

As a final validation of AIREBO-M, we compute the Hugoniot equation of state for orthorhombic PE.^{20} This state function, given an initial state, encodes the accessible pressure states a material can be driven to by density shock-waves of various speeds. As PE becomes more widely used for protective armors and other impact-oriented material applications, models that accurately describe its shock properties will grow in utility. To generate the T = 300 K Hugoniot curve, we perform a series of NVT and NPT simulations, utilizing the quasistatic methodology of Erpenbeck,^{25} as performed by Chantawansri *et al.*^{26}

For our standard of comparison, we follow other recent studies of shock in PE^{4,26} and use the accepted shock Hugoniot curve for crystalline PE derived by Pastine.^{40} Work by Carter and Marsh has shown Pastine’s crystalline Hugoniot curve to be an accurate ordered limit of the available experimental data on semicrystalline samples.^{41}

Fig. 6 compares our computed Hugoniot curve for AIREBO-M to Pastine’s curve. For comparison, we plot sev eral other Hugoniots computed by Mattsson *et al.*^{4} with various models, including AIREBO, the all-atom optimized potential for liquid simulations (OPLS-AA),^{42} and the exponential-6 model (EXP-6) of Borodin and Smith.^{43} AIREBO-M (red circles) produces a Hugoniot curve for the orthorhombic system that closely follows the accepted crystalline curve, which we interpret to mean our model effectively captures the shock properties of the crystal system. In contrast, AIREBO (blue squares) and OPLS (magenta triangles) greatly over-predict the stiffness of the shock response, while REAXFF (green triangles) gives a Hugoniot that is slightly softer than Pastine’s curve. The EXP-6 (cyan pentagons) shows similar perormance to AIREBO-M. This is unsurprising, since both AIREBO-M and EXP-6 utilize exponential repulsive forces.

## V. SUMMARY AND CONCLUSIONS

We have developed a new potential, AIREBO-M, that describes the high pressure response of hydrocarbons. In the original AIREBO,^{14} intermolecular interactions were modeled with Lennard-Jones potentials that rise too steeply at small separations.^{4} These are replaced by Morse potentials in AIREBO-M. The extra degree of freedom in the Morse potential is used to fit high pressure data while retaining most equilibrium properties of AIREBO. We plan to provide AIREBO-M as part of the AIREBO model driver being developed by S. J. Stuart for the Open Knowledge-base of Interatomic Models (Open KIM).^{44}

The C–C interaction was fit to experimental data for the layer spacing of graphite from 0 to 14 GPa.^{35} To match the experimental results for zero pressure, the equilibrium C–C spacing must be reduced by about 3% from that in AIREBO. This change has little effect on hydrocarbons, where C–C interactions are normally screened by hydrogen. The fit was then tested against DFT-D calculations for graphite bilayers.^{36} Excellent agreement was obtained down to separations corresponding to significantly larger pressures (∼18 GPa) than were achieved in experiments.

The repulsive region of C–H and H–H interactions was fit to dimer interaction energies calculated with Møller-Plesset second-order perturbation theory, while retaining the interaction energy and equilibrium spacing of AIREBO. The resulting potential closely matches the entire attractive range of AIREBO, while softening the hard core repulsion. Thus, the high pressure response is corrected without sacrificing the quality of AIREBO’s fits to ambient pressure data. AIREBO-M reproduces the equilibrium pressure-volume isotherm and Hugoniot equation of state for crystalline polyethylene up to 40 GPa.

The Morse potentials for intermolecular interactions developed here can be paired with different versions of the REBO2 for covalent interactions. The calculations shown above were performed with the REBO2 potential of Ref. 16 as implemented in LAMMPS. This has been shown to overestimate the forces needed to break carbon bonds under tension.^{45} Better agreement can be obtained by screening covalent interactions and extending their range.^{45,46}

There are active efforts to extend REBO2 to other elements, including sulphur, oxygen, and nitrogen. We expect that Morse potentials will also be useful in describing their intermolecular interactions. Fig. 3 suggests that an approximate translation between Lennard-Jones and Morse potentials should use Eq. (3), and then multiply *α* by 0.877. Another useful extension would be to consider configuration dependent intermolecular interactions. For example, carbons with sp2 and sp3 coordination may have different dispersion forces.^{47}

## Acknowledgments

This research was performed within the Center for Materials in Extreme Dynamic Environments (CMEDE) under the Hopkins Extreme Materials Institute at Johns Hopkins University. Financial support was provided by the Army Research Laboratory under the MEDE Collaborative Research Alliance, through Grant No. W911NF-12-2-0022 and by the Modeling Complex Systems IGERT, through Grant No. DGE0801471. We thank Lars Pastewka, Steven J. Stuart, and Bill and Elaine Pappas for useful discussions.