Modeling intermolecular and intramolecular modes of liquid water using multiple heat baths: Machine learning approach

Seiji Ueno [email protected] [ Yoshitaka Tanimura [email protected]

Abstract

The vibrational motion of molecules in dissipative environments, such as solvation and protein molecules, is composed of contributions from both intermolecular and intramolecular modes. The existence of these collective modes introduces difficulty into quantum simulations of chemical and biological processes. In order to describe the complex molecular motion of the environment in a simple manner, we introduce a system-bath model in which the intramolecular modes with anharmonic mode-mode couplings are described by a system Hamiltonian, while the other degrees of freedom, arising from the environmental molecules, are described by heat bath. Employing a machine-learning based approach, we determine not only the system parameters of the intramolecular modes but also the spectral distribution of the system-bath coupling to describe the intermolecular modes, using the atomic trajectories obtained from molecular dynamics (MD) simulations. The capabilities of the present approach are demonstrated for liquid water using MD trajectories calculated from the SPC/E model and the polarizable water model for intramolecular and intermolecular vibrational spectroscopies (POLI2VS) by determining the system parameters describing the symmetric-stretch, asymmetric-stretch and bend modes with intramolecular interactions and the bath spectral distribution functions for each intramolecular mode representing the interaction with the intra-molecular modes. From these results, we were able to elucidate the energy relaxation pathway between the intramolecular modes and the intermolecular modes in a non-intuitive manner

keywords:

American Chemical Society, LaTeX

HPCS]HPC Systems Inc., Japan \altaffiliationDepartment of Chemistry, Kyoto University, Kyoto, Japan \altaffiliationDepartment of Chemistry, Kyoto University, Kyoto, Japan \abbreviations

1 1. INTRODUCTION

Elucidating the effects of molecular environments on chemical and biological processes in both classical and quantum regimes has been an important problem in chemical physics over the last four decades. ^{1, 2, 3, 4, 5} Theories of quantum open systems have been used to construct models of practical interest, ^{6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 25, 26, 27, 28, 29, 19, 20, 21, 22, 23, 24} in particular to account for line shapes in NMR^{25, 26} and linear and nonlinear laser spectra.^{27, 28, 29} The key feature of the environment is that it gives rise to irreversible dynamics through which the system evolves toward the thermal equilibrium state at finite temperature.⁵ The environmental effects that arise from solid state materials, solvation, and protein molecules are generally described by an interaction between a primary molecular system coupled to a harmonic oscillator bath (HOB). Two commonly used models of this kind are the Caldeira-Leggett (Brownian) model², Brownian Oscillator model (BOM) and the spin-Boson model.² The HOB model, whose distribution takes a Gaussian form, exhibits wide applicability, despite its simplicity. This is because the influence of the environment can in many cases be approximated by a Gaussian process, due to the cumulative effect of the large number of weak environmental interactions. In such a situation, the ordinary central limit theorem is applicable, and hence the Gaussian distribution function is appropriate.⁵ The distinctive features of the HOB model are determined by the spectral distribution function (SDF) of the coupling strength, $J(\omega)=\sum c_{j}^{2}\delta(\omega-\omega_{j})$ , where $c_{j}$ is the coupling strength between the system and the $j$ th bath oscillator, with frequency $\omega_{j}$ . By properly choosing the form of the SDF, the properties of the bath can be adjusted to represent a variety of environments consisting of, for example, solvates and protein molecules.

Several analytical and numerical approaches have been developed for such models to treat the quantum and classical motion of molecules in condensed phases. These include the generating functional approach^{30, 31, 32}, the Redfield equation approach,^{33, 34, 35} the quasi-adiabatic propagator path integral (QUAPI) approach, ^{36, 37} the multi-configurational time-dependent Hartree (MCTDH) approach, ^{38, 39, 40, 41} and the reduced hierarchy equations of motion (HEOM) approach.^{42, 43, 44, 45, 46, 47, 48, 49, 50, 51} These approaches have been applied to problems of practical interest, in particular to the investigation of chemical reaction processes, ^{6, 7, 8, 9} non-adiabatic transitions,^{10, 11, 12, 13} quantum device systems,^{14, 15} ratchet rectification,^{16, 17, 18} exciton transfer,^{19, 20, 21, 22, 23, 24} and the analysis of the linear and nonlinear laser spectra.^{27, 28, 29}

Because the key feature of environments is determined by the system-bath coupling and the SDF, a methodology to determine these is significant. Typically, the SDF is estimated from linear and nonlinear infrared and Raman spectra, both experimentally^{52, 53, 54, 55, 56, 57} and numerically.^{61, 62, 63, 58, 59, 60, 21, 22, 23} The SDFs of entire molecular systems have also been evaluated directly from molecular dynamics (MD) trajectories on the basis of normal modes^{64, 65, 66, 68, 67, 69, 70, 71, 72, 73, 74, 75}, the velocity-velocity autocorrelation function,^{65, 69} and Fourier transformations.⁷⁶ While the method based on optical spectra does not allow estimation of optically inactive system modes, the approaches based on MD trajectories have the capability to describe all intermolecular and intramolecular vibrational motions. However, in such trajectory-based methods, it is not easy to separate the system and the heat-bath, in particular in the case of the intermolecular modes. This is because these methods are based on naive statistics of the complex and irregular motion of molecular trajectories, and thus, with these methods, it is not easy to distinguish the bath modes from the background noise. Moreover, because the SDF generally depends on the system parameters, for example depending on the intermolecular mode-mode coupling strengths, with the existing methods, it is not easy to obtain the system parameters and the bath SDF simultaneously.

In this study, we employ a data mining approach, which allows us to obtain variable information from massive data in physics and chemistry^{77, 78, 79} on the basis of the statistical analysis.^{80, 81} We then develop an efficient algorithm to optimize both the system parameters and the SDF through analysis of MD trajectories. This algorithm allows us to obtain not only the SDF of each intramolecular mode, but also the variance of the bath-correlation function that describes the vibrational dephasing of molecular modes. We demonstrate this approach for liquid water using MD trajectories obtained from SPC/E model^{82, 83} and POLI2VS model.^{58, 59}

This paper is organized as follows. In Sec. 2, we introduce a multi-mode system-bath model and describe the data-mining approach that we use to determine the system parameters, the system-bath interaction, and the SDFs. In Sec. 3, we present model potentials for intramolecular modes of liquid water in the framework of BOM. In Sec. 4, we present the details of the calculation. In Sec.5, we present the results for liquid water obtained through analysis of the MD trajectories obtained using various force fields. Section 6 is devoted to concluding remarks.

2 2. THEORY

2.1 A model with multiple heat baths

In order to describe the vibrational modes of molecular liquids, we consider a model that consists of primary intramolecular modes coupled to intermolecular environmental modes, which are regarded as bath systems. These bath systems are represented by ensembles of harmonic oscillators. ²⁷ The model is constructed by extending a Brownian (or Caldeira-Leggett) Hamiltonian to include a nonlinear system-bath interaction,^{84, 85, 86, 87, 88, 89, 90, 91, 92} which causes the frequency and amplitude of the intramolecular modes to vary in time or to be inhomogeneously distributed. We can describe both situations within a unified framework by adjusting the nonlinear system-bath coupling strength.^{5, 29} The total Hamiltonian is expressed as

\displaystyle H=\sum_{s}\left(H^{(s)}_{\mathrm{S}}+H^{(s)}_{\mathrm{B}}+H^{(s)}_{\mathrm{I}}\right)+\sum_{s\neq s^{\prime}}U_{s,s^{\prime}}(q_{s},q_{s^{\prime}}),

(1)

where

\displaystyle H^{(s)}_{\mathrm{S}}=\frac{p_{s}^{2}}{2m_{s}}+U_{s}(q_{s})

(2)

is the Hamiltonian for the $s$ th mode, with mass $m_{s}$ , coordinate ${q_{s}}$ , momentum ${p_{s}}$ , and potential $U_{s}(q_{s})$ . The interaction between the modes $s$ and $s^{\prime}$ is given by

U_{s,s^{\prime}}(q_{s},q_{s^{\prime}})=g_{11}q_{s}q_{s^{\prime}}+\frac{1}{2}\left(g_{21}q_{s}^{2}q_{s^{\prime}}+g_{12}q_{s}q^{2}_{s^{\prime}}\right).

(3)

The forms of the potentials and interactions should be chosen to accurately describe the system dynamics in a simple manner. The bath Hamiltonian for the $s$ th mode is expressed as

\displaystyle H^{(s)}_{\mathrm{B}}=\sum_{j_{s}}\left(\frac{p_{j_{s}}^{2}}{2m_{j_{s}}}+\frac{m_{j_{s}}\omega_{j_{s}}^{2}x_{j_{s}}^{2}}{2}\right)+\sum_{j_{s}}\left(\frac{\alpha_{j_{s}}^{2}V_{s}^{2}({q_{s}})}{2m_{j_{s}}\omega_{j_{s}}^{2}}\right),

(4)

where the momentum, coordinate, mass, and frequency of the $j_{s}$ th bath oscillator are given by $p_{j_{s}}$ , $x_{j_{s}}$ , $m_{j_{s}}$ and $\omega_{{j_{s}}}$ , respectively. The last term in the above is a counter term, which maintains the translational symmetry of the system in the case ${U}_{s}(q_{s})={U}_{s,s^{\prime}}(q_{s},q_{s^{\prime}})=0$ .¹ The system-bath interaction is given by

H^{(s)}_{\mathrm{I}}=-V_{s}(q_{s})\sum_{j_{s}}\alpha_{j_{s}}x_{j_{s}},

(5)

which consists of linear-linear (LL) and square-linear (SL) interactions,

V_{s}(q_{s})\equiv V^{(s)}_{\mathrm{LL}}q_{s}+\frac{1}{2}V^{(s)}_{\mathrm{SL}}q_{s}^{2},

(6)

with coupling strengths $V^{(s)}_{\mathrm{LL}}$ , $V^{(s)}_{\mathrm{SL}}$ , and $\alpha_{j_{s}}$ .⁵ As shown in Refs. [ 5, 29], while the LL interaction contributes mainly to energy relaxation, the SL system-bath interaction leads to vibrational dephasing in the slow modulation case, due to the frequency fluctuation of the system vibrations.^{88, 89, 90, 91, 93, 92} Then, combining Eq.(4) and Eq.(5), we obtain

H^{(s)}_{\mathrm{B}}+H^{(s)}_{\mathrm{I}}=\sum_{j_{s}}\left(\frac{p_{j_{s}}^{2}}{2m_{j_{s}}}+\frac{m_{j_{s}}\omega_{j_{s}}^{2}\tilde{x}_{j_{s}}^{2}}{2}\right),

(7)

where $\tilde{x}_{j_{s}}=x_{j_{s}}-({\alpha_{j_{s}}V_{s}({q_{s}})})/({2m_{j_{s}}\omega^{2}_{j_{s}}})$ is the re-oriented bath coordinate. This model has been used to derive predictions for the study of single-mode systems employing 2D Raman^{88, 89, 90, 91}, 2D THz-Raman,⁶² and 2D IR signals,^{93, 92} and for the study of two-mode systems employing 2D IR^{94, 95} and 2D THz-IR signals.⁶³ Here, we apply this model to describe the intermolecular modes of molecular liquids. We assume that the influences of the fluctuation and dissipation on individual modes are all independent and that the correlations of the fluctuations among different modes can be ignored.^{94, 95} Then, the SDF is defined as

J_{s}(\omega)\equiv\hbar\sum_{j_{s}}\frac{\alpha_{j_{s}}^{2}}{2m_{j_{s}}\omega_{j_{s}}}\delta(\omega-\omega_{j_{s}}),

(8)

This definition characterizes the nature of the bath.

2.2 Optimizing the likelihood probability distribution

Let us onsider the trajectory of the $k$ th molecule in phase space, expressed as $\left(\mathbf{q}_{k}(t),\mathbf{p}_{k}(t)\right)$ , describing the intramolecular motion of interest obtained from the MD simulation. We compute the set of trajectories $\left(\mathbf{q}_{k}(t_{i}+i\Delta t),\mathbf{p}_{k}(t_{i}+i\Delta t)\right)$ with time step $\Delta t$ for all integer $i$ satisfying $0\leq i\leq N-1$ , where $N$ is the total number of time steps. Using a data mining approach, we attempt to reproduce the MD trajectories for the intramolecular modes by adjusting the system parameters in Eq.(1) and the SDF for each system mode. The trajectory of the $j_{s}$ th bath oscillator for the $s$ th system mode in Eqs. (4) and (5) is assumed to take the form

\tilde{x}_{j_{s}}(t)=A_{j_{s}}\sin(\omega_{j_{s}}t+\phi_{j_{s}}),

(9)

where $A_{j_{s}}$ and $\phi_{j_{s}}$ are the amplitude and phase of the ${j_{s}}$ th bath oscillator. The phase $\phi_{j_{s}}$ is chosen randomly to avoid recursive motion of the oscillator. This implies that the bath oscillators described by Eq.(7) are harmonic for any form of the system-bath coupling.

In the LL coupling case, the bath parameters and the system-bath interactions are expressed as a set of latent variables in the machine learning context, defined as

\mathbf{z}_{k}=\left(\left\{c^{k}_{j_{1}}\right\},\left\{c^{k}_{j_{2}}\right\},\cdots,\left\{c^{k}_{j_{N}}\right\}\right)

(10)

where $\left\{c^{k}_{j_{s}}\right\}$ is the set of bath coupling parameters with the element,

c^{k}_{j_{s}}=V_{\mathrm{LL}}^{(s)}\alpha^{k}_{j_{s}}A^{k}_{j_{s}},

(11)

for the mode $s$ .

A set consisting of the system potential parameters and the SDFs is denoted by $\mathbf{\Sigma}$ . In the present approach, these are optimized in the same manner. The trajectory at time $t_{i}+i\Delta t$ for integer $i$ ( $1\leq i\leq N$ ) that is calculated using the HOB model on the basis of the MD trajectories at $t_{i}+(i-1)\Delta t$ is then expressed as

(\tilde{\mathbf{q}}_{k}(t_{i}+i\Delta t),\tilde{\mathbf{p}}_{k}(t_{i}+i\Delta t))=\hat{L}(\Delta t;\mathbf{z}_{k},\mathbf{\Sigma})(\mathbf{q}_{k}(t_{i}+(i-1)\Delta t),\mathbf{p}_{k}(t_{i}+(i-1)\Delta t)),

(12)

where $\hat{L}(\Delta t;\mathbf{z}_{k},\mathbf{\Sigma})$ is the Liouville operator of $\mathbf{q}_{k}$ and $\mathbf{p}_{k}$ for the model Hamiltonian (1) with the bath functions given in Eq. (9). The differences between the coordinates and momenta for a given MD trajectory and the corresponding HOB trajectory are defined as $\Delta\mathbf{q}_{k}(t_{i}+i\Delta t)\equiv\left(\tilde{\mathbf{q}}_{k}(t_{i}+i\Delta t)-\mathbf{q}_{k}(t_{i}+i\Delta t)\right)$ and $\Delta\mathbf{p}_{k}(t_{i}+i\Delta t)\equiv\left(\tilde{\mathbf{p}}_{k}(t_{i}+i\Delta t)-\mathbf{p}_{k}(t_{i}+i\Delta t)\right)$ .

It should be noted that the influence of the environmental molecules can in many cases be approximated by a Gaussian process, because it can be treated as the cumulative effect of a large number of weak environmental interactions, in which case the ordinary central limit theorem is applicable. In this case, the distribution function of the HOB model, which is Gaussian, is appropriate. Thus, the difference between the actual influence of the environmental molecules and that described by the HOB model is due to the non-Gaussian aspect of the molecular motion.

Moreover, if we choose the time step, $\Delta t$ , to be sufficiently large in comparison with the time scale of the frequency fluctuations of the intramolecular modes, we can regard the dynamics to be Markovian, in which case $\Delta\mathbf{q}_{k}(t_{i}+i\Delta t)$ and $\Delta\mathbf{p}_{k}(t_{i}+i\Delta t)$ for different $i$ are not correlated. Here, we indeed consider this situation. Thus, we assume that the time evolution of these variables obeys a Gaussian-Markovian process, while the intramolecular modes themselves are described by the system Hamiltonian.

The mean squares of the differences given above are is expressed as

\sigma_{qq}\delta\mathbf{q}_{k}^{2}(\Delta t)=\frac{1}{N}\sum_{i=1}^{N}\Delta\mathbf{q}_{k}^{2}(t_{i}+i\Delta t)

(13)

and

\sigma_{pp}\delta\mathbf{p}_{k}^{2}(\Delta t)=\frac{1}{N}\sum_{i=1}^{N}\Delta\mathbf{p}_{k}^{2}(t_{i}+i\Delta t),

(14)

where $\sigma_{qq}$ and $\sigma_{pp}$ represent the deviations between the BOM and MD trajectories. We next introduce a joint probability distribution function $\delta P\left(\mathbf{q}_{k}(\Delta t),\mathbf{p}_{k}(\Delta t)\mid\mathbf{q}_{k},\mathbf{p}_{k};\mathbf{z}_{k},\mathbf{\Sigma}\right)$ , which describes the time evolution of the probability distribution from the state at time $t$ to that at time $t+\Delta t$ for a given set of system parameters, $\mathbf{\Sigma}$ , and bath parameters, $\mathbf{z}_{k}$ :

\delta P\left(\mathbf{q}_{k}(\Delta t),\mathbf{p}_{k}(\Delta t)\mid\mathbf{q}_{k},\mathbf{p}_{k};\mathbf{z}_{k},\mathbf{\Sigma}\right)=\exp\left[-\sigma_{qq}\delta\mathbf{q}_{k}^{2}(\Delta t)-\sigma_{pp}\delta\mathbf{p}_{k}^{2}(\Delta t)\right].

(15)

The parameters $\sigma_{qq}$ and $\sigma_{pp}$ indicate the accuracy of the BOM. In ordinary situations, these take finite values, because the description of molecular trajectories based on the simple BOM is limited, and because the choice of the parameters $\mathbf{z}_{k}$ may not yield the best fit. We determine these parameters using the maximum likelihood series estimating method (MLSEM) and the gradient estimating method. To use MLSEM, we maximize the logarithm likelihood ratio, defined by

G(\mathbf{z}_{k})=\log\left[\delta P\left(\mathbf{q}_{k}(\Delta t),\mathbf{p}_{k}(\Delta t)\mid\mathbf{q}_{k},\mathbf{p}_{k};\mathbf{z}_{k},\mathbf{\Sigma}\right)\right].

(16)

Then, using Eqs. (13)-(16), we obtain

G(\mathbf{z}_{k})=\frac{1}{N}\sum_{i=1}^{N}\left[\Delta\mathbf{q}^{2}_{k}(t_{i}+i\Delta t)+\Delta\mathbf{p}^{2}_{k}(t_{i}+i\Delta t)\right],

(17)

which can be interpreted as the mean squared error of the BOM trajectories from the MD trajectories.

We then estimate $\mathbf{z}_{k}$ using the gradient estimating method for a given trajectory for the $k$ th molecule,

\mathbf{z}_{k}=\mathop{\rm arg~{}min}\limits_{\mathbf{z}_{k}}G(\mathbf{z}_{k}).

(18)

We repeat this procedure for all system molecules, and compute the average over these molecules in order to obtain the optimized set of $\mathbf{z}$ .

2.3 The bath spectral distribution functions

The total energy of the $j_{s}$ th oscillator is given by

E_{j_{s}}=m_{j_{s}}\omega_{j_{s}}^{2}(A^{k}_{j_{s}})^{2}.

(19)

We assume that the probability distribution for each harmonic oscillator is described by the canonical ensemble with energy $E_{j_{s}}$ . The expectation values of the amplitude $A^{k}_{j_{s}}$ are then evaluated as

\left<A^{k}_{j_{s}}\right>=\frac{1}{\sqrt{\pi\beta m_{j_{s}}\omega_{j_{s}}^{2}}},

(20)

where $\beta=1/k_{\mathrm{B}}T$ is the inverse temperature of the system. Substituting $\alpha^{k}_{j_{s}}$ into Eq. (8), and using Eqs. (11), and using (20), we obtain

J^{k}(\omega_{j_{s}})\sim\frac{\pi^{2}}{2V_{\mathrm{LL}}^{(s)}}\beta\omega_{j_{s}}(c^{k}_{j_{s}})^{2}.

(21)

The SDF for the $s$ th mode is evaluated as the average of $J^{k}(\omega_{j_{s}})$ over the samples $k$ .

3 3. INTRAMOLECULAR SYSTEM MODEL FOR LIQUID WATER

Here, we demonstrate the data-mining approach for liquid water to construct the system-bath Hamiltonian. The system Hamiltonian for the intramolecular vibrational modes of a water molecule consists of the symmetric-stretch, asymmetric-stretch and bending modes, which we assume to take the form

\sum_{s}H^{(s)}_{\mathrm{S}}=\sum_{\alpha\in\{\mathrm{O,H_{1},H_{2}}\}}\left(\frac{1}{2}m_{\alpha}\vec{v}_{\alpha}^{2}\right)+U(r_{1},r_{2},\theta),

(22)

where $m_{\alpha}$ and $\vec{v}_{\alpha}$ represent the mass and velocity of the O, H1 and H2 atoms, and $U(r_{1},r_{2},\theta)$ is the water model potential.

We employ two models, a harmonic oscillator (HO) based model, defined as

U_{\mathrm{HO}}(r_{1},r_{2},\theta)=k_{r}\left(r_{1}-r_{0}\right)^{2}+k_{r}\left(r_{2}-r_{0}\right)^{2}+k_{t}\left(\theta_{1}-\theta_{0}\right)^{2},

(23)

where $k_{r}$ and $k_{t}$ are the force constants for the OH stretching and HOH bending motion, and a Morse oscillator (MO) based model, defined as

	$\displaystyle U_{\mathrm{MO}}(r_{1},r_{2},\theta)$		$\displaystyle=D\left(e^{-2a(r_{1}-r_{0})}-2e^{-a(r_{1}-r_{0})}\right)$		(24)
			$\displaystyle+D\left(e^{-2a(r_{2}-r_{0})}-2e^{-a(r_{2}-r_{0})}\right)+k_{t}\left(\theta_{1}-\theta_{0}\right)^{2},$		(24)

where $D$ is the dissociation energy, and $a$ is the width parameter. In the present study, the other intermolecular vibrational modes are treated as the bath oscillators, described by $H_{\mathrm{I}}^{(s)}$ and $H_{\mathrm{B}}^{(s)}$ .^{88, 89, 90, 91, 62, 93, 92, 94, 95, 63} When the intramolecular and intermolecular modes are not well separated, the profiles of SDFs become very sensitive for a choice of system potential. We found that the optimized SDFs for the HO and MO models differ significantly in the POLI2VS case, while they are similar in the SPC/E flexible case. For this reason, here we employ the HO model in the SPC/E flexible case and the MO model in the POLI2VS case, setting $D=444.188$ kJ/mol. The equilibrium values of $r_{0}$ and $\theta_{0}$ for these models are evaluated through use of $\left<r\right>$ and $\left<\theta\right>$ . The trajectories obtained with this procedure are described in the $xyz$ coordinates for each atom. The loss function given in Eq. (16) is evaluated by comparing the predicted $xyz$ trajectories with the MD $xyz$ trajectories.

We describe the intra-vibrational modes in terms of the two HO bond lengths and the HOH bond angle of a water molecule as

	$\displaystyle r_{1}=\left\|\vec{x}_{\mathrm{O}}-\vec{x}_{\mathrm{H1}}\right\|,$		(25)
	$\displaystyle r_{2}=\left\|\vec{x}_{\mathrm{O}}-\vec{x}_{\mathrm{H2}}\right\|,$		(26)
	$\displaystyle\theta=\arccos\left(\frac{\left(\vec{x}_{\mathrm{O}}-\vec{x}_{\mathrm{H1}}\right)\cdot\left(\vec{x}_{\mathrm{O}}-\vec{x}_{\mathrm{H2}}\right)}{r_{1}r_{2}}\right),$		(27)

where $\vec{x}_{\mathrm{O}}$ and $\vec{x}_{\mathrm{H}_{k}}$ are the positions of the oxygen atom and the $k$ th hydrogen atom.

The system coordinates for the symmetric-stretch, asymmetric-stretch, and bend modes are expressed as

	$\displaystyle q_{\mathrm{sym}}=\frac{1}{2}\left(r_{1}+r_{2}-2r_{0}\right),$		(28)
	$\displaystyle q_{\mathrm{asym}}=\frac{1}{2}\left(r_{1}-r_{2}\right),$		(29)
	$\displaystyle q_{\mathrm{bend}}=\theta-\theta_{0},$		(30)

where $r_{0}$ is the equilibrium length of the OH bond. The system and system-bath interaction, $H^{(s)}_{\mathrm{B}}$ and $H^{(s)}_{\mathrm{I}}$ , are described with these system coordinates. The intramolecular mode-mode interaction, $U_{s,s^{\prime}}(q_{s},q_{s^{\prime}})$ , in the system Hamiltonian can be obtained by rewriting the intramolecular potentials appearing in Eqs.(23) and (24) in terms of Eqs.(28)-(30), respectively. In the present approach, the system parameters, such as $k_{r}$ and $k_{t}$ , may vary during optimization. For this reason, we also optimize the mode-mode coupling constants $g_{ij}$ given in Eq. (3).

4 4. CALCULATION DETAILS

The MD trajectories were obtained by carrying out the simulations for liquid water with 216 molecules in a cubic box with periodic boundary conditions. We employed GROMACS 2018 for the flexible SPC/E model^{82, 83} and the code developed by Hasegawa for the POLI2VS model.⁵⁸ The equations of motion were integrated using the velocity-Verlet algorithm with ${\Delta}t=0.1$ fs. The MD trajectories were constructed for each 0.2 fs up to 300 fs at the temperature about $300$ K.

Because we model the fast intermolecular modes, which arise from short-range intermolecular interactions, it is not necessary to carry out large-scale simulations with many molecules.⁹⁶ The external motion of molecules, i.e., the translational and rotational motion, was then eliminated from the obtained trajectories. The time step $\Delta t$ for the BOM Hamiltonian, given in Eq.(1), was set to 0.1 fs. The initial phase, $\phi_{j_{s}}^{k}$ , was generated from a random number uniformly distributed over $[0,\pi)$ at every iteration step. We chose the number of bath oscillator for the $s$ th mode as $N_{s}=1500$ . The frequency of oscillators $\omega_{j_{s}}$ was set to $j_{s}\Delta\omega$ with $\Delta\omega=1\ \mathrm{ps}^{-1}\approx 5.31\ \mathrm{cm}^{-1}$ ; while almost all features of spectroscopic signals of water can be described within this resolution, the highest frequency of bath modes is $8000$ cm^-1 that covers the overtone of OH stretching mode.

We developed the Python code using the TensorFlow library⁹⁸ to employ the Adam algorithm⁹⁷ for minimization of the loss function Eq. (16). The learning rate of the Adam algorithm was set to 0.01. The training data were obtained from the classical MD simulations: the 3888 trajectories for POLI2VS and the 216 trajectories for SPC/E were provided, respectively. The initial values for $\mathbf{z}_{k}$ was set to 0. The optimizations for the system parameters and $\mathbf{z}_{k}$ were then carried out. We chose 22 and 108 trajectories as a minibath for the SPC/E flexible case and the POLI2VS case, respectively. The same trajectories were used for each optimization epoch. The iterations were repeated until the optimization parameters of the potentials converged.

Table 1: Equilibrium OH bond length and HOH angle for SPC/E flexible and POLI2VS force field models.

Model	$r_{0}$	$\theta_{0}$
SPC/E Flexible	0.10282 nm	1.8177 rad
POLI2VS	0.09779 nm	1.8379 rad

5 5. RESULTS AND DISCUSSION

Refer to caption — Figure 1: The learning curve of the mean squared error for $q$ between the predicted trajectories and the actual MD trajectories for the SPC/E flexible (blue curve) and the POLI2VS force fields (orange curve).

The optimized parameters for the potential, the intramolecular mode-mode coupling, and the linear and nonlinear system-bath couplings are listed in Table 2. We employed approximately 50000 epochs for the optimization of the potential parameters, while we needed around 10000 epochs for the optimization of $J_{s}(\omega)$ . We display the learning curves for the mean square of the coordinate, $\Delta q$ , and the mode-mode coupling parameters, $g_{ij}$ , in Figs. 1 and 2. The frequencies of the modes determined from the force constants and the parameters in the HO and MO models are found to be smaller than the values determined from the force field model used in the MD simulations. This is because the LL and SL interactions enhance the system frequency, as shown by analytical and numerical analyses in the harmonic case.^{87, 88, 89, 90, 91, 92} The anharmonic mode-mode coupling constants differ significantly in the HO and MO cases. This is because some of the contributions to the anharmonic mode-mode coupling have already been taken into account by the anhrmonicity of the MO potential, as explained in Sec. 3. The obtained $V_{\mathrm{SL}}/V_{\mathrm{LL}}$ values indicate that the SL interaction plays a larger role in the symmetric-stretch case than in the asymmetric-stretch case, while the effect of the SL interaction is negligibly small in the bend case. This is due to the formation of the hydrogen bonding (HB) between the system molecule and the surrounding molecules, which changes the vibrational frequencies of the O-H1 and O-H2 bonds. The formation of the HB bonds for both the O-H1 and O-H2 bonds occurs in a correlated manner in the symmetric case, while it occurs in an anti-correlated manner in the asymmetric case. Because the two formations occur simultaneously in the correlated case, the symmetric-stretch mode exhibits a larger SL effect. Contrastingly, because the O-H bond lengths do not change significantly in the bending case, the effect of SL coupling is very small.

In Fig. 3, we display the results of the SDF for the symmetric-stretch (orange curve), asymmetric-stretch (blue curve), and bend (black curve) of H₂O molecules evaluated for (a) the SPC/E flexible and (b) POLI2VS cases. It should be noted that the present results for each intramolecular mode with the LL+SL interaction are similar to those obtained in the case with the LL interaction only (not presented). This indicates that the SL interaction is important in the system dynamics, while it plays a minor role for the bath modes.

In each case, the intramolecular symmetric-stretch, asymmetric-stretch, and bend peaks are observed near 3100cm^-1, 3200cm^-1, and 1600cm^-1, respectively, while the intermolecular hydrogen bonded (HB) librational and HB translational peaks are observed in the range 800-400 cm^-1 and near 300 cm^-1, respectively. In each intramolecular mode, we observe the peak at a frequency similar to that of its own mode, because each intramolecular mode can interact with the modes of the surrounding water molecules as an intermolecular interaction.

For the asymmetric-stretch mode (blue curve), the peak near 2200 cm^-1 arises from the coupling between the HB translation and the bend modes, whereas, for the bend mode (black curve), the peak near 3800 cm^-1 and 2600 cm^-1 arise from the coupling between the HB liberation and symmetric-stretch modes and the HB translation and bending modes, respectively. These combination bands indicate the existence of energy transfer pathways among the modes: They are predicted on the basis of analysis of the energy transfer process.⁶⁹ In the case of POLI2VS, we further observed combination bands between the bend and asymmetric-stretch modes and the symmetric- and asymmetric-stretch modes in the ranges 5000-5500 cm^-1 and 7000-7500 cm^-1, respectively. (See Appendix). The assigned intermolecular modes and combination bands are considered in Table 3. In the POLI2VS case, we found that the combination band peak of the HB transitional mode and bending mode is merged with the symmetric-stretch mode peak, and together they appear as a broadened peak. Such features that describe the coupling between the inter- and intra-molecular vibrational modes have not been thoroughly explored on the basis of the infrared and Raman spectroscopic analysis, because they are optically inactive. The present machine learning approach has the capability to estimate optically inactive combination bands that may be different for different force fields and different simulation approaches.

The intermolecular mode frequencies obtained using POLI2VS are higher than those obtained using SPC/E, because POLI2VS was constructed for quantum simulations, while here we used it for classical simulations. For this reason, the appearance of the combinational bands is also different in the POLI2VS case. In order to obtain an accurate description, we must employ quantum trajectories.⁵⁹

Table 2: Optimized potential parameters.

		SPC/E flexible	POLI2VS
	$k_{t}$ [kJ/mol/rad²]	$382.168$	$364.963$
	$k_{r}$ [kJ/mol/nm²]	$311708$	-
	$a$ [nm^-1]	-	$21.1293$
	Bend. [rad^-1]	$-0.0485218$	$0.17656853$
$V_{\mathrm{SL}}/V_{\mathrm{LL}}$	Sym. [nm^-1]	$10.8475$	$18.41797858$
	Asym. [nm^-1]	$-1.73348$	$-0.219222965$
	$q_{\mathrm{bend}}q_{\mathrm{sym}}$ [kJ/mol/rad/nm]	$-217.080$	$-468.367$
$g_{ss^{\prime}}$	$q_{\mathrm{bend}}q_{\mathrm{asym}}$ [kJ/mol/rad/nm]	$-2.03171$	$0.06029$
	$q_{\mathrm{sym}}q_{\mathrm{asym}}$ [kJ/mol/nm²]	$-103.004$	$-516.444$
	$q_{\mathrm{bend}}^{2}q_{\mathrm{sym}}$ [kJ/mol/rad²/nm]	$-212.825$	$-3019.51$
$g_{sss^{\prime}}$	$q_{\mathrm{bend}}^{2}q_{\mathrm{asym}}$ [kJ/mol/rad²/nm]	$101.093$	$100.569$
	$q_{\mathrm{sym}}^{2}q_{\mathrm{asym}}$ [kJ/mol/nm³]	$250917$	$-4511.28$
	$q_{\mathrm{bend}}q_{\mathrm{sym}}^{2}$ [kJ/mol/rad/nm²]	$4123.39$	$64048.4$
$g_{ss^{\prime}s^{\prime}}$	$q_{\mathrm{bend}}q_{\mathrm{asym}}^{2}$ [kJ/mol/rad/nm²]	$-624.620$	$-14091.7$
	$q_{\mathrm{sym}}q_{\mathrm{asym}}^{2}$ [kJ/mol/nm³]	$241830$	$345203$

Table 3: Peak positions for the assigned intermolecular modes and convention bands observed in the bend mode. Results for the intramolecular bending and stretching modes are not presented.

Mode	SPC/E Flexible	POLI2VS
HB stretching mode	85 cm^-1	69 cm^-1
Rotation	557 cm^-1	509 cm^-1
HB trans. + bend.	2623 cm^-1	-
HB lib. + sym.	3759 cm^-1	4170 cm^-1

In the present system-bath model, it is not necessary that the frequency of the system mode be the same as the frequency of the molecular mode in the MD simulations, because the system-bath interaction modifies the central frequency of the system oscillator, as has been found in investigations employing Brownian motion theory³ and optical spectroscopic theory with the LL and SL interactions.^{87, 88, 89, 90, 91, 92} Although the peak heights in the evaluated spectral distribution may depend on the conditions of the optimization, for example, on the value of $\Delta t$ in the HOB simulations, the overall features of the results are insensitive to such details. Verifying the validity of the system Hamiltonian and the system-bath interaction must be done in a case-by-case manner by, for example, evaluating the reaction rates and vibrational spectrum in any given application of interest

6 6. CONCLUSION

We studied a system-bath model to simulate the dynamics of molecular liquids characterized by molecular trajectories obtained from MD simulations. We demonstrated that the machine learning approach is a versatile tool for estimating the model parameters for the system potential, mode-mode interactions, the system-bath coupling, and, most importantly, the SDFs of the bath. An optimized system-bath Hamiltonian evaluated with the present approach allows us to carry out numerically expensive simulations that are not feasible in the case of the full quantum MD simulations, for example for calculations of 2D infrared^{92, 93, 94, 95, 99, 100, 101} and 2D infrared-Raman spectroscopies,^{63, 102, 103} where the quantum dynamics of intramolecular modes play a significant role.

The obtained SDFs can also be used for characterization of the dynamics that play an important role in vibrational spectroscopy and energy translation by taking account of the various combination bands that cannot be explored by spectroscopic means because they are optically inactive. In particular, we discovered the existence of intermolecular interactions between the stretching modes of the system and the bath modes that consist of the stretching modes of the surrounding water molecules.

The difference between the model parameters and SDFs for different MD force fields indicates that we can use such information to verify the description of the MD force fields in a critical manner. Moreover, if we use trajectories obtained from quantum dynamical simulations instead of classical simulations, such as ab initio MD simulations^{104, 105, 106, 107} or quantum MD simulations using a force field designed for quantum simulations^{59, 108, 109, 110, 111}, we should be able to obtain a system-bath Hamiltonian more suitable for quantum dynamical simulations.

Here, we employed the HO and MO Hamiltonians with an anharmonic mode-mode interaction as a trial system. Instead of these, we could also consider a system potential with the same form as that in the force field model used in the MD simulations. Although such a system-bath Hamiltonian is much more difficult to study on the basis of open quantum dynamics theories, we found that such a model can describe the MD results better than the simple HO and MO models (results not presented). In order to model a complex system more accurately, we should include physical observables, such as the optical dipole and polarization, in order to obtain the optimized model Methodologies and model functions developed in various machine learning approches, for example an approach developed for the physics-informed neural network (PINN),^{112, 113, 114} may be useful.

The key feature of the present study is that it treats a system-bath model, in which the primary molecular dynamics are described by a system Hamiltonian, while the other degrees of freedom are described by a harmonic heat bath that is characterized by SDFs. This approach can be applied not only to systems described by molecular coordinates, such as chemical reaction systems^{6, 7, 8, 9} and vibrational systems^{27, 28, 29}, but also systems described by electronic states, such as systems exhibiting non-adiabatic transitions^{10, 11, 12, 13}, quantum devices^{14, 15}, ratchet rectification^{16, 17, 18}, and exciton-transfer systems^{19, 20, 21, 22, 23, 24} in situation that we have microscopic data describing not only molecular coordinates but also electronic states.¹¹⁵ As a future investigation, we plan to extend the present study in such a direction.

{acknowledgement}

The authors would like to thank Dr. Taisuke Hasegawa for providing the source code for the POLI2VS potential model. YT is thankful to professor Shinji Saito for helpful discussions concerning the assignment of the combination bands appearing in the spectral distribution functions. Financial support from HPC Systems Inc. is acknowledged.

Appendix A Appendix: Overtone of OH stretching in the POLI2VS model

Figure 4 illustrates the high-frequency overtone area of $J_{s}(\omega)$ obtained from POLI2VS. (The low-frequency part is presented in Fig. 2(b).) In the SPC/E flexible case, we could not observe these peaks. In the symmetric and asymmetric modes, small peaks appear as overtones of the OH stretching mode. The peak intensities of these overtones become larger if we employ the harmonic potential for the OH bonds instead of the Morse potential. This indicates that the anharmonicity of the MD potential, which cannot be taken into account by the system potential,is described by the overtone peak of $J(\omega)$ .

References

1 Caldeira, A. O.; Leggett, A. J. Path Integral Approach to Quantum Brownian Motion. Physica A: Statistical Mechanics and its Applications 1983, 121 (3), 587–616.
2 Leggett, A. J.; Chakravarty, S.; Dorsey, A. T.; Fisher, M. P. A.; Garg, A.; Zwerger, W. Dynamics of the Dissipative Two-State System. Rev. Mod. Phys. 1987, 59 (1), 1–85.
3 Ulrich Weiss. Quantum Dissipative Systems, 4th ed.; World Scientific: Singapore, 2012
4 Heinz-Peter Breuer; Francesco Petruccione. The Theory of Open Quantum Systems; Oxford University Press: New York, 2002.
5 Tanimura, Y. Stochastic Liouville, Langevin, Fokker-Planck, and Master Equation Approaches to Quantum Dissipative Systems. J. Phys. Soc. Jpn. 2006, 75 (8), 082001.
6 Wolynes, P. G. Quantum Theory of Activated Events in Condensed Phases. Phys. Rev. Lett. 1981, 47 (13), 968-971.
7 Voth, G. A.; Chandler, D.; Miller, W. H. Rigorous Formulation of Quantum Transition State Theory and Its Dynamical Corrections. J. Chem. Phys. 1989, 91 (12), 7749-7760.
8 Hänggi, P.; Talkner, P.; Borkovec, M. Reaction-Rate Theory: Fifty Years after Kramers. Rev. Mod. Phys. 1990, 62 (2), 251-341.
9 Y. Tanimura and P. G. Wolynes, J. Chem. Phys. 96, 8485 (1992).
10 Garg, A.; José Nelson Onuchic; Vinay Ambegaokar. Effect of Friction on Electron Transfer in Biomolecules. J. Chem. Phys. 1985, 83 (9), 4491.
11 Peter G Wolynes. Dissipation, Tunneling, and Adiabaticity Criteria for Curve Crossing Problems in the Condensed Phase. J. Chem. Phys. 1987, 86 (4), 1957.
12 Ikeda, T.; Tanimura, Y. Phase-Space Wavepacket Dynamics of Internal Conversion via Conical Intersection: Multi-State Quantum Fokker-Planck Equation Approach. Chem. Phys. 2018, 515, 203–213.
13 Ikeda, T.; Tanimura, Y. Low-Temperature Quantum Fokker–Planck and Smoluchowski Equations and Their Extension to Multistate Systems. J. Chem. Theory Comput. 2019, 15 (4), 2517–2534.
14 Mason, B. A.; Hess, K. Quantum Monte Carlo Calculations of Electron Dynamics in Dissipative Solid-State Systems Using Real-Time Path Integrals. Phys. Rev. B 1989, 39 (8), 5051-5069.
15 Sakurai, A.; Tanimura, Y. Self-Excited Current Oscillations in a Resonant Tunneling Diode Described by a Model Based on the Caldeira-Leggett Hamiltonian. New J. Phys. 2014, 16 (1), 015002.
16 Peter Reimann; Milena Grifoni; Peter Hänggi. Quantum Ratchets. Phys. Rev. Lett. 1997, 79 (1), 10.
17 Hänggi, P.; Marchesoni, F. Artificial Brownian Motors: Controlling Transport on the Nanoscale. Rev. Mod. Phys. 2009, 81 (1), 387-442.
18 Kato, A.; Tanimura, Y. Quantum Suppression of Ratchet Rectification in a Brownian System Driven by a Biharmonic Force. J. Phys. Chem. B 2013, 117 (42), 13132-13144.
19 Ishizaki, A.; Fleming, G. R. Theoretical Examination of Quantum Coherence in a Photosynthetic System at Physiological Temperature. Proceedings of the National Academy of Sciences 2009, 106 (41), 17255-17260.
20 Strümpfer, J.; Schulten, K. Open Quantum Dynamics Calculations with the Hierarchy Equations of Motion on Parallel Computers. J. Chem. Theory Comput. 2012, 8 (8), 2808-2816.
21 Adolphs, J.; Renger, T. How Proteins Trigger Excitation Energy Transfer in the FMO Complex of Green Sulfur Bacteria. Biophysical Journal 2006, 91 (8), 2778–2797.
22 Kreisbeck, C.; Kramer, T.; Aspuru-Guzik, A. Disentangling Electronic and Vibronic Coherences in Two-Dimensional Echo Spectra. J. Phys. Chem. B 2013, 117 (32), 9380-9385.
23 Kreisbeck, C.; Kramer, T.; Aspuru-Guzik, A. Scalable High-Performance Algorithm for the Simulation of Exciton Dynamics. Application to the Light-Harvesting Complex II in the Presence of Resonant Vibrational Modes. J. Chem. Theory Comput. 2014, 10 (9), 4045-4054.
24 Lee, M. K.; Coker, D. F. Modeling Electronic-Nuclear Interactions for Excitation Energy Transfer Processes in Light-Harvesting Complexes. J. Phys. Chem. Lett. 2016, 7 (16), 3171–3178.
25 Malcolm H. Levitt. Spin Dynamics; John Wiley and Sons, Ltd: New York, 2001.
26 N.G. Van Kampen. Stochastic Processes in Physics and Chemistry; Elsevier: Amsterdam, 1981.
27 Shaul Mukamel. Principles of Nonlinear Optical Spectroscopy; Shaul Mukamel: New York, 1995.
28 Volkhard May; Oliver Kühn. Charge and Energy Transfer Dynamics in Molecular Systems; Wiley-VCH: Weinheim, 2003.
29 Tanimura, Y.; Ishizaki, A. Modeling, Calculating, and Analyzing Multidimensional Vibrational Spectroscopies. Acc. Chem. Res. 2009, 42 (9), 1270-1279.
30 Tanimura, Y.; Mukamel, S. Real-Time Path-Integral Approach to Quantum Coherence and Dephasing in Nonadiabatic Transitions and Nonlinear Optical Response. Phys. Rev. E 1993, 47 (1), 118-136.
31 Okumura, K.; Tanimura, Y. Unified Time-Path Approach to the Generating Functional of the Brownian Oscillator System: The Bilinearly Corrected Feynman Rule for Nonequilibrium Processes. Phys. Rev. E 1996, 53 (1), 214-227.
32 Tanimura, Y.; Okumura, K. First-, Third-, and Fifth-Order Resonant Spectroscopy of an Anharmonic Displaced Oscillators System in the Condensed Phase. J. Chem. Phys. 1997, 106 (6), 2078-2095.
33 Schröder, M.; Kleinekathöfer, U.; Schreiber, M. Calculation of Absorption Spectra for Light-Harvesting Systems Using Non-Markovian Approaches as Well as Modified Redfield Theory. J. Chem. Phys. 2006, 124 (8), 084903.
34 Egorova, D.; Gelin, M. F.; Thoss, M.; Wang, H.; Domcke, W. Effects of Intense Femtosecond Pumping on Ultrafast Electronic-Vibrational Dynamics in Molecular Systems with Relaxation. J. Chem. Phys. 2008, 129 (21), 214303.
35 Montoya-Castillo, A.; Berkelbach, T. C.; Reichman, D. R. Extending the Applicability of Redfield Theories into Highly Non-Markovian Regimes. J. Chem. Phys. 2015, 143 (19), 194108.
36 Makri, N.; Makarov, D. E. Tensor Propagator for Iterative Quantum Time Evolution of Reduced Density Matrices. I. Theory. J. Chem. Phys. 1995, 102 (11), 4600-4610.
37 Makri, N. Blip Decomposition of the Path Integral: Exponential Acceleration of Real-Time Calculations on Quantum Dissipative Systems. J. Chem. Phys. 2014, 141 (13), 134117.
38 Meyer, H.-D.; Manthe, U.; Cederbaum, L. S. The Multi-Configurational Time-Dependent Hartree Approach. Chem. Phys. Lett. 1990, 165 (1), 73-78.
39 Manthe, U.; Meyer, H. -D.; Cederbaum, L. S. Wave-packet Dynamics within the Multiconfiguration Hartree Framework: General Aspects and Application to NOCl. J. Chem. Phys. 1992, 97 (5), 3199-3213.
40 Wang, H.; Thoss, M. Multilayer Formulation of the Multiconfiguration Time-Dependent Hartree Theory. J. Chem. Phys. 2003, 119, 1289.
41 Pullerits, T.; Kühn, O. Exciton-Vibrational Coupling in the Dynamics and Spectroscopy of Frenkel Excitons in Molecular Aggregates. Phys. Rep. 2015, 567, 1-78.
42 Yoshitaka Tanimura; Ryogo Kubo. Time Evolution of a Quantum System in Contact with a Nearly Gaussian-Markoffian Noise Bath. J. Phys. Soc. Jpn. 1989, 58 (1), 101-114.
43 Tanimura, Y. Nonperturbative Expansion Method for a Quantum System Coupled to a Harmonic-Oscillator Bath. Phys. Rev. A 1990, 41 (12), 6676-6687.
44 Ishizaki, A.; Tanimura, Y. Quantum Dynamics of System Strongly Coupled to Low-Temperature Colored Noise Bath: Reduced Hierarchy Equations Approach. J. Phys. Soc. Jpn. 2005, 74 (12), 3131-3134.
45 Tanimura, Y. Reduced Hierarchical Equations of Motion in Real and Imaginary Time: Correlated Initial States and Thermodynamic Quantities. J. Chem. Phys. 2014, 141 (4), 044114.
46 Tanimura, Y. Real-Time and Imaginary-Time Quantum Hierarchal Fokker-Planck Equations. J. Chem. Phys. 2015, 142 (14), 144110.
47 Xu, R.-X.; Cui, P.; Li, X.-Q.; Mo, Y.; Yan, Y. Exact Quantum Master Equation via the Calculus on Path Integrals. J. Chem. Phys. 2005, 122 (4), 041103.
48 Han, P.; Xu, R.-X.; Li, B.; Xu, J.; Cui, P.; Mo, Y.; Yan. Kinetics and Thermodynamics of Electron Transfer in Debye Solvents: An Analytical and Nonperturbative Reduced Density Matrix Theory. J. Phys. Chem. B 2006, 110 (23), 11438-11443.
49 Tang, Z.; Ouyang, X.; Gong, Z.; Wang, H.; Wu, J. Extended Hierarchy Equation of Motion for the Spin-Boson Model. J. Chem. Phys. 2015, 143 (22), 224112.
50 Popescu, B.; Rahman, H.; Kleinekathöfer, U. Chebyshev Expansion Applied to Dissipative Quantum Systems. J. Phys. Chem. A 2016, 120 (19), 3270–3277.
51 Nakamura, K.; Tanimura, Y. Hierarchical Schrödinger Equations of Motion for Open Quantum Dynamics. Phys. Rev. A 2018, 98 (1), 012109.
52 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 H₂O(l) at 25°C between 15,000 and 1 cm^-1. Appl Spectrosc 1996, 50 (8), 1047-1057.
53 Palese, S.; Mukamel, S.; Miller, R. J. D.; Lotshaw, W. T. Interrogation of Vibrational Structure and Line Broadening of Liquid Water by Raman-Induced Kerr Effect Measurements within the Multimode Brownian Oscillator Model ^†. J. Phys. Chem. 1996, 100 (24), 10380-10388.
54 Ji, X.; Ahlborn, H.; Space, B.; Moore, P. B. A Theoretical Investigation of the Temperature Dependence of the Optical Kerr Effect and Raman Spectroscopy of Liquid CS₂. J. Chem. Phys. 2000, 113, 8693.
55 Castner, E. W.; Chang, Y. J.; Chu, Y. C.; Walrafen, G. E. The Intermolecular Dynamics of Liquid Water. J. Chem. Phys. 1995, 102 (2), 653-659.
56 Fecko, C. J.; Eaves, J. D.; Tokmakoff, A. Isotropic and Anisotropic Raman Scattering from Molecular Liquids Measured by Spatially Masked Optical Kerr Effect Spectroscopy. J. Chem. Phys. 2002, 117 (3), 1139-1154.
57 Markus Wendling; Tõnu Pullerits; Milosz A. Przyjalgowski; Simone I. E. Vulto; Thijs J. Aartsma; Rienk van Grondelle; Herbert van Amerongen. Electron-Vibrational Coupling in the Fenna-Matthews-Olson Complex of Prosthecochloris Aestuarii Determined by Temperature-Dependent Absorption and Fluorescence Line-Narrowing Measurements. J. Phys. Chem. B 2000, 104 (24), 5825-5831.
58 Hasegawa, T.; Tanimura, Y. A Polarizable Water Model for Intramolecular and Intermolecular Vibrational Spectroscopies. J. Phys. Chem. B 2011, 115 (18), 5545-5553.
59 Liu, X.; Liu, J. Critical Role of Quantum Dynamical Effects in the Raman Spectroscopy of Liquid Water. Molecular Physics 2018, 116 (7-8), 755-779.
60 Fanourgakis, G. S.; Xantheas, S. S. Development of Transferable Interaction Potentials for Water. V. Extension of the Flexible, Polarizable, Thole-Type Model Potential ( TTM3-F , v. 3.0) to Describe the Vibrational Spectra of Water Clusters and Liquid Water. J. Chem. Phys. 2008, 128 (7), 074506.
61 Hasegawa, T.; Tanimura, Y. Calculating Fifth-Order Raman Signals for Various Molecular Liquids by Equilibrium and Nonequilibrium Hybrid Molecular Dynamics Simulation Algorithms. J. Chem. Phys. 2006, 125 (7), 074512.
62 Ikeda, T.; Ito, H.; Tanimura, Y. Analysis of 2D THz-Raman Spectroscopy Using a Non-Markovian Brownian Oscillator Model with Nonlinear System-Bath Interactions. J. Chem. Phys. 2015, 142 (21), 212421.
63 Ito, H.; Tanimura, Y. Simulating Two-Dimensional Infrared-Raman and Raman Spectroscopies for Intermolecular and Intramolecular Modes of Liquid Water. J. Chem. Phys. 2016,
64 Cho, M.; Fleming, G. R.; Saito, S.; Ohmine, I.; Stratt, R. M. Instantaneous Normal Mode Analysis of Liquid Water. J. Chem. Phys. 1994, 100 (9), 6672-6683.
65 Saito, S.; Ohmine, I. Off-Resonant Fifth-Order Nonlinear Response of Water and CS2: Analysis Based on Normal Modes. J. Chem. Phys. 1998, 108 (1), 240-251.
66 Ohmine, I.; Saito, S. Water Dynamics: Fluctuation, Relaxation, and Chemical Reactions in Hydrogen Bond Network Rearrangement. Acc. Chem. Res. 1999, 32 (9), 741–749.
67 Imoto, S.; Xantheas, S. S.; Saito, S. Molecular Origin of the Difference in the HOH Bend of the IR Spectra between Liquid Water and Ice. J. Chem. Phys. 2013, 138 (5), 054506.
68 Imoto, S.; Xantheas, S. S.; Saito, S. Ultrafast Dynamics of Liquid Water: Frequency Fluctuations of the OH Stretch and the HOH Bend. J. Chem. Phys. 2013, 139 (4), 044503.
69 Yagasaki, T.; Saito, S. A Novel Method for Analyzing Energy Relaxation in Condensed Phases Using Nonequilibrium Molecular Dynamics Simulations: Application to the Energy Relaxation of Intermolecular Motions in Liquid Water. J. Chem. Phys. 2011, 134 (18), 184503.
70 Imoto, S.; Xantheas, S. S.; Saito, S. Ultrafast Dynamics of Liquid Water: Energy Relaxation and Transfer Processes of the OH Stretch and the HOH Bend. J. Phys. Chem. B 2015, 119 (34), 11068-11078.
71 Jeon, J.; Cho, M. Direct Quantum Mechanical/Molecular Mechanical Simulations of Two-Dimensional Vibrational Responses: N -Methylacetamide in Water. New J. Phys. 2010, 12 (6), 065001.
72 Jeon, J.; Cho, M. An Accurate Classical Simulation of a Two-Dimensional Vibrational Spectrum: OD Stretch Spectrum of a Hydrated HOD Molecule. J. Phys. Chem. B 2014, 118 (28), 8148-8161.
73 Shiga, M.; Okazaki, S. An Influence Functional Theory of Multiphonon Processes in Molecular Vibrational Energy Relaxation. J. Chem. Phys. 1998, 109 (9), 3542-3552.
74 Shiga, M.; Okazaki, S. Molecular Dynamics Study of Vibrational Energy Relaxation of CN^- in H₂O and D₂O Solutions: An Application of Path Integral Influence Functional Theory to Multiphonon Processes. J. Chem. Phys. 1999, 111 (12), 5390–5401.
75 Mikami, T.; Okazaki, S. An Analysis of Molecular Origin of Vibrational Energy Transfer from Solute to Solvent Based upon Path Integral Influence Functional Theory. J. Chem. Phys. 2002, 119 (9), 4790-4797.
76 Gottwald, F.; Ivanov, S. D.; Kühn, O. On Computing Spectral Densities from Classical, Semiclassical, and Quantum Simulations. J. Chem. Phys. 2019, 150 (8), 084109.
77 Curtarolo, S.; Morgan, D.; Persson, K.; Rodgers, J.; Ceder, G. Predicting Crystal Structures with Data Mining of Quantum Calculations. Phys. Rev. Lett. 2003, 91 (13), 135503.
78 Dam, H. C.; Ho, Tu Bao; Sugiyama, Ayumu. Simulation-Based Data Mining Solution to the Structure of Water Surrounding Proteins. IJCAI-2011 2011, 2424.
79 Chen, C.; Lu, Z.; Ciucci, F. Data Mining of Molecular Dynamics Data Reveals Li Diffusion Characteristics in Garnet Li₇La₃Zr₂O₁₂. Sci. Rep. 2017, 7 (1), 40769.
80 Ian H. Witten, Eibe Frank, Mark A. Hall and Christopher J. Pal, Data Mining; Practical Machine Learning Tools and Techniques (Morgan Kaufmann, 2016).
81 Bishop, C. M. Pattern Recognition and Machine Learning; Information science and statistics; Springer: New York, 2006.
82 Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. GROMACS: A Message-Passing Parallel Molecular Dynamics Implementation. Comp. Phys. Comm. 1996, 91 (1–3), 43–56.
83 Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1–2, 19–25.
84 Madden, P. A.; Lynden-Bell, R. M. Theory of Vibrational Linewidths. Chemical Physics Letters 1976, 38 (1), 163-165.
85 David W. Oxtoby. Dephasing of Molecular Vibrations in Liquids. Adv. Chem. Phys. 1979, 40, 1.
86 Bader, J. S.; Berne, B. J. Quantum and Classical Relaxation Rates from Classical Simulations. J. Chem. Phys. 1994, 100 (11), 8359-8366.
87 Okumura, K.; Tanimura, Y. Two-Time Correlation Functions of a Harmonic System Nonbilinearly Coupled to a Heat Bath: Spontaneous Raman Spectroscopy. Phys. Rev. E 1997, 56 (3), 2747-2750.
88 Steffen, T.; Tanimura, Y. Two-Dimensional Spectroscopy for Harmonic Vibrational Modes with Nonlinear System-Bath Interactions. I. Gaussian-White Case. J. Phys. Soc. Jpn. 2000, 69 (9), 3115–3132.
89 Tanimura, Y.; Steffen, T. Two-Dimensional Spectroscopy for Harmonic Vibrational Modes with Nonlinear System-Bath Interactions. II. Gaussian-Markovian Case. J. Phys. Soc. Jpn. 2000, 69 (12), 4095–4106.
90 Kato, T.; Tanimura, Y. Vibrational Spectroscopy of a Harmonic Oscillator System Nonlinearly Coupled to a Heat Bath. J. Chem. Phys. 2002, 117 (13), 6221–6234.
91 Kato, T.; Tanimura, Y. Two-Dimensional Raman and Infrared Vibrational Spectroscopy for a Harmonic Oscillator System Nonlinearly Coupled with a Colored Noise Bath. J. Chem. Phys. 2004, 120 (1), 260–271.
92 Atsunori Sakurai; Yoshitaka Tanimura. Does ℏ Play a Role in Multidimensional Spectroscopy Reduced Hierarchy Equations of Motion Approach to Molecular Vibrations. J. Phys. Chem. A 2011, 115 (16), 4009–4022.
93 Ishizaki, A.; Tanimura, Y. Modeling Vibrational Dephasing and Energy Relaxation of Intramolecular Anharmonic Modes for Multidimensional Infrared Spectroscopies. J. Chem. Phys. 2006, 125 (8), 084501.
94 Ishizaki, A.; Tanimura, Y. Dynamics of a Multimode System Coupled to Multiple Heat Baths Probed by Two-Dimensional Infrared Spectroscopy ^†. J. Phys. Chem. A 2007, 111 (38), 9269-9276.
95 Tanimura, Y.; Ishizaki, A. Modeling, Calculating, and Analyzing Multidimensional Vibrational Spectroscopies. Acc. Chem. Res. 2009, 42 (9), 1270-1279.
96 Jo, J.-Y.; Ito, H.; Tanimura, Y. Full Molecular Dynamics Simulations of Liquid Water and Carbon Tetrachloride for Two-Dimensional Raman Spectroscopy in the Frequency Domain. Chem. Phys. 2016, 481, 245–249.
97 Kingma, D. P.; Ba, J. Adam: A Method for Stochastic Optimization. arXiv:1412.6980 [cs] 2014.
98 Abadi, M.; Agarwal, A.; Barham, P.; Brevdo, E.; Chen, Z.; Citro, C.; Corrado, G. S.; Davis, A.; Dean, J.; Devin, M.; Ghemawat, S.; Goodfellow, I.; Harp, A.; Irving, G.; Isard, M.; Jia, Y.; Jozefowicz, R.; Kaiser, L.; Kudlur, M.; Levenberg, J.; Mane, D.; Monga, R.; Moore, S.; Murray, D.; Olah, C.; Schuster, M.; Shlens, J.; Steiner, B.; Sutskever, I.; Talwar, K.; Tucker, P.; Vanhoucke, V.; Vasudevan, V.; Viegas, F.; Vinyals, O.; Warden, P.; Wattenberg, M.; Wicke, M.; Yu, Y.; Zheng, X. TensorFlow: Large-Scale Machine Learning on Heterogeneous Distributed Systems arXiv:1603.04467 [cs] 2016.
99 Fecko, C. J. Ultrafast Hydrogen-Bond Dynamics in the Infrared Spectroscopy of Water. Science 2003, 301 (5640), 1698–1702.
100 Auer, B.; Kumar, R.; Schmidt, J. R.; Skinner, J. L. Hydrogen Bonding and Raman, IR, and 2D-IR Spectroscopy of Dilute HOD in Liquid D₂O. PNAS 2007, 104 (36), 14215–14220.
101 Cowan, M. L.; Bruner, B. D.; Huse, N.; Dwyer, J. R.; Chugh, B.; Nibbering, E. T. J.; Elsaesser, T.; Miller, R. J. D. Ultrafast Memory Loss and Energy Redistribution in the Hydrogen Bond Network of Liquid H2O. Nature 2005, 434 (7030), 199–202.
102 Grechko, M.; Hasegawa, T.; D’Angelo, F.; Ito, H.; Turchinovich, D.; Nagata, Y.; Bonn, M. Coupling between Intra- and Intermolecular Motions in Liquid Water Revealed by Two-Dimensional Terahertz-Infrared-Visible Spectroscopy. Nat. Commun. 2018, 9 (1), 885.
103 Pastorczak, M.; Nejbauer, M.; Radzewicz, C. Femtosecond Infrared Pump–Stimulated Raman Probe Spectroscopy: The First Application of the Method to Studies of Vibrational Relaxation Pathways in the Liquid HDO/D 2 O System. Phys. Chem. Chem. Phys. 2019, 21 (31), 16895–16904.
104 Fanourgakis, G. S.; Xantheas, S. S. Development of Transferable Interaction Potentials for Water. V. Extension of the Flexible, Polarizable, Thole-Type Model Potential ( TTM3-F , v. 3.0) to Describe the Vibrational Spectra of Water Clusters and Liquid Water. J. Chem. Phys. 2008, 128 (7), 074506.
105 Reddy, S. K.; Moberg, D. R.; Straight, S. C.; Paesani, F. Temperature-Dependent Vibrational Spectra and Structure of Liquid Water from Classical and Quantum Simulations with the MB-Pol Potential Energy Function. J. Chem. Phys. 2017, 147 (24), 244504.
106 Schmidt, M.; Roy, P.-N. Path Integral Molecular Dynamic Simulation of Flexible Molecular Systems in Their Ground State: Application to the Water Dimer. J. Chem. Phys. 2018, 148 (12), 124116.
107 Dawson, W.; Gygi, F. Equilibration and Analysis of First-Principles Molecular Dynamics Simulations of Water. J. Chem. Phys. 2018, 148 (12), 124501.
108 Sauceda, H. E.; Chmiela, S.; Poltavsky, I.; Müller, K.-R.; Tkatchenko, A. Molecular Force Fields with Gradient-Domain Machine Learning: Construction and Application to Dynamics of Small Molecules with Coupled Cluster Forces. J. Chem. Phys. 2019, 150 (11), 114102.
109 Sifain, A. E.; Lubbers, N.; Nebgen, B. T.; Smith, J. S.; Lokhov, A. Y.; Isayev, O.; Roitberg, A. E.; Barros, K.; Tretiak, S. Discovering a Transferable Charge Assignment Model Using Machine Learning. J. Phys. Chem. Lett. 2018, 9 (16), 4495-4501.
110 Chmiela, S.; Tkatchenko, A.; Sauceda, H. E.; Poltavsky, I.; Schütt, K. T.; Müller, K.-R. Machine Learning of Accurate Energy-Conserving Molecular Force Fields. Sci. Adv. 2017, 3 (5), e1603015.
111 Stefan Chmiela; Huziel E. Sauceda; Klaus-Robert Müller; Alexandre Tkatchenko. Towards Exact Molecular Dynamics Simulations with Machine-Learned Force Fields. Nat. Commun. 2018, 9 (3887), 1-10.
112 Raissi, M.; Perdikaris, P.; Karniadakis, G. E. Physics-Informed Neural Networks: A Deep Learning Framework for Solving Forward and Inverse Problems Involving Nonlinear Partial Differential Equations. J. Comput. Phys. 2019, 378, 686-707.
113 Tartakovsky, A. M.; Marrero, C. O.; Perdikaris, P.; Tartakovsky, G. D.; Barajas-Solano, D. Learning Parameters and Constitutive Relationships with Physics Informed Deep Neural Networks. arXiv:1808.03398 [physics] 2018.
114 Lu, L.; Meng, X.; Mao, Z.; Karniadakis, G. E. DeepXDE: A Deep Learning Library for Solving Differential Equations. arXiv:1907.04502 [physics, stat] 2019.
115 Pavlo O. Dral; Mario Barbatti; Walter Thiel. Nonadiabatic Excited-State Dynamics with Machine Learning. J. Phys. Chem. Lett. 2018, 9, 5660-5663.