Nuclear quantum effects in liquid water have profound implications for several of its macroscopic properties related to the structure, dynamics, spectroscopy, and transport. Although several of water’s macroscopic properties can be reproduced by classical descriptions of the nuclei using interaction potentials effectively parameterized for a narrow range of its phase diagram, a proper account of the nuclear quantum effects is required to ensure that the underlying molecular interactions are transferable across a wide temperature range covering different regions of that diagram. When performing an analysis of the hydrogen-bonded structural networks in liquid water resulting from the classical (class) and quantum (qm) descriptions of the nuclei with two interaction potentials that are at the two opposite ends of the range in describing quantum effects, namely the flexible, pair-wise additive q-TIP4P/F, and the flexible, polarizable TTM3-F, we found that the (class) and (qm) results can be superimposed over the temperature range T = 250-350 K using a surprisingly simple, linear scaling of the two temperatures according to T(qm) = α T(class) + ΔT, where α = 0.99 and ΔT = −6 K for q-TIP4P/F and α = 1.24 and ΔT = −64 K for TTM3-F. This simple relationship suggests that the structural networks resulting from the quantum and classical treatment of the nuclei with those two very different interaction potentials are essentially similar to each other over this extended temperature range once a model-dependent linear temperature scaling law is applied.

Owing to its importance in life emanating from its function as a solvent and the principal medium for biological activity, water has received a lot of attention aimed at the understanding of the structure and dynamics of its fleeting hydrogen-bonded structural networks to the exploration of its phase diagram and its hydration functions for several chemical and biological species. Central to the desired unified description of water’s properties across its phase diagram is the ability to describe the underlying intermolecular interactions and collective phenomena on equal footing in a wide range of macroscopic conditions (such as temperature and pressure) that are associated with dissimilar structural motifs. Of particular importance is the identification of the structural patterns and corresponding dynamics of its transitional hydrogen-bonded structural network that is responsible for many of water’s anomalous properties. Since the first simple model for water introduced in 1933 by Bernal and Fowler1 and the first computer simulations of the liquid by Barker and Watts2 in 1969 and by Rahman and Stillinger3 in 1971, there has been a plethora of models used to describe the underlying interactions with increasing degree of sophistication; see a recent review4 for the range of interaction potentials developed until 2002.

Quantum effects have been shown to affect the macroscopic properties of liquid water. It can be argued in very general terms that any difference in a thermodynamic property of H2O vs. D2O is due to nuclear quantum effects, since in classical mechanics the nuclear masses kinetic energy contribution to the partition function separates from that of the potential energy, and hence can be integrated out.5 Most obvious in this regard is the experimentally observed higher melting temperature of D2O (3.8 °C), and even more so of T2O (4.5 °C), both of which are closer to classical water in their thermodynamic properties when compared to H2O. In very simple words, the shift in the melting temperature can be explained by zero-point energy, which results in a kinetic energy larger than kBT/2 per degree of freedom and consequently in stronger hydrogen bonds than H2O, so one needs a larger temperature in the classical case to match the properties obtained from quantum simulations (it is however established by now that this picture is oversimplified).6 

The majority of existing interaction potentials for water have been parameterized to describe a set of macroscopic properties in a narrow range of the phase diagram via classical (Newtonian) simulations, and as such they can be viewed as models that effectively include quantum effects in the potential. Most of them fail to describe parts of the phase diagram that are different from the ones they were parameterized.7 In contrast, when the underlying interactions are obtained from first principles (i.e., in an ab initio manner without any adjustable parameters), it is expected that these would be transferable across different environments, provided that the level of electronic structure theory used is capable of describing the subtle changes in the various components of the interaction energy that occur upon the change of the underlying hydrogen-bonding network.8 Moreover, in order for the ab initio based potentials to be transferable across different environments, explicit account of the zero-point energy should be considered, usually included via nuclear statistical mechanical simulations, viz. Path Integral Molecular Dynamics (PIMD) or Centroid Molecular Dynamics (CMD) approaches.9–19 Nuclear quantum statistical simulations with electronic structure (DFT-based) potentials20 have already been reported for liquid water and ice;21,22 however, due to their computational cost, they cannot currently be adopted to explore larger parts of the phase diagram of water or employ more accurate but computationally expensive functionals. To this end, alternative approaches, aiming at combining the merits of the two practices, have been recently introduced.23 Among the transferable classical potentials for water, the development of the family of Thole-Type Model (TTM2.1-F and TTM3-F) interaction potentials24–26 has relied on a different philosophy than the rest of the classical potentials: highly accurate electronic structure binding energies of small clusters27–29 were used to fit the interactions between fragments, which are themselves described by an ab initio derived analytical monomer potential and dipole moment surface (DMS),30 appropriately modified to describe the condensed phase;31 as such they are appropriate for nuclear quantum PIMD and CMD, rather than classical (Newtonian) statistical simulations. In conjunction with quantum simulation protocols, they have been previously used to obtain structural, dynamical, spectral, thermodynamic, and transport properties of liquid water20,31–35 with remarkable accuracy. This approach is distinctively different than others prior to that effort, which, for instance, have relied on reproducing salient features such as the rotational barriers that have been experimentally measured for small water clusters (dimer and trimer).36,37 An extension of this approach fitting to several tens of thousands of water dimer and trimer configurations has recently resulted in the development of the MB-pol potential.38–40 

On the opposite side of the range of water models associated with the description of quantum effects, several pair-wise additive potentials have been re-parametrized (starting from TIP4P/200541) to reproduce macroscopic properties of liquid water when including nuclear quantum effects explicitly in CMD or PIMD simulations. For example, the rigid, pair-wise additive TIP4PQ/2005 interaction potential treating long range Coulombic interactions with either the reaction field or Ewald summation methods was found to reproduce several macroscopic properties of liquid water.42 The q-TIP4P/F potential is a flexible pairwise-additive model that also has been parametrized to effectively describe quantum effects in liquid water.6 CMD simulations have been reported with the TIP4P and SPC/E interaction potentials as well;43 however, since these potentials already include quantum effects in an effective manner in connection with classical simulations, the quantum simulations double-count quantum effects.

The early work of Rossky and co-workers9 suggested that quantum effects induce changes in the liquid’s structural properties of the order of an increase of ∼50 °C, whereas Stern and Berne10 reported similar results, viz. a less structured liquid and a correction to the binding enthalpy by −1.5 kcal/mol when compared to the classical results. Fanourgakis et al.32 were the first ones to estimate the magnitude of quantum effects in liquid water at T = 298 K with the TTM2.1-F potential, whereas Paesani et al.33 subsequently reported quantitative comparisons between classical and quantum simulations for the structure, enthalpy, and diffusion coefficient with the same force field over an extended (∼100 K) temperature range. They suggested a value of 25–30 K for this effect based on simulations with the flexible, polarizable TTM2.1-F potential. They also demonstrated that this is not a constant temperature correction, since quantum effects vary in a non-linear fashion as a function of temperature. The same authors have also reported the impact of quantum effects in the pre-melting of the surface of ice,44 whereas subsequent versions of the potential (TTM3-F)31 have been used in conjunction with quantum45 and semi-classical35 simulations to model the infrared (IR) spectra in liquid water. In addition, quantum fluctuations have been reported to either promote or inhibit glass formation46 and as such can play a role in the area of water’s phase diagram defined by the glass transition and homogeneous nucleation temperatures (also known as “no-man’s land”).33,47 Although previous studies have indicated that quantum effects have a ∼30–50 K effect on the melting temperature with some models such as ST2 or q-SPC/Fw, Manolopoulos and co-workers6 have found a difference of just 8 K between the classical and quantum melting points of the q-TIP4P/F potential. They also reported that quantum effects change the diffusion coefficient of the q-TIP4P/F model by a factor of 1.1 and the diffusion coefficient of the TTM3-F model by a factor of 0.95. In this respect, the TTM3-F potential behaves quite differently than the TTM2.1-F (i.e., quantum effects for the diffusion coefficient vary from a quantum/classical ratio of 0.95 for TTM3-F to 1.5× for TTM2.1-F). In both cases (melting point and diffusion coefficient), the small(er) net quantum effects with the q-TIP4P/F and TTM3-F potentials have been attributed by Manolopoulos and co-workers6 to the competition between opposing inter- and intra-molecular effects. In simple words, the inter-molecular effect originates from quantum delocalization, tentatively lowering the structuring of water, while the intra-molecular effect originates from the anharmonicity of the OH-stretch bond that increases its average dipole in the quantum case and thus the strength of hydrogen bonds. Markland and Berne, in their recent study of quantum mechanical effects in water using isotopic fractionation,48 have cautioned that these competing quantum effects may be overpredicted by the TTM3-F potential. A recent publication49 has investigated the shift of the phase diagram of the various forms of crystalline ice with the TIP4PQ/2005 potential using classical and quantum Path Integral Monte Carlo (PIMC) simulations. The overall structure of the phase diagram computed with the two protocols was very similar, with a shift of 15-20 K to lower temperatures.

The concept of using a temperature scale relative to the melting temperature (Tm), consistent with the temperature shifts used to superimpose structural features between classical and quantum simulations first introduced by Kuharski and Rossky,9 has been previously invoked34 when comparing the experimentally measured anisotropy of liquid water to the one computed via classical and quantum simulations with the TTM3-F potential. However, the melting temperature is just one point in the phase diagram and it might not necessarily be the most decisive one in elucidating the structural properties of liquid water; for instance, there is no apparent correlation between Tm and the density maximum for various water models since temperature differences vary from 11 to 37 K.50 Building upon the previously identified temperature shifts and/or scaling between the classical and quantum results, we seek to further quantify these effects and validate earlier insights over a larger temperature range. To that end, here we rely on tools of complex network analysis to compare the structural properties around a given water molecule obtained from classical and quantum simulations of liquid water over a wide temperature range (270-350 K). We consider the TTM3-F and q-TIP4P/F interaction potentials for water, which can be thought of lying at the opposite ends of the range of potentials describing quantum effects, as they exhibit competing quantum effects with different strengths.6 That is, in the case of TTM3-F, the inter-molecular quantum effect overcompensates the intra-molecular contribution, in contrast to q-TIP4P/F, for which the competition in the intra-molecular effect still dominates in the relevant temperature range.48 We characterize the local structure of water through hydrogen-bond patterns, schematically depicted in Fig. 1, which are identified along the simulation trajectories by a clustering algorithm. We have recently introduced that approach to explore the free-energy surface of simpler water models in connection with purely classical simulations.51,52

FIG. 1.

Hydrogen-bonded structures that act as attractors of the seven most populated free energy basins (gradient clusters).

FIG. 1.

Hydrogen-bonded structures that act as attractors of the seven most populated free energy basins (gradient clusters).

Close modal

Classical and Path Integral Molecular Dynamics (PIMD) simulations at constant temperature and pressure conditions (NPT ensemble) with the q-TIP4P/F6 and TTM3-F31 potentials were initially performed to determine the liquid density for the temperature range T = 250-350 K and pressure P = 1 atm. For the path integral simulations, a total of n = 32 ring polymer beads were used to ensure converged results. The Andersen thermostat53,54 and Berendsen barostat54 were employed to maintain the temperature and pressure, respectively. Cubic cells consisting of 128 water molecules with standard periodic boundary conditions in the three directions were employed for the description of the liquid phase. Real space electrostatic and van der Waals interactions were considered up to a cut-off distance (Rc) equal to half the box size, while for interactions beyond Rc, the standard long-range corrections for the van der Waals55 and the Ewald summation technique56 for the electrostatic interactions were employed. A multiple time-step scheme combined with the velocity Verlet algorithm was employed for the integration of the classical and quantum trajectories. The intramolecular forces were computed every 0.25 fs, while the intermolecular forces were updated every 0.5 fs. The quantum simulations were significantly sped up by using the recently developed57 fast path integral method for polarizable force fields. Simulations at constant volume (NVT ensemble) were performed at the density determined before from the (NPT) simulations. The atomic coordinates obtained during classical trajectories or from the individual beads of the ring-polymer during the path integral simulation were saved to the disk for further analysis every 4 fs. The total integration time for the simulation at each temperature was 1 ns or more. All simulations were performed using an in-house code that has been efficiently parallelized.58 

Our approach is based on building a conformation space network (as described in the next paragraph), which requires a definition for the formation of a hydrogen bond. We used the hydrogen-bond criterion introduced by Skinner and co-workers,59 which is based on a molecular orbital point of view since it defines the strength of a hydrogen bond via the amount of electron transfer (called the “occupancy N” in the following) from the lone pair of the oxygen of the hydrogen-bond accepting water molecules into the antibonding σOH orbital of the hydrogen-bond donating water molecule. It has been shown in Ref. 59 that the occupancy N correlates extremely well with

Nr,ψ=er0.343Å(7.10.05ψ+0.00021ψ2),
(1)

where r is the intermolecular O⋯H distance and ψ the angle between the intermolecular O⋯H vector and the normal to the intramolecular plane of the hydrogen-bond accepting water molecule. Unless otherwise noted, we considered the coordinates of individual beads (instead of the centroid) from the PIMD to determine the occupancy N from Eq. (1). The two water models considered in this study (q-TIP4P/F and TTM3-F) are flexible, in contrast to the ones used in Ref. 59. Before determining the hydrogen bonds, each water molecule from the MD trajectory is therefore replaced by a rigid counterpart with rOH = 0.983 Å and bond angle 106°, aligning the bisector between both OH groups, the molecule’s planes and the center of mass.

The distribution of occupancies is bimodal (see Fig. S1 of the supplementary material), with a high occupancy indicating a hydrogen bond and a low occupancy indicating a broken hydrogen bond, respectively. Hence, the minimum of that distribution may serve as a cut-off value to distinguish an intact hydrogen bond from a broken one. The precise position of this minimum depends on both the water model and the temperature. For consistency, we chose the same value Ncut = 0.0085 as in Ref. 59, which was used to derive the cut-off value for the SPC/E water model at 300 K, despite the fact that the minimum is found at somewhat lower occupancy values for the TTM3-F and q-TIP4P/F (N ≈ 0.0071 at 300 K for both models). In either case, the minimum is not very pronounced, indicating that the distributions of hydrogen-bonded versus broken water molecules strongly overlap with respect to the occupancy N. As such, any choice of a cut-off is to a certain extent arbitrary. Nevertheless, for a direct comparison of the amount of hydrogen bonding in the various models and at different temperatures, it is appropriate to use the same cut-off to analyze the results of all simulations.

The main idea behind a conformation space network is to map the dynamics obtained from a MD trajectory onto a conformation space network, where nodes and links represent microstates (i.e., local structures) and the transitions between them, respectively. The conformation space network is constructed in essentially the same way as in Refs. 51 and 52. In brief, we first define microstates that represent the local structure around a given water molecule by evaluating the connectivity to neighboring water molecules via hydrogen-bonds, using the hydrogen bond definition described in the previous paragraph.59 Water molecules up to the second hydration layer are considered. Figure 1 shows representative examples of such hydrogen-bonded structural networks, but structures can be much more complicated, containing, for example, loops (see Fig. 1 of Ref. 51). We assume that each water molecule has a maximum of four binding sites. In relatively rare occasions (ca. 2%) when two water molecules hydrogen-bond to one of the hydrogens, or three water molecules hydrogen-bond to the oxygen, we selected those with the stronger hydrogen bond. In contrast to Refs. 51 and 52, we no longer distinguish the two hydrogen sites of the central water molecule; all symmetries due to the exchange of identical particles (i.e., water molecules and the two hydrogens of each individual water molecule) are resolved.

The time evolution of these microstates reveals the conformation space network [see Figs. 2(a) and 2(b)]. This is quantified by counting how often each microstate (node) is visited along a sufficiently long MD trajectory, indicated in Figs. 2(a) and 2(b) by the area of the circles. Furthermore, when two microstates appear at subsequent time steps in the MD trajectory for a given water molecule, they are connected by a link [represented by the lines in Figs. 2(a) and 2(b)]. The thickness of these lines is proportional to the number of times such a transition occurs. The conformation space network is often referred to as a representation of the transition probability matrix in a Markov-State-Model (assuming that the Markov property is fulfilled).60 The analysis of the MD trajectories at different temperatures is based on the same (finite) set of possible microstates; however, the statistical weights of both the microstates (nodes) and the links between them will be different [compare Fig. 2(a)vs. Fig. 2(b)].

FIG. 2.

Schematic representation of the “full” conformation space network [panels (a) and (b)] and the reduced conformation network [“gradient clusters,” panels (c) and (d)] that remains after deleting links in the gradient clustering algorithm (see text and Ref. 51 for details). The area of the nodes represents the statistical weights of the microstates, and the thickness of the lines represents those of the links. Two examples are shown: (a) and (c) at a low temperature, where the fully hydrogen-bonded structure #1 (see Fig. 1) acts as the only attractor on a funnel-like free energy surface, and (b) and (d) at a higher temperature, when an entropically stabilized side-minimum evolves, since the statistical weights of both nodes and links change as a function of temperature. In the second case, the gradient-clustering fragments the full conformation space network into two basins. Panels (e) and (f) show the corresponding free energy landscapes along an abstract (i.e., not necessarily physical) reaction coordinate, i.e., the x-axis in panels (a) and (d).

FIG. 2.

Schematic representation of the “full” conformation space network [panels (a) and (b)] and the reduced conformation network [“gradient clusters,” panels (c) and (d)] that remains after deleting links in the gradient clustering algorithm (see text and Ref. 51 for details). The area of the nodes represents the statistical weights of the microstates, and the thickness of the lines represents those of the links. Two examples are shown: (a) and (c) at a low temperature, where the fully hydrogen-bonded structure #1 (see Fig. 1) acts as the only attractor on a funnel-like free energy surface, and (b) and (d) at a higher temperature, when an entropically stabilized side-minimum evolves, since the statistical weights of both nodes and links change as a function of temperature. In the second case, the gradient-clustering fragments the full conformation space network into two basins. Panels (e) and (f) show the corresponding free energy landscapes along an abstract (i.e., not necessarily physical) reaction coordinate, i.e., the x-axis in panels (a) and (d).

Close modal

The conformation space network constructed in that way is subsequently analyzed by a method termed “gradient-clustering,”51 which is a simpler version of the stochastic steepest descent algorithm introduced earlier.61 Gradient-clustering partitions the free-energy landscape into basins by building a reduced conformation space network [Figs. 2(c) and 2(d)] from the original (full) conformation space network schematically shown in Figs. 2(a) and 2(b). For each node, only the one link with the largest statistical weight is kept. Deleting all other links partitions the network into one or more fragments that are connected in a star-like manner [Figs. 2(c) and 2(d)], which are considered to be basins of the free-energy surface. The physical motivation for doing so is that nodes in the neighborhood of a barrier between two basins will maintain only the link to the valley, to which they relax the fastest. Following the most probable link results in a “steepest descent” on the free-energy surface, where microstates in the neighborhood of a barrier will connect to either one basin or the other.

Panels (e) and (f) of Fig. 2 trace the corresponding free energy surface for the two examples at low and high temperatures, which changes from one with only one minimum into one that evolves a side minimum. The split-off temperature is characterized by the limiting case between Figs. 2(e) and 2(f), i.e., by a free energy surface with an inflection point, and hence is sharply defined. It is important to point out that this sharply defined split-off temperature certainly does not reflect any phase transition, since the molecular system under study is finite, i.e., a hydrogen-bonded water cluster with two hydration layers and with a maximum size of 17 water molecules. The height of the barrier just above the split-off temperature is marginal, while it would be infinite for a phase transition in an infinite system.

We have tested that the results from the conformation space network are very robust against the details of the analysis. For example, we qualitatively obtain the same result for Figs. 3 and 4 (vide infra) when using either the positions of individual beads of the PIMD to evaluate whether or not there is a hydrogen bond via Eq. (1) or the centroid coordinates, despite the fact that the distribution of occupancies N deviates significantly in the two cases, see Figs. S1 and S2 (supplementary material). This is since the network analysis evaluates kinetic connectivity of microstates and as such averages out much of the thermal noise in the motion of individual beads. To this end, robustness is a strength of complex networks.62 

FIG. 3.

Population of the seven most populated gradient clusters for classical (solid lines) and quantum (dashed lines) simulations with (a) q-TIP4P/F and (c) TTM3-F water potentials as a function of temperature. Each gradient cluster can be characterized by one of the structures indicated in Fig. 1, which acts as an attractor of the corresponding free energy minimum (gradient cluster #1 in red, #2 in cyan, #3 in green, #4 in dark green, #5 in blue, #6 in violet, and #7 in magenta). Panels (b) and (d) show the same populations when adjusting the temperature of the classical simulation (given at the bottom) to that of the quantum simulation (given at the top) according to Eq. (2) with α = 0.99 and ΔT = −6 K for q-TIP4P/F and α = 1.24 and ΔT = −64 K for TTM3-F.

FIG. 3.

Population of the seven most populated gradient clusters for classical (solid lines) and quantum (dashed lines) simulations with (a) q-TIP4P/F and (c) TTM3-F water potentials as a function of temperature. Each gradient cluster can be characterized by one of the structures indicated in Fig. 1, which acts as an attractor of the corresponding free energy minimum (gradient cluster #1 in red, #2 in cyan, #3 in green, #4 in dark green, #5 in blue, #6 in violet, and #7 in magenta). Panels (b) and (d) show the same populations when adjusting the temperature of the classical simulation (given at the bottom) to that of the quantum simulation (given at the top) according to Eq. (2) with α = 0.99 and ΔT = −6 K for q-TIP4P/F and α = 1.24 and ΔT = −64 K for TTM3-F.

Close modal
FIG. 4.

Correlation between the temperatures at which the gradient clusters split-off (squares), obtained from classical vs. quantum simulations with the q-TIP4P/F (blue) and TTM3-F (red) interaction potentials. The switching points have been determined with an accuracy of ±1.25 K. The solid lines represent linear fits that establish the scaling law of Eq. (2). The numbers refer to the structures shown in Fig. 1. The dotted circle marks the point of equal stability between the classical and quantum descriptions of liquid water with the TTM3-F model.

FIG. 4.

Correlation between the temperatures at which the gradient clusters split-off (squares), obtained from classical vs. quantum simulations with the q-TIP4P/F (blue) and TTM3-F (red) interaction potentials. The switching points have been determined with an accuracy of ±1.25 K. The solid lines represent linear fits that establish the scaling law of Eq. (2). The numbers refer to the structures shown in Fig. 1. The dotted circle marks the point of equal stability between the classical and quantum descriptions of liquid water with the TTM3-F model.

Close modal

Panels (a) and (c) in Fig. 3 show the populations of the various gradient clusters obtained following the procedure outlined in Sec. II C for the q-TIP4P/F and TTM3-F interaction potentials as a function of the absolute temperature, overlaying the results from the classical (solid lines) and the quantum (dashed lines) simulations. In either case, only gradient cluster #1 is present with ∼80% probability at low temperatures; this cluster is represented by the fully hydrogen-bonded structure #1 (see Fig. 1), which acts as an attractor of a funnel-like free-energy surface [see Figs. 2(c) and 2(e)].52 As the temperature increases, other gradient clusters split off from gradient cluster #1 in a stepwise manner with various patterns of missing hydrogen bonds. While those split-off temperatures are in principle sharply defined [see panels (e) and (f) in Fig. 2], we find that due to statistical noise, the algorithm jumps back and forth in a narrow temperature range of typically ±0.3 K. In total seven gradient clusters can be identified at the highest temperature considered in this study (370 K), each characterized by one of the hydrogen-bonded structures shown in Fig. 1. These structures are “attractors” of the corresponding free energy basins. Note that many more microstates may exist in each of these gradient clusters, i.e., structures containing loops or more voids, which, however, are kinetically close to the representative structures. While the population of gradient cluster #1 continuously decreases with temperature, the populations of all other gradient clusters, once they are formed, increase (except when further split-offs occur, e.g., gradient cluster #6 splitting off from gradient cluster #5). Hence, the fully hydrogen-bonded gradient cluster #1 is enthalpically favored, while all other gradient clusters are entropically stabilized at higher temperatures. This overall behavior is identical to that previously observed during classical simulations for simpler water potentials such as SPC, TIP4P/2005, and TIP3P.51,52

The overall pattern of split-offs as a function of temperature is qualitatively similar for q-TIP4P/F vs. TTM3-F and also quite similar for the analysis of the classical vs. quantum trajectories for a given water model. Nevertheless, the particular temperatures, at which the split-offs occur, do differ, as seen from Fig. 4, which correlates the switch-off points of the classical vs. quantum simulation of the two water models. In either case, a high degree of correlation is found that is very close to the linear in the temperature range considered in this study. The correlation of switching points can be fit to the following temperature scaling law (Fig. 4, solid lines):

T(qm)=αT(class)+ΔT,
(2)

with α = 0.99 and ΔT = −6 K for q-TIP4P/F and α = 1.24 and ΔT = −64 K for TTM3-F. The correlation coefficients are very high with R2 = 0.998 for q-TIP4P/F and R2 = 0.996 for TTM3-F. When overlaying the results obtained during the classical vs. quantum simulations shown in panels (b) and (d) of Fig. 3, the temperature axes (bottom and top) are scaled according to Eq. (2), revealing an excellent agreement of the corresponding conformation space networks.

There is a general consensus that “classical” water is more stable than “quantum” water,9,10,48 as evidenced for example by the experimentally observed higher melting temperature of D2O (3.8 °C) and T2O (4.5 °C) as compared to H2O. In that regard, it is important to realize that the network analysis measures the amount of structuring with respect to hydrogen-bonded water network, and not stability per se, since the energy of a configuration is not evaluated by the algorithm. There is nevertheless a loose connection between the amount of structuring and stability, in the sense that structuring competes with thermal energy. That is, if a particular water model at a given temperature is more structured than another one, then this implies that it is also more stable, since it requires a higher temperature to destroy the structure.

The lower height of the first peak of the oxygen–oxygen radial distribution function (RDF) is consistent with the commonly accepted lower degree of water structuring in the quantum case. Figure 5 shows that this is indeed the case for both water models in the temperature range considered in this study; this is true both when plotting it in an absolute temperature (left panels) or when plotting it in the rescaled temperature according to Eq. (2). The general resemblance of the RDFs with the two models may imply that the properties of the two water models are similar in that temperature range.

FIG. 5.

Height of the first peak of the oxygen–oxygen radial distribution function (RDF) as a function of temperature for the q-TIP4P/F (top) and TTM3-F (bottom) potentials. Left/right panels show the results in an absolute/scaled temperature range [according to Eq. (2) with α = 0.99 and ΔT = −6 K for q-TIP4P/F and α = 1.24 and ΔT = −64 K for TTM3-F]. The red lines trace the results of the classical simulations, and the blue lines trace the ones of the quantum simulations.

FIG. 5.

Height of the first peak of the oxygen–oxygen radial distribution function (RDF) as a function of temperature for the q-TIP4P/F (top) and TTM3-F (bottom) potentials. Left/right panels show the results in an absolute/scaled temperature range [according to Eq. (2) with α = 0.99 and ΔT = −6 K for q-TIP4P/F and α = 1.24 and ΔT = −64 K for TTM3-F]. The red lines trace the results of the classical simulations, and the blue lines trace the ones of the quantum simulations.

Close modal

Figure 4 nevertheless suggests that according to the network analysis, the two water models behave qualitatively differently with respect to quantum effects. First, the scaling laws lie on either side of the diagonal for the most part of the temperature range considered in this study. In that plot, the diagonal traces equal amount of structuring between the classical and quantum water (assuming we can equate stability with structuring). Furthermore, there is only a constant temperature offset with the slope α being practically 1 in the case of q-TIP4P/F, while the slope α significantly deviates from 1 for TTM3-F (α = 1.24). Consequently, in the case of TTM3-F, there is a point of equal stability at T ≈ 267 K, marked by the dotted circle in Fig. 4. Above that temperature (which covers most of the liquid phase) quantum TTM3-F water is more structured than classical TTM3-F water in contrast to the common expectation. Classical q-TIP4P/F water, on the other hand, is always more structured than quantum q-TIP4P/F. It is suggestive that the (essentially constant) temperature difference between the quantum and classical q-TIP4P/F water (∼9 °C at 300 K) agrees extremely well with the difference of 8 K calculated for their melting points.6 However, a more systematic study, considering more than one water model, would be needed to confirm that there is indeed a correlation between the melting point, which is a true phase transition, and the split-off points found in the network analysis.

Some discussion regarding the behavior of the TTM3-F interaction potential is in order. Unlike its previous version (v. 2.1), which predicted a contribution of ∼1.0 kcal/mol to the enthalpy of liquid water at T = 298 K due to quantum effects,32 version 3.0 in contrast predicts a negligible correction of ∼0.1 kcal/mol.31 We have previously discussed31 this difference in terms of the different intramolecular potential energy surface between the two models, whereas Manolopoulos and co-workers6 have also elaborated on the factors controlling this difference. There are two competing factors that determine the properties of liquid water during a classical simulation and a quantum simulation: one is the quantum fluctuations that lead to faster diffusion, orientation dynamics, and less structured RDFs, and the other is the magnitude of the various interactions (mainly electrostatic), which during quantum simulations are stronger and have exactly the opposite trend for the previous properties. The latter is due to the longer intramolecular OH bonds during a quantum simulation, a fact that leads to a higher dipole moment and therefore stronger electrostatic interactions compared to the ones during a classical simulation. Assuming that the magnitude of the quantum fluctuations is similar for all models, the electrostatic interactions have a large effect on the magnitude of quantum corrections. Our simulations with the TTM3-F potential suggest that at low temperatures the individual monomer dipole moments are 2.767 D (classical at T = 287.5 K) and 2.876 D (quantum at T = 285 K), whereas at higher temperatures they are 2.710 D (classical at T = 350 K) and 2.826 D (quantum at T = 350 K). Therefore, over a wide temperature range, the difference between the fragment dipole moments obtained via quantum and classical simulations remains practically the same (within ∼0.11 D). This value is quite high and is the result of the non-linear dipole moment surface (DMS) that the TTM3-F model incorporates, according to which the charge on the hydrogen atom increases as the OH bond elongates. In contrast, for a typical model that instead uses a linear DMS (i.e., one with constant charges) such as TIP4P, qTIP4P/F, SPC and others, the difference in the fragment dipole moment obtained via quantum and classical simulations is much smaller (0.037 D for the qTIP4P/F).6 Therefore, the TTM3-F potential yields larger differences in the electrostatic interactions between quantum and classical simulations when compared to any other typical water model with constant charges (i.e., a linear monomer DMS). The difference in the electrostatic interactions manifests itself in the diffusion coefficient with estimates for the ratio Dquant./Dclass. being 1.15 for the qTIP4P/F potential33 and 0.95 for the TTM3-F potential.6 Version 2.1 of the TTM model also incorporates a non-linear DMS in which, in contrast to the DMS of version 3.0, the charge on the hydrogen atom decreases as the OH bond elongates and predicts, as expected, a higher value for Dquant./Dclass. = 1.5.33 

In a recent paper, Markland and Berne48 investigated the isotopic fractionation between liquid and vapor phase, a measure that can determine quantum effects experimentally at any temperature (in contrast to the freezing point). Experimentally, the “more classical” HOD molecule in H2O prefers the liquid phase with a certain free energy difference at temperatures below 500 K, reflecting the more stable hydrogen bonding of HOD in H2O.63 This result can almost quantitatively be reproduced with the q-TIP4P/F model, while the TTM3-F potential predicts the opposite at 300 K,48 presumably since the nonlinear DMS overcompensates quantum delocalization and thus leads to more stable hydrogen bonds in the quantum case.

Therefore, the previously reported results of the isotopic fractionation and the diffusion constant suggest that the q-TIP4P/F and TTM3-F interaction potentials behave qualitatively differently with respect to quantum effects. The qualitative agreement of the conclusions that can be drawn from these two properties, viz., the opposite behavior between q-TIP4P/F and TTM3-F, suggests that the analysis of the conformation space network, which we introduce here as a new method to study the structure of water as a function of temperature, is indeed a valid approach.

We have shown that the hydrogen bond structural networks in liquid water obtained from the quantum and the classical descriptions of the nuclei with two different interaction potentials remain very similar throughout the whole temperature range investigated (T = 250-370 K), when rescaling the temperature according to Eq. (2). In contrast to the general consensus, the analysis of the structural networks suggests that a classical description of water is not necessarily more structured than the quantum one, depending on the force field; rather an inversion of regimes may be found above a temperature of equal stability. That inversion of regimes is attributed to a counterbalance of quantum effects. Quantum fluctuations render quantum water more disordered, while on the other hand the on-the-average larger OH bond length increases the water dipole and thus induces order. The TTM3-F potential overestimates the intra-molecular component and hence the inversion of regimes is found at just above freezing point, while experimentally an inversion of regimes is found only around 500 K, as judged from the isotopic fractionation.48,63 To this end, our results can provide important guidelines for the development of interaction potentials to capture quantum effects correctly.

The present study nonetheless shows that the hydrogen bonded networks obtained via classical and quantum simulations are structurally very similar over the examined temperature range even for interaction potentials that are based on different philosophies and have been reported to describe quantum effects quite differently. Once a relationship similar to Eq. (2) is established for a given force field (the parameters being model-dependent), quantum effects, which inevitably need to be added to any classical, ab initio based water potential as well as DFT-based MD simulations, can be incorporated by such a simple rescaling law in a computationally extremely efficient manner. Such rescaling approaches have of course been discussed before, as extensively described in the Introduction, but the method proposed here, which resolves the free-energy surface of water by identifying complex hydrogen-bond patterns, can provide a refined and structurally better-founded scaling law. This can ultimately lay the ground to accurately account for larger parts of the phase diagram of liquid water using more sophisticated and realistic water interaction potentials.

See supplementary material for Fig. S1, comparing the distribution of occupancy N calculated from beads vs. the centroid in the PIMD, as well as Fig. S2, comparing the resulting correlations between classical and quantum water.

We thank Francesco Rao for his contributions at an early stage of the project. Part of this work was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences (for S.S.X.), and by the Swiss National Science Foundation (SNF) through the NCCR MUST (for P.H.). Pacific Northwest National Laboratory (PNNL) is a multiprogram national laboratory operated for DOE by Battelle. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

1.
J. D.
Bernal
and
R. H.
Fowler
,
J. Chem. Phys.
1
,
515
(
1933
).
2.
J. A.
Barker
and
R. O.
Watts
,
Chem. Phys. Lett.
3
,
144
(
1969
).
3.
A.
Rahman
and
F. H.
Stillinger
,
J. Chem. Phys.
55
,
3336
(
1971
).
5.
D.
Chandler
,
Introduction to Modern Statistical Mechanics
(
Oxford University Press
,
Oxford
,
1987
).
6.
S.
Habershon
,
T. E.
Markland
, and
D. E.
Manolopoulos
,
J. Chem. Phys.
131
,
024501
(
2009
).
7.
J. L. F.
Abascal
,
E.
Sanz
,
R.
García Fernández
, and
C.
Vega
,
J. Chem. Phys.
122
,
234511
(
2005
).
8.
C.
Millot
,
J.-C.
Soetens
,
M. T. C. M.
Costa
,
M. P.
Hodges
, and
A. J.
Stone
,
J. Phys. Chem. A
102
,
754
(
1998
).
9.
R. A.
Kuharski
and
P. J.
Rossky
,
J. Chem. Phys.
82
,
5164
(
1985
).
10.
H. A.
Stern
,
F.
Rittner
,
B. J.
Berne
, and
R. A.
Friesner
,
J. Chem. Phys.
115
,
2237
(
2001
).
11.
M. W.
Mahoney
and
W. L.
Jorgensen
,
J. Chem. Phys.
115
,
10758
(
2001
).
12.
J.
Cao
and
G. A.
Voth
,
J. Chem. Phys.
99
,
10070
(
1993
).
13.
J.
Cao
and
G. A.
Voth
,
J. Chem. Phys.
100
,
5106
(
1994
).
14.
J.
Cao
and
G. A.
Voth
,
J. Chem. Phys.
101
,
6168
(
1994
).
15.
J.
Cao
and
G. A.
Voth
,
J. Chem. Phys.
101
,
6157
(
1994
).
16.
J.
Cao
and
G. J.
Martyna
,
J. Chem. Phys.
104
,
2028
(
1996
).
17.
G. J.
Martyna
,
M. L.
Klein
, and
M. E.
Tuckerman
,
J. Chem. Phys.
97
,
2635
(
1992
).
18.
T. F.
Miller
and
D. E.
Manolopoulos
,
J. Chem. Phys.
123
,
184503
(
2005
).
19.
S.
Habershon
,
D. E.
Manolopoulos
,
T. E.
Markland
, and
T. F.
Miller
,
Annu. Rev. Phys. Chem.
64
,
387
(
2013
).
20.
K.
Laasonen
,
M.
Sprik
, and
M.
Parrinello
,
J. Chem. Phys.
99
,
9080
(
1993
).
21.
B.
Chen
,
I.
Ivanov
,
M. L.
Klein
, and
M.
Parrinello
,
Phys. Rev. Lett.
91
,
215503
(
2003
).
22.
J. A.
Morrone
and
R.
Car
,
Phys. Rev. Lett.
101
,
017801
(
2008
).
23.
X.
He
,
O.
Sode
,
S. S.
Xantheas
, and
S.
Hirata
,
J. Chem. Phys.
137
,
204505
(
2012
).
24.
C. J.
Burnham
,
J.
Li
,
S. S.
Xantheas
, and
M.
Leslie
,
J. Chem. Phys.
110
,
4566
(
1999
).
25.
C. J.
Burnham
and
S. S.
Xantheas
,
J. Chem. Phys.
116
,
1479
(
2002
).
26.
G. S.
Fanourgakis
and
S. S.
Xantheas
,
J. Phys. Chem. A
110
,
4100
(
2006
).
27.
S. S.
Xantheas
,
C. J.
Burnham
, and
R. J.
Harrison
,
J. Chem. Phys.
116
,
1493
(
2002
).
28.
S.
Bulusu
,
S.
Yoo
,
E.
Aprà
,
S.
Xantheas
, and
X. C.
Zeng
,
J. Phys. Chem. A
110
,
11781
(
2006
).
29.
G. S.
Fanourgakis
,
E.
Aprà
, and
S. S.
Xantheas
,
J. Chem. Phys.
121
,
2655
(
2004
).
30.
H.
Partridge
and
D. W.
Schwenke
,
J. Chem. Phys.
106
,
4618
(
1997
).
31.
G. S.
Fanourgakis
and
S. S.
Xantheas
,
J. Chem. Phys.
128
,
074506
(
2008
).
32.
G. S.
Fanourgakis
,
G. K.
Schenter
, and
S. S.
Xantheas
,
J. Chem. Phys.
125
,
141102
(
2006
).
33.
F.
Paesani
,
S.
Iuchi
, and
G. A.
Voth
,
J. Chem. Phys.
127
,
074506
(
2007
).
34.
F.
Paesani
,
S.
Yoo
,
H. J.
Bakker
, and
S. S.
Xantheas
,
J. Phys. Chem. Lett.
1
,
2316
(
2010
).
35.
J.
Liu
,
W. H.
Miller
,
G. S.
Fanourgakis
,
S. S.
Xantheas
,
S.
Imoto
, and
S.
Saito
,
J. Chem. Phys.
135
,
244503
(
2011
).
36.
R. S.
Fellers
,
C.
Leforestier
,
L. B.
Braly
,
M. G.
Brown
, and
R. J.
Saykally
,
Science
284
,
945
(
1999
).
37.
C.
Leforestier
,
F.
Gatti
,
R. S.
Fellers
, and
R. J.
Saykally
,
J. Chem. Phys.
117
,
8710
(
2002
).
38.
G. R.
Medders
,
V.
Babin
, and
F.
Paesani
,
J. Chem. Theory Comput.
9
,
1103
(
2013
).
39.
V.
Babin
,
G. R.
Medders
, and
F.
Paesani
,
J. Phys. Chem. Lett.
3
,
3765
(
2012
).
40.
V.
Babin
,
C.
Leforestier
, and
F.
Paesani
,
J. Chem. Theory Comput.
9
,
5395
(
2013
).
41.
J. L. F.
Abascal
and
C.
Vega
,
J. Chem. Phys.
123
,
234505
(
2005
).
42.
E. G.
Noya
,
C.
Vega
,
L. M.
Sesé
, and
R.
Ramírez
,
J. Chem. Phys.
131
,
124518
(
2009
).
43.
L.
Hernández de la Peña
and
P. G.
Kusalik
,
J. Chem. Phys.
125
,
054512
(
2006
).
44.
F.
Paesani
and
G. A.
Voth
,
J. Phys. Chem. C
112
,
324
(
2008
).
45.
F.
Paesani
,
S. S.
Xantheas
, and
G. A.
Voth
,
J. Phys. Chem. B
113
,
13118
(
2009
).
46.
T. E.
Markland
,
J. A.
Morrone
,
B. J.
Berne
,
K.
Miyazaki
,
E.
Rabani
, and
D. R.
Reichman
,
Nat. Phys.
7
,
134
(
2011
).
47.
O.
Mishima
and
H. E.
Stanley
,
Nature
396
,
329
(
1998
).
48.
T. E.
Markland
and
B. J.
Berne
,
Proc. Natl. Acad. Sci. U. S. A.
109
,
7988
(
2012
).
49.
C.
McBride
,
E. G.
Noya
,
J. L.
Aragones
,
M. M.
Conde
, and
C.
Vega
,
Phys. Chem. Chem. Phys.
14
,
10140
(
2012
).
50.
C.
Vega
and
J. L. F.
Abascal
,
J. Chem. Phys.
123
,
144504
(
2005
).
51.
F.
Rao
,
S.
Garrett-Roe
, and
P.
Hamm
,
J. Phys. Chem. B
114
,
15598
(
2010
).
52.
D.
Prada-Gracia
,
R.
Shevchuk
,
P.
Hamm
, and
F.
Rao
,
J. Chem. Phys.
137
,
144504
(
2012
).
53.
H. C.
Andersen
,
J. Chem. Phys.
72
,
2384
(
1980
).
54.
H. J. C.
Berendsen
,
J. P. M.
Postma
,
W. F.
van Gunsteren
,
A.
DiNola
, and
J. R.
Haak
,
J. Chem. Phys.
81
,
3684
(
1984
).
55.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulations of Liquids
(
Oxford University Press
,
Oxford
,
1987
).
56.
T. M.
Nymand
and
P.
Linse
,
J. Chem. Phys.
112
,
6152
(
2000
).
57.
G. S.
Fanourgakis
,
T. E.
Markland
, and
D. E.
Manolopoulos
,
J. Chem. Phys.
131
,
094102
(
2009
).
58.
G. S.
Fanourgakis
,
V.
Tipparaju
,
J.
Nieplocha
, and
S. S.
Xantheas
,
Theor. Chem. Acc.
117
,
73
(
2007
).
59.
R.
Kumar
,
J. R.
Schmidt
, and
J. L.
Skinner
,
J. Chem. Phys.
126
,
204107
(
2007
).
60.
J.-H.
Prinz
,
H.
Wu
,
M.
Sarich
,
B.
Keller
,
M.
Senne
,
M.
Held
,
J. D.
Chodera
,
C.
Schütte
, and
F.
Noé
,
J. Chem. Phys.
134
,
174105
(
2011
).
61.
D.
Prada-Gracia
,
J.
Gómez-Gardeñes
,
P.
Echenique
, and
F.
Falo
,
PLoS Comput. Biol.
5
,
e1000415
(
2009
).
62.
R.
Albert
,
H.
Jeong
, and
A.-L.
Barabasi
,
Nature
406
,
378
(
2000
).
63.
J.
Horita
and
D. J.
Wesolowski
,
Geochim. Cosmochim. Acta
58
,
3425
(
1994
).

Supplementary Material