Density anomaly of water at negative pressures from first principles

Using molecular dynamics simulations based on ab initio trained high-dimensional neural network potentials, we study the equation of state of liquid water at negative pressures. From density isobars computed for various pressures down to p  =  −230 MPa we determine the line of density maxima for two potentials based on the BLYP and the RPBE functionals, respectively. In both cases, dispersion corrections are included to account for non-local long-range correlations that give rise to van der Waals forces. We have followed the density maximum down to negative pressures close to the spinodal instability. For both functionals, the temperature of maximum density increases with decreasing pressure under moderate stretching, but changes slope at MPa and MPa for BLYP and RPBE, respectively. Our calculations confirm qualitatively the retracing shape of the line of density maxima found for empirical water models, indicating that the spinodal line maintains a positive slope even at strongly negative pressures.


Introduction
In contrast to most other liquids, water expands on freezing. At atmospheric pressure, the density of water reaches a maximum at a temperature of 4 °C due to a reorganization of hydrogen bonds compensating for the thermal expansion usually occurring upon heating [1][2][3]. With increasing pres sure, the temperature of maximum density decreases, leading to a line of density maxima (LDM) with negative slope. In an attempt to explain this density anomaly and other anoma lies, such as the strong increases of the isobaric heat capacity and the isothermal compressibility in the supercooled regime, it has been proposed that a liquid-liquid transition exists in supercooled water between two liquids of different densities [4,5]. In this scenario, the apparent divergence of response functions is related to a second critical point terminating the liquid-liquid coexistence line. While this picture has been questioned [6][7][8][9], a recent analysis indicates that the discrep ancy at the origin of this debate might be due to issues with the simulations [10]. The second critical point picture was preceeded by the stability limit conjecture [11,12], in which the apparent singularities are attributed to a retracing spinodal instability line. Other scenarios to explain the anomalies of water include the singularity free scenario [13][14][15], in which pronounced maxima in the response functions are caused by the existence of the LDM, and weak crystallization theory [15,16], in which fluctuations close to the stability limit of the liquid are responsible for the behavior of the response functions. An overview of the various scenarios put forward as explanation for the anomalous properties of supercooled water can be found in [17] and [18].
The general shape of the LDM bears some significance for the behavior of the spinodal emanating from the vaporliquid critical point. If the LDM maintains a negative slope for negative pressures, the LDM is bound to intersect the liquidvapor spinodal as envisaged in the stability limit conjecture. It was noted by Speedy [11] that at the intersection point with the LDM, the spinodal line needs to have a horizontal tan gent implying a reentrant spinodal, which could trace back to positive pressures and give rise to the anomalies observed in the response functions. Later it has been pointed out [19] that the stability limit conjecture requires another liquidvapor critical point at low temperatures, somewhat reducing the interest in this scenario. Note, however, that if the insta bility is not towards the vapor phase, there is no need for an additional vapor-liquid critical point at low temperatures. Other scenarios for the explanation of the anomalies of water do not consider a reentrant spinodal such that in these cases the LDM must bend back and acquire a positive slope at suf ficiently negative pressures in order to avoid intersecting the spinodal [20]. Since the LDM allows to distinguish between different scenarios, efforts have been made to determine the LDM of water deep into the negative pressure regime.
While the experimental determination of the equation of state of water at positive pressures is relatively straightfor ward and the LDM has been located over a wide pressure range [21][22][23], experiments at negative pressures are more challenging because the stretched liquid is only metastable and easily decays via cavitation [24,25]. However, if hetero geneous cavitation is avoided, water can sustain negative pres sures in excess of −140 MPa due to the strong cohesive forces acting between the molecules [26]. In recent experiments, Caupin and collaborators [27,28] have determined the equa tion of state of water down to −137 MPa by measuring the sound velocity with Brillouin light scattering in a microscopic Berthelot tube [24,26]. The LDM obtained from their data bends to more negative slopes as the pressure is decreased. However, the LDM does not reach a turning point in the pres sure range studied in the experiments.
Complementary information on the behavior of the LDM is provided by molecular computer simulations, which can be carried out at even lower pressures than the experiments because cavitation can easily be avoided by simulating sys tems at short length and time scales [4,[29][30][31][32][33][34]. Besides allowing to study the anomalies of water at more extreme pressures than in experiments, molecular computer simula tions yield insights into their microscopic origin. Also, com parison of the LDM for a particular water model with the experimental LDM over a wide pressure range is an important test for the quality of the model. The LDM has been deter mined for various simple water models such as ST2 [4,29], SPC/E [30,33] and TIP4P/2005 [30][31][32], TIP5P [33], and mW [34]. The LDMs computed for these water models differ considerably, but in all cases they show a similar behavior in the (T, p)plane. Starting from high positive pressures and low temperatures the negative slope increases until the LDMs ultimately bend back towards lower temperatures at suffi ciently negative pressures in contrast to the perspective of the stability limit conjecture.
Since empirical water models are typically parametrized to reproduce experimental data determined for positive pres sures, one might question their applicability to the negative pressure regime. In principle, ab initio simulations offer the possibility to carry out atomistic simulations of liquid water over a wide range of pressures and temperatures in a predic tive and unbiased way without suffering from this limitation [16,35]. In practice, however, the high computational cost of such calculations severely limits accessible system sizes and simulation times, making a comprehensive study of the den sity anomaly of water from first principles impractical. High dimensional artificial neural networks have been suggested as a solution to this limitation [36][37][38]. The general idea of this approach is to train a neural network (NN) to reproduce ener gies and forces learned from a high quality set of reference data, for instance obtained with density functional theory. The neural network, which is nothing else than a very flexible fit ting function, is then used to carry out molecular dynamics simulations at the same accuracy but at a much lower cost than by applying the first principle method directly. Recently, we have developed a neural network potential for ice and water based on different dispersion corrected functionals and have shown that this model reproduces several properties of liquid water including its melting point and density maximum at atmospheric pressure [3].
In this work, we use this neural network potential to deter mine the LDM of water from first principles over a wide range of pressures. As the training data set of the neural network potential included stretched configurations typical for nega tive pressures, it reproduces the energetics of water even for strongly negative pressures close to the spinodal stability limit. Calculations are carried out for two parameteriza tions of the neural network potential trained on energies and forces calculated with the BLYP [39,40] and the RPBE [41] functionals, respectively. As such functionals based on the generalized gradient approximation (GGA) fail to properly describe dispersion interactions, Grimme's D3 corrections are added to the training data [42]. As we have shown recently [3], such van der Waals corrections are strictly necessary to give hydrogen bonds just the right geometry and flexibility required for developing a density maximum. Our simulations demonstrate that the LDM sensitively depends on the under lying functional. While the computed LDMs deviate to some extent from the experimental LDM [21-23, 27, 28, 43], they both reach a turning point in agreement with the behavior observed from empirical water models.
The remainder of the article is organized as follows. In sec tion 2 we give a brief overview of the computational methods used in this work. In particular, we discuss the basic idea of the neural network approach and review the properties of the neural network potential parameterized to reproduce the prop erties of liquid water and ice. Results for the LDM are pre sented and discussed in section 3 and some conclusions are provided in section 4.

Neural network potential
As reported in detail recently [3], we have developed an atomistic neural network potential capable of describing the structural and dynamical properties of liquid water and ice across the phase diagram. In this approach, devised by Behler and Parrinello [36], the total potential energy of a system is written as a sum of local energies depending on the chemical environ ment up to a cutoff radius around individual atoms. The energy associated with each atom is expressed by an arti ficial neural network that takes socalled symmetry functions as input. The manybody symmetry functions, which are struc tural fingerprints characterizing the chemical environment of each atom, depend on the distances and angles between the atoms within a certain cutoff (here, a cutoff radius of 6.35 Å is used) and they are constructed to be invariant with respect to global rotations and translations as well as the exchange of atoms of the same species. Note that electrostatic interactions are reproduced only locally and longrange electrostatic inter actions are neglected in the present potential.
For water, two separate neural networks are needed, one for the energy contributions of hydrogen atoms and one for contributions of oxygen atoms. No information about bonding or connectivity of the atoms is required by the neural network, and a specification of the atomic positions is sufficient. In fact, neural networks of the type used in this work are capable of describing the formation and breaking of chemical bonds such that they can be used for the simulation of chemical reactivity.
Our neural network potential for water and ice is based on feedforward neural networks with 2 hidden layers containing 25 nodes each. There are 27 and 30 nodes in the input layer of the neural network for the hydrogen and oxygen atoms, respectively, corresponding to the number of symmetry func tions used to characterize the atomic environment of these two elements. The information entering the input layer is processed by the neural network and the atomic energy contribution is returned via a single output neuron, following the common scheme for data processing in feedforward neural networks (see e.g. [37] for a detailed discussion). As activation function, the hyperbolic tangent function is used here for the hidden layers and the linear function for the output layer to allow for large variations in energy. Hence, there are 1376 parameters in the neural network for hydrogen and 1451 parameters in that of oxygen, leading to a total of 2827 parameters for the entire neural network potential.
These parameters were adjusted in the training phase in order to minimize the deviation between energies and forces predicted by the neural network and a set of reference data computed with density functional theory for a large set of configurations taken from the liquid phase and various solid phases of ice (ice Ih, XI, IX, II, XIV, XV, VIII, and X). Initial configurations for the training set were generated using molecular dynamics with the flexible nonpolarizable SPC/ Fw water model [44] and by artificially distorting the perfect lattice structures of the ice polymorphs. More configurations for the training set were generated with a first version of the neural network potential, which was then further refined using the larger training set. The final training set contained about 7000 configurations. As described in [3], separate neural net works were trained to reproduce energies computed with two popular GGA density functionals: BLYP [39,40] and RPBE [41], each with and without the D3 dispersion corrections of Grimme [42], leading to four different neural network poten tials for liquid water and ice. In the present work, we only use the potentials including dispersion corrections, which we denote as BLYPD3 and RPBED3. The predictions of the network were validated using test data not included in the training set. The root mean square error of these neural net work potentials is about 2 meV per water molecule for the energies and about 70 meV Å −1 for the forces. While energy (or force) calculations performed with the neural network potential are far less costly than ab initio simulations, they are still considerably more expensive than empirical water models based on a small number of LennardJones interac tions centers and Coulombic charges, such as TIP4P or SPC.
As has been known for a while, BLYP water without dis persion corrections is overstructured and displays glassy behavior with a very low water mobility [16,45]. Including van der Waals corrections softens the liquid and yields to a dif fusivity and viscosity much closer to the experimental value [3]. For RPBE water the opposite effect occurs. While without dispersion corrections the liquid is too mobile, inclusion of these corrections reduces the diffusion constant and viscosity leading to good agreement with experiment. For both func tionals, dispersion corrections also have a strong effect on the density of the liquid. In both cases, the density of liquid water at atmospheric pressure is strongly underestimated and no density maximum occurs. With dispersion corrections, on the other hand, a density maximum appears at approximately the right temperature for both models [3]. Dispersion correc tions also improve the estimates for the melting temperature (T m = 283 K for BLYPD3 and T m = 274 K for RPBED3) and yield an ice density that is lower than that of the liquid in contrast to the case without added van der Waals interactions [3]. Both effects are due to a hydrogen bonding network with just the right level of flexibility, which allows second shell water molecules to penetrate the first shell leading to the den sity maximum at about 4 °C.

Molecular dynamics simulations
In this work we carry out extended molecular dynamics (MD) simulations to determine the LDM at negative pres sures for the two neural networks trained to reproduce the properties of the BLYPD3 and RPBED3 functionals. All simulations are carried out using the molecular dynamics program LAMMPS [46] into which we have incorporated the neural network potential. To determine the temperature of the density maximum we performed MD simulations for various pressures and temperatures in the NpT ensemble for a system consisting of 360 water molecules using the equa tions of motion of Shinoda et al [47] integrated with a time step of 0.5 fs. The lengths of the simulations were between 2 ns and 37 ns depending on the temperature. In total, we car ried out 91 simulations for the BLYPD3 potential and 93 for the RPBED3 potential, corresponding to a total simulation time of approximately 1.2 μs, of which roughly 370 ns were used for equilibration. Temperatures were mostly chosen in 10 K intervals in the range from 230 to 360 K. First, we cal culated the density ρ(T, p) as a function of temperature T at pressures p = 0.1, −100 −160 MPa for both functionals and later further pressure values were added to follow the shape of the LDM. For the sake of simplicity, the simulations on the initial isobars were started from only two distinct water configurations (one for each of the functionals). These were obtained by melting ice structures within 2 ns followed by an equilibration period of 10 ns at 300 K and 0.1 MPa. Upon adding more isobars, initial configurations for new trajectories were chosen from already existing simulation data close by in the (T, p)plane. For the raw density time series we carried out an analysis of the autocorrelation function to determine statis tical errors with algorithms provided in the pymbar package [48][49][50]. From the resulting uncorrelated samples we calcu lated averages and error bars.

Fitting the equation of state
In an effort to create a comprehensive view on the equation of state for the two neural network water models, we followed the approach of Poole et al [29] and fitted a 2D model function to the densities ρ(T, p) obtained from the molecular dynamics simulations. Here, C i is the Chebyshev polynomial of the first kind of order i and the a ij are fitting parameters. The sum mation limits i max and j max define the highest occurring poly nomial order along the T and p axis, respectively. We found that (i max , j max ) = (4, 3) provides a sufficiently complex model function to fit the simulation data. Since the Chebyshev polynomials map arguments in the interval [ for the BLYPD3 neural network potential. Note that the model parameters a ij have units of kg m −3 .
Once the parameters a ij are optimized, the corresponding model LDM can be easily derived from (1) via which allows us to interpolate between simulation points in the (T, p)plane.

Spinodal line
While direct access to the spinodal line is not possible from our MD simulations, we estimate its location following an approach suggested by Poole et al [51] theoretically rooted in meanfield theory [52]. In this approach, the pressureversus density isotherms are fitted by a power law, where A is a constant and δ is a critical exponent. For each isotherm a separate fit is performed, resulting in predictions for the spinodal pressure p s (T) and density ρ s (T) as a func tion of T. In all our fits, we use the critical exponent δ = 2 predicted by mean field theory [52] and fit only the remaining parameters (i.e. A, p S and ρ s ) to the simulation data. Results of the fitting procedure are discussed in the next section.

Results and discussion
Density isobars obtained from the molecular dynamics simu lations for various densities are shown in figure 1. The sym bols with error bars denote simulation results while the colored solid lines are obtained from the fitted model function of (1) for the same pressures at which the simulations where car ried out. The solid black lines, which are also obtained from the fitted model function, connect the maxima of the density isobars. As can be inferred from the figure, the model function ρ m (T, p) reproduces the simulation results very well over the entire temperature and pressure ranges explored in the simu lations. Note that particularly for the BLYPD3 case the error bars of the density are larger for lower temperatures due to the slower decay of correlations. Both for the BLYPD3 and the RPBED3 functional, the density isobars show a pronounced density maximum down to the most negative pressures for which we have carried out simulations. In experiments it is very challenging to determine the equation of state of water at strongly negative pressures because under these conditions the liquid is only marginally  metastable and very prone to cavitation [27,28]. In simula tions, however, cavitation does not occur easily due to small systems sizes and short simulation times. For both models, the LDM bends back to lower temper atures for sufficiently negative pressures. This change of slope occurs at pressures of p ≈ −200 MPa and p ≈ −20 MPa for the BLYPD3 and RPBED3 functional, respectively. Note, however, that in the BLYPD3 case the effect is less pronounced and hardly noticeable in the simulation results directly. The backbending can be clearly identified in the LDM obtained from the fitted equation of state, which combines information from all simulations. The retracing character of the LDM is also clearly visible in figure 2, in which the density ρ m is shown as a colorcoded map in the (T, p)plane. Here, the boundaries between dif ferent color shades are isochores and the lines of maximum density, shown as solid black lines, connect the minima of these curves. Data for the figure have been obtained from the model function (1) fitted to the simulation data. Note that for a given temperature and pressure, the density for BLYPD3 water is typically higher than that of RPBED3 water by more than 10%.
Using the the fitting procedure described in section 2.4, we have estimated the location of the spinodal in the (T, p)plane. Results of such fits are shown in figure 3 for both the BLYPD3 and the RPBED3 functional. Note that for the data shown, the functional form of (5), suggested by mean field theory [52], fits the simulation results very well for an exponent of δ = 2. A reasonable fit, however, could not be achieved for all temperatures of figure 1 due to insufficient amount of sampling points or noisy data. The spinodal pres sures p s (T) obtained from the fits are shown in figure 4 for both models.
In figure 4, the LDMs calculated for the two NN water models are compared to LDMs obtained previously for various empirical water models including ST2, SPC/E, TIP4P/2005, TIP5P and mW water. Also shown in the figure is the LDM determined experimentally for positive [21][22][23] and negative pressures [27,28,43]. As in the case of the two models studied in this work, all empirical water models yield retracing LDMs. Their specific position in the (T, p)plane is, however, quite different. The experimental results available at this point, shown in figure 4 as filled black triangles and circles, show a downward curvature hinting at retracing behavior. However, the data are not available at sufficiently negative pressures to answer this question conclusively. Of all empirical water models, TIP4P/2005 water reproduces the experimental LDM best. At atmospheric pressure, the temperature of maximum density of RPBED3 water is very close to the experimental value, but at lower and higher pressures considerable devia tions occur and at strongly negative pressures BLYPD3 water agrees more closely with the experimental behavior.

Conclusion
Using a high dimensional neural network trained with ab initio reference data for two popular density functionals (BLYPD3 and RPBED3) including dispersion corrections, we have determined the equation of state of liquid water under negative pressure, paying particular attention to the shape of the LDM. We find that the density maximum, deter mined for these two models previously at atmospheric pres sure [3], persists down to very negative pressures close to the spinodal instability. The specific location and shape of the LDM in the (T, p)plane differs in the two models, but in both cases the LDM bends back for sufficiently negative pressures, confirming findings obtained with several different empirical water models. Hence, our simulations, which are based on a description of liquid water derived from first principles without experimental information, indicate that a retracing LDM is a robust property that does not depend on model details, making it plausible to expect this behavior also for real water.
The substantial discrepancies between the LDMs obtained for the two ab initio based models and the experimental curve, both at negative and positive pressures, point to the deficien cies of GGA density functionals in capturing the properties of water even if dispersion corrections are included. An improved description of liquid water and ice, including an accurate pre diction of the equation of state and the LDM at negative pres sures, may be achieved by training highdimensional neural network potentials with energies and forces determined with first principle approaches beyond density functional theory, such as Møller-Plesset perturbation theory [53], coupled cluster theory [54] or the random phase approximation [55]. Lines of density maxima for the BLYPD3 functional (red line) and for the RPBED3 functional (blue line) compared with LDMs of several popular empirical water models (ST2 [29], SPC/E [30], TIP4P/2005 [31], TIP5P [33], and mW [34]) taken from the literature (grey lines). The single symbols mark the temperatures of maximum density obtained from simulations with the NN potential for the BLYPD3 functional (empty red circle) and the RPBED3 functional (empty blue circle) as reported previously [3]. Also shown are experimental data for negative pressures (black circles) [27,28,43] and for positive pressures (black triangles) [21][22][23]. The position of the spinodal lines is indicated by blue (RPBED3) and red (BLYPD3) squares connected via dashed lines.