New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology
Liquid water is a dynamic polydisperse branched polymer
Contributed by William A. Goddard III, December 17, 2018 (sent for review October 9, 2018; reviewed by Charles L. Brooks III and Michael L. Klein)
This article has a Letter. Please see:
- Water is not a dynamic polydisperse branched polymer - June 25, 2019
See related content:

Significance
In contrast to ice, in which each water makes strong hydrogen bonds (SHBs) to four neighbors, we show that upon melting the number of SHBs drops quickly to two in liquid water. These two SHBs couple into chains containing ∼150 waters resembling a branched polymer. The lifetime of each SHB at 298 K is 90.3 fs (11 OH vibrational periods), so the polymer branches evolve dynamically. This dynamics-branched polymer paradigm may explain long-standing puzzles of water, such as the critical point at 227 K in supercooled water (which may correspond to a glass transition caused by an increase in the SHB lifetime). It may explain the observed angular correlations in water that persist for 20 nm.
Abstract
We developed the RexPoN force field for water based entirely on quantum mechanics. It predicts the properties of water extremely accurately, with Tmelt = 273.3 K (273.15 K) and properties at 298 K: ΔHvap = 10.36 kcal/mol (10.52), density = 0.9965 g/cm3 (0.9965), entropy = 68.4 J/mol/K (69.9), and dielectric constant = 76.1 (78.4), where experimental values are in parentheses. Upon heating from 0.0 K (ice) to 273.0 K (still ice), the average number of strong hydrogen bonds (SHBs, rOO ≤ 2.93 Å) decreases from 4.0 to 3.3, but upon melting at 273.5 K, the number of SHBs drops suddenly to 2.3, decreasing slowly to 2.1 at 298 K and 1.6 at 400 K. The lifetime of the SHBs is 90.3 fs at 298 K, increasing monotonically for lower temperature. These SHBs connect to form multibranched polymer chains (151 H2O per chain at 298 K), where branch points have 3 SHBs and termination points have 1 SHB. This dynamic fluctuating branched polymer view of water provides a dramatically modified paradigm for understanding the properties of water. It may explain the 20-nm angular correlation lengths at 298 K and the critical point at 227 K in supercooled water. Indeed, the 15% jump in the SHB lifetime at 227 K suggests that the supercooled critical point may correspond to a phase transition temperature of the dynamic polymer structure. This paradigm for water could have a significant impact on the properties for protein, DNA, and other materials in aqueous media.
- water structure
- molecular dynamics
- liquid–liquid critical point
- radial distribution function
- density-functional theory
Due to its importance for life and the environment, the properties of liquid water have been studied thoroughly by numerous experimental and computational techniques. Even so, bulk water exhibits many anomalies that challenge understanding (1, 2). For example, supercooled water exhibits a critical point at 227 K (3) and polarization-resolved second-harmonic scattering experiments (4) show that water orientations are correlated over length scales of 20 nanometers (nm).
Our understanding of ice is excellent, with the many observed phases quite consistent with each H2O making strong hydrogen bonds to each of its four neighbors [O–O distance in ice Ih ∼ 2.75 Å at 1 atm (5)]. For liquid water the only direct experimental structural information is the oxygen–oxygen radical distribution function, gOO, from neutron scattering, shown in Fig. 1. Here the first peak from 2.5 to 3.4 Å describes the first shell of waters around each particular water molecule.
Comparison of first peak in the oxygen–oxygen radial distribution function (gOO) predicted from RexPoN (black) with experiments and other models at T = 298 K and 1-atm pressure. The gOO for RexPoN was computed by averaging the trajectory every 1 ps during 1-ns MD-NVT simulations for a periodic cell containing 216 water molecules. (A) Comparison of the gOO of RexPoN with the neutron scattering experiment (19) (red). RexPoN leads to a first O–O peak at R1 = 2.84 Å and gOO(R1) = 2.34, in excellent agreement with neutron at R1 = 2.86 Å and gOO(R1) = 2.50. The computed number of SHBs is ∑SHB = 2.14 for rOO ≤ 2.93 Å (dark-blue area) with an average lifetime (τ) of 90.3 fs. RexPoN leads to a first minimum of gOO at 3.43 Å with a total coordination number Noo = 4.7, in good agreement with experiment (∼4.5). (B) Comparison of experimental gOO curves obtained by combining of X-ray and neutron scattering experiments [X-ray1 (20), X-ray2 (21), and X-ray3 (22)] with the neutron gOO. To remove the contribution from electrons (i.e., O–H interaction and oxygen lone-pair distortions), these X-ray results use various empirical FFs. Therefore, the most reliable gOO curve is neutron where there is no inference from electrons. (C) gOO curves obtained by QM MD with PBE (23), PBE0+Tkatchenko and Scheffler (TS)-vdW (24), PBE-D3 (25), and SCAN (26) flavors of DFT. The first peak and minimum of the DFT differ dramatically from experiment, explaining why they lead to very inaccurate structural properties for water. (D) gOO computed using the CC-pol (27) and MB-pol (28) FF based on QM calculations. These models lead to sharp and high first peak, indicating that they accommodate too many water molecules in the first shell, leading to too many SHBs. Thus, NOO = 5.6 for CC-pol and 6.2 for MB-pol. (E) gOO from the popular empirical TIP4P-2005f (29), SPC/E (30), and ST2 (31) FFs. These models also lead to a very sharp and high first peak compared with neutron data. All empirical FFs fitted to the gOO assume ΣSHB = 4, leading to a very sharp peak at ∼2.77 Å with a maximum as high as 3.3 instead of 2.5.
Thus, all of our detailed structural information about liquid water is derived from computer simulations. Here one must postulate a model potential or force field (FF) that is then used for molecular-dynamics (MD) simulations to predict the properties. The parameters in these FFs are generally adjusted to reproduce two or three of such experimental properties as density, heat of vaporization, diffusion, and/or melting point. However, Fig. 1 shows that all previous FFs and quantum mechanics (QM)-derived MD (QM MD) fail to reproduce the experimental description of the first shell of waters. Thus, all rise too fast at small R, go too high at the first peak at too short a distance, and then fall too fast to the first minimum. This suggests that all previous FF are missing something important about the H2O–H2O interactions in the first shell of neighbors.
This dramatic disagreement motivated us to develop the RexPoN FF for water, based entirely on QM with no empirical information. There are four components to RexPoN aimed at independently describing the short- and long-range interactions:
i) Nonbond (NB) or van der Waals (vdW) (Pauli repulsion and London dispersion) based on fitting the density-functional theory (DFT) equation of state (EOS) for H2 and O2, which after including phonon corrections matches well the experimental EOS. This ensures an accurate description of long-range vdW interactions.
ii) Polarization and charge optimization based on the polarizable charge equilibration (PQEq) (6, 7) model with finite-sized charges and shells (not point charges). This has been validated to reproduce the QM polarization of various molecules as point dipoles are brought in from various directions.
iii) Hydrogen bond (HB) corrections with explicit lone pairs based on extensive coupled-cluster single-double triple, CCSD(T), ab initio QM on water dimers (8).
iv) Reactive HO bond and HOH angles based on CCSD(T) ab initio QM to describe exactly the bond breaking (8).
We denote this FF as RexPoN to indicate the use of Reax, PQEq, and NB. The mathematical formulation of the FF and details of the parameters for each energy term are provided elsewhere (9).
Results
Structural and Thermodynamics Properties.
Table 1 compares the predicted properties using RexPoN with experiment and with QM MD and other FFs. RexPoN predicts:
i) A melting point (Tmelt) between 273 and 273.5 K, very close to the experimental value of 273.15 K (10). This can be compared with Tmelt = 146 K for transferable intermolecular potential with three points (TIP3P), 215 K for extended simple point charge model (SPC/E), and 420 K for Perdew–Burke–Ernzerhof (PBE).
ii) A density of 0.9965 g/cm3 at 298 K, agreeing exactly with experiment (it also captures the maximum density of 1.000 g/cm3 at 4 °C, where underline indicates some uncertainty). This can be compared with 0.944 g/cm3 for PBE DFT and 1.05 for strongly constrained and appropriately normed (SCAN) DFT.
iii) A heat of vaporization of 10.36 kcal/mol, low by 0.16 kcal/mol from experiment. This can be compared with 6.2 kcal/mol for PBE DFT.
iv) An entropy at 298 K of 68.4, low by 2%. This can be compared with 51.32 for PBE DFT.
v) A dielectric constant at 298 K of 76.1, low by 3%. This can be compared with 112 for PBE DFT and 68.4 for MB-pol (based on QM).
Of most interest is the gOO from RexPoN shown in Fig. 1. The predicted gOO for the first peak (up to 3.43 Å) is in excellent agreement with experiment with a slight shift in the peak position from 2.86 to 2.84 Å and the peak height from 2.50 to 2.34. In addition, the computed coordination number, NOO = 4.7, is in good agreement with experiment, NOO = 4.5. Fig. 1 shows that no other FF or QM description captures correctly the shape and height of this first peak in the gOO.
The accurate dielectric constant from RexPoN shows that the polarization effects and dipole moments for various molecular configuration in bulk water are captured correctly. This is crucial for describing the arrangements and alignments of water molecules in the first shell that in turn determine the strong hydrogen bonds (see below).
We consider that the spectacular agreement with experiment of this purely QM description (no empirical data) indicates that RexPoN captures accurately the properties of water.
We consider these excellent predictions to be a strong validation of the accuracy of RexPoN FF for describing the properties of ice and water.
Melting Simulations of Ice Ih.
To understand how ice evolves to water from 273.0 to 273.5 K, we examined the melting of ice, starting with ice Ih (96 molecules per cell). In this analysis we heated the system slowly from 0 to 273.5 K while keeping a pressure of 1 atm and followed the initial steps of melting. Since H2O dimer leads to a strong hydrogen bond (SHB) with a well depth of 4.98 kcal/mol at rOO = 2.93 Å, we define a SHB in bulk ice or water as having an OO distance ≤2.93 Å.
For ice, the sum of the SHBs to each water is ΣSHB = 4.0 at 0 K, decreasing slowly to 3.3 by 273 K, due to vibrations.
Fig. 2 shows that upon heating to 273.5, the ΣSHB drops quickly to 3.3 during the first 10 ps. The ΣSHB stays around 3.3 for 150 ps of MD at the density of the ice (MD-NVT). However, at 142 ps the initial steps of melting initiate and melting is completed in just 33 ps (at 175 ps) of additional MD simulations allowing the density to optimize (MD-NPT). First (step 1), the blue and green waters molecules break partly out of the same ice hexagon, reducing ΣSHB to ∼2. Then (step 2), for ∼20 ps these two molecules go in and out of their hexagon with ΣSHB ranging from 1.75 to 3.3. This is followed (step 3) by weakening of the initial hexagon and nearby ones for ∼10 ps, which by 175 ps (step 4) leads to full melting and ΣSHB = 2.3. Movie S2 shows this process.
Illustration of the initial steps of melting of Ice-Ih at 273.5 K. The system with a cell of 96 water molecules was heated from 0 to 273.5 K over 10 ps followed by MD-NVT equilibration for 150 ps. Then the volume was allowed to adjust to 1 atm during 100 ps of MD-NPT simulations. (A) The change of the ΣSHB versus time for all 96 molecules (red). Also shown are ΣSHB for water molecules number 70 (green) and 71 (blue). All data were averaged over 1-ps time intervals for better visibility. The initiation of melting is dominated by the blue and green waters in the same hexagonal ring of ice, involving four distinct steps (S1 to S4): at 142 ps (S1) the blue and green neighboring molecules (71 and 70) start to break out of the same hexagonal ring, moving toward the center of the ring with their ΣSHB both dropping from 3.3 to 2.6 at 152 ps. By 158 ps, ΣSHB has dropped to 1.75 for blue while for green it increased back to 3.1. From 142 to 162 ps (S2) the blue and green molecules move in and out of the hexagonal ring repeatedly. Then, from 162 to 175 ps (S3) the local hexagons near these two molecules collapse. Finally, the melting propagates quickly with the whole structure melting to ΣSHB ∼ 2.25 by 175 ps (S4). SI Appendix, Fig. S14 shows the MD trajectory for the case of T = 273.0 K, where we see that there is no melting over this same period of 260 ps. (B) Snapshots of the system corresponding to the four steps shown in A, with molecules 71 (blue) and 70 (green) highlighted. The remaining molecules are gray.
Structure of Water.
For 298 K, we equilibrated the bulk system in a cubic cell (a = 18.65 Å) containing 216 water molecules at 298 K and 1 atm (see ref. 11 for the details). This led to a density of 0.9965 g/cm3, in exact agreement with experiment and an average ΣSHB = 2.14. We were shocked by these unexpected results of only two strong HBs at 298 K.
Since each H2O on the average makes 2 SHBs, we next traced out the connections between the SHBs. We find that connecting just the SHBs leads to branched polymer chains with a largest cluster of ∼151 waters at 298 K, containing a main chain of 39 H2O and 15 side chains ranging from 1 to 22 H2Os long. This largest average cluster decreases to 72 at 350 K and to 36 at 400 K. For lower temperatures the largest cluster size increases to 168 at 277 K (the density maximum), 177 at 273.5 K (the initial melt), and to 216 at 150 K. Our calculations using 216 molecules per periodic cell might underestimate the true polymer sizes. Thus, at 298 K and below these clusters may bond to their images in adjacent cells, leading to infinite sizes (∞) in our simulations.
To study the structure of liquid water at other temperatures, we started with the equilibrated system at 298 K. Then, we heated the system separately to 350 and 400 K and cooled separately to various temperatures from 298 to 150 K during 10-ps MD-NVT simulations. Each system was then relaxed using constant temperature–constant pressure MD-NPT simulations at 1 atm for 400 ps followed by 100-ps MD-NVT equilibrations. We used the trajectories from the last 10 ps to analyze the structures.
Table 2 shows that ΣSHB drops dramatically from ΣSHB = 4.0 of ice at 0 K to ΣSHB = 2.28 at 273.5 K (liquid), and then decreases slowly to ΣSHB = 2.14 at 298 K and ΣSHB = 1.6 at 350 K. At 298 K we find one strong donor (SD) and one strong acceptor, and similar to ΣSHB the SD number drops at higher temperatures (0.71 at 400 K) while increasing at lower temperatures (1.55 at 150 K).
Summary of polymer properties as a function of temperature for the bulk system with 216 H2O molecules
Table 2 lists the average lifetime of each SHB, which is 90.3 fs at 298 K, increasing to 175.1 fs at 150 K and decreasing to 67.9 fs at 400 K. This can be compared with an OH vibrational period of 8.4 fs, so that the SHB should be measurable with femtosecond and attosecond spectroscopy.
The percent populations of ΣSHB (1⇓⇓⇓–5) as a function of temperature are shown in Fig. 3C. Thus, liquid water at 298, 277, and 273 K is a polydisperse dynamic multibranched polymer in which most H2Os form SHBs to just two others, with occasional branch points at waters bonded to 3 H2Os terminating at H2Os with ΣSHB = 1 at the end of the branches, as shown in Fig. 3 A and B. The molecules that are not part of the largest cluster are mainly in small clusters with sizes ranging from 1 to 10 H2Os. The average number of H2Os in these small clusters is 65 (30%) at 298 K decreasing to 39 (18%) at 273.5 K. At 350 K the big polymer chains split into smaller ones, with the number of H2Os bonded to two other molecules decreasing to 35% while the number of H2Os bonded to only one water molecule is 33% and 17% have ΣSHB = 3 (Fig. 3C). At 350 and 400 K, many more clusters have sizes up to 10 H2Os. There are an average of 144 and 180 water molecules in these small clusters at 350 and 400 K, respectively. At 150 K, the connected polymer includes all water molecules, which forms a highly branched polymer with mostly ΣSHB = 3 and 4.
Structure and HB network of water. (A) Topology of the largest H2O polymer for a single timeframe at 298 K and 1 atm. The average of all trajectories over 10 ps leads to a size of 151 for the largest water polymer cluster at 298 K. This structure was selected from the last 10-ps trajectories of MD-NVT simulations of a bulk cell containing 216 waters. It contains 151 water molecules connected via SHBs. The arrows point from the donor to the acceptor. The polymer backbone (longest chain, shown as red) contains 39 molecules. The side branches include: 1 chain with 22 molecules (silver), 1 chain with 20 molecules (yellow), 1 chain with 18 molecules (purple), 1 chain with 14 molecules (green), 1 chain with 11 molecules (orange), 1 chain with 6 molecules (blue), 1 chain with 4 molecules (tan), 3 chains with 3 molecules (pink), 3 chains with 2 molecules (magenta), and 2 chains with 1 molecules (cyan). Movie S1 shows this polymer from a variety of angles. SI Appendix, Figs. S8 and S9 show the same polymer chain after an additional 1, 10, 100, and 1,000 fs to provide a sense of the dynamical structure. Table 2 provides the polymer size at other temperatures. (The system was initially relaxed and equilibrated during 0.5-ns MD-NPT and -NVT simulations) (B) Same structure as in A from the same view but showing the OH bonds. The black lines show the SHBs. (C) The percent population for molecules with various ΣSHB as a function of temperature. The maximum population at 150 K (41%) is for ΣSHB ∼ 3.5, which drops to ΣSHB ∼ 2 at 273.5 K (36%). For temperatures of 277 and 298 K the maximum stays at ΣSHB ∼ 2 with about 36.0% population. At 350 K, the maximum population is ΣSHB ∼ 1.8 (36%). By 400 K the maximum is for ΣSHB ∼ 1.3. For a branched polymer chain each molecule with ΣSHB = 3 is a branch point, which generally terminates with ΣSHB = 1 so that these numbers are comparable at 298 K. Also shown are results at 298 K for a 4-nm-thick periodic slab (432 molecules) showing a maximum at ΣSHB ∼ 2 similar to bulk at 298 K. For each temperature the system was heated or cooled from the 298 K to the desired temperature over 10 ps. The system was then relaxed using MD-NPT simulations for 400 ps followed by 100-ps MD-NVT equilibrations. The ΣSHB were averaged over the last 10 ps of simulations for each temperature. The data for the slab were taken from the last 10 ps of MD-NVT simulations for 200 ps at 298 K. SI Appendix, Table S2 provides the slab ΣSHB populations at other temperatures.
Surface Evaporation.
We also constructed a 4 nanometer (nm) thick 2D slab with the same a and b parameters as the cubic cell (a = b = 18.64 Å) but doubled its size to 37.28 Å in the nonperiodic direction to form a system with 432 water molecules. The slab was heated from 0 K to various temperatures during 10 ps followed by 200 ps relaxation at the final temperatures using MD-NVT simulations. This leads to cluster sizes with ΣSHB only slightly smaller than the bulk at the same temperatures (see Fig. 3C and SI Appendix, Fig. S11 and Table S2 for 298 K). To examine evaporation, we tracked the molecules that came off the surface of this finite slab at various temperatures. We did not see any evaporation until 373.5 K (boiling point of water). We found that over a period of 200 ps, 11 H2O evaporate into the vacuum. All 11 were termini (i.e., ΣSHB = 1) of short polymer chains ranging in size from 2 to 5 H2Os. No waters evaporated from the larger clusters. Snapshots of one case is shown in SI Appendix, Fig. S13 (here the evaporated molecule was an acceptor).
Definition of the Strong Hydrogen Bond.
We defined a strong hydrogen bond (SHB) as ≤2.93 Å, the equilibrium distance for water dimer, which is bound by 4.98 kcal/mol. This can be compared with the 10.52 kcal/mole heat of vaporization of water at 298 K. However, the interpretation of water as a dynamic branched polymer chain is not sensitive to the definition of a SHB. Thus, section A.6 of the SI Appendix reanalyzes the dynamics using an OO distance ≤3.08 Å (allowing an extra 0.15 Å for zero-point motion from 2.93 Å) as the criterion, which leads to essentially the same results.
Discussion
An important question here is why we and all other theorists missed discovering the polymer nature of water earlier. The answer is that we all assumed that water makes four hydrogen bonds just as in ice, so our FFs were all biased toward this concept. The hint that there was a problem has, in fact, been apparent since the neutron scattering experiments in 1982. Previous FFs have not been able to capture the shape and critical parameters of the first peak in the neutron-derived gOO. RexPoN reproduces this shape fully.
Thus, Fig. 1 shows the gOO of the empirical FFs rise far too sharply at small R, go far too high with a peak magnitude up to 3.3 compared with 2.5 from experiment, leading to a peak slightly too short (∼2.77 Å, 0.1 Å too short), then fall far too quickly. This results in too many oxygens within the first neighbor shell.
Indeed, Fig. 1 and SI Appendix, section A.1 show that the same problem exists with all previous QM MD and QM-based FFs. The QM MD methods find a first peak that is dramatically different from experiment. The QM-based FFs also result in a sharp high peak, which increases NOO within the first neighbor shell [NOO = 5.6 for pair polarizable analytical potential (CC-pol) and NOO = 6.2 for many body potential with polarization (MB-pol) compared with NOO = 4.5 from experiment]. This means many more SHBs.
Although the gOO was not a criterion in developing the RexPoN FF, it is the first FF to capture the shape, height, and position of the first neighbor peak in the neutron-derived gOO. This excellent agreement in the shape and the height and position of the first peak in the radial distribution function for RexPoN compared with neutron scattering experiments shows that it is the reduced number of SHBs that leads to this more gradual slope in gOO at short R, compared with previous FF. It is this reduced number of SHB in water that leads to the polymeric structure.
The next question is, what is the physics driving this dynamic polymer nature of water. Along each polymer chain, we find that most waters have their SHBs align so that the OH bonds point in the same direction. This allows polarization of the charge to stabilize the polymer chain. Thus, the polymer stabilization arises from induced electrostatic attraction and polarization.
The next question is why the RexPoN FF captures this polymeric nature of water. We believe that this is because RexPoN is an FF that allows the Gaussian charges and polarization to adjust dynamically (at the 1-fs timescale), permitting the electron densities at short distances to respond to the vibrations. In addition, the long-range vdW interactions of RexPoN retain the DFT accuracy for the solid crystals that matches the experimental equation of state.
An additional question is why the QM MD calculations (PBE and SCAN) did not discover the polymer nature. This may be the timescale: We used up to 1 ns of MD to equilibrate the system with 96 or 216 waters (288 or 648 atoms), which is very far from practical with QM. The SCAN MD included 64 H2O for 20 ps of MD. It also could be because our short-range interactions are based on high-level ab initio wave functions (CCSDT) that are likely much more accurate than DFT.
Comparison with Nilsson and Saykally Experiments.
After discovering that water at 298 K has an average of ∼2 SHBs, we found a paper (Nilsson, 2004) (11) that in analyzing experimental X-ray absorption spectroscopy and X-ray Raman scattering data concluded that H2O at 298 K has 2.2 ± 0.5 SHBs, in agreement with our results! Quickly many papers appeared dismissing this interpretation. Thus, within months the Saykally group published a paper showing that the total electron yield in the near-edge X-ray absorption fine-structure spectrum is consistent with “a locally symmetric, strongly H-bonded (ice-like) domain in both solid and liquid water” (12). They concluded that only 1.5 ± 0.5 kcal/mol in hydrogen bonding is lost in going from ice at 251 K to water at 288 K, which they considered to invalidate the Nilsson interpretation. In fact, RexPoN finds a decrease of 1.7 kcal/mol over this temperature range, in agreement with Saykally. Thus, the RexPoN calculations confirm both the Saykally analysis of Nilsson’s experiments and the Nilsson (2004) interpretation of X-ray experiments.
Potential Impacts.
This revelation concerning the polymeric nature of liquid water may have a dramatic impact on our perceptions about water, with possible implications on such physical properties as viscosity, diffusion, and solvation of molecules, ions, membranes, and proteins.
One puzzling problem with water that has withstood three decades of study without a convincing explanation is the liquid–liquid critical point (LLCP) observed near 227 K. Although some consider the LLCP to occur at high pressures (13, 14), recent experimental works (using femtosecond X-ray pulses) observed the LLCP at 227 K and at ∼0 bar3.
Fig. 4 shows that the SHB lifetime in supercooled water increases monotonically as the temperature decreases until 227 K, at which point it jumps by 15% and then continuously increases monotonically as the temperature decreases further. This suggests that the supercooled critical point may correspond to the glass-transition temperature of the dynamic polymer water structure, below which the SHB lifetimes are too long to accommodate to the lower temperature, This would likely lead to a decrease from the expected density and to a maximum in the compressibility, as observed experimentally near the supercooled critical point at 227 K. To prove that RexPoN fully explains the critical-point behavior of water requires far more extensive simulations, but these results suggests the potential impact of the paradigm that water is a dynamic polydisperse branched polymer might have on physical and chemical phenomena involving water.
The predicted lifetime of the SHBs in water from RexPoN MD simulations with a periodic box of 216 H2O. Each calculation was based on a fully equilibrated structure at 298 K that was quenched to the above temperatures in 10 ps, and then equilibrated at 1 atm for 0.5 ns. The SHB lifetime increases continuously as the temperature is decreased to 227 K, then jumps by 15% (from 139 fs above 227 K to 121 fs below 227 K), and then continues to increase continuously as the temperature is decreased further. We suggest that dramatic increase in SHB lifetime at 227 K may indicate a glass transition in water polymer at 227 K that may underlie the critical point observed experimentally in supercooled water at 227 K.
Another puzzle with water is that polarization-resolved second-harmonic scattering (SHS) experiments (4) show that the orientations of water molecules are correlated to distances of 20 nm (200 Å) at 298 K whereas previous FF would have difficulty in explaining correlations longer than 1 nm. We suspect that our dynamic polydisperse polymer paradigm for water might explain such long-range correlations.
A possibly important point in explaining the critical behavior of supercooled water and the long-range angular correlations is that our simulations are on pure water, with no salt or dissolved substances. Certainly, we can expect that the presence of salt might dramatically alter the polymer nature of water. Indeed, the SHS studies mentioned above (4) show a dramatic decrease in the long-range orientational correlations of water as salt increases above some threshold.
Another indication of possible polymer character for water is that recent QM metadynamics calculations of electrocatalytic reactions with five layers of explicit solvent discovered that Grotthuss chains of up to six waters dominate the transition-state structures for the proton transfer events involved in the oxygen evolution reaction on Pt (15) and for the CO reduction reaction on Cu (16). This somewhat surprising result may arise from the polymer chains in the water solvent.
We suspect that quickly quenching the plethora of small clusters at 373 K to low temperatures might form ice much more quickly than quenching the large polymers at 298 K, which could explain the Mpemba effect (17).
Also, stabilization of these H2O polymers by antifreeze proteins might play a role in how they depress the freezing point (18).
All previous FFs for H2O were based on ΣSHB ∼ 4 at 298 K, which might dramatically affect their accuracy in predicting the properties of proteins, DNA, membranes, and materials in water solvent. It is not obvious how this paradigm that water is a dynamic polydisperse polymer on timescales of hundreds of femtoseconds might affect these results and interpretations. It will be most interesting to use RexPoN to reexamine previous MD on water-containing systems.
In summary, we developed the RexPoN FF for water based purely on accurate QM and validated that it leads to extremely accurate values for melting temperature, density, heat of vaporization, standard molar entropy, and static dielectric constant. Then we applied RexPoN to liquid water to find the very surprising result that the average water makes just two SHBs having OO distances of 2.93 or smaller. Then we analyzed the network resulting from connecting the SHBs and discovered that the nature of water is that of a dynamic multibranched polydisperse polymer network. This paradigm for water could have a significant impact on our understanding of anomalies of water and nearly any phenomena that include water.
Methods
RexPoN partitions the total energy of the system into NB and covalent energy terms. The NB part is the sum of electrostatic (Eelectro), vdW (EvdW), and HB (EHB) energies. The covalent part includes the covalent O–H bond (Ebond) and the H–O–H angle (Eangle) energy terms.
The mathematical formulation of RexPoN is described in a separate paper (9) where the training set, optimization process for each term, and complete list of parameters are provided.
Acknowledgments
We thank Dr. Sergey Zybin and Prof. Tod Pascal for helpful discussions. S.N. was supported by the Joint Center for Artificial Photosynthesis, a Department of Energy (DOE) Energy Innovation Hub, supported through the Office of Science of the US DOE under Award DE-SC0004993. W.A.G. was supported by the Computational Materials Sciences Program funded by the US DOE, Office of Science, Basic Energy Sciences, under Award DE-SC00014607. The calculations were carried out on the Extreme Science and Engineering Discovery Environment, which is supported by National Science Foundation Grant ACI-1548562.
Footnotes
- ↵1To whom correspondence should be addressed. Email: wag{at}caltech.edu.
Author contributions: S.N. and W.A.G. designed research; S.N. performed research; S.N. and W.A.G. analyzed data; and S.N. and W.A.G. wrote the paper.
Reviewers: C.L.B., University of Michigan; and M.L.K., Temple University.
The authors declare no conflict of interest.
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1817383116/-/DCSupplemental.
Published under the PNAS license.
References
- ↵
- ↵
- Ouyang JF,
- Bettens RP
- ↵
- Kim KH, et al.
- ↵
- Duboisset J,
- Brevet P-F
- ↵
- ↵
- Naserifar S,
- Brooks DJ,
- Goddard WA 3rd,
- Cvicek V
- ↵
- Oppenheim JJ,
- Naserifar S,
- Goddard WA III
- ↵
- ↵
- Naserifar S,
- Goddard WA 3rd
- ↵
- Haynes WM
- ↵
- Wernet P, et al.
- ↵
- Smith JD, et al.
- ↵
- ↵
- ↵
- Cheng T, et al.
- ↵
- Cheng T,
- Xiao H,
- Goddard WA 3rd
- ↵
- Jin J,
- Goddard WA III
- ↵
- Wen X, et al.
- ↵
- Narten AH,
- Thiessen WE,
- Blum L
- ↵
- ↵
- ↵
- ↵
- Pham TA,
- Ogitsu T,
- Lau EY,
- Schwegler E
- ↵
- ↵
- Bankura A,
- Karmakar A,
- Carnevale V,
- Chandra A,
- Klein ML
- ↵
- Wiktor J,
- Ambrosio F,
- Pasquarello A
- ↵
- Bukowski R,
- Szalewicz K,
- Groenenboom GC,
- van der Avoird A
- ↵
- Reddy SK, et al.
- ↵
- ↵
- ↵
- ↵
- Pascal TA,
- Schärf D,
- Jung Y,
- Kühne TD
- ↵
- Chen M, et al.
- Zhang C,
- Hutter J,
- Sprik M
Citation Manager Formats
Sign up for Article Alerts
Article Classifications
- Physical Sciences
- Physics
- Biological Sciences
- Biophysics and Computational Biology