Time-Dependent Friction Effects on Vibrational Infrared Frequencies and Line Shapes of Liquid Water.

From ab initio simulations of liquid water, the time-dependent friction functions and time-averaged nonlinear effective bond potentials for the OH stretch and HOH bend vibrations are extracted. The obtained friction exhibits not only adiabatic contributions at and below the vibrational time scales but also much slower nonadiabatic contributions, reflecting homogeneous and inhomogeneous line broadening mechanisms, respectively. Intermolecular interactions in liquid water soften both stretch and bend potentials compared to the gas phase, which by itself would lead to a red-shift of the corresponding vibrational bands. In contrast, nonadiabatic friction contributions cause a spectral blue shift. For the stretch mode, the potential effect dominates, and thus, a significant red shift when going from gas to the liquid phase results. For the bend mode, potential and nonadiabatic friction effects are of comparable magnitude, so that a slight blue shift results, in agreement with well-known but puzzling experimental findings. The observed line broadening is shown to be roughly equally caused by adiabatic and nonadiabatic friction contributions for both the stretch and bend modes in liquid water. Thus, the quantitative analysis of the time-dependent friction that acts on vibrational modes in liquids advances the understanding of infrared vibrational frequencies and line shapes.

Introduction

The OH stretch band in liquid water is significantly red-shifted and broadened compared to the gas-phase spectrum, while the HOH bend frequency is in fact slightly blue-shifted when going from gas to the liquid phase.[1] The broadening of the OH stretch band in liquid water is typically rationalized by a combination of homogeneous and inhomogeneous effects.[2,3] Inhomogeneous line broadening is associated with different hydrogen-bonding environments of individual OH bonds, which in the limit when the hydrogen-bonding pattern changes more slowly than the OH vibrational period and in the presence of nonlinearities in the OH bond potential produce vibrational frequencies that vary over time.[4−6] Homogeneous line broadening reflects the fast coupling of OH bonds to their neighboring water molecules, mostly via hydrogen bonding, which reduces the vibrational lifetime because the vibrational energy is quickly transported to neighboring molecules and thus dissipated into collective modes.[7] Indeed, the vibrational lifetime of the OH stretch is very short (of the order of 190 fs[8,9]) and thus only 19 times longer than the OH-stretch vibrational period itself (of the order of 10 fs). The experimentally observed red shift of the OH stretch band is usually rationalized by strong hydrogen bonding in liquid water, which extends and thereby softens the OH bond[2,3] (in fact, the relationship between the hydrogen-bond strength, the OH bond length and the red shift of the stretch band has been amply and partly controversially discussed in literature[5,10−12]). According to such reasoning, the rather small frequency shift of the water bending mode when going from gas to liquid water could be argued to imply that the bond angle potential is only weakly perturbed by the liquid water environment and thus that the coupling of bend vibrations to the hydration environment is weak. This interpretation is puzzling though, because the vibrational lifetime of the water bending mode in liquid water is rather short (around 170 fs[13,14]) and thus only 8.5 times longer than the vibrational period of 20 fs, an even smaller ratio than for the stretch mode. The short bend vibrational lifetime reflects quick energy dissipation into librational modes,[15−17] which in turn can be rationalized by efficient multiphonon energy relaxation based on the excitation of librational overtones in liquids.[18]

Time-dependent or, equivalently, frequency-dependent friction arises whenever the dynamics of a many-particle system is described in a low-dimensional reaction-coordinate space,[19−25] and its relevance for infrared (IR) spectra was clearly demonstrated in the past.[26−31] All friction contributions that decay faster or similarly as the vibrational period stem from adiabatic solvent degrees of freedom and account for dissipation into intra- and intermolecular degrees of freedom (including vibrational overtones);[32−35] these friction contributions dominate the vibrational energy relaxation and lead to homogeneous line broadening. Friction contributions that decay much slower than the vibrational time scale describe the slowly changing nonadiabatic hydration environment and in conjunction with nonlinear bond potentials induce inhomogeneous line broadening, as our results explicitly demonstrate. Of course, there is no clear-cut separation between adiabatic and nonadiabatic solvent relaxation modes,[36−39] prompting for a time-scale bridging framework to treat the dynamic coupling of molecular vibration modes and their environment. In fact, the frequency-dependent friction function, which appears in the generalized Langevin equation (GLE), is the appropriate framework to account for all these effects, with the only drawback that nuclear quantum effects can at the current level of the formalism not be included without making additional approximations. Only recently developed extraction methods that account for nonlinearities[40] allow obtaining these time-dependent friction functions from ab initio molecular dynamics (aiMD) simulations and with high enough accuracy; this we self-consistently demonstrate by deriving vibrational spectra from numerical simulations of the GLE that are virtually indistinguishable from the vibrational spectra directly obtained from aiMD simulation trajectories.

In this paper, we address the puzzle posed by the different line shifts of the water stretch and bend modes by analyzing the vibrational water dynamics in terms of the time-averaged nonlinear bond potentials (as a function of the bond length for the OH stretch and the bond angle for the HOH bend) and the corresponding time-dependent friction functions, which are extracted from extensive aiMD simulations for 256 H2O molecules. In particular, we show that the slight blue shift of the water bend mode when going from gas to the liquid phase is not caused by a stiffening of the bend potential, which would explain the blue shift,[1] but rather by the time dependence of the friction acting on bending vibrations. We find that the liquid environment in fact significantly softens the time-averaged bond potentials, and it does so quite similarly for the stretch and bend modes. Neglecting the frequency dependence of the friction, both stretch and bend bands would thus be expected to be red-shifted by comparable amounts when going from gas to the liquid phase, in stark contrast to the experimental finding.[1] It turns out that nonharmonic bond-potential effects are rather unimportant for the band position and thus cannot explain this puzzling finding. Likewise, frequency-independent friction shifts the bands insignificantly and only increases the line width, in agreement with expectations.[4] In contrast, the frequency dependence of the friction is crucial and not only leads, in conjunction with nonlinearities in the bond potentials, to inhomogeneous line broadening but also gives rise to pronounced blue shifts for both stretch and bend bands. The mechanism for this blue shift is very general,[26] as we analytically demonstrate. The compensation of the potential red shift and the friction blue shift is incomplete for the stretch band but almost perfect for the bend band, so the stretch band exhibits a significant net red shift from gas to liquid, while the bend band shows only a slight blue shift in both experiments and simulations. The absence of a significant frequency shift of the bend mode does by no means imply that bend vibrations couple less to their environment than stretch vibrations (as has been demonstrated previously[15−17]); rather, it is the subtle balance of the potential and friction contributions to the line shift, which both are caused by interactions with the liquid environment, that is different for the stretch and bend bands. We conclude that the coupling of water stretch and bend vibrations to other intra- and intermolecular degrees of freedom, as quantified by the time-averaged bond potentials and friction functions, is of similar strength, which explains their similar vibrational life times, although their frequency shifts are rather different, which we rationalize by a subtle difference of the compensatory potential and friction effects. The spectral blue shift due to frequency-dependent friction is a very general mechanism; it transpires that the concept of frequency-dependent friction is important for advancing the understanding of vibrational spectroscopy.

System, Spectra, and Model

We primarily analyze aiMD simulations of 256 H2O (and D2O for comparison) molecules in the liquid phase at 300 K that neglect nuclear quantum effects. Figure A compares the trajectories of the mean OH bond length of a single H2O molecule in liquid H2O (blue line) and in the gas phase (gray line), both at 300 K (see Methods for simulation details). The increase of the mean and the variance of the bond length in the liquid phase compared to the gas phase is clearly visible, which reflects the shift and softening of the OH bond potential due to hydrogen bonding in the liquid phase. The slow fluctuations of the oscillation amplitude reflect vibrational energy relaxations that occur over about 100 fs in the liquid phase (pure dephasing due to fluctuations of the vibrational frequency[41,42] is not easily visible in the time domain). Similarly, the bond-length velocity autocorrelation function (VACF) in Figure B demonstrates a significantly faster decay and thus a decreased vibrational lifetime in the liquid phase. Although the OH-stretch absorption spectrum is (apart from electronic and collective effects) straightforwardly related to the OH bond-length VACF via Fourier transformation, it turns out that a careful analysis of the molecular vibrations in terms of the GLE reveals interesting information on the mechanisms that determine the vibrational frequencies and line shapes.

Figure 1

(A) Trajectory of the OH bond length, averaged over both OH bonds in a single water molecule, from ab initio molecular dynamics (aiMD) simulations of one H2O in the gas phase (gray line) and for 256 H2O molecules in the liquid phase (blue line), both at 300 K. (B) Corresponding velocity autocorrelation functions.

Linear IR spectroscopy experiments measure the absorbed power of light at angular frequency ω = 2πf, which is proportional to the imaginary part of the dielectric susceptibility χ̃″(ω). Linear-response theory relates χ̃″(ω) to the total dipole-moment autocorrelation (see section I of the Supporting Information), allowing IR spectra to be calculated from equilibrium simulations.[45,46]Figure A compares the IR absorption spectrum from aiMD simulations of liquid H2O (gray solid line) and D2O (blue solid line) with corresponding experimental data (gray and blue broken lines, respectively). One discerns the stretch band (around 3300 cm–1 for H2O and 2400 cm–1 for D2O in the aiMD results) and the bend band (at 1650 cm–1 for H2O and 1200 cm–1 for D2O). The librational absorption band is produced by a large number of different intermolecular vibrational modes[47] that are dominated by rotational vibrations of water molecules in their hydrogen-bond environment (around 700 cm–1 for H2O and 550 cm–1 for D2O) and by translational vibrations of water molecules against each other around 200 cm–1 for both H2O and D2O. The agreement between the absorption spectra from aiMD simulations, which fully account for electronic and nuclear polarizations, and from experiments is good, which suggests that the chosen simulation method is well-suited for modeling IR spectra, although the agreement is known to be partly due to a cancellation of approximations in the employed density functional theory (DFT) and the neglect of nuclear quantum effects.[48,49] Molecular simulations of liquids including nuclear quantum effects have rather recently become feasible, mostly via centroid and ring-polymer molecular dynamics (MD) techniques.[50] Unfortunately, the projection techniques we apply on the classical nuclear trajectories from our aiMD simulations are not available on the quantum level; therefore, the analysis we do in this paper can currently not be done for quantum systems without additional drastic and uncontrolled approximations. Nevertheless, the interplay of potential and frequency-dependent friction effects we explore in this paper presumably is not modified by nuclear quantum effects in a fundamental way, so that the conclusions we draw with regard to the importance of frequency-dependent friction effects should remain valid even beyond our classical treatment of the nuclei.

Figure 2

Absorption spectra of aiMD simulations at 300 K of 256 H2O molecules are shown as gray solid lines in panels A–C. The spectra up to 2000 cm–1 are multiplied by a factor three. (A) Comparison to aiMD spectra for liquid D2O (blue solid line) and experimental data (obtained for 298 K), shown as a gray broken line for H2O[43] and a blue broken line for D2O.[44] (B) Comparison to aiMD simulations of a single H2O (blue solid line). The normal-mode frequencies of a single H2O are shown as vertical dotted lines. (C) Comparison to power spectra of the symmetric stretch, dOHs (blue solid line); antisymmetric stretch, dOHa (blue broken line); and bend mode vibrations, ϕHOH (blue dotted line), which are averaged over all molecules and rescaled to match the absorption spectrum.

Figure B compares simulated liquid H2O (gray) and single H2O (blue solid line) spectra at 300 K. The single-water spectrum shows sharp peaks which perfectly coincide with the normal-mode frequencies of a single water molecule (vertical dotted lines, computed on the same DFT level as the aiMD simulations) at 1607, 3675, and 3772 cm–1, which are within 20 cm–1 of the experimental values 1594.7, 3657.1, and 3755.9 cm–1.[51,52] Note that the OH-stretch band consists of two modes, namely, the low-frequency symmetric mode, where both OH bonds vibrate in phase, and the high-frequency antisymmetric mode, where the OH bonds vibrate out of phase, which do not clearly separate in the liquid spectrum. The symmetric stretch mode in the gas phase shows a much smaller intensity than the antisymmetric stretch mode, in agreement with experiment,[52] which is caused by electronic polarization effects. The OH-stretch peak in the liquid is significantly red-shifted and enhanced compared to gas phase, which is typically rationalized by the softening of the OH bond potential and the constructive collectivity of OH-stretching vibrations in the liquid (see section II in the Supporting Information);[2,3,53] the significant enhancement is noteworthy, because one could expect the friction acting on the OH bond to be much stronger in the liquid and thus to reduce the vibrational amplitude. In contrast, the HOH-bending mode in the liquid is slightly blue-shifted and not enhanced, which can be rationalized by collective effects that are slightly destructive (see section II in the Supporting Information). All these effects are fully accounted for by the frequency-dependent friction acting on the different vibrational modes, as explained below.

The vibrational modes of a water molecule can be described by the bond angle ϕHOH and the symmetric and antisymmetric stretch distances, dOHs = (dOH1 + dOH2)/2 and dOHa = (dOH1 – dOH2)/2, where the two OH bond distances in a water molecule are denoted as dOH1 and dOH2, which are all based on the nuclear positions in the aiMD simulations, as illustrated in the inset in Figure C. The power spectra of these three vibrational modes, averaged over all water molecules in the liquid, are shown in Figure C (ϕHOH as dotted, dOHs as solid, and dOHa as broken blue lines) and compared to the absorption spectrum from the total dipole moment. The agreement of the line frequencies and shapes is quite good, except that the absorption spectrum is red-shifted compared to the dOHs and dOHa vibrational spectra. This red shift is due to dipolar correlations between neighboring water molecules (as mentioned above and discussed in detail in section II in the Supporting Information) and also due to electronic polarization effects, as shown in section III in the Supporting Information. The dOHs and dOHa spectra overlap significantly, with a small red shift of the dOHs spectrum relative to the dOHa spectrum, in accordance with previous observations.[54] The vibrational spectrum of the ϕHOH mode overlaps perfectly with the spectrum from the total (nuclear and electronic) dipole moment, which is due to the fact that the bending angle vibrations of neighboring water molecules are only weakly (and in fact anti-) correlated, as shown in section II in the Supporting Information.[53] We conclude that the absorption spectrum calculated from the total system polarization (including nuclear and electronic polarization from all water molecules and their correlations) match the power spectra based on the single-water nuclear-coordinate-based vibrational modes rather faithfully. This good agreement lies at the heart of the common interpretation of IR absorption spectra in terms of molecular vibrations; it also validates our approach, because it means that the conclusions from our GLE analysis, which in the present formulation can be applied only on one-dimensional reaction coordinates that are derived from nuclear positions, can also be used to interpret simulated and experimental absorption spectra.

In the following, we will analyze the dynamics of different water vibrational modes based on the one-dimensional GLE[55,56]which contains an in-general nonharmonic time-independent potential U(x) that corresponds to a free energy as it results from integrating out all other degrees of freedom except x(t). The memory kernel Γ(t) describes the time-dependent friction acting on the fluctuating variable x(t), which can be the bond angle (ϕHOH), the symmetric stretch distance (dOHs), or the antisymmetric stretch distance (dOHa). The random force FR(t) has zero mean ⟨FR(t)⟩ = 0 and fulfills the fluctuation–dissipation relation ⟨FR(t)FR(t′)⟩ = kBTΓ(t – t′). The GLE approach as introduced by Mori and Zwanzig in the 1960s is an exact projection of the full dynamics of a multiparticle system onto a reduced set of coordinates. Given a one-dimensional trajectory x(t) from the aiMD simulations, the effective mass m, the potential U(x), and the friction function Γ(t) are uniquely determined and can be accurately extracted,[40] as described in sections IV and V in the Supporting Information. It follows that all intra- and intermolecular interaction effects on the molecular vibration dynamics are accurately taken into account: as a crucial test of the validity of the GLE in the formulation of eq , of our extraction methods, and of our simulation methods of the GLE, we will further demonstrate below that the GLE accurately reproduces the vibrational mode spectra calculated directly from the aiMD simulations. Therefore, our description of the vibrational dynamics of a water molecule in the liquid phase via the GLE is accurate, or, more strictly speaking, it accurately reproduces the vibrational dynamics obtained in our ab initio simulations.

For a harmonic potential, U(x) = kx2/2, the vibrational power spectrum that follows from the GLE eq can be given in closed form as (see section VI in the Supporting Information)where the frequency-dependent friction is obtained by a single-sided Fourier transform Γ̃(ω) = ∫0∞ dt eiωΓ(t). In the limit of frequency-independent friction Γ̃(ω) = γ, this yields the standard Lorentzian line shape[57] (see section VII in the Supporting Information)which will be shown below to give only a poor account of our simulated vibrational spectra. Nonharmonic potentials are parametrized aswhere x0 is the position of the minimum of U(x). Vibrational spectra in the presence of nonharmonic potentials are obtained from numerical simulations of the GLE using a parametrized friction function of the form[58−60]consisting of n exponentially decaying components with time scales τe and friction coefficients γ as well as l oscillating and decaying components with amplitudes a, oscillation frequencies ω, and decay time scales τo (see sections VIII–X in the Supporting Information for details).

Methods

All Born–Oppenheimer aiMD simulations are performed with the CP2K 4.1 software package using a contracted double-ζ basis set for the valence electrons, optimized for small molecules and short ranges (DZVP-MOLOPT-SR-GTH), dual-space pseudopotentials, the BLYP exchange–correlation functional, D3 dispersion correction, and a cutoff for the plane-wave representation set to 400 Ry.[61−63] A time step of 0.5 fs is used under NVT conditions at 300 K by coupling all atoms to a CSVR thermostat with a time constant of 100 fs.[64] The bulk systems contain 256 molecules subject to periodic boundary conditions in a cubic cell of size (1.9734 nm)3, corresponding to densities of 996.4 kg/m3 for H2O and 1107.8 kg/m3 for D2O. The total trajectory lengths of the liquid systems are 230 ps for H2O and 130 ps for D2O. Simulations of a single H2O molecule, representing the gas phase data, are performed in the NVE ensemble with 47 initial configurations sampled from a 25 ps NVT simulation using an individual thermostat with a time constant of 10 fs for each atom. The NVE simulations are each run for 10 ps with a time step of 0.25 fs. The distributions of their initial configurations sample well the equilibrium distributions as shown in section XI of the Supporting Information.

Linear response theory relates the dielectric susceptibility χ(t) to the equilibrium autocorrelation of the dipole moment C(t) = ⟨(t)(0)⟩, reading in Fourier space[65] (see section I in the Supporting Information)with system volume V, thermal energy kBT, and vacuum permittivity ϵ0. IR absorption spectra can therefore be calculated straight-forwardly from sufficiently long trajectories from aiMD simulation data using eq and the Wiener–Khintchine relation,[53,66] derived in section XII in the Supporting Information. Quantum corrections have previously been addressed, but are not applied here.[67] The molecular dipole moments are obtained after Wannier-center localization of the electron density at a time resolution of 2 fs. The Wannier centers are assigned to the molecule of the nearest oxygen, which always results in exactly four Wannier centers per water molecule. A charge of −2e is assigned to each Wannier center, which together with the nuclear charges, reduced by the electronic charges of the inner shells, allows for the calculation of the dipole moment. The power spectra are smoothed using a Gaussian kernel with width 10 cm–1. The normal-mode analysis is performed for an energetically minimal configuration of a single H2O using the implementation in CP2K 4.1 and the same ab initio model as for the aiMD simulation.

Results and Discussion

We start with a discussion of the symmetric stretch mode dOHs. The potential (which actually corresponds to a free energy) U(dOHs) from the aiMD simulations for liquid water (gray solid line) is in Figure A compared with a nonharmonic fit according to eq (blue broken line); the harmonic contribution is shown as a gray dotted line. The comparison of the liquid and gas-phase bond potentials in Figure B shows that the minimum of the potential (i.e., the most probable OH bond length) increases from x0 = 97.50 pm in the gas phase to x0 = 99.25 pm in the liquid; at the same time the harmonic force constant decreases from k/kBT = 0.404 pm–2 in the gas phase to k/kBT = 0.274 pm–2 in the liquid. This softening of the potential is due to elongation of the bond, caused by hydrogen bonding in the liquid, and will in the absence of frequency-dependent friction effects be shown to induce a pronounced spectral red shift. Furthermore, the potential nonharmonicity increases, as can be seen by comparing the reduced cubic potential coefficient in the liquid phase k̃3 = k3/kBT(k/kBT)−3/2 = −0.0840 with the value in the gas phase k̃3 = −0.0485.

Figure 3

Results for the symmetric stretch coordinate dOHs from aiMD simulations. (A) Potential U(dOHs) for 256 H2O in the liquid phase (gray solid line) compared to the nonharmonic fit according to eq (blue broken line) and the harmonic part (gray dotted line). (B) Potential U(dOHs) for a single H2O molecule in the gas phase (green solid line) compared with the nonharmonic fit according to eq (blue broken line) and the harmonic part (gray dotted line). The liquid-phase potential (gray solid line) is shown for comparison. (C and D) Friction as a function of time and frequency (gray lines) compared with the fit according to eq (blue lines). Real and imaginary parts in panel D are shown as solid and broken lines, respectively; the spectrum on top is the full absorption spectrum from aiMD. The blue dotted line in panel C shows a single exponential with decay time τ = 10 fs; the dotted horizontal line in panel D shows the constant real friction γOHs = Γ̃′(fOHs) evaluated at the symmetric OH stretch vibrational frequency fOHs = 3390 cm–1. The gray circle denotes the static friction Γ̃′(0). (E) Vibrational power spectrum ωχ̃″ (gray solid line) compared to models of varying complexity: normal mode of single H2O (broken vertical line), Lorentzian with harmonic potential and constant friction γOHs (gray dotted line), nonharmonic potential and constant friction γOHs (blue broken line), harmonic potential and frequency-dependent friction Γ̃(f) (purple solid line), and nonharmonic potential and frequency-dependent friction Γ̃(f) (blue solid line). (F) Vibrational power spectrum ωχ̃″ using the nonharmonic potential and different values of the constant friction γ, where γOHs = Γ̃′(fOHs) is the friction evaluated at the symmetric OH stretch vibrational frequency. The gray solid line is the spectrum from aiMD simulations.

The time-dependent friction function for the symmetric stretch mode Γ(t) extracted from aiMD simulations (gray line in Figure C) shows multiexponential decay characterized by relaxation times from a few femtoseconds to many picoseconds, which is appreciated by comparison with a single-exponential function with decay time τ = 10 fs (dotted blue line, the logarithmic time axis should be noted). This in particular means that Γ(t) accounts for solvent relaxations that are equally fast (adiabatic) and slower (nonadiabatic) compared to the OH vibrational period of about 10 fs; one thus expects homogeneous as well as inhomogeneous line broadening to occur, as indeed borne out by our analysis below. The oscillations that appear in Γ(t) at around 10–250 fs reflect the dissipative coupling of symmetric stretch vibrations to antisymmetric stretch as well as higher-harmonic bend and librational modes. This is illustrated by the real and imaginary frequency-dependent friction components Γ̃′(ω) + iΓ̃″(ω) = ∫0∞ dt eiω Γ(t) in Figure D (solid and broken gray lines, respectively), which exhibit maxima at the OH-stretching and HOH-bending frequencies and also at their higher harmonics.

The friction function thus accounts for the frequency-dependent vibrational energy dissipation of a given vibrational mode within a water molecule as well as into the surrounding water and in particular accounts for resonances between different vibrational modes and their overtones. The resonances contained in the friction function thus are equivalent to Fermi resonances,[32−35] which typically arise in a quantum formulation and in our classical description are caused by nonlinear intra- and intermolecular couplings in the multidimensional potential landscape that describes the nuclear vibrations. Also non-Condon effects, which arise because of modifications of the transition dipole moment of a vibrational mode due to time-dependent changes of the solvation environment of a molecule,[36] are included via the interplay of the potential U(dOHs) and the time-dependent friction function Γ(t). Interestingly, the symmetric stretch shows a much stronger frictional damping at the characteristic frequency of the bending mode than the antisymmetric stretch mode, shown in section XIII of the Supporting Information, which points to a stronger dissipative coupling of bending vibrations with symmetric than with antisymmetric stretch vibrations.

For simulations of the GLE, which are necessary for the analysis of the coupling between nonlinearities in the potential and frequency-dependent friction, we fit Γ̃′(ω) by the expression eq with a sum of three exponential and six oscillating functions, see section VIII of the Supporting Information for details. The fit shown in blue in Figure C,D describes the simulated friction function equally well in the time as well as in the frequency domain.

The vibrational spectrum of the dOHs mode directly extracted from aiMD simulations is shown in Figure E as a gray solid line. The simplest possible model for a vibrational line shape is the Lorentzian model eq for a harmonic potential and a constant, frequency-independent friction. Using k/kBT = 0.274 pm–2 from the harmonic fit in Figure A and the friction γOHs = Γ̃′(fOHs) in Figure D at the stretch vibrational frequency fOHs = 3390 cm–1, we obtain the gray dotted line in Figure E. Compared to the normal-mode frequency of the gas phase, denoted by a vertical green broken line, the Lorentzian is significantly red-shifted by about 500 cm–1; the width of the Lorentzian reflects homogeneous line broadening due to adiabatic solvent friction that is described by the frequency-independent constant γOHs. Note that the Lorentzian is considerably red-shifted and narrower compared to the spectrum extracted from the aiMD simulation (gray line). Interestingly, the friction γOHs that acts at the vibration frequency fOHs is about 2 orders of magnitude smaller than the friction in the static limit f = 0, as seen in Figure D, which explains why the stretch vibrational dynamics shown in Figure A is rather weakly damped. The vibrational power spectrum in the presence of the full nonharmonic potential U(dOHs) and constant friction γOHs, obtained from numerical simulations of the memoryless Langevin equation (blue broken line, see section IX in the Supporting Information for details), is only slightly red-shifted with respect to the Lorentzian obtained for a harmonic potential, which is expected based on perturbation theory.[4] We conclude that nonlinearities in the potential have for constant friction only an insignificant influence on the line frequency and shape. The peak frequency of a Lorentzian does not depend on the value of the constant friction γ (see section XIV in the Supporting Information), which is approximately true also in the presence of the nonharmonic potential U(dOHs), as demonstrated in Figure F where spectra from numerical simulations for varying γ are compared. We next check for the influence of time-dependent friction on the spectrum. For a harmonic potential and for frequency-dependent friction, the spectrum is determined analytically by eq and shown in Figure E as a purple solid line. A significant blue shift compared to the results for constant friction is obtained, so that the position of the spectrum agrees very well with the simulated spectrum, while the line shape is too narrow. The blue shift can be understood based on simple and rather general analytic arguments, as shown below. The spectrum obtained from the GLE in the presence of the nonharmonic potential U(dOHs) and time-dependent friction Γ(t), shown by the blue solid line in Figure E (here numerical simulations are employed), is significantly broadened compared to the results obtained for a harmonic potential and time-dependent friction Γ(t) (purple line). This reflects the effects of inhomogeneous line broadening,[4] i.e., the effects of a slowly varying hydration environment of a vibrating bond that elongates or compresses the bond length in conjunction with a nonlinear bond potential, and reproduces the spectrum extracted from the aiMD simulations (gray line) almost perfectly; in fact, inhomogeneous line broadening is quite substantial and accounts for 52% of the total line broadening. This means that the GLE, when used in conjunction with the properly extracted nonharmonic time-averaged potential U(dOHs) and time-dependent friction Γ(t), reproduces the system dynamics very well, which is not guaranteed in general because the projection onto the GLE neglects nonlinear friction effects.[31]

The ϕHOH water bending coordinate is analyzed analogously: The bend angle potential U(ϕHOH) in Figure A extracted from aiMD simulations (gray line) includes significant nonquadratic contributions as appreciated by a comparison of the nonharmonic fit (blue broken line) with the harmonic part (dotted line) and as witnessed by the magnitude of the reduced cubic and quartic fit parameters k̃3 = −0.0596 and k̃4 = k4/kBT(k/kBT)−2 = −0.00323. Different from the situation for the stretch potential, the liquid environment shifts the most probable bending angle only very slightly. Nevertheless, the potential is softened considerably, as is seen by a comparison of the shape and fit parameters of the gas and liquid-phase potentials U(ϕHOH) in Figure B, which can be rationalized by the fact that attractive electrostatic interactions, which are predominant for strongly correlated polar liquids such as water, exhibit negative curvature throughout their entire interaction range. The time-dependent friction Γ(t) extracted from the simulations in Figure C (gray line) shows a broad decay but more pronounced oscillations compared to the stretch vibrations in Figure C. The fit (blue solid line) to the simulated real frequency-dependent friction Γ̃′(f) (gray solid line) in Figure D requires two exponential and 14 oscillatory functions to describe the simulated data satisfactorily, see section VIII in the Supporting Information for details. The dissipative damping is significantly more pronounced at stretch frequencies around 3400 cm–1 and at the overtones of the bending around 3300 cm–1 and around 4950 cm–1 than at the bending fundamental around 1650 cm–1 itself, indicative of the nonlinear coupling between different modes and overtones (where it should be noted that coupling of bend vibrations to higher-frequency modes and overtones are reduced when quantum effects are properly included[16]).

Figure 4

Results for the bend coordinate ϕHOH from aiMD simulations. (A) Potential U(ϕHOH) for 256 H2O molecules in the liquid phase (gray solid line) compared to the nonharmonic fit according to eq (blue broken line) and the harmonic part (gray dotted line). (B) Potential U(ϕHOH) for a single H2O molecule in the gas phase (green solid line) compared with the nonharmonic fit according to eq (blue broken line) and the harmonic part (gray dotted line). The liquid-phase potential (gray solid line) is shown for comparison. (C and D) Friction as a function of time and frequency (gray lines) compared with the fit according to eq (blue lines). Real and imaginary parts in panel D are shown as solid and broken lines; dotted lines denote negative values of the imaginary part; the spectrum on top is the full absorption spectrum from aiMD. The dotted horizontal line in panel D shows the constant real friction γHOH = Γ̃′(fHOH) evaluated at the bend vibrational frequency fHOH = 1650 cm–1. The gray circle denotes the static friction Γ̃′(0). (E) Vibrational power spectrum ωχ̃″ (gray solid line) compared to models of varying complexity: normal mode of single H2O (broken vertical line), Lorentzian with harmonic potential and constant friction γHOH (gray dotted line), nonharmonic potential and constant friction γHOH (blue broken line), harmonic potential and frequency-dependent friction Γ̃(f) (purple solid line), and nonharmonic potential and frequency-dependent friction Γ̃(f) (blue solid line). (F) Vibrational power spectrum ωχ̃″ using the harmonic potential part and the constant friction γHOH (gray dotted line), the real frequency-dependent friction only Γ̃′(f) (blue broken line), and the real and imaginary frequency-dependent friction Γ̃′(f) + iΓ̃″(f) (blue solid line). The gray solid line is the spectrum from aiMD simulations.

The vibrational spectrum of the ϕHOH coordinate from the aiMD simulations is shown in Figure E as a gray solid line and is weakly blue-shifted from the gas-phase normal mode (vertical green broken line), which is a surprising fact and will be explained now by compensatory potential and friction effects. The spectrum from the Lorentzian model eq (gray dotted line) using only the harmonic potential part of U(ϕHOH) and the frequency-independent friction γHOH = Γ̃′(fHOH), obtained at the bending peak at fHOH = 1650 cm–1 (horizontal broken line in Figure D), is significantly red-shifted and is not modified much by including the nonharmonic potential contributions (blue broken line). Upon including the complex frequency-dependent friction Γ̃(f) but only the harmonic part of U(ϕHOH), the purple line is obtained, which is blue-shifted with respect to the constant-friction case and reaches the frequency of the simulated curve but is too narrow. Including the complex frequency-dependent friction Γ̃(f) and also the full nonharmonic potential U(ϕHOH), the GLE (indicated by the blue line) rather accurately reproduces the position and width of the simulated spectrum. In agreement with our stretch-vibration results in Figure E, we detect considerable inhomogeneous line broadening (amounting to 47% of the total line broadening) from the comparison of the results with and without nonharmonic potential contributions in the presence of frequency-dependent friction. In contrast to the stretch–vibration results, we see that the blue shift induced by including the frequency dependence of the friction almost exactly cancels the red shift due to the softening of the bond potential in the liquid phase, which means that the frequency dependence of Γ̃(f) close to the characteristic bend-mode frequency is more pronounced compared to the stretch mode. It transpires that the fine details of the frequency dependence of the friction at the vibrational frequency determine vibrational line shape and position, which we now analyze in more detail.

It turns out that the imaginary and real parts of the frequency-dependent friction influence the line position and shape quite differently.[26] This is illustrated in Figure F by comparing spectra using only the harmonic part of the potential for constant friction (gray dotted line), for purely real frequency-dependent friction Γ̃′(f) (blue broken line), and for friction that contains both real and imaginary frequency-dependent parts Γ̃′(f) + iΓ̃″(f) (blue solid line), note that for purely imaginary friction the spectrum according to eq exhibits a singularity and thus is not shown. It is in fact the imaginary part Γ̃″(f) that gives rise to the blue shift, as is now explained by a simple analytical argument.

For this we consider a single-exponential memory function Γ(t) = γτ–1 exp(−t/τ). The single-sided Fourier transform is given as Γ̃(ω) = ∫0∞ dt eiωΓ(t) = γ/(1 – iτω) with the asymptotic limits Γ̃(ω) ≃ γ(1 + iωτ) for small ω and Γ̃(ω) ≃ iγ/(ωτ) for large ω; both deviations from the zero-frequency limit Γ̃(ω → 0) ≃ γ turn out to be imaginary, which already hints at why the imaginary part of the friction determines the line position, as demonstrated in Figure F. A general form that contains both asymptotic limits is given by Γ̃(ω) ≃ γ + iaω + ib/ω, where a = γτ and b = 0 for small ω and a = 0 and b = γ/τ for large ω. By inserting this asymptotic form into eq , the Lorentzian line shape eq is recovered but with an effective mass meff = m – a and an effective potential curvature keff = k + b. The vibrational frequency turns out to beand in fact increases both in the small and large frequency limits, because a and b are positive constants for single-exponential memory. Thus, a blue shift of the vibrational frequency is very generally expected for frequencies where the frequency-dependent friction is described by the asymptotic form Γ̃(ω) ≃ γ + iaω + ib/ω with positive a and b. In fact, this functional form is able to describe the stretch and band friction functions rather accurately around the stretch and band frequencies, respectively, as inspection of Figures D and 4D shows. In other words, frequency-dependent friction can lead to a blue shift of a vibrational band irrespective of whether the imaginary friction function decreases or increases in the vicinity of the vibrational band.

The full width at half-maximum of a Lorentzian is given as γ/meff ≃ γ/(m – a); thus, the line width is, within the harmonic approximation, predicted to slightly increase for the stretch band (because Γ̃″(f) slightly increases at the stretch vibrational frequency in Figure D and thus a is positive) but to stay rather constant for the bend band (because Γ̃″(f) slightly decreases at the bend vibrational frequency Figure D and thus a presumably is small and dominated by b). These predictions are in good agreement with the results shown in Figures E and 4E for the scenario of a harmonic potential and friction-dependent friction (purple lines). Clearly, the exact line shape and position are determined by the interplay of nonlinearities of the potential and frequency-dependent friction, which can be accessed only by simulations or perturbation theory, but the simple harmonic model discussed here allows appreciation of part of the mechanisms at play.

Conclusions

Frequency-independent friction, which reflects the fast adiabatic dissipative channels available for a specific vibrational mode, modifies the vibrational spectral line width via homogeneous line broadening but not the line position. This is strictly true for a harmonic potential but holds approximately even in the presence of nonharmonic potential contributions, as we demonstrate in Figure F. On the other hand, the full frequency dependence of the friction, which in particular accounts for the slower solvent relaxation processes, gives rise to a spectral blue shift and additional line broadening. The latter reflects what is typically called inhomogeneous line broadening. In contrast, softening of the bond potential in the liquid environment, which is due to hydrogen bonding and hydration interactions, gives rise to a red shift. Therefore, we find that the line shapes and positions of the bend and stretch bands in liquid water can be interpreted in terms of the compensatory effects of frequency-dependent friction and harmonic as well as nonharmonic potential contributions. For stretch vibrations, the bond softening dominates and therefore the stretch vibration is red-shifted when going from gas to liquid water. For bend vibrations the potential-induced red shift and the friction-induced blue shift almost exactly compensate. This of course does not imply that the coupling of bend vibrations to the hydrating liquid environment is weaker than for stretch vibrations, as one might naively guess from only looking at the frequency shifts. Rather, the contrary is true. It turns out that it is the imaginary part of the frequency-dependent friction that gives rise to the blue shift, in line with previous arguments.[26] The situation is rather complex, though, because the effects due to the frequency-dependency of the friction and due to nonlinearities in the potential do not decouple.

Our methodology is different from previous approaches to describe the infrared line shapes of water[6,68] because we deploy the time-averaged bond potential as it naturally emerges via the projection formalism used to derive the GLE. This in particular means that in our approach, inhomogeneous line broadening enters via the time-dependent friction function, not via a time-dependent bond potential, as in previous theories.

The GLE framework we use to derive the time- or frequency-dependent friction function is not constrained to nuclear reaction coordinates that describe the vibrations of a single molecule, which forms the topic of this paper. Rather, electronic polarization degrees of freedom can be included as well and also collective effects that stem from dipolar correlations between neighboring molecules can be accounted for by suitable reaction coordinates. Likewise, it would be interesting to model Raman spectra, which reveal a perspective on the vibrational molecular modes that is very different from IR spectroscopy.[69]

As mentioned before, our ab initio simulations neglect nuclear quantum effects, owing to the fact that methods to extract friction functions from path-integral simulations are not yet available. This approximation presumably is permissible in the present context, as we target the general compensatory effects the liquid environment has on bond potentials and the bond friction function, which should not be fundamentally changed by nuclear quantum effects. This is corroborated by results by Marsalek and Markland,[49] who reported a red shift of the OH peak by about 200 cm–1 and of the HOH bend peak by about 50 cm–1 upon inclusion of nuclear quantum effects in their simulations, which are somewhat smaller than the shifts due to potential and frequency-dependent friction effects we find here. Nonetheless, in the future, it would be highly desirable to develop techniques that would allow extracting GLE parameters from path integral simulations[48,49] and from mixed quantum/classical approaches.[32,34,70]