Benchmark of ReaxFF force ﬁeld for subcritical and supercritical water

Water in the subcritical and supercritical states has remarkable properties that make it an excellent solvent for oxidation of hazardous chemicals, waste separation, and green synthesis. Molecular simulations are a valuable complement to experiments in order to understand and improve the relevant sub- and super-critical reaction mechanisms. Since water molecules under these conditions can act not only as a solvent but also as a reactant, dissociative force ﬁelds are especially interesting to investigate these processes. In this work, we evaluate the capacity of the ReaxFF force ﬁeld to reproduce the microstructure, hydrogen bonding, dielectric constant, diffusion, and proton transfer of sub- and super-critical water. Our results indicate that ReaxFF is able to simulate water properties in these states in very good quantitative agreement with the existing experimental data, with the exception of the static dielectric constant that is reproduced only qualitatively. Published by AIP Publishing. https://doi.org/10.1063/1.5031489

Water in the subcritical and supercritical states has remarkable properties that make it an excellent solvent for oxidation of hazardous chemicals, waste separation, and green synthesis. Molecular simulations are a valuable complement to experiments in order to understand and improve the relevant sub-and super-critical reaction mechanisms. Since water molecules under these conditions can act not only as a solvent but also as a reactant, dissociative force fields are especially interesting to investigate these processes. In this work, we evaluate the capacity of the ReaxFF force field to reproduce the microstructure, hydrogen bonding, dielectric constant, diffusion, and proton transfer of sub-and super-critical water. Our results indicate that ReaxFF is able to simulate water properties in these states in very good quantitative agreement with the existing experimental data, with the exception of the static dielectric constant that is reproduced only qualitatively. Published by AIP Publishing. https://doi.org/10.1063/1.5031489

I . INTRODUCTION
Pressure stabilises liquid water above the usual room temperature boiling point of 373.2 K, in a state called superheated or subcritical. At even higher temperatures, above 647.1 K and pressure of 22.06 MPa, water enters in the supercritical regime, a fluid state where the liquid and gas phases cannot be distinguished, and the density evolves smoothly as the external conditions move along an isotherm or an isobar. In nature, sub-and super-critical water (both called SC herein) appears as a free fluid or embedded into rocks and minerals in earth's mantle, playing an important role in tectonics as a lubricant and in rock melting. [1][2][3] From a technological point of view, SC water is used as a solvent for the synthesis of nanoparticles hazardous waste oxidation and chemical separation. [4][5][6][7][8][9] All these applications are possible due to the unique properties of water in the SC state. The high temperature and intermediate density between the room temperature liquid and the steam modify the microstructure of water. The mediumrange order found at room conditions is lost, 10 the coordination of water molecules changes from tetrahedral to trigonal, 11 and the number of hydrogen bonds (HBs) between water molecules is considerably reduced. 12 The structural changes are accompanied by a decrease in the dielectric constant, 13 an a) Electronic mail: hegoi.manzano@ehu.eus increased diffusivity, and lower viscosity. 14 As a consequence, SC water has a lower solvation capacity of hydrophilic species like dissolved ions, which tend to form ionic pairs, and the reactivity increases, reducing synthesis times from hours to seconds. 4 Also, SC water becomes a good solvent for organic solutes. 9,15 Even more important, all these properties can be easily controlled tuning the temperature and pressure, achieving the desired conditions for a given synthesis, oxidation, or extraction.
As in many other fields, atomistic simulations could be a valuable complement to experiments to fully understand SC processes. Ab initio methods are theoretically the most appropriate methodology, yet usually of limited applicability due to the size and time scales needed to study the reactivity of species in solution. Reactive empirical potentials [16][17][18] are capable of simulating reactions in systems with thousands of atoms during tens of nanoseconds, making them an excellent alternative. In addition, the fast reaction kinetics under SC conditions make the gap between experiments and simulations smaller.
In this work, we focus on the ReaxFF force field. 18 ReaxFF has been used to simulate supercritical state processes such as wastewater purification and hydrogen production. [19][20][21][22] However, the validity of the force field to reproduce the basic properties of SC water has never been systematically studied. To cover this blank space, we performed molecular dynamics (MD) simulations with ReaxFF of several supercritical and subsupercritical states. For each thermodynamic state, the most relevant characteristics of the supercritical state were tested: density, dipole moment and static dielectric constant, microstructure, hydrogen bonding, diffusivity, and proton transfer. All the results are compared with the existing experimental data or with results from other water models.

II . COMPUTATIONAL METHODS
All the simulations were done using the ReaxFF reax/c implementation 23 on the LAMMPS simulation code. 24 Initially, 309 water molecules were randomly packed in a 2.1 nm side cubic simulation box, 25 and the atomic positions were relaxed by conjugated gradient minimisation to a local energy minimum. Then, random velocities were assigned to the atoms from a gaussian velocity distribution at the desired temperature. The molecular dynamics simulations were done in the isobaric-isothermal ensemble, with the Nose-Hoover style thermostat and barostat and relaxation times of 25 and 100 steps, respectively. 26 The simulation boxes were kept orthogonal and the deformation was coupled in the 3 dimensions. A Verlet integration scheme 27 was used, with a time step of 0.2 fs. The simulations were run for 2.5·10 6 steps, reaching a total time of 500 ps. The equilibrium was typically achieved in less than 10 ps in all the cases, and half of the simulation time, 200 ps, was used for production. The trajectory steps for property analyses were collected every 0.2 ps. For the proton hopping analysis, extra 50 ps runs were done after adding excess hydrogen atoms, and finer trajectories were collected to capture the fast hops, every 0.05 ps for a total time of 20 ps.
The calculations were performed at 6 temperatures (400 K, 500 K, 600 K, 700 K, 800 K, and 900 K) along 4 isobars (25 MPa, 50 MPa, 75 MPa, and 100 MPa). These thermodynamic states include supercritical and subcritical or superheated water. In addition, bulk water in ambient conditions (293 K and 0.1 MPa) was computed for comparison of structure and dynamics.

A . ReaxFF water model
There are uncountable number of specific purpose water models, 28,29 each with its strengths and weaknesses. First, it must be mentioned that ReaxFF O/H parameters are not specific for water. The force field formalism 30 requires that each chemical element is described by a single species and the same parameters describing H 2 O must be able to describe, for instance, hydroxyl bonds in surfaces 23,[31][32][33][34] or O-H interactions in bulk materials. [35][36][37][38] Nevertheless, within the typical classification of water models, 28 ReaxFF could be labeled as a three-center, flexible, polarisable, and dissociative model.
Most water models describe the attraction and repulsion forces between water molecules by a classic Coulombic function and the short-range repulsion via a Lennard-Jones potential centred at the oxygen. In addition to the intermolecular interactions, flexible models include functions to account for the vibrational internal degrees of freedom, usually harmonic potentials for both the bond stretching and angle bending. Typically, in ReaxFF, the description of inter-and intramolecular potential energy includes the following terms: (1) The first four terms (U bond , U angle , U torsion , and U over ) describe the short-range interactions between bonded atoms. They correspond to the interatomic two-body bonding energy, the three-body angular energy, the four-body torsion energy, and a penalty energy to avoid over coordination. Obviously, the torsional term is not needed in the case of water. These short-range terms are bond order-dependent. The bond order is computed at each integration step using an empirical formula. 18,30 If the bond order between a given pair of atoms becomes zero, they do these four interaction energies and the chemical bond is effectively broken.
Next, U vdW and U Coulomb describe long-range interactions. They are bond order-independent, and unlike common models where these contributions are excluded from bonded atoms, they are computed for every pair of atoms in the system. Polarisability has been suggested important to describe water properties and especially for transferability to investigate interfaces or solutions. 28 Flexible models have implicit polarisability due to changes in the molecular geometry. In addition, there are various methods to include polarisability explicitly, 28,29 such as adding extra sites carrying a charge, linked to water oxygen by a Drude oscillator, 39 or by the interaction between atomic multipoles. 40 In ReaxFF, polarisability is taken into account by the atomic charge fluctuations. At each integration step, a charge equilibration method is used to compute on the fly the atomic charges. Thus, variations in the molecular geometry and interaction with environment (surfaces, dissolved ions, other solvent molecules) imply variations in the molecule dipole moment. Typically, ReaxFF uses the electronegativity equalisation method (EEM) 41 to compute the atomic charges, yet the more sophisticated ACKS2 method 42 has been recently implemented.
Finally, the U specific takes into account system-specific contributions to the energy and are included only if needed. For water, specific terms include hydrogen bonding and lone pairs, both being bond order-dependent. For the detailed expressions and functional forms, see Ref. 43.
The first ReaxFF water parameters were developed in the context of hydrocarbon oxidation. 43 These parameters described gas phase water properties but failed to reproduce liquid conditions. To avoid reparametrisation of the whole hydrocarbon set, a new set of water parameters were then developed for aqueous systems, 34 giving rise to the division of ReaxFF in "combustion" and "aqueous" branches. 30 Recently, the aqueous water parameters were modified for a better description of bulk water density and proton transfer reactions. 44 The new set of parameters, named ReaxFF2017, have been further modified to improve the weak interactions between hydrocarbons and water, refitting the parameters related to the vdW force, hydrogen bonding, and Coulomb terms. The force field set is described in detail in Ref. 45 and will be herein named "ReaxFF2017-weak."

A . Density
The ReaxFF2010 water parametrization 34 underestimated bulk water density by almost 10%. 44 This can be considered as a poor performance since most empirical potentials are capable of good agreement in water density for a wide range of temperatures and pressures. 28,29,46,47 The ReaxFF2017-weak version corrects such disagreement for bulk water under ambient conditions and according to our simulations, 45 so it does for the range of SC conditions simulated in this study. In Fig. 1, the density at different temperatures along four isobars is presented. The lines correspond to the experimental data 48 and the points correspond to the simulation results. As the temperature increases along an isobar, the density decreases following a sigmoidal trend. Below the supercritical temperatures (400 K, 500 K, and 600 K), subcritical water can be considered as a superheated liquid, as reflected by their high densities ρ ≥ 0.7 g cm −1 . For the lower pressure, 25 MPa, the system is close to the critical point and there is an important change in density close to the critical T c = 647 K. It is possible to distinguish the liquid-like and gas-like zones. As the pressure increases, we move into the supercritical region and the density change along an isobar is smoother. All these features are captured by the ReaxFF2017weak, and since most of supercritical water properties depend to a great extent on its density, the capacity to reproduce it accurately is the first good symptom of the force field quality.

B . Internal structure
SC conditions have little impact on the internal structure of water molecules. Figure 2 shows a characteristic bond length-bond angle distribution. We did not find any trend in the bond distance variation with temperature or density. The most probable O-H distance averaged from all the pressure and temperatures is r OH = 0.981 ± 0.012 Å, in agreement with the accepted experimental value of r OH = 0.97 Å. 49 It must be noted that the standard deviation is smaller than that from the fluctuations in a single simulation (±0.07 Å) and therefore without statistical significance. The φ angle is 104.93 • for the most dense conditions, again in good agreement with experiments, yet closer to the gas phase (104.52 • ) than to the liquid phase (106.06 • ). 49 It decreases linearly with density by just 1 • . Despite the trend, the reduction is again smaller than the standard deviation for a single thermodynamics state (±5 • ) and almost negligible, so we can consider that overall SC conditions have no impact on the internal structure of water molecules.

C . Atomic charges, dipole moment, and dielectric constant
While the internal structure remains in practice constant in the studied range, the atomic charges follow a clear trend with density. As the density decreases, so do the atomic charges. Figure 3 shows the normalised charge distributions for the highest (400 K, 100 MPa) and lowest (1000 K, 25 MPa) density cases. For the lowest density, the most probable charge on oxygen atoms is q O = −0.614 and on hydrogen atoms is q H = 0.307, the values are in good agreement with charges derived from MP2 electrostatic potential for a gas phase (isolated) water molecule. 50 These values increase up to q O = −0.719 and q H = 0.360 in high-density states due to stronger intermolecular interactions. The increase in the charge separation is consequently translated into the dipole moment distribution. At low density, the dipole moment µ = 1.77 D is slightly smaller than the MP2 value 50 and the experimental value of gas phase µ = 1.85 D. 51 As the density increases, so does the charge on atoms, and evenly the dipole moment reaches a maximum value of µ = 2.03 D. This maximum value should compare with the dipole moment of bulk water under room conditions µ = 2.8 D, but it is about 25% lower. Nevertheless, the increase in dipole moment with temperature and pressure is captured, and so it is the small but non-negligible distribution wideness that is predicted by Density Functional Theory (DFT) simulations. 52 Microscale changes in the molecular dipolar moment are translated in macroscopic changes in the dielectric constant. In general, the dielectric constant of SC water is reduced with respect to the bulk value, a relevant effect that favours kinetics during SC synthesis and solvation of organic compounds in supercritical oxidation. 4 where 0 is the vacuum permittivity, V is the average system volume, k b is the Boltzmann constant, and T is the average temperature. Since the water dipole moment computed with ReaxFF is lower than the experimental value, we could expect an underestimation of the dielectric constant. However, for bulk water and room conditions, the dielectric constant computed from ReaxFF dipole fluctuations and using Eq. (2) is 84.8 (not shown), ≈8% larger than the experimental value (78.5). This can be considered good agreement since the relative error of obtained from different force fields 28,29 is commonly higher than 50% and even different studies with the same force field 47 report values with more than 30% difference. Such a result is likely due to a high correlation of the dipole moments that overestimates the value of the dipole fluctuations. Or in other words, the Kirkwood g-factor is probably overestimated by the force field. The static dielectric constant obtained from ReaxFF as a function of the density for different temperatures is showed in Fig. 3. At the lowest temperature (400 K), and thus highest densities, the static dielectric constant fluctuates between 50 and 60 depending on the pressure, a value in good agreement with the experimental data. 13 For higher temperatures, ReaxFF overestimates the static dielectric constant, yet overall the experimental trend is captured well, with an asymptotic decrease in as the density goes to zero. Again, this might be due to the mentioned overestimation of the Kirkwood g-factor that persists over temperature.

D . Water microstructure
The microstructure of water is usually analyzed by looking at the site-site pair distribution functions [radial distribution function (RDF)], g(r). At room temperature, the ReaxFF2017-weak gives an accurate description of the sitesite RDFs, 45 comparable with other common models. Only a small underestimation of the g O-O (r) fist peak intensity can be noted. The extension of the good agreement at room conditions to SC conditions must be evaluated. In Fig. 4, we show the intermolecular RDFs for various site combinations at different temperatures, together with those for the ambient conditions reference state. As in Ref. 54, we also found that pressure has little impact on the RDFs since plots at the same temperature but different pressures do not show important differences. In other words, the temperature is the variable influencing most the microstructure of water. In general, the RDFs show the desired trends: there is an increasing broadening of the short-range bands as the temperature increases, and the long-range structure is lost. These changes are clear in the g O-O (r), where the long-range structure hallmark of water tetrahedral arrangement 55 is lost. However, the band at lower distances in the g O-H (r) indicates that hydrogen bonds are still present, yet with a less defined structure as temperature increases. 56 Finally, the two well-defined peaks of the g H-H (r) progressively merge as the temperature increases. 55 For further analysis, Fig. 4 also shows a comparison between simulated RDFs and those obtained from Neutron Diffraction with Isotopic Substitution (NDIS) experiments. 10 The NDIS-97 dataset 57 for ρ = 0.72 g cm −1 and T = 573 K was compared with our simulations at ρ = 0.70 g cm −1 and T = 600. The pressure in the simulated states is 25 MPa, while in the experiment it should be ≈20 MPa according to the water phase diagram. 48 Although the thermodynamic states are not exactly the same, they should be similar enough for a qualitative comparison. In general, the ReaxFF2017-weak RDFs match very well with those derived from NDIS experiments. The good match on the onset of the RDFs indicates that the exclusion volume is perfectly captured and the intensities and coordination numbers in simulation are in good agreement with experiments. The position of the first peak on the g O-O (r) and its decay is well captured, except for the small bump at 3.75 Å, which is not present in simulations, suggesting that the exclusion volume from the first coordination shell is slightly less structured in ReaxFF. The same conclusion can be extracted from the decay of the second peak of the g O-H (r) RDF. Nevertheless, the short-range peaks originated from the hydrogen bonds are perfectly reproduced. Finally, despite the little correlation in the g H-H (r), the simulation matches well with the experiment, with the exception of a sharper first peak that has the maximum at ≈0.4 Å shorter distance. Overall, ReaxFF reproduces the RDFs and hence SC water microstructure, with good accuracy.

E . Hydrogen bonding
In the SC regime, water's hydrogen bond (HB) network is substantially altered. It is generally accepted that hydrogen bonds persist to temperatures higher than 800 K, 12 although their number decreases substantially as temperature increases. 12,62-65 We used a simple geometric criterion to define hydrogen bonds based on the bond length and angles. 67,68 It is considered that a hydrogen bond is formed if the distance between the oxygen from the HB-donor molecule and the oxygen HB-acceptor molecule is smaller than 3.5 Å and the angle by the vectors formed by the mentioned oxygen atoms and the O-H bond of the HB-donor water molecule is smaller than 30 • ; see Fig. 5. It has been demonstrated that the geometric criterion fails at high temperatures and pressures. Due to the high density induced by pressure, in a single trajectory snapshot, two water molecules could be close enough to fulfil the geometric criterion, yet be energetically repulsive. 60 In such a scenario, a combined geometric plus energetic criterion has proven to have more physical insight. 69 In any case, both the geometric and the combined criteria predict a very close number of hydrogen bonds up to a limiting value of n HB = 2, which in terms of pressure and temperature corresponds to 300 MPa and 773 K, respectively. Furthermore, the performance of different criteria might depend also on the water model choice since other studies found that energy-based HB definition underestimated the n HB = 2 and the geometrybased definition was more appropriate. 59 For simplicity, we chose a combined length and angle probability analysis. In Fig. 5, we show the results for 3 supercritical states and the room temperature and pressure reference state for comparison. In the latter, the hydrogen bonds are well defined. The relevant probability area is essentially confined to 3 Å and 25 • , well beyond the cutoffs. For a system at 400 K and 25/100 MPa, with a density comparable with the room conditions state (see Fig. 1), the HB region is also well defined, with slightly large distances and angles. The low probability areas beyond cutoffs start to gain importance and the most probable areas are more diffused. The system at 800 K and 100 MPa has the highest possible density at the highest studied temperature, and therefore the most susceptible case to errors. The most probable combination of bond lengths and bond angles is still within the defined range, but the structure of the hydrogen bonds is less defined, and the probability tends to a constant value for all the configurational space without prohibited areas. This graph suggests that the geometry-based criterion might start to overestimate the number of hydrogen bonds at the highest studied temperatures.
In Fig. 5(b), we plot the average number of hydrogen bonds per water molecule n HB as a function of density for the studied pressures. The n HB decreases nearly linearly with the density from the room temperature value, 3.73 with ReaxFF. It is interesting to note that all the isobars collapse into the same line. The results match very well with previous simulation studies using different water models and methods: molecular dynamics with the SPCE model, 58 the SPC model, 59 and the RPOL model 61 as well as Monte Carlo with the TIP4P model. 60 The largest deviation corresponds to a SPCE study that 58 overestimates n HB for the whole density range, probably because a geometric criterion based only on distances was used. These results suggest that ReaxFF describes hydrogen bonding with the same accuracy as other models and serves as a posteriori validation of the used geometric criterion. J. Chem. Phys. 148, 234503 (2018)  58 squares-SPC, 59 rhombus-TIP4P, 60 and circles-RPOL. 61 (c) Average number of hydrogen bonds per molecule along the studied isobars as a function of temperature. The open symbols correspond to experimental data: rhombus, 62 circles, 12 squares, 63 triangles, 64 and inverted triangles. 65 Dashed lines in (b) and (c) indicate the percolation HB network limit. 66 Figure 5(c) also shows a comparison of the n HB with experimental results as a function of the temperature for the studied pressures. We chose this plot instead of rewriting the experimental data in terms of pressure for a good comparison with previous well-known plots from the literature. 60,70 As the temperature increases, both simulations and experiments show a progressive decrease in the n HB , going to zero for very high temperatures. There is good agreement between experiments and simulations in the supercritical regime, above 600 K. The apparent disagreement below that point is due to the different thermodynamics states since experiments correspond to room pressure conditions and their density is smaller. In this case, the n HB from different isobars do not collapse, which indicates again that the number of hydrogen bonds depends on the density.
It is generally accepted that supercritical water forms a heterogenous fluid with dense patches of molecules connected through HB and low-density zones with isolated molecules or small clusters. From a geometrical point of view, it is described as a transient gel, until the n HB = 1.58 when percolation is lost. 66 That threshold is showed as a dashed line in Fig. 5. We performed a cluster analysis of the simulated trajectories using OVITO code 71 to check if the percolation threshold is captured by ReaxFF. In Fig. 6, snapshots of the analysis are presented for the 75 MPa isobar, using the HB criterion to define if two molecules are linked. At a low temperature, 500 K, the largest cluster, showed in black, spans through the whole unit cell. As described before, there are few small clusters of 2-4 or isolated molecules in low-density regions. At 600 K, the number of small clusters increases, yet the system is still percolated. In agreement with the n HB threshold, the HB network percolation is lost between 600 K and 700 K. The largest cluster is composed by a considerable number of molecules, about 20% of the total, but percolation is clearly lost. At 800 K, even the largest clusters are formed by few molecules, less than 10, and the small clusters prevail.

F . Diffusion
High temperature and low density imply a larger diffusion of water in the SC regime. We computed the selfdiffusion coefficients (D) from the mean square displacement of the molecular dynamics trajectories by the following relationship: The obtained values of D are plotted in Fig. 7 as a function of the density, together with experimental results from nuclear magnetic resonance (NMR) 72,74 and quasi-elastic neutron scattering. 14,73 The self-diffusion coefficient increases exponentially as the density decreases. At densities close to 1 g cm −1 , the value is in the range of 8 cm 2 s −1 ·10 −5 , slightly larger than the bulk self-diffusion coefficient. 75 As for the average number of hydrogen bonds, the diffusivity values of the four isobars collapse in a single curve, indicating that density is the main factor governing diffusion. The agreement between ReaxFF and experimental values is good in the whole range of densities. Most force fields are also capable of reproducing the self-diffusion coefficient at medium-high densities 59,61,70,76,77 within the typical 10% error expected both from the simulations and the experimental data. 72 At very low densities, ρ ≤ 0.2 g cm −1 , the relative errors can be higher, yet the trend is generally captured.

G . Proton transfer
Finally, we study water reactivity in the SC conditions. At saturated water pressure, water self-dissociation increases with temperature up to a maximum at 523 K and then decreases again. Pressure increases self-dissociation and shifts the maximum temperature to higher values. Within the range explored in this work, the maximum expected concentration of H 3 O/OH ionic pairs is about 4.5·10 −11 molecules nm −1 at 600 K FIG. 7. Self-diffusion coefficient of supercritical water obtained from ReaxFF molecular dynamics. Open symbols represent experimental data: triangles, 72 rhombus, 14 inverted triangles, 73 and squares. 74 and 100 MPa, 78 and consequently, we did not observe any self-dissociation in our simulations.
Hence, we turn to evaluate the proton transfer in systems with excess protons. For each thermodynamic state, an extra hydrogen atom was added to the final configuration of each trajectory in a random position, and after equilibration, the number of proton hops taking place during 10ps were counted. In Fig. 8(a), we show a characteristic hopping scheme plotting the id of the hydronium ion oxygen atom (H 3 O+) evolution as a function of time for 2 different temperatures, 600 K and 900 K, at 50 MPa. The hopping profiles show two types of events. On the one hand, there are resonant periods in which the H atoms hop very fast between two water molecules. Those proton transfer events do not translate into effective proton diffusion since the excess H is actually shared between the two water molecules forming what is known as "special-bond." 79 Such a picture has been attributed to the delocalisation of the proton in transient H 5 O 2 + or H 9 O 4 + complexes. 79 On the other hand, there are effective diffusive proton transfer events when an excess proton moves from the special bond regime toward a neighbouring molecule. 80,81 Both temperature profiles are in general similar, but at high temperature there are less hops of both types. Also, resonant periods are shorter and there are intervals in which the excess hydrogen is localised in a single water molecule for ps forming the hydronium ion (H 3 O+). The localisation is due to the lack of neighbouring water molecules to form complexes at very low densities. Overall, the picture is consistent with previous studies. 81 To evaluate the effect of SC conditions on proton transfer, we plotted in Fig. 8(b) the average hopping rate (R hop ) obtained from five simulations as a function of temperature. The computed rates include both types of proton transfer events, the resonant between a pair of molecules and the diffusive ones. Counting this way, the obtained hopping rates are in the range between 5 and 20 ps −1 , in good agreement with Car-Parrinello MD simulations. 79,83 The hopping rate increases to a maximum at low temperatures, 500 K or 600 K, depending on the pressure and then decreases considerably. The non-monotonic behavior originates from the competition between the two opposite contributions of temperature. Temperature facilitates overcoming the energy barrier for proton transfer. However, as temperature increases, the connectivity of water molecules decreases, and since proton transfer takes place along a hydrogen bond, 81-84 the hopping becomes less probable. Hence, there is a maximum at the highest possible temperature for which the hydrogen bond network is percolated. At higher temperatures, the excess proton can be trapped into an isolated water molecule or a patch of water molecules, decreasing its hopping chances. It is interesting to note that the highest proton hopping rate takes place at the conditions for maximum self-dissociation, 600 K and 100 MPa, 78 which could be expected due to the similar nature of both processes. 82,83 To finish, we must mention that these proton transfer rates must be necessarily qualitative since proton transfer has a small contribution from tunneling, which cannot be captured with empirical models. 85

IV . CONCLUSIONS
The ability of ReaxFF to reproduce chemical reactions makes it a promising tool to investigate sub-and super-critical (SC) synthesis and oxidation processes, yet the validity of the force field under high temperature and pressure conditions had to be evaluated. In this work, we have tested the description of SC water properties by the new parametrisation of water in ReaxFF. 44 We computed some of the most relevant properties that make SC water an important solvent from a technological point of view and compare the results with the available experimental data. The results can be summarised in the following items: • Density: ReaxFF describes accurately the density, with relative errors lower than 5% in the whole range of studied temperatures and pressures. A bad description of density was arguably the most problematic issue of the first generation of water ReaxFF parameters, and the second generation clearly solves it. • Dipole moment: The dipole moment, µ, is solely influenced by the changes in the atomic charges since the internal structure of the molecules remains almost constant in SC conditions. It increases with the density, from µ = 1.77 D at low densities, in close agreement with experiments, to µ = 2.08 D at high densities, which is smaller than the expected value. • Static dielectric constant: The static dielectric constant ( ) computed with ReaxFF decreases toward zero as density decreases. The trend is in qualitative agreement with experimental observations, yet the results are apparently little sensitive to temperature, while experiments predict, for a given density, lower as temperature increases. • Microstructure: ReaxFF reproduces the microstructure of water as depicted by neutron diffraction scattering experiments. The exclusion volume, distance, and intensity of the first coordination shell are captured in excellent agreement and just minor differences can be observed at longer distances. • Hydrogen bond network: We use a simple geometrybased criterion to count the number of hydrogen bonds per water molecule in the SC regime. The results are in very good agreement between experimental values and other models. Furthermore, ReaxFF reproduces the percolation limit found in previous simulations. • Self-diffusion: The self-diffusion coefficient of water is captured correctly within the typical experimental error. The diffusion coefficient spans over 3 orders of magnitude within the studied range, increasing asymptotically as density goes to zero. • Proton transfer: The general picture of proton transfer is consistent with high-level ab initio simulations and experiments. The hopping rate of an excess proton increases with temperature and pressure to a maximum close to the critical point and then reduces due to the lower connectivity between water molecules.
Overall, we can conclude that the description of SC water structure and dynamics by ReaxFF is satisfactory, being quantitative in most cases. This opens the possibility of molecular scale simulations of chemical reactions in the SC regime with water as a solvent, which can help in the understanding and design of SC technological processes.