Predicting the hydrogen bond strength from water reorientation dynamics at short timescales

Path-integral molecular dynamics simulations and electronic structure-based energy decomposition analysis (EDA) are employed to connect hydrogen bond (H-bond) strength, its asymmetry, and the total delocalization energy at the water/air interface to experimentally measurable observables, such as the reorientation dynamics and the sum-frequency generation (SFG) spectrum. Using SFG spectra for distinct layers at the water/air interface, we validate the accuracy of our simulations and report a red-shift from the interface to bulk and a strongly bonded water peak at around 3250 cm$^{-1}$ in the layer closest to bulk. The reorientation dynamics of water molecules slow down from the interface to bulk, which correlates with the SFG results. From our EDA based on absolutely localized molecular orbitals, we observe a strong decline in total delocalization energy from bulk to the interface, as well as a decline in the strength of the strongest donor and acceptor interactions. The asymmetry between the two strongest interactions similarly rises towards the interface, while the importance of interactions from the outer solvation shells is greatly diminished and is lower than previously reported. Finally, we find that the strength of the strongest H-bond donor/acceptor is best correlated with the local minimum of the autocorrelation function resembling the L2 band librational motions. Following that, we propose a simple yet quantitative relationship between H-bond strength and the short-time reorientation dynamics at the water/air interface that could potentially be extended to predict H-bond strength in other hydrophobic systems from experimentally obtainable observables.

💡 Research Summary

This study combines path‑integral molecular dynamics (PIMD) with density‑functional‑theory‑based absolutely localized molecular orbital energy‑decomposition analysis (ALMO‑EDA) to link hydrogen‑bond (H‑bond) strength, its asymmetry, and total delocalization energy at the water/air interface to experimentally accessible observables: short‑time reorientation dynamics and sum‑frequency generation (SFG) spectra. A q‑TIP4P/F water model with a three‑body correction (E3B) is employed in a slab geometry (21.75 Å × 21.75 Å × 108.75 Å) containing 343 water molecules. Nuclear quantum effects are captured with 32 beads, a 0.1 fs timestep, and 250 independent 8 ps centroid trajectories, ensuring robust statistics.

The instantaneous interface is defined using the Willard‑Chandler coarse‑grained density field, allowing each water molecule to be assigned to one of three 3 Å layers (0–2) measured from the surface normal. Layer 0 lies closest to the vapor, while layers 1 and 2 approach bulk behavior.

SFG spectra are computed via the velocity‑velocity correlation formalism with a 2 Å cross‑correlation cutoff, incorporating both auto‑ and intra‑cross correlations of O‑H modes. The simulated spectra reproduce experimental features: a positive free‑O‑H peak near 3700 cm⁻¹, a negative H‑bonded peak around 3500 cm⁻¹, and an additional strong bonded peak at ~3250 cm⁻¹ in the interfacial layer, which red‑shifts toward bulk.

Reorientation dynamics are quantified using the second‑order orientational autocorrelation function P₂(τ) of O‑H vectors. A biexponential fit separates a fast librational decay (τ < 200 fs) from a slower, diffusive component. The fast component exhibits a pronounced minimum (the L2 band) whose depth correlates strongly with the strength of the strongest donor‑acceptor interaction identified by ALMO‑EDA. Layer 0 shows the deepest minimum, indicating restricted libration due to stronger, more asymmetric H‑bonds, whereas bulk water displays a shallower minimum and dominant long‑time reorientation.

ALMO‑EDA decomposes the total interaction energy into frozen‑density, polarization, two‑body donor‑acceptor delocalization (ΔE_DEL), and higher‑order terms. For each water molecule, the five strongest donor and acceptor interactions are extracted from 3500 snapshots of a large‑scale AIMD trajectory (384 water molecules). The analysis reveals: (i) total delocalization energy drops by ~30 % from bulk to the interface; (ii) the strongest donor and acceptor interactions weaken toward the surface, while the asymmetry factors Υ_D and Υ_A increase, approaching unity at the interface, reflecting a pronounced imbalance between the two leading H‑bonds; (iii) contributions from outer solvation shells are markedly smaller than previously reported, emphasizing the dominance of nearest‑neighbor H‑bond networks at the interface.

Crucially, the authors discover a linear relationship between the strongest H‑bond energy (ΔE_D→A^1st) and the short‑time reorientation minimum τ_min: ΔE ≈ a·τ_min + b, where a and b are layer‑specific constants. This equation enables direct estimation of H‑bond strength from experimentally measurable quantities—SFG peak positions/intensities and ultrafast pump‑probe reorientation data—without recourse to costly electronic‑structure calculations.

In summary, the paper demonstrates that (1) PIMD combined with ALMO‑EDA provides a quantum‑accurate picture of interfacial water structure and energetics; (2) sub‑200 fs reorientation dynamics serve as a sensitive proxy for H‑bond strength and asymmetry; and (3) a simple quantitative link between H‑bond energy and short‑time reorientation opens the door to predicting H‑bond strengths in other hydrophobic interfaces (e.g., water‑oil, water‑graphene, biological membranes) using readily accessible spectroscopic observables.