Mass density fluctuations in quantum and classical descriptions of liquid water

Understanding the benefits and limitations of ab initio approaches based on quantum density functional theory (DFT) for describing aqueous phase processes in bulk and in the vicinity of interfaces continues to be an active area of research. Many studies regarding the accuracy of DFT to describe both bulk and interfacial properties of neat water have been performed and focus on how the role of simulation protocol affects the computable observables.^1–25 More recently, the efficacy of DFT-based methods to describe water in environments ranging from the gas phase to the condensed phase has been called into question.^26,27 One solution to the problem is to use a sophisticated classical empirical interaction potential based on a fit to the energetics of configurations obtained with high-level wavefunction methods.^28–31 These recent studies have produced excellent agreement with structural and spectroscopic properties of water and are designed to be correctly coupled with path integral calculations to explore the role of nuclear quantum effects (NQEs).

The advantages of using an empirical representation interaction over DFT-based methods are clear from the point of efficiency. Until recently, the phase behaviour of DFT-based methods has been informed by relatively short simulation times and small system sizes.^1–3 The results of these studies produced interesting results pertaining to the melting points and boiling points of popular DFT functionals.^1–3 It should be noted that earlier studies of these thermodynamic properties were performed with exchange-correlation (XC) functionals that did not contain a correction for long-range dispersion interactions that are absent from DFT.¹ 

The recent empirical corrections due to Grimme^32,33 have greatly enhanced the agreement with experiment over a range of structural, dynamical, and thermodynamic properties.^2,4–9 One of the most important thermodynamic properties of DFT water that was markedly improved was the mass density at ambient conditions.^4,6–8 This improvement in the mass density using the empirical corrections to the dispersion interaction has allowed for rapid progress to be made in the understanding of ions and reactivity in the vicinity of the air-water interface.³⁴ 

In similar spirit to the fitting empirical potentials to high-level wavefunction methods discussed above, empirical interaction potentials using DFT-based levels of electronic structure with and without dispersion have been constructed.¹⁰ These potentials have afforded the opportunity to perform simulations for relevant time scales and system sizes to converge the properties of bulk liquid water.¹⁰ The results of this study further corroborate some past careful studies using DFT interaction potentials^2,4–9 and clear up many inconsistencies regarding the thermodynamic properties of DFT water. This aforementioned study also suggests a picture where the revised functional of Perdew-Burke-Ernzerhof (revPBE)³⁵ in conjunction with Grimme’s third generation of dispersion (D3) produces an effective description of liquid water over a range of condensed phase environments.

One reason to consider an alternative to parameterized empirical potentials is to understand processes that involve the response of liquid water to a notional interface. Here, we desire to exploit the flexibility of DFT-based interaction potentials to correctly describe the short-range response to an arbitrary perturbation from the bulk liquid, namely, solutes or macroscopic interfaces. To this end, the short-range response to hard sphere cavities of various sizes obtained with DFT was directly compared to two popular fixed charged empirical potentials.¹¹ This study suggests that the quantitative differences observed in the short-range response between different water models lead to questions about the quality of interaction potentials needed to obtain solvation free energies of ions. Indeed, an earlier study on the local structure of ions as determined by the extended x-ray absorption fine structure (EXAFS) technique found that DFT-based interaction potentials were required in order to reproduce the accurately measured short-range structure.³⁶ Although progress is being made toward high-quality empirical force fields for ions based on fits to high-level wavefunction methods,^37,38 to the extent that empirical potentials can reproduce the details of local solvent response to interfaces remains important research.

The main focus of previously detailed studies of DFT-based methods has been on the equilibrium structural and dynamical properties.¹⁰ Herein, we compare and contrast empirical potentials against DFT for phenomena that are germane to computing solvation free energies, namely, mass density fluctuations. The choice of empirical potentials for this study is the SPC/E and MB-pol models of water. The former is chosen because of both its popularity and use in the study of hydrophobicity; the latter is chosen because of its demonstrated accuracy in producing the correct potential energy surfaces as benchmarked by high-level wavefunction methods. This will require that we establish the DFT simulation protocol to quantify the role of mass density fluctuations under both isothermal (NVT) and isobaric (NpT) ensembles for system sizes that are relevant to DFT studies. The importance of capturing the mass density fluctuations at short and long length scales forms the corner stone of the theory of hydrophobicity and solvation. Furthermore, the examination of fluctuations provides an additional self-consistent check on the thermodynamic properties of surface tension and isothermal compressibility.^39,40 Going beyond traditional probes of aqueous structure, we contrast the local structure of ambient water by examining the distribution of the 5th nearest neighbor distance (d₅). This order parameter was found to be relevant for describing the experimental structure of water under pressure and possibly a diagnostic for providing signatures of differences between empirical and DFT models of liquid water.⁴¹ The goal of this study is to provide a clear comparison of mass density fluctuations between different representations of interaction. This will require the development of DFT simulation protocol that provides a robust and consistent picture of structure and their fluctuations, thus further advancing our understanding of the utility of using quantum descriptions of interaction based on DFT to inform our understanding of complex phenomena in the condensed phase.

All the simulations presented here have been carried out using the CP2K program within the Born-Oppenheimer approximation, i.e., the wavefunction was optimized to the ground state at each time step. The QUICKSTEP module within CP2K was used to employ the Gaussian and plane wave (GPW) method.^42,43 In this GPW method, both the Gaussian and plane wave basis sets are used to linearly expand molecular orbitals and electronic density, respectively. Our model system consisted of 64 water molecules in a cubic simulation box under periodic boundary conditions. All NpT simulations were carried out at the ambient thermodynamic conditions, namely, the temperature was set to 300 K and the pressure was set to 1 bar using the reversible algorithm due to Tuckerman and co-workers.⁴⁴ The time step was maintained to be 0.5 fs. Nose-Hoover thermostats were employed to all degrees of freedom using the “massive” thermostatting. The time constant of the thermostat and the barostat was set to be 11.12 fs (corresponding to 3000 cm⁻¹) and 300 fs, respectively. All the NVT simulations were carried out using 256 water molecules in a cubic box of side length of 19.7319 Å providing a density of 0.997 g/cm³ at a temperature of 300 K. Both revPBE³⁵ and BLYP^45,46 functionals were used with the Grimme dispersion correction^32,33 known as D3 and D2, respectively. The core electrons were replaced by the norm-conserving pseudopotentials of Goedecker and co-workers⁴⁷ to carry out the simulations efficiently. Two types of basis set were used: a triple-$ζ$ valence Gaussian basis set augmented with two sets of d-type or p-type polarization functions (TZV2P) and the molecularly optimized double-$ζ$ basis set (MOLOPT-DZVP-SR-GTH which we will refer to as MOLOPT in the remaining text).⁴⁸ Both of these basis sets were previously successfully used with these functionals in the NVT simulations of bulk water at ambient, high pressure, and high temperature conditions.^41,49 In a NpT molecular dynamics run, longer simulation times are required to obtain the equilibration and to sample the fluctuation. We ran the simulations to produce a 100 ps long trajectory, from which the last 50 ps was used to gather statistics. The NpT dynamics were carried out using a larger reference simulation cell to ensure a constant number of grid points and provide a lower bound on the electron density cutoff. The reference simulation cell used was 19% larger than the original simulation cell and corresponding to the density of 0.59 g/cm³.

It has been established that the cutoff of 400 Ry for the expansion of electron density in the plane wave basis produces converged results in the NVT ensemble. However, in the case of NpT ensemble, a much larger cutoff is needed to produce the converged virial. It has been reported by McGrath et al.⁵⁰ that an NpT Monte Carlo simulation with a cutoff of 1200 Ry produced 10% lower density than that with a cutoff of 280 Ry. Another more recent NpT Monte Carlo simulation by Del Ben et al. used a cutoff of 800 Ry.⁷ They confirmed that changing cutoff from 800 to 1200 Ry did not affect the density. In the original NpT MD simulation, Schmidt et al. found that increasing cutoff from 600 to 1200 Ry did not change the density.⁶ 

However, there are many options for simulations in the NpT ensemble within CP2K and it is instructive to provide useful information regarding how simulation protocol can affect the outcome. A summary of these options in addition to convergence tests is detailed in  Appendix A. By using the standard Fourier interpolation technique, the total pressure (as defined by $13TrΠ$⁠, where $Π$ is defined in Ref. 6) was sufficiently converged to at a cutoff of 800 Ry to reproduce a mass density in agreement with previous studies (see  Appendix A).

Having established the simulation protocol used in this study, we can turn to the calculation of the mass density of DFT-based interaction potentials. The mass density can be calculated from an NpT run, using the aforementioned protocol, by taking an average of the instantaneous fluctuating volume over the simulation time. Figure 1 shows the variation of instantaneous mass density and the corresponding running average with simulation time for the 64 water box using revPBE-D3 functionals with TZV2P and MOLOPT basis sets. The calculated average value and the root mean square deviation are given in Table I. Our estimates provide a picture where the revPBE-D3 functional is providing a density of 0.962 g/cm³ and 0.988 g/cm³ with TZV2P and MOLOPT basis sets, respectively. These values are in good agreement with the experimental density of 0.997 g/cm³. The difference of 0.01-0.03 g/cm³ does not account for more than 1% in the lattice constant making up the simulation supercell.

Our calculated value of 0.96 g/cm³ is consistent with the previously reported value of 0.96 for revPBE-D2 by Lin et al.⁴ However, in the cited study, they did not calculate the density directly from an NpT ensemble. Instead, they used an indirect method where the total energy was calculated as a function of the scaled lattice constants for a given snapshot obtained with a Car-Parrinello molecular dynamics (CPMD) trajectory. The equilibrium mass density was obtained from the minimum of the interpolated energy. Our values are less than the previously reported value of 1.02 g/cm³ by Wang et al.¹² for revPBE using the nonlocal van der Waals (vdW) correlation functional proposed by Dion et al.¹³ However, these calculations are also not obtained from a traditional NpT ensemble. Rather, this study calculated the equilibrium density from the pressure-density curve obtained from NVT simulations at different volumes. To the best of our knowledge, our results report the first NpT simulations of revPBE-D3 water and its equilibrium density at ambient conditions. Like other popular gradient corrected (GGA) functionals (e.g., PBE and BLYP), in the case of revPBE-D3, the density has been significantly improved (from 0.69 to 0.96) towards the experimental value with the inclusion of dispersion correction (Grimme D3). This is consistent with the consensus that GGA functionals require the dispersion corrections to obtain a physically reasonable description of liquid water.

Figure 2(a) depicts the oxygen-oxygen radial distribution functions (RDFs) from our NpT simulations using both TZV2P and MOLOPT basis sets along with the experimental radial distribution functions previously published in Ref. 51 by Skinner et al. Our calculated RDF using the TZV2P basis set shows an excellent agreement with the experimental data. Most importantly, the position of the first peak is in the correct position, and the first minimum contains the correct amount of disorder as compared to the experiment. This suggests that revPBE-D3 water has the potential to display better diffusivity at 300 K as compared to other popular GGA functionals. Previous Monte Carlo simulations in the NpT ensemble for BLYP-D3, PBE0-ADMM-D3, and MP2 have predicted the first minimum in the RDF to be significantly more shallow than the experiment although a good mass density is reproduced.⁷ The diffusion constant for revPBE-D3 water has recently been calculated by Marsalek and Markland.⁵² An 800 ps of classical revPBE-D3 simulation provided a system size corrected diffusion coefficient of $2.22±0.05×10−9$ m² s⁻¹, within statistical error bars of the experimental value of $2.41±0.15×10−9m2s−1$⁠.⁵² The only significant deviation from the results herein is the height of the first peak that is higher by 0.2 when compared to the experiment. It should be noted that we did not include the nuclear quantum effect (NQE) into our simulations. It was previously reported that the inclusion of NQE might influence the height of the first peak towards the experimental value.^53–57 

RDFs calculated using the short-range molecular optimized basis set (MOLOPT)⁴⁸ at the double-$ζ$ level are also in good agreement with the experiment. Interestingly, MOLOPT produces the correct height of the first peak but is slightly shifted to larger distances than the experiment. Moreover, the first minimum suggests less structuring as compared to both TZV2P and the experimental results. To understand the origins of the difference between the two basis sets, we compared the RDF calculated with and without the D3 dispersion correction (see Fig. 15 in  Appendix C for the corresponding RDFs). Our calculations indicate that in the absence of dispersion correction, both the basis sets give similar RDFs; however, in the presence of D3 dispersion correction, MOLOPT results deviate from that of TZV2P. This indicates that the origin of the difference is due to the matching of the basis set with the Grimme dispersion correction scheme. Since the original D3 parameters were optimized with TZV2P basis sets, they do not work as well with the MOLOPT basis sets. A comparison of RDF with D3 dispersion correction to that without (see Fig. 15 in  Appendix C) clearly shows that the structure of water becomes softer in the presence of dispersion correction for both basis sets. This is consistent with a previously observed phenomenon that van der Waals interactions in simulations of water alter structure from mainly tetrahedral to high-density-like.^12,58

As another self-consistent check of our NpT protocol, we compare the RDF calculated from our NpT ensembles to those calculated from NVT ensembles. Theoretically, the NVT and NpT approaches should yield the same results if the protocol in both approaches is converged. Indeed, our RDFs from both NpT and NVT simulations are similar as shown in Fig. 2(b). This is a clear improvement in our understanding between the different approaches to simulation. Previous results for MP2 water using NpT Monte Carlo simulations provided a very different RDF from that obtained by simulations using other ensembles.^7,59

Additional comparisons between revPBE-D3 and BLYP-D2 were carried out in the NpT ensemble. BLYP-D2 has been a popular choice for numerous past studies of water and is known to produce satisfactory results regarding the mass density.^6,49,60 Our simulations suggest that the mass density obtained using BLYP-D2 at 300 K is 1.04 g/cm³ [see Fig. 3(a)]. This is slightly higher than the reported value of 0.992 g/cm³ (⁠$±$0.036) by Schmidt et al.⁶ This deviation might be attributed to the difference in the temperature (330 K used by Schmidt et al.). The slightly higher density obtained here with BLYP-D2 (1.04 g/cm³) is also comparable to the BLYP-D3 density reported by Del Ben et al. (1.066 g/cm³)⁷ and Ma et al. (1.07 g/cm³).⁸ 

Figure 3(b) compares the oxygen-oxygen RDFs for revPBE-D3 and BLYP-D2, both using TZV2P basis sets, along with the experimental RDF. The height of the first peak for BLYP-D2 is too pronounced and the first minimum is also significantly deeper when compared to the experimental results. Moreover, the position of the first peak is also slightly at a lower distance compared to the experiment. Overall, our research suggests that revPBE-D3 is producing a better overall mass density and liquid structure as determined by the experimental oxygen-oxygen RDF than BLYP-D2.

Going beyond the RDF to provide a more detailed description of the local structure of water is one way to differentiate between different representations of interaction for water. Understanding the local structure is crucial to understand the anomalous thermodynamic and kinetic behavior of water. Elucidating whether the local structure of water is just a random collection of states generated by hydrogen bond fluctuations or the competition between different specific locally favored structures in the free energy landscape could play a key role in advancing our understanding of the bulk homogeneous phase of water. To characterize the local structure of water, a wide range of different order parameters have been used.⁶¹ The most widely used order parameter is the tetrahedral order parameter (q) that is focused only on the first shell water. Previous studies have shown that order parameters that describe the second shell order may play a role in our understanding of the bulk homogeneous phase of water.⁶¹ In a recent work,⁴¹ it was observed that the behaviour of the 5th nearest neighbour water molecule is crucial to understand the change in local structure of water under pressure (from ambient to 360 MPa pressure). The so-called d₅ order parameter has been previously used to investigate the local structure of supercooled water.^62,63 Here, we have focused on the distance of the 5th water from the central water molecule as a suitable order parameter to analyze the local structure of water at ambient conditions.

To this end, we have analyzed the revPBE-D3, SPC/E, and MB-pol (in the NVT ensemble under bulk periodic boundary conditions in a supercell containing 256 water molecules) in terms of this d₅ order parameter. We have calculated d₅ as follows: For a water molecule i, we ordered all the other water molecules in the simulation according to the increasing radial distance to that water oxygen from the ith water oxygen (d_ji). Then the order parameter d₅ is simply the distance between the ith water oxygen and its 5th water oxygen (d_5i). Figure 4(a) displays the probability distribution of d₅ over all the water molecules in the simulation box for the revPBE-D3, MB-pol and SPC/E waters. The average values of d₅ for the revPBE-D3, MB-pol and SPC/E waters at ambient conditions are 3.49 Å, 3.36 Å, and 3.38 Å, respectively. Although all models show quantitative differences, SPC/E and revPBE-D3 are in better agreement than with MB-pol.

We can take this analysis a step further and consider the influence of the hydrogen bonding between the 5th water molecule to any of the four water molecules comprising first solvation shell. Specifically, we calculated d₅ for all water molecules in the simulation and divided them into two groups. The first group represents configurations in which the 5th water forms a hydrogen bond to any of the first shell waters. The second group represent the 5th that is not hydrogen bonded to molecules comprising the first shell. We use the standard hydrogen bond criteria of the distance between two oxygens being less than 3.5 Å and O–H–O angle being less than 30°. The results are shown in Figs. 4(b)–4(d) where the SPC/E and MB-pol models show distinct behavior when compared to revPBE-D3. However, when examining the hydrogen bonding distributions of the 5th water, MB-pol and revPBE-D3 seem to be in better qualitative agreement. The computed distances for the hydrogen bonded 5th water have the average value of 3.51, 3.43, and 3.40 Å and for the non-hydrogen bonded 5th water are 3.46, 3.33, and 3.36 Å for the DFT, MB-pol, and SPC/E water, respectively. Although all absolute distances are different, there seems to be the largest difference between hydrogen bonding and non-hydrogen bonding in the MB-pol representation of interaction.

Finally, we can look at the average value of d₅ for each individual water molecule. Figure 5 displays the distribution of the mean value of d₅ for each individual water for all three models. The salient point of this analysis is that it clearly demonstrates the inflexibility of the SPC/E water model predicting a very narrow distribution of the average d₅. This is an indication that every SPC/E water molecule has nearly the identical average environment. On the other hand, DFT water produces a wide distribution in d₅ suggesting that, on average, water explores a wide range of local environments even under bulk homogeneous conditions. The MB-pol model seems to capture this local heterogeneity and has a broader distribution than SPC/E but remains significantly more restricted than DFT. Capturing the flexibility of d₅ under ambient conditions seems to be an indication of the ability of a water model to describe the correct structure under different environments.⁴¹ To the extent that this is a relevant distribution under bulk homogeneous conditions at ambient conditions is yet to be determined experimentally. Nevertheless, this order parameter that is presented here is able to discern between different representations of interaction.

Now that the properties of two DFT-based interaction potentials have been established with respect to experimental radial distribution functions and mass density, we can further push our understanding of water, and we examine and compare the quality of the mass density fluctuations that are related to thermodynamic variables. To start, we examine the isothermal compressibility, $κT$ using the instantaneous volume fluctuations in the NpT ensemble using the following formula:

where V is the instantaneous volume of the system, k_B is the Boltzmann constant, and T is the simulation temperature. The volume fluctuations were averaged over a trajectory. Figure 6 depicts the variation of the isothermal compressibility as a function of simulation time. An average over the 50 ps trajectory provided a value of 42 Mbar⁻¹ for revPBE-D3 and 35 Mbar⁻¹ for BLYP-D2. The revPBE-D3 value, as expected, is very close to the experimental value of 45 Mbar⁻¹ and is in much better agreement than the previously computed value using PBE-D3 (21 Mbar⁻¹) and for PBE0-D3 (32 Mbar⁻¹) by Gaiduk et al.⁶⁴ Interestingly, both PBE-D3 and PBE0-D3 were reported to produce a density of water very close to the experimental density (1.02 g/cm³ and 0.96 g/cm³, respectively). revPBE-D3 water also produces a compressibility in better agreement with the experiment than vdW-DF, vdW-DF^PBE and VV10 levels of theory (18.2, 32.2 and 59.0 Mbar⁻¹, respectively) as computed by Corsetti et al.⁶⁵ It is interesting to note that among the four popular classical empirical water models (i.e., SPC/E, TIP4P, TIP4P/2005, and TIP5P), SPC/E and TIP4P/2005 are found by Helena et al. to produce the best agreement with the experimental compressibility at ambient conditions.⁶⁶ However, SPC/E fails to produce the change of compressibility with temperature. TIP4P/2005 is the only point charge model parameterized to produce the temperature dependence of the isothermal compressibility for liquid water.⁶⁶ The MB-pol model produced an isothermal compressibility of 45.9 Mbar⁻¹ at ambient conditions which is in excellent agreement with the experiment.⁶⁷ The MB-pol model was also found to produce the experimental pattern of changing isothermal compressibility with temperature.⁶⁷ The iAMOEBA model of flexible water that used three body terms for electronic polarizability was also found to provide an excellent pattern of changing isothermal compressibility with temperature.⁶⁸ As an additional self-consistent check to the estimates of the isothermal compressibility that were provided above, we examine the impact of both the thermostat and barostat frequencies on our results.

To this end, we have conducted separate simulations of 20 ps in length using a significantly higher time constant (i.e., lower frequency) for both, namely, 1 ps. Figure 7 shows the computed instantaneous mass density and the isothermal compressibility for these revised simulations. As shown in Fig. 7, the revPBE-D3 density is not significantly affected by the lower thermostat and barostat frequency. Interestingly, the BLYP-D2 mass density decreased from $∼$1.02 g/cm³ [see Fig. 3(a)] to $∼$0.99 g/cm³. Moreover, there seems to be a nontrivial but small dependence on the choice of thermostat frequency on the resulting mass density of DFT-based water. It should be pointed out that the reason for the higher barostat frequency was to be able to sample volume fluctuation over the significantly shorter simulation time afforded by DFT.⁶ A more reasonable barostat time constant of 1 ps is generally used in conjunction with classical empirical potentials. Not surprisingly, the isothermal compressibility displays a significant dependence on the barostat and thermostat values producing lower values for both functionals studied herein [see Figs. 6 and 7(b)]. To further our understanding, we examined the isothermal compressibility of the SPC/E water using a barostat of both time constants 300 fs and 2 ps. The values for $κT$ were 44.7 and 45.2 Mbar⁻¹ for the lower and higher frequencies, respectively. This is not surprising since the SPC/E model is rigid and likely has no significant coupling to either the low or high frequency barostat. Nevertheless, a more systematic study on the effects of the coupling between the representation of interaction and the barostat and thermostat frequencies is needed.

Having explored the dependence of molecular interaction on mass density fluctuations in bulk, we now move to the vicinity of the air-water interface. It is the mass density fluctuations in the vicinity of the air-water interface that provides the direct connection to the surface tension. Only recently, has the surface tension of water been computed using an energy based methodology.⁶⁹ The advantage of examining the fluctuations directly is that we are not concerned with the convergence of the pressure as was discussed earlier. We begin by computing the mean density profile as a function of the distance from the instantaneous interface. The air-water interface was simulated using a 20 Å × 20 Å × 50 Å slab. Three different water models were used: the DFT-based ab initio model (revPBE-D3), the empirical many body potential model (MB-pol), and the classical potential model (SPC/E). In order to check the quality of DFT slab simulations, we computed the net dipolar moment of the slab confirming that it vanishes within the statistical uncertainty determined by the bulk fluctuations. The instantaneous interface was calculated using the method proposed by Willard and Chandler and has been shown to be a superior coordinate for studying interfaces over the widely used Gibbs dividing surface.¹⁷ Following Willard and Chandler, a coarse-grained time dependent density field was defined as

where $ζ$ is the coarse graining length and r_i is the position of the ith particle at time t. Considering the molecular correlation length of water, the value of $ζ$ is usually chosen to be 2.4 Å. The instantaneous surface is defined by the isosurface h(x, y) having a density equal to half of the bulk density. Once the interface is identified, we then calculate the distance of each water molecule from the instantaneous interface for each configuration (a_i) which is as follows:

where s_i(t) is h(x, y) for the corresponding r(x, y, z) for the ith configuration and n(t) is the surface normal vector at h(x, y). The mean mass density profile (⁠$ρ$⁠) is then calculated as a function of the distance from the instantaneous interface (z) using the following equation:

where L is the length of the simulation cell and $δ$ is the Dirac’s delta function.

Figure 8 depicts the mean density as a function of distance from the instantaneous interface for the water models studied here. In general, the density profiles for all three models (i.e., DFT, SPC/E, and MB-pol) show the well-defined peak with clear minima at the interfacial region indicating that water molecules are more structured in the vicinity of the interface. Our work is consistent with previous studies using both classical and ab initio potentials.^14,17 Among the three models, revPBE-D3 produces a less structured water in the vicinity of the interface than either MB-pol or SPC/E that are both in near quantitative agreement. Both SPC/E and MB-pol have similar density for the interfacial layer of water as reflected in the height of the first peak (i.e., $≈$1.7 gm/cc) in Fig. 8, whereas revPBE provided a lower value of density (i.e., $≈$1.35 gm/cc). The results in Fig. 8 provide an important self-consistent check on the mass density presented in Sec. III A. A good estimate of the mass density can be gleaned from Fig. 8 as it converges to the value of 1 g/cm³ for all interaction potentials in this study. BLYP-D2 simulations that are not shown here have been performed and have shown a similar excellent agreement with the value of 1 g/cm³.^14–16 Our results suggest that all protocols, namely, NpT and NVT in slab geometry, will converge as longer NpT simulations can be performed using barostat time constants on the orders of picoseconds.

Given that there are some quantitative differences in the mass density profile shown in Fig. 8, it would be useful to provide some measure to how the structural averages affect the thermodynamic property of surface tension through an analysis of the fluctuations in height of the instantaneous interface, namely, h(x, y). Figure 9 displays the power spectrum of the instantaneous water-vapor interface. The Fourier transform [$h̃(k)$] of the instantaneous interface, h(x, y) is related to the surface tension through macroscopic capillary-wave theory for wave vectors less than $≈$2$π$/9 Å.⁷⁰ As one can glean from Fig. 9, all models studied herein are nearly indistinguishable in terms of their height fluctuations even though revPBE-D3 produced both a slightly wider and an understructured interface, as was shown in Fig. 8. As a guide, we have plotted the linear response curve that is consistent with the experimental surface tension, $γ$ of 72.0 mJ/m² in the range 0.01 Å⁻¹ $<k<$ 0.7 Å⁻¹. One clearly sees the deviations from linear response as would be expected for short distances (large k). A comparison within the linear regime having smaller k (i.e., 0.01 Å⁻¹ $<k<$ 0.7 Å⁻¹) suggests a qualitative picture where all of the models presented herein provide satisfactory agreement with the experimental surface tension of water.⁷¹ In order to be quantitative, one would have to simulate much larger surfaces in order to have a significant linear region to extract a precise surface tension. It should be noted that the surface tensions of SPC/E and MB-pol are 63.6⁷² and 66.8⁶⁷ mJ/m², respectively. However, we are encouraged by our results presented in this study as they suggest that the fluctuations at scales relevant to modern DFT simulations are accurately represented.

The aforementioned results on surface tension are probing the response of water models to large hydrophobic interfaces at length scales where capillary waves dominate. The opposite limit of small interfaces is equally important and getting the balance between small length scale response and large length scale response correct forms the basis of describing the hydrophobic effect.⁷³ The free energy of forming a large macroscopic interface is given by the surface tension multiplied by the surface area of the interface. This relationship breaks down for very small interfaces at the molecular scale.⁷³ Beyond the hydrophobic effect, these molecular scale interfaces are important for estimating solvation free energies of small molecules and for the hydrophobic interaction between small molecules in solution. For small cavities, the free energy of cavity formation energy can be estimated with the Widom particle insertion formula,

where U_cav is a hard sphere repulsion that acts only on the oxygen atoms out to a radius of R_cav, and $β=1∕kBT$⁠. This expression can be rewritten in terms of the probability of observing a cavity of a given size in pure water,

We can estimate this energy by monitoring the probability of observing a cavity of a given size in a pure water simulation. Figure 10 compares this quantity for the revPBE-D3, MB-pol, and SPC/E models. It is calculated using the slab simulation described above. We see that MB-pol and revPBE-D3 with the MOLOPT basis set agree over the whole size range studied. SPC/E has a comparatively lower cavity formation energy for the larger sizes. A comparison of the cavity formation free energy between different DFT functionals and protocols is discussed in  Appendix B. It suffices to say that differences between basis sets and functionals appear at the larger cavity radii (e.g., $>2.5$ Å) where enhanced sampling methods are needed for proper convergence.⁷⁴ 

Finally, Fig. 11 depicts the changes in the cavity formation energy as a function of z for the revPBE-D3/MOLOPT case in the slab geometry. This is equivalent to the potential of mean force on moving a hard sphere solute across the air-water interface. This is an important quantity for building improved simple models of the distribution of solutes at the air-water interface.^75–77 Again, to probe larger cavity sizes, it will be necessary to use a biasing potential to improve the sampling.⁷⁴ But overall, for smaller cavity radii (e.g., $<$2.5 Å), there is a good agreement between all models indicating that the molecular scale response is robust to all methods studied here.

In conclusion, we have established the simulation protocol for DFT calculations using popular GGA functionals to accurately study the structure and mass density fluctuations of water across different scales. Our results showed that NpT simulations using the standard Fourier interpolation technique at a cutoff of 800 Ry were sufficient to reproduce a converged mass density in agreement with previous studies for the popular revPBE functional along with TZV2P basis set and Grimme D3 dispersion correction. Moreover, we have demonstrated that the quality of the structure and mass density fluctuations is robust across a variety of ensembles for liquid water at 300 K using DFT. Specifically, the simulated structure of water obtained from the NpT ensembles was shown to be consistent with that obtained from the NVT ensembles in both bulk and slab simulation geometries for all DFT functionals. The computed density for all DFT functionals was in the neighborhood of 1 g/cm³ with a slight dependence on protocol. Our research suggests that revPBE-D3 provides an excellent description of water at ambient conditions in agreement with a recent study that used a sophisticated fitting scheme to derive an empirical potential based on revPBE-D3.¹⁰ However, it was recently shown that revPBE-D3 water with the inclusion of NQE was found to worsen the agreement with a variety of experiments.⁵² Interestingly, using the more accurate hybrid density functionals in conjunction with NQE provides an excellent agreement with structural, dynamical, and spectral properties of bulk liquid water.⁵² This finding is consistent with a recent calculation of water clusters up to pentamer revealing that revPBE-D3 benefits from a subtle cancellation of error and that more accurate meta-GGA functionals in conjunction with NQE may provide the correct description of liquid water.⁷⁸ 

In order to ascertain differences between the empirical and the DFT-based interaction potentials for water, we have investigated the local structure in simulations of ambient water by looking into the distribution of the d₅ order parameter that represents the distance of the 5th nearest neighbor from a tagged water molecule. We demonstrated that this order parameter probes the local heterogeneity of water and demonstrates that the DFT-based potentials exhibit a broad range of local environments, in contrast to the empirical models.

In addition to the mass density and local structure, we examined mass density fluctuations and the response to molecular scale and macroscopic (e.g., air-water) interfaces. All empirical models studied produce an isothermal compressibility in agreement with the experimental results. However, the DFT results showed rather large discrepancies depending on the simulation protocol used in this study. Interestingly, when the free energy of forming a molecular sized cavity in water was computed, there was a striking agreement between all representations of interaction. It is interesting to see such good agreement in a free energy when stark differences are present in both local structure and compressibility between the quantum and classical representations of interaction. Examining the response of water to the air-water interface also produced an excellent agreement between all representations of interaction for the system sizes studied herein. Specifically, for the system sizes studied, all three models (revPBE, MB-pol, and SPC/E) were found to qualitatively reproduce the experimental surface tension within the framework of capillary wave theory.

We gratefully acknowledge John Fulton and Sotiris Xantheas for useful discussions and Professor Francesco Paesani for his guidance regarding the studies using MB-pol. We also gratefully acknowledge Professor Joost VandeVondele for discussions regarding the system setup. This work was supported by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences Grant No. BES DE-FG02-09ER46650, which supported MD simulations, data analysis, and manuscript preparation. PNNL is a multiprogram national laboratory operated by Battelle for the U.S. Department of Energy. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under Contract No. DE-AC02-05CH11231. Additional computing resources were generously allocated by PNNL’s Institutional Computing program.

Here we examine the effects of the electron density cutoff in conjunction with the use of grid interpolation techniques that are present in the CP2K code. We performed simulations using the 64 water box with a cutoff of 1000, 2000, and 3000 Ry for 30 ps each. We use the grid interpolation method, i.e., the electron density is calculated using a smoothing protocol⁴³ using the keyword options XC_SMOOTH_RHO NN10 and XC_DERIV SPLINE2_SMOOTH. Figure 12 shows the instantaneous density and the running average with simulation time using the aforementioned cutoffs and grid interpolation technique. All simulations were energy conserving and produce a good liquid structure; however, both the quality of the density fluctuations and the slow convergence of the mass density on the plane wave cutoff are observed in contrast to simulations highlighted in Fig. 1. From the examination of Fig. 12, it is clear that the mass density has not satisfactorily converged even at 3000 Ry.

This can be understood by examining Fig. 13 that displays the variation of the pressure at different electron density cutoffs using both the smoothing and Fourier interpolation techniques. It is clear that convergence of the pressure is achieved at $∼$800 Ry when the Fourier interpolation, namely, using no smoothing protocol, is used. Because grid interpolation requires an abnormally high electron density cutoff, we choose to perform all simulations using the Fourier interpolation. This affords a set of reproducible results as a function of system size and across ensembles where we can confidently focus on the quality of fluctuations that are important to ascertain differences between the descriptions of molecular interaction.

Figure 14 is a comparison between different DFT functionals. We can see that the MOLOPT in the NpT ensemble calculation agrees with the MOLOPT slab calculation. This is to be expected as both of these simulations allow the water cell to fluctuate in size and are at their natural density. On the other hand, the NVT calculation has significantly higher cavity formation energies as the simulation cell cannot fluctuate to compensate for the cavity. The BLYP-D2 and revPBE-D3 using TZV2P basis set results show some differences compared to the revPBE-D3 MOLOPT results for the larger cavity sizes indicating that there is a degree of basis set and functional dependence to this quantity. This study suggests that for cavities with a small radius (⁠$≤2Å$⁠), the NVT ensemble produces converged cavity free energies.

Figure 15 is a comparison for the oxygen-oxygen RDFs between the revPBE and revPBE-D3 for two types of basis sets, TZV2P and MOLOPT. We can see that without the dispersion correction, both the basis sets provide the same RDFs. However, when the dispersion correction is applied, MOLOPT basis set produces a less structured RDF than the TZV2P. This indicates that the origin of the difference between these two basis sets in RDF is due to the matching of the basis set with the D3 dispersion correction. The Grimme’s D3 parameters were originally optimized with the TZV2P basis set, which needs to be re-optimized with MOLOPT to get the consistent RDF from both basis sets.