{'fraction_non_alphanumeric': 0.02638030542928648, 'fraction_numerical': 0.01688033096685224, 'mean_word_length': 4.7317037470726, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Raman Spectrum and Polarizability of Liquid Water from Deep Neural Networks

We introduce a scheme based on machine learning and deep neural networks to model the environmental dependence of the electronic polarizability in insulating materials. Application to liquid water shows that training the network with a relatively small number of molecular configurations is sufficient to predict the polarizability of arbitrary liquid configurations in close agreement with ab initio density functional theory calculations. In combination with a neural network representation of the interatomic potential energy surface,the scheme allows us to calculate the Raman spectra along 2-nanosecond classical trajectories at different temperatures for H2O and D2O. The vast gains in efficiency provided by the machine learning approach enable longer trajectories and larger system sizes relative to ab initio methods, reducing the statistical error and improving the resolution of the low-frequency Raman spectra. Decomposing the spectra into intramolecular and intermolecular contributions elucidates the mechanisms behind the temperature dependence of the low-frequency and stretch modes.


Introduction
Raman scattering has been widely used to study rotational and vibrational spectra of gases and condensed-phase systems 1 . The Raman effect arises from the inelastic scattering of visible light with matter, in which incident radiation is shifted in the frequency domain due to vibrations and rotations of the scatterer 2 . Raman spectroscopy thus probes the same frequency region as infrared (IR) spectroscopy, but different selection rules apply for each technique, making Raman and IR complementary tools to investigate the rotational and vibrational signatures of condensed-phase systems.
The Raman spectrum of a system can be modeled from the polarizability time-correlation function obtained from molecular simulations. The sensitivity of the polarizability to the environmental dependence of the electronic structure demands consistent quantum mechanical approaches to model the potential and the polarizability surfaces. Traditional numerical simulations employ ab initio molecular dynamics (AIMD), in which electronic interactions are computed on th 13 e fly using density functional theory (DFT) and the system polarization described through the modern theory of polarization 3,4 . This approach was used to analyze Raman spectra of heavy water, with polarizabilities evaluated at every AIMD step using Density Functional Perturbation Theory (DFPT) 5,6 . However, while AIMD methods with the appropriate functional bring much-needed predictive accuracy, their computational intensity forbids their application on large systems size and time scales 7 . Empirical potentials provide a less expensive alternative, but are less robust and generalizable. One concern is that the optimal parameters for the potential energy surface do not reproduce the correct polarizability surface 8 ; thus, success-ful polarizable models must take care to parameterize the polarizability surface. One such potential, POLI2VS, closely matches the observed IR spectrum, but the low-frequency and librational modes of the Raman spectrum are inconsistent with experiment 9 . Moreover, non-reactive potentials are unable to model systems with mutating chemical environments, such as acidic or alkaline aqueous solutions in which proton transfer occurs on a picosecond timescale 10 . An alternative to empirical potentials is MB-pol, a many-body potential including up to 3-body terms plus induction, parameterized with high-level quantum mechanical calculations on small molecular clusters. This is quite accurate for liquid water and more promising than empirical potentials, as it models the potential and polarization/polarizability surfaces within a self-consistent framework, but, being limited to molecular systems, this method cannot treat dissociation 11,12 . A second consideration in modeling the Raman spectrum is the treatment of nuclear quantum effects (NQEs). Recent studies have approached the quantum TCF in several ways: approximate path integral methods such as centroid molecular dynamics (CMD) and (thermostatted) ring polymer molecular dynamics (TRPMD) 12,13 ; the local monomer approximation and other mixed quantum-classical methods, which treat exactly a small subset of vibrational modes 14 ; and the linearized semiclassical initial value representation (LSC-IVR) 15 . However, each of these approximations has its drawbacks 16,17 , and it is difficult to distinguish effects due to approximate quantum dynamics from those due to the choice for the potential energy surface (PES). Our goal in this paper is to demonstrate the extent to which a classical-nuclei approach using a consistent neural-network-based PES and polarizability surface successfully reproduces experimental results.
AIMD. Some of the authors have also recently used DNNs to learn the electric polarization as a sum of local atomic contributions 27 . This method allows the construction of ab initio-level polarization surfaces, which were used to compute the IR spectra of liquid and crystalline water under different pressure conditions. Alternative machine learning models based on kernel instead of DNN representations have been recently reported in the literature to predict tensorial properties extracted from ab initio calculations [28][29][30] . These models have been used to describe the polarizability and the Raman spectra of isolated molecules and dispersion bound molecular crystals 30 .
The strategy we adopt here is to use models based on DNNs to represent the quantities required by Raman spectra calculations: the interatomic potential energy and force, as well as the polarizability tensor. Upon training with ab initio data, DNN-based simulations reproduce AIMD results at orders of magnitude lower computational cost, enabling accurate simulations of large-scale systems at timescales prohibitively long for AIMD. We use our DNN-predicted interatomic potential and polarizabilities to compute the Raman spectrum of liquid water at different temperatures.
Due to the wide relevance of liquid water for science and technology, experimental studies over the past several decades have investigated Raman spectra of water at different thermodynamic conditions [31][32][33][34][35] . Although most experiments agree on the position of primary and secondary peaks, their interpretation is still not free of controversy. For instance, it is well understood that a decrease in water temperature redshifts and broaden the primary OH stretch peak, a fact attributed to larger stability of hydrogen bonds at lower temperatures. But interpretations diverge on whether and how specific features can be assigned to unique hydrogen bond configurations. Therefore, theoretical modeling can point the way forward in identifying the physical processes at play.

Methods
In this section, we first introduce the theory of Raman spectra, relating the Raman line shape to the autocorrelator of the system polarizability. We then discuss the DNN representation of the Wannier centroid polarizabilities, before providing detail on the DFT calculations used to train the network.

Raman Line Shape from Ab initio Molecular Dynamics
The differential cross section of Raman scattering can be written in terms of the Fourier transform of the time autocorrelation function of the electronic polarizability of the system according to 2 : Here ω is the Raman frequency shift, α α α is the polarizability tensor of the sample,ε s andε i are, respectively, the polarization directions of the scattered and incident light, λ s is the wavelength of the scattered light. The angular brackets denote ensemble average, and the integration is over the time t. Factoring out the dependence on λ s , experiments typically report reduced line shapes R(ω) containing an arbitrary constant factor 36 : (2) Here the Bose Einstein (BE) factor n BE (ω) = 1 − exp(−βhω) is introduced when studying the low-frequency features of the spectrum, which would otherwise be obscured by the Rayleigh line 37 . For fluid systems, it is convenient to decompose the polarizability tensor into a spherical part α = 1 3 Trα α α and a traceless anisotropic tensor β β β = α α α − 1α, yielding the isotropic and anisotropic components of the line shape: In the above formulae, the polarizability depends on the nuclear coordinates, which are quantum mechanical operators but are treated here classically to compute the equilibrium time correlation functions via a molecular dynamics simulation. This amounts to neglecting nuclear quantum effects (NQEs) in the dynamics of the nuclei. In liquid water, NQEs are small but not negligible. Their influence on static equilibrium properties has been quantified in experiments and simulations based on Feynman path integrals 38 , but it is difficult to predict the effect on dynamic properties using statistical simulation methods. Leaving this issue aside, we approximate the environmental dependence of the polarizability with classical mechanics.
Nevertheless, the electronic polarizability itself can only be derived from quantum mechanics. Within the Born-Oppenheimer approximation, the electronic polarizability at time t measures the response of the instantaneous polarization (µ µ µ) of the sample at time t to an infinitesimally small uniform electric field while the nuclear positions are held fixed: We adopt first-principles density functional theory (DFT) to describe the electronic ground state of the system and use the modern theory of polarization to compute µ µ µ. The derivative in Eq. 4 can be expressed analytically with density functional perturbation theory (DFPT) 39 , requiring the solution of the self-consistent response equations for the electrons, or it can be calculated numerically within the electric enthalpy framework 40,41 by applying small but finite electric fields ±ε ε ε to the sample: In practice, the two formulations are equivalent. We adopt here the one based on Eq. 5, which does not require a specialized DFPT code but only a DFT minimization code. The macroscopic polarization µ µ µ of a bulk periodic system is conveniently expressed, modulo a quantum, in terms of the position vectors of the nuclei (r r r i ) and the maximally localized Wannier centers (w w w k ), where e is the unit electronic charge, Z i are atomic numbers, and we have assumed a spin-saturated system 42 . The Wannier centers are obtained from a unitary transformation that minimizes the spatial spread in the occupied orbital subspace 43 . We use a valence-only pseudopotential approach so that the nuclear charges eZ i are the charges of the ions consisting of the nuclei and the frozen core electrons, and the Wannier centers correspond to the valence electrons. Specializing to water, which contains oxygen (r r r O i ) and hydrogen (r r r H m ) ions, the polarization vector is: In water, four Wannier centers can be uniquely associated to their nearest oxygen ion. These four Wannier centers remain close to the same oxygen during dynamical evolution, even when the water molecule to which the oxygen belongs dissociates leading to formation of hydronium and hydroxyl complexes. It is convenient to define a Wannier centroid by the average position of the four Wannier centers associated to oxygen i: Then, the electronic polarizability of the liquid water sample is just the sum of the centroid polarizabilities α α α i : In absence of molecular dissociation, the centroid polarizabilities can be viewed as effective molecular polarizabilities. These are useful for interpreting the spectrum, but only their sum is a physical observable. Its time correlation function yields the Raman line shape through Eq. 3. The polarizability α α α(t) is accessible on the fly in AIMD trajectories, which thus provide a way to compute the Raman spectra using the same DFT approximation for the spectral calculations and for modeling the potential interactions that generate the atomic trajectories. Computational cost restricts these calculations to relatively short trajectories (≈ 100 ps) and small simulation boxes (≈ 100 molecules). The ensuing statistical errors limit the accuracy of Raman spectra calculations for water, particularly at low frequency where the Raman signal is very weak.
Recent progress with machine learning (ML) techniques applied to molecular simulation greatly alleviates these difficulties, as ML interatomic potentials can reproduce accurately AIMD trajectories at a cost that is several orders of magnitude lower than AIMD and scales linearly with system size. Our group has developed a versatile deep neural network (DNN) representation of the potential energy surface of multi-atomic systems called deep potential (DP) that has been used in several applications 18,26 . Recently, the DP representation was generalized to describe the environmental dependence of the polarization µ µ µ 27 . In the next section we discuss how this approach can be extended to the electronic polarizability α α α.

Deep Neural Network for the Environmental Depen-
dence of the Polarizability Tensor Let X (r r r 1 , r r r 2 , ..., r r r i , ..., r r r N ) be an extensive physical property, such as the potential energy U, the polarization µ µ µ, or the polarizability α α α, which depends on the atomic positions. We consider systems for which X can be decomposed into a sum of local components X i that depend on the coordinates of all atoms inside a finite neighborhood N i of the atom at r r r i , i.e.: where r i j = ||r r r i − r r r j || is the distance between i and j, and r c is a predefined cut-off distance. There is no restriction, besides finite range, on the functional form of the environmental dependence in Eq. 9. In a condensed phase, only the global property X is a meaningful observable while the local quantities X i are effective properties that depend on the adopted decomposition. As discussed in the previous section, in water the electronic contribution to the polarization and the electronic polarizability are conveniently given by sums of Wannier centroid contributions. Since each centroid is uniquely associated to an oxygen atom, only the oxygen neighborhoods enter the sum in Eq. 9. By contrast, all atomic neighborhoods, oxygen and hydrogen, are included in the sum when X is the potential energy U. The local quantity is an effective property whose precise value depends on the adopted decomposition. In the case of potential energy, X i denotes the atomic energy U i , whose summation gives rise to the total potential energy U of the system in the DNN model. The analytical negative gradients of U with respect to atomic positions define the interatomic forces. In the case of polarizability, for the liquid water system we consider here, X i becomes the molecular polarizability α α α i , where we only consider i to be oxygen.
The environmental dependence of X i can be accurately represented by a DNN-parametrized function X γ γ γ i , where γ γ γ denotes the parameters of the DNN model. Due to the local dependence of X γ γ γ i on the neighborhood of i, the DNN model is scalable by construction. X γ γ γ i should also satisfy some additional criteria. First, it should depend continuously on the atomic positions and, in terms of efficiency, be orders of magnitude faster than ab initio models. Moreover, it is crucial to preserve the translational, rotational, and permutational symmetry of the quantity that is learned. The potential energy is a scalar quantity and should be invariant upon translation and rotation of the system and identical particle permutation. The polarization and the polarizability are vectorial and tensorial quantities, respectively, and should be translationally and permutationally invariant, but rotationally covariant. The above requirements are achieved by means of two DNNs, an embedding DNN and a fitting DNN. The embedding DNN maps the positions of the atoms belonging to the ith neighborhoods to features that are invariant under symmetry operations, while the fitting DNN maps these features in a way that is symmetry invariant for the potential energy, or covariant for the polarization and the polarizability. The number of hidden layers and outputs is refined in the training procedure.
In detail, we make a local frame transformation to the primed | 3 coordinates, which are defined relative to r r r i : Translational symmetry is preserved by construction in the local frame. It is convenient to use generalized coordinates q q q k , which weight atoms according to their distance from site i, r k ≡ (r r r k · r r r k ) 1/2 , and provide continuous evolution as atoms enter/exit the neighborhood. For this purpose we introduce a smooth weight function s(r ) equal to 1/r k at short distances and decaying smoothly to zero as r k approaches r c , the radius of the neighborhood. The four-component vector q q q k = (q 1 k , q 2 k , q 3 k , q 4 k ) is then given by (s(r k ), s(r k )x k /r k , s(r k )y k /r k , s(r k )z k /r k ), in terms of the Cartesian components of r r r k . We use the matrix Q Q Q = (Q ki ) = (q i k ), which has N i rows and 4 columns, to represent the set of generalized coordinates {q q q k } in a neighborhood.
The embedding DNN is the matrix E E E = (E ik ) = E i (s(r k )) with M rows and N i columns, whose elements are found by training, which maps each element in the set {s(r k ), r k ∈ N i } onto M outputs. Multiplication of E E E by Q Q Q gives the matrix T T T = E E EQ Q Q with M rows and 4 columns, whose generic element is: In Eq. 11, the permutationally invariant sum over the atoms is a smooth function of N i . The last three columns of T T T ( j = 2, 3, 4) transform covariantly under rotation because (Q k2 , Q k3 , Q k4 ) transforms like r r r k . Let S S S be the matrix formed by the first M The elements of D D D are invariant under rotation and permutation. D D D captures the local features of the neighboring pattern of i in a faithful and adaptive way. While a fitting network from D D D to a scalar would properly define U i , one needs different procedures for vectorial and tensorial quantities. The case of the polarization vector was discussed in Ref. 27. Here we specialize to the polarizability tensor. In this case, the output of the fitting network is an M × M diagonal matrix F F F = (F jk ), which is mapped onto α α α i , in a permutationally invariant and rotationally covariant way by right-and left-multiplications with the last three columns of T T T and its transpose: Finally, the parameters γ γ γ are determined by training, i.e., an optimization process that minimizes a loss function, which is here the mean square difference between the DNN prediction and the training data. The Adam stochastic gradient descent method 44 is adopted for the optimization.
The DNN for the polarizability should be combined with a DNN for the PES to study the evolution of the polarizability along MD trajectories. For consistency, the two networks should be trained with electronic structure data at the same level of theory. In this paper, we use a DNN for the polarizability and a DP representation of the PES based on the same DFT data generated with the SCAN functional approximation. Since the ab initio electronic structure data are expensive, efficient learning strategies are crucial. To collect a minimal set of ab initio data for training, we use the iterative learning scheme Deep Potential Generator (DP-GEN) of Ref. 45. In this approach, an ensemble of DNN models, initially trained with a limited set of ab initio data but with different initializations of the network parameters, are used to efficiently explore the configuration space. A small subset of the visited configurations is selected with an error indicator, defined as the variance of the predictions within the ensemble DNNs. The protocol is repeated until all the explored configurations are described with satisfactory accuracy. The error indicator, here chosen to be the maximum standard deviation of atomic forces, exploits the highly non-linear dependence of the DNN models on the network parameters. As a consequence, different initializations of the parameters lead to different local minima in the landscape of the loss function, originating an ensemble of minimizing DNNs. In our experience, good DNN models constructed with the above procedure require significantly less ab initio data in the target thermodynamic range than learning approaches based on independent AIMD sampling data.

DFT Data and Deep Neural Networks for the Potential and Polarizability of Water
The ground state electronic structures of the equilibrium configurations of liquid water within DFT were calculated using the SCAN functional approximation 46 . SCAN predicts with sufficient accuracy the molecular structure of liquid water at equilibrium 47 as well as the dipole moment (SCAN:1.84 D, experiment 48 : 1.86 D) and isotropic polarizability (SCAN: 1.41 Å 3 ,experiment 49 : 1.47Å 3 ) of the water molecule in gas phase. We used the CP code of the Quantum ESPRESSO package 50,51 to compute the electronic ground state, the potential energy, and the forces on the atoms at selected molecular configurations. The same code was also used to compute the Wannier center coordinates and the electronic polarizability at fixed nuclear positions via the electric enthalpy method 40,52 . Norm-conserving pseudo-potentials of Troullier-Martins type 53 were used for both oxygen and hydrogen atoms, and the wavefunctions and charge density were plane-wave expanded with an energy cutoff of 110 and 440 Ry, respectively. Total energy was converged to 10 −8 Hartree or lower. The polarizability was estimated numerically with Eq. 5 using ε = 0.001 a.u., which falls well within the linear response regime. An average error of 0.0005 Å 3 for the DFT polarizabilities was estimated from a higher-order finite difference method based on a 5-point stencil. All Raman response calculations were performed at the experimental equilibrium density of water.
With the DP-GEN scheme we explored a pressure range of 1-10 5 bar with liquid water (192 atoms/cell), ice Ic (192 atoms/cell) and ice Ih (288 atoms/cell). The temperature range of the exploration ranged from 270-370 K for liquid water and 50-270 K for ices Ih and Ic. At the end of the iterative training procedure, our training data contained a set of 2056, 2015 and 1708 con- figurations of liquid water, ice Ih and ice Ic, respectively. The referred training data set included only atomic forces and energy for each atomic configuration, the data needed to train the DNN potential energy surface (DP). The polarizability DNN, on the other hand, was trained only with effective molecular polarizabilities evaluated for the entire DP training set of liquid water. The resulting DP predicts water density (ρ) and diffusion coefficient (D) of 1.07 ± 0.02 g/cm 3 and 0.17 ± 0.01 Å 2 /ps, respectively, in close agreement with SCAN-AIMD (ρ = 1.05 ± 0.03 g/cm 3 and D = 0.19 ± 0.02 Å 2 /ps) 47 .
DP-based Molecular dynamics (DPMD) simulations of liquid water were carried out with 512 water molecules in a periodically repeated cubic cell of 24.9 Å size (H 2 O density of 1.0 g/cm 3 ). The system was initially equilibrated at constant volume coupled to a single Nosé-Hoover thermostat 54,55 for 200 ps. The simulations proceeded at constant volume and energy for 2 ns, the only section of the simulation used to compute the spectra. The classical equations of motion were integrated with the velocity-Verlet algorithm with a time step of 0.25 and 0.5 fs for H 2 O and D 2 O, respectively. All simulations were performed with the Lammps 56 package interfaced with the DeepMD-Kit 57 . The DP-GEN package 58 was used to realize the iterative learning scheme.

Numerical Modeling of Raman Spectra
Given the polarizability tensor as a function of time, we numerically evaluate the classical TCF of the cell polarizability and take its discrete Fourier transform (Eq. 3). Furthermore, we also obtain the contribution of intermolecular coupling to the Raman spectra by decomposing the TCF of the cell polarizability into in-tramolecular and intermolecular terms.
The decomposition of the system polarizability into effective molecular polarizabilities enables us to distinguish spectral features due to autocorrelations (within the same effective molecule) from those due to intermolecular coupling, which weaken with increasing temperature. Analogously to Wan et al. 6 , in evaluating Eq. 3, we first calculated the intramolecular and intermolecular TCF. For the isotropic TCF, these are defined as: where N = 512, and with analogous definitions for the anisotropic TCF.
The decomposition of the cell polarizability into effective molecular polarizabilities also enables us to isolate sources of noise in the intermolecular spectrum. The sum in Eq. 14 was taken only over pairs {i, j} within the first two shells of neighboring molecules (6 Å), as determined by plotting the radial distribution function g OO (r). Including interactions from the third shell and beyond only adds to the zero-frequency component of the Fourier spectrum. The cutoffs were enforced on the coordinates of the oxygen atoms at each initial time used in the computation of the TCF.
Next, the spectrum was obtained by taking a discrete Fourier transform of the TCF and multiplying by the appropriate prefactors: e −iω k t m C(t m )∆t (15) where T = 2.5 ps is the length of the TCF and the times t m were discretized into intervals of ∆t = 0.25 fs for H 2 O and 0.5 fs for D 2 O. The frequencies ω k = 2πk/T , where k is an integer, run from 0 to the Nyquist frequency, discretized into units of 13.33 cm −1 . This provides a bound on the resolution of the spectrum, so we report the locations of peaks to the nearest 10 cm −1 .

Results and discussion
In this section, we present the isotropic and anisotropic Raman spectra computed from DPMD, compare to experiment, and discuss the temperature dependence of the low-frequency and stretch modes. We first show in Fig. 1 the ability of our neural network model to predict effective molecular polarizabilities from ab initio data. The molecular polarizabilities were obtained from a set of 416 liquid water atomic configurations not included in the training data. Our neural network model predicts similar polarizabilities distributions to DFT, with a better agreement for the isotropic polarizability than for the off-diagonal components of the polarizability tensor.

Comparison with Experiment
Our key results are shown in Fig. 2, which compares the DPMD spectra to experimental data obtained from Brooker et al 59   Scherer et al 60 . In order to properly compare the intensities from simulation and experiment, we set the integral of the DPMD isotropic spectrum of water between 2700 and 4000 cm −1 equal to the integral of the same region in the experimental spectrum. The intensities reported for all the other spectra in Fig. 2 (except the experimental low-frequency spectra) are relative to the units defined in Fig. 2a. The experimental intensities of low-frequency region were estimated by comparing the amplified and unamplified intensities of the HOH (or DOD) bending peaks reported by Brooker et al. 59 . Magnification factors of 50× and 40× for the amplified low frequency spectra relative to their unamplified counterparts were crudely estimated for H 2 O and D 2 O, respectively, but this crude estimate has an error of at least 30%, and therefore the experimental intensities shown in Fig. 2c should be taken with reservation. To make a direct comparison, we plot the temperature-reduced spectra at the OH stretch region, R(ω), and the frequency-reduced spectra, R(ω)ω, in the low-frequency region as done in Brooker et al. 59 For H 2 O, we successfully captured the peaks in the isotropic and anisotropic spectra, at 3250 cm −1 and 3470 cm −1 respectively. We are also able to reproduce the accurate location of the peaks in the D 2 O spectrum, at ≈ 2350 cm −1 for the isotropic spectrum and ≈ 2530 cm −1 for the anisotropic spectrum. The quoted wavenumbers of the peaks are defined as the positions, to the nearest 10 cm −1 , of the maxima in R(ω) in the interval 2800 to 4000 cm −1 for H 2 O, and 2000 to 3000 cm −1 for D 2 O. We also observe a shoulder at ≈ 3750 cm −1 of the H 2 O spectrum, blueshifted relative to experiment, which originates form the OH stretch of transiently broken H-bonds. The main shortcoming of our results is the absence of a shoulder at 3400 cm −1 in the H 2 O spectrum, previously assigned to the asymmetric stretch of water. With the exception of LSC-IVR 15 , this shoulder was likewise absent from the theoretical Raman spectra obtained with other methods 6,13 .
Another metric for comparison to experiment is the full width at half maximum (FWHM) of the peaks in the OH and OD stretch, which is sensitive to the choice of functional and the method of approximate quantum dynamics. Several theoretical studies have predicted an stretch band that is either too broad (such as spectra computed using the PBE functional 6,13 ) or too narrow (such as the H 2 O spectra computed using a semiclassical approach 14 and MB-pol 12 ). We find a FWHM of 430 cm −1 (270 cm −1 ) for the isotropic OH (OD) stretch, greater than the experimental width of 418 cm −1 (236 cm −1 ) determined from the temperaturereduced spectrum in Brooker et al. In the anisotropic spectrum, we find widths of 370 cm −1 (286 cm −1 ) for OH (OD), compared to 328 cm −1 (306 cm −1 ) from Brooker's experiment and 296 cm −1 (225 cm −1 ) from Scherer et al.. Possible causes for this artificial broadness are the spurious electronic self-interaction present in the SCAN functional, the missing quantum fluctuations in our treatment of molecular dynamics, or a combination of both effects.
With DPMD, we are able to produce accurate results not only in the stretch band, but across the entire spectrum. We qualitatively reproduce an array of features present in the anisotropic spectra, with approximately correct intensities, over an order of magnitude less than the stretching band. These features include the network modes below 300 cm −1 , examined in greater detail below; the libration mode from ≈ 400 cm −1 to ≈ 800 cm −1 ; the bending peaks at 1670 cm −1 and 1210 cm −1 in H 2 O and D 2 O respectively; and the combination bands at 1640 cm −1 in D 2 O and 2250 cm −1 in H 2 O. The most significant discrepancy between the DPMD-predicted and experimental data is the librational peak, which extends approximately 100 cm −1 higher in wavenumber compared to experiment. This indicates that the libration potential is stiffer than in experiment, which could result from overbonding. Moreover, our results do not capture the significant suppression of the libration peak in D 2 O compared to H 2 O seen in the experimental data. Previous simulations, both those using the harmonic approximation 6 and those employing approximate quantum dynamics 12,13 , have also produced an incorrect intensity in this region.

Temperature Dependence of the OH/OD Stretch
A novel contribution of this paper is that, by running the DPMD simulations at six temperatures from approximately 280 to 370 K, we can observe the temperature dependence of OH/OD stretch. Further, we decompose the spectra into intermolecular and intramolecular contributions, which provides insights into the mechanisms driving the temperature dependence.
The theoretical spectra in Figs. 3a and 3b exhibit a blueshift in the isotropic stretch bands due to weakening hydrogen bonds as the temperature increases. Experimentally, this blueshift manifests as a change in the relative intensities of the symmetric (≈ 3200 cm −1 ) and asymmetric (≈ 3400 cm −1 ) OH stretch peaks. The shoulder is not visible in the DPMD-computed spectra, but we do accurately reproduce the magnitude of the blueshift. Qualitatively, the experimentally observed decrease in maximum intensity is captured for both H 2 O and D 2 O.
While several experimental studies have proposed Gaussian spectral decomposition of the stretch band as a mean for understanding the role of different hydrogen bond configurations, there is no consensus on the number of components 32 , or even whether the multistructure model is the correct starting point 31 . Ab initio modeling offers an alternative route to decomposing the spectrum: by splitting into the intramolecular and intermolecular contributions, we can understand how the relative intensity of couplings between and within molecules varies as the temperature increases. This is shown for the isotropic spectrum of H 2 O in Fig. 3c.
There are several features worth noting in Fig. 3c, which together are responsible for the temperature dependence of the total spectrum. First, the intermolecular spectrum exhibits a maximum in the range of 3190 to 3280 cm −1 , followed by a minimum in the range of 3450 to 3570 cm −1 , both of approximately equal magnitude. At a given temperature, both extrema tend to redshift the frequency at which the total spectrum has a maximum. Second, as the temperature increases, both extrema are blueshifted and reduced in intensity. The same behavior occurs in the intramolecular spectrum, but to a lesser extent: the peak is blueshifted from 3430 to 3520 cm −1 with a slight decrease in intensity. Thus, we can attribute the overall blueshift in the spectrum not only to the separate blueshifting of the two components, but perhaps more importantly, to the reduced role of intermolecular coupling as the temperature rises. We also note that at all temperatures, the intramolecular contribution dominates.
The extent of intermolecular coupling in the Raman spectra of liquid water can be controlled experimentally through isotopic substitution. Raman spectra of dilute HOD in D 2 O, for instance, is able to probe the OH stretch almost uncoupled from the surrounding oscillators. The temperature dependence of both experiment and our DPMD simulations of HOD in D 2 O (Fig. 4) Fig. 4 we also observe a good agreement between DPMDpredicted and the experimental spectra, although DPMD predicts broader peaks. The experimental spectra also contains a shoulder at ≈ 3600 cm −1 , usually attributed to the oscillation of non-Hbonded OH groups. This shoulder is blueshifted in the DPMDpredicted spectrum at 370 K.

Temperature Dependence of the Low Frequency Spectrum
Turning now to the low-frequency regime, the anisotropic spectrum of H 2 O below 300 cm −1 is shown in Fig. 5. At 300 K, we find the nominal 60 cm −1 and 180 cm −1 peaks whose presence has been firmly established in experiment 34,61 . Walrafen et al. further observe that as the temperature increases, the 180 cm −1 peak decreases in intensity, disappearing entirely in the gas phase. This temperature dependence is evident, but slightly less pronounced in the DNN-predicted spectra relative to experiment. . Experimental spectra were multiplied by a constant factor in order to match maximum intensities of lowest temperature DPMD and experimental spectra. Experimental data taken from Walrafen et al. 34 .
While the two low-frequency peaks are experimentally attested, their origin is not entirely clear. Previous ab initio simulations of the IR spectrum of water and ice have assigned these peaks to H bond network modes: bending of H bonded oxygen atoms (≈ 70 cm −1 ) and stretching of H bonded oxygen atoms (≈ 200 cm −1 ) 62 . Another school of thought assigns both peaks to restricted translational modes, modeling the 180 cm −1 mode as a harmonic oscillator which changes from underdamped to overdamped as the temperature increases 63 .
In an AIMD study of the low-frequency D 2 O spectrum, Wan et al. attribute the 60 cm −1 peak to intramolecular dipole induceddipole modes and the 200 cm −1 peak to intermolecular charge fluctuations 6 . The DNN-predicted spectra in Fig. 6, obtained from a more accurate functional, larger system size, and longer trajectories than Wan et al, are consistent with these conclusions. At 300 K, the nominal 60 cm −1 peak is dominated by the maximum in the intramolecular spectrum at 70 cm −1 , with a much weaker, redshifted contribution from the intermolecular spectrum at 30 cm −1 . In addition, while the intramolecular spectrum exhibits a shoulder at 190 cm −1 , a more distinct peak appears in the intermolecular spectrum, at 200 cm −1 .
The division into intramolecular and intermolecular components also sheds some light on the temperature dependence of the peaks. As the temperature rises, the 60 cm −1 peak in the intramolecular spectrum is slightly redshifted, and the shoulder at 210 cm −1 redshifts and disappears. As in the OH stretch, the intermolecular contribution is far more sensitive to temperature; the nominal 180 cm −1 peak redshifts from 210 cm −1 at 281 K to 160 cm −1 at 363 K, dropping in intensity by a factor of 0.6. By accessing long time scales and controlling noise, DPMD enables us to accurately model this intermolecular contribution to the spectrum, which is crucial in explaining the temperature dependence at both low and high frequencies.

Conclusion
In this paper, we have constructed a DNN representation of the physical properties of a molecular system, focusing specifically on the effective molecular polarizability. By doing so, we were able to model the temperature dependence of the Raman spec-trum, a useful tool for examining the local structure of condensed phase systems. Liquid water offers a particularly interesting testing ground for the DNN framework, as its complex hydrogen bond networks lead to unique vibrational signatures which have thus far defied full theoretical explanation. The DPMD-predicted trends in the spectra for H 2 O and D 2 O at temperatures ranging from 280 to 360 K are in good agreement with experiment, although nuclear quantum effects must be included to complete the picture. We emphasize, however, that while the Raman spectra presented here are for water, the DNN framework is fully generalizable to other systems of interest in chemical physics. Since DPMD greatly reduces the uncertainty due to limited statistics in AIMD, it offers a useful starting point for assessing the systematic errors of DFT functionals and approximate quantum corrections relative to experiment. Although the SCAN functional describes H-bonds and intermediate van der Waals forces significantly better than GGA functionals, it is still affected by self-interaction errors. Going forward, a study of the effect of using different functional approximations in the Raman spectra will be useful. Of particular interest is the effect of functional approximations, which reduce the self-interaction errors. It will also be worthwhile to study how quantum corrections affect the Raman spectra. The ordersof-magnitude efficiency gains of DPMD over AIMD should enable various semi-classical methods and even permit the analytic continuation of imaginary time data.

Conflicts of interest
There are no conflicts to declare.