# Time-dependent friction effects on vibrational infrared frequencies and line
shapes of liquid water

Florian N. Brünig  Freie Universität Berlin, Germany  Otto Geburtig  Freie
Universität Berlin, Germany  Alexander von Canal  Freie Universität Berlin,
Germany  Julian Kappler  Freie Universität Berlin, Germany  Roland R. Netz
<EMAIL_ADDRESS>Freie Universität Berlin, Germany

aiMD

    ab initio Molecular Dynamics
ATR

    attenuated total reflection
DFT

    density functional theory
FPT

    first-passage time
GLE

    generalized Langevin equation
GH

    Grote-Hynes
FTIR

    Fourier-transform infrared spectroscopy
FWHM

    full width at half maximum
IR

    infrared
LE

    Langevin equation
MD

    molecular dynamics
MFPT

    mean first-passage time
MFP

    mean first-passage
MSD

    mean squared displacement
NM

    normal-mode
PTP

    $p(\text{TP}|q)$
PME

    particle mesh Ewald [pronk2013gromacs]
PMF

    potential of mean force
PGH

    Pollak-Grabert-Hanggi
RC

    reaction coordinate
RDF

    radial distribution function
RTT

    round-trip time
TP

    transition path

##  I Infrared power spectra from linear-response theory

Assuming linear response of an observable $x(t)$ with respect to a force that
couples to an observable $y(t)$, the response function $\chi_{xy}(t)$ is
related to the correlation function $C_{xy}(t^{\prime})=\langle
x(t+t^{\prime})y(t)\rangle$ for $t\geq 0$ [1]

$\displaystyle\chi_{xy}(t)=-\frac{1}{k_{B}T}\frac{d}{dt}C_{xy}(t),$ (S1)

where $k_{B}T$ is the thermal energy. Realizing that $\chi(t)$ is single-
sided, i.e. $\chi(t)=0$ for $t<0$, the Fourier transform is calculated as

$\displaystyle\widetilde{\chi}_{xy}(\omega)$
$\displaystyle=-\frac{1}{k_{B}T}\int_{-\infty}^{\infty}dt\ e^{i\omega
t}\frac{d}{dt}C_{xy}(t)$ $\displaystyle\widetilde{\chi}_{xy}(\omega)$
$\displaystyle=-\frac{1}{k_{B}T}\left(C_{xy}(0)-i\omega\int_{0}^{\infty}dt\
e^{i\omega t}C_{xy}(t)\right)$ $\displaystyle\widetilde{\chi}_{xy}(\omega)$
$\displaystyle=-\frac{1}{k_{B}T}\left(C_{xy}(0)-i\omega\widetilde{C}_{xy}^{+}(\omega)\right),$
(S2)

where the superscript + denotes a single-sided Fourier transform. In case of
$x=y$, $C_{xx}(t)$ is an autocorrelation function, which is real and
symmetric, therefore it follows for the imaginary part of the response
function in Fourier space

$\displaystyle\widetilde{\chi}_{xx}^{\prime\prime}(\omega)$
$\displaystyle=\frac{1}{k_{B}T}\omega\operatorname{Re}(\widetilde{C}^{+}_{xx}(\omega))$
(S3)
$\displaystyle=\frac{1}{k_{B}T}\frac{\omega}{2}\widetilde{C}_{xx}(\omega).$
(S4)

When computing the power spectra of a stochastic process $x(t)$, limited to
the time domain $[0,L_{t}]$, the Wiener-Khintchine theorem, eq. (S97) in
section XIV, can be used to express $\widetilde{C}_{xx}(\omega)$ in terms of
$\tilde{x}(\omega)$, turning eq. (S4) into

$\displaystyle\widetilde{\chi}_{xx}^{\prime\prime}(\omega)=\frac{\omega}{2k_{B}TL_{t}}|\tilde{x}(\omega)|^{2}.$
(S5)

When computing the power spectra of the observable $x(t)$ from the ensemble
average of equilibrium trajectories, a decomposition of $x(t)$ into two parts
$x(t)=x_{1}(t)+x_{2}(t)$ gives rise to three contributions in the total power
spectrum

$\displaystyle\omega\widetilde{\chi}_{xx}^{\prime\prime}(\omega)$
$\displaystyle=\frac{\omega^{2}}{2k_{B}T}\left[\widetilde{C}_{1}(\omega)+\widetilde{C}_{2}(\omega)+2\widetilde{C}_{1,2}(\omega)\right]$
$\displaystyle=\omega\left[\widetilde{\chi}^{\prime\prime}_{1}(\omega)+\widetilde{\chi}^{\prime\prime}_{2}(\omega)+\widetilde{\chi}^{\prime\prime}_{1,2}(\omega)\right],$
(S6)

where the cross-correlation contribution $\chi^{\prime\prime}_{1,2}(\omega)$
was defined such that it equals the difference spectrum

$\displaystyle\widetilde{\chi}^{\prime\prime}_{\mathrm{diff}}(\omega)$
$\displaystyle=\widetilde{\chi}_{xx}^{\prime\prime}-\widetilde{\chi}^{\prime\prime}_{1}(\omega)-\widetilde{\chi}^{\prime\prime}_{2}(\omega)=\widetilde{\chi}^{\prime\prime}_{1,2}(\omega)$
$\displaystyle=\frac{\omega}{k_{B}T}\widetilde{C}_{1,2}(\omega).$ (S7)

A positive cross-correlation spectrum hints to in-phase motion, a negative
cross-correlation spectrum to out-of-phase motion of $x_{1}(t+t^{\prime})$ and
$x_{2}(t)$ at a given frequency.

In case of $x(t)$ being the polarization $\bm{p}(t)$ of the system, coupling
to an external electric field $\bm{E}(t)$, the dimensionless dielectric
susceptibility $\chi(t)$ is given by

$\displaystyle\widetilde{\chi}(\omega)=\frac{1}{V\epsilon_{0}D}\langle\widetilde{\chi}_{\bm{p}\bm{p}}(\omega)\rangle,$
(S8)

where $\epsilon_{0}$ is the vacuum permittivity, $V$ is the system volume and
an average is performed over the $D$ dimensions of $\bm{p}$.

##  II Decomposition of water infrared spectra into single-molecular and
collective components

\begin{overpic}[width=212.47617pt]{{pc_spectra/spec_dft_nvt_256}.pdf}
\put(2.0,55.0){\large\bf A}
\end{overpic}\begin{overpic}[width=212.47617pt]{{pc_spectra/spec_dft_nvt_256_nuc}.pdf}
\put(2.0,55.0){\large\bf B} \end{overpic}

Figure S1:  Decomposition of the IR spectra into total (grey solid lines),
single-molecular (blue solid lines) and collective contributions (purple solid
lines), obtained from the total dipole-moment trajectory from the aiMD
simulation of 256 H2O molecules including nuclear and electronic charges after
Wannier localization (A) and an approximation using partial charges on the
nuclear dynamics of the same aiMD simulation (B, see the following section III
for details).

Fig. S1A and B shows a complete decomposition of the total (nuclear and
electronic) IR spectra obtained from aiMD simulations of 256 H2O molecules
(grey solid lines), as introduced in the main text, into single-molecular
(blue solid lines) and collective components (purple solid lines), which
follow from the molecular dipole moments $\bm{p}_{i}(t)$ of all the molecules
in the bulk as [2]

$\displaystyle\widetilde{\chi}^{\prime\prime}_{\rm tot}(\omega)$
$\displaystyle=\widetilde{\chi}^{\prime\prime}_{\rm
mol}(\omega)+\widetilde{\chi}^{\prime\prime}_{\rm coll}(\omega)$ (S9)
$\displaystyle\sim\langle\sum_{i}\bm{p}_{i}(0)\sum_{j}\bm{p}_{j}(t)\rangle$
$\displaystyle=\sum_{i}\langle\bm{p}_{i}(0)\bm{p}_{i}(t)\rangle+\sum_{i}\langle\bm{p}_{i}(0)\sum_{j\neq
i}\bm{p}_{j}(t)\rangle.$ (S10)

The data clearly shows, that whereas in the OH stretching regime the
collective effects are constructive and lead to an amplification of the total
infrared spectrum with respect to the single-molecular spectrum, in the HOH
bending regime the collective effects are destructive and lead to a slight
decrease of the amplitude of the total IR spectrum with respect to the single-
molecular spectrum. In the regime of the librations between
$300\text{\,}\mathrm{c}\mathrm{m}^{-1}800\text{\,}\mathrm{c}\mathrm{m}^{-1}$
the collective effects contribute constructive as well as destructive. At the
small signature around $200\text{\,}\mathrm{c}\mathrm{m}^{-1}$, associated
with translational vibrations of water molecules against each other, the
collective effects contribute constructive.

##  III Infrared spectra from nuclear coordinates

\begin{overpic}[width=216.81pt]{{spectra_h2o_classical}.pdf} \end{overpic}
Figure S2: Comparison of IR spectra obtained from the total dipole-moment
trajectory including nuclear and electronic charges from the aiMD simulation
(grey solid line) and an approximation using partial charges on the nuclear
dynamics of the same aiMD simulation (blue solid line).

Fig. S2 shows a comparison of the IR spectrum of the aiMD simulation, using
the total dipole moment including nuclear and electronic charges and a nuclei-
only power spectrum, which is obtained by assigning the partial charges
$1\text{\,}\mathrm{e}$ to the hydrogen atoms and $-2\text{\,}\mathrm{e}$ to
the oxygen atoms of the aiMD simulation in post processing. Both spectra are
qualitatively similar but the nuclei-only spectrum is significantly increased
at lower frequencies, due to the neglect of electronic polarization effects.
In the OH-stretching regime at around $3300\text{\,}\mathrm{c}\mathrm{m}^{-1}$
both spectra coincidentally have a similar amplitude. Interestingly the
nuclei-only spectrum is blue-shifted with respect to the spectrum from nuclear
and electronic charges, indicating a slow down of the total dipole-moment
dynamics in this regime due to the electronic degrees of freedom.

##  IV Extraction and parametrization of the potentials

From the entire set of trajectories of a vibrational coordinate a histogram is
created with 50 equidistant bins centered at $x_{i}$. From this histogram
$h(x_{i})$ the potential (or free energy) is calculated as

$\displaystyle U(x_{i})=-\log(h(x_{i}))k_{B}T,$ (S11)

where the potential if shifted vertically so that $\min(U[x_{i}])=0$. The
largest and smallest value of $x_{i}$ where $U(x_{i})<8\,k_{B}T$ are the edges
for a new histogram with 99 equidistant bins. From this histogram the
potential is again calculated as shown above. This potential is then fitted to
a 4th-order polynom, $U(x)=a+bx+cx^{2}+dx^{3}+ex^{4}$, using the Levenberg-
Marquardt algorithm. The position of the minimum of the fit function, $x_{0}$,
is subsequently obtained using Newton’s method, both implemented in scipy
v1.5. Eventually, $U(x)$ is rewritten as

$\displaystyle U(x)=$ $\displaystyle
k_{0}+\frac{k}{2}(x-x_{0})^{2}+\frac{k_{3}}{3}(x-x_{0})^{3}+\frac{k_{4}}{4}(x-x_{0})^{4}$
(S12) $\displaystyle k_{4}=$ $\displaystyle 4e$ (S13) $\displaystyle k_{3}=$
$\displaystyle 3(d+k_{4}x_{0})$ (S14) $\displaystyle k=$ $\displaystyle
2(c+k_{3}x_{0}-\frac{6}{4}k_{4}x_{0}^{2})$ (S15) $\displaystyle k_{0}=$
$\displaystyle
a-\frac{k_{2}}{2}x_{0}^{2}+\frac{k_{3}}{3}x_{0}^{3}-\frac{k_{4}}{4}x_{0}^{4}.$
(S16)

Note that $k_{0}$ is small but can be non-zero because of discretization
effects. The results of the fits for the different vibrational coordinates are
presented in fig. S3A–D, by comparing the second derivatives of the fit
functions, $U^{\prime\prime}(x)=k+2k_{3}(x-x_{0})+3k_{4}(x-x_{0})^{2}$, (blue
broken lines) with second derivatives obtained numerically from the data,
$U^{\prime\prime}(x_{i})=(U(x_{i+1})-2U(x_{i+1})+U(x_{i-1}))/(x_{i}-x_{i-1})^{2}$,
using histograms with 300 equidistant bins. The numerically obtained derivates
are smoothed (grey and purple dots) by iterative convolution with a flat
window function, which has a width of three bins. The plots illustrate that
for all vibrational coordinates the fit functions and particularly the
harmonic fit parameter $k$ (shown as horizontal grey dotted lines) match well
the numerically obtained curvature of the potential around the position of the
minima at $x_{0}$ (shown as vertical grey dotted lines).

\begin{overpic}[width=212.47617pt]{{./f/dft_nvt_256_f_ohs_2ndDFit}.pdf}
\put(2.0,55.0){\large\bf A}
\end{overpic}\begin{overpic}[width=212.47617pt]{{./f/dft_nvt_256_f_hoh_2ndDFit}.pdf}
\put(2.0,55.0){\large\bf B}
\end{overpic}\begin{overpic}[width=212.47617pt]{{./f/dft_nvt_256_f_oha_2ndDFit}.pdf}
\put(2.0,55.0){\large\bf C}
\end{overpic}\begin{overpic}[width=212.47617pt]{{./f/dft_nvt_256_f_oh_2ndDFit}.pdf}
\put(2.0,55.0){\large\bf D} \end{overpic}

Figure S3: Second derivatives of the potential, obtained numerically and
iteratively smoothed (grey and purple dots), of the various vibrational
coordinates (A: $d_{\mathrm{OHs}}$ coordinate, B: $\phi_{\mathrm{HOH}}$
coordinate, C: $d_{\mathrm{OHa}}$ coordinate and D: $d_{\mathrm{OH}}$
coordinate.). The second derivatives of the fit functions are shown as blue
broken lines. The vertical grey dotted lines denote the position of the minima
at $x_{0}$ and the horizontal grey dotted lines denote the values of the
harmonic fit parameters $k$. The insets show a wider range of data, for which
$U(x_{i})<8\,k_{B}T$.

##  V Extraction of the time-dependent friction kernels

The extraction of time-dependent friction kernels from trajectories with an
anharmonic system potential is performed using a modification of a recently
published approach [3] and based on previous work by Harp and Berne [4]. The
derivation starts from the generalized Langevin equation (GLE)

$\displaystyle
m\ddot{x}(t)=-\int_{0}^{t}\Gamma(t-t^{\prime})\dot{x}(t^{\prime})dt^{\prime}-\nabla
U[x(t)]+F_{R}(t),$ (S17)

which is multiplied by the initial velocity $\dot{x}(0)$ and ensemble averaged

$\displaystyle\begin{split}m\left<\dot{x}(0)\ddot{x}(t)\right>=-\int_{0}^{t}dt^{\prime}\,\Gamma(t^{\prime})&\left<\dot{x}(0)\dot{x}(t-t^{\prime})\right>-\left<\dot{x}(0)\nabla
U[x(t)]\right>,\end{split}$ (S18)

assuming the random force $F_{R}(t)$ to be uncorrelated with $\dot{x}(0)$.
Defining the correlation functions as

$\displaystyle C_{vv}(t)$ $\displaystyle=\left<\dot{x}(0)\dot{x}(t)\right>,$
(S19) $\displaystyle C_{v\nabla U}(t)$ $\displaystyle=\left<\dot{x}(0)\nabla
U[x(t)]\right>,$ (S20) $\displaystyle C_{x\nabla U}(t)$
$\displaystyle=\left<x(0)\nabla U[x(t)]\right>,$ (S21)

eq. (S18) is written as

$\displaystyle
m\frac{d}{dt}C_{vv}(t)=-\int_{0}^{t}dt^{\prime}\,\Gamma(t^{\prime})$
$\displaystyle C_{vv}(t-t^{\prime})-C_{v\nabla U}(t),$ (S22)

and integrated over time

$\displaystyle
mC_{vv}(t)-mC_{vv}(0)=-\int_{0}^{t}dt^{\prime\prime}\int_{0}^{t^{\prime\prime}}dt^{\prime}\,\Gamma(t^{\prime})$
$\displaystyle
C_{vv}(t^{\prime\prime}-t^{\prime})-\int_{0}^{t}dt^{\prime\prime}\ C_{v\nabla
U}(t^{\prime\prime}).$ (S23)

For the terms on the right hand side, one finds

$\displaystyle\int_{0}^{t}dt^{\prime\prime}\,C_{v\nabla U}(t^{\prime\prime})$
$\displaystyle=\int_{0}^{t}dt^{\prime\prime}\,\left<\dot{x}(0)\nabla
U[x(t^{\prime\prime})]\right>$
$\displaystyle=\int_{0}^{t}dt^{\prime\prime}\,\left<\dot{x}(-t^{\prime\prime})\nabla
U[x(0)]\right>$ $\displaystyle=\left<\nabla
U[x(0)]\int_{0}^{t}dt^{\prime\prime}\,\dot{x}(-t^{\prime\prime})\right>$
$\displaystyle=\left<\nabla
U[x(0)]\int_{-t}^{0}dt^{\prime\prime}\,\dot{x}(t^{\prime\prime})\right>$
$\displaystyle=\left<\nabla U[x(0)]x(0)\right>-\left<\nabla
U[x(0)]x(-t)\right>$ $\displaystyle=C_{x\nabla U}(0)-C_{x\nabla U}(t),$ (S24)

and

$\displaystyle\int_{0}^{t}dt^{\prime\prime}\int_{0}^{t^{\prime\prime}}dt^{\prime}\,\Gamma(t^{\prime})C_{vv}(t^{\prime\prime}-t^{\prime})$
$\displaystyle=\int_{0}^{t}dt^{\prime\prime}\int_{0}^{t^{\prime\prime}}dt^{\prime\prime\prime}\,\Gamma(t^{\prime\prime}-t^{\prime\prime\prime})C_{vv}(t^{\prime\prime\prime})$
$\displaystyle=\int_{0}^{t}dt^{\prime\prime}\int_{0}^{t^{\prime\prime}}dt^{\prime\prime\prime}\,\Gamma(t^{\prime\prime}-t^{\prime\prime\prime})C_{vv}(t^{\prime\prime\prime})$
$\displaystyle=\int_{0}^{t}dt^{\prime\prime\prime}\,C_{vv}(t^{\prime\prime\prime})\int_{t^{\prime\prime\prime}}^{t}dt^{\prime\prime}\,\Gamma(t^{\prime\prime}-t^{\prime\prime\prime})$
$\displaystyle=\int_{0}^{t}dt^{\prime\prime\prime}\,C_{vv}(t^{\prime\prime\prime})\int_{0}^{t-t^{\prime\prime\prime}}dt^{\prime\prime}\,\Gamma(t^{\prime\prime})$
$\displaystyle=\int_{0}^{t}dt^{\prime\prime\prime}\,C_{vv}(t^{\prime\prime\prime})G(t-t^{\prime\prime\prime}),$
(S25)

where

$\displaystyle G(t)=\int_{0}^{t}dt^{\prime}\,\Gamma(t^{\prime}),$ (S26)

is the integral over the friction kernel. Inserting eq. (S24) and eq. (S25)
into eq. (S23), one obtains

$\displaystyle
mC_{vv}(t)-mC_{vv}(0)=-\int_{0}^{t}dt^{\prime}\,C_{vv}(t^{\prime})G(t-t^{\prime})-C_{x\nabla
U}(0)+C_{x\nabla U}(t).$ (S27)

The mass is found to be

$\displaystyle m=\frac{C_{x\nabla U}(0)}{C_{vv}(0)},$ (S28)

which follows from multiplying the GLE eq. (S17) by the initial position
$x(0)$, averaging over the ensemble and evaluating at $t=0$

$\displaystyle m\left<x(0)\ddot{x}(0)\right>$ $\displaystyle=-\left<x(0)\nabla
U[x(0)]\right>$ (S29) $\displaystyle\quad-m\left<\dot{x}(0)\dot{x}(0)\right>$
$\displaystyle=-\left<x(0)\nabla U[x(0)]\right>.$ (S30)

Using eq. (S28) in eq. (S27), one finds

$\displaystyle C_{vv}(t)\frac{C_{x\nabla
U}(0)}{C_{vv}(0)}=-\int_{0}^{t}dt^{\prime}\,C_{vv}(t^{\prime})G(t-t^{\prime})+C_{x\nabla
U}(t).$ (S31)

Eq. (S31) is a Volterra equation of first kind which can be discretized in
time, $t=i\Delta t$, and solved numerically,

$\displaystyle C_{vv}^{i}\frac{C_{x\nabla U}^{0}}{C_{vv}^{0}}$
$\displaystyle=-\sum_{j=0}^{i}w_{j,i}\Delta t\,C_{vv}^{i}G^{i-j}+C_{x\nabla
U}^{i}$ (S32) $\displaystyle=-G^{i}w_{i,i}\Delta
t\,C_{vv}^{0}-\sum_{j=0}^{i-1}w_{j,i}\Delta t\,C_{vv}^{i}G^{i-j}+C_{x\nabla
U}^{i},$ (S33)

where $w_{j,i}$ are integration weights and the corresponding iteration
relation reads

$\displaystyle G^{i}=-\frac{1}{\omega_{i,i}\Delta
tC_{vv}^{0}}\left(\sum_{j=0}^{i-1}\omega_{i,j}\Delta
tC_{vv}^{i}G^{i-j}+\frac{C_{x\nabla U}^{0}}{C_{vv}^{0}}C_{vv}^{i}-C_{x\nabla
U}^{i}\right).$ (S34)

Thus by calculating the necessary correlation functions,
$C_{vv}^{i}=C_{vv}(i\Delta t)$ and $C_{x\nabla U}^{i}=C_{x\nabla U}(i\Delta
t)$, from equilibrium trajectories, the integral of the time-dependent
friction kernel, $G^{i}=G(i\Delta t)$, can be obtained by iteration from the
initial value $G^{0}=0$. The friction kernel $\Gamma^{i}=\Gamma(i\Delta t)$ is
subsequently obtained by a numerical derivative.

##  VI Power spectrum of the Generalized Langevin Equation

A random process $x(t)$ with mass $m$ and subject to time-dependent (or
frequency-dependent) friction $\Gamma(t)$ in a harmonic potential
$U(x)=\frac{k}{2}x^{2}$ can be described by the generalized Langevin equation

$\displaystyle
m\ddot{x}(t)=-\int_{0}^{t}dt^{\prime}\,\Gamma(t-t^{\prime})\dot{x}(t^{\prime})-kx(t)+F_{R}(t),$
(S35)

obeying for the random force $F_{R}(t)$ the fluctuation-dissipation relation
$\langle F_{R}(t)F_{R}(t^{\prime})\rangle=k_{B}T\Gamma(t-t^{\prime})$. A
Fourier transform gives

$\displaystyle-\omega^{2}m\tilde{x}(\omega)=i\omega\widetilde{\Gamma}^{+}(\omega)\tilde{x}(\omega)-k\tilde{x}(\omega)+\widetilde{F}_{R}(\omega),$
(S36)

as well as

$\displaystyle\langle\widetilde{F}_{R}(\omega)\widetilde{F}_{R}(\omega^{\prime})\rangle$
$\displaystyle=k_{B}T\int_{-\infty}^{\infty}dt\ e^{i\omega
t}\int_{-\infty}^{\infty}dt^{\prime}\
e^{i\omega^{\prime}t^{\prime}}\Gamma(t-t^{\prime})$ (S37)
$\displaystyle=k_{B}T\int_{-\infty}^{\infty}dt\ e^{i\omega
t}\int_{-\infty}^{\infty}dt^{\prime}\
e^{i\omega^{\prime}t^{\prime}}\Gamma(t-t^{\prime})$
$\displaystyle=k_{B}T\int_{-\infty}^{\infty}dt^{\prime}\
e^{i(\omega^{\prime}+\omega)t^{\prime}}\int_{-\infty}^{\infty}dt\
e^{i\omega(t-t^{\prime})}\Gamma(t-t^{\prime})$ $\displaystyle=2\pi
k_{B}T\delta(\omega^{\prime}+\omega)\widetilde{\Gamma}(\omega)$ (S38)

The absorbed power $\omega\widetilde{\chi}^{\prime\prime}(\omega)$ is derived
in Fourier space from linear response to an external force $F_{\rm ext}$,
which is calculated according to

$\displaystyle\widetilde{\chi}(\omega)$
$\displaystyle=\frac{\tilde{x}(\omega)}{\widetilde{F}_{\mathrm{ext}}(\omega)}$
(S39)
$\displaystyle=(k-m\omega^{2}-i\widetilde{\Gamma}^{+}(\omega)\omega)^{-1},$
(S40)

and using eqs. (S38) and (S40) the autocorrelation in Fourier space is given
by

$\displaystyle\widetilde{C}_{xx}(\omega)$
$\displaystyle=\int_{-\infty}^{\infty}dt\ e^{i\omega t}\langle
x(t)x(0)\rangle$ (S41)
$\displaystyle=\frac{1}{4\pi^{2}}\int_{-\infty}^{\infty}dt\ e^{i\omega
t}\int_{-\infty}^{\infty}d\omega^{\prime}\
e^{-i\omega^{\prime}t}\int_{-\infty}^{\infty}d\omega^{\prime\prime}e^{-i\omega^{\prime\prime}0}\langle\tilde{x}(\omega^{\prime})\tilde{x}(\omega^{\prime\prime})\rangle$
$\displaystyle=\frac{1}{4\pi^{2}}\int_{-\infty}^{\infty}d\omega^{\prime}\
2\pi\delta(\omega-\omega^{\prime})\int_{-\infty}^{\infty}d\omega^{\prime\prime}\langle\tilde{x}(\omega^{\prime})\tilde{x}(\omega^{\prime\prime})\rangle$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}d\omega^{\prime\prime}\langle\tilde{x}(\omega)\tilde{x}(\omega^{\prime\prime})\rangle$
(S42)
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}d\omega^{\prime\prime}\langle\widetilde{\chi}(\omega)\widetilde{F}_{R}(\omega)\widetilde{\chi}(\omega^{\prime\prime})\widetilde{F}_{R}(\omega^{\prime\prime})\rangle$
$\displaystyle=k_{B}T\int_{-\infty}^{\infty}d\omega^{\prime\prime}\widetilde{\chi}(\omega)\widetilde{\chi}(\omega^{\prime\prime})\delta(\omega^{\prime\prime}+\omega)\widetilde{\Gamma}(\omega)$
$\displaystyle=k_{B}T\widetilde{\chi}(\omega)\widetilde{\chi}(-\omega)\widetilde{\Gamma}(\omega),$
(S43)

which is further simplified by realizing that the Fourier transform of the
purely real function $\chi(t)$ is even
$\widetilde{\chi}(-\omega)=\widetilde{\chi}^{*}(\omega)$

$\displaystyle\widetilde{C}_{xx}(\omega)$
$\displaystyle=k_{B}T\widetilde{\chi}(\omega)\widetilde{\chi}^{*}(\omega)\widetilde{\Gamma}(\omega)$
$\displaystyle=k_{B}T\left|\widetilde{\chi}(\omega)\right|^{2}\widetilde{\Gamma}(\omega)$
$\displaystyle=\frac{k_{B}T\widetilde{\Gamma}(\omega)}{\left|k-m\omega^{2}-i\widetilde{\Gamma}^{+}(\omega)\omega\right|^{2}}.$
(S44)

From eq. (S4) the power spectrum follows

$\displaystyle\omega\widetilde{\chi}^{\prime\prime}(\omega)=\frac{\omega^{2}\
\text{Re}\widetilde{\Gamma}^{+}(\omega)}{\left|k-m\omega^{2}-i\widetilde{\Gamma}^{+}(\omega)\omega\right|^{2}}.$
(S45)

##  VII Power spectrum of the damped harmonic oscillator

The absorbed power $\omega\widetilde{\chi}^{\prime\prime}(\omega)$ of the
damped harmonic oscillator described by the differential equation

$\displaystyle m\ddot{x}(t)=-\gamma\dot{x}(t)-kx(t)+F_{\mathrm{ext}}(t),$
(S46)

is computed from the linear response in Fourier space

$\displaystyle\widetilde{\chi}(\omega)$
$\displaystyle=\frac{\tilde{x}(\omega)}{\tilde{F}_{\mathrm{ext}}(\omega)}$
(S47) $\displaystyle=(k-m\omega^{2}-i\gamma\omega)^{-1}$ (S48)
$\displaystyle=\frac{k-m\omega^{2}+i\gamma\omega}{(k-m\omega^{2})^{2}+\gamma^{2}\omega^{2}},$
(S49)

where $\tilde{x}(\omega)$ is the oscillating variable, $m$ is the mass,
$\gamma$ the friction coefficient, $k$ the spring constant of the harmonic
potential and $\tilde{F}_{\mathrm{ext}}(\omega)$ an external force. For the
power spectrum follows

$\displaystyle\omega\widetilde{\chi}^{\prime\prime}(\omega)=\frac{\gamma\omega^{2}}{(k-m\omega^{2})^{2}+\gamma^{2}\omega^{2}},$
(S50)

which by introducing the time scales $\tau=2\gamma/k$,
$\tau_{\omega}=\sqrt{m/k}$ and length scale $L$ with $L^{2}=2k_{B}T/k$
converts to

$\displaystyle\omega\widetilde{\chi}^{\prime\prime}(\omega)=\frac{L^{2}}{k_{B}T}\frac{\tau\omega^{2}}{4(1-\tau_{\omega}^{2}\omega^{2})^{2}+\tau^{2}\omega^{2}}.$
(S51)

In spectroscopy this is known as a Lorentzian line shape [5], which in the
overdamped case, $\tau_{\omega}\to 0$, reads

$\displaystyle\omega\widetilde{\chi}^{\prime\prime}(\omega)$
$\displaystyle=\frac{\gamma\omega^{2}}{k^{2}+\gamma^{2}\omega^{2}}$ (S52)
$\displaystyle=\frac{L^{2}}{k_{B}T}\frac{\tau\omega^{2}}{4+\tau^{2}\omega^{2}}.$
(S53)

Eq. (S52) is also known as the Debye line shape.

##  VIII Parametrization of the extracted memory kernels

The memory kernels $\Gamma(t)$, extracted from the aiMD simulation data as
described in section V, are truncated at $50\text{\,}\mathrm{p}\mathrm{s}$,
over-sampled at a time resolution of $\Delta
t=$0.025\text{\,}\mathrm{f}\mathrm{s}$$ and subsequently Fourier transformed
using the FFT algorithm implemented in numpy v1.19. The real part of the
Fourier-transformed memory kernel is then fitted to a combination of $n$
exponential and $l$ oscillating memory kernels according to eq.
(LABEL:eq:kernel_param) in the main text. The Fourier-transformed expressions
of the fundamental kernels are found to be

$\displaystyle\widetilde{\Gamma}_{\rm
osc}^{+}(\omega,a_{i},\tau_{i},\omega_{i})$
$\displaystyle=\frac{a_{i}}{2}\Bigl{(}\frac{1+\frac{i}{\omega_{i}\tau_{i}}}{\frac{1}{\tau_{i}}+i\omega_{i}-i\omega}+\frac{1-\frac{i}{\omega_{i}\tau_{i}}}{\frac{1}{\tau_{i}}-i\omega_{i}-i\omega}\Bigr{)}$
(S54) $\displaystyle\widetilde{\Gamma}_{\rm
exp}^{+}(\omega,\gamma_{i},\tau_{i})$
$\displaystyle=\gamma_{i}\dfrac{1}{1-i\omega\tau_{i}}$ (S55)
$\displaystyle\tilde{\Gamma}^{+}(\omega)$
$\displaystyle=\sum\limits_{i=1}^{n}\widetilde{\Gamma}_{exp}^{+}(\omega,\gamma_{i},\tau^{e}_{i})+\sum\limits_{i=1}^{l}\widetilde{\Gamma}_{osc}^{+}(\omega,a_{i},\tau^{o}_{i},\omega_{i})$
(S56)

The real part of the Fourier-transformed memory kernel
$\widetilde{\Gamma}^{\prime+}(\omega)$ is fitted in the $(2n+3l)$-dimensional
parameter space using the Levenberg-Marquardt algorithm implemented in scipy
v1.5. $n$ and $l$ are iteratively increased until the fit quality does not
improve significantly. The initial values for all $\gamma_{i}$,
$\tau^{e}_{i}$, $a_{i}$, $\tau^{o}_{i}$ and $\omega_{i}$ are chosen suitably.
After fitting in Fourier space, another fit is performed in the time domain,
by subtracting the $l$ oscillating fit functions from the memory kernels in
the time domain and again fitting the $n$ exponential functions to the
remainder. This significantly improves the fits for the long-time tails of the
memory kernels. The hereby obtained fit parameters for the friction memory
kernels for the different vibrational coordinates in the bulk water system are
summarized in tabs. VIII–VIII.

Table S1: Fit parameters according to eq. (LABEL:eq:kernel_param) in the main
text of the friction memory kernel of the $d_{\mathrm{OHs}}$ coordinate.

$\gamma_{i}$ [u/ps] |  $\tau_{i}^{e}$ [ps] |
---|---|---
2758 | 0.115 |
25941 | 1.28 |
54005 | 5.72 |
$a_{i}$ [u/ps2] |  $\tau_{i}^{o}$ [ps] |  $\omega_{i}$ [THz]
8468 | 0.00944 | 170
3525 | 0.0146 | 1290
5985 | 0.0256 | 614
1095 | 0.0599 | 675
10560 | 0.0807 | 28.7
6318 | 0.122 | 313
Table S2: Fit parameters according to eq. (LABEL:eq:kernel_param) in the main
text of the friction memory kernel of the $\phi_{\mathrm{HOH}}$ coordinate.
Table S3: Fit parameters according to eq. (LABEL:eq:kernel_param) in the main
text of the friction memory kernel of the $d_{\mathrm{OHa}}$ coordinate. Table
S4: Fit parameters according to eq. (LABEL:eq:kernel_param) in the main text
of the friction memory kernel of the $d_{\mathrm{OH}}$ coordinate.

##  IX Simulation of the generalized Langevin equation

The GLE

$\displaystyle
m\ddot{x}(t)=-\int_{0}^{t}\Gamma(t-t^{\prime})\dot{x}(t^{\prime})dt^{\prime}-\nabla
U[x(t)]+F_{R}(t),$ (S57)

with a sum of $n$ exponential and $l$ oscillating memory kernels, analogous to
eq. (LABEL:eq:kernel_param) in the main text,

$\displaystyle\Gamma(t)$ $\displaystyle=\sum\limits_{i=1}^{n}a_{i}\
e^{-t/\tau_{i}}+\sum\limits_{i=1}^{l}a_{i}\
e^{-t/\tau_{i}}\Bigl{(}\cos(\omega_{i}t)+\dfrac{1}{\tau_{i}\omega_{i}}\sin(\omega_{i}t)\Bigr{)},$
(S58)

can be efficiently simulated using a Markovian embedding (which for the
oscillating components is derived in detail in section X)

$\displaystyle\dot{x}(t)$ $\displaystyle=v(t)$ (S59) $\displaystyle
m_{x}\dot{v}(t)$ $\displaystyle=-\nabla
U[x(t)]+\sum\limits_{i=1}^{n}a_{i}[y_{i}(t)-x(t)]+\sum\limits_{i=1}^{l}a_{i}[z_{i}(t)-x(t)]$
(S60) $\displaystyle\gamma_{i}\dot{y}_{i}(t)$
$\displaystyle=a_{i}[x(t)-y_{i}(t)]+F_{i}(t)$ for $n$ exp. components (S61)
$\displaystyle\dot{z}_{i}(t)$ $\displaystyle=w_{i}(t)$ (S62) $\displaystyle
m_{i}\dot{w}_{i}(t)$
$\displaystyle=-\gamma_{i}w_{i}(t)+a_{i}[x(t)-z_{i}(t)]+F_{i}(t)$ for $l$ osc.
components (S63)

For the $n$ exponential components we obtain
$a_{i}=\frac{\gamma_{i}}{\tau_{i}^{e}}$ from the fit parameters $\gamma_{i}$
and $\tau_{i}^{e}$. To get $m_{i}$ and $\gamma_{i}$ for the $l$ oscillating
components from the fit parameters $a_{i}$, $\tau^{o}_{i}$ and $\omega_{i}$ we
use

$\displaystyle m_{i}$
$\displaystyle=a_{i}((\tau^{o}_{i})^{-2}+\omega_{i}^{2})^{-1}$ (S64)
$\displaystyle\gamma_{i}$ $\displaystyle=2\frac{m_{i}}{\tau^{o}_{i}}$ (S65)

The random force $F_{i}$ present in eqs. (S61) and (S63) is

$\displaystyle F_{i}$ $\displaystyle=\sqrt{2k_{B}T\gamma_{i}\delta t^{-1}}\Xi$
(S66)

Where $\Xi$ is a Gaussian random distribution with zero mean and a standard
deviation of one. In the numerical simulation we use the numpy.random.normal
function to generate these random numbers. Eqs. (S60)–(S63) are numerically
solved using a 4th-order Runge-Kutta scheme to get the trajectory of $x$. The
timestep of the simulation is $1\text{\times}{10}^{-6}\text{\,}\mathrm{ps}$
and only every 500th value is stored. Thus the trajectory is evaluated at a
time step of $0.5\text{\,}\mathrm{f}\mathrm{s}$. The trajectory contains
$2.5\text{\times}{10}^{7}$ values which corresponds to a simulation time of
$1.25\text{\times}{10}^{4}\text{\,}\mathrm{ps}$. Each trajectory is divided
into 50 parts of $250\text{\,}\mathrm{ps}$ and their spectra are calculated
using the Wiener-Kintchine relation introduced in section XIV. All 50 spectra
are then averaged to give the final spectrum. Starting velocity and position
were set to be zero.

The memoryless Langevin equation, with a friction constant $\gamma$ replacing
the integral over the memory kernel

$\displaystyle m\ddot{x}(t)=-\gamma\dot{x}(t^{\prime})-\nabla
U[x(t)]+F_{R}(t),$ (S67)

is likewise simulated using a 4th-order Runge-Kutta scheme. The timestep of
these simulations is $0.1\text{\,}\mathrm{f}\mathrm{s}$ for a total of
$10^{7}$ steps, which corresponds to a simulation time of
$1000\text{\,}\mathrm{ps}$. Each trajectory is divided into 200 parts of
$5\text{\,}\mathrm{ps}$ and their spectra are calculated using the Wiener-
Kintchine relation introduced in section XIV. All 200 spectra are then
averaged to give the final spectrum.

##  X Oscillating memory kernels

###  X.1 Oscillating memory kernels from Hamiltonian

Oscillating memory kernels arise from inertial dynamics in orthogonal degrees
of freedom at time scales of the order of the primary reaction coordinate. The
relation of the time scales of the orthogonal coordinate and the memory kernel
is understood from an explicit derivation of the dynamics of the primary
coordinate $x$ with velocity $v$, linearly coupled to an orthogonal coordinate
$y$ with velocity $w$ and coupling constant $k$. The Hamiltonian of the system
is given as

$\displaystyle H=\frac{m}{2}v^{2}+\frac{m_{y}}{2}w^{2}+\frac{k}{2}(x-y)^{2},$
(S68)

defining the Hamilton equations

$\displaystyle\dot{x}(t)$ $\displaystyle=v$ (S69) $\displaystyle m\dot{v}(t)$
$\displaystyle=k\left[x(t)-y(t)\right]$ (S70) $\displaystyle\dot{y}(t)$
$\displaystyle=w(t)$ (S71) $\displaystyle m_{y}\dot{w}(t)$
$\displaystyle=-k\left[x(t)-y(t)\right].$ (S72)

These Newtonian equations of motion can straight-forwardly be extended to
Langevin equations by introducing random forces $F_{x}(t),F_{y}(t)$ and
friction terms with coefficients $\gamma_{x},\gamma_{y}$

$\displaystyle\dot{x}(t)$ $\displaystyle=v$ (S73) $\displaystyle m\dot{v}(t)$
$\displaystyle=-\gamma_{x}v(t)+k\left[x(t)-y(t)\right]+F_{x}(t)$ (S74)
$\displaystyle\dot{y}(t)$ $\displaystyle=w(t)$ (S75) $\displaystyle
m_{y}\dot{w}(t)$
$\displaystyle=-\gamma_{y}w(t)-k\left[x(t)-y(t)\right]+F_{y}(t).$ (S76)

The random forces are designed to introduce dynamics on the order of $k_{B}T$
but no net energy transfer into the system. Therefore they have zero mean,
$\langle F_{i}(t)\rangle=0$, and their strength on the order of the friction
forces $\langle
F_{i}(t)F_{i}(t^{\prime})\rangle=2\gamma_{i}k_{B}T\delta(t-t^{\prime})$,
fulfilling the fluctuation-dissipation relation.

To obtain the GLE including the memory kernel we will now find a solution for
$y$ using eq. (S76) and insert this solution in eq. (S74). Eq. (S76) is
essentially a driven damped harmonic oscillator, which can be solved in a
general fashion using matrix exponentials and introducing the inertial time
scale $\tau_{m_{y}}=m_{y}/\gamma_{y}$ and the memory time scale
$\tau_{k_{y}}=\gamma_{y}/k$. First we write the coupled first order ordinary
differential eqs. (S75) and (S76) in the matrix form

$\displaystyle\dot{\vec{y}}(t)$ $\displaystyle=A\vec{y}(t),$ (S77)

where

$\displaystyle\vec{y}(t)$ $\displaystyle=\begin{pmatrix}y\\\
\dot{y}\end{pmatrix}$ (S78) $\displaystyle A$
$\displaystyle=\begin{pmatrix}0&1\\\
-(\tau_{m_{y}}\tau_{k_{y}})^{-1}&-\tau_{m_{y}}^{-1}\end{pmatrix},$ (S79)
$\displaystyle A^{-1}$
$\displaystyle=\begin{pmatrix}-\tau_{m_{y}}^{-1}&(\tau_{m_{y}}\tau_{k_{y}})^{-1}\\\
-1&0\end{pmatrix}.$ (S80)

The solution is given as

$\displaystyle\begin{split}\vec{y}(t)&=\exp{\left(A[t-t_{0}]\right)}\vec{y}_{0}\\\
&+\int_{t_{0}}^{t}dt^{\prime}\exp{\left(A[t-t^{\prime}]\right)}\begin{pmatrix}0\\\
-(\tau_{m_{y}}\tau_{k_{y}})^{-1}x(t^{\prime})+m_{y}^{-1}F_{y}(t^{\prime})\end{pmatrix},\end{split}$
(S81)

Next a partial integration over the integral containing $x(t)$ is performed

$\displaystyle\begin{split}\vec{y}(t)&=\exp{\left(A[t-t_{0}]\right)}\vec{y}_{0}\\\
&-\left[A^{-1}\exp{\left(A[t-t^{\prime}]\right)}\begin{pmatrix}0\\\
-(\tau_{m_{y}}\tau_{k_{y}})^{-1}x(t)\end{pmatrix}\right]_{t_{0}}^{t}\\\
&+A^{-1}\int_{t_{0}}^{t}dt^{\prime}\exp{\left(A[t-t^{\prime}]\right)}\begin{pmatrix}0\\\
-(\tau_{m_{y}}\tau_{k_{y}})^{-1}v(t^{\prime})\end{pmatrix}\\\
&+\int_{t_{0}}^{t}dt^{\prime}\exp{\left(A[t-t^{\prime}]\right)}\begin{pmatrix}0\\\
m_{y}^{-1}F_{y}(t^{\prime})\end{pmatrix}.\end{split}$ (S82)

The matrix exponential can be evaluated by diagonalizing the matrix $A$. After
introducing
$\omega_{0}=\sqrt{(2\tau_{m_{y}})^{-2}-(\tau_{m_{y}}\tau_{k_{y}})^{-1}}$ the
Eigenvalues $\lambda_{1,2}$ of $A$ are given as

$\displaystyle\lambda_{1,2}=-(2\tau_{m_{y}})^{-1}\pm\omega_{0},$ (S83)

and the matrix exponential reads

$\displaystyle\begin{split}\exp(At)&=\frac{1}{\lambda_{2}-\lambda_{1}}\begin{pmatrix}\lambda_{2}e^{\lambda_{1}t}-\lambda_{1}e^{\lambda_{2}t}&e^{\lambda_{2}t}-e^{\lambda_{1}t}\\\
\lambda_{1}\lambda_{2}\left(e^{\lambda_{1}t}-e^{\lambda_{2}t}\right)&\lambda_{2}e^{\lambda_{2}t}-\lambda_{1}e^{\lambda_{1}t}\end{pmatrix}\end{split}$
(S84)
$\displaystyle\begin{split}&=e^{-t/(2\tau_{m_{y}})}\begin{pmatrix}\cosh(\omega_{0}t)+\frac{\sinh(\omega_{0}t)}{2\tau_{m_{y}}\omega_{0}}&-\frac{\sinh(\omega_{0}t)}{\omega_{0}}\\\
\frac{\sinh(\omega_{0}t)}{\tau_{m_{y}}\tau_{k_{y}}\omega_{0}}&\cosh(\omega_{0}t)-\frac{\sinh(\omega_{0}t)}{2\tau_{m_{y}}\omega_{0}}\end{pmatrix},\end{split}$
(S85)

giving

$\displaystyle\begin{split}&A^{-1}\exp(At)=\tau_{m_{y}}\tau_{k_{y}}e^{-t/(2\tau_{m_{y}})}\\\
&\begin{pmatrix}-\frac{1}{\tau_{m_{y}}}\cosh(\omega_{0}t)+\frac{1}{2\tau_{m_{y}}^{2}\omega_{0}}\sinh(\omega_{0}t)+\frac{1}{\tau_{m_{y}}\tau_{k_{y}}\omega_{0}}\sinh(\omega_{0}t)&-\frac{1}{2\tau_{m_{y}}\omega_{0}}\sinh(\omega_{0}t)-\cosh(\omega_{0}t)\\\
\frac{1}{\tau_{m_{y}}\tau_{k_{y}}}\cosh(\omega_{0}t)+\frac{1}{\tau_{m_{y}}\tau_{k_{y}}\omega_{0}}\sinh(\omega_{0}t)&\frac{1}{\tau_{m_{y}}\tau_{k_{y}}\omega_{0}}\sinh(\omega_{0}t)\end{pmatrix}.\end{split}$
(S86)

The solution for $y(t)=\vec{y}_{1}(t)$ is now obtained from eqs. (S82)–(S86)
and by sending $t_{0}\to-\infty$, i.e. assuming an equilibrated system

$\displaystyle\begin{split}y(t)&=x(t)-\int_{-\infty}^{t}dt^{\prime}v(t^{\prime})e^{-(t-t^{\prime})/(2\tau_{m_{y}})}\\\
&\left[\cosh(\omega_{0}(t-t^{\prime}))+\frac{1}{2\tau_{m_{y}}\omega_{0}}\sinh(\omega_{0}(t-t^{\prime}))\right]\\\
&-\int_{-\infty}^{t}dt^{\prime}m_{y}^{-1}F_{y}(t^{\prime})e^{-(t-t^{\prime})/(2\tau_{m_{y}})}\frac{1}{\omega_{0}}\sinh(\omega_{0}(t-t^{\prime})).\end{split}$
(S87)

To obtain the generalized Langevin equation for $x(t)$ containing the
oscillating memory kernel

$\displaystyle\Gamma(t)=ke^{-t/(2\tau_{m_{y}})}\left[\cosh(\omega_{0}t)+\frac{1}{2\tau_{m_{y}}\omega_{0}}\sinh(\omega_{0}t)\right],$
(S88)

eq. (S87) is inserted in eq. (S74)

$\displaystyle\begin{split}m\dot{v}(t)&=-\gamma_{x}v(t)+F_{x}(t)+\int_{-\infty}^{t}dt^{\prime}v(t^{\prime})\Gamma(t-t^{\prime})+F^{\prime}_{y}(t).\end{split}$
(S89)

###  X.2 Underdamped limit: oscillating kernel

The kernel of eq. (S88) can be analysed for two limiting cases. For
Im$(\omega_{0})>(2\tau_{m_{y}})^{-1}$, which is equivalent to
$4\tau_{m_{y}}>\tau_{k_{y}}$, the underdamped case, $\omega_{0}$ is purely
imaginary and the kernel is

$\displaystyle\Gamma(t)=ke^{-t/(2\tau_{m_{y}})}\cos(\mathrm{Im}(\omega_{0})t)+\frac{1}{2\tau_{m_{y}}\mathrm{Im}(\omega_{0})}\sin(\mathrm{Im}(\omega_{0})t).$
(S90)

The integral is evaluated to

$\displaystyle\gamma=\int_{0}^{\infty}dt\Gamma(t)=k\frac{4\tau_{m_{y}}}{1-4\tau_{m_{y}}^{2}\mathrm{Im}(\omega_{0})^{2}}=k\tau_{k_{y}}=\gamma_{y},$
(S91)

###  X.3 Overdamped limit: exponential kernel

For Re$(\omega_{0})<(4\tau_{m_{y}})^{-1}$ and $\tau_{k_{y}}>(4\tau_{m_{y}})$,
the overdamped case, it follows that $\omega_{0}$ is purely real. The integral
equally gives

$\displaystyle\gamma=\int_{0}^{\infty}dt\Gamma(t)=k\frac{4\tau_{m_{y}}}{1-4\tau_{m_{y}}^{2}\omega_{0}^{2}}=k\tau_{k_{y}}=\gamma_{y}.$
(S92)

##  XI Decomposition of the OH stretch mode

\begin{overpic}[width=216.81pt]{{spec_ohDecomp_dft_nvt_256}.png} \end{overpic}
Figure S4: Comparison of IR spectra in the OH-stretching regime obtained from
the total dipole-moment trajectory including nuclear and electronic charges
from the aiMD simulation of 256 H2O molecules (grey solid line), an
approximation using partial charges on the nuclear dynamics of the same aiMD
simulation (blue sold line, introduced in section III), experimental FTIR data
(grey broken line) [6] and rescaled power spectra of the $d_{\mathrm{OH}}$
(green solid line), $d_{\mathrm{OHs}}$ (purple solid line) and
$d_{\mathrm{OHa}}$ (purple broken line) vibrational coordinates.

As introduced in the main text, vibrational coordinates for the OH stretch
mode are given by
$d_{\mathrm{OHs}}=({d_{\mathrm{OH}}}_{1}+{d_{\mathrm{OH}}}_{2})/2$ for the
symmetric stretching-vibration and
$d_{\mathrm{OHa}}=({d_{\mathrm{OH}}}_{1}-{d_{\mathrm{OH}}}_{2})/2$ for the
anti-symmetric stretching vibration. Fig. S4 compares the power spectra of
these coordinates, obtained from averaging over all 256 H2O molecules in the
aiMD simulation, to power spectra of the $d_{\mathrm{OH}}$ coordinate,
experimental FTIR data and the IR spectra obtained from the total dipole-
moment trajectory including nuclear and electronic charges from the aiMD
simulation and an approximation using partial charges on the nuclear dynamics
of the aiMD simulation, as introduced in section III. The power spectrum of
the $d_{\mathrm{OH}}$ coordinate is blue-shifted with respect to the power
spectra considering the total dipole moment, indicating a slow down of the
collective nuclear and electronic dynamics with respect to the single-molecule
dynamics, as discussed in section II. This could be similar to the effect of
the slow down observed in the IR spectra including electronic degrees of
freedom compared to the nuclei-only IR spectra, discussed in section III. As
expected for the decomposition, the power spectra of $d_{\mathrm{OHs}}$ and
$d_{\mathrm{OHa}}$ exactly sum up to the $d_{\mathrm{OH}}$ power spectrum,
there exists no remaining cross-correlation spectrum, so they form orthogonal
coordinates.

The effective masses of the coordinates are determined by the equipartition
theorem $m_{\mathrm{OH-
sy.}}=k_{B}T/\langle\dot{d}_{\mathrm{OHs}}^{2}\rangle=$1.901\,04\text{\,}\mathrm{u}$$,
$m_{\mathrm{OH-
as.}}=k_{B}T/\langle\dot{d}_{\mathrm{OHa}}^{2}\rangle=$1.845\,71\text{\,}\mathrm{u}$$
and
$m_{\mathrm{OH}}=k_{B}T/\langle\dot{d}_{\mathrm{OH}}^{2}\rangle=$0.936\,477\text{\,}\mathrm{u}$$.
Analytic estimates are determined by the reduced mass of the $d_{\mathrm{OH}}$
coordinate $m_{\rm OH}=(m_{\rm O}m_{\rm H})/(m_{\rm O}+m_{\rm
H})=$0.947\,617\text{\,}\mathrm{u}$$. It follows for the masses of the
$d_{\mathrm{OHs}}$ and $d_{\mathrm{OHa}}$ coordinates $m_{\mathrm{OH-
as.}}=m_{\mathrm{OH-sy.}}=2m_{\mathrm{OH}}=$1.895\,23\text{\,}\mathrm{U}$$.

\begin{overpic}[width=212.47617pt]{{./f/dft_nvt_256_f_oh_afit}.pdf}
\put(2.0,55.0){\large\bf A}
\end{overpic}\begin{overpic}[width=212.47617pt]{{f/f_oh_1vsA_afit}.pdf}
\put(2.0,56.0){\large\bf B} \end{overpic}

\begin{overpic}[width=212.47617pt]{{oh_kernel/kernel_rs3Exp8Osc_dft_nvt_256_noInt}.png}
\put(2.0,56.0){\large\bf C} \end{overpic}

\begin{overpic}[width=212.47617pt]{{oh_kernel/kernel_ft1wSpec3Exp8OscDelta_dft_nvt_256}.pdf}
\put(2.0,56.0){\large\bf D}
\end{overpic}\begin{overpic}[width=212.47617pt]{{./oh_spectra/spec_ftAFitDelta3Exp8Osc_dft_nvt_256}.pdf}
\put(2.0,55.0){\large\bf E} \end{overpic}

Figure S5: Results for the stretch coordinate $d_{\mathrm{OH}}$ from aiMD
simulations. A, B: Potential $U(d_{\mathrm{OH}})$ for 256 H2O in the liquid
phase (grey solid line) and for a single H2O (green solid line in B), both
compared to the non-harmonic fit according to eq. (LABEL:eq:pot_param) in the
main text (blue broken line) and the harmonic part (grey dotted line). C, D:
Friction as a function of time and frequency (grey lines) compared with the
fit according to eq. (LABEL:eq:kernel_param) in the main text (blue lines).
Real and imaginary parts in (D) are shown as solid and broken lines, the
spectrum on top is the full absorption spectrum from aiMD. The dotted
horizontal line in (D) shows the constant real friction $\gamma_{\rm
OH}=\widetilde{\Gamma}^{\prime}(f_{\rm OH})$ evaluated at the OH stretch
vibrational frequency $f_{\rm OH}=$3400\text{\,}\mathrm{c}\mathrm{m}^{-1}$$.
The grey circle denotes the static friction $\widetilde{\Gamma}^{\prime}(0)$.
E: Power spectrum $\omega\widetilde{\chi}^{\prime\prime}$ (grey solid line)
compared to models of varying complexity: Lorentzian with harmonic potential
and constant friction $\gamma_{\rm OH}$ (grey dotted line), non-harmonic
potential and constant friction $\gamma_{\rm OHs}$ (blue broken line),
harmonic potential and frequency-dependent friction $\widetilde{\Gamma}(f)$
(purple solid line), non-harmonic potential and frequency-dependent friction
$\widetilde{\Gamma}(f)$ (blue solid line).

\begin{overpic}[width=212.47617pt]{{./f/dft_nvt_256_f_oha_afit}.pdf}
\put(2.0,55.0){\large\bf A}
\end{overpic}\begin{overpic}[width=212.47617pt]{{f/f_oha_1vsA_afit}.pdf}
\put(2.0,56.0){\large\bf B} \end{overpic}

\begin{overpic}[width=212.47617pt]{{oha_kernel/kernel_rs2Exp9Osc_dft_nvt_256_noInt}.png}
\put(2.0,56.0){\large\bf C} \end{overpic}

\begin{overpic}[width=212.47617pt]{{oha_kernel/kernel_ft1wSpec2Exp9OscDelta_dft_nvt_256}.pdf}
\put(2.0,56.0){\large\bf D}
\end{overpic}\begin{overpic}[width=212.47617pt]{{./oha_spectra/spec_ftAFitDelta2Exp9Osc_dft_nvt_256}.pdf}
\put(2.0,55.0){\large\bf E} \end{overpic}

Figure S6: Results for the anti-symmetric stretch coordinate
$d_{\mathrm{OHa}}$ from aiMD simulations. A, B: Potential
$U(d_{\mathrm{OHa}})$ for 256 H2O in the liquid phase (grey solid line) and
for a single H2O (green solid line in B), both compared to the non-harmonic
fit according to eq. (LABEL:eq:pot_param) in the main text (blue broken line)
and the harmonic part (grey dotted line). C, D: Friction as a function of time
and frequency (grey lines) compared with the fit according to eq.
(LABEL:eq:kernel_param) in the main text (blue lines). Real and imaginary
parts in (D) are shown as solid and broken lines, the spectrum on top is the
full absorption spectrum from aiMD. The dotted horizontal line in (D) shows
the constant real friction $\gamma_{\rm
OHa}=\widetilde{\Gamma}^{\prime}(f_{\rm OHa})$ evaluated at the OH stretch
vibrational frequency $f_{\rm OHa}=$3440\text{\,}\mathrm{c}\mathrm{m}^{-1}$$.
The grey circle denotes the static friction $\widetilde{\Gamma}^{\prime}(0)$.
E: Power spectrum $\omega\widetilde{\chi}^{\prime\prime}$ (grey solid line)
compared to models of varying complexity: normal mode of single H2O (broken
vertical line), Lorentzian with harmonic potential and constant friction
$\gamma_{\rm OHa}$ (grey dotted line), non-harmonic potential and constant
friction $\gamma_{\rm OHa}$ (blue broken line), harmonic potential and
frequency-dependent friction $\widetilde{\Gamma}(f)$ (purple solid line), non-
harmonic potential and frequency-dependent friction $\widetilde{\Gamma}(f)$
(blue solid line).

The potentials, power spectra and analysis in terms of frequency-dependent
friction of the $d_{\mathrm{OHs}}$ and $\phi_{\mathrm{HOH}}$ coordinates are
shown in the main text. The remaining discussion is shown here for the
$d_{\mathrm{OH}}$ coordinate in fig. S5 and for the $d_{\mathrm{OHa}}$
coordinate in fig. S6. The potential of the $d_{\mathrm{OH}}$ coordinate shows
strong non-harmonic contributions in fig. S5A, seen from the reduced potential
coefficients $\tilde{k}_{3}=k_{3}/k_{B}T(k/k_{B}T)^{-3/2}=-0.128$ and
$\tilde{k}_{4}=k_{4}/k_{B}T(k/k_{B}T)^{-2}=0.00503$. The potential of the
$d_{\mathrm{OHa}}$ coordinate is symmetric by definition and therefore has a
negligible cubic and only a small quartic contribution $\tilde{k}_{4}=-0.0141$
in fig. S6A, when comparing to the potential of the $d_{\mathrm{OHs}}$
coordinate, shown in fig. LABEL:ohs_frictionA in the main text to have a
significant cubic contribution. When compared to the potential of a single
water molecule in figs. S5B and S6B, the $d_{\mathrm{OH}}$ coordinate shows a
shift of the minimum from $97.37\text{\,}\mathrm{p}\mathrm{m}$ to
$99.06\text{\,}\mathrm{p}\mathrm{m}$, i.e. elongation of the bond length due
to hydrogen bonding. As expected the potential of $d_{\mathrm{OHa}}$
coordinate is centered around zero for both systems. Both coordinates show a
significant potential softening in the liquid phase, which follows from
comparison of the dominant harmonic contributions of
$k/(k_{B}T)=$0.139\text{\,}\mathrm{p}\mathrm{m}^{-2}$$ for $d_{\mathrm{OH}}$
and $k/(k_{B}T)=$0.272\text{\,}\mathrm{p}\mathrm{m}^{-2}$$ for
$d_{\mathrm{OHa}}$ with the respective gas-phase values of
$k/(k_{B}T)=$0.183\text{\,}\mathrm{p}\mathrm{m}^{-2}$$ for $d_{\mathrm{OH}}$
and $k/(k_{B}T)=$0.352\text{\,}\mathrm{p}\mathrm{m}^{-2}$$ for
$d_{\mathrm{OHa}}$.

The memory kernels in the time domain, shown in figs. S5C for the
$d_{\mathrm{OH}}$ coordinate and S6C for the $d_{\mathrm{OHa}}$ coordinate,
are discussed further below in fig. S7. The frequency-dependent friction,
shown in figs. LABEL:ohs_frictionD in the main text for the $d_{\mathrm{OHs}}$
coordinate, S5D for the $d_{\mathrm{OH}}$ coordinate and S6D for the
$d_{\mathrm{OHa}}$ coordinate show remarkable differences. The frequency-
dependent friction of the $d_{\mathrm{OHs}}$ coordinate shows a much stronger
contribution in the regime of the HOH bending band at around
$1650\text{\,}\mathrm{c}\mathrm{m}^{-1}$, indicating that it couples more
intensely to the HOH bending mode than the $d_{\mathrm{OHa}}$ mode.

Figs. S5E and S6E show the power spectra of the two coordinates, compared to
the different GLE models. As shown in the main text already for the
$d_{\mathrm{OHs}}$ and $\phi_{\mathrm{HOH}}$ coordinates, the frequency-
dependent-friction model including an non-harmonic potential outperforms the
other models for all coordinates. In case of the $d_{\mathrm{OHa}}$
coordinate, shown in fig. S6E, the frequency-dependent-friction model with the
harmonic potential also performs very well, due to the negligible non-harmonic
contribution to the potential in fig. S6A.

\begin{overpic}[width=212.47617pt]{{./kernels_dft_nvt_256_compare_noInt}.pdf}
\end{overpic} Figure S7: Comparison of time-dependent friction kernels
extracted from the aiMD simulations for the $d_{\mathrm{OH}}$ (grey solid
line), $d_{\mathrm{OHs}}$ (blue solid line) and $d_{\mathrm{OHa}}$ (purple
solid line) coordinates.

The time-dependent friction kernels of the $d_{\mathrm{OH}}$,
$d_{\mathrm{OHs}}$ and $d_{\mathrm{OHa}}$ coordinates are compared in fig. S7,
scaled by the mass of the respective coordinate. The time-dependent friction
kernel of the $d_{\mathrm{OHa}}$ coordinate differs significantly from the one
of the $d_{\mathrm{OHs}}$ coordinate. However, all friction kernels show
various oscillating decay time scales between
$10\text{\,}\mathrm{f}\mathrm{s}$ and $5\text{\,}\mathrm{p}\mathrm{s}$, which
is discussed in detail the main text for the example of $d_{\mathrm{OHs}}$.

##  XII Constant-friction line shape

The dependence of the Lorentzian line shape function eq. (S51) on the friction
coefficient $\gamma$ for the $d_{\mathrm{OHs}}$ coordinate is illustrated in
fig. S8A and for the $\phi_{\mathrm{HOH}}$ coordinate in fig. S8B. A variation
of $\gamma$ does not shift the peak position but significantly changes the
width. This follows analytically from the maximum of eq. (S51)

$\displaystyle\frac{\partial}{\partial\omega}\left[\omega\widetilde{\chi}^{\prime\prime}(\omega)\right]$
$\displaystyle=\frac{L^{2}}{k_{B}T}\left[-\frac{8\tau\omega\left(\omega^{4}\tau_{\omega}^{4}-1\right)}{\left(\omega^{2}\left(\tau^{2}-8\tau_{\omega}^{2}\right)+4\omega^{4}\tau_{\omega}^{4}+4\right){}^{2}}\right]$
(S93)
$\displaystyle\overset{!}{=}0\quad\text{for}\quad\omega=\tau_{\omega}^{-1},$
(S94)

which does not depend on $\gamma$. When considering non-harmonic effects in
the potential, simulated power spectra still show no significant shifts for
variation of the friction coefficient, as shown in figs. S8C and D.

\begin{overpic}[width=212.47617pt]{{ohs_spectra/spec_DeltaFitMaxVar_dft_nvt_256}.pdf}
\put(2.0,55.0){\large\bf A}
\end{overpic}\begin{overpic}[width=212.47617pt]{{hoh_spectra/spec_DeltaFitMaxVar_dft_nvt_256}.pdf}
\put(2.0,55.0){\large\bf B}
\end{overpic}\begin{overpic}[width=212.47617pt]{{ohs_spectra/spec_DeltaAFitMaxVar_dft_nvt_256}.pdf}
\put(2.0,55.0){\large\bf C}
\end{overpic}\begin{overpic}[width=212.47617pt]{{hoh_spectra/spec_DeltaAFitMaxVar_dft_nvt_256}.pdf}
\put(2.0,55.0){\large\bf D} \end{overpic}

Figure S8: Power spectra, obtained from the aiMD simulations (grey solid
lines) and compared to spectra from a constant friction model with the
harmonic fit to the potential (A, B) and considering non-harmonic effects (C,
D). The spectra are shown for variations of the friction coefficient
$\gamma_{\rm OHS}=\Gamma^{\prime}(f_{\rm OHS})$ and $\gamma_{\rm
HOH}=\Gamma^{\prime}(f_{\rm OHS})$, where $f_{\rm
OHS}=$3390\text{\,}\mathrm{c}\mathrm{m}^{-1}$$ and $f_{\rm
HOH}=$1650\text{\,}\mathrm{c}\mathrm{m}^{-1}$$ are at the maxima of the
respective power spectra from the simulations. A, C: $d_{\mathrm{OHs}}$
coordinate and B, D: $\phi_{\mathrm{HOH}}$ coordinate.

##  XIII Simulations of single H2O in the NVE ensemble

For the aiMD simulations of single H2O molecules, representing the gas phase,
47 initial configurations were sampled from a
$25\text{\,}\mathrm{p}\mathrm{s}$ NVT simulation and subsequently simulated in
the NVE ensemble. The NVT simulation was temperature-controlled using an
individual thermostat with a time constant of
$10\text{\,}\mathrm{f}\mathrm{s}$ for each atom. The NVE simulations were each
run for $10\text{\,}\mathrm{p}\mathrm{s}$ with a time step of
$0.25\text{\,}\mathrm{f}\mathrm{s}$. The distributions of their initial
configurations are shown in fig. S9 to sample well the equilibrium
distributions from the NVT trajectory.

\begin{overpic}[width=212.47617pt]{{f/1dist_ohs_nve}.png}
\put(2.0,55.0){\large\bf A} \end{overpic}

\begin{overpic}[width=212.47617pt]{{f/1dist_hoh_nve}.png}
\put(2.0,55.0){\large\bf B} \end{overpic}

\begin{overpic}[width=212.47617pt]{{f/1dist_oha_nve}.png}
\put(2.0,55.0){\large\bf C} \end{overpic}

\begin{overpic}[width=212.47617pt]{{f/1dist_oh_nve}.png}
\put(2.0,55.0){\large\bf D} \end{overpic}

Figure S9: Comparison of the distributions for the different coordinates
sampled from the aiMD simulations of a single H2O molecule under NVT
conditions (blue solid lines) and 47 samples taken from the NVT trajectory as
initial conditions for equivalent simulations under NVE conditions (orange
dots). A: $d_{\mathrm{OHs}}$ coordinate, B: $\phi_{\mathrm{HOH}}$ coordinate,
C: $d_{\mathrm{OHa}}$ coordinate and D: $d_{\mathrm{OH}}$ coordinate.

##  XIV Wiener-Khintchine theorem

The correlation function $C_{xy}(t)$ of two stochastic processes $x(t)$ and
$y(t)$ limited to the interval $[0,L_{t}]$ is efficiently computed from the
Fourier-transformed expressions $\tilde{x}(\omega)$ and $\tilde{y}(\omega)$
according to

$\displaystyle
C_{xy}(t)=\frac{1}{2\pi(L_{t}-t)}\int_{-\infty}^{\infty}d\omega\ e^{-i\omega
t}\tilde{x}(\omega)\tilde{y}^{*}(\omega),$ (S95)

where the asterisk denotes the conjugate form. This is known as the Wiener-
Khintchine theorem [7]. Both sides of eq. (S95) are Fourier-transformed to
give

$\displaystyle\int_{-\infty}^{\infty}dt\ e^{i\omega t}\ 2\pi
L_{t}\left(1-\frac{t}{L_{t}}\right)C_{xy}(t)$
$\displaystyle=\tilde{x}(\omega)\tilde{y}^{*}(\omega),$ (S96)

which in the limit of large $L_{t}$ reduces to

$\displaystyle\widetilde{C}_{xy}(\omega)$
$\displaystyle=L_{t}^{-1}\tilde{x}(\omega)\tilde{y}^{*}(\omega).$ (S97)

Eq. (S95) can be derived starting off with the definition of the correlation
function

$\displaystyle C_{xy}(t)=\frac{1}{L_{t}-t}\int_{0}^{L_{t}-t}dt^{\prime}\
x(t^{\prime}+t)y(t^{\prime}),$ (S98)

and making use of the convolution theorem

$\displaystyle C_{xy}(t)$
$\displaystyle=\frac{1}{4\pi^{2}(L_{t}-t)}\int_{0}^{L_{t}-t}dt^{\prime}$
$\displaystyle\ \ \int_{-\infty}^{\infty}d\omega\
e^{-i\omega(t+t^{\prime})}\tilde{x}(\omega)\int_{-\infty}^{\infty}d\omega^{\prime}\
e^{-i\omega^{\prime}t^{\prime}}\tilde{y}(\omega^{\prime})$
$\displaystyle=\frac{1}{4\pi^{2}(L_{t}-t)}\int_{-\infty}^{\infty}d\omega\
e^{-i\omega t}\tilde{x}(\omega)\int_{-\infty}^{\infty}d\omega^{\prime}\
\tilde{y}(\omega^{\prime})$ $\displaystyle\ \
\int_{0}^{L_{t}-t}dt^{\prime}e^{-it^{\prime}(\omega+\omega^{\prime})}$
$\displaystyle=\frac{1}{4\pi^{2}(L_{t}-t)}\int_{-\infty}^{\infty}d\omega\
e^{-i\omega t}\tilde{x}(\omega)\int_{-\infty}^{\infty}d\omega^{\prime}\
\tilde{y}(\omega^{\prime})$ $\displaystyle\ \
2\pi\delta(\omega+\omega^{\prime})$
$\displaystyle=\frac{1}{2\pi(L_{t}-t)}\int_{-\infty}^{\infty}d\omega\
e^{-i\omega t}\tilde{x}(\omega)\tilde{y}(-\omega),$ (S99)

noting that $\tilde{y}(-\omega)=\tilde{y}^{*}(\omega)$ for a real function
$y(t)$ in order to obtain eq. (S95).

## References

  * Kohler et al. [1972] Kohler, F.; Findenegg, G.; Fischer, J.; Posch, H. _The Liquid State_ ; Verlag Chemie Weinheim, 1972.
  * Carlson et al. [2020] Carlson, S.; Brünig, F. N.; Loche, P.; Bonthuis, D. J.; Netz, R. R. Exploring the Absorption Spectrum of Simulated Water from MHz to Infrared. _J. Phys. Chem. A_ 2020, _124_ , 5599–5605.
  * Daldrop et al. [2018] Daldrop, J. O.; Kappler, J.; Brünig, F. N.; Netz, R. R. Butane dihedral angle dynamics in water is dominated by internal friction. _Proc. Natl. Acad. Sci._ 2018, _115_ , 5169–5174.
  * Harp and Berne [1970] Harp, G.; Berne, B. Time-Correlation Functions, Memory Functions, and Molecular Dynamics. _Phys. Rev. A_ 1970, _2_ , 975–996.
  * Schrader [1995] Schrader, B., Ed. _Infrared and Raman Spectroscopy: Methods and Applications_ ; Wiley-VCH, 1995.
  * Bertie and Lan [1996] Bertie, J. E.; Lan, Z. Infrared Intensities of Liquids XX: The Intensity of the OH Stretching Band of Liquid Water Revisited, and the Best Current Values of the Optical Constants of H2O(l) at 25∘C between 15,000 and 1 cm-1. _Appl. Spectrosc._ 1996, _50_ , 1047–1057.
  * Wiener [1930] Wiener, N. Generalized Harmonic Analysis. _Acta Math._ 1930, _55_ , 117–258.

  *[IR]: infrared
  *[aiMD]: ab initio Molecular Dynamics
  *[GLE]: generalized Langevin equation
  *[FTIR]: Fourier-transform infrared spectroscopy