Water’s dual nature and its continuously changing hydrogen bonds

A model is proposed for liquid water that is a continuum between the ordered state with predominantly tetrahedral coordination, linear hydrogen bonds and activated dynamics and a disordered state with a continuous distribution of multiple coordinations, multiple types of hydrogen bond, and diffusive dynamics, similar to that of normal liquids. Central to water’s heterogeneous structure is the ability of hydrogen to donate to either one acceptor in a conventional linear hydrogen bond or to multiple acceptors as a furcated hydrogen. Linear hydrogen bonds are marked by slow, activated kinetics for hydrogen-bond switching to more crowded acceptors and sharp first peaks in the hydrogen-oxygen radial distribution function. Furcated hydrogens, equivalent to free, broken, dangling or distorted hydrogens, have barrierless, rapid kinetics and poorly defined first peaks in their hydrogen-oxygen radial distribution function. They involve the weakest donor in a local excess of donors, such that barrierless whole-molecule vibration rapidly swaps them between the linear and furcated forms. Despite the low number of furcated hydrogens and their transient existence, they are readily created in a single hydrogen-bond switch and free up the dynamics of numerous surrounding molecules, bringing about the disordered state. Hydrogens in the ordered state switch with activated dynamics to make the non-tetrahedral coordinations of the disordered state, which can also combine to make the ordered state. Consequently, the ordered and disordered states are both connected by diffusive dynamics and differentiated by activated dynamics, bringing about water’s continuous heterogeneity.

Evidently, water displays characteristics of both homogeneous and two-state models and is therefore intermediate between these two extremes and seemingly contradictory. The purpose of this study is to resolve this dichotomy and find a model consistent with both of them. This will be done by examining how different coordinations affect the free energy barriers and rates of HB switching between different coordinations. It builds on an earlier study which investigated the nature of non-tetrahedral water [88]. Using a trajectory from an all-atom molecular dynamics simulation of ambient water, we examine the potentials of mean force, radial distribution functions, and switching kinetics for every hydrogen as a function of three factors: the number of acceptor HBs of the hydrogen's nearest and next-nearest acceptors and the hydrogen's rank to the nearest acceptor. One important problem to address is how to define HBs. Many factors are known to contribute to the strength of HBs such as whole molecule vibration [88,89], bifurcation of a donor to more than one acceptor [17-21, 44, 88, 90-92], strain induced by small, closed HB networks such as cyclic dimers [25,89] and cyclic trimers [88,93,94], competition between donors to a single acceptor [25,88,95], secondshell effects [83,96,97] and polarization [98][99][100][101]. It has been proposed that it may not even be possible to determine water's HBs instantaneously [4,65,72,94,102]. Nonetheless, almost all HB definitions use a fixed, mean-field, two-body deviation from the idealized linear geometry that is independent of other neighbouring molecules. Hydrogens with no acceptor within the cut-off are said to have a broken HB even though acceptors are usually not far beyond. The nearest-acceptor HB definition [25], on the other hand, whereby a donor donates to its nearest acceptor, considers each hydrogen's surrounding environment on a case-by-case basis according to the relative sizes of OH distances and involves no fixed cut-off parameters. It relies on the property that a hydrogen usually has one much nearer non-covalently bound oxygen and assumes that the transition state of HB switching has the hydrogen midway between the nearest two acceptors. Provision can be made for having fewer [25] or more accepting oxygens [88]. While the overall results here should not depend too much on the HB definition, we use the nearest-acceptor HB definition because the predicted HBs closely match those in quenched inherent structures [25]. A second problem in the analysis of HB kinetics is the need to extract rate constants for each type of HB switch. Rather than deriving them from the logarithm of the decay of the survival functions of each HB, a much simpler equilibrium kinetics analysis is implemented here based on the rates of HB switching for each type of switch. These rates are constant because the fractions of each HB coordination are constant at equilibrium, with the constant of proportionality in each case being the rate constant.

Molecular dynamics simulations
A cubic periodic box of 375 TIP4P/2005 water molecules [103] and volume approximately 2.2 2.2 2.2 × × nm 3 is simulated using the sander module of the AMBER 14 software package [104]. 500 steps of steepest-descent minimization are followed by a molecular dynamics simulation in the NVT ensemble at 298 K using a Langevin thermostat with collision frequency of 5 ps −1 for 100 ps, and in the NPT ensemble using a Berendsen barostat at 1 bar and relaxation time 2 ps for a further 2 ns. Data collection proceeds over a further 1 ns under the same NPT conditions with SHAKE applied to all bonds, Figure 1. Schematic of the distances and angles for the nearest and next-nearest acceptors relative to the hydrogen (H) of interest. r HO,1 is the distance of a hydrogen to its nearest non-bonded oxygen (O) and r OO,1 is the distance from the hydrogen's bonded oxygen to its nearest non-bonded oxygen. θ HOO,1 is the angle between the hydrogen's bonded HO vector and the OO vector between a hydrogen's bonded oxygen and its nearest non-bonded acceptor. The quantities denoted by 2 instead of 1 are for the hydrogen's next-nearest oxygen. a 2 fs timestep, periodic boundary conditions, an 8 Å cutoff and particle mesh Ewald with default AMBER parameters. Coordinates are saved every 2 fs for analysis, giving 500 000 structures in total. The calculated density under these conditions is 0.997 0.008 ± g cm −3 .

Hydrogen bonds (HBs) analysis
HBs are defined using the nearest-acceptor definition [25] whereby a hydrogen donates to its nearest non-covalently bonded oxgyen. To keep the number of coordination states to a manageable number, broken HBs are not considered in the analysis, whereby only the strongest donor-acceptor interaction between a pair of water molecules would be allowed as an HB [25]. The number of HBs a molecule accepts is denoted a and the fraction of molecules accepting a HBs is denoted x(a). To every hydrogen is assigned the acceptor number of its nearest acceptor a 1 and next-nearest acceptor a 2 . Fractions x(a 1 ) and x(a 2 ) are recorded for each molecule with nearest and next-nearest acceptors a 1 and a 2 , respectively. Each nearest acceptor and next-nearest acceptor pair of a hydrogen is denoted (a a , 1 2 ) and its fraction x a a , 1 2 ( ). Snapshots are analysed every 2 fs. Errors are determined from the standard deviation of averages from five 200 ps blocks.

Potentials of mean force (PMFs)
PMFs are used to plot the free energy as a function of reaction coordinate for the switch of a donor from its nearest acceptor to its next-nearest acceptor. Similar to previous analysis [88], the reaction coordinate is r r r HO HO,1 HO,2 ∆ = − , where r HO,1 is the HO distance of a hydrogen to its nearest oxygen and r HO,2 is the distance to its next-nearest oxygen, both excluding the covalently bonded oxygen (see figure 1). For TIP4P/2005 specifically, the distances are taken between the hydrogen and the site of the oxygen's charge. A histogram N r HO ( ) ∆ for the count of r HO ∆ with a 0.01 nm bin-width is constructed and the PMF of the switch is calculated using the equation where R is the universal gas constant and T is the temperature. Normalization of N r HO ( ) ∆ , which shifts the PMF by a constant, is unnecessary when only differences are examined.
The equation for a hydrogen's switch from the nearestacceptor with a 1 to the next-nearest acceptor with a 2 is written as a a a a , 1 , 1 where the order of molecules is the same for both the reactant and product states, even though the second molecule in the pair is the nearest acceptor in the product state. The trans ition state of a HB switch, where the donor switches from its nearest-acceptor to its next-nearest acceptor, is assumed to be where r 0 HO ∆ = , an assumption that will be examined later. Compared to the previous study which only showed the reactant half of the PMF up to the transition state [88], this study plots the full PMFs from reactants to products. When averaged over all switches, the product PMF with r 0 HO ∆ > is the mirror image of the reactant PMF with ∆ < r 0 HO by symmetry. PMFs are also constructed as a function of a 1 or (a a , 1 2 ). The reactant PMF, defined as ∆ < r 0 HO , is constructed directly from r HO ∆ . The product PMF with r 0 HO ∆ > for a reaction starting from (a a , 1 2 ) is constructed from the reactant PMF of the switch a a a a 1, 1 , , the reverse of equation (2). Further PMFs are constructed as a function of the ranking of the donor to its nearest acceptor a 1 1 … together with (a a , 1 2 ). To examine an angular reaction coordination instead of the one in distance, PMFs are constructed using a 5° bin-width as a function of HOO HOO,1 HOO,2 θ θ θ ∆ = − , where HOO,1 θ is the angle between the hydrogen's bonded HO vector and the OO vector connecting a hydrogen's bonded oxygen and nearest acceptor (see figure 1). HOO,2 θ is defined similarly except that the OO vector involves the hydrogen's nextnearest acceptor. To show the PMFs in terms of the explicit values of r OH,1 and r OH,2 rather than only their difference, two-dimensional PMFs are built for r HO,1 and r HO,2 using a 0.01 nm bin-width. This is done for all hydrogens or for only the outermost hydrogens. Snapshots are analysed every 2 fs   for all PMFs except the 2D PMFs for which the snapshot frequency is every 100 fs.

Radial distribution functions
g r HO ( ) and g r OO ( ) are constructed using a 0.01 nm bin-width where r HO is the HO distance as before and r OO is the OO distance, taken from the Lennard-Jones centres of the TIP4P/2005 oxygen. g r HO ( ) is plotted according to the value of the contributing hydrogen's value of a 1 , (a 1 , a 2 ) or (a a , 1 2 ) including rank. g r OO ( ) is plotted according to the values of a 1 , (a 1 , a 2 ) or (a a , 1 2 ) including rank for each hydrogen bonded to the oxygen. Snapshots are analysed every 10 fs.

HB switch percentages and rate constants
A HB switch is defined at the point when r HO ∆ changes from negative to positive. It is referred to here as a donor switch to distinguish it from the two other kinds of switches considered below, which are termed nearest-acceptor and nextnearest-acceptor switches. The HB switch percentages for a given acceptor pair are calculated by dividing the number of switches of (a a , 1 2 ) by the total number of switches. Rate constants k a a , 1 2 ( ) for a HB switch (a a , 1 2 ) are calculated using the rate equation of first-order kinetics x a a t k a a x a a d , d , , At equilibrium, x a a , 1 2 ( ) is constant, allowing equation (3) to be solved for k a a , 1 2 ( ). Rate constants in units of ps −1 are calculated from the total number of switches divided by the number of molecules and by the length of the simulation in ps. Switch percentages and rate constants are also accrued for a given a 1 , a 2 or for all switches.
Switch percentages and rate constants are calculated for a second kind of switch, termed a nearest-acceptor switch, for how often a hydrogen's a 1 changes because another donor switches to or from that acceptor rather than the hydrogen switching itself. The possible values of a 1 ∆ considered are ±1 because instances where a 1 changed by two after a single timestep were found to be negligible. These rate constants are evaluated as a function of a 1 but not of a 2 because they were found to depend weakly on a 2 . Switch percentages and rate constants are also calculated for how often a 2 changes without the hydrogen of interest switching, termed a next-nearest-acceptor switch. This happens either because another donor to the next-nearest acceptor changes or because the next-nearest acceptor itself is replaced by another molecule which was previously more distant. Because the second mechanism involves a change of molecule, a 2 ∆ may now be ±1, ±2 or ±3. If a next-nearest neighbour is replaced by another molecule with the same a such that a 0 ∆ = , no change is recorded. These rate constants are evaluated as a function of a 2 but not of a 1 because little dependence was found on a 1 . Overall, donor switches, nearestacceptor switches and next-nearest-acceptor switches represent the three possible fates of a hydrogen's acceptor pair. For all rate constants, snapshots are analysed every 2 fs. Errors are determined from the standard deviation of averages from five 200 ps blocks.

Hydrogen bonds (HBs)
The percentages of x(a) are listed in table 1 and those of x(a 1 ), x(a 2 ) and x a a , 1 2 ( ) are listed in table 2. Their standard deviations are mostly negligible. By definition, there are no hydrogens with a 1 = 0 because every hydrogen has a nearest oxygen and so this column is omitted. The proportions of x(a) are similar to those reported earlier [25] for TIP4P/2005 but without the small difference due to broken HBs. The values of x(a 1 ) differ from x(a) because x(a 1 ) is weighted by the number of hydrogens donating to acceptor a 1 . Interestingly, x(a 2 ) also differs from x(a 1 ) in that molecules with lower a 2 are more likely. This reflects local anti-correlation between acceptors of high and low values because every donor is matched by an acceptor. This anti-correlation also manifests in the non-ideal crossterms x a a , 1 2 ( ) for acceptors of high and low values which otherwise differ little from the uncorrelated values x a x a

Potentials of mean force (PMFs)
The PMF averaged over all switches as a function of r HO ∆ is shown in figure 2. Its free energy barrier is 4 kJ mol −1 . When the PMFs are plotted as a function a 1 (figure 3), they now differ. The higher the value of a 1 , the lower the barrier, and the more favourable the products versus the reactants. Evidently, the more donors competing for the same acceptor, the weaker their HBs and the more readily they move to another acceptor. For a going from 1 to 4, the respective barriers are 9, 6, 2 and <1 kJ mol −1 while the free energy differences of product minus reactant range from 9, 4, −4, and −8 kJ mol −1 .
When the PMFs are plotted as a function of acceptor pairs (a a , 1 2 ) (figure 4), they differ even more. The larger a 1 and the smaller a 2 , the lower is the barrier and the more stable the products than the reactants. The PMFs on the diagonal with a a 1 1 2 = + are symmetrical with the hydrogen having an equal preference for its nearest two acceptors and only a small barrier between them. The (2, 2) hydrogen is stable and has a 7 kJ mol −1 barrier to climb to reach the barely stable (1, 3) pair. Hydrogens with a a 1 2 < are even stabler than the tetrahedral (2, 2) case while hydrogens with a a 1 1 2 > + have no barrier whatsoever. By looking up the acceptor fractions in table 4, it can be seen that 67% of hydrogens are stably bonded to their nearest acceptor but 33% are not: 27% have an equal preference for their nearest two acceptors while 6% are metastably bonded and would prefer to donate to a 2 . The angular PMFs as a function of When the PMFs in r HO ∆ are plotted according to the rank of a hydrogen donating to its nearest acceptor as well as (a a , 1 2 ) (figure 6), the PMFs differ even further. Here we consider the outermost hydrogen, shown in black, which is the one most likely to switch. Of the four most prevalent cases (table 2), the outermost hydrogen of (2, 2) has a large positive free energy difference to shift from one acceptor to the other while that of (3, 1) has no barrier and a large negative free energy difference. For (3,2) there is no barrier at all and for (2, 1) the barrier is less than 2 kJ mol −1 . The outermost hydrogens for cases with a greater excess of a 1 over a 2 are even less stable.
To see more clearly how the switch depends on the individual values of r HO,1 and r HO,2 , 2D PMFs are plotted for different (a a , 1 2 ) for all hydrogens in figure 7 and for just the outermost hydrogens in figure 8. These plots mirror those in figures 4 and 6 but reveal that HB switching takes place via two distinct mechanisms, either directly via a single transition state for shorter HBs, or via two transition states and the broad plateau of a weakly stable intermediate with longer HBs. The path via the single transition state has r 0.23 HO ∼ nm, indicative of a bifurcated hydrogen, while the path via the intermediate has r 0.3 HO ∼ nm. Given that the average OH distance to a hydrogen's third and fourth-nearest acceptors are calculated from the simulation to be 0.31 and 0.33 nm, which is also close to the second peak in g r HO ( ), this intermediate represents what could be interpreted as either a free hydrogen or a hydrogen furcated to more than two acceptors. The similar occurrence of this hydrogen for all plots is consistent with the nearest two acceptors not being so influential.
PMFs were also considered as a function of the acceptor values for the nearest-three molecules (a a a , , 1 2 3 ) where a 3 is the number of HBs accepted by the switching hydrogen's  third-nearest acceptor. However, they showed little dependence on a 3 , with plots very similar to those in figure 4 (data not shown). Similarly, PMFs as a function of (a a a , , 1 2 ) had little dependence on a (data not shown). The acceptor numbers of other water molecules, being more distant, are unlikely to have a significant effect on a hydrogen's PMF of HB switching.

Radial distribution functions
The HO and OO radial distribution functions g r HO ( ) and g r OO ( ) as a function of a 1 are shown in figure 9. The larger a 1 , the less defined is the first peak, which matches the lower barrier of HB switching. The differences are stronger for the plots versus r HO , first because the first shell depends almost exclusively on a 1 , whereas the first shell of the r OO plots includes OO distances of molecules with a values other than a 1 , and second because the rankings are based on r HO rather than r OO . An important distinction should be made with the g(r) plots made in earlier work [25] in that those plots are a function of a, the number of hydrogen bonds accepted by the water in question. They show the opposite trend that larger a leads to the larger first peak in g(r), reflecting the larger number of hydrogens donating to that acceptor. When the g r HO ( ) plots are shown as a function of a 2 and the ranking to the nearest acceptor, as in figure 10, the delocalization of the hydrogen is even more striking. For a 1 equal to 3 or 4, the first peak is essentially absent for the outermost hydrogen and there is little dependence on a 2 . The peaks for the other hydrogens are sharper but become less defined for lower values of a 2 . The g r OO ( ) plots in figure 11 show little dependence on a 2 for the reasons already discussed.

Switch percentages
The percentages of donor, nearest-acceptor and next-nearestacceptor switches are listed in table 3 according to (a a , 1 2 ), a 1 and a 2 . The total number of switches recorded in 1 ns is 1 519 346. Four acceptor pairs dominate the donor switches: (3,2) at 29%, (2, 2) and (3, 1) at 21%, and (2, 1) at 16%. This trend is consistent with an earlier study [25], which obtained slightly different values because it considered multiple, simultaneous switches as distinct compound switches. The number of nearest-acceptor switches in the simulation is 4 865 617 and the number of next-nearest-acceptor switches is 15 258 001. Evidently, the number of next-nearest-acceptor switches is 3.2 times the number of nearest-acceptor switches, which in turn is 3.2 times the number of donor switches. Because the number of hydrogens grows with distance, these longerrange switches occur more often, especially the next-nearestacceptor switches which may only involve whole-molecule vibration exchanging second and third-nearest acceptors rather than a HB switch. The most common nearest-acceptor switch  at 33% alternates between acceptors with a equal to 2 and 3. Almost all next-nearest-acceptor switches have ∆ = ± a 1 2 (92%) while the most common value of a 2 is 2 (46%), reflecting the predominance of the tetrahedral coordination. Unlike for the nearest-acceptor switches, the contrib ution of a 2 = 3 is comparable to that of a 2 = 1. The prevalence of a switch depends on both the fraction of the acceptor pair and its stability. The contribution of fraction can be factored out by considering the rate constants, which are examined next. Table 4 shows the rate constants of HB switching determined by equation (3) as a function of (a a , 1 2 ), as well as running totals as a function of a 1 , a 2 and the grand total. The trends of the rate constants are consistent with the PMFs in that larger values of a 1 and lower values of a 2 increase the rate constant. This matches the trends in the sizes of the barriers in the PMF. Also shown in table 4 are the rate constants of the nearestacceptor switches listed according to a 1 and a 1 1 ∆ = ± , and the rate constants of the next-nearest-acceptor switches according to a 2 and a 2 ∆ . The greater frequency of the longer-range switches is mirrored in the larger overall rate constants. The rate constant of next-nearest-acceptor switches, 20.3 ps −1 , is 3.2 times that of the nearest-acceptor switches, 6.4 ps −1 , which in turn is 3.2 times that of the donor switches, 2.0 ps −1 . For comparison and not given in table 4, the rate constant for a change of rank ±1 within the set of donors for a given acceptor is ∼12 ps −1 . This rate constant reflects whole-molecule vibration and is found to be largely independent of the acceptor pair. As expected, the rate constants for a 1 1 ∆ = + become smaller with increasing a 1 because it is harder to switch to a more crowded acceptor. The opposite is seen for a 1 1 ∆ = − . The combined effect is that molecules with a 1 deviating from two are more likely to switch, consistent with the stability of the tetrahedral coordination. For the next-nearest-acceptor switches, rate constants are larger for smaller a 2 ∆ and the further a 2 deviates from 2.

Discussion
The results presented here reveal a strong dependence of a hydrogen's stability and switching mechanism on the surrounding water molecules. The PMFs may be either favourable or unfavourable, and with negligible or substantial barriers relative to the thermal energy. The most important and representative PMFs are the four PMFs starting from the (2, 1), (2, 2), (3, 1) and (3, 2) acceptor pairs. Of these, there are effectively only two kinds of PMF because (2, 2) and (3,1) are the mirror-image of each other while (2, 1) and (3, 2) both have similar shapes. Hence the discussion is focussed on these acceptor pairs. The PMFs for other acceptor pairs occur less commonly and are generally more extreme forms of the same trends.
A key feature associated with the non-tetrahedral coordinations is that a hydrogen's linear HB is destabilized when the hydrogen's nearest acceptor has more HBs than its nextnearest acceptor and when the hydrogen is the outermost donor. If the difference a a 1 2 − is only one, then the outermost hydrogen has an equal preference for its nearest and nextnearest acceptors. For the (3, 2) pair the hydrogen is either bifurcated between the two nearest acceptors or furcated between even more acceptors. The (2, 1) pair is similar but marginally more stable in the linear form. The g r HO ( ) plots are consistent with this and indicate that the first peak for the outermost hydrogen is absent for (3, 2) and substantially weakened for (2, 1). There is some ambiguity about how one might describe these hydrogens. The bifurcated form has been proposed before as important in water [17-21, 44, 88, 90-92]. As noted previously [88], their 11% prevalence closely matches the 10% of hydrogens defined as broken HBs of many models [27,35,49,51,54,55] or cut-off-based HB definitions [14,105,106]. These hydrogens could also be interpreted as free, dangling, distorted, multifurcated, non-linear or delocalized. Arguments in favour of terms which suggest them having multiple weak interactions are the negligible amount of nonred-shifted gas-phase OH stretches in the liquid [2,72,107,108] the ability of these hydrogens to interact with both the first and second shells, collapsing the second shell and bringing about the higher density and larger coordination of less tetrahedral regions [31,37,40,46,47,58,59]. Recognising these subtle differences, for the remainder of this text we use the term 'furcated' hydrogen to reflect their interactions, albeit weak, with multiple acceptors.
As to whether HB switching via the furcated hydrogens of the (3, 2) and (2, 1) pairs is activated or barrierless, figure 4 indicates that a hydrogen moving from one acceptor to another experiences a small barrier. However, figure 6 shows that at any one time the outermost hydrogen for a (3,2) pair experiences no barrier. The ranking of hydrogens continually interchanges due to whole molecule vibration in a barrierless process involving thermal energies. This means that the barrier for the whole process is less than the thermal energy and is essentially equivalent to there being no barrier. A similar picture holds for the outermost hydrogen for a (2, 1) pair, except that there is still a small barrier, making these hydrogens switch about half as often as outermost (3,2) hydrogens (tables 3 and 4). We distinguish these two     on the switching of neighbouring molecules, both of which are more rapid on average than any single switch (table 4). The other two main kinds of switch are between the (2, 2) and (3, 1) pairs. The switch from the (2, 2) pair is activated and slower than the furcated switches with a rate constant of 0.9 ps −1 . The slowness is because the hydrogen is moving to a more crowded acceptor and is consistent with the slower dynamics of tetrahedral coordinations [9,20,26,33,38,109,110]. Despite this, it comprises 21% of the switches because of the high fraction of (2, 2) pairs. It indicates that hydrogens do not always switch from over-coordinated to under-coordinated molecules. The resulting (3, 1) pair represents an unstable hydrogen. Such hydrogens occur when the difference a a 1 2 − exceeds one, such that the outermost hydrogen's interaction with its nearest acceptor is unstable relative to its next-nearest acceptor and it would prefer to switch to its next-nearest acceptor. The (3, 1) pair makes up 5% of hydrogens, with the other such pairs making up 1%. The hydrogen with a (3, 1) acceptor pair has three main fates, all of them largely barrierless. It may switch to its next-nearest acceptor to make the more stable (2, 2) pair with a rate constant of 8.2 ps −1 (table 4) in a barrierless switch. With almost equal probability, a nearest-acceptor switch may lower the hydrogen's rank and change the pair to (2, 1) with a rate constant of 10.1 ps −1 . Alternatively, the next-nearest-acceptor switch may change the pair to (3,2) with a rate constant roughly three times larger at 25.8 ps −1 . These last two mechanisms are the route by which furcated hydrogens are created from molecules with tetrahedral coordination. For the over-and under-coordinated furcated hydrogens to find each other to form an unstable hydrogen, there is a mild entropic barrier of 3.4 kJ mol −1 , estimated from table 2 as − = − ( ( )/ ( )) ( / ) RT x x RT ln 3, 1 3 ln 5.24 21.12 . Based on all these findings, the different mechanisms of HB switching with their barrierless or activated dynamics are consistent with the multitude of data that water has heterogeneous structure and dynamics. One state, comprising about a third of the hydrogens under ambient conditions, is characterized by a diversity of coordinations, both tetrahedral and  non-tetrahedral, and a diversity of hydrogens, whether furcated, unstable or as linear HBs, which continuously interchange with rapid HB switching and diffusive dynamics. The remaining two-thirds of hydrogens have more linear HBs, predominantly tetrahedral coordination and slower, activated dynamics. These states could be mapped onto the traditional two-state models for water with a stabler, lower-density tetrahedral state and a dynamic, disordered, high-density state, the latter more characteristic of a normal liquid. However, the classification used here is artificially discrete, HB-definitiondependent, and averages over the more continuous shifts between different HBs and PMFs as barriers grow and fall, as molecules vibrate and exchange their nearest and next-nearest neighbours and for individual molecules whose two hydrogens straddle more and less ordered regions. Importantly, the ordered and dis ordered states interconvert by two mechanisms: either non-tetrahedral coordinations encroach on tetrahedral coordinations via barrierless HB switches, or tetrahedral coordinations convert to non-tetrahedral coordinations in activated HB switches or do the reverse in entropically activated HB switches. Via the diffusive mechanism, water's structure continuously shifts between the two states, such that neither state is pure nor separated by sharp boundaries. Via the activated mechanism, the two states are separated not by a boundary between different types of water molecule in space but by a HB switch or by entropic dilution. Hence water is best described as a continuum of ordered and disordered states. This is intermediate between the two main kinds of models proposed for water, namely mixture and continuum models. At one extreme, two-state models that refer to water as a 'mixture' imply that the states are distinct entities separated by sharp boundaries in space. However, the lack of evidence for such boundaries has been a continual point of criticism of such models [1][2][3][4][5][6][7], even though water's tetrahedral component is known to partially self-aggregate [39,40,55,59,97,109,[111][112][113][114][115][116]. At the other extreme, treating water as a single homogeneous continuum, which combines the ordered and disordered states, convolutes the homogeneity of  the dominant tetrahedral coordination with the homogeneity of diffusive dynamics characteristic of a normal liquid. These homogeneities are quite different, with the crystal and gas being the respective limiting forms. Their duality lies at the centre of a related long-standing question of whether a liquid is more like a crystal, more like a gas or shares a continuum of both [28,55,63,[117][118][119][120]. The homogeneous crystal-like or gas-like models both fail at the opposite extremes, and only heterogeneous crystal-gas continuum models succeed over the full liquid-range. We examine in more detail the features which make the heterogeneous behaviour prevalent in water compared to normal liquids. In normal liquids the weaker and less-directional interactions bring about close-packing and higher coordination numbers than in water. The average coordination is lower than the stable crystalline coordination, which is relatively rare in the liquid [121][122][123]. This brings about a continuum of smaller molecular motions, a smoother energy landscape and diffusive dynamics. Local dynamics increases only gradually with decreasing coordination state accompanied by decreasing density. There is little of the activated dynamics and heterogeneous nature of water except near the glass transition [124][125][126]. By way of contrast, water's most notable feature is its well-known ability to form strong, linear HBs. Their highly directional nature, water having only two hydrogens, and the hydrogens subtending a near-tetrahedral bond angle, make water's most stable packing the unusually low tetrahedral coordination with its notably low density. This means that molecules in the liquid can access coordination states that are either higher, lower or the same as the stable tetrahedral coordination, with little scope for fluctuations more than ±1. Therefore, a single coordination, the stable tetrahedral coordination, dominates the liquid phase. Directional HBs bring about strong correlations between the orientations of neighbouring water molecule. These correlations extend over many neighbours [109,127,128], making the dynamics of water molecules highly dependent on not just their first shell but also the second and third shells. Both the low tetrahedral coordination and the strong correlations are conducive to activated, jump-like dynamics for a number of reasons. First, strong, directional HBs are obviously stabilising. Second, a donor has few nearby acceptors to which to switch. Third, it is difficult for a donor to switch to many of these acceptors either because they may already be engaged in HBs to the donor's molecule, because they are in an unfavourable orientation, because they are relatively distant owing to the open structure, or because the next-nearest acceptor is equally or more crowded with other donors than the donor's nearest acceptor. Water still accesses the disordered, diffusive state by way of the furcated hydrogens. Their number may be low but this is compensated by their ability to be readily created from a single HB switch from the ordered state, predominantly the activated (2,2) to (1,3) switch. They can rapidly Table 3. Switch percentages a for donors and nearest and next-nearest acceptors.  exchange in a barrierless fashion coupled to whole-molecule vibration, reducing the stability of multiple surrounding molecules, permitting diffusive dynamics and making water a liquid. For water at higher temperatures [129] or pressures [34], the disordered state dominates as in a normal liquid, and the tetrahedral model clearly fails. With decreasing temper ature, at low pressure the ordered state dominates and converges to the low-density liquid, amorphous ice and ice Ih phases, while at higher pressure the disordered state dominates and converges to the high-density liquid, amorphous and crystalline ice phases [29,43,84,86,130]; as for any liquid near its glass transition, dynamics is highly heterogeneous. Whether water's two states ever truly separate at a liquid-liquid transition with an associated critical point remains a tantalizing conjecture, given the insufficient stability of supercooled water to crystallization [85,87,[130][131][132]. There is no requirement for them to do so in the proposed model. Regarding more technical aspects of this study, the results make clear that the stability of each hydrogen is highly cooperative and depends on the HBs of its neighbours and cannot be determined from a donor-acceptor interaction in isolation as quantified by a typical HB definition or by the nearest-acceptor definition. Nor can it be determined from the average PMF, which hides heterogeneity and by itself would suggest that water's dynamics is moderately activated. Nor can it be assessed by looking only at the coordination of individual molecules in isolation. In the proposed model, tetrahedral and non-tetrahedral coordinations appear in both the ordered and disordered states, and correlations out to the third shell need to be considered in assessing a molecule's stability, if not further. Such analyses can be quite challenging. A brute-force second-shell approach involving 16 neighbours is difficult enough [83,97], assuming molecules have four neighbours. A brute-force thirdshell analysis would have to examine on the order of a daunting 4 4 4 64 × × = neighbouring molecules. The third-shell approach here is more efficient by considering only a subset of 4 + 4 = 8 neighbours on average per hydrogen, giving 16 neighbours in total per molecule, based on the assumption that the main determinant of stability is the number of acceptors.
An issue for any kinetics study is the choice of reaction coordinate. PMFs inevitably average over all other degrees of freedom and may miss important mechanisms. This is made very clear from the qualitatively different PMFs as a function of a 1 , a 2 , and the donor's ranking to the nearest acceptor. How meaningful the sizes of the free energy barriers are depends on how closely the transition states align for the individual PMFs being averaged over. The nearest-acceptor HB definition [25] performs reasonably well at identifying the transition states of HB switching for most types of switch, based on how closely the maximum in the PMFs match the point where the r 0 OH ∆ = . This also validates the assignment of HBs. However, it misses the absence of a transition state for   the outermost furcated or unstable hydrogens. Including the broken HB criterion based on competiting HO interactions [25] should help address this. The overall rate constant of HB switching at 2.0 ps −1 gives a HB decay time of 0.5 ps which is on the low side compared to 1 ps values by other measures [77,133]. This discrepancy is well-recognized and occurs because transient fluctuations away from a linear HB are not ignored [13,134,135]. Furcated and unstable hydrogens contribute disproportionately to HB breakages because of their frequent recrossing of the point r 0 HO ∆ = . Such breakages could be excluded as HB switches [14,33,70,71,99,105,[134][135][136] and would certainly lead to longer decay times and smaller rate constants. However, this would favour the activated switches between stable start and end states at the expense of the diffusive switches. The kinetics analysis used here, based on constant rates at equilibrium rather than non-equilibrium decay, is straightforward and convenient to implement but does not permit an examination of any deviations from first-order exponential decay, deviations that are more likely for diffusive dynamics. Nonetheless, these issues do not affect the proposed structure because the PMFs, on which the proposed structure is based, are calculated independently of the rate constants.
The final issue concerns the use of a classical all-atom force field. It remains to be seen if the inclusion of other effects may alter the picture such as polarization [49,[98][99][100][101], as captured in polarisable force fields or ab initio methods. The balance between the two kinds of structure relates to the level of ordering indicated in standard measures such as the radial distribution function. More disorder means fewer linear HBs and more furcation. Given that the TIP4P/2005 water model is known to generate a more-ordered structure than real water [103], the proportion of the disordered state is likely to be even larger. Nuclear quantum effects [137,138], which delocalize protons between their surrounding acceptors, may further shift the balance in non-obvious ways.

Conclusions
This study presents a coherent model for the structure and dynamics of liquid water that combines both the competing continuum and mixture models. Water is made up of a continuum of two states that interchange via either diffusive or activated dynamics. Consistent with two-state models, the ordered state has the dominant tetrahedral coordination and is characterized by linear HBs and activated dynamics. The disordered state has a continuous distribution of coordination states, furcated and unstable hydrogens, linear HBs, and diffusive dynamics, similar to normal liquids. PMFs and rate constants of HB switching, supported by radial distribution functions, provide a way to assign HBs and explain their stability and structural origin. Furcated hydrogens are equivalent to free, broken, distorted or dangling hydrogens and involve a hydrogen donating to more than one acceptor. They occur when a hydrogen's nearest acceptor has one more donor than its next-nearest acceptor. They are highly mobile because they spread by barrierless, whole-molecule vibration, bringing about continuous exchange between the ordered and disordered states. They are formed from the non-tetrahedral coordinations brought about by activated HB switches in the ordered state and are consumed by the reverse process. The anomalous low density of water's ordered state and water's ability to alternate between the ordered and disordered states in only a single HB switch, whether diffusive or activated, are central to its continuous, heterogeneous nature.