Density isobar of water and melting temperature of ice: Assessing common density functionals

Water is ubiquitous on our planet, and life on Earth is unthinkable without it. Therefore, it is not surprising that water has received a great deal of attention from the scientific community. It was also one of the first liquids to be intensively studied by computer simulations.^1,2 At present, carefully constructed force fields play a key role in the study of water.^3–33 The best force fields are constructed to give an excellent description of all relevant properties of water (density, viscosity, interfacial free energy, phase transitions, etc.) while also giving outstanding agreement with high-level quantum chemistry calculations for small clusters. More recently, there has been a gradual shift toward supplementing or even substituting force fields with first-principles electronic structure calculations. These calculations predominantly rely on density functional theory (DFT), with the initial studies emerging in the early 1990s.^34–37 As computing power advanced, the calculations evolved to encompass more meaningful timescales and systems.^38–51 The advent of machine learning potentials has led to a further leap in the number of properties that can be predicted and the timescales that can be simulated.^52–65 In particular, notable is the capability to simulate density isobars and melting temperatures within a matter of days or even hours, rendering the comparison of these properties with experimental data straightforward.

Indeed, there is little doubt that the density isobars predicted by density functional theory are not yet on par with those obtained from carefully constructed semi-empirical models.^66–68 So why bother with the costly density functional theory calculations? The explanation is straightforward: as our microscopic understanding of water has progressed over time, the scientific community has redirected its focus toward more intricate questions. These encompass the interactions between water and surfaces, the nucleation of water at nanoparticles, and processes such as the electrolysis-induced hydrogen production or the hydrogen-to-water conversion in fuel cells.⁶⁵ Understanding the interface between water and surfaces, frequently metallic, is imperative in these cases. In particular, the simulations of electrolysis and fuel cells, where bonds are continuously broken and reformed, necessitate the application of first-principles calculations. However, conducting high-level quantum chemical calculations, for which some force fields have been calibrated, remains challenging for surfaces or practically impossible for metallic ones. Consequently, density functional theory is the go-to approach. Employing density functional theory potentially adversely impacts the properties of water, as already discussed above. So ideally, one seeks a density functional that accurately describes the liquid state of water, encompassing at least the density and melting temperature. It is, however, understood that this alone will not suffice to describe other properties such as the solvation energy of ions, transport properties of protons or ions, or interfacial energies, but simple thermodynamic properties are a first stepping stone. Overall, identifying a density functional that reasonably describes the density isobars and melting point of water is a pressing scientific issue that will strongly impact future water studies. This is the core goal of the present work, where we focus on testing commonly applied density functionals via machine learning potentials. Ionic quantum effects are neglected in the present study but are expected to be small (1% and 2% for the density and temperatures, respectively).⁶⁰ 

Our methodology closely follows our previous work,⁶⁹ where we have constructed machine-learning potentials (MLPs) for water using the RPBE-D3/zd density functional. The adopted strategy uses first on-the-fly training along the isobar using NPT simulations. Simulations cover the melting of ice and water in a temperature range between the supercooled state and slightly below the evaporation point. Further structures are added from parallel-tempering simulations. The final dataset consists of about 1500 structures with typically 64 molecules and corresponds to dataset 1 in Ref. 69. This will be referred to as the “base” dataset in the following. To obtain a highly accurate database, we have opted to use the hardest and most accurate potentials available for VASP^70,71 (O_h and H_h). Following Ref. 69, the Brillouin zone was sampled at the Γ-point only, and the final datasets correspond to a 2000 eV energy cutoff to avoid errors in the pressure. For the present work, we use the MLP implementation available in VASP. The on-the-fly training relies on Gaussian process regression using Bayesian regression.^72–74 For the production runs, we rely on kernel-ridge regression to solve the least squares problem using a singular value decomposition.^75,76

In a first attempt, we re-calculated the entire dataset 1 from Ref. 69 using all considered semi-local functionals at 2000 eV. We expect this strategy to work well for functionals that yield a volume within 5% of the base functional (RPBE-D3/zd), as the volume fluctuations during training are typically 5%. For SCAN (very different melting temperature and density) and optPBE (significantly larger density), we additionally performed a complete on-the-fly training using the same strategy as outlined in Ref. 69. In Table I, this strategy is named “on-the-fly.” For the hybrid functional revPBE0, we first upgraded a smaller on-the-fly RPBE-D3/zd database containing 662 structures (dataset 2 in Ref. 69) to revPBE0-D3/zd using Δ-learning, then continued training starting from this force field and refining it using on-the-fly training. This on-the-fly training was done using the hybrid functional, an energy cutoff of 1100 eV, and cooling liquid water twice from 350 to 320 K, respectively, to 250 K, as well as melting cubic ice. The final dataset comprises 1010 structures and covers a much larger set of volumes than the original database. We term the structures in this database “hybrid” in the following.

To obtain the final MLPs for hybrid functionals, baseline RPBE-D3/zd calculations were conducted at 2000 eV for all structures in the original dataset (either the base or the hybrid dataset). Subsequently, Δ-learning was employed to “upgrade” this dataset to different hybrid functionals. This was achieved by machine learning the energy, force, and stress differences between the desired hybrid functional and the RPBE-D3/zd functional at 1100 eV for 300–600 structures. When learning the difference between the two functionals, we found that using only 5 and 64 structures resulted in test set errors of 6 and 2.5 meV/Å, respectively. This does not add a significant error to the baseline error of about 26 meV/Å for semi-local functionals (error propagation follows the Gaussian law). The Δ-force field obtained in this manner is calculated for all structures in the original RPBE-D3/zd dataset, and the predicted energy, force, and stress differences are added to the baseline dataset. This upgraded dataset is then passed to the machine-learning algorithm and subsequently used in the production runs. Δ-learning was used to upgrade the “base” dataset for revPBE0-D3/BJ and the “hybrid” dataset for all hybrid functionals.

The Δ-learning approach taken here does not introduce sizable additional errors. Indeed, comparison with a test set calculated at 2000 eV using the hybrid functional gives a force error of 28 meV/Å, which is only slightly larger than the errors we found for semi-local density functionals (26 meV/Å). The main source for the slight error increase is that we added the difference between the hybrid functional and RPBE-D3/zd evaluated at 1100 eV instead of 2000 eV.

Nevertheless, it is important to be cautious when employing Δ-learning. The primary sources of errors stem not from capturing the difference between the two functionals but rather from potential gaps in the training dataset, which may fail to cover the necessary range of configurations. Consequently, predictions for the baseline semi-local functional can become relatively inaccurate, in the present case, if there is a significant difference between the volumes predicted by the hybrid and semi-local functional. The hybrid dataset remedies this issue since it spans a volume range covering ±10% around the experimental density.

Figure 1 shows the density–temperature isobars obtained using a Langevin thermostat and a Parrinello–Rahman barostat (friction coefficients in both cases being 2 ps⁻¹). Parallel tempering is performed for 128 molecules (for further details, we refer to Ref. 69). As shown in Ref. 88, parallel tempering notably improves the sampling of liquid water densities. In most simulations, we see that images at low temperatures tend to stop exchanging with images at higher temperatures. Although the structure of ordered ice is not recovered, the density around 210 K is always very close to the density of hexagonal ice reported in Table II. This can also be seen in Fig. 2, where we show the difference in density between ice and the parallel tempering runs along the isobar. Therefore, these structures are best categorized as low-density amorphous ice.^89–91 

The first important observation is that none of the functionals gives a satisfactory density for liquid water in the full range. At standard conditions, it should be around 1 g/cm³, but the best-performing functionals are about 4% away from experiments, only SCAN being in good agreement. However, such an agreement only holds at that particular thermodynamic state where it accidentally crosses the experimental isobar. In addition, the experimental density difference between liquid water and ice is not well reproduced; experimentally, ice is about 90% less dense than water.

We start the discussion with the semi-local revised Perdew–Burke–Ernzerhof (RPBE⁷⁷ and revPBE⁸⁰) functionals, all of which consistently give too low densities (8% lower than experiments, in the best case). The results also depend very much on the choice of damping, i.e., how the attractive van der Waals (vdW) terms are “damped” to zero as the constituents approach each other. The damping is necessary because the approximate semi-local DFT functionals already account for some of the interactions associated with correlation effects. Observing such a strong dependence on the choice of damping is disturbing. In particular, only the zero damping (zd), i.e., the original functional form proposed by Grimme et al.,⁷⁸ yields a higher density for water than for ice, while the Becke–Johnson damping (BJ)⁷⁹ incorrectly predicts ice to be denser than water so that the slope of the melting line is likely to be inverted.

The baseline DFT functionals also have some influence on the results. The RPBE and revPBE functionals are closely related. The revPBE functional is a revised version of the original PBE functional. It includes modifications to the exchange part of the functional aimed at improving its accuracy in describing properties such as the reaction energies of molecules. RPBE was primarily aimed at improving adsorption energies at surfaces and, as opposed to revPBE, locally observes the Lieb–Oxford criterion.⁷⁷ Overall, revPBE-D3/zd gives the best agreement with the experimental density, although it underestimates the density of ice and, even more so, water. This underestimation is also true for the difference in densities between water and ice, as illustrated in Fig. 2. The ice melting points are presented in Table II. They were determined by performing direct coexistence simulations in anisotropic NpT calculations on 864 molecules (half being ice and half being liquid, divided in two rectangular prisms separated by two planar interfaces). Pressure is kept at standard conditions, and various temperatures are covered using a 2 K grid. The temperatures reported are usually the average of the lowest temperature at which melting occurred and the highest temperature at which freezing occurred.

Considering the intrinsic error of direct coexistence simulations and our 2 K grid, we expect the errors to be of about ±2 K. We note that the agreement with our previously reported melting temperature for larger ensembles for RPBE-D3/zd⁶⁹ as well as with previous work⁵⁴ is excellent. Table II shows that most semi-local functionals overestimate the experimental value of 273.15 K. Among these, the -zd ones align reasonably well (5 and 10 K higher for RPBE-D3/zd and revPBE-D3/zd, respectively), whereas the -BJ versions notably deviate (37 K higher in the best case), even exceeding the temperature of maximum density.

Now we discuss the results for hybrid functionals as well as optPBE. For optPBE and revPBE0-D3/BJ (“base”), we illustrate results for the original database (open symbols), as well as on-the-fly results for optPBE and results using the hybrid database covering a larger density variation for revPBE0-D3/BJ (filled symbols). The predicted densities are about 1% smaller using the databases constructed from RPBE-D3/zd structures, and the error increases with the density difference to that density functional. The error becomes even larger for revPBE0-D3/zd (not shown). This is a sign of the somewhat limited extrapolation capabilities of the MLPs applied in the present work, but we note that despite a density difference of 10%–20%, density errors remain smaller than 1% for optPBE and revPBE0-D3/BJ.

Let us begin our more detailed discussion with optPBE, a functional that introduces vdW interactions through terms dependent on the density at two distinct points in space. This departs from conventional semi-local density functionals, which typically only consider semi-local information at a single point in space. As already noted above, the isobars for the two training sets are very similar. Both simulations lack a clear density maximum. Instead, the density increases almost linearly with decreasing temperature, with an abrupt jump to a low-density amorphous ice phase around 210 K. This result is reproducible, as confirmed by the two independent MLPs and, thus, independent of database details. A very similar unsatisfactory result is observed for revPBE0-D3/zd. Again, there is no discernible density maximum, and an abrupt transition to low-density amorphous ice structure occurs around 210 K. Both functionals, optPBE and revPBE0-D3/zd, overestimate the density difference between water and ice, which appears to be the primary cause of the lack of a discernible density maximum. In addition, their melting points occur at a very low temperature (around 240 K). On the other hand, the revPBE0-D3/BJ functional predicts densities that are about 5% lower than in experiments, somewhat similar but slightly better than for the semi-local functionals with vdW corrections. Similar to revPBE0-D3/zd, PBE0-D3/zd yields a larger density than in experiments, but now a density maximum is discernible in the isobar, whose shape agrees quite well with experiments although shifted about 30 K to higher temperatures. Its difference in density between water and ice gives the best agreement with experiments but shifted by about 20 K to higher temperatures. For PBE0-D3/BJ, the density isobar is very similar in shape to revPBE0-D3/BJ but shifted to higher densities and temperatures. In fact, the melting points for all the functionals under consideration are, in general, significantly too high (by ∼30 K), except for optPBE and revPBE0-D3/zd, which are too low (by ∼30 K).

It is worth noting that there are certain similarities between semi-local and hybrid functionals regarding the choice of the damping term. In particular, the use of Becke–Johnson (BJ) damping tends to increase the density of ice compared to water. While this adversely affects the performance of semi-local functionals, it has a somewhat positive effect on hybrid functionals, restoring the density maximum. However, now the density difference between ice and water is again too small, so the agreement with the experiment is still far from good.

Empirically mixing the machine-learning potentials for revPBE0-D3 with zero damping and Becke–Johnson damping yields the most favorable outcome (revPBE0-D3/mix). To obtain the revPBE0-mix results, we averaged the energies, forces, and stresses of the revPBE0-D3/BJ and revPBE0-D3/zd datasets with equal weighting and retrained a force-field on that dataset. This is equivalent to averaging the damping factor before the D3 dispersion corrections are multiplied by the damping function and added to the DFT energies and forces. The approach demonstrates excellent agreement in predicting the density of water, as well as its difference with ice, and yields a well-defined density maximum, although at temperatures that are ∼10 K lower than the experimental ones (with the maximum in density occurring slightly below the melting point).

Finally, we comment on the results using the SCAN functional. We observe a difference in the water and ice density, as well as a water density isobar that closely matches the experimental one but shifted to higher temperatures by about 50 K. Overall, the density of water is about 5% too large, but that seems acceptable considering the larger errors of other functionals. The density maximum is also at a temperature (323 K) that is significantly too high, whereas the melting temperature is too low compared to the temperature of the density maximum but still higher than the experimental value (by ∼ 30). Hence, if one desires to simulate liquid water using SCAN, the temperature should be set well above 330 K.

The density isobars and melting temperatures have been calculated before for several functionals, including RPBE-D3/zd, SCAN, and revPBE0-D3/zd. For RPBE-D3/zd, we have already documented good agreement (density and melting temperatures within 1% and 2%, respectively) with FHI-AIMS⁵⁴ in our previous work.⁶⁹ We also summarize these results in Table II and show the isobars in the supplementary material.

For SCAN, about 10 K higher melting temperatures and temperatures of the density maxima⁹² and 1.5% smaller densities have been reported,⁹⁵ overall in good agreement with the present work. The discrepancies might arise from a different plane wave energy cutoff (1500 eV in Refs. 92 and 95 and 2000 eV here), which we found quite relevant for water simulations.⁶⁹ 

For revPBE0-D3/zd, the disagreement with previous work^60,93 is more troublesome, but previous publications also differ significantly from each other. Most of the discrepancies seem to be a consequence of the higher water density in the present work.

To confirm the water densities predicted via the MLPs, we estimated the DFT volume using histogram reweighting for RPBE-D3/zd, SCAN, and revPBE0-D3/zd. We performed 1.5 ns simulations (1 mio time steps, H mass set to 8 a.u.) for 64 molecules at 273 K (RPBE-D3/zd and revPBE0-D3/zd) and 325 K (SCAN), and selected 1000 equally spaced configurations for the three cases. These were re-weighted using the energy difference between the desired functional and the MLP. For revPBE0-D3/zd, we chose an energy cutoff of 1100 eV; for RPBE-D3/zd and SCAN, 1100 and 2000 eV were used (both cutoffs yield identical results for the corrected volume). In all cases, the reweighted average DFT density is slightly smaller than the MLP density (0.1%–0.2%); however, the statistical uncertainties are about 0.5%. In particular, we note that the volume and the energy difference between DFT and the MLPs are not statistically significantly correlated. This evaluation confirms that the MLP-predicted volume is consistent with the VASP revPBE0-D3/zd ground truth.

We speculate that the shortcoming of previous calculations is related to specifics of the codes (see the supplementary material) or holes in the MLP database: when we set out from the initial RPBE-D3/zd database, the density isobars were similar to previous calculations, but as we added more high-density data from the on-the-fly learning, the isobars approached the present results. Furthermore, we show in the supplementary material that a cutoff of 6 Å on the exact exchange—as used in previous calculations—favors lower densities.

It is not easy to predict or disentangle the trends associated with different density functionals, but some relationships can be identified. It appears that the equilibrium density of water and ice is related to the strength of the two-body interactions predicted by the functional. In particular, it has been found that SCAN predicts a significantly too large two-body interaction,^96–98 while RPBE-D3 predicts too weak two-body interactions.⁹⁸ In line with these considerations, the BLYP-D3 functional predicts quite accurate two-body interactions and very good densities.^53,98

However, the shape of the density isobars and melting curves are difficult to relate to simple properties. The likely explanation is that the maximum in the density isobar relates to a subtle balance between three-body and two-body interactions, for which the present density functionals show a relatively large scatter. The only consistent observation that emerges from our study is that an accurate reproduction of the density difference between ice and water seems to be crucial for a reasonable isobar shape. Functionals that overestimate the density difference fail to produce a density maximum similar to the experiment. Instead, they predict an abrupt transition to a low-density amorphous ice phase (optPBE and revPBE0-D3/zd). Conversely, functionals with too small a density difference, or one in the right range, tend to produce more reasonable density isobar results, although potentially shifted in temperature and density. Notably, the SCAN functional falls into the latter category, with a density that is too high (5%) and a temperature deviation of almost 20%.

Regarding the damping of the vdW interactions, we can rationalize our observations as follows: VdW interactions only affect the two-body terms (3-body vdW terms were not considered in the present work). For water, the vdW corrections are typically 25% larger for zero damping than for Becke–Johnson damping for all functionals considered. This implies that the density is consistently greater for zero damping (compare Fig. 1). Furthermore, since ice is less dense and the inter-molecular distances are larger, vdW corrections tend to be smaller for ice-like structures. This means that ice is less affected by vdW corrections than water, especially for Becke–Johnson damping. In other words, increasing the strength of the vdW corrections (BJ → zd) also increases the density difference between water and ice, which can be advantageous or disadvantageous depending on the baseline functional.

In summary, we have assessed the density isobar of water and the melting temperature of ice for semi-local and hybrid functionals with vdW corrections across a range of commonly used density functionals. The outcomes, unfortunately, are not particularly pleasing, as none of the investigated functionals show good agreement with experimental density isobars. Although the shape of the isobars may resemble the experimental one for some functionals, a quantitative agreement would require a shift in density and temperature. It is challenging to determine whether the inadequacy lies in the vdW correction or the specifics of the density functionals, but both factors likely contribute. In addition, it is important to recognize that correlation effects are inherently non-separable, as they operate at all distances and wavelengths, and electronic excitations at one length-scale mix with those at other length scales. This makes it challenging to isolate vdW interactions from other contributions. Employing a damping function to truncate vdW interactions at a certain range is necessarily heuristic, and even parameter fitting to an extensive dataset does not ensure accurate predictions beyond the test set or for specific problems such as water and ice.

Among the semi-local functionals, we recommend either revPBE-D3/zd, as it provides the best overall density and accurate temperatures for the density maximum and melting, or RPBE-D3/zd, which is only slightly inferior for the density, although better in the melting point, both giving similar differences in the density between water and ice. However, a common issue with all semi-local functionals based on revised versions of PBE is an underestimation of the density of ice and water. This inadequacy is likely related to too weak two-body interactions.⁹⁸ 

In the case of hybrid functionals, there are often only small improvements in terms of densities and melting points, although previous publications report notably better performance for other crucial properties such as proton transport and reactivity.^99,100 Zero damping results are generally unsatisfactory, as the density difference between ice and water is too large, precluding the observation of an isobar density maximum resembling experimental results. Although BJ damping yields more reasonable results, they are not markedly superior to those using the semi-local functionals with zero damping for the considered properties. We only achieve excellent agreement with the experiment when combining an MLP with zero damping and BJ damping. We conclude that the state of density functionals for water is unsatisfactory, and more development effort will be required to resolve this issue. A simple fix suggested by the present work is to adjust the damping parameters in the vdW functionals specifically to water. However, this will adversely affect the performance of the vdW corrections for other systems. For revPBE0, we have empirically observed that averaging the standard (zero-damping) and Becke–Johnson damping functions yields a very reasonable density isobar.

In the supplementary material, the density isobars of this work are compared to those of previous studies, and we attempt to clarify discrepancies.

This research was funded in part by the Austrian Science Fund (FWF) through the SFB TACO No. 10.55776/F81 and in part by 10.55776/I6467. For open access purposes, the author has applied a CC BY public copyright license to any author-accepted manuscript version arising from this submission. Computer resources and technical assistance were provided by the Vienna Scientific Cluster (VSC).

The authors have no conflicts to disclose.

Pablo Montero de Hijes: Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Christoph Dellago: Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Writing – review & editing (equal). Ryosuke Jinnouchi: Methodology (equal); Software (equal); Writing – review & editing (equal). Georg Kresse: Conceptualization (equal); Data curation (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.