The nature of hydrogen bonding (H-bonding) manifests the delicateness of molecular physics and is a cornerstone to implement various chemical processes and biological functions.^[1–3] Advancing explorations are unveiling this mystery with state-specific spectroscopy,^[4] real-space identification^[5] and direct observation of a single hydrogen bond.^[6–8] Illustrated phenomenon including dimer dissociation, tunneling and quantum nuclear effect (QNE) further emphasizes the dynamic aspect of H-bonding, which is intrinsic due to the zero-point vibration (ZPV) and is the underlying behaviors of all H-bonded networks. While experiments give the final picture containing all effects and inevitable turbulence from environment and observation, theoretical simulations facilitate the construction of a step-by-step roadmap towards the dynamic nature of H-bonding. Particular physical effects could be abstracted using ideally isolated models through hypothetically classified theoretical levels, enable interpretations like the vibrational origin of QNE in deuterated H-bond,^[7,9,10] the possible instantaneous H-bonding under ZPV,^[11] and the ZPV-allowed dissociation path for water dimer cations.^[12] However, as the base bearing effects at upper theoretical levels and attractive processes in complex cases, the quasi-classical ZPV dynamic effect on an elementary ground-state H-bond model is still unclear. Consequently, in ZPV fluctuation, whose normal modes dominate its geometrical variance, how its electronic structure characteristics would respond and even in what range the interaction energy might change are still problematic, leaving incomplete theoretical templet for understanding ground-state nature of all H-bond systems.

As one of the most elementary H-bond archetypes,^[13–16] water dimer has been systematically studied to provide an ideal single-bond reference for the geometric parameter, electronic structure, and interaction strength.^[17–20] Experimentally, geometric structure of a gas-phase H-bonded water dimer was first characterized by distance between two oxygen atoms (R_O…O), which is 2.976 Å including an estimated +0.03 Å contribution assigned to the classical vibrational averaging effect.^[21] Theoretically, high-precision ab-initio calculations predicted equilibrium R_O…O to be around 2.91 Å^[22,23] and the remaining inconsistence about 0.03 Å between observation and prediction could be complemented by adding NQE to dynamic simulations.^[9] While this NQE elongation has been attributed to motions in H-bond bending and covalent bond stretching, the infrastructural quasi-classical understanding on intrinsic ZPV effects is unexpectedly lacking, where bending, stretching and all other motions could be intuitively analyzed within the normal-mode discourse. Such a discontinuity in fundamental models would hinder comprehensive understandings towards practical dynamic behaviours of water systems, where ultrafast vibrational processes and a mode-specific resolution could be crucial.^[24,25] Similar deficiencies also present for understanding the H-bonding electronic structure and interaction energy, where further insight may reside. On one hand, many electronic-structure probes have been established to illustrate the characteristic of H-bond including the covalent-like feature,^[26–33] but only few criteria like NMR shielding parameters^[10,32] considered the dynamic influence at a general view under finite temperatures and bulk phase. On the other hand, theoretical energy decomposition analysis (EDA) was often conducted on a quasi-static model of water dimer, either in equilibrium structure or during an artificial manipulation of chosen freedoms.^[34–43] These strategies overlook the fact that a ground-state water dimer intrinsically vibrates with ZPVs, thus its more practical bonding nature would correspond to a vibrationally varying electronic structure as well as interaction energy and physical components within. Besides the significance in understanding an elementary H-bond, water dimer is also a key cluster to explain the dynamic behaviours of water systems such as interface diffusion^[44,45] and a benchmark model for developing theoretical methodology.^[46,47] Therefore, a demonstration and analysis of its intrinsic dynamic properties might also contribute to multidisciplinary advancements.

In this paper, we quantify the ZPV-induced fluctuation in geometric parameters, electronic-structure features and interaction energy of the H-bond in a ground-state water dimer through mode-specific ab initio molecular dynamics (AIMD) simulations, accompanied by electronic-structure analysis and EDA. Within the H-bond length (R_O…H) fluctuation up to 0.54 Å, two low-frequency intermolecular modes dominate its long-time trend with a spatial contribution about 0.30 Å and a temporal period about 220 fs. Two O–H covalent bond stretching modes are responsible for short-time R_O…H fluctuation peaks, giving minor contribution about 0.23 Å and shorter periods within 10 fs. The former intermolecular modes are responsible for the time-averaged R_O…H elongation about 0.03 Å that converges over 1 picosecond, while the latter ones are for the time-averaged elongation of that O–H covalent bond length (R_O–H) about 0.02 Å, converging within 100 fs. Remaining contribution from ∠ O…O–H variance is hard to assign to specific modes and generally upshifts upper bound of instantaneous fluctuation by about 0.10 Å. The ZPV-induced electronic-structure fluctuation in a water dimer crosses the border dividing into partially covalent and noncovalent H-boning for static models, thus being called an extension of fundamental understandings towards intrinsic and more practical dynamic scenario. Similar extension is further emphasized by EDA, whose temporal evolution conforms with a superficially random H-bonding energy around 2.63 kcal/mol to 6.04 kcal/mol and reveals the spontaneous alternation of the underlying bonding characteristics. Furthermore, spatial projection of EDA claims that the ground-state dynamics of H-bonding energy is beyond the capability of quasi-static dragging models. Distribution of total interaction energy to R_O…H is even qualitatively distinct from the conventional U-shape profile, due to intrinsic fluctuation under ZPV involving all geometric freedoms instead of only artificially chosen ones. The significance of intermolecular modes is further supported by both electronic-structure and EDA analyses, with intermolecular modes driving other modes in real space to fulfil the fluctuation of electronic-structure functions and manipulate the general trend of relative weights from different physical effects in H-bonding energy.

Calculations and simulations are carried out by adopting the second-order Møller–Plesset perturbation theory (MP2)^[48,49] with the augmented correlation-consistent polarized valence quadruple-ζ aug-cc-pVQZ basis sets. Equilibrium structure of water dimer is fully optimized and vibration frequencies are then calculated. The ZPV energy of each modeis converted from the obtained harmonic frequencies with the ZPV energy scaling factor of 0.9959 and fundamental frequency for referencing adopts a scaling factor of 0.9601 for modes higher than 1000 cm⁻¹ and 1.0698 for modes lower.^[50] Simulations are initialized on the equilibrium water dimer with all (specific) normal modes excited by ZPV energy for all-ZPV (single- and multi-mode) cases, replacing the thermal sampling with initial vibrational energy, and evolve under NVE conservation and Born–Oppenheimer molecular dynamics. Simulation timescale is around 350 fs, which could cover one complete vibration period of each mode and ensure that the effect of ZPVE flowing among different modes on the exploring of H-bonding nature in ground-state water dimer under the ZPV is slight. Simulation adopts Hessian based predictor-corrector method and corresponding step size relates to displacement under mass-weighted coordinates.^[51,52] Step size for all-ZPV case is 0.05 amu^{1 / 2} ⋅ Bohr, which gives an averaged temporal step size around 0.25 fs. Equal and smaller displacement steps with lower bound of 0.01 amu^{1 / 2} ⋅ Bohr are used for single-mode case. Fourier fitting curves of R_O…H and R_O–H are conducted to implement the extrapolation for time-averaging effect under modulated initial vibration phase. About 300 snapshots are almost evenly taken from the simulation trajectory to calculate and analyze the electronic-structure and conduct the EDA. Electronic-structure topological analysis in the quantum theory of atoms in molecules (QTAIM)^[26] is conducted to obtain the line critical point (LCP). The electron density at LCP (ρ_C) and the total electron energy density at LCP (H_C)^[31] are calculated. Symmetry-adapted perturbation theory (SAPT)^[53] under the second-order truncation with charge transfer analysis^[54] (SAPT2-ct) adopting the aug-cc-pVQZ basis sets is used to decompose the H-bonding energy, whose theoretical frame and interaction energy evaluation are both close to MP2.^[55] An MP2 molecular dynamics simulation is performed by using the smaller triple-ζaug-cc-pVTZ basis sets for conducting the estimation under longer simulation time. Optimization, simulations and electronic structure calculations are performed in Gaussain 09 D.01 package.^[56] The SAPT calculations are conducted in PSI4 package.^[57] The AIM electronic-structure analyses are accomplished with the help of Multiwfn code.^[58]

Water dimer is fully optimized and its 12 normal modes are depicted in the left panel of Fig. 1. Both calculated fundamental frequencies and approximate temporal periods reproduced from dynamic simulations excited by corresponding ZPV energies (ZPVE) are marked. While the frequency values are more for the parallelism with previously reported calculations, the reproduced periods are in better consistence with experimental observations^[59,60] thus are used to number modes from slower to faster ones. The right panel shows the plots of the dynamic variance of R_O…H when all modes are vibrating with ZPV energy on a temporal scale long enough to cover complete periods of all vibrations. Contributions from each mode when the mode is solely vibrating for about one period are also depicted. The dimer scheme decomposes R_O…H into three geometric determinants as the intermolecular O…O distance, R_O…O, the deviation angle of the H-bonded O–H covalent bond from the O…O axis, ∠O⋯O–H, and the length of that covalent bond, R_O–H. Colour codes of red, green, and blue are respectively assigned. Vibration modes seeming to dominate three determinants are shaded by the same colour code, which are respectively in low, middle and high frequency ranges. Optimized equilibrium geometry (Table S1 in Supporting information) is in good consistence with high-precision ab initio predictions^[22,23] while the R_O…O difference is about 0.007Å.

In general, the zero-point vibration of all modes tends to dynamically lengthen the H-bond. The fluctuation is in the shown temporal range from 1.752 Å to 2.293 Å, giving that a 0.541 Å amplitude from −0.190 Å to 0.351 Å is compared with the equilibrium length of 1.942 Å. A preliminary comparison with single-mode simulations could roughly attribute this fluctuation to three groups of mode contributions. Specifically, two low-frequency intermolecular modes as modes 2 and 4 might dominate the long-time trend. The red horizontal and vertical bar in the right panel of Fig. 1 approximately mark their maximal temporal period about 220 fs and the remarkable spatial bias tends towards elongation. And the trend of R_O…H curve by exciting all normal modes under ZPV energy coincides with that by exciting the mode 4. Middle-frequency modes with maximal period around 100 fs may also assist this elongation with a smaller amplitude as shown by green bars. Modes 5, 6 and 8 seem to contribute more in this group than the remaining mode 7. Two high-frequency covalent bond stretching modes on the proton-donor molecule, i.e., modes 9 and 11, shape the short-time fluctuation peaks. Their periods are within 10 fs but might leads the total bias to shorten the H-bond that is opposite to previous two groups, as shown by blue bars. Further rationalization of different spatial effects of three groups could be achieved by relating them to three geometric determinants of the H-bond. Intuitively shown by vibration schemes in left panels, modes 2 and 4 may dominate the variance of R_O…O, while modes 5, 6, and 8 may have remarkable influence on ∠O⋯O–H. Both the anharmonic lengthening of R_O…O and deviation of ∠O⋯O–H would elongate R_O…H, while modes 9 and 11 are likely to be responsible for anharmonically lengthening the R_O–H rather than shortening the H-bond. Schemes of these modes in Fig. 1 are temperately framed by colour codes. However, further quantification considering the blending effect of modes within each group and the time-averaging effect is urged to possess more practical significance.

Figure 2 demonstrates the fluctuations of three geometric determinants and their net influences on R_O…H, with comparisons between the all-ZPV cases and the blended contributions from corresponding modes in each group. A zoom-in strategy is presented by in sequence focusing on R_O…O, ∠O⋯O–H and R_O–H variances, whose fluctuation time scales are from the longest to the shortest. The net R_O…H influence of each determinant is derived by variable separation, following the simple trigonometric equation using R_O…O, R_O–H, and ∠O⋯O–H to construct R_O…H (Eq. (S1) in Supporting information). When deriving the net R_O…H variance due to the R_O…O fluctuation, ∠O⋯O–H and R_O–H are frozen at their equilibrium values. For net influence from ∠O⋯O–H variance, only R_O–H is frozen but the dynamic contribution from R_O…O in the former step is subtracted out to leave only the angle effect. For R_O–H, both net influences in former two steps are subtracted out.

It is verified that the group of modes 2 and 4 indeed dominates the total fluctuation of R_O…O, whose variance is consistent with the all-ZPV trend. The geometric projection conveys this domination to the long-time rising and the H-bond length decreasing, since the derived line penetrates the grey background curve depicting all-ZPV R_O…H. However, figure 2(b) illustrates that modes 5, 6, 7, and 8 cannot qualitatively reproduce the all-ZPV ∠O⋯O–H fluctuation. The inconsistency is remarkable no matter whether angle changes of four modes are simply summed up when they are solely vibrating or taken from another simulation with four modes simultaneously excited, similar failure is also projected to their net influence on R_O…H. In contrast, all-ZPV projection still reasonably directs the fluctuation trend of the remaining R_O…H after subtracting the long-time components from R_O…O variance, which verifies the validity of the variable separation strategy. Figure 2(c) is for R_O–H and verifies that the blending of modes 9 and 11 govern the fluctuation of this determinant as well as the corresponding short-time R_O…H peaks. Quantitatively, the net effect on variance of H-bond length from ΔR_O…O implemented by modes 2 and 4 is the largest, ranging from −0.114 Å to 0.190 Å and giving an obvious bias towards H-bond elongation. The variance of H-bond length from R_O–H by modes 9 and 11 ranges from −0.132 Å to 0.102 Å, inclining to H-bond shortening. The net effect from ∠O⋯O–H is difficult to assign to a few middle-frequency modes but indeed assists the H-bond elongation as preliminarily analyzed above, whose all-ZPV effect is from −0.005 Å to 0.096 Å with the smallest amplitude among all three determinants.

As a characteristic geometric parameter for water dimer, R_O…O in Fig. 2(a) fluctuates in the range of 2.825 Å–3.115 Å with a amplitude about 0.29 Å in all-ZPV vibration, which is in good agreement with the previous AIMD reported value 2.863 Å–3.102 Å obtained by using the same method but smaller basis sets.^[12] If the report focuses on the dissociation process of a dimer cation, no further mode-specific information about ground-state ZPVs can be revealed. Actually, the basic determination of low-frequency intermolecular vibrations in water dimer is still in process with recent experimental development.^[19,59,60] Vibration eigenvectors under the harmonic approximation as depicted in Fig. 2(a) are usually used for intuitively illustrating the real-space nuclear motion corresponding to observed bands, where modes 2 and 4 are also respectively referred to as the proton-acceptor wagging mode and the O–O stretching mode.^[48,49]

In Fig. 2(a), the blending effect of modes 2 and 4 are represented by simply summing up two solely vibrating cases. Such a fully-decoupled procedure successfully reproduces the total R_O…O thus suggesting that there is no obvious coupling between two modes in the all-ZPV simulation. This fortunately avoids the well-known deficiency of quasi-classical simulations in studying the ZPV effects that the coupling between ZPVs, namely, the energy leak from one ZPV mode to other ones might happen but should be conceptually forbidden due to quantum principles.^[61] Seeming to have no remarkable physical error, in this work, we could keep faithful to the outcome of classical nuclear motion on the quantum potential-energy-surface, which could better connect with upper-building studies focusing on particular differences between quasi-classical and fully quantum simulations (i.e., NQE),^[7,9,10] instead of taking forcible adjustments on kinetic parameters to imitate quantum nuclear behaviour, which measure has been stated to be inappropriate for the evolution of multiple ZPVs.^[61] For detailed verification, a simulation simultaneously exciting mode 3 and 4 but no other modes gives a delayed R_O…O curve compared with all-ZPV result and deviates especially after 200 fs (Fig. S1 in Supporting information), maybe suffering the artificially enhanced ZPV coupling after the unrealistic isolation of these two modes from the all-ZPV environment. For R_O–H, two kinds of curves from modes 9 and 11 almost coincide with each other and both qualitatively reproduce the total results (Fig. S1 in Supporting information) thus in neither choice, the physical insight or result validity discussed above would be violated, and only the collective curve from simultaneous excitation is plotted in Fig. 2 for clarity.

The failure in reproducing ∠O⋯O–H using both two blend measures seems more complicated. After quantitatively checking angle variances from each of all 12 ZPVs (Fig. S2), the fluctuation of this geometric determinant is found to be delocalized with respect to ZPV modes, where only four high-frequency O–H stretching modes actually have negligible effects. Simultaneously exciting eight lower modes indeed better reproduces the total ∠O⋯O–H trend (Fig. S3) but is incapable of providing further significant understandings. In detail, ∠O⋯O–H changing has two directions, namely, the in-plane direction (the plane is determined by two oxygen and the H-bonded hydrogen atoms) performed by modes 2, 4, 5, 7, and 8, and the out-of-plane direction by modes 1, 3, and 6, as shown in Fig. 1. Consequently, non-collinear cooperation among these modes makes it hard to intuitively assign the scalar variance of ∠O⋯O–H to individual ZPV modes. Meanwhile, this non-collinearity irrationalizes the decoupled sum measure then invalidate the comparison between two blending curves for verifying the possible ZPV coupling as discussed above. Considering less reliability in technical and possibly theoretical aspect, no additional mode-specific analysis is conducted on this angle determinant in this work. Studies exceeding the scope of this research reveal that the angle variance is closely related to the vibrational-rotational-tunnelling in water dimer, giving isomers far from the ground-state structure focused here, and the NQE needs to be taken into considertion.^[19,20]

After verifying the instantaneous domination of the above mode groups on R_O…O and R_O–H variance, we further investigate the more practical cumulative time averaging effects, where different initial vibration phases of two modes in the same group would induce different convergence behaviours. Technically, single-mode curve is fitted by Fourier series then analytically extrapolated with manipulated initial phases, i.e., 0, π/2, π, and 3π/2. The extrapolation of two modes belonging to the same group are to sum up and to cumulatively average over simulation time. Finally, 16 phase combinations for two groups are drawn in Fig. 2(d), with curves respectively coloured in red and blue gradation. It is illustrated that modes 2 and 4 influence the fluctuation of cumulatively averaged R_O…O in an effective time window around 10 fs–1000 fs. Before this window, curves are almost straight lines due to the relatively long vibration period. Initial R_O…O values before 10 fs are remarkably modulated by different initial phases with an amplitude over 0.40 Å. After this window, results from different phase combinations almost coverage to the same value. Similar time window is much earlier and narrower for R_O–H due to the faster vibration of modes 9 and 11, which is beneath 100 fs and the lower bound is shorter than 1 fs with initial R_O–H modulation amplitude around 0.24 Å. Consequently, the averaged H-bond length shows separated temporal regions that are in sequence sensitive to accumulations of R_O–H and R_O…O variance as highlighted by coloured shades in the middle panel of Fig. 2(d). The all-ZPV R_O…H beneath 10 fs is almost a mirror image of R_O–H in the same range, with the geometrical out-of-phase relation discussed above, and that above 10 fs obviously takes the trend of R_O…O curves.

The characteristic R_O…O in Fig. 2(d) is converged around 2.94 Å, specifically, 2.937 Å by mode extrapolation and 2.941 Å by all-ZPV estimation. Estimation involves another more efficient AIMD simulation adopting smaller basis sets (Fig. S4). This value is in good consistence with previous reports predicting the equilibrium distance around 2.91 Å^[22,23] and an anharmonic dynamic effect around +0.03 Å estimated with classical treatment of nuclear motions.^[21] The experimental observation^[21] such as 2.976 Å is alsomarked in the figure by using orange dotted line, while the remaining 0.03 Å difference could be attributed to secondary effect from NQE dynamic correction on the lightest hydrogen nuclei.^[9] It is worth noting that this dotted line is buried in the early-time brush of diverse red curves approximately before 100 fs–200 fs, which covers R_O…O around 2.8 Å–3.1 Å. The contact of specific fluctuating curves with the dotted line could even be found around 300 fs. Therefore, a long enough time scale is suggested with the quasi-classical frame to obtain reliable averaged R_O…O results, not only for theoretical simulations but also for experiments that may collect data after specific excitations through ultrafast techniques like pump-probe strategy.^[24,25] The converged value for R_O–H gives good agreement between the mode extrapolation (0.984 Å) and all-ZPV simulation (0.980 Å). However, the all-ZPV R_O…H curve does not show a convergence trend till the end of demonstrated time scale in Fig. 2(d). In fact, it would rise higher in longer simulation time due to the exchange of donor and acceptor role of two molecules through relative rotation (Fig. S5). Similar exchange has been reported in previous quasi-classical simulation^[4] and further physical description of the exchange process is beyond the scope of this work for involving the vibrational-rotational-tunnelling and QNE^[19,20] as mentioned above for the angle determinant.

As the origin that in-principle determines ground-state properties of a system, electronic structure has been long expected to provide an insight into the geometric appearances about the H-bonding. For an equilibrium system, criteria indicating the covalent-like characteristic and the strength of H-bonding have been established with the help of the quantum theory of atoms in molecules (QTAIM).^[26] At the bond critical points (BCP) that are uniquely determined by QTAIM topological analysis in the H-bonding region, sign for Laplacian of electron density (∇²ρ_C) and total electron energy density (H_C) are widely used to classify the bonding nature into covalent (∇²ρ_C < 0 and H_C < 0), partially covalent (∇²ρ_C > 0 and H_C < 0) and noncovalent (∇²ρ_C > 0 and H_C > 0), while the value of electron density (ρ_C) is usually in proportion to bonding strength.^[27,30,31] However, the ZPV dynamic fluctuation would diversify the static equilibrium electronic structure into a vibrationally evolving entity, thus presenting a diffuse distribution in real-space for critical points and in value for all corresponding electronic-structure functions. With more than 300 snapshots nearly-evenly taken from the all-ZPV simulation within 350 fs, figure 3 demonstrates such a dynamic nature for ground-state water dimer and the underlying relation to the variance of three geometric determinants as well as corresponding ZPV modes. Methodologically, QTAIM analyses could be used in a non-equilibrium system to reveal the dependence of interaction on internuclear distance^[26,62,63] and the effect of nuclear vibration has also been discussed when developing the connotation of the theory.^[64] Including non-equilibrium cases, the more generalized terminology of atomic interaction line (AIL) is used to denote the maximum-electron-density path connecting two interacting atoms which is originally termed the bond path in the equilibrium case,^[26,28] and corresponding dynamic counterparts of BCP would refer to critical points in AILs (LCPs). In this work, ∇²ρ_C, H_C, and ρ_C below are used for investigating the properties of all LCPs including the equilibrium BCP (denoted as BCP_eq).

It is shown in the upper panel that the ZPV fluctuation qualitatively changes the sign of H_C, while ∇²ρ_C > 0 is preserved in the whole process. Thus, the H-bond in the ground-state water dimer dynamically crosses territories of partially covalent and noncovalent H-bonding, judged by experiences from equilibrium systems. On the other hand, the electron density distribution plotted in the lower panel remarkably varies from 0.011 a.u. to 0.037 a.u., giving a fluctuation amplitude around 50% compared with the equilibrium value of 0.025 a.u, which implies a possibly remarkable dynamic variance of the interaction strength. The real-space distribution of H-bond LCPs gradually forms a pancake-like cluster with a maximum radius around 0.30 Å and maximum thickness around 0.21 Å as shown in the inset of Fig. 3. Unlike the unitary critical point in equilibrium H-bond model (BCP_Eq in Fig. 3), this extends the basic picture when understanding the ground-state water dimer from a real-space electronic-structure angle.

For further clarifying the inclination of electronic-structure fluctuation and the mode contribution, distribution regions for criterion functions in Fig. 3 are divided and all data points are coloured. With positions of equilibrium H_C and ρ_C being denoted by crossings of dashed lines, the regions with more negative H_C and ρ_C larger meanwhile shorter R_O…H are shaded in light yellow colour, indicating more covalent-like and stronger bonding characteristics. Regions with opposite conditions for less covalent and weakened ones are shaded in light grey colour. The colouring of data points is achieved by decomposing ΔR_O…H of each snapshot into variances of three geometric determinants. Then, contributions from determinant variances could be conveyed to corresponding ZPV mode groups with the above clarified geometric correspondence. Specifically, variable-separation strategy used in Fig. 2 is adopted and absolute values of net effects from three determinants are divided by their sum to obtain the normalized contributions. The colour code used above is also inherited when variance of three determinants respectively gives dominant contributions, while the gradation among three vertex colours is for intermediate snapshots as shown by the colour triangle in Fig. 3.

Generally, the ZPV fluctuation drives the H-bond electronic-structure towards weakened and less covalent-like bonding characteristics, given more H_C and ρ_C data points belonging to grey regions. The points in the colour themes of red, blue and green colour mainly in sequence concentrate in grey regions, yellow regions and around equilibrium values with an inclination to grey ones. Thus, when the other two determinants are coincidentally around equilibrium values, R_O…O variance implemented by modes 2 and 4 dynamically drives the electronic-structure of H-bond towards noncovalent and weakened bonding characteristics, while that of R_O–H by modes 9 and 11 is towards the opposite direction, and effects of the angle variance is less remarkable or biased, being in consistence with the above geometric regulations. Colouring of real-space LCPs more intuitively explains such contributions. For instance, while the oxygen atoms of the left molecule in the inserted scheme are fixed for all snapshots, the red LCPs are mainly present on the right of BCP_Eq, corresponding to the elongated H-bond lengths by intermolecular departure from modes 2 and 4. Thus, the weaker electronic structure contact between two water molecules is reasonably manifested, as revealed by H_C and ρ_C values.

More interesting cases are given by the cooperation between different mode groups. There are purple points at the left end of H_C and ρ_C distribution, which indicates that the superposition between R_O…O shortening and R_O–H elongation, achieved by the constructive cooperation between corresponding two mode groups, is responsible for the most covalent-like and strongest bonding characteristics. Besides, the destructive cancelation between these two groups also presents purple points in the middle of both patterns, reasonably giving moderate fluctuations around equilibrium values. At the opposite ends, olive points indicate that the weakest H-bonding is more attributed to R_O…O elongation and ∠O⋯O–H deviation. The LCPs real-space locations corresponding to these three cases are respectively purple ones on the left of BCP_Eq, purples ones around the centre of LCPs cluster and olive ones in the top-right and bottom-right parts of the cluster. It is noteworthy that R_O…O variance participates in all three cases and penetrate the whole H_C and ρ_C distribution pattern. Therefore, though modes 2 and 4 cannot produce the whole LCPs cluster in real space, they serve as the driving modes carrying other modes to fulfil the ZPV fluctuation on CP criterion functions. Such a role matches their domination in the general trend of geometric fluctuations revealed above.

It should be mentioned that the of the equilibrium water dimer in Fig. 3 is small but negative (about −0.001 a.u.), while previous reports gave positive value (about 0.002 a.u.).^[31,65,66] Though in either case its magnitude is small and satisfies the H-bonding criterion of |H_c|,^[31,67] the sign-changing of implies the inconsistence in classifying the H-bonding into partially-covalent and noncovalent interaction strictly according to criteria introduced above. This discrepancy might be mainly attributed to different sizes of basis sets in producing the electron density. Equilibrium becomes negative when large basis sets like aug-cc-pVQZ in this work are used, and becomes positive based on smaller ones like 6-311++G(d,p) that are more widely adopted in previous reports. Besides, more rigorous calculation methods with larger basis sets similar to those used for electron-density benchmarking^[68] also give negative values when performing on a more precise equilibrium structure (Table S2). Specific discussion of equilibrium result and related technical details would not influence above conclusions, which focus on the remarkable electronic-structure fluctuation induced by the ZPV and extension of the H-bond understanding beyond static view.

The significance of interaction energy for characterizing bonding is indisputable. Theoretical decomposition of total interaction energy into physically meaningful components permits further revealing the bonding nature of water dimer.^[35–43] Considering ZPV motion is implemented by vibrations with distinct temporal periods, resulting energy outcomes may convey characteristics of different modes in temporal fluctuation. For projection on the particular spatial variable, the interaction energy may differ from conventional impressions of U-shape profiles around the equilibrium geometry, considering that all geometric freedoms now vary along the eigenvectors of 12 normal modes instead of being mostly frozen in artificial manipulation. In Fig. 4, interaction energy and decomposed terms of above snapshots are analyzed. Technically, the EDA adopts the symmetry-adapted perturbation theory (SAPT),^[53] giving E_int = E_elst + E_exch + E_ind + E_disp, where E_int, E_elst, E_exch, E_ind, and E_disp are respectively the term for total interaction, electrostatic, first-order exchange, induction and dispersion energy. For the intuitive explanation and visualization convenience, first two terms are summed up to represent steric effect, i.e., E_ster = E_elst + E_exch, and similar measure has been taken in previous SAPT analysis on water dimer.^[43] Charge transfer energy (CT) is also presented which is covered in the induction energy. The difference between induction energy values calculated in the dimer basis and the monomer basis can therefore be reasonably taken as a measure of the CT.^[54] With clear physical backgrounds,^[53] the SAPT results has been used to construct force fields describing the water dimer potential.^[20,69] The adopted theory level (see Methods for details) shows good consistence in both form^[55] and value with the above ab initio method, giving equilibrium interaction energy values: about 4.93 kcal/mol obtained by using SAPT and 4.90 kcal/mol in above calculation after conducting the basis sets superposition correction (BSSE).^[70]

The temporal evolutions of interaction energy and decomposed terms are demonstrated in the left panel of Fig. 4, with the upper part showing energy values and the lower part referring to percentages of E_exch, E_ind, and E_disp in E_int. Generally, no obvious temporal regularity is found in evolution of E_int, which seems to randomly fluctuate and return to its equilibrium value. However, decomposed energy terms present clear periodic patterns which prove slight influence of energy flowing among vibrational modes. Their general rise and fall show approximately two periods, whose behaviour reminds us about the long-time trend of ΔR_O…H within the same timescale shown in Fig. 2(a), thus might be attributed to the change of intermolecular distance governed by modes 2 and 4. Meanwhile, all decomposed curves carry synchronized shock peaks with periods stabilized around 10 fs. Apparently, they are energy responses to the ZPV stretching of H-bonded O–H covalent bond implemented by modes 9 and 11, as analyzed in Fig. 2(c). Therefore, the mode dependence of ZPV effect in H-bonding manifests in its energy infrastructure. Besides, there is irregular turbulence seeming to emerge every few shock peaks, which interferes the peak-valley pattern of decomposed terms then induces erratic fusing of E_int peaks. This relatively dim feature might relate to the angle variance, where related modes have diffuse fluctuation periods being all several times longer than stretching modes.

More significantly, the percentage evolution reveals the spontaneous periodic change of H-bonding nature in ground-state water dimer beneath the superficially random E_int. This insight is quantitatively highlighted by distinct term weights marked at three different moments but all correspond to the same equilibrium E_int, as marked by dashed and dotted lines in left panels. Specially, E_ind and E_disp dominate the stabilizing effect at the equilibrium moment with weights respectively around 50.57% and 44.07%, then could be further enhanced to about 70.20% and 53.19% due to R_O–H stretching when two molecules are still near around 30 fs, while that minor effect from E_ster is flipped to destabilization from 5.36% to −23.39%. Finally, weights of all three terms are homogenized to about 33.26%, after the long-time evolution around 130 fs that corresponds to the intermolecular departure process. Clearly, such a dynamic picture extends the static view where only the equilibrium result is used to represent the H-bond characteristic. The detailed explanation about the temporal features involves the correspondence between ZPV induced geometric variance and energy responses in different terms, which could be more naturally discussed in spatial projection of EDA below. The EDA result of equilibrium state is in good consistence with the previous report.^[43]

Energy values with respect to R_O…H are plotted in the right panel of Fig. 4, where referencing profiles obtained by pushing and pulling two fixed monomers with fixed H-bond angle are also drawn in dashed curves. Remarkably, ZPV E_int distribution violates the U-shape profile with an opposite curvature in the most part and reverses back in a bounce at the tail for longer R_O…H. Quantitatively, its fluctuation range is unexpectedly large compared with results in artificial manipulation and generally rises. The ZPV E_int now fluctuates around −2.63 kcal/mol to −6.04 kcal/mol, presenting over 370% of the amplitude of the U-shape profile in dragging for the same ΔR_O…H effect. The latter range is around −4.01 kcal/mol to −4.93 kcal/mol, similar to the range reported in Refs. [42,43]. Moreover, the ZPV E_int in the far-left region is even lower than the equilibrium E_int, whose value is conventionally the minimum point to anchor dragging curves using fixed monomers. Such numerical anomalies suggest that the ZPV motion allows not only widened dynamic redundance but also enhanced instantaneous strength for intermolecular interaction, possibly at the cost of monomeric deformation that is difficult to capture by static models. These remarkable distinctions may dramatically extend our understanding of the H-bonding nature in ground-state water dimer and are clearly explained by combining EDA and geometry analysis below.

The E_int divergence from dragging curve could be mainly attributed to behaviours of induction and steric terms. In phenomenological assignment, the reversed E_int curvature seems to be due to the collective effect of much steepened E_ind when R_O…H is shorter than the equilibrium distance and the remarkably shifted E_ster that almost fills up the original concavity of the corresponding curve. Physically, induction term represents the interaction between induced multipole moments of one monomer and the permanent multipole moments of the other monomer,^[38] thus further E_ind dive implies that the extra intramolecular deformation might effectively enhance related monomeric electronic-structure properties. On the other hand, the general rise of total interaction energy is apparently dominated by the E_ster shift, as well as its tail bounce that is almost reproduced by E_ster pattern. Intuitively, the steric term is sensitive to the pair-wise interatomic distance between two monomers and their relative orientation. Its components of exchange repulsion grow fast when intermolecular electronic-structure overlap is stronger in closer distance, while the other component of electrostatic attraction increases slower in the approaching to the charged entities and weakens if they deviate from the electrostatically preferred relative orientation. Consequently, the rise and bounce of E_ster might be related to occasions where ZPV deformation effectively draws near molecular constituents of two molecules under the same R_O…H and disturbs their relative orientation from E_elst preferred one. Thereafter, that the interaction energy reduces to what is corresponding to the geometric difference is further explained by the ZPV in specific regions and whether they could implement the above inferred physical effects.

The extreme case in E_ind deepening is an ideal staring point towards the complete dynamic picture. Geometrically, to reach the far-left end for the shortest R_O…H, the ZPV fluctuation requires R_O–H to be maximally stretched on the basis of linearized ∠O⋯O–H and exhausted R_O…O approaching redundance. Compared with quasi-static model under the same R_O…H, this stretch effectively further pushes 0.13 Å away from the molecular continum of the donor molecule behind that denoted hydrogen atom, whose value can be obtained from Fig. 2(c) considering the nearly linear ∠O⋯O–H. Consequently, the lengthened O–H bond in donor molecule loosens its molecular continents and presents less compact electronic structure, especially along the H-bonding direction, thus might be easier to polarize by the acceptor molecule and produce larger induced multipole moments to enhance E_ind. Extending to dynamic scenario, ideally reversed effect with weakened E_ind arises when R_O…H shrinks compared with the equilibrium length, providing ∠O⋯O–H and R_O…O hardly vary. Such occasions are temporally guaranteed by the relatively fast vibration of stretching modes 9 and 11, which allows R_O–H to stretch and shrink several times before those slower modes could obviously modulate the other two geometric determinates. The R_O…H shortening effect solely from O–H bond stretching would be returned, thus snapshots with such effects might be 0.13 Å away from the left end. Delightfully, E_ind points in Fig. 4 indeed relatively rise until R_O…H is beyond a specific scope close to this evaluated value, thus the E_ind steepening across rising-to-lowering points in this local region seems to be explained. The detailed temporal evolution (Fig. S6) verifies the R_O–H dependant polarizability of the donor molecule, which is precisely out-of-phase synchronized with E_ind shock peaks and fluctuates around the equilibrium value, especially for its diagonal components corresponding to the H-bonding direction. The permanent dipole moment of the donor molecule, which can offer an approach to elucidating fundamental questions of H-bonds and whose fluctuating values are in a similar range to that in previous study^[71] and lower than that in the other study,^[72] is also checked but found to be more dependent on the variance period of ∠O⋯O–H instead of R_O–H, thus has slight influence in this linear-angle preferring region but could support the electrostatic explanation for Ester below. And accordingly the Raman spectra of frequency modes (9 and 11) are marked as shown in Fig. S7, which is similar to previous ones. Besides, fluctuations of polarizability and permanent dipole moments of the acceptor molecule lose the reasonable temporal and phase synchronicity with E_ind peaks, hence are not major physical origins here.

Similar explanation strategy is applied to E_ster. Before analyzing the general rise and bounce tail, its minor feature in the far-left region discussed above could be apparently rationalized. Namely, relative further distance between molecular continents weakens the component of E_exch and E_elst (Fig. S8) while E_exch decreases more remarkably in such a close range. Consequently, their sum E_ster is lowered. This explains that only in this region of Fig. 4, the general rise of E_ster is obviously suppressed and some scatter points are even lower than the dragging curve, assisting the E_int lowering as a minor effect besides E_int drops considerably. In contrast, E_ster bounce around the opposite end with R_O…H approximately longer than 2.2 Å corresponds to completely reversed occasions, which gives real-space olive LCPs in the top-right and bottom-right parts of the cluster in Fig. 3. Therein, the H-donating O–H covalent bond squishes moreover has to be swung away, given that the exhausted R_O…O lengthening only pushes R_O…H to around 2.13 Å and that plus R_O–H shortening is still within 2.23 Å as shown in Fig. 2. Effectively, with respect to the same R_O…H, the molecular continum of donor molecule is now drawn near meanwhile swung forward, thus is closer to the acceptor molecule. Consequently, E_exch repulsion is simply enhanced. However, the weaker enhancement of E_elst attraction is further compensated for by the disorientation of permeant multipole moments from their electrostatically preferred linear ∠O⋯O–H, whose angle preference has been previously reported.^[43] Without effective balancing from E_elst, the E_ster in this region is almost dominated by E_exch thus remarkably rises and finally gives a bounce towards the longest R_O…H (Fig. S8), with the fiercer R_O–H shrink and ∠O⋯O–H deviation for further longer R_O…H. The E_ster turning neck around 2.2 Å in Fig. 4 is close to the above-mentioned upper bound of 2.23 Å, beyond which the angle swing is indispensable and the balancing from ΔE_elst to ΔE_exch starts to obviously fail.

Far from extreme cases at both ends, the middle part of E_ster pattern corresponds to the mixture of different geometric variances and mainly implements the general rise. In this region, E_exch snapshots just diffuse on both sides of the referencing curve due to ZPV fluctuations having no longer obvious bias. However, E_elst still one-side rises (Fig. S8) and drives E_ster, since the ZPV disorientation effect is active until the far-left region and always inclines to electrostatically less preferred nonlinear angle. So far, main features of ZPV extension to H-bonding energy in a spatially projected picture have been well explained. The E_disp variance is also in consistent with the above analyses, which slightly lowers due to the shorter effective distance between monomers at left end then slightly rises due to reversed effect at the opposite end. The charge transfer term is well conserved to the dragging case and has small weight in the induction term. Coloring assignment of E_int snapshots according to ΔR_O…H contributions from three geometric determinants, and the ∠O⋯O–H similar to Fig. 3 is also conducted, respectively (Figs. S9 and S10) In each small angle region (about 5° as shown in the right figures, when R_O…H is shorter than the equilibrium distance, the E_int curves all remarkably rise and when it is larger than the equilibrium distance, the curve level is off. The tail bounce of E_int curve is due to the increase of R_O…H caused by angle, which is also shown in Fig. S9. Being in consistence with electronic-structure analyses, assignment supports the throughout participation of R_O…O across the fluctuation process, the remarkable contribution from ∠O⋯O–H variance at the bounce tail and in the middle part, as well as the ΔR_O–H effect towards the shortest R_O…H.

After explaining the correspondence between geometry term and energy term, the details of temporal evolutions in left panels of Fig. 4 could now be easily understood. The narrowing of energy curves from the beginning moment to around 130 fs indicates the intermolecular departure by modes 2 and 4 and suppresses the shock-peak fluctuation from stretching modes 9 and 11, which is reasonable to consider that the energy responses of all intermolecular physical effects to monomeric fluctuation weaken when two molecules separate far. On the other hand, weight homogenization in this process is induced by different R_O…O sensitivities of these effects. Specifically, the induction decreases faster than dispersion though it initially contributes more to the stabilization, and the flipping of steric term is due to the slower response of stabilizing electrostatic attraction than destabilizing exchange repulsion towards longer R_O…O. Besides, the phase cooperation between E_ster and other terms switches from out-of-phase to in-phase around 70 fs within the first EDA period. This is also attributed to the dominating term in E_ster changing from E_exch to E_elst, since the former gives strengthened repulsion in R_O–H lengthening while the latter and other terms become more attractive. The pattern irregularity assigned to angle-varying modes could be further explained as the disorientation of multipole moments, which intuitively shallows E_ind valleys every a few stretching periods, and influencing patterns of E_ster as well as the total interaction energy as discussed above.

In this work, the intrinsic dynamic nature of the elementary H-bond in a ground-state water dimer is theoretically revealed within the quasi-classical frame, giving mode-specific explanation of ZPV induced geometric, electronic-structure and energy fluctuations. Two intermolecular modes generally weaken the H-bond and act as the driving modes in all ZPVs, which leads to long-time R_O…H elongation, real-space electronic-structure excursion towards non-covalent H-bond and H-bonding energy characteristic towards the same direction. Two stretching modes are for minor opposite effects such as short-time R_O…H shortening and electronic-structure features as well as EDA shock peaks towards stronger and covalent-like interaction. The overall dynamic nature extends conventional H-bond understandings based on equilibrium and quasi-static dragging models in both temporal aspect and spatial aspect, weaving the lining for the thorough comprehension and fine manipulation of H-bonded systems.