Liquid water is a dynamic polydisperse branched polymer

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 (T_melt) between 273 and 273.5 K, very close to the experimental value of 273.15 K (10). This can be compared with T_melt = 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/cm³ at 298 K, agreeing exactly with experiment (it also captures the maximum density of 1.000 g/cm³ at 4 °C, where underline indicates some uncertainty). This can be compared with 0.944 g/cm³ 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 g_OO from RexPoN shown in Fig. 1. The predicted g_OO 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, N_OO = 4.7, is in good agreement with experiment, N_OO = 4.5. Fig. 1 shows that no other FF or QM description captures correctly the shape and height of this first peak in the g_OO.

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 H₂O dimer leads to a strong hydrogen bond (SHB) with a well depth of 4.98 kcal/mol at r_OO = 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.

• 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/cm³, 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 H₂O 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 H₂O and 15 side chains ranging from 1 to 22 H₂Os 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).

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 H₂Os form SHBs to just two others, with occasional branch points at waters bonded to 3 H₂Os terminating at H₂Os 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 H₂Os. The average number of H₂Os 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 H₂Os bonded to two other molecules decreasing to 35% while the number of H₂Os 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 H₂Os. 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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Structure and HB network of water. (A) Topology of the largest H₂O 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 H₂O evaporate into the vacuum. All 11 were termini (i.e., ΣSHB = 1) of short polymer chains ranging in size from 2 to 5 H₂Os. 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.