Skip to main content
  • Submit
  • About
    • Editorial Board
    • PNAS Staff
    • FAQ
    • Accessibility Statement
    • Rights and Permissions
    • Site Map
  • Contact
  • Journal Club
  • Subscribe
    • Subscription Rates
    • Subscriptions FAQ
    • Open Access
    • Recommend PNAS to Your Librarian
  • Log in
  • My Cart

Main menu

  • Home
  • Articles
    • Current
    • Special Feature Articles - Most Recent
    • Special Features
    • Colloquia
    • Collected Articles
    • PNAS Classics
    • List of Issues
  • Front Matter
  • News
    • For the Press
    • This Week In PNAS
    • PNAS in the News
  • Podcasts
  • Authors
    • Information for Authors
    • Editorial and Journal Policies
    • Submission Procedures
    • Fees and Licenses
  • Submit
  • About
    • Editorial Board
    • PNAS Staff
    • FAQ
    • Accessibility Statement
    • Rights and Permissions
    • Site Map
  • Contact
  • Journal Club
  • Subscribe
    • Subscription Rates
    • Subscriptions FAQ
    • Open Access
    • Recommend PNAS to Your Librarian

User menu

  • Log in
  • My Cart

Search

  • Advanced search
Home
Home

Advanced Search

  • Home
  • Articles
    • Current
    • Special Feature Articles - Most Recent
    • Special Features
    • Colloquia
    • Collected Articles
    • PNAS Classics
    • List of Issues
  • Front Matter
  • News
    • For the Press
    • This Week In PNAS
    • PNAS in the News
  • Podcasts
  • Authors
    • Information for Authors
    • Editorial and Journal Policies
    • Submission Procedures
    • Fees and Licenses

New Research In

Physical Sciences

Featured Portals

  • Physics
  • Chemistry
  • Sustainability Science

Articles by Topic

  • Applied Mathematics
  • Applied Physical Sciences
  • Astronomy
  • Computer Sciences
  • Earth, Atmospheric, and Planetary Sciences
  • Engineering
  • Environmental Sciences
  • Mathematics
  • Statistics

Social Sciences

Featured Portals

  • Anthropology
  • Sustainability Science

Articles by Topic

  • Economic Sciences
  • Environmental Sciences
  • Political Sciences
  • Psychological and Cognitive Sciences
  • Social Sciences

Biological Sciences

Featured Portals

  • Sustainability Science

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
Research Article

Liquid water is a dynamic polydisperse branched polymer

View ORCID ProfileSaber Naserifar and View ORCID ProfileWilliam A. Goddard III
PNAS February 5, 2019 116 (6) 1998-2003; first published January 24, 2019; https://doi.org/10.1073/pnas.1817383116
Saber Naserifar
aMaterials and Process Simulation Center, California Institute of Technology, Pasadena, CA 91125
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
  • ORCID record for Saber Naserifar
William A. Goddard III
aMaterials and Process Simulation Center, California Institute of Technology, Pasadena, CA 91125
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
  • ORCID record for William A. Goddard III
  • For correspondence: wag@caltech.edu
  1. 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:

  • Reply to Head-Gordon and Paesani: Liquid water, a branched polymer with ∼100-fs short-lived heterogeneous hydrogen bonds
    - Sep 10, 2019
  • Article
  • Figures & SI
  • Info & Metrics
  • PDF
Loading

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.

Fig. 1.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. 1.

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

View this table:
  • View inline
  • View popup
Table 1.

Summary of predicted properties from the RexPoN FF compared with experimental data, empirical FF (28, 32, 33), and QM MD [temperature Tmelt (K) at 1-atm pressure]

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.

Fig. 2.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. 2.

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

View this table:
  • View inline
  • View popup
Table 2.

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.

Fig. 3.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. 3.

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.

Fig. 4.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. 4.

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.Etotal=ENB+Ecovalent,[1]ENB=Eelectro+EvdW+EHB,[2]Ecovalent=Ebond+Eangle.[3]

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.

View Abstract

References

  1. ↵
    1. Guillot B
    (2002) A reappraisal of what we have learnt during three decades of computer simulations on water. J Mol Liq 101:219–260.
    OpenUrlCrossRef
  2. ↵
    1. Ouyang JF,
    2. Bettens RP
    (2015) Modelling water: A lifetime enigma. Chimia (Aarau) 69:104–111.
    OpenUrl
  3. ↵
    1. Kim KH, et al.
    (2017) Maxima in the thermodynamic response and correlation functions of deeply supercooled water. Science 358:1589–1593.
    OpenUrlAbstract/FREE Full Text
  4. ↵
    1. Duboisset J,
    2. Brevet P-F
    (2018) Salt-induced long-to-short range orientational transition in water. Phys Rev Lett 120:263001.
    OpenUrl
  5. ↵
    1. Bergmann U, et al.
    (2007) Nearest-neighbor oxygen distances in liquid water and ice observed by x-ray Raman based extended x-ray absorption fine structure. J Chem Phys 127:174504.
    OpenUrlCrossRefPubMed
  6. ↵
    1. Naserifar S,
    2. Brooks DJ,
    3. Goddard WA 3rd,
    4. Cvicek V
    (2017) Polarizable charge equilibration model for predicting accurate electrostatic interactions in molecules and solids. J Chem Phys 146:124117.
    OpenUrl
  7. ↵
    1. Oppenheim JJ,
    2. Naserifar S,
    3. Goddard WA III
    (2018) Extension of the polarizable -charge equilibration model to higher oxidation states with applications to Ge, As, Se, Br, Sn, Sb, Te, I, Pb, Bi, Po, and At elements. J Phys Chem A 122:639–645.
    OpenUrl
  8. ↵
    1. Shank A,
    2. Wang Y,
    3. Kaledin A,
    4. Braams BJ,
    5. Bowman JM
    (2009) Accurate ab initio and “hybrid” potential energy surfaces, intramolecular vibrational energies, and classical ir spectrum of the water dimer. J Chem Phys 130:144314.
    OpenUrlPubMed
  9. ↵
    1. Naserifar S,
    2. Goddard WA 3rd
    (2018) The quantum mechanics-based polarizable force field for water simulations. J Chem Phys 149:174502.
    OpenUrl
  10. ↵
    1. Haynes WM
    (2014) CRC Handbook of Chemistry and Physics (CRC Press, Boca Raton, FL).
  11. ↵
    1. Wernet P, et al.
    (2004) The structure of the first coordination shell in liquid water. Science 304:995–999.
    OpenUrlAbstract/FREE Full Text
  12. ↵
    1. Smith JD, et al.
    (2004) Energetics of hydrogen bond network rearrangements in liquid water. Science 306:851–853.
    OpenUrlAbstract/FREE Full Text
  13. ↵
    1. Palmer JC, et al.
    (2014) Metastable liquid-liquid transition in a molecular model of water. Nature 510:385–388.
    OpenUrlCrossRefPubMed
  14. ↵
    1. Poole PH,
    2. Sciortino F,
    3. Essmann U,
    4. Stanley HE
    (1992) Phase behaviour of metastable water. Nature 360:324–328.
    OpenUrlCrossRef
  15. ↵
    1. Cheng T, et al.
    (2017) Mechanism and kinetics of the electrocatalytic reaction responsible for the high cost of hydrogen fuel cells. Phys Chem Chem Phys 19:2666–2673.
    OpenUrl
  16. ↵
    1. Cheng T,
    2. Xiao H,
    3. Goddard WA 3rd
    (2015) Free-energy barriers and reaction mechanisms for the electrochemical reduction of CO on the Cu (100) surface, including multiple layers of explicit solvent at pH 0. J Phys Chem Lett 6:4767–4773.
    OpenUrl
  17. ↵
    1. Jin J,
    2. Goddard WA III
    (2015) Mechanisms underlying the Mpemba effect in water from molecular dynamics simulations. J Phys Chem C 119:2622–2629.
    OpenUrl
  18. ↵
    1. Wen X, et al.
    (2016) Antifreeze proteins govern the precipitation of trehalose in a freezing-avoiding insect at low temperature. Proc Natl Acad Sci USA 113:6683–6688.
    OpenUrlAbstract/FREE Full Text
  19. ↵
    1. Narten AH,
    2. Thiessen WE,
    3. Blum L
    (1982) Atom pair distribution functions of liquid water at 25 C from neutron diffraction. Science 217:1033–1034.
    OpenUrlAbstract/FREE Full Text
  20. ↵
    1. Soper AK
    (2007) Joint structure refinement of x-ray and neutron diffraction data on disordered materials: Application to liquid water. J Phys Condens Matter 19:335206.
    OpenUrlCrossRefPubMed
  21. ↵
    1. Skinner LB, et al.
    (2013) Benchmark oxygen-oxygen pair-distribution function of ambient water from x-ray diffraction measurements with a wide Q-range. J Chem Phys 138:074506.
    OpenUrlCrossRefPubMed
  22. ↵
    1. Soper A
    (2000) The radial distribution functions of water and ice from 220 to 673 K and at pressures up to 400 MPa. Chem Phys 258:121–137.
    OpenUrlCrossRef
  23. ↵
    1. Pham TA,
    2. Ogitsu T,
    3. Lau EY,
    4. Schwegler E
    (2016) Structure and dynamics of aqueous solutions from PBE-based first-principles molecular dynamics simulations. J Chem Phys 145:154501.
    OpenUrl
  24. ↵
    1. DiStasio RA Jr,
    2. Santra B,
    3. Li Z,
    4. Wu X,
    5. Car R
    (2014) The individual and collective effects of exact exchange and dispersion interactions on the ab initio structure of liquid water. J Chem Phys 141:084502.
    OpenUrlCrossRefPubMed
  25. ↵
    1. Bankura A,
    2. Karmakar A,
    3. Carnevale V,
    4. Chandra A,
    5. Klein ML
    (2014) Structure, dynamics, and spectral diffusion of water from first-principles molecular dynamics. J Phys Chem C 118:29401–29411.
    OpenUrl
  26. ↵
    1. Wiktor J,
    2. Ambrosio F,
    3. Pasquarello A
    (2017) Note: Assessment of the SCAN+rVV10 functional for the structure of liquid water. J Chem Phys 147:216101.
    OpenUrl
  27. ↵
    1. Bukowski R,
    2. Szalewicz K,
    3. Groenenboom GC,
    4. van der Avoird A
    (2007) Predictions of the properties of water from first principles. Science 315:1249–1252.
    OpenUrlAbstract/FREE Full Text
  28. ↵
    1. Reddy SK, et al.
    (2016) 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.
    OpenUrl
  29. ↵
    1. González MA,
    2. Abascal JL
    (2011) A flexible model for water based on TIP4P/2005. J Chem Phys 135:224516.
    OpenUrlPubMed
  30. ↵
    1. Berendsen H,
    2. Grigera J,
    3. Straatsma T
    (1987) The missing term in effective pair potentials. J Phys Chem 91:6269–6271.
    OpenUrlCrossRef
  31. ↵
    1. Mason PE,
    2. Brady JW
    (2007) “Tetrahedrality” and the relationship between collective structure and radial distribution functions in liquid water. J Phys Chem B 111:5669–5679.
    OpenUrlPubMed
  32. ↵
    1. Pascal TA,
    2. Schärf D,
    3. Jung Y,
    4. Kühne TD
    (2012) On the absolute thermodynamics of water from computer simulations: A comparison of first-principles molecular dynamics, reactive and empirical force fields. J Chem Phys 137:244507.
    OpenUrl
  33. ↵
    1. Vega C,
    2. Abascal JL
    (2011) Simulating water with rigid non-polarizable models: A general perspective. Phys Chem Chem Phys 13:19663–19688.
    OpenUrlCrossRefPubMed
    1. Chen M, et al.
    (2017) Ab initio theory and modeling of water. Proc Natl Acad Sci USA 114:10846–10851.
    OpenUrlAbstract/FREE Full Text
    1. Yoo S,
    2. Zeng XC,
    3. Xantheas SS
    (2009) On the phase diagram of water with density functional theory potentials: The melting temperature of ice I(h) with the Perdew-Burke-Ernzerhof and Becke-Lee-Yang-Parr functionals. J Chem Phys 130:221102–221104.
    OpenUrlCrossRefPubMed
    1. Schmidt J, et al.
    (2009) 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–11964.
    OpenUrlCrossRefPubMed
    1. Zhang C,
    2. Hutter J,
    3. Sprik M
    (2016) Computing the Kirkwood G-factor by combining constant Maxwell electric field and electric displacement simulations: Application to the dielectric constant of liquid water. J Phys Chem Lett 7:2696–2701.
    OpenUrl
    1. Grossman JC,
    2. Mitas L
    (2005) Efficient quantum Monte Carlo energies for molecular dynamics simulations. Phys Rev Lett 94:056403.
    OpenUrlPubMed
    1. Mark P,
    2. Nilsson L
    (2001) Structure and dynamics of the TIP3P, SPC, and SPC/E water models at 298 K. J Phys Chem A 105:9954–9960.
    OpenUrlCrossRef
PreviousNext
Back to top
Article Alerts
Email Article

Thank you for your interest in spreading the word on PNAS.

NOTE: We only request your email address so that the person you are recommending the page to knows that you wanted them to see it, and that it is not junk mail. We do not capture any email address.

Enter multiple addresses on separate lines or separate them with commas.
Liquid water is a dynamic polydisperse branched polymer
(Your Name) has sent you a message from PNAS
(Your Name) thought you would like to see the PNAS web site.
CAPTCHA
This question is for testing whether or not you are a human visitor and to prevent automated spam submissions.
Citation Tools
Liquid water is a dynamic polydisperse branched polymer
Saber Naserifar, William A. Goddard
Proceedings of the National Academy of Sciences Feb 2019, 116 (6) 1998-2003; DOI: 10.1073/pnas.1817383116

Citation Manager Formats

  • BibTeX
  • Bookends
  • EasyBib
  • EndNote (tagged)
  • EndNote 8 (xml)
  • Medlars
  • Mendeley
  • Papers
  • RefWorks Tagged
  • Ref Manager
  • RIS
  • Zotero
Request Permissions
Share
Liquid water is a dynamic polydisperse branched polymer
Saber Naserifar, William A. Goddard
Proceedings of the National Academy of Sciences Feb 2019, 116 (6) 1998-2003; DOI: 10.1073/pnas.1817383116
Digg logo Reddit logo Twitter logo Facebook logo Google logo Mendeley logo
  • Tweet Widget
  • Facebook Like
  • Mendeley logo Mendeley
Proceedings of the National Academy of Sciences: 116 (6)
Table of Contents

Submit

Sign up for Article Alerts

Article Classifications

  • Physical Sciences
  • Physics
  • Biological Sciences
  • Biophysics and Computational Biology

Jump to section

  • Article
    • Abstract
    • Results
    • Discussion
    • Methods
    • Acknowledgments
    • Footnotes
    • References
  • Figures & SI
  • Info & Metrics
  • PDF

You May Also be Interested in

Abstract depiction of a guitar and musical note
Science & Culture: At the nexus of music and medicine, some see disease treatments
Although the evidence is still limited, a growing body of research suggests music may have beneficial effects for diseases such as Parkinson’s.
Image credit: Shutterstock/agsandrew.
Large piece of gold
News Feature: Tracing gold's cosmic origins
Astronomers thought they’d finally figured out where gold and other heavy elements in the universe came from. In light of recent results, they’re not so sure.
Image credit: Science Source/Tom McHugh.
Dancers in red dresses
Journal Club: Friends appear to share patterns of brain activity
Researchers are still trying to understand what causes this strong correlation between neural and social networks.
Image credit: Shutterstock/Yeongsik Im.
Yellow emoticons
Learning the language of facial expressions
Aleix Martinez explains why facial expressions often are not accurate indicators of emotion.
Listen
Past PodcastsSubscribe
Goats standing in a pin
Transplantation of sperm-producing stem cells
CRISPR-Cas9 gene editing can improve the effectiveness of spermatogonial stem cell transplantation in mice and livestock, a study finds.
Image credit: Jon M. Oatley.

Similar Articles

Site Logo
Powered by HighWire
  • Submit Manuscript
  • Twitter
  • Facebook
  • RSS Feeds
  • Email Alerts

Articles

  • Current Issue
  • Special Feature Articles – Most Recent
  • List of Issues

PNAS Portals

  • Anthropology
  • Chemistry
  • Classics
  • Front Matter
  • Physics
  • Sustainability Science
  • Teaching Resources

Information

  • Authors
  • Editorial Board
  • Reviewers
  • Librarians
  • Press
  • Site Map
  • PNAS Updates

Feedback    Privacy/Legal

Copyright © 2021 National Academy of Sciences. Online ISSN 1091-6490