Using Force Matching To Determine Reactive Force Fields for Water

Dec 8, 2016 - We present a method for the creation of classical force fields for water under dissociative thermodynamic conditions by force matching t...
2 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/JCTC

Using Force Matching To Determine Reactive Force Fields for Water under Extreme Thermodynamic Conditions Lucas Koziol,† Laurence E. Fried, and Nir Goldman* Physical and Life Sciences Directorate, Lawrence Livermore National Laboratory, Livermore, California 94550, United States S Supporting Information *

ABSTRACT: We present a method for the creation of classical force fields for water under dissociative thermodynamic conditions by force matching to molecular dynamics trajectories from Kohn−Sham density functional theory (DFT). We apply our method to liquid water under dissociative conditions, where molecular lifetimes are less than 1 ps, and superionic water, where hydrogen ions diffuse at liquid-like rates through an oxygen lattice. We find that, in general, our new models are capable of accurately reproducing the structural and dynamic properties computed from DFT, as well as the molecular concentrations and lifetimes. Overall, our force-matching approach presents a relatively simple way to create classical reactive force fields for a single thermodynamic state point that largely retains the accuracy of DFT while having the potential to access experimental time and length scales.

■

INTRODUCTION Understanding chemistry under high pressures and temperatures in the laboratory can require detailed atomistic and kinetic knowledge about reactants and products under these conditions. However, in many cases insufficient data exist for the equation of state and chemical reactivity of these materials under the extreme conditions attained during experiments.1 For example, studies on carbon-rich materials under pressures up to 30 GPa (1 GPa = 10 kbar) suggest that slow chemical kinetics can predominate beyond the time scales of nanosecond laserdriven compression experiments, even at temperatures of thousands of Kelvin.2 Computer simulations such as molecular dynamics (MD) hold promise as an independent route to determining the equation of state and chemical reactivity under these conditions.3−9 Such studies can provide simple chemical pictures of ionized intermediates and reaction mechanisms and can help identify atomic-scale properties that determine observed macroscopic kinetics. These types of results have proven to make experiments more tractable by aiding in their interpretation and helping to narrow the number of different materials and reactive conditions to investigate.10 Accurate modeling of the breaking and forming of bonds in condensed phases frequently requires use of quantum theories such as Kohn−Sham density functional theory (DFT). DFT © XXXX American Chemical Society

remains one of the most popular and widely used modeling methods in condensed matter physics, computational chemistry, and materials science for prediction of material properties under extreme conditions.11−15 DFT has been shown to accurately reproduce the phase boundaries and thermal decomposition of many materials,16,17 particularly at higher thermodynamic conditions like planetary interiors,12,18−20 where choice of a specific exchange−correlation functional is less important.21 Characterizing the fundamental properties and reactivity of water at extreme temperature and pressure is still a formidable challenge and is important for developing accurate equation of state models and for understanding the reaction pathways of energetic molecules. Recent DFT-MD simulations of an organic explosive found that water exhibited catalytic function and was not a stable byproduct, as previously thought.22 Combined experimental and DFT simulation efforts have proposed that water adopts a superionic phase above ∼30 GPa and 2000 K,23,24 with oxygen atoms frozen in a body-centeredcubic lattice and hydrogen atoms with liquid-like mobility. However, in slightly lower pressure regimes (e.g., 10−30 GPa), Received: July 14, 2016 Published: December 8, 2016 A

DOI: 10.1021/acs.jctc.6b00707 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

conditions33,48 and to improve the accuracy of semiempirical quantum methods.21,49 In this work, we develop new, efficient force-matched force fields for MD simulations of reactive water under extreme thermodynamic conditions, using the simplest functional form viable (two body interactions plus overcoordination). In contrast to previous work, we explore a much simpler force field determined through force matching and ascertain its applicability to a relatively narrow range of conditions where dissociation is a frequent event. Doing so allows us to focus on chemical reactivity specifically, rather than pressure and temperature-induced phase transitions. Due to its rapid chemistry under extreme conditions, water serves as an excellent model system for the development of reactive force fields for systems with several element types, where more complicated force fields might be prohibitively expensive. We find that our force field approach can be rapidly and efficiently optimized to DFT-MD trajectories under reactive conditions. It can thus be used to extend MD simulations to much longer time scales than possible with DFT-MD while potentially retaining much of its accuracy. The accuracy of our forcematched force fields is assessed through comparison to chemical, structural, and dynamical properties from DFT-MD. Additional comparison is made to results from a ReaxFF force field developed for water under extreme conditions, and results from Pinilla et al., where applicable.

water exists in a molecular state, though with the absence of a persistent hydrogen bond network.6,25,26 Furthermore, frequent OH bond breaking leads to extremely short-lived dissociation products H+, H3O+, and OH− as part of rapid intermolecular proton transfer.13,26,27 DFT-MD simulations, though, require immense computational effort per simulation time step, and consequently, are usually limited to picosecond time scales and nanometer system sizes. In contrast, many chemical events occur over nanosecond time scales or longer,28 and experiments can probe micron length scales or beyond.29−31 Difficulty thus generally arises in determining MD models for chemical bonding that are both accurate and computationally efficient. Classical force fields treat interatomic interactions by analytic functional forms and gain several orders of magnitude efficiency over quantummechanical methods, allowing for time and length scales that more closely resemble those of experiments. These methods are widely used in the context of weak (noncovalent) interactions.32−35 Covalent bond breaking and formation, however, is more challenging to describe using simplified potentials,36−39 in part due to the dynamic evolution of the potential energy surface as compounds are destroyed and new ones are formed. Pinilla et al. successfully determined a relatively complex reactive water model based on iterative matching to DFT forces and equation of state data from ambient conditions and discussed the ice VIII−ice X transition and the superionic phase.39 However, their force field was generated exclusively from DFT simulation of water under ambient, nondissociative conditions and underpredicts the transition temperature for solid ice to superionic water. The ReaxFF force field has emerged as a good balance between accuracy and computational cost for reactive materials and conditions38 and has been applied to a diversity of reactive systems, including fullerene combustion,40 detonation of energetic molecules,41,42 and catalysis on metal surfaces.43 Nevertheless, the high complexity of the ReaxFF force field is a disadvantage, in terms of the ability to modify the potential as needed for new systems and conditions. In general, force fields consist of simplified analytical expressions covering all interatomic interactions. These expressions contain a set of parameters, {p}, that must be optimized to reproduce high-quality reference data at similar chemical environments and thermodynamic conditions to which the force field will be used. The accuracy of force fields is known to be limited to conditions close to their reference data, and in general, this lack of transferability is a major limitation compared to quantum-mechanical methods. Historically, reference data have often included experimental or theoretical structural information (lattice constants, bond lengths) or energetic information (binding energies, heats of formation). More recently, the force-matching method has become a common way to optimize force fields. Force matching utilizes atomic forces obtained from a high-quality reference MD simulation (e.g., DFT) to minimize the leastsquares error between it and a given force field model, e.g., χ2 =

∑ (Fi ,DFT − Fi ,{p})2 i

■

METHOD AND COMPUTATIONAL DETAILS Initial DFT-MD calculations were performed on liquid water at two different extreme thermodynamic states: 1.50 g/mL and 2000 K (10.7 GPa), corresponding to moderately dissociative conditions, and 2.25 g/mL and 2000 K (41.3 GPa), corresponding to highly dissociative conditions. These state points are close to the planetary isentropes of Neptune and Uranus, where the high degree of dissociation could contribute to the large magnetic field measurements from these planets.23 An additional force-matched potential was created for superionic water at 2.95 g/mL and 2000 K to match one set of conditions from Pinilla et al., as well. Hereafter, we refer to these conditions by their density, only. All liquid water DFT calculations were performed with the OpenMX code,50 using norm-conserving pseudopotentials51 and the Perdew−Burke− Ernzerhof generalized gradient approximation functional (PBE).52 (NB: Details of our DFT simulations of superionic water are discussed below.) Starting configurations for our MD simulations were taken from previous DFT NVT simulations with PBE that were run for up to 10 ps. DFT has several documented shortcomings for water under unreactive conditions where it remains a molecular liquid.53,54 However, simulations with PBE have been shown to yield accurate physical and chemical results for water under the highly dissociative conditions of our study.55,56 In particular, use of different exchange-correlation functionals yields little to no change in the observed chemistry of materials under extreme conditions.13,14,21 For liquid water, we have employed a double-ζ polarizable (DZP) localized atomic basis set, optimized to be nonzero inside a 5.0 bohr radius.57 Dispersion interactions were included using the DFT-D2 approach.58 All MD simulations discussed in this work were run in the canonical (NVT) ensemble with global a Nosé−Hoover thermostat,59,60 where the system is thermostated to the total system temperature (vs a local or “massive” thermostat, which applies an independent

(1)

The index i runs over atomic coordinates (including x, y, z) in the reference MD trajectory. Force matching was first introduced by Ercolessi and co-workers44 and has been utilized in several applications45−47 including for water under ambient B

DOI: 10.1021/acs.jctc.6b00707 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation thermostat to each particle61). MD simulations using a planewave code62−64 with a cutoff of 1000 eV, the PBE functional, and projector-augmented wave pseudopotentials65,66 yielded molecular concentrations and lifetimes for H2O, H+, OH−, and H3O+ that all agreed with the localized basis set calculations to within 1%. The system at 1.50 g/mL contained 32 water molecules with cell lattice dimensions of 8.609 × 8.609 × 8.609 Å3, whereas the system at 2.25 g/mL contained 54 water molecules with cell dimensions of 8.953 × 8.953 × 8.953 Å3. These dimensions have been shown to be large enough to permit Γ-point sampling only of the Brillouin zone.13 Each DFT calculation was run for up to 10 ps using a time step of 0.2 fs, although force matching to the first 5 ps only was sufficient to yield an optimal potential energy surface. A typical snapshot of each MD simulation is shown in Figure 1. Here, we observe relatively frequent dissociation at 1.50 g/mL, though the simulation mainly consists of H2O molecules. In contrast,

the simulation at 2.25 g/mL yields an even higher degree of dissociation and larger concentrations of H+ and OH−. The force field model introduced here contains three components: short-ranged central forces, electrostatic interactions due to fixed point charges, and a short-ranged manybody overbonding term. The short-ranged central force component was motivated by the central force water model of Stillinger and co-workers.67 In this case, these forces account for both intra- and intermolecular interactions, in contrast to previous models for ambient water, where the intramolecular interactions were treated separately with monomer potentials (e.g., Partridge−Schwenke, or harmonic bonds/angles).48 Separation of intra- and intermolecular potentials requires maintaining a bond topology at every MD step, which precludes straightforward treatment of bond breaking. Our use of a central-force model to include intramolecular interactions is inherently applicable at extreme conditions, where the high internal energy samples vibrationally hot regions of configuration space, largely reducing bond directionality. Our long-range electrostatic force terms contain nonpolarizable atomic charges for oxygen and hydrogen, qO and qH. Unconstrained optimization of these partial charges could yield unphysical results, e.g., qO and qH ≈ 0. Hence, for the work discussed here we have chosen to fix the atomic partial charges to values computed from DFT for water under similar conditions, where qH = +0.365 and qO = −0.710.13,14 We note that computed values of qO and qH vary by only 5−15% over pressures from 8−68 GPa and temperatures from 400−3600 K. Force matching of the central force terms was performed from our DFT computed total atomic forces to linear combinations of Chebyshev polynomials of the first kind, Tn(x): Ecf =

∑ ∑ ∑ Cna a Tnij(x) Fc(rij) i j

i

j>i

n

(2)

Here, Ecf is the central force potential energy neglecting the Coulomb interaction. The indices i and j correspond to each atom and ai and aj to their specific atom types (e.g., oxygen or hydrogen). The Chebyshev polynomials can be taken to arbitrary order n, and the variable x is subject to the constraint −1 ≤ x ≤ +1. Chebyshev polynomials were chosen rather than simple powers of x, because the oscillatory nature of the Chebyshev polynomial leads to better behavior of the polynomial coefficients. Fc(rij) is a smooth cutoff function which is set to zero beyond a set cutoff, e.g., for rij < rmax aiaj : Fc(rij) = (1 − rij/ramax )3 i aj

(3)

The Cn parameters are the linear coefficients to be determined from the force-matching process. Permutational invariance of the potential with respect to atom labeling requires that Cani aj = Canj ai. Force matching then determines a set of n polynomial coefficients for each atomic pair (i.e., O−O, O−H, and H−H). The internuclear distance rij was transformed to the Morse variable, xij = exp(−rij/λaiaj), where λ is the Morse variable range parameter,68,69 defined in this case individually for each type of atomic pair interaction. The Morse variable tends toward zero in the limit of large rij, permitting a physically correct description of the short-range potential in regions of weak or dissociative interactions.70 The Chebyshev polynomial variable x was then set to be within the range [−1, +1] by performing the operation x = (xij − xavg)/xdiff. In this case, xavg and xdiff were

Figure 1. MD snapshots of liquid water under moderate (1.50 g/mL) and highly dissociative (2.25 g/mL) conditions. Here, H2O molecules are color coded blue, OH− are silver, H3O+ are red, and H+ are yellow. C

DOI: 10.1021/acs.jctc.6b00707 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation defined by minimum and maximum allowed rij values according to the equations xavg = 0.5(xmax + xmin)

Analysis of the force autocorrelation function for each DFT simulation revealed that the forces were completely decorrelated by ∼80−100 fs. For the liquid densities, the DFT forces were sampled every ∼160 fs, yielding a total of 31 sets of forces for the 1.50 g/mL system and 28 for the 2.25 g/mL system to use for the force-matching step. Doubling the DFT-force data set by taking configurations every ∼80 fs yielded virtually identical results. In each case, preprocessing of the forces was performed where the Coulombic forces due to the fixed values of qO and qH (computed with Ewald sums) were subtracted from the DFT forces before performing the force matching. A major limitation of the central force approximation is the neglect of many-body effects on covalent bonding. Hybridization leads to the effect that bond-dissociation energies decrease with increasing coordination number of bonded atoms. This limitation was observed in our molecular dynamics simulations using only the central force terms: we observed long-lived H3O+ species in the simulation that were not present in the DFT-MD (discussed below). Hence, to account for this qualitative deficiency, a second set of force-matching calculations were performed for each system with an additional many-body potential energy term to penalize overcoordinated species. This correction is identical in form to the overcoordination term in ReaxFF.38 Using the original ReaxFF notation, this is written as

(4)

xdiff = 0.5(xmax − xmin)

(5)

exp(−rmin aiaj /λaiaj) and xmin min raiaj by the approximate

exp(−rmax aiaj /λaiaj).

where xmax = = We minimum internuclear have defined distance sampled by our DFT-MD simulation for a specific atomic pair. Hence, the Chebyshev variable x was set to a value max of +1 if rij < rmin aiaj and a value of −1 if rij > raiaj . The values used min max for λaiaj, raiaj , and raiaj are shown in Table 1. In this work, the Table 1. Morse Variable Parameters Used for Our ForceMatched Potentials at 1.50 and 2.25 g/mLa density (g/mL)

pair type

λaiaj

rmin aiaj (Å)

rmax aiaj (Å)

1.50

O−O O−H H−H O−O O−H H−H

2.20 2.20 2.20 2.50 2.50 2.50

1.90 0.70 0.90 1.90 0.75 0.85

6.00 6.00 6.00 6.00 6.00 6.00

2.25

a

The Chebyshev polynomial order was set to 12 for both liquid densities studied here.

⎛ ⎞ 1 ⎟⎟ Ej ,over = pover *Δj *⎜⎜ ⎝ 1 + exp(λ6*Δj ) ⎠

Chebyshev order was set to 12 for all O−O, O−H, and H−H interactions for the liquid densities, and to eight for superionic water. Increasing the Chebyshev order to significantly higher values (e.g., 48) yielded only minimal improvement in accuracy. The linear coefficients Caniaj were readily obtained by linear least-squares fitting using the singular value decomposition (SVD) method. SVD allows for the optimal coefficient values to be found while limiting overparametrization through the use of a lower cutoff bound on the magnitude of the singular values. Because the potential energy is a linear function of all fitting parameters, we have removed the need for nonlinear optimization techniques, such as iterative direct gradient minimization approaches that can become trapped in a local minimum (e.g., Levenberg−Marquardt), or global minimum searches (e.g., simulated annealing), which are computationally expensive. An additional close-ranged repulsive energy function was added to the central-force interactions to prevent sampling of interatomic distances smaller than rmin aiaj during MD simulation, where the force-matched potential could possibly yield unphysical results. Here, we used a penalty function Vp in terms of a penalty distance rp, defined as follows: Vp = ps rp3

(6)

min ⎧ rij − pd < ramin ⎪ ra a + p − rij d i j i aj rp = ⎨ ⎪ ⎩0 otherwise

(7)

(8)

The parameters pover and λ6 are to be fitted, and Δj is the difference between an oxygen atom’s bond order and its valency,

Δj ,O =

−2 ∑ BOOH ij (9)

i ,H

BOOH ij

In ReaxFF, the bond order is a function of interatomic distance only. We limit our force field to simulations with maximum OH bond order of 1 (e.g., a single bond). The ReaxFF bond order is then given as ⎛ ⎛ rij ⎞ pbo,2 ⎞ ⎜ BOOH exp = p · ij ⎜ bo,1 ⎜ r ⎟ ⎟⎟ ⎝ 0⎠ ⎠ ⎝

(10)

Equation 10 describes an S-shaped curve (pbo,1 is negative) decreasing from 1 to 0 around r0, which should be close to the equilibrium bond length for an OH bond. The form of eqs 8−10 is such that Ej,over is positive when the total molecular bond order is greater than oxygen’s valency of 2 and rapidly decreases to 0 for a total bond order less than 2. A detailed description of the ReaxFF functional forms can be found in ref 38. Inclusion of eq 10 in our force-matching calculations required the use of a nonlinear least-squares optimization approach. This was done in a two-step, self-consistent fashion where first, a single iteration of the Powell method was used to choose new values for the nonlinear parameters in the overbonding potential (λ6, pbo,1, pbo,2, and r0). Immediately following the Powell iteration, the remaining linear coefficients (pover from eq 8 and the Cn coefficients for the Chebyshev polynomials) were determined though an SVD calculation before proceeding with the next Powell optimization step. Additional manual tuning of the pover parameter was performed at the conclusion of the Powell optimization to yield better

The parameter ps is a penalty function scaling factor, set to a value of 106, and pd is the penalty distance, set to a value of 10−2 for this work. This prevents simulations from occasionally achieving intermolecular separations less than those sampled in the DFT MD simulations. Testing showed that MD simulation results were insensitive to specific parameter values. D

DOI: 10.1021/acs.jctc.6b00707 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation agreement with H2O molecular lifetimes and oxygen and hydrogen diffusion constants from DFT.

■

RESULTS AND DISCUSSION We present results from potentials created by force matching with Chebyshev polynomials (labeled FM), as well as results from potentials that also include the many-body overbonding function (labeled FM-OB). Tables of parameter values for each central force and overbonding potential are included in Supporting Information. We first discuss results from 1.5 and 2.25 g/mL, only. The root-mean-square (RMS) error in the forces these four potentials (FM and FM-OB at each liquid density) and system pressures from NVT simulation are shown in Table 2. Our expectation is that each force-matched potential Table 2. Root Mean Square Error and System Pressures from the FM and FM-OB Force Fields for Each Densitya density (g/mL)

force field

RMS error (kcal/mol Å)

pressure (GPa)

DFT pressure (GPa)

1.50

FM FM-OB FM FM-OB

14.19 13.86 17.33 17.36

−1.3 5.9 7.6 9.5

10.7

2.25 a

41.3

Pressures from DFT-MD simulation are shown in the final column.

applies to its fitted thermodynamic state point, only. We note that force matching can result in incorrect pressures for a periodic system due to the volume dependence of the DFT energy that is generally neglected with this method.21 This is particularly true for water under extreme, dissociative conditions, where it can exhibit metallic properties.55 Improved accuracy for the pressure can be achieved by matching to components of the stress tensor during the fitting process.21 However, in this work we choose to focus on the chemical and dynamic properties from the FM and FM-OB potentials, exclusively. Inclusion of the overbonding potential with manual tuning of pover had a minimal effect on reducing the overall RMS error, although the effect on the resulting chemical reactivity was more pronounced. The distribution of forces on oxygen and hydrogen atoms for both FM and FM-OB potentials are shown in Figures 2 and 3. At 1.50 g/mL, both FM and FM-OB potentials agree well with results from DFT-MD, although each yielded slightly less accurate results for low-magnitude forces acting on the hydrogens. At 2.25 g/mL, we observe similar agreement, with both potentials again performing slightly more poorly for the low-magnitude hydrogen forces. Regardless, the distribution of forces appears to be well represented by both FM and FM-OB approaches at both sets of conditions. Analysis of the shortranged central-force part of our models for both densities (Figure 4) indicates that for the FM-OB force field, inclusion of the overbonding potential resulted in a deeper potential well depth for O−H bond dissociation. In addition, the repulsive wall for the HH interaction at 2.25 g/mL for the FM-OB model is significantly stiffer than that from the FM model. MD simulations were performed at each density using both FM and FM-OB potentials. Simulations were performed using the same system sizes as our DFT-MD calculations and an identical time step, for up to 100 ps. Doubling of the system size in each Cartesian direction (yielding 256 H2O molecules at 1.50 g/mL and 432 H2O molecules at 2.25 g/mL, respectively) produced nearly indistinguishable results for all properties

Figure 2. Distribution of total forces acting hydrogen (a) and oxygen (b) at 1.50 g/mL. The solid black lines correspond to results from DFT, the dotted red line to force-matching pairwise interactions only (FM), and the solid red line to force matching with additional manual tuning of overbonded terms (FM-OB).

discussed here. In each case, comparison is made to a ReaxFF force field developed for pressures and temperatures similar to those studied from our study,71−73 simulated with the same system size for 100 ps. Radial distribution functions (RDF or g(R)) at 1.50 g/mL (Figure 5) indicate that overall the comparison between all four calculations (DFT-MD, FM, FM-OB, and ReaxFF) for the OO distribution is favorable. All four simulation types yield very similar positions of the first, second, and third peaks, with only a slight difference between peak heights and values of the minima between peaks. In contrast, for the OH RDF, the FM and FM-OB calculations yield much closer agreement with the height of the first peak from DFT-MD as well as the nonzero height of the first minimum, which is indicative of molecular dissociation. The larger first peak from ReaxFF including the small value for the height of the first minimum is indicative of a smaller degree of O−H bond dissociation. The HH RDF from the FM and FM-OB simulations show some discrepancies with results from DFT-MD, with first peak positions at ∼1.9 Å (compared to ∼1.6 Å from DFT-MD) and peaks that are somewhat more broad. In addition, the FM and FM-OB simulations do not yield the slight second and third peaks at ∼2.5 and ∼3.4 that are present in the DFT-MD simulation. It is possible that these discrepancies are due to the errors in low E

DOI: 10.1021/acs.jctc.6b00707 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

lying forces that are present in our force fields and that a higher fidelity force-matching result such as through explicit calculation of three-body forces would rectify these errors. The ReaxFF simulation is able to get the correct position of the first peak in the g(RHH), though the peak is too high, and the remaining part of the RDF also does not include the slight second and third peaks seen in the DFT-MD result. In general, we observe strong agreement between the RDF’s for 2.25 g/mL from the FM, FM-OB, and DFT-MD simulations, and somewhat worse agreement from ReaxFF (Figure 6). For the OO RDF, all four simulation methods yield similar peak positions, with the FM-OB result showing the closest agreement with DFT-MD. The results from FM indicate slightly lower peak heights, though the result from ReaxFF shows somewhat higher peak heights for the first solvation shell in particular. The FM and FM-OB simulations exhibit close agreement with DFT-MD for the OH RDF, though the FMOB result indicates a small peak at distances slightly less than 2.0 Å that is absent from both FM and DFT-MD. In contrast, the results from ReaxFF show a higher first peak intensity and a significantly lower first minimum, indicative again of a smaller amount of molecular dissociation relative to DFT-MD. Finally, the FM and FM-OB simulations show excellent agreement with DFT-MD for the HH RDF at this density, whereas, ReaxFF yields a peak height for the first solvation shell that is too high and is at a slightly lower distance relative to DFT-MD. Similarly to earlier work,20,74,75 we have determined the specific chemical products formed in our simulations using a pre-established methodology of optimal bond cutoff distances and lifetimes.16 In this work, molecules were identified by having all bonds existing for a minimum lifetime of 20 fs. This time period corresponds to about two periods of a fast vibration (OH stretch). The bond distance criteria were based on the first minimum of the corresponding RDF. With these bonding criteria, specific molecular species were then defined by recursively creating a data tree of all atomic bonds branching from the original bonded pair. At 1.50 g/mL, the mole fractions

Figure 3. Distribution of total forces acting hydrogen (a) and oxygen (b) at 2.25 g/mL. The solid black lines correspond to results from DFT, the dotted red line to FM, and the solid red line to FM-OB.

Figure 4. Short-ranged central force contribution to the FM (dotted red line) and FM-OB (solid red line) force fields at 1.5 and 2.25 g/mL. F

DOI: 10.1021/acs.jctc.6b00707 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 5. Radial distribution functions for water at 1.50 g/mL and 2000 K. The same legend applies as from Figure 3, with the gray line corresponding to ReaxFF.

conditions, each simulation yielded a highly reactive mixture, with molecular lifetimes for all chemical species well below 0.100 ps. Consequently, the overall chemistry is extremely dissociative and exact comparisons between predicted concentrations and lifetimes are arguably less significant. DFT-MD calculations yielded a mole fraction and lifetime for H2O of 0.710 and 0.056 ps, respectively, whereas the FM calculations yielded values of 0.399 and 0.026 ps, and the FM-OB calculations yielded values of 0.555 and 0.043 ps. The predicted concentrations of H+ and OH− from the FM and FM-OB calculations are somewhat higher than those from DFT-MD (e.g., 0.237 and 0.191 vs 0.139 for the OH− mole fraction for FM and FM-OB vs DFT), though all three sets of simulations yield lifetimes for H+ and OH− of ∼0.020 ps. However, the FM simulation yielded a mole fraction of 0.109 and lifetime of 0.021 ps for H3O+ whereas FM-OB and DFT-MD yielded similar mole fractions of 0.019 and 0.034, and lifetimes of 0.012 and 0.013 ps, respectively, for this species. In contrast, ReaxFF did yield a mole fraction and lifetime for H2O that compared more closely to results from DFT-MD. However, ReaxFF also yielded exceedingly small amounts of OH− (mole fraction of 0.006) as well as a number of chemical species that were not present in any of the other simulations, such as an abundance of H2O2 (mole fraction of 0.242). In addition, the ReaxFF simulations did not yield any H+ or H3O+, as was found in our other simulations.

of all chemical species were renormalized to exclude extremely short-lived water dimers and higher order clusters that generally occurred at very low concentration (