First principles molecular dynamics simulation protocol is established using revised functional of Perdew-Burke-Ernzerhof (revPBE) in conjunction with Grimme’s third generation of dispersion (D3) correction to describe the properties of water at ambient conditions. This study also demonstrates the consistency of the structure of water across both isobaric (NpT) and isothermal (NVT) ensembles. Going beyond the standard structural benchmarks for liquid water, we compute properties that are connected to both local structure and mass density fluctuations that are related to concepts of solvation and hydrophobicity. We directly compare our revPBE results to the Becke-Lee-Yang-Parr (BLYP) plus Grimme dispersion corrections (D2) and both the empirical fixed charged model (SPC/E) and many body interaction potential model (MB-pol) to further our understanding of how the computed properties herein depend on the form of the interaction potential.

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.1 

The recent empirical corrections due to Grimme32,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.34 

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.10 These potentials have afforded the opportunity to perform simulations for relevant time scales and system sizes to converge the properties of bulk liquid water.10 The results of this study further corroborate some past careful studies using DFT interaction potentials2,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)35 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.11 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.36 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.10 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 (d5). 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.41 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.44 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−1) 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/cm3 at a temperature of 300 K. Both revPBE35 and BLYP45,46 functionals were used with the Grimme dispersion correction32,33 known as D3 and D2, respectively. The core electrons were replaced by the norm-conserving pseudopotentials of Goedecker and co-workers47 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).48 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/cm3.

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.50 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.7 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.6 

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).

1. Mass density and radial distribution functions

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/cm3 and 0.988 g/cm3 with TZV2P and MOLOPT basis sets, respectively. These values are in good agreement with the experimental density of 0.997 g/cm3. The difference of 0.01-0.03 g/cm3 does not account for more than 1% in the lattice constant making up the simulation supercell.

FIG. 1.

The instantaneous density fluctuation and its running average as a function of simulation time from the NpT simulation at revPBE-D3/TZV2P (black) and revPBE-D3/MOLOPT (red) levels of theory.

FIG. 1.

The instantaneous density fluctuation and its running average as a function of simulation time from the NpT simulation at revPBE-D3/TZV2P (black) and revPBE-D3/MOLOPT (red) levels of theory.

Close modal
TABLE I.

Density, compressibility, and structural data obtained from NpT simulations at various levels of theories for bulk water at ambient conditions.

PropertyrevPBE-D3/TZV2PrevPBE-D3/MOLOPTBLYP-D2/TZV2PExpt.
ρ (g/cm30.962 ± 0.029 0.988 ± 0.040 1.04 ± 0.026 0.997 
κT (Mbar−142 … 35 45 
1st max r (Å) 2.80 2.82 2.75 2.80 
1st max gOO(r) 2.74 2.50 3.24 2.57 
1st min r (Å) 3.45 3.66 3.35 3.45 
1st min gOO(r) 0.82 0.91 0.72 0.84 
PropertyrevPBE-D3/TZV2PrevPBE-D3/MOLOPTBLYP-D2/TZV2PExpt.
ρ (g/cm30.962 ± 0.029 0.988 ± 0.040 1.04 ± 0.026 0.997 
κT (Mbar−142 … 35 45 
1st max r (Å) 2.80 2.82 2.75 2.80 
1st max gOO(r) 2.74 2.50 3.24 2.57 
1st min r (Å) 3.45 3.66 3.35 3.45 
1st min gOO(r) 0.82 0.91 0.72 0.84 

Our calculated value of 0.96 g/cm3 is consistent with the previously reported value of 0.96 for revPBE-D2 by Lin et al.4 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/cm3 by Wang et al.12 for revPBE using the nonlocal van der Waals (vdW) correlation functional proposed by Dion et al.13 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.7 The diffusion constant for revPBE-D3 water has recently been calculated by Marsalek and Markland.52 An 800 ps of classical revPBE-D3 simulation provided a system size corrected diffusion coefficient of 2.22±0.05×109 m2 s−1, within statistical error bars of the experimental value of 2.41±0.15×109m2s1.52 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 

FIG. 2.

RDFs for oxygen-oxygen distances: (a) RDFs obtained from the NpT simulations at revPBE-D3/TZV2P (black) and revPBE-D3/MOLOPT (red) basis sets, compared to the experimental RDF (blue dashed) obtained from XRD (taken from Ref. 51) and (b) RDFs obtained from NpT (solid lines) and NVT (dashed lines) ensembles for revPBE-D3/TZV2P (black) and revPBE-D3/MOLOPT (red) levels of theory.

FIG. 2.

RDFs for oxygen-oxygen distances: (a) RDFs obtained from the NpT simulations at revPBE-D3/TZV2P (black) and revPBE-D3/MOLOPT (red) basis sets, compared to the experimental RDF (blue dashed) obtained from XRD (taken from Ref. 51) and (b) RDFs obtained from NpT (solid lines) and NVT (dashed lines) ensembles for revPBE-D3/TZV2P (black) and revPBE-D3/MOLOPT (red) levels of theory.

Close modal

RDFs calculated using the short-range molecular optimized basis set (MOLOPT)48 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/cm3 [see Fig. 3(a)]. This is slightly higher than the reported value of 0.992 g/cm3 (±0.036) by Schmidt et al.6 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/cm3) is also comparable to the BLYP-D3 density reported by Del Ben et al. (1.066 g/cm3)7 and Ma et al. (1.07 g/cm3).8 

FIG. 3.

(a) The instantaneous density fluctuation and its running average as a function of simulation time from the NpT simulation at the BLYP-D2/TZV2P level of theory. The calculated average value and the root mean square deviation are given in Table I; (b) RDFs for oxygen-oxygen distances at the BLYP-D2/TZV2P (black) level of theory, compared to that at revPBE-D3/TZV2P (red) level of theory, and the experimental RDF (blue dashed) obtained from XRD (taken from Ref. 51).

FIG. 3.

(a) The instantaneous density fluctuation and its running average as a function of simulation time from the NpT simulation at the BLYP-D2/TZV2P level of theory. The calculated average value and the root mean square deviation are given in Table I; (b) RDFs for oxygen-oxygen distances at the BLYP-D2/TZV2P (black) level of theory, compared to that at revPBE-D3/TZV2P (red) level of theory, and the experimental RDF (blue dashed) obtained from XRD (taken from Ref. 51).

Close modal

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.

2. Local structure

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.61 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.61 In a recent work,41 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 d5 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 d5 order parameter. We have calculated d5 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 (dji). Then the order parameter d5 is simply the distance between the ith water oxygen and its 5th water oxygen (d5i). Figure 4(a) displays the probability distribution of d5 over all the water molecules in the simulation box for the revPBE-D3, MB-pol and SPC/E waters. The average values of d5 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.

FIG. 4.

Probability distribution of (a) d5 order parameter in bulk ambient water obtained from the NVT simulations at the revPBE-D3/TZV2P (black), MB-pol (red) and SPC/E (blue) levels of theory. Distribution of d5 based on when the 5th water is hydrogen bonded (solid line) to any of the first shell waters vs. that when it is not hydrogen bonded (dashed line) for (b) SPC/E, (c) revPBE-D3, and (d) MB-pol.

FIG. 4.

Probability distribution of (a) d5 order parameter in bulk ambient water obtained from the NVT simulations at the revPBE-D3/TZV2P (black), MB-pol (red) and SPC/E (blue) levels of theory. Distribution of d5 based on when the 5th water is hydrogen bonded (solid line) to any of the first shell waters vs. that when it is not hydrogen bonded (dashed line) for (b) SPC/E, (c) revPBE-D3, and (d) MB-pol.

Close modal

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 d5 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 d5 for each individual water molecule. Figure 5 displays the distribution of the mean value of d5 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 d5. 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 d5 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 d5 under ambient conditions seems to be an indication of the ability of a water model to describe the correct structure under different environments.41 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.

FIG. 5.

Probability distribution of the average value of d5 order parameter for each individual water in the simulation box obtained from the NVT simulations of revPBE-D3/TZV2P (black), MB-pol (red), and SPC/E (blue) levels of theory.

FIG. 5.

Probability distribution of the average value of d5 order parameter for each individual water in the simulation box obtained from the NVT simulations of revPBE-D3/TZV2P (black), MB-pol (red), and SPC/E (blue) levels of theory.

Close modal

1. Isothermal compressibility

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:

κT=V2V2kBTV,
(1)

where V is the instantaneous volume of the system, kB 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−1 for revPBE-D3 and 35 Mbar−1 for BLYP-D2. The revPBE-D3 value, as expected, is very close to the experimental value of 45 Mbar−1 and is in much better agreement than the previously computed value using PBE-D3 (21 Mbar−1) and for PBE0-D3 (32 Mbar−1) by Gaiduk et al.64 Interestingly, both PBE-D3 and PBE0-D3 were reported to produce a density of water very close to the experimental density (1.02 g/cm3 and 0.96 g/cm3, respectively). revPBE-D3 water also produces a compressibility in better agreement with the experiment than vdW-DF, vdW-DFPBE and VV10 levels of theory (18.2, 32.2 and 59.0 Mbar−1, respectively) as computed by Corsetti et al.65 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.66 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.66 The MB-pol model produced an isothermal compressibility of 45.9 Mbar−1 at ambient conditions which is in excellent agreement with the experiment.67 The MB-pol model was also found to produce the experimental pattern of changing isothermal compressibility with temperature.67 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.68 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.

FIG. 6.

The running average of the isothermal compressibility calculated from the NpT simulation at the revPBE-D3/TZV2P (black) and BLYP-D2/TZV2P (red) levels of theory.

FIG. 6.

The running average of the isothermal compressibility calculated from the NpT simulation at the revPBE-D3/TZV2P (black) and BLYP-D2/TZV2P (red) levels of theory.

Close modal

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/cm3 [see Fig. 3(a)] to 0.99 g/cm3. 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.6 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−1 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.

FIG. 7.

(a) The instantaneous density fluctuation and its running average as a function of simulation time and (b) the isothermal compressibility as a function of simulation time. Both were calculated from the NpT simulation at the revPBE-D3/TZV2P (black) and BLYP-D2/TZV2P (red) levels of theory with a thermostat and a barostat frequency of 1 ps.

FIG. 7.

(a) The instantaneous density fluctuation and its running average as a function of simulation time and (b) the isothermal compressibility as a function of simulation time. Both were calculated from the NpT simulation at the revPBE-D3/TZV2P (black) and BLYP-D2/TZV2P (red) levels of theory with a thermostat and a barostat frequency of 1 ps.

Close modal

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.69 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.17 Following Willard and Chandler, a coarse-grained time dependent density field was defined as

ρ(r,t)=i=1N(2πζ2)32exp12|rri(t)|()ζ2,
(2)

where ζ is the coarse graining length and ri 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 (ai) which is as follows:

ai={[si(t)ri(t)]n(t)},
(3)

where si(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:

ρ(z)=1L2i=1Nδ(aiz),
(4)

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/cm3 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/cm3.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.

FIG. 8.

Mean density profile as a function of distance from the instantaneous interface for revPBE-D3 (black), MB-pol (red), and SPC/E (blue) slab geometries.

FIG. 8.

Mean density profile as a function of distance from the instantaneous interface for revPBE-D3 (black), MB-pol (red), and SPC/E (blue) slab geometries.

Close modal

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 Å.70 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/m2 in the range 0.01 Å−1<k< 0.7 Å−1. 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 Å−1<k< 0.7 Å−1) suggests a qualitative picture where all of the models presented herein provide satisfactory agreement with the experimental surface tension of water.71 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.672 and 66.867 mJ/m2, 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.

FIG. 9.

Fourier transform of the instantaneous interface configuration h(x, y) [h̃(k)] for (a) revPBE-D3, (b) MB-pol, and (c) SPC/E. The dashed line is the fit for γ=72.0mJm2 to the capillary wave theory (h̃(k)21βγk2).

FIG. 9.

Fourier transform of the instantaneous interface configuration h(x, y) [h̃(k)] for (a) revPBE-D3, (b) MB-pol, and (c) SPC/E. The dashed line is the fit for γ=72.0mJm2 to the capillary wave theory (h̃(k)21βγk2).

Close modal

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.73 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.73 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,

ΔμXcav=kBTlnexpβUcav0,
(5)

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

ΔμXcav=kBTlnp0(Rcav).
(6)

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.74 

FIG. 10.

Cavity formation free energy for three water models as a function of cavity size calculated using the slab geometry revPBE-D3 (MOLOPT) (black solid line), MB-pol (red dashed line), and SPC/E (blue solid line).

FIG. 10.

Cavity formation free energy for three water models as a function of cavity size calculated using the slab geometry revPBE-D3 (MOLOPT) (black solid line), MB-pol (red dashed line), and SPC/E (blue solid line).

Close modal

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.74 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.

FIG. 11.

Cavity formation free energy with revPBE-D3 (MOLOPT) for several cavity sizes (1 Å, black solid line; 1.5 Å, red solid line; 2.0 Å, green solid line; 2.5 Å, blue solid line; 3.0 Å, black dashed line; and 3.5 Å, red dashed line) as a function of position within the slab. This is equivalent to the potential of mean force for a hard sphere crossing the air-water interface.

FIG. 11.

Cavity formation free energy with revPBE-D3 (MOLOPT) for several cavity sizes (1 Å, black solid line; 1.5 Å, red solid line; 2.0 Å, green solid line; 2.5 Å, blue solid line; 3.0 Å, black dashed line; and 3.5 Å, red dashed line) as a function of position within the slab. This is equivalent to the potential of mean force for a hard sphere crossing the air-water interface.

Close modal

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/cm3 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.10 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.52 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.52 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.78 

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 d5 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 protocol43 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.

FIG. 12.

The instantaneous density fluctuation and its running average as a function of simulation time from the NpT simulation at various density cutoffs (1000 Ry, black; 2000 Ry, red; and 3000 Ry, blue) at the revPBE-D3/TZV2P level of theory with grid interpolation.

FIG. 12.

The instantaneous density fluctuation and its running average as a function of simulation time from the NpT simulation at various density cutoffs (1000 Ry, black; 2000 Ry, red; and 3000 Ry, blue) at the revPBE-D3/TZV2P level of theory with grid interpolation.

Close modal

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.

FIG. 13.

The effect of density cutoff on pressure for Fourier (black circles) and grid (red diamonds) interpolation methods from the NpT simulation at the revPBE-D3/TZV2P level of theory. The inset shows the low cutoff region for Fourier interpolation.

FIG. 13.

The effect of density cutoff on pressure for Fourier (black circles) and grid (red diamonds) interpolation methods from the NpT simulation at the revPBE-D3/TZV2P level of theory. The inset shows the low cutoff region for Fourier interpolation.

Close modal

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.

FIG. 14.

Cavity formation free energy for the different DFT water models as a function of cavity size. The NpT revPBE-D3/MOLOPT (red dashed), revPBE-D3/TZV2P (blue), BLYP-D2/TZV2P (violet dashed), NVT revPBE-D3/MOLOPT (yellow dash-dot), and slab revPBE-D3/MOLOPT (black) calculations agree. There is a non-trivial basis set and functional dependence for the larger cavity sizes.

FIG. 14.

Cavity formation free energy for the different DFT water models as a function of cavity size. The NpT revPBE-D3/MOLOPT (red dashed), revPBE-D3/TZV2P (blue), BLYP-D2/TZV2P (violet dashed), NVT revPBE-D3/MOLOPT (yellow dash-dot), and slab revPBE-D3/MOLOPT (black) calculations agree. There is a non-trivial basis set and functional dependence for the larger cavity sizes.

Close modal

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.

FIG. 15.

RDFs for oxygen-oxygen distances: (a) RDFs obtained from NVT simulations at revPBE/TZV2P (black solid line) and revPBE/MOLOPT (red dashed line) basis sets. (b) RDFs obtained from NVT simulations at revPBE-D3/TZV2P (black solid line) and revPBE-D3/MOLOPT (red dashed line) basis sets.

FIG. 15.

RDFs for oxygen-oxygen distances: (a) RDFs obtained from NVT simulations at revPBE/TZV2P (black solid line) and revPBE/MOLOPT (red dashed line) basis sets. (b) RDFs obtained from NVT simulations at revPBE-D3/TZV2P (black solid line) and revPBE-D3/MOLOPT (red dashed line) basis sets.

Close modal
1.
S.
Yoo
,
X. C.
Zeng
, and
S. S.
Xantheas
, “
On the phase diagram of water with density functional theory potentials: The melting temperature of ice Ih with the Perdew–Burke–Ernzerhof and Becke–Lee–Yang–Parr functionals
,”
J. Chem. Phys.
130
,
221102
(
2009
).
2.
S.
Yoo
and
S. S.
Xantheas
, “
Communication: The effect of dispersion corrections on the melting temperature of liquid water
,”
J. Chem. Phys.
134
,
121105
(
2011
).
3.
M. J.
McGrath
,
J. I.
Siepmann
,
I.-F. W.
Kuo
,
C. J.
Mundy
,
J.
VandeVondele
,
J.
Hutter
,
F.
Mohamed
, and
M.
Krack
, “
Simulating fluid-phase equilibria of water from first principles
,”
J. Phys. Chem. A
110
,
640
646
(
2006
).
4.
I.-C.
Lin
,
A. P.
Seitsonen
,
I.
Tavernelli
, and
U.
Rothlisberger
, “
Structure and dynamics of liquid water from ab initio molecular dynamics: Comparison of BLYP, PBE, and revPBE density functionals with and without van der Waals corrections
,”
J. Chem. Theory Comput.
8
,
3902
(
2012
).
5.
O.
Akin-Ojo
and
F.
Wang
, “
Effects of the dispersion interaction in liquid water
,”
Chem. Phys. Lett.
513
,
59
62
(
2011
).
6.
J.
Schmidt
,
J.
VandeVondele
,
I.-F. W.
Kuo
,
D.
Sebastiani
,
J. I.
Siepmann
,
J.
Hutter
, and
C. J.
Mundy
, “
Isobaric-isothermal molecular dynamics simulations utilizing density functional theory: An assessment of the structure and density of water at near-ambient conditions
,”
J. Phys. Chem. B
113
,
11959
(
2009
).
7.
M.
Del Ben
,
M.
Schönherr
,
J.
Hutter
, and
J.
VandeVondele
, “
Bulk liquid water at ambient temperature and pressure from MP2 theory
,”
J. Phys. Chem. Lett.
4
,
3753
3759
(
2013
).
8.
Z.
Ma
,
Y.
Zhang
, and
M. E.
Tuckerman
, “
Ab initio molecular dynamics study of water at constant pressure using converged basis sets and empirical dispersion corrections
,”
J. Chem. Phys.
137
,
044506
(
2012
).
9.
R.
Jonchiere
,
A. P.
Seitsonen
,
G.
Ferlat
,
A. M.
Saitta
, and
R.
Vuilleumier
, “
Van der Waals effects in ab initio water at ambient and supercritical conditions
,”
J. Chem. Phys.
135
,
154503
(
2011
).
10.
T.
Morawietz
,
A.
Singraber
,
C.
Dellago
, and
J.
Behler
, “
How van der Waals interactions determine the unique properties of water
,”
Proc. Natl. Acad. Sci. U. S. A.
113
,
8368
8373
(
2016
).
11.
R. C.
Remsing
,
M. D.
Baer
,
G. K.
Schenter
,
C. J.
Mundy
, and
J. D.
Weeks
, “
The role of broken symmetry in solvation of a spherical cavity in classical and quantum water models
,”
J. Phys. Chem. Lett.
5
,
2767
2774
(
2014
).
12.
J.
Wang
,
G.
Romań-Ṕerez
,
J. M.
Soler
,
E.
Artacho
, and
M.-V.
Fernández-Serra
, “
Density, structure, and dynamics of water: The effect of van der Waals interactions
,”
J. Chem. Phys.
134
,
024516
(
2011
).
13.
M.
Dion
,
H.
Rydberg
,
E.
Schröder
,
D. C.
Langreth
, and
B. I.
Lundqvist
, “
Van der Waals density functional for general geometries
,”
Phys. Rev. Lett.
92
,
246401
(
2004
).
14.
J.
Kessler
,
H.
Elgabarty
,
T.
Spura
,
K.
Karhan
,
P.
Partovi-Azar
,
A. A.
Hassanali
, and
T. D.
Kühne
, “
Structure and dynamics of the instantaneous water/vapor interface revisited by path-integral and ab initio molecular dynamics simulations
,”
J. Phys. Chem. B
119
,
10079
10086
(
2015
).
15.
T. D.
Kühne
,
T. A.
Pascal
,
E.
Kaxiras
, and
Y.
Jung
, “
New insights into the structure of the vapor/water interface from large-scale first-principles simulations
,”
J. Phys. Chem. Lett.
2
,
105
113
(
2010
).
16.
M. D.
Baer
,
C. J.
Mundy
,
M. J.
McGrath
,
I.-F. W.
Kuo
,
J. I.
Siepmann
, and
D. J.
Tobias
, “
Re-examining the properties of the aqueous vapor–liquid interface using dispersion corrected density functional theory
,”
J. Chem. Phys.
135
,
124712
(
2011
).
17.
A. P.
Willard
and
D.
Chandler
, “
Instantaneous liquid interfaces
,”
J. Phys. Chem. B
114
,
1954
1958
(
2010
).
18.
G.
Miceli
,
S.
de Gironcoli
, and
A.
Pasquarello
, “
Isobaric first-principles molecular dynamics of liquid water with nonlocal van der Waals interactions
,”
J. Chem. Phys.
142
,
034501
(
2015
).
19.
I.-C.
Lin
,
A. P.
Seitsonen
,
M. D.
Coutinho-Neto
,
I.
Tavernelli
, and
U.
Rothlisberger
, “
Importance of van der Waals interactions in liquid water
,”
J. Phys. Chem. B
113
,
1127
1131
(
2009
).
20.
K.
Forster-Tonigold
and
A.
Groß
, “
Dispersion corrected rpbe studies of liquid water
,”
J. Chem. Phys.
141
,
064501
(
2014
).
21.
M. J.
Gillan
,
D.
Alfè
, and
A.
Michaelides
, “
Perspective: How good is DFT for water?
,”
J. Chem. Phys.
144
,
130901
(
2016
).
22.
R. A.
DiStasio
, Jr.
,
B.
Santra
,
Z.
Li
,
X.
Wu
, and
R.
Car
, “
The individual and collective effects of exact exchange and dispersion interactions on the ab initio structure of liquid water
,”
J. Chem. Phys.
141
,
084502
(
2014
).
23.
R. Z.
Khaliullin
and
T. D.
Kühne
, “
Microscopic properties of liquid water from combined ab initio molecular dynamics and energy decomposition studies
,”
Phys. Chem. Chem. Phys.
15
,
15746
15766
(
2013
).
24.
T. D.
Kühne
,
M.
Krack
, and
M.
Parrinello
, “
Static and dynamical properties of liquid water from first principles by a novel Car–Parrinello-like approach
,”
J. Chem. Theory Comput.
5
,
235
241
(
2009
).
25.
H.-S.
Lee
and
M. E.
Tuckerman
, “
Structure of liquid water at ambient temperature from ab initio molecular dynamics performed in the complete basis set limit
,”
J. Chem. Phys.
125
,
154507
(
2006
).
26.
G. R.
Medders
,
A. W.
Gotz
,
M. A.
Morales
,
P.
Bajaj
, and
F.
Paesani
, “
On the representation of many-body interactions in water
,”
J. Chem. Phys.
143
,
104102
(
2015
).
27.
G. A.
Cisneros
,
K. T.
Wikfeldt
,
L.
Ojamae
,
J.
Lu
,
Y.
Xu
,
H.
Torabifard
,
A. P.
Bartok
,
G.
Csanyi
,
V.
Molinero
, and
F.
Paesani
, “
Modeling molecular interactions in water: From pairwise to many-body potential energy functions
,”
Chem. Rev.
116
,
7501
7528
(
2016
).
28.
V.
Babin
,
G. R.
Medders
, and
F.
Paesani
, “
Toward a universal water model: First principles simulations from the dimer to the liquid phase
,”
J. Phys. Chem. Lett.
3
,
3765
3769
(
2012
).
29.
G. R.
Medders
,
V.
Babin
, and
F.
Paesani
, “
A critical assessment of two-body and three-body interactions in water
,”
J. Chem. Theory Comput.
9
,
1103
1114
(
2013
).
30.
G. S.
Fanourgakis
and
S. S.
Xantheas
, “
Development of transferable interaction potentials for water. V. Extension of the flexible, polarizable, Thole-type model potential (TTM3–F, v. 3.0) to describe the vibrational spectra of water clusters and liquid water
,”
J. Chem. Phys.
128
,
074506
(
2008
).
31.
C.
Burnham
,
D.
Anick
,
P.
Mankoo
, and
G.
Reiter
, “
The vibrational proton potential in bulk liquid water and ice
,”
J. Chem. Phys.
128
,
154519
(
2008
).
32.
S.
Grimme
, “
Semiempirical GGA-type density functional constructed with a long-range dispersion correction
,”
J. Comput. Chem.
27
,
1787
(
2006
).
33.
S.
Grimme
,
J.
Antony
,
S.
Ehrlich
, and
H.
Krieg
, “
A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H–Pu
,”
J. Chem. Phys.
132
,
154104
(
2010
).
34.
M. D.
Baer
,
W. I.-F.
Kuo
,
D. J.
Tobias
, and
C. J.
Mundy
, “
Toward a unified picture of the water self-ions at the air–water interface: A density functional theory perspective
,”
J. Phys. Chem. B
118
,
8364
8372
(
2014
).
35.
Y.
Zhang
and
W.
Yang
, “
Comment on ‘generalized gradient approximation made simple
,’”
Phys. Rev. Lett.
80
,
890
(
1998
).
36.
M. D.
Baer
,
V.-T.
Pham
,
J. L.
Fulton
,
G. K.
Schenter
,
M.
Balasubramanian
, and
C. J.
Mundy
, “
Is iodate a strongly hydrated cation?
,”
J. Phys. Chem. Lett.
2
,
2650
2654
(
2011
).
37.
D. J.
Arismendi-Arrieta
,
M.
Riera
,
P.
Bajaj
,
R.
Prosmiti
, and
F.
Paesani
, “
i-TTM model for ab initio-based ion−water interaction potentials. 1. Halide−water potential energy functions
,”
J. Phys. Chem. B
120
,
1822
1832
(
2016
).
38.
M.
Riera
,
A. W.
Gotz
, and
F.
Paesani
, “
The i-TTM model for ab initio-based ion–water interaction potentials. II. Alkali metal ion–water potential energy functions
,”
Phys. Chem. Chem. Phys.
18
,
30334
30343
(
2016
).
39.
D. G.
Triezenberg
and
R.
Zwanzig
, “
Fluctuation theory of surface tension
,”
Phys. Rev. Lett.
28
,
1183
(
1972
).
40.
A.
Nilsson
,
C.
Huang
, and
L. G.
Pettersson
, “
Fluctuations in ambient water
,”
J. Mol. Liq.
176
,
2
16
(
2012
).
41.
L. B.
Skinner
,
M.
Galib
,
J. L.
Fulton
,
C. J.
Mundy
,
J. B.
Parise
,
V.-T.
Pham
,
G. K.
Schenter
, and
C. J.
Benmore
, “
The structure of liquid water up to 360 MPa from x-ray diffraction measurements using a high Q-range and from molecular simulation
,”
J. Chem. Phys.
144
,
134504
(
2016
).
42.
G.
Lippert
,
J.
Hutter
, and
M.
Parrinello
, “
A hybrid Gaussian and plane wave density functional scheme
,”
Mol. Phys.
92
,
477
(
1997
).
43.
J.
VandeVondele
,
M.
Krack
,
F.
Mohamed
,
M.
Parrinello
,
T.
Chassaing
, and
J.
Hutter
, “
Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach
,”
Comput. Phys. Commun.
167
,
103
(
2005
).
44.
M. E.
Tuckerman
,
J.
Alejandre
,
R.
López-Rendón
,
A. L.
Jochim
, and
G. J.
Martyna
, “
A Liouville-operator derived measure-preserving integrator for molecular dynamics simulations in the isothermal–isobaric ensemble
,”
J. Phys. A: Math. Gen.
39
,
5629
(
2006
).
45.
A. D.
Becke
, “
Density-functional exchange-energy approximation with correct asymptotic behavior
,”
Phys. Rev. A
38
,
3098
(
1988
).
46.
C.
Lee
,
W.
Yang
, and
R. G.
Parr
, “
Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density
,”
Phys. Rev. B
37
,
785
(
1988
).
47.
S.
Goedecker
,
M.
Teter
, and
J.
Hutter
, “
Separable dual-space Gaussian pseudopotentials
,”
Phys. Rev. B
54
,
1703
(
1996
).
48.
J.
VandeVondele
and
J.
Hutter
, “
Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases
,”
J. Chem. Phys.
127
,
114105
(
2007
).
49.
A.
Bankura
,
A.
Karmakar
,
V.
Carnevale
,
A.
Chandra
, and
M. L.
Klein
, “
Structure, dynamics, and spectral diffusion of water from first-principles molecular dynamics
,”
J. Phys. Chem. C
118
,
29401
29411
(
2014
).
50.
M. J.
McGrath
,
J. I.
Siepmann
,
I.-F. W.
Kuo
,
C. J.
Mundy
,
J.
VandeVondele
,
J.
Hutter
,
F.
Mohamed
, and
M.
Krack
, “
Isobaric-isothermal Monte Carlo simulations from first principles: Application to liquid water at ambient conditions
,”
ChemPhysChem
6
,
1894
(
2005
).
51.
L. B.
Skinner
,
C.
Huang
,
D.
Schlesinger
,
L. G. M.
Pettersson
,
A.
Nilsson
, and
C. J.
Benmore
, “
Benchmark oxygen-oxygen pair-distribution function of ambient water from x-ray diffraction measurements with a wide Q-range
,”
J. Chem. Phys.
138
,
074506
(
2013
).
52.
O.
Marsalek
and
T. E.
Markland
, “
Quantum dynamics and spectroscopy of ab initio liquid water: The interplay of nuclear and electronic quantum effects
,”
J. Phys. Chem. Lett.
8
,
1545
1551
(
2017
).
53.
R. A.
Kuharski
and
P. J.
Rossky
, “
A quantum mechanical study of structure in liquid H2O and D2O
,”
J. Chem. Phys.
82
,
5164
(
1985
).
54.
J.
Lobaugh
and
G. A.
Voth
, “
A quantum model for water: Equilibrium and dynamical properties
,”
J. Chem. Phys.
106
,
2400
(
1997
).
55.
M. W.
Mahoney
and
W. L.
Jorgensen
, “
Quantum, intramolecular flexibility, and polarizability effects on the reproduction of the density anomaly of liquid water by simple potential functions
,”
J. Chem. Phys.
115
,
10758
(
2001
).
56.
J. A.
Morrone
and
R.
Car
, “
Nuclear quantum effects in water
,”
Phys. Rev. Lett.
101
,
017801
(
2008
).
57.
M.
Ceriotti
,
W.
Fang
,
P. G.
Kusalik
,
R. H.
McKenzie
,
A.
Michaelides
,
M. A.
Morales
, and
T. E.
Markland
, “
Nuclear quantum effects in water and aqueous systems: Experiment, theory, and current challenges
,”
Chem. Rev.
116
,
7529
7550
(
2016
).
58.
A.
Møgelhøj
,
A. K.
Kelkkanen
,
K. T.
Wikfeldt
,
J.
Schiøtz
,
J. J.
Mortensen
,
L. G.
Pettersson
,
B. I.
Lundqvist
,
K. W.
Jacobsen
,
A.
Nilsson
, and
J. K.
Nørskov
, “
Ab initio van der Waals interactions in simulations of water alter structure from mainly tetrahedral to high−density−like
,”
J. Phys. Chem. B
115
,
14149
14160
(
2011
).
59.
J.
VandeVondele
,
F.
Mohamed
,
M.
Krack
,
J.
Hutter
,
M.
Sprik
, and
M.
Parrinello
, “
The influence of temperature and density functional models in ab initio molecular dynamics simulation of liquid water
,”
J. Chem. Phys.
122
,
014515
(
2005
).
60.
M. J.
McGrath
,
I.-F. W.
Kuo
, and
J. I.
Siepmann
, “
Liquid structures of water, methanol, and hydrogen fluoride at ambient conditions from first principles molecular dynamics simulations with a dispersion corrected density functional
,”
Phys. Chem. Chem. Phys.
13
,
19943
19950
(
2011
).
61.
E.
Duboué-Dijon
and
D.
Laage
, “
Characterization of the local structure in liquid water by various order parameters
,”
J. Phys. Chem. B
119
,
8406
8418
(
2015
).
62.
M. J.
Cuthbertson
and
P. H.
Poole
, “
Mixturelike behavior near a liquid-liquid phase transition in simulations of supercooled water
,”
Phys. Rev. Lett.
106
,
115706
(
2011
).
63.
R. S.
Singh
,
J. W.
Biddle
,
P. G.
Debenedetti
, and
M. A.
Anisimov
, “
Two-state thermodynamics and the possibility of a liquid-liquid phase transition in supercooled TIP4P/2005 water
,”
J. Chem. Phys.
144
,
144504
(
2016
).
64.
A. P.
Gaiduk
,
F.
Gygi
, and
G.
Galli
, “
Density and compressibility of liquid water and ice from first-principles simulations with hybrid functionals
,”
J. Phys. Chem. Lett.
6
,
2902
2908
(
2015
).
65.
F.
Corsetti
,
E.
Artacho
,
J. M.
Soler
,
S. S.
Alexandre
, and
M.-V.
Fernández-Serra
, “
Room temperature compressibility and diffusivity of liquid water from first principles
,”
J. Chem. Phys.
139
,
194502
(
2013
).
66.
L. P.
Helena
,
J. L.
Aragones
,
C.
Vega
,
E. G.
Noya
,
J. L. F.
Abascal
,
M. A.
Gonzalez
, and
C.
McBride
, “
Anomalies in water as obtained from computer simulations of the TIP4P/2005 model: Density maxima, and density, isothermal compressibility and heat capacity minima
,”
Mol. Phys.
107
,
365
374
(
2009
).
67.
S. K.
Reddy
,
S. C.
Straight
,
P.
Bajaj
,
C. H.
Pham
,
M.
Riera
,
D. R.
Moberg
,
M. A.
Morales
,
C.
Knight
,
A. W.
Gotz
, and
F.
Paesani
, “
On the accuracy of the MB-pol many-body potential for water: Interaction energies, vibrational frequencies, and classical thermodynamic and dynamical properties from clusters to liquid water and ice
,”
J. Chem. Phys.
145
,
194504
(
2016
).
68.
H.
Pathak
,
J.
Palmer
,
D.
Schlesinger
,
K. T.
Wikfeldt
,
J. A.
Sellberg
,
L. G.
Pettersson
, and
A.
Nilsson
, “
The structural validity of various thermodynamical models of supercooled water
,”
J. Chem. Phys.
145
,
134507
(
2016
).
69.
Y.
Nagata
,
T.
Ohto
,
M.
Bonn
, and
T. D.
Kühne
, “
Surface tension of ab initio liquid water at the water-air interface
,”
J. Chem. Phys.
144
,
204705
(
2016
).
70.
F.
Sedlmeier
,
D.
Horinek
, and
R. R.
Netz
, “
Nanoroughness, intrinsic density profile, and rigidity of the air-water interface
,”
Phys. Rev. Lett.
103
,
136102
(
2009
).
71.
A. J.
Patel
,
P.
Varilly
,
S. N.
Jamadagni
,
H.
Acharya
, and
S. G. D.
Chandler
, “
Extended surfaces modulate hydrophobic interactions of neighboring solutes
,”
Proc. Natl. Acad. Sci. U. S. A.
108
,
17678
17683
(
2011
).
72.
C.
Vega
and
E.
de Miguel
, “
Surface tension of the most popular models of water by using the test-area simulation method
,”
J. Chem. Phys.
126
,
154707
(
2007
).
73.
D.
Chandler
, “
Interfaces and the driving force of hydrophobic assembly
,”
Nature
437
,
640
647
(
2005
).
74.
E.
Xi
,
R. C.
Remsing
, and
A. J.
Patel
, “
Sparse sampling of water density fluctuations in interfacial environments
,”
J. Chem. Theory Comput.
12
,
706
713
(
2016
).
75.
Y.
Levin
,
A. P.
dos Santos
, and
A.
Diehl
, “
Ions at the air−water interface: An end to a hundred−year−old mystery?
,”
Phys. Rev. Lett.
103
,
257802
(
2009
).
76.
T. T.
Duignan
,
D. F.
Parsons
, and
B. W.
Ninham
, “
Ion interactions with the air-water interface using a continuum solvent model
,”
J. Phys. Chem. B
118
,
8700
8710
(
2014
).
77.
S.
Vaikuntanathan
,
G.
Rotskoff
,
A.
Hudson
, and
P. L.
Geissler
, “
Necessity of capillary modes in a minimal model of nanoscale hydrophobic solvation
,”
Proc. Natl. Acad. Sci. U. S. A.
113
,
E2224
E2230
(
2016
).
78.
L. R.
Pestana
,
N.
Mardirossian
,
M.
Head-Gordon
, and
T.
Head-Gordon
, “
Ab initio molecular dynamics simulations of liquid water using high quality meta−GGA functionals
,”
Chem. Sci.
8
,
3554
3565
(
2017
).