Tumbling with a Limp: Local Asymmetry in Water's Hydrogen Bond Network and its Consequences

Ab initio molecular dynamics simulations of liquid water under equilibrium ambient conditions, together with a novel energy decomposition analysis, have recently shown that a substantial fraction of water molecules exhibit a significant asymmetry between the strengths of the two donor and/or the two acceptor interactions. We refer to this recently unraveled aspect as the"local asymmetry in the hydrogen bond network". We discuss how this novel aspect was first revealed, and provide metrics that can be consistently employed on simulated water trajectories to quantify this local heterogeneity in the hydrogen bond network and its dynamics. We then discuss the static aspects of the asymmetry, pertaining to the frozen geometry of liquid water at any given instant of time and the distribution of hydrogen bond strengths therein, and also its dynamic characteristics pertaining to how fast this asymmetry decays and the kinds of molecular motions responsible for this decay. Following this we discuss the spectroscopic manifestations of this asymmetry, from ultrafast X-ray absorption spectra to infrared spectroscopy and down to the much slower terahertz regime. Finally, we discuss the implications of these findings in a broad context and its relation to the current notions about the structure and dynamics of liquid water.


Introduction
The many remarkable properties exhibited by liquid water remain subject of intensive on going research. While the exceptional nature of many of these properties warrant investigation in their own right, it is also clear that a better understanding of the molecular basis of the properties of liquid water is of paramount relevance to biology, chemistry, materials science, and geology, just to name a few disciplines. 1,2 It is currently recognized that water is nolonger merely the background against which many reactions, and most biochemical reactions, are taking place. Water is just as important as any of the key (bio)chemical players on any stage. 3,4 Shortcomings in our understanding of liquid water are, ipso facto, limitations to our understanding of the chemistry in aqueous or humid environments. Despite of intensive efforts, many fundamental details are unknown and many questions remain heavily debated. It is remarkable that we basically do not understand exactly how microwaves heat water, 5 nor how heat dispersion through the HB network actually works. Similarly, the microscopic mechanism underlying water's dielectric function in the low frequency region have been under constant debate, so is the basic tetrahedral structure of water, which has been challenged based on some interpretations of water's thermodynamic and spectroscopic characteristics. 6 At the heart of the properties of liquid water is its dynamic hydrogen bond (HB) network with its fluctuating local tetrahedral geometry and collective behaviour. In fact, the remarkable properties of water and also its so-called anomalies can ultimately be traced back in one way or another to this dynamic network when combined with the small molecular size and molecular polarity. 7,8 Because of this, characterizing individual HB dynamics and localized motional modes is far from sufficient for a full understanding of the complexities of the underlying dynamical processes. The collective and emergent nature of the ensuing properties simply prohibits such a naive reductionist view. As a drastic example of this, it is enough to recognize that while it is true that water without HBs would have been a gas at ambient conditions, it is also true that without the HB cooperativity, water would have been a very viscous liquid that would hardly be similar to the matrix of life that we know. This is because the average HB energy in liquid water is an order of magnitude higher than thermal fluctuations at room temperature 9 and it is only through the collective motion that water flows as it does. It is also because of these cooperative modes that the dynamics of water remains a difficult problem for analytical theory that cannot be adequately handled by molecular hydrodynamic theory or continuum-model-based theories, which is why atomistic computer simulations are usually required to interpret experimental findings. Indeed, simulations have provided insights that were crucial for the interpretation of some of the most important spectroscopic investigations of water. [10][11][12] In this article, we are going to review some of our efforts in unraveling a particular aspect of water's HB network, which is the local asymmetry in the strengths of HBs donated or accepted from/to a water molecule in bulk liquid water under ambient conditions. 9,13,14,21,40,43,48,61 While it is well known that the HB strength in water or in ice exhibit a broad distribution trivially due to thermal fluctuations, the novel and non-trivial aspect that was revealed is that in liquid water (but not in hexagonal ice), thermal motion leads to a significant population of molecules, where the strength of the two donor or acceptor interactions have an extreme disparity with as much as 25% of the molecules having one donor (acceptor) interaction that is more than six times stronger than the other donor (acceptor) interaction, and a significant pop-ulation having both kinds of interactions highly asymmetric. This interesting feature of water's HB network was revealed using ab initio molecular dynamics (AIMD) simulations ? ? combined with a condensed phase energy decomposition analysis based on absolutely localized molecular orbitals (ALMO-EDA). 13,21 The strength of this method is that it permits a rigorous quantification of the amount of charge-transfer by locating the variationally lowest energy state without charge-transfer. 9 . Thereby, the issues of under/overestimation of charge-transfer and contamination due to charge overlap effects are circumvented. 15 The decomposition of the collective interaction energy in the condensed phase into physically meaningful components revealed this significant instantaneous asymmetry within the strength of the local donorâȂŞacceptor contacts.
In the following sections of this review we will discuss how this novel aspect was first revealed, its static aspects pertaining to the frozen geometry of liquid water at any given instant of time, and its dynamic characteristics pertaining to how fast this asymmetry decays and the kinds of molecular motions that brings about such decay. Following this we discuss the consequences of this asymmetry, in particular how it shows up or influences different kinds of spectroscopy from ultrafast X-ray absorption (XA) spectra to line shapes in infrared (IR) spectroscopy and down to the much slower terahertz (THz) regime. Finally, we discuss the implications of these findings in a broad context and in particular, its relation to the current notions about the structure and dynamics of liquid water.

ALMO-EDA and the Discovery of the Instantaneous Asymmetry in Liquid Water's Hydrogen Bond Network
2.1 ALMO-EDA: Energy decomposition analysis based on absolutely localized molecular orbitals Intermolecular bonding is the end result of a complicated interplay of multipolar electrostatic interactions, polarization effects, geometric distortions, Pauli repulsion, and charge-transfer interactions (also known as covalent, delocalization, or donor-acceptor orbital interactions). 15 The goal of an electron-decomposition analysis (EDA) is to isolate these physically relevant components from the total energy of an interacting system composed of several subunits. In case of liquid water the subunits are of course the individual water molecules, and the total energy is, for example, that obtained from an ab initio calculation (wave function-based or density functional-based) for a system composed of many interacting water molecules. In this case, what is of interest to us is bulk liquid water as simulated using periodic boundary conditions. The EDA method based on ALMOs, has the key advantage that it permits the calculation of the energy lowering associated with the electron transfer due to hydrogen bonding, using a self-consistent approach that includes cooperativity effects, which are essential for a correct description of the HB network. The ALMO-EDA method works by separating the total interaction energy (∆E T OT ) into molecular frozen density (∆E FRZ ), molecular polarization (∆E POL ), and pair-wise charge-transfer (∆E CT ) terms, i.e.
Therein, (∆E i FRZ ) is the energy change in each molecule i within the system, which is brought about by bringing infinitely separated distorted molecules (at their given geometries in bulk water) into the bulk without any relaxation of the electronic structure of the molecules. This roughly parallels the idea of atomization energy, except that the building blocks in this case are not individual atoms, but molecules with all their internal degrees of freedom (nuclear positions and electron density) totally frozen. Thus, (∆E i POL ) is due to the intramolecular relaxation of each molecule's electrons in the field of all other molecules in the system. These two terms are molecular and additive, i.e. the total contribution to the energy of the system is simply a sum over all the molecular contributions. The ALMO method is designed such that this relaxation constrains the molecular orbitals to remain strictly localized on their respective molecules, hence the name absolutely localized molecular orbitals. 15 Finally, the contribution from parwise two-body charge-transfer interactions (∆E CT ) is computed as the energy difference between the relaxed (i.e. polarized) ALMOs and the fully optimized and delocalized molecular orbitals, which are obtained by a full self-consistent field calculation of the whole system. Hence, ∆E i, j CT accounts for the energy lowering due to electron transfer from the occupied ALMOs on one molecule i to the virtual orbitals of another molecule j. In the current implementation this term is computed using a single noniterative Rayleigh-Schrödinger (RS) perturbative correction ∆E RS CT starting from the converged ALMO solution, i.e. neglecting a much smaller second-order energy change due to induction that accompanies such occupiedvirtual mixing. Detailed descriptions of the ALMO-EDA terms with their mathematical details and algorithmic implementations have been given by Khaliullin et al. 9,13,[15][16][17]21 ) Most importantly for the present context is that the RS perturbative occupied-virtual electron transfer can be further decomposed into forward and back-donation components for each pair (i, j) of molecules: Such a transfer of electrons from an acceptor to a donor is characteristic for hydrogen bonding, where a HB acceptor mostly acts as an electron donor, and vice versa. 19 The charge-transfer term is also, by definition, directly related to the covalency of a hydrogen bond. 15,20,21 Moreover, we have also shown that the chargetransfer term can be accurately estimated using proton nuclear magnetic resonance spectroscopy ( fig. 1), providing a means to experimentally quantify this important aspect of intermolecular bonding. 21 The elegance of ∆E CT (and also ∆Q CT , the fractional charge-transfer that is responsible for the ∆E CT term) is that these metrics are solely based on the two-body delocalization-energy and charge-transfer components, which are readily computable for any given molecular configuration, and as such can be directly used to study static HB networks, as well as the HB rearrangement dynamics without incorporating any (arbitrary) geometry-based criterion of a HB.

Instantaneous asymmetry in the first coordination shell of liquid water
When the ALMO-EDA technique was used to analyze snapshots extracted from AIMD trajectories of liquid water, the fractional electron transfer through a HB turned out to be a few milli-electrons only, but nevertheless this contributes significantly to the total HB energy and is roughly equal to the contribution of the polarization term. 9,17 Figure 2 depicts the average of the strongest five acceptor and donor contributions to the average delocalization energy of a molecule (∆E CT ). Examination of the depicted pattern of chargetransfer interactions reveals that electron delocalization is domi- nated by two strong donor (accpetor) interactions that together are responsible for 93% of the delocalization energy of a single molecule. The third and the fourth strongest donor (acceptor) interactions contribute only 5% and correspond to back donation of electrons to (from) the remaining two first-shell neighbours. This is to say that a HB donor still weakly donates electrons to the unoccupied orbitals of the HB acceptor, and vice versa. The remaining 2% correspond to the delocalization energy to (from) the second and more distant coordination shells.
Comparison of the strengths of the first and second strongest donorâȂŞacceptor interactions ( ∆E CT ∼ 25 and ∼ 12 kJ mol −1 , respectively) with that of the water dimer (∼ 9 kJ mol −1 ) immediately suggests that each water molecule can be considered to form on average two donor and two acceptor bonds, in agreement with the usual tetrahedral structural picture of the coordination shell in liquid water. What is striking in fig. 2, however, is the substantial difference in the strengths of the first and second strongest interactions, which implies that a large fraction of water molecules experience a significant asymmetry in their local environment. The same pattern is also clear when the energies of the two strongest donors or acceptors are plotted together, which is depicted in the bottom two panels of fig. 6. In a simple-minded picture, where the strength of a HB directly correlates with its length, this instantaneous asymmetry would correspond to a structural picture where, due to thermal agitation, many water molecules are in distorted local tetrahedral environments. We will shortly scrutinize this simple geometric picture in more quantitative rigor, and anticipating the discussion, we assert that things are not that simple regarding the relation between HB strength and the geometry in the first coordination shell.
To characterize the aforementioned asymmetry of the two strongest donor or acceptor contacts of a molecule, a dimensionless asymmetry parameter was introduced, where γ A is used to denote the asymmetry within the acceptor and γ D for the donor interactions. The asymmetry parameters are obviously zero when the two strongest contacts are equally strong (perfectly symmetric) and equals to one when the second contact ceases to exist. The joint probability distribution of the molecules according to their γ asymmetry parameters is depicted in fig. 3 together with the lines indicating the quartiles separating the molecules into four groups of equal sizes. The shape of the distribution demonstrates that most molecules form highly asymmetric donor or acceptor contacts at any instance of time. For example, the line at γ ≈ 1 2 indicates that for 75% of the molecules either γ A or γ D is more than 0.5, which means that for the majority of molecules the strongest donor or acceptor contact is at least two times stronger than the second strongest. Furthermore, the peak in the region of high γ in fig. 3 indicates the presence of molecules with significant asymmetry in HB strengths, both for donated and accepted HBs simultaneously. This is also clear from the line at γ ≈ 5 6 , indicating that 25% of the molecules have the strongest interaction six times stronger than the second-strongest.
To understand the origin of this asymmetry, the geometry of donorâȂŞacceptor pairs involved in the first and second strongest interactions was compared. It was found that the strength of the interaction is greatly affected by the intermolecular HB distance R = d(O D − O A ) and to a weaker extent by the HB angle fig. 4 for a depiction of R and β ), while the other geometric parameters have only minor influence on ∆E CT . The histograms of ∆E CT , R, and β for the two strongest donor interactions are shown in fig. 5. The strong overlap between all the histograms suggests that some second strongest interactions have the same energetic and geometric characteristics as the strongest contacts. This implies that the observed electronic asymmetry cannot be attributed to the presence of two distinct types of HBs âȂŤ weak and strong! It is rather a result of continuous deformations of a typical HB. Another important conclusion that can be made from the distributions in fig. 5 is that relatively small variations of the HB distance ∆R ∼ 0.2 Å and angle ∆θ ∼ 5 • to 10 • entails remarkable changes in the strength and electronic structure of HBs. Specifically, the average intermolecular oxygen-oxygen distances The normalized probability density function of the asymmetry parameters γ A and γ D , respectively. The probability of finding a molecule in a bin can be found by dividing the corresponding density value by the number of bins (that is, 400). The dashed black lines at γ ≈ 1 2 , 2 3 , 5 6 partition all molecules into four groups of equal sizes. Reproduced from Ref. 13. for the strongest and second strongest interactions differ only by ∆R ∼ 0.2 Å and are not large enough to justify asymmetric models of liquid water. This is a sober reminder that with the highly nonlinear relation between the HB strength and HB geometry, simpleminded projections of the features of one onto the another can be very misleading, and a high degree of asymmetry in the former does not imply a dramatic distortion in the latter. Furthermore, analysis of the structure of the molecular chains defined by the strongest bonds (that is, one donor and one acceptor for each molecule) shows that their directions are random and without any long-range order (e.g. rings, spirals or zig-zags) on the length scale of the simulation box (∼ 15 Å).
Because of the slow decay of the distribution tails in fig. 5, a quantification of the concentration of single-donor and singleacceptor molecules was not attempted. Defining such configurations using a distance, angle or energy cutoff is an unavoidably arbitrary procedure. A quantitative analysis of the network, which was performed for the same molecular configurations than in Kühne et al. shows that according to the commonly used geometric definitions of hydrogen bonds, 23-25 the structure of water is distorted tetrahedral with only a small fraction of broken bonds.
Finally, with an asymmetry in liquid water that is due to thermal distortions, it is very interesting to compare this with the results of ALMO-EDA in hexagonal ice, which is generally regarded to be a solid with symmetric HBs. Such a comparison has been performed and indeed does provide further insight into the ori-  gin of the asymmetry in liquid water and its relation to the HB network structure. 14 Figure 6 shows the joint distribution of the strengths of the two pairs of donor and acceptor HBs in ice. The ensuing two-dimensional distribution is characterized by the peak centered at −18.7 kJ mol −1 , large deviation from the average values (σ = 7.1 kJ mol −1 ), and a small correlation coefficient of 0.1. Such a tiny correlation coefficient indicates that the two HBs are essentially independent from each other and the asymmetry in ice (points in the distribution that are significantly displaced from the diagonal line) arises trivially from the very broad distribution of HB strengths. As in the case of liquid water, this asymmetry is a result of thermal fluctuations around the average symmetric structure. The distribution of the strength of the donor interactions in liquid water exhibits a drastically different pattern with two pronounced features. In addition to a broad peak resembling the one for ice, there is a sharp peak in the region of high γ D . The center of the first peak is shifted to lower energies (−12 kJ mol −1 ) and is somewhat broader than that of ice. The second peak indicates the presence of molecules with one intact and one broken donor HB (∆E C→N > −1 kJ mol −1 ). To estimate the fraction of molecules responsible for the sharp asymmetric feature, we draw a somewhat arbitrary boundary at γ D = 0.8 (dashed line in Figure 6), which divides the distribution into the regions of ice-like configurations and highly asymmetric configurations. Figure 6 shows that 26% (18%) of molecules in the liquid phase are characterized by γ D > 0.8 (γ A > 0.8) compared to 2% in ice.
In hexagonal ice, the uncorrelated thermal motion of molecules around their crystallographic sites broadens the range of HB ener-gies, and thus creates a noticeable asymmetry in the donor and acceptor contacts of each water molecule. This is trivially a broad distribution without any particular correlation between the strengths of the pair of donor (or acceptor) interactions at any center. In water, while the majority of molecules still exhibit HB patterns similar to those in ice and retains a four-fold coordination with only moderately distorted tetrahedral configurations, there is the drastic difference of a large fraction of molecules with very weak HBs, which are elongated by as much as 0.5 ÃĚ and exhibits a wide range of angular distortions. A detailed analysis is given in Ref. 14. These results imply that the traditional view of water as a four-fold coordinated nearly tetrahedral liquid is more appropriate than the recently proposed asymmetric model. 6,23,26 However, the substantial fraction of molecules with broken HBs undoubtedly affects the physical properties and chemical behavior of liquid water. Starting from Section 3 we investigate some of the consequences of this aspect of liquid water.

Contrasting the asymmetry in donor versus acceptor interactions
Figure 6 reveals an interesting difference between donor and acceptor interactions in liquid water. While the corresponding plots in ice are identical, they are rather different in liquid water. Thermal disorder affects the strengths of donor and acceptors interactions in a quantitatively different manner. First, while also significant, the fraction of molecules with broken acceptor bonds (18% with γ D > 0.8), is lower than the corresponding fraction in γ A (26%). Furthermore, the distribution of the acceptor interactions does not exhibit a high-γ A peak, which would match the high-γ D peak. This difference indicates that within the broken HBs, only the donor of electrons remains under-coordinated, while the acceptor (i.e., hydrogen atom) forms a HB with another donor that becomes overcoordinated. The existence of a significant fraction of overcoordinated donors is supported by the relatively large contribution of the third interaction shown in Figure 2 and is consistent with the well-known fact that the distribution of electron acceptors around a water molecule is more disordered than that of the donors. 27 This phenomenon is attributed to the existence of the so-called âȂIJnegativity trackâȂİ between the lone pairs of a water molecule, which facilitates the disordered motion of electron acceptors around the central donor. 28,29 This asymmetry in donor versus acceptor interactions is also reminiscent of the asymmetry in solvation energy of anions and cations, ? and a possible relation between both phenomena

Dynamics of the asymmetry
The overlapping distributions in Figure 5 suggest that, despite the difference in the strength of the donorâȂŞacceptor contacts, their nature is similar and that the strongest interacting pair can become the second strongest in the process of thermal motion and vice versa. To estimate the time scale of this process, it is necessary to examine how the average energies of the first two strongest interactions fluctuate in time. The instantaneous values at time t (solid lines in fig. 7a) were calculated using the ALMO EDA terms for 3501 snapshots separated by 20 fs (448128 local configurations) by averaging over time origins τ separated by 100 fs and over all surviving triples: where M(τ,t) is the number of triples that survived from time τ to τ + t. A triple is considered to survive a specified time interval if the central molecule has the same two strongest-interacting neighbours in all snapshots within this interval. Figure 7a shows that the strength of the first two strongest interactions oscillates rapidly and that ∼ 80 fs after an arbitrarily chosen time origin, the strongest interaction becomes slightly weaker than the second strongest (note that first and second refer to their order at t = 0). The amplitude of the oscillations decreases and the strengths of both interactions approaches the average value of ∼ 20 kJ mol −1 on the time scale of hundreds of femtoseconds. The decay of the oscillations indicates fast decorrelation of the timeseparated instantaneous values because of the strong coupling of a selected pair of molecules with its surroundings. In other words, while the energy of a particular HB fluctuates around its average value, this bond has approximately equal chances of becoming weak or strong after a certain period of time independently of its strength at t = 0. This effect is due to the noise introduced by the environment and can be observed in ultrafast IR spectroscopy experiments. 30 The time averages shown in fig. 7 are physically meaningful and can be calculated accurately only for the time intervals that are shorter than the average lifetime of a HB. 31 The small residual asymmetry that is still present after 500 fs ( fig. 7a) is an indication of the slow non-exponential relaxation behaviour that characterizes the kinetics of many processes in liquid water. 24 Specifically, the non-exponential relaxation of the HB lifetime is due to the coupling between HB fluctuations and diffusion, 24 and the non-exponential tail in fig. 7a a likely manifestation of this process. At variance, the fast relaxation component in fig. 7a is more correlated with HB stretch motions and closely matches the time scale of 170fs found from the spectral diffusion of the OH stretch peak in IR spectroscopy. 12,30,33 In addition to the instantaneous values of ∆E C→N (t) at time t, the dashed lines in fig. 7a show the corresponding averages over time interval of length t. These values were calculated by averaging over time origins τ, all snapshots lying in the time interval from τ to τ + t and over all surviving triples: The dashed lines in fig. 7a show that any neighbour-induced asymmetry in the electronic structure of a water molecule can be observed only with an experimental probe with a time-resolution of tens of femtoseconds or less. On longer time scales, the asymmetry is destroyed by the thermal motion of molecules and only the average symmetric structures can be observed in experiments with low temporal resolution. An examination of the time dependence of all two-body and some three-body geometric parameters that characterize the relative motion of molecules reveals the mechanism of the relaxation. The curves in fig. 7a and b show that the asymmetry is almost completely lost on the time scale of a single cycle of the HB stretch motion (∼ 200 cm −1 , which corresponds to ∼ 170 fs). It is worth noting here that the decay of the time correlation function of the instantaneous fluctuations in energy δ E 1 (0) δ E 2 (t) follows the same time scale. Relaxation of the asymmetry is thus primarily caused by low-frequency vibrations of the molecules relative to each other. The minor differences in the behaviour of the curves, in particular at 80 fs, indicate that the relaxation of the asymmetry is possibly also influenced by some faster degrees of freedom. The temporal changes in the HB angles towards the average value ( fig. 7c) show that librations of molecules play a minor role in the relaxation process.
The kinetics and mechanism of the asymmetry relaxation presented here are supported by data from ultrafast IR spectroscopy, which can directly observe the spectral diffusion of the OH stretch peak with a fast component that monoexponentially decays at 170 fs. 12,32,33 In conclusion, it is important to point out that the dynamics of the asymmetry closely parallels the dynamics of the spectral diffusion of the OH stretch peak, with a fast component that is mostly associated with the intermolecular O-O stretch motion, and a longer time dynamics that is associated with the collective HB network restructuring. 34 3 Electronic Signature of the Asymmetry: XA Spectroscopy The time behaviour described above implies that the instantaneous asymmetry can, in principle, be detected by X-ray spectroscopy, which has a temporal resolution of several femtoseconds and is highly sensitive to perturbations in the electronic structure of molecules. 23 To identify possible relationships between the spectroscopic features and asymmetry, the X-ray absorption (XA) spectrum of liquid water was calculated at the oxygen Kedge using the half-core-hole transition potential formalism within all-electron density functional theory. ? ? ? Although the employed computational approach overestimates intensities in the post-edge part of the spectrum and underestimates the pre-edge peak and overall spectral width, it provides an accurate description of the core-level excitation processes and semi-quantitatively reproduces the main features of the experimentally measured spectra. The localized nature of the 1s core orbitals allows to disentangle the spectral contributions from molecules with different asymmetry.
To this end, all molecules were separated into four groups according to the asymmetry of their donor and acceptor environments, as shown in fig. 3. As already mentioned, choosing boundaries at γ = 1 2 , 2 3 , 5 6 distributes all molecules into four groups of approximately equal sizes (that is, 25 ± 2%). Figure 8 shows four XA spectra obtained by averaging the individual contributions of molecules in each group. It reveals that molecules in the symmetric environment exhibits pronounced post-edge peaks, whereas molecules with high asymmetry of their environment are characterized by the amplified absorption in the pre-edge region. Furthermore, the relationship between the asymmetry and absorption intensity is non-uniform: the pre-edge peak is dramatically increased in the spectrum for the 25% of molecules in the most asymmetric group (γ = 5 6 ), for which the strongest interaction is more than six times stronger than the second strongest. As a consequence, the preedge feature of the calculated XA is dominated by the contribution of molecules in the highly asymmetric environments (Figure 8c). The pronounced pre-edge peak in the experimentally measured XA spectrum of liquid water has been interpreted as evidence for the so-called "rings and chains" structure, where ∼ 80% of molecules have two broken HBs. 23,26 These results however suggest that this feature of the XA spectrum can be explained by the presence of a smaller fraction of water molecules with high instantaneous asymmetry. Although the employed XA modelling methodology does not allow a precise estimate of the size of this fraction, the result is consistent with that of recent theoretical studies at an even higher level of theory, which have demonstrated that the main features of the experimental XA spectra can be reproduced in simulations based on conventional nearly tetrahedral models. 35,36 Thus, the investigation of the relation between local inhomogeneities in the HB network and the XA spectrum revealed an interesting and important connection between relatively small geometric perturbations in the HB network and the large asymmetry in the electronic ground state, as well as the XA spectral signatures of the core-excitation processes. This helps to reconcile the two existing and seemingly opposite models of liquid water -the traditional symmetric and the more recently proposed asymmetric model.

Vibrational Signature of the Asymmetry: Inhomogeneous broadening of the O-H stretch peak
On the one hand, the structural picture obtained from XA spectroscopy essentially corresponds to snapshots of the system that are frozen on the time scales of nuclear motion. Raman and IR spectroscopies, on the other hand, directly probe the resonance frequencies of the vibrational motions of the ionic cores of the atoms. Moreover, the intrinsic time resolution of time-resolved IR spectroscopy is about 50 fs. 37 It is also well known that the strength of a HB is correlated with the frequency of the covalent OH bond  stretch motion, with stronger HBs being associated with a red-shift in the frequency. 38,39 In fact, it is mainly because of the disorderinduced inhomogeneities in HB strength that the line width of the OH stretch peak in vibrational spectra of liquid water are much broader than that in the gas-phase, i.e. they are inhomogeneously broadened with some contribution from homogeneous broadening as well. 10,12,38 Furthermore, there seems not to be a substantial extent of motional narrowing and thus the total linewidth is close to the inhomogeneous limit. 38 Taken these features of the OH stretch peak together, one can immediately conclude that the instantaneous asymmetry might exhibit some observable effects within vibrational spectroscopy. In fact, it is tempting to say that vibrational spectroscopy, and its lower frequency sibling, THz spectroscopy, are the natural techniques to probe any structural or dynamic manifestations of the local HB asymmetry. As already alluded to above, ALMO-EDA offers a powerful and consistent way to quantify the strength of HB interactions in liquid water, and hence to quantitatively establish the aforementioned relationships between the structure and dynamics of the HB network and its characteristics on the OH IR and Raman peaks, as obtained from linear, non-linear, and time-resolved techniques. 10 Along these lines of reasoning, ALMO-EDA was used to investigate the relationship between the vibrational fluctuations of the OH stretching modes and the strength of the HB, as quantified using the ∆E CT ( fig. 9). 40? The instantaneous fluctuations in frequency of the OH modes were calculated using a wavelet-based timeseries analysis. 41,42 The clear correlation in fig. 9 indeed shows that ∆E CT , taken as a measure of HB strength, linearly correlates with the frequency of the covalent OH stretch vibration. A linear regression yields which can be used to estimate the average charge-transferassociated stabilization of a HB, corresponding to a given OH frequency of the electron-accepting water molecule. The root mean square error of the fit is 7.27 kJ mol −1 . Such a large dispersion of the HB strength in relation to the OH vibrational frequency (and vice versa) agrees well with previous findings regarding the relation between O-O HB length and the OH frequency. 34 It reaffirms our previous remark that the relation between the HB strength and the local environment cannot be trivially inferred from a single geometric or dynamical variable. It is also important to note that for similar reasons, the dispersion of the observed data points varies with the frequency/HB strength, stronger HBs exhibit a higher variation in the associated OH frequencies.
Most importantly, the linear trend in fig. 9 suggests a mean to unravel the consequences of asymmetry on the inhomogeneouslybroadened OH stretch peak. In order to see how this is possible, let us now describe the two O-H stretch vibrations on a single water molecule as two coupled harmonic oscillators which we denote as O-H1 and O-H2, with the same mass m and restoring force constants k 1 and k 2 that are coupled by the intramolecular harmonic coupling constant k (see fig. 10). If k 1 = k 2 , as would be the case for water in the gas phase, or for a symmetric HB environment in condensed phase, then the coupled normal modes of the system ν 1 = [1, 1] and ν 3 = [−1, 1] are the familiar symmetric and asymmetric stretch modes. The energy splitting between the two eigenmodes is given by is used to denote that k 1 = k 3 in this case). We can see in this case that both modes equally contribute to the two coupled normal modes, i.e. the two O-H stretch vibrations are perfectly and symmetrically coupled, and the hydrogen atoms move together with the same amplitude. Harmonic mode analysis on a water monomer in vacuum using DFT gives a splitting ∆ω 13 = 104 cm −1 at 0 K ( fig. 11A). 43 For comparison, CCSDT(T) calculations at the complete basis set limit yield ∆ω 13 = 109.8 cm −1 . 44 However, when k 1 = k 2 , the two normal modes are ν 1,3 = [δ ± √ 1 + δ 2 , 1], with a splitting ∆ω 13 ∝ (1 + δ 2 ) where δ = (k 2 − k 1 )/2k . 43 The asymmetry in the restoring force leads to a decoupling of the two O-H stretch vibrations, and the degree of decoupling depends on the difference between k 1 and k 2 . In fact, ab initio calculations reveal that the decoupling of the two O-H stretch modes in the HB donor of the water dimer is almost complete, where each mode being composed roughly of 85% from one O-H stretch, with only a 15% contribution from the other O-H being mixed in. 45 In this case normal mode analysis shows, as expected, that the frequency splitting increases to 221 cm −1 (fig. 11C).
The simple model of two coupled harmonic oscillators shows that one consequence of the local asymmetry in the HB environment is to decouple the two intramolecular O-H bond vibrations, as well as to increase the frequency splitting between the resultant symmetric and the asymmetric stretch modes. Indeed, it becomes possible to quantify the HB asymmetry with the degree of decoupling f d between the two O-H bond vibrations where C KO−H1 and C KO−H2 are the contributions of the O-H1 and O-H2 bonds to the normal mode k ∈ {ν 1 , ν 2 }, respectively. When f d is equal to 0, the two O-H bonds are completely coupled in their motions. An increase in f d , just like an increase in the frequency splitting, is a signature of water molecules in asymmetric HB environments and, according to the previous findings, should be manifest in liquid water. This line of reasoning has been tested with AIMD simulations of liquid water. 43 To this end, a pre-requisite is to extract the localized normal modes (ν 1 and ν 3 ) from the total vibrational density of states of the bulk liquid, because the latter is the quantity that is directly accessible from MD simulations. One way to achieve this is by requiring maximal localization of the power spectra of the two local modes in the frequency domain. If we take the localization criterion as a minimization of the quantity ω 2n − ω n , where n is an integer constant, then this leads to the requirement of minimizing the following functional 46 : with respect to ν k , where β = 1/k B T , and Pν k is the vibrational density of states of ν k . It can be shown that for n = 2, this method is equivalent to a normal mode analysis performed on the thermally averaged Hessian matrix and generalized to anharmonic systems at finite temperatures, and that the zero-temperature limit yields the usual normal modes. 46 With this procedure, the two normal modes of each water molecule, ν 1 and ν 3 , were extracted from MD trajectories at finite temperatures. A splitting of 103 cm −1 was obtained from MD trajectories of a single water molecule in vacuum at 300 K ( fig. 11B). The essentially identical results obtained at 0 K (using the familiar normal mode analysis) and 300 K (from MD trajectories) not only validate this approach, but they further indicate that temperature and anharmonicity effects alone have negligible influence on the mode coupling. The average value of the parameter f d turns out to be 0.07, in clear contrast to the water dimer case with f d = 0.92, thus indicating an almost full decoupling between the two O-H stretch vibrations of the HB donor molecule, as already pointed out. Finally, in the case of liquid water at 300 K, f d turns out to be 0.82, with an average splitting of 137 cm −1 (fig. 11D). In this case, the modes resembles those of the water dimer and, interestingly, also share similarities with the instantaneous normal modes of water molecules at water/vapor interfaces. 47 Interestingly, it was found that for f d ≈ 1 and ∆ω 13 > 400 cm −1 , the dipole moment of water molecules shifts to a lower value by 0.25 D with respect to the fully symmetric case in bulk water. This resembles interfacial water molecules in water/vapor systems. 48 By empirically fitting the data obtained from MD, it was also shown that the relation between f d and the frequency splitting is closely approximated by ∆ω 13 , where C is an empirical fitting parameter. This relation holds in the range f d = 0 until f d ≈ 0.8, after which, the highly asymmetric configurations of water exhibit a broader distribution in ∆ω 13 and a skew towards higher frequency splitting. For the configurations, where a HB is totally broken, ∆ω 13 increases up to more than 400 cm −1 .
Although the asymmetry in HB strength seems to be the strongest factor in determining the frequency splitting, and hence, the inhomogeneous broadening of the O-H stretch peak in liquid water, other factors exhibit a role as well. In a later investigation the role of the intramolecular coupling between the two O-H stretch modes (k ) was scrutinized. 48 It was found that k does modulate the frequency splitting, in a reverse fashion to f d , with the overall effect that the observed frequency splitting in liquid water is less than the value that is predicted solely based on the effect of f d . As a consequence, the instantaneous asymmetry of the local HB environment around a water molecule in bulk water was shown to account to the observed inhomogeneous broadening of the O-H stretch peak. While the role played by the finitetemperature distribution of HB strength in the O-H stretch, inhomogeneous broadening in neither surprising, nor a novel finding (see Ref. 39 for instance), instead, the new aspect here is that the distribution of HB strength is associated with the peculiar feature of mostly coming in pairs of strong-weak HBs, which are localized on the water molecules. A feature that is most clearly seen on the right side of fig. 12 as f d approaches 1.

Enhancement of Asymmetry Under External Electric Fields
The study of water under electric fields is immensely interesting for several reasons. First, water does exist under moderately to substantially strong electric fields in a variety of natural settings, for example within biological membrane channels, in the vicinity of electrodes and ions in solution, and in cracks at crystal surfaces. The electric field intensity in some of these cases can be of the order of 10 9 V m −1 . [49][50][51] Another reason of interest for studying water under electric fields, is that the electric field-induced anisotropy can give rise to new interesting features 52? ? ? -54 and can also enhance the equilibrium structural/dynamical features in a liquid, facilitating their study and hence providing new ways to understand the complexities of liquid state kinetics and thermodynamics. 55 A good case for the latter is the long-established study of dielectric relaxation both experimentally and theoretically, 56,57 and for the former, the discovery of non-vanishing rotational-translational cross-correlations in water under electric fields. 52,58? External electric fields, both static and fluctuating, are known to induce a variety of field-induced anisotropies in liquid water. 59,60 In the context of HB asymmetry, the significance of application of an external electric field lies in the possibility of exploiting these fieldinduced anisotropies to enhance the asymmetry, thus for instance to facilitate its experimental investigation. Naturally, in this case, the origin of the -possibly enhanced -anisotropy would not be the thermal fluctuations, but rather the external field itself.
The effect of an electric field on the local HB asymmetry in water was investigated with AIMD using an intense electric field square pulse with an intensity of 4.3 × 10 9 V m −1 . 61 The pulse strength was chosen such that it induces an ultrafast re-orientation of the water molecules on a sub-picosecond time scale. Under these conditions, it was shown that within 300 fs, the water molecules reach a steady-state average orientation angle of 37 degrees ( cos(θ ) ≈ 0.8), where the angle is calculated between the molecular bisector and the direction of the external field. This ballistic re-orientation of the water molecules substantially increases the temperature of the system up to 650 K. Once the pulse is switched off, this orientational anisotropy decays exponentially and vanishes within 750 fs. The influence of the pulse on the asymmetry of the HB network is shown in Figure 13, which depicts the joint probability distribution of the asymmetry parameters γ A and γ D at various times after the pulse. To distinguish the electric field induced effects from effects that are only due to the high temperature of the system (in particular the drop in the HB density), the joint probability distribution was compared to that found in a field-free microcanonical trajectory simulated at an average temperature of 650 K (rightmost plot of Figure 13). We see in Figure 13 that immediately following the pulse, the probability distribution has its peak at the top right corner of the plot, where the molecules exhibit a high level of asymmetry simultaneously in the two asymmetry parameters. The asymmetry patterns in fig. 13 are very distinct from the situation in liquid water under ambient conditions, where the largest population of molecules exhibits high asymmetry in one, but not in both asymmetry parameters. This is also distinctively different from the field-free situation in hexagonal ice, where the asymmetries in both parameters are uncorrelated. 14 Comparison to field-free conditions shows that the electric field appreciably enhances the asymmetry. This enhanced asymmetry then gradually decays once the field is switched off, so that after 1 ps, the joint distribution has almost fully relaxed to the situation found in the high-temperature trajectory. The enhanced HB asymmetry under an electric field can be explained by the field-induced anisotropy in HB strength. It was found that the pulse leads to a strong anisotropy in the HB orientation along with the field-induced molecular re-orientation. Moreover, the better a HB is aligned with the field, the shorter it is. This arrangement makes charge-transfer through the HB more favorable as it corresponds to a lowering in the electric potential along the HB, and thus leads to an overall strengthening of the HB. These findings agree with the directional pair density distributions that were recently reported from simulations of water under electric fields. ? ? ? After the pulse, the HB orientational anisotropy decays on the same time scale as the molecular one. Consequently, it was found that the strongest acceptor or donor interaction is typically pointing along the field axis, while the weaker one is more-orless in the orthogonal plane. The molecules that are simultaneously engaged in two HBs and exhibit a high degree of asymmetry (γ > 0.8), are those simultaneously donating (or accepting) one HB in parallel and another one in an orientation that is more-or-less orthogonal to the field.

Water under a single THz Pulse: Influence of asymmetry on molecular re-orientation
For a long time, one obstacle that has hindered a better understanding of water's HB network has been the absence of an experimental technique that directly probes this network. The pico-to sub-picosecond lifetimes of HBs 31 are too short for the NMR and dielectric spectroscopy time window and is only indirectly accessible by time-resolved IR spectroscopy. ? The recent technological feasibility of intense THz laser pulses has lately changed this state of affairs, opening new possibilities for the direct coherent excitation and control of the intermolecular and collective HB modes in water. ? ? ? ?
Leveraging these experimental advances, we have recently investigated the (sub)picosecond pathways of energy transfer from a THz pulse to water, using a novel THz experimental setup ( fig. 14) in combination with MD simulations. ? So far, the physical nature of such fast relaxation processes within the HB network of water has been poorly understood and is heavily debated. 57 Our study has elucidated that the energy of a single THz pulse with a frequency of 1 THz couples mostly to the rotational dynamics of the water molecules, with 85% of the energy going directly to rotational motions and only 15% of the energy going into restricted molecular translational motion. The highly efficient rotranslational coupling in liquid water, coupled with the rigidity of the liberational potential (15 -20 THz), then leads to an increase of the molecular translational kinetic energy of water molecules after the pulse, which lasts for ∼ 1 ps and decays on the same time scale as the observed THz Kerr effect. Interestingly, THz spectroscopy can also shed some light on the local asymmetry of the HB network. When we examined the relation between fieldinduced molecular re-orientation, and the asymmetry parameters γ, we found that the molecules with high asymmetry (γ A/D > 0.8) are orientationally more labile than the molecules with low γ (γ A/D < 0.2). The average angle between the molecular bisectors and the field axis is depicted in Figure 15 for high-and low-γ, where γ was calculated at the time corresponding to the peak THz field intensity (10 8 V m −1 ), denoted as t = 0 on the x-axis. As expected from the lifetime of the asymmetry, the impact on the orientational dynamics can be observed for a few hundred femtoseconds before and after the point in time where γ is calculated. We believe that focusing on HB strengths and asymmetries in these strengths, can be a more fruitful pathway than trying to identify geometric defects in the structure of water, and can help us understanding the relationship between orientational relaxation of a single molecule and the collective relaxation, a problem which remains at the frontier in studies of orientational relaxation. 57,63 Fig. 15 Dielectric alignment of water molecules with high (black) and low (red) asymmetry parameter under a single THz pulse. High-γ is defined as γ > 0.8, whereas low-γ corresponds to γ < 0.2. The average angle between the molecular bisector and the field axis is denoted as cos θ . The vertical dotted line marks the reference point in time at which γ was calculated and corresponds to the peak field in the THz pulse. The THz field profile is shown in the top panel.
small population of mobile fast re-orienting water molecules to explain dielectric relaxation, 57,67? and in various hypotheses of local structural defects and their propagation. 68,69 In all these cases, the core assumptions are fundamentally about the structure and dynamics of the HB network on different time scales, its ability to spatially and temporally sustain local defects, the life times of these defects, and their spatial propagation through the network. This is to say that, as a means towards fruitful progress, a possible falsification of any of these hypotheses must be ultimately based on a better understanding of the dynamics of the HB network. In this regard, it is unfortunate that in many cases the indiscriminate use of descriptive words like "structure" and "heterogeneity" can lead to much confusion if not accompanied by clear operational definitions and a context under which these definitions are valid. The literature on water structure and dynamics is full of descriptions such as "enhanced structure", "structure formation/breaking", and "dynamic polymer", just to name a few. Such explanations can be useful and insightful, but they also have the potential to quickly become mere issues of semantics and nothing more. We agree with previous criticisms that the indiscriminate use of descriptive words in the water literature has only lead to much confusion. 70 In this article, we have discussed an aspect of the HB network that has been recently unraveled using a combination of AIMD together with ALMO-EDA. 9 We have shown that in equilibrium liquid water, at any instance of time, a large population of water molecules have a high disparity between the strengths of the two HBs they are donating or the two HBs they are accepting, an aspect which we refer to as "local asymmetry in the HB network". Here, the strength of the HB is quantified by the amount of charge-transfer from the HB-acceptor to the HB-donor and the associated energy lowering (∆E CT ), both quantities are obtained from ALMO-EDA. We have also shown that the extent of asymmetry can be quantified using the dimensionless asymmetry parameters γ A (electron acceptor interaction) and γ D (electron donor interactions), where γ runs from zero (fully symmetric HB envi-ronment with two equally strong HB interactions) to one (one HB is totally broken). Thus, we provide operationally defined metrics that can be consistently employed on simulated water trajectories to quantify the heterogeneity in the HB network and the dynamics of this heterogeneity.
The idea that there is a broad distribution of HB strengths in liquid water is neither novel nor surprising, and indeed our quantitative metrics of HB strength confirms this picture, but furthermore, the local asymmetry presents a stronger statement. It is neither necessary nor self-evident that a broad distribution of HB strengths implies a significant population of water molecules with a strong asymmetry in both asymmetry parameters simultaneously (75% of molecules with γ A and/or γ D greater than 0.5). This is particularly clear when comparing the picture in liquid water to that in hexagonal ice. But again, it is important to emphasize here that the picture we find is that of a continuous distribution of HB strengths rather than any kind of two-state picture. When we discuss the behavior of high-versus low-γ water molecules, it is meant to emphasize how the extremes of the distribution are behaving, an aspect which is particularly important given the significant abundance of molecules with high γ, but this is not to imply that water is a mixture of two kinds of molecules or HBs. The comparison with ice is once again illuminating in that it further corroborates that the high asymmetry in water is a consequence of an interplay of thermal disorder, the high connectivity and cooperativity of the HB network, and the diffusive aspect of HB dynamics.
Regarding the traditional tetrahedral structure of liquid water, it is trivial to realize that in a liquid system, what might be disordered on some short time scale, becomes ordered when viewed over longer time scales, possibly just as it is no more -or lessprofound to realize that a liquid below some relaxation time can be viewed as a solid glass and can similarly support transverse and longitudinal phonon modes. [71][72][73] Keeping in mind the dynamics of the asymmetry and how it decays, the strength of donorâȂŞacceptor interactions we find, suggest that each molecule in liquid water at ambient conditions forms, on average, two donor and two acceptor bonds. It is because of the very strong dependence of the HB energy on the local geometry (in particular the exponential dependence on the HB length) that even small thermal distortions in the tetrahedral HB network induce a significant asymmetry in the strength of the contacts causing one of the two donor (acceptor) interactions to become, at any instance of time, substantially stronger than the other. Thus, the instantaneous structure of water is strongly asymmetric only according to the electronic criteria, not the geometric one. Overall, the picture provided by ALMO-EDA does not warrant any departure from the traditional tetrahedral structure. With respect to the dynamics of the asymmetry, we have shown that it decays on a time scale of several hundred femtoseconds. Intermolecular vibrations (O-O stretch) and librations of OH groups of HBs are primarily responsible for the relaxation of the instantaneous asymmetry. The time scales, which we find, closely match those obtained from studying the time-dependence of the OH spectral diffusion. The long-time non-exponential tail of the relaxation seems to be related to the non-exponential behavior of HB kinetics, which can be traced back to the translational diffusional aspect of the kinetics. 24 So what are the consequences of this asymmetry on the spectroscopic observables of liquid water? The challenge in figuring out the answer is that the causal relation between the HB energy (and its asymmetry) on the one hand, and the structure and dynamics of water (e.g. order parameters, power spectra, HB lengths and angles and their distributions), on the other hand, is very far from trivial. We know for instance that the removal of some of the HB partners in liquid water has the potential to slow down water rotations, as it breaks cooperativity, but the removal of all HBs, possibly by dispersion in non-polar solvents, can greatly speed up the rotation of water molecules. 70,74? ? So what happens when a water molecule is simultaneously donating a strong and a weak HB? This is somewhat like a water molecule that is tumbling on one strong and one weak leg. It turns out that in XA spectroscopy, this results in a population of water molecules that gives an amplified response in the pre-edge peaks. In case of the vibrational O-H stretch peak, it turns out that the ensuing decoupling between the symmetric and antisymmetric normal modes can largely explain the observed inhomogeneous broadening, whereas in pulsed THz spectroscopy, it turns out that the molecules with high asymmetry are orientationally more "labile" on a time scale close to half a picosecond. Thus, we see that the manifestations are present on different time scales spanning several orders of magnitude. Hence, trying to answer the question posed in the beginning of this paragraph, this can only be considered as work in progress. There are many more properties of water where asymmetry can play an important role, which remains to be investigated. Examples of the latter are attempts to explain the nature of energy dissipation processes, ? and more generally, of fast relaxation processes in water, aspects that are still poorly understood and intensely debated. 57,67? Recent findings have also pointed towards an important role of intermolecular modes in vibrational energy relaxation, a role particularly played by strong HBs. 33 It would also be very interesting to see how asymmetry might dictate certain energy dissipation pathways in this case.
We owe much of what we know about liquid water to interpretations of experimental findings that were founded on force field MD simulations. 75? ? ? ? ? ? , 76 It is an intriguing question whether any of the existing interpretations would be significantly modified if a force field that captures the asymmetry aspect were employed. Recently, Naserifar and Goddard 77 have published a work, where they describe the results of MD simulations of water using a new force field parametrized using a combination of DFT and CCSD(T) benchmark calculations. Interestingly, they find that each water molecule on average has two strong and two weak HBs (based on a distance criterion), and that the relaxation time between a strong and a week hydrogen bond is ∼ 100 fs. Despite of our agreement with the criticism that has been raised against this work 78 , these finding are intriguing and motivate further careful investigation and comparison with previous experimental and theoretical work.

Conflicts of interest
There are no conflicts to declare.