Tensorial Properties via the Neuroevolution Potential Framework: Fast Simulation of Infrared and Raman Spectra

Infrared and Raman spectroscopy are widely used for the characterization of gases, liquids, and solids, as the spectra contain a wealth of information concerning, in particular, the dynamics of these systems. Atomic scale simulations can be used to predict such spectra but are often severely limited due to high computational cost or the need for strong approximations that limit the application range and reliability. Here, we introduce a machine learning (ML) accelerated approach that addresses these shortcomings and provides a significant performance boost in terms of data and computational efficiency compared with earlier ML schemes. To this end, we generalize the neuroevolution potential approach to enable the prediction of rank one and two tensors to obtain the tensorial neuroevolution potential (TNEP) scheme. We apply the resulting framework to construct models for the dipole moment, polarizability, and susceptibility of molecules, liquids, and solids and show that our approach compares favorably with several ML models from the literature with respect to accuracy and computational efficiency. Finally, we demonstrate the application of the TNEP approach to the prediction of infrared and Raman spectra of liquid water, a molecule (PTAF–), and a prototypical perovskite with strong anharmonicity (BaZrO3). The TNEP approach is implemented in the free and open source software package gpumd, which makes this methodology readily available to the scientific community.


I. INTRODUCTION
Infrared (IR) and Raman spectroscopy are widely used techniques for the non-destructive characterization of the dynamics and to some extent chemistry of materials spanning the entire range from the gas phase to condensed matter [1][2][3].Over the years, various theoretical approaches have been developed for simulating IR and Raman spectra, including in particular methods based on ab-initio molecular dynamics (MD) simulations [4][5][6][7][8].While these approaches are capable of reproducing experimental IR and Raman spectra of gases, liquids and solids [5,[7][8][9], they are severely limited with respect to the system sizes and time scales attainable for two main reasons [5,10]: Firstly, ab-initio MD simulations rely on computationally demanding electronic structure calculations that scale strongly with system size in order to obtain energy and forces at every time step.Secondly, similarly expensive calculations of dipole moment (µ), polarizability (α) or electric susceptibility (χ) are required for at least many thousand configurations to achieve numerical convergence of the underlying correlation functions [5].
MD simulations can be accelerated by using classical force fields [11][12][13] or empirical interatomic potentials [14,15], which approximate the potential energy surface (PES) with physically motivated yet constrained functions and few fitted parameters.The accuracy of such approaches for general materials is, however, often limited, negatively affecting the prediction of IR and Raman spectra [16].Machine learning (ML) potentials are well suited to address this challenge as they bridge between the accuracy of quantum mechanical methods and the computational efficiency of classical force fields or empirical interatomic potentials [17][18][19][20][21].The power of this approach, in particular for capturing vibrational properties of materials has been shown repeatedly, see, e.g., Refs.22-26.The calculation of µ, α or χ can be accelerated using parametric models in similar fashion.Considering only static charges, the dipole moment is given by µ = N i=1 Q i r i , where Q i and r i are the charge and position of atom i.Many classical force fields [11][12][13] assign fixed charges to atoms and thereby provide a convenient approach for calculating µ.Such fixed-charge models neglect, however, polarization effects, which can lead to large errors [27].While this situation can in principle be ameliorated by fluctuating-charge models [28,29], the latter tend to lack robustness and can be difficult to generalize [10,30].
Both α and χ describe the dielectric response to an applied electric field.For α or χ, the bond polarizability model is one of the most frequently used parametric ones, and has for example been applied to alkanes [31,32], zeolites [33] as well as carbon nanotubes [34].However, this simple model often suffers from unsatisfactory transferability when used in different environments [35].POLI2VS [36] and MB-pol [37] are two other parametric models that can be used for predicting µ and α, but are limited to molecular systems such as water [10].
The successful applications of ML potentials have inspired the development of ML dipole, polarizability, and susceptibility models [22,[38][39][40][41].For µ, a rank-1 tensor, both partial-charge and the partial-dipole ML models have been developed [30].The objective of the partialcharge models is to assign proper partial charges for atoms in order to fit the total dipole moment [22,30,42].
Here, one concern is the balance between the fitting quality of µ and the reproducibility of total charges [22,30].By contrast partial-dipole models such as symmetryadapted Gaussian process regression (SA-GPR) [38], tensorial embedded-atom neural network (T-EANN) [39], and deep potential (DP) [40] treat µ as a sum of vectors [30,38] that can be determined from atom-centered chemical environments.
While this approach works for µ, which is a rank-1 tensor, it does not transfer to the construction of ML models for α or χ, which are rank-2 tensors.This has motivated the pioneering development of the SA-GPR method for tensorial properties [38] as well as later the T-EANN [39,43] and DP models [10].
The combination of ML potentials with ML models for µ, α or χ enables the simulations of IR and Raman spectra.This approach has been used to predict, e.g., the IR spectra of methanol, n-alkanes, and a peptide [22], IR and Raman spectra of liquid water [10,20,39,44] or the Raman spectra of various solid materials [45].While these earlier studies have established the usefulness of ML models for predicting IR and Raman spectra, there is still ample room for improvement of current models for µ, α or χ in terms of computational and data efficiency [30,39] as well as the accessibility of these techniques in order to lower the threshold for the widespread adoption of such approaches.This situation motivates the present work, in which we introduce accurate as well as computationally and data efficient ML models for rank-1 and rank-2 tensors based on the NEP framework [21,46,47].We demonstrate the efficacy and efficiency of the resulting TNEP approach by training models for µ, α, and χ, and combining these with models for the PES to predict IR and Raman spectra for a molecule (PTAF -), a liquid (water), and a solid (BaZrO 3 ; Fig. 1).We make this methodology available via the gpumd package [47], enabling comprehensive simulations of high-quality IR and Raman spectra with limited user effort.

MD trajectory
Frame 1 Frame 2 Frame N

Autocorrelation functions
Infrared spectrum e.g., liquid water in Figure 4 Time Raman spectrum e.g., BaZrO 3 in Figure 7 Structures

MD simulation
< l a t e x i t s h a 1 _ b a s e 6 4 = " 2 T 9 f 0 s

II. METHODOLOGY
A. NEP models for the PES Since the ML models for µ and α that we introduce below are based on the NEP framework for modeling PESs [21,46,47], we first provide a brief review of the latter.Originally NEPs are ML potentials that model the highdimensional PES of finite or extended systems, in the spirit of the neural network potential model proposed by Behler and Parrinello [48].In this formalism, the total energy of the system is given by the sum of atomic site energies U = i U i .The site energy U i for a given atom i depends on the local environment of the atom, which is represented by an abstract vector q ν i with a number of components indexed by ν.The function mapping from the descriptor to the site energy is represented by a feedforward neural network (also known as a multilayer perceptron) with typically a single hidden layer.The input layer of the neural network is thus the descriptor vector and the output layer consists of a single node whose value is the site energy U i of the considered atom i, which can be formally expressed as From the energy, we can derive the rank-2 virial ten-sor that serves as the foundation for the dipole and polarizability models developed in the present work.For a given structure with N atoms, the virial tensor can be expressed as [47] where r υ ij is the υ-component of the vector r ij ≡ r j − r i , and r i is the position of atom i.We refer to the term ∂U i /∂r ν ij as the partial force, explicit expressions for which have been presented in the original works developing the NEP approach [21,47].

B. TNEP rank-1 tensor models
To develop a ML model for predicting µ, we first note that it is a rank-1 tensor commonly expressed as a vector, in contrast to the energy, which is a rank-0 tensor (i.e., a scalar).The partial force in Eq. ( 2) is a vector, but the summation of it over the whole structure would be zero as a result of Newton's third law.To obtain a vector representation that does not vanish for a general structure, we note that the quantity defined in Eq. ( 2) is a rank-2 tensor that can adopt both positive and negative values (as it is the virial tensor in the context of PES models).We can thus obtain an expression for a vector quantity by contracting this rank-2 tensor with a vector.A natural choice for the vector to be contracted is r ij , which yields the following expression for rank-1 tensors such as the dipole moment where r 2 ij = r ij •r ij is the distance squared between atoms i and j.We note that U i here should have the dimension of charge instead of energy.Crucially this goes to show that the NEP formalism for PESs can be directly used to construct a ML model for rank-1 tensors such as the dipole moment.Below we refer to Eq. ( 3) as the TNEP dipole model.

C. TNEP rank-2 tensor models
To develop ML models for predicting α or χ, we first note that these are rank-2 tensors.Clearly, the quantity defined in Eq. ( 2) is an ideal candidate.However, using only Eq. ( 2) to represent α or χ does not lead to high regression accuracy because the diagonal terms of α or χ are usually much larger than the off-diagonal ones.We therefore represent α (and equivalently χ) as a combination of Eqs. ( 1) and ( 2) as follows where δ υν is the Kronecker delta.Note that both the first and second term on the right-hand side contribute to the diagonal elements of α υν , but only the second term contributes to the off-diagonal elements.U i here has the dimension of polarizability instead of energy, yet the entire NEP formalism can be reused.Below we refer to Eq. ( 4) as the TNEP polarizability or susceptibility model.

D. Loss functions
The NEP approach is named after the underlying ML model (a neural network) and the separable natural evolution strategy used as the training algorithm [49].The latter is a principled real-valued black-box optimization method that is very well suited for training the weight and bias parameters in the neural network, of which there are typically a few thousand.The optimization is driven by the minimization of a loss function that is given by the weighted sum of the root-mean-square errors (RMSEs) of physical quantities as well as L 1 and L 2 regularization terms.For the construction of PES models, the physical quantities included in the loss function are the energies, forces, and virial tensors of the structures in the training set, where ∆U (z), ∆F (z), and ∆W (z) are the RMSEs of energies, forces, and virials calculated using a set of trainable parameters z, and λ e , λ f , and λ v are the corresponding relative weights.Explicit expressions for the regularization terms can be found in Ref. 47.For the construction of dipole TNEP models, the loss function is defined in terms of the RMSE of the dipole ∆µ(z) For the construction of polarizability TNEP models, the loss function is defined in terms of the RMSE of the polarizability ∆α(z)

E. Dielectric response
It is instructive to recall some relations that describe the response of finite (such as molecules) and extended systems (such as solids and liquids) to an applied electric field.
If a molecule is subjected to an electric field E the resulting displacement of nuclei and electrons induces a dipole, which is given by [50] µ ind = αE, where α is the molecular polarizability.
For an extended system such as a solid or a liquid, one considers equivalently the dipole moment per unit volume, i.e., the polarization where χ is the electric susceptibility.In the context of bulk liquids the latter has also been referred to as the bulk polarizability.For clarity in the following, we use the term polarizability only to refer to the molecular polarizability.There are different conventions for expressing µ, α, and χ leading to different units (Sect.S7).Here, we use e • bohr for µ and bohr 3 for α whereas χ is unitless.
We note that under certain conditions, one can approximately connect the molecular polarizability and the electric susceptibility via the Clausius-Mossotti relation, which is based on a mean-field treatment of local field effects (see Sect.S8 in the Supporting Information).

F. The IR intensity
The IR absorption cross section is given by [50] where n is the refractive index of the material, c the speed of light, β = 1/k B T and M (ω) is the absorption lineshape given by the Fourier transform of the autocorrelation function (ACF) of the (total) dipole moment µ, where ⟨• • • ⟩ indicates the average over time origins and ε is the polarization of the light [50].For an isotropic sample, the time correlation should be averaged over the three directions, i.e., the lineshape reduces to one third of the trace of the dipole time correlation.Since the lineshape is sampled classically, we make a classical approximation for the prefactor by expanding the Boltzmann factor to first order, which gives G. The Raman intensity The differential Raman cross-section for Stokes scattering is given by [50-52] where n is the polarization of observed light, ε is the polarization of the incoming light, and Ω is a solid angle.
Here, it is assumed that the frequency of the incoming light ω in is significantly larger than the Raman shift ω, and significantly smaller than the band gap, i.e., far from any electronic excitations.L(ω) is the Raman lineshape given by the Fourier transform of the time-dependent polarizability α(t) (finite systems) or susceptibility χ(t) (extended systems), e.g., in the case of the former Note that the elements of the polarizability (or susceptibility) tensor are selected by the polarization of the incoming and outgoing light as indicated in Eq. (10).
Polarized Raman measurements can be directly related to Eq. ( 10) by combinations of the Raman lineshape L(ω).One can also calculate an average spectrum for isotropic samples [50].The polarizability tensor (and equivalently the susceptibility tensor) can also be written as α = γI + β where γ = Tr(α)/3 and β is a traceless tensor to obtain the isotropic (polarized) and anisotropic (depolarized) spectrum.This leads to the decomposition The electric susceptibility (Sect.II E) can be separated into an electronic and an ionic contribution where the general frequency dependence of these terms is emphasized.For the prediction of Raman spectra we only need to consider the electronic contribution χ e (ω).Furthermore, we limit ourselves to non-resonant Raman spectroscopy.This means that we require the electric susceptibility in the ion-clamped static limit, i.e., χ e (0), and do not have to consider the frequency dependence of χ e (ω), which arises from electronic transitions.

H. Workflow for simulations of IR and Raman spectra
By combining a NEP model for the PES with TNEP models for dipole, polarizability or susceptibility, one obtains a simple yet general workflow for the computation of IR and Raman spectra (Fig. 1).Starting from a NEP PES model, large-scale MD simulations are performed to sample the PES via the gpumd package, typically for a few hundred picoseconds.TNEP dipole, polarizability or susceptibility models are then employed to predict µ(t), α(t) or χ(t) along the trajectory.Finally, IR or Raman spectra are obtained via Fourier transformation of the respective ACFs via Eqs.( 9) or (10).[39] was reported for the averaged molecular polarizability obtained via the Clausius-Mossotti relation (Sect.S8).The validation RRMSEs for χ/ρ should be somewhat higher than that for the averaged molecular polarizability (also see Table S5).

III. PERFORMANCE EVALUATION
In this section, we evaluate the performance of TNEP dipole, polarizability, and susceptibility models in comparison with models from the literature with respect to both regression accuracy and computational speed.The comparison includes the molecules H 2 O, (H 2 O) 2 , and H 5 O 2 + (the Zundel cation), as well as a set of configurations representing liquid water.Structures with dipole, polarizability, and/or susceptibility data were retrieved from the repository maintained by the developers of the SA-GPR models [38,53] (see Sect.S1 in the Supporting Information for details).The data set for each of these systems comprises 1000 configurations, half of which were use for training, while the other half were used for validation.The hyperparameters used in the training of the TNEP models are presented in Tables S1 and S2.In the case of the SA-GPR method, the results for liquid water were computed using a publicly available model [54] while the models for the molecules were trained by us (see Sect.S3 for details).In the case of the T-EANN method, we only use those data available in the literature [39].

A. Dipole moment
The TNEP dipole models can achieve very high precision when predicting µ for both molecules and liquid water with very low RMSEs (Table I) and coefficients of determination (R 2 ) very close to one (Fig. S2).
As a further, more intuitive measure, one can also consider the root-mean-square-error relative to standard deviation (RRMSE) [39], defined as the RMSE divided by the standard deviation of the reference data (Fig. 2a).For the water monomer (H 2 O) all three methods yield extremely small RRMSEs below 0.1%.For the other three systems, including liquid water, the TNEP and SA-GPR models achieve comparable accuracy while the T-EANN models perform systematically worse.This behavior is particularly pronounced for liquid water and might arise since the T-EANN model uses the positions relative to the center of mass as input, which are not well defined in periodic systems [55,56].
Neutral molecules.The µ of neutral molecules such as H 2 O or (H 2 O) 2 is uniquely defined.In the TNEP approach µ is calculated by summing over atomic contributions which, by contrast to, e.g., the T-EANN approach, does not require choosing a reference point.Therefore, the TNEP dipole models are naturally suitable for neutral molecules.
In this context, we note that we also trained and validated a model for the QM7B data set containing thousands of neutral organic molecules [57,58], for which we make similar observations (Sect.S4).The TNEP model yields a very low RMSE for the validation set of 1.80 × 10 −3 e • bohr/atom and a very high R 2 score for the validation set of about 0.998.
Charged molecules.The µ of charged molecules is non-unique and depends on the choice of the reference point [4,59].For charged molecules, one should therefore employ the relative permanent dipole µ r defined with respect to the center of mass, when training TNEP dipole models.The reference µ in the H 5 O 2 + data set [38,53] have already been transformed to µ r .Therefore, the absolute dipole moment of H 5 O 2 + including the movement of the center of mass should then be µ = µ r + e • r COM .The same procedure has been applied to the PTAF -molecule below (Sect.IV C).
Periodic systems.Traditional methods for calculating µ cannot be applied to periodic systems since the position operator is not uniquely defined [56,60].This issue is overcome via the modern theory of polarization [38,60,61], which provides a rigorous definition for the polarization of periodic systems and established a methodology for calculating µ.It was therefore used in the present work to obtain µ for periodic systems includ-

B. Polarizability and susceptibility
The RMSEs for the diagonal and off-diagonal elements of α of (H 2 O), (H 2 O) 2 and H 5 O 2 + are quite small (Table II), indicating the high accuracy of the TNEP polarizability model.The coefficients of determination are larger than 0.98 mirroring this trend (Fig. S7 and Fig. S8).For liquid water, we consider χ/ρ, which has the same unit as the polarizability per atom.The RM-SEs for χ/ρ are on the same order of magnitude as the RMSEs for α (Table II).
The NEP models achieve an accuracy that is comparable to the T-EANN and SA-GPR models for the polarizability of (H 2 O) 2 and H 5 O 2 + as well as the susceptibility of liquid water (Fig. 2b).While the performance for the water monomer H 2 O is worse, the TNEP model still yields a validation RRMSE of less than 1%.
As a further test we constructed a TNEP polarizability model for the QM7B data set (Sect.S4).The RMSE values for the validation set are 4.64×10 −2 bohr 3 atom −1 and 2.58 × 10 −2 bohr 3 atom −1 for the diagonal and offdiagonal elements of α, respectively.For comparison, Wilkins et al. [62] reported a higher RMSE value of 5.50 × 10 −2 bohr 3 atom −1 over both the diagonal and offdiagonal elements of α using a SA-GPR model.

C. Computational speed
It is now instructive to evaluate the computational performance of TNEP models in comparison with publicly available SA-GPR models [53,54].To this end, we consider liquid water systems with varying numbers of atoms.Starting from a cell containing 96 atoms, larger samples with up to 69 984 atoms were created by replication.
The SA-GPR models can only be run serially on a central processing unit (CPU).In contrast, the TNEP model can be run on CPUs using NEP CPU [63], e.g., via the interface provided by the calorine package [64], or on Here, the SA-GPR results were obtained using the TENSOAP-FAST implementation [54].
graphics processing units (GPUs) by using the gpumd package.The SA-GPR and TNEP (CPU) models were tested on a server containing two Intel XEON Platinum 8275CL processors with a system memory of 256 GB, while the TNEP (GPU) models were tested on a heterogeneous server containing two Intel XEON Gold 6148 processors and an Nvidia GeForce RTX 4090 card with a graphics memory of 24 GB.
The comparisons show that for system sizes ≳ 1000 atoms the TNEP CPU models are at least one order of magnitude faster than the SA-GPR models on CPUs for both dipole and polarizability (Fig. 3).On CPUs the TNEP models exhibit nearly perfect weak scaling over the system sizes considered here.In contrast, the SA-GPR models show a notable decrease in speed as the system size increases.Running the TNEP models on GPUs enables an additional speed up by an order of magnitude or more.For very small systems the GPU implementation is limited by IO.In addition we note that gpumd allows one to evaluate TNEP models on-the-fly during MD simulations for prediction of tensorial properties with a small impact on simulation speed (Sect.S10).

IV. APPLICATIONS
Having established the accuracy and computational performance of the TNEP approach by comparison with reference data sets, we now demonstrate the application of NEP and TNEP models in combination for predicting IR and Raman spectra of molecules, liquids, and solids.To this end, we employ the correlation function approach outlined above (Sect.II H and Fig. 1).

A. IR spectrum of water
Firstly, we developed a NEP PES model for liquid water using energy, atomic forces, and virial data from density functional theory (DFT) calculations (Sect.S2).
Next, a system of 216 water molecules was equilibrated in the NPT ensemble for 100 ps using the trained PES model at 298 K and 1 bar, followed by a further equilibration run in the NVT ensemble for another 100 ps.Three production runs were carried out in the NVE ensemble for a duration of 200 ps.A time step of 0.5 fs was used throughout.We note that quantum effects can be actually rather pronounced in water as has been shown by path integral MD simulations in, e.g., Refs.65-67.Here, we, however, decided to carry out classical MD simulations in order to enable a one-to-one comparison with the results of earlier studies.
The time dependence of the dipole (µ(t)) was computed for the production trajectories with a spacing of 1 fs using the TNEP dipole model for liquid water described above (Sect.III A).The IR spectrum was then obtained by Fourier transforming the dipole moment ACF via Eq.( 9).The final IR spectrum was obtained by averaging the IR spectra from the production runs.
For comparison, we also ran a 200 ps MD simulation with the TIP3P force field [68] via the CP2K software package [69], where the TIP3P force field uses charges of −0.834 e and 0.417 e for oxygen and hydrogen, respectively.
The NEP-TNEP method yields an IR spectrum that is in very good agreement with experimental data [70,71] over the entire frequency range from 0 to 4000 cm −1 (Fig. 4a).This includes the hydrogen-bond stretching band [10] between 160 and 250 cm −1 , the libration band [10] from 400 to 800 cm −1 associated with hindered molecule rotations [37], the bending modes [37,72] at about 1650 cm −1 as well as the OH stretching band [37,72] from 2800 to 4000 cm −1 .The NEP and TNEP models for PES and µ in conjunction with the underlying exchange-correlation functional thus succeed in capturing the entire range stretching from the soft intermolecular to the stiff intramolecular modes.This performance is also observed for the DP model (Fig. 4a).
By comparison classical models produce rather large errors for the location of several features in the IR spectrum of water.MD simulations with classical force fields [68,73] such as TIP3P (Fig. 4a) and SPC/E tend to predict a blue-shifting of the bending modes by roughly 100 to 200 cm −1 .A similar tendency was also observed for the POLI2VS model [36].The results from the MB-pol model on the other hand exhibit a blue-shift of the OH stretching band by about 50 cm −1 [37].
The width of the OH stretching band has been proven to be quite difficult to predict due to the anharmonicity of the OH stretch mode [37].The NEP-TNEP approach yields a value of 380 cm −1 for the full width at half maximum of this band, which is in good agreement with experimental estimates of about 350 cm −1 from Downing's experiment [70].Both NEP-TNEP and DP predictions exhibit a slight high-frequency tail for this band, which is not visible in the experimental spectra.This small difference could originate from the strongly constrained and appropriately normed (SCAN) functional [76] that was used for generating the PES training data [10,77] and/or the absence of quantum effects in the (classical) MD simulations [10,37].Brooker et al. [74], and Morawietz et al. [75].Simulated spectra from T-EANN [39], MB-pol [37], and DP [10,44] models were adapted from the literature.In (a) and (b) the spectra were normalized by the integral between 80 and 2500 cm −1 , while in (c) they were normalized by the integral between 1000 and 2500 cm −1 .

B. Raman spectra of water
To obtain the Raman spectra of liquid water we sampled the time dependence of χ(t) using the TNEP susceptibility model and subsequently computed the ACFs for the same trajectories used for the prediction of the IR spectra.The full spectrum given by Eq. ( 10) and averaged over the available trajectories was then split into isotropic (polarized) and anisotropic (depolarized) contributions via Eq.(12).
The anisotropic spectrum predicted by the NEP-TNEP approach is overall in very good agreement with experimental data (Fig. 4b) [74,75].The locations of peaks and relative intensities of the stretching, bending, and librational modes in the simulated anisotropic Raman spectra are all well produced.It is noteworthy that in the low frequency region below approximately 1000 cm −1 , the variation between the experimental spectra is larger than the variation between the ML models and the experimental data.This could be related to difficulties associated with processing the experimental raw data in this frequency region.
The T-EANN and DP models yield similar results as the NEP-TNEP approach in the region up to about 1900 cm −1 .On the other hand, all ML models underestimate the intensity of the association band between 1900 and 2500 cm −1 , which is arising from the combination of librational and bending modes [37,75].Here, the NEP-TNEP prediction is actually still the one that comes closest to the experimental spectra.
The broad high-frequency peak above 3000 cm −1 , which is associated with the OH stretch mode, is notably blue-shifted and broadened for the T-EANN model, while the DP model strongly underestimates the intensity of this peak.In contrast, the NEP-TNEP combination predicts this feature in good agreement with the experimental data.
Finally, the parametric MB-pol model yields the worst agreement with experiment, for example, strongly overestimating the intensity of the bending band while underestimating the libration band.
With regard to the isotropic Raman spectrum (Fig. 4c), one should first note the variation among the experimental data.In particular in the region below 1000 cm −1 , the resulting uncertainty is comparable or even larger than the deviation between the NEP-TNEP prediction and the experimental data, while the position of the libration band predicted by T-EANN appears redshifted.With regard to the higher frequency region both NEP-TNEP and T-EANN reproduce the bending band well.In the case of NEP-TNEP this also applies for the OH stretch band, whereas in the case of T-EANN a blueshift can be observed similar to the anisotropic spectrum (Fig. 4b).

C. IR spectrum of PTAF -
The NEP-TNEP method for predicting IR spectra can be easily adopted for other molecular systems as long as the underlying observables to be learned are available.Naturally, this includes the molecular configurations along a chemical reaction, such that experimentally observable spectral changes can be connected to metastable complexes.One such complex is PTAF -(see inset in Fig. 5), the intermediate reaction minimum in the deprotection reaction 1-phenyl-2-trimethylsilylacetylene (PTA) with tetra-n-butylammonium fluoride [78][79][80][81].
To train NEP and TNEP models, we obtained PES and µ data for a set of 20 170 structures via DFT calculations using the ORCA code [82], the PBE functional [83], and a def2-TZVP basis set [84] while enforcing tight convergence of the self-consistent field cycles.Subsequently, MD simulations at various temperatures were performed in the NVE ensemble using a timestep of 0.1 fs for 1 ns, during which µ(t) was recorded with a time resolution of 0.5 fs.
The IR spectra obtained via the analysis of the ACF of µ show a pronounced temperature dependence in particular of the linewidths (Fig. 5).The molecule supports several soft modes with frequencies in the region below 250 cm −1 , which are associated with bending of and rotation about the ethynyl linker.These modes in particular FIG. 5. IR spectra for the metastable PTAF -complex (see inset) at various temperatures.The gray dashed line represents the broadened integrated absorption coefficients of the harmonic spectrum obtained directly from DFT calculations.The overall agreement is good considering the lack of anharmonic corrections (intensity and vibrational frequencies) and temperature sensitivity of the spectrum obtained from DFT calculations.
lead to strong mode coupling (i.e., anharmonicity), which underlies the changes in linewidth and the redistribution of the dipole strength across the spectrum.Here, the computational efficiency of the NEP-TNEP implementation in gpumd was crucial to resolve these features, as it enabled sampling on the nanosecond time scale, which would be prohibitive for a DFT-MD simulations and computationally very expensive for a CPU implementation.

D. Raman spectra of BaZrO3
BaZrO 3 is a perovskite that is being investigated, e.g., as a proton conductor for applications in fuel cells.It has also been the subject of various fundamental studies, as it is a prototypical antiferroelectric perovskite [86][87][88].It features soft and strongly temperature-dependent phonon modes [89,90], which have been carefully analyzed with Raman spectroscopy [85], rendering BaZrO 3 an ideal application for the present approach.
For benchmarking, we constructed models for χ using both the TNEP and SA-GPR approaches.The reference data set comprised cubic and tetragonal supercells with up to 40 atoms.The training structures were taken from MD simulations at different temperatures and pressures, generated using a NEP PES model constructed in an earlier study [90].In total the reference data set contained 940 structures.140 structures were randomly placed in a hold-out set for validation, while training sets were com-  piled by the shuffle-split method (random selection with replacement) with 200 to 800 structures and five data sets per training set size.A comparison of models generated using different choices for the size of the neural network as well as the descriptor demonstrates that viable models can be obtained for a wide range of parameters, and that even small models with as few as 1500 or so parameters can yield very good results (Fig. S12).Yet fine-tuning of these parameters as well as the regularization parameters (Fig. S13) allows one to maximize model performance.
The convergence of RMSEs and R 2 scores with training set size is similar for TNEP and SA-GPR with a slightly better performance for TNEP (Fig. 6).In both cases, training sets of about 400 structures already yield very good models, demonstrating the data efficiency of these approaches.This behavior has also been observed in the construction of models for amino acids [91].
Next MD simulations were carried out using 12×12×12 supercells (8640 atoms) and a timestep of 1 fs using the NEP model for the PES.Following equilibration at 300 K and 0 GPa in the NPT ensemble, the time-dependent sus-ceptibility χ(t) was recorded for 500 ps using a time resolution of 5 fs.For production, we used a TNEP model for χ trained against the full data set but we found that models based on at least approximately 400 structures to yield results that are practically indistinguishable within the statistical uncertainty.The Raman lineshape was subsequently obtained via the ACF of χ according to Eq. (11).We then computed the Raman spectra for parallel (Fig. 7a,b) and crossed polarization (Fig. 7c,d), which in Porto notation correspond to Z(XX) Z and Z(XY ) Z, respectively, where X and Y are arbitrary crystal axes.The final spectra were obtained by averaging over 20 independent MD trajectories.
The results are overall in very good agreement with experiment, especially considering the very strong anharmonicity of this material and the strong temperature dependence of the vibrational spectrum [89].The main difference with respect to the position of the peaks is a slight red-shift in the predicted spectra in the region above 600 cm −1 .This overly soft response can be attributed to the underlying exchange-correlation functional (vdW-DF-cx, Refs.92,93), which the NEP model truthfully reproduces.One can also observe an inversion in the intensity of the low and high energy features.This effect is almost certainly due to the classical sampling used here.It is rather common to correct for quantum effects in IR and first order Raman spectra by including a factor similar to the prefactor in Eq. ( 9).In the case of BaZrO 3 the room-temperature Raman spectrum arises, however, due to second-order scattering, i.e., due to combinations of modes.In that case, the application of the commonly used correction factor is no longer valid.Here, we therefore omit such corrections entirely.
The Raman spectra depend on the crystal orientation with respect to the excitation laser.The present approach allows one to readily map out this dependence via Eqs.(10) and (11) (Fig. 7b,d).While we are unaware of experimental measurements of the polarization dependence for BaZrO 3 , we note that such experiments have been carried out for, e.g., NaCl [94].As demonstrated in the former study, such measurements can provide valuable additional information.

V. CONCLUSIONS
In this contribution, we have introduced an extension of the NEP approach to tensors, resulting in the TNEP scheme.This was achieved by constructing expressions for rank-1 and rank-2 tensors based on the expression for the virial, which is a rank-2 tensor that arises naturally from derivatives of the energy (a rank-0 tensor) with respect to the atomic distances.This approach, which can be extended to tensors of higher rank, thus allows one to easily construct models that are equivariant.
We demonstrated the accuracy of this approach and its computational efficiency by constructing models for the dipole moment µ, the molecular polarizability α, and the electric susceptibility χ for several molecules, a liquid as well as two crystalline materials.In particular, the computational speed of the current method and its implementation in the gpumd package provide a significant advantage both in terms of the time scales and system sizes that can be sampled.
Finally, we applied the approach to predict IR and Raman spectra of liquid water, the molecule PTAF -, and the perovskite BaZrO 3 in very good agreement with available experimental data, illustrating the range of systems that can be readily addressed using the TNEP methodology introduced here.

FIG. 1 .
FIG.1.Workflow for simulations of IR and Raman spectra using NEP models for the PES and TNEP models for the dipole moment µ, the polarizability α or the susceptibility χ.

FIG. 2 .
FIG.2.RRMSEs for the validation sets according to TNEP, T-EANN, and SA-GPR models for water systems for (a) µ as well as (b) α and χ/ρ.Validation RRMSEs for liquid water from T-EANN[39] was reported for the averaged molecular polarizability obtained via the Clausius-Mossotti relation (Sect.S8).The validation RRMSEs for χ/ρ should be somewhat higher than that for the averaged molecular polarizability (also see TableS5).

FIG. 3 .
FIG.3. of computational speed of SA-GPR and TNEP models for dipole (µ) and susceptibility (χ) of liquid water.Here, the SA-GPR results were obtained using the TENSOAP-FAST implementation[54].
FIG. 6. Variation of (a) RMSEs and (b) R 2 scores with training set size for TNEP and SA-GPR models based on five training sets per size generated by shuffle-split.In the case of TNEP, we used Nneu = 20, n R max = n A max = 4, and λ1 = λ2 = 2 × 10 −3 (compare Figs. S12 and S13).

FIG. 7 .
FIG.7.Raman spectra of BaZrO3 for (a,b) parallel and (c,d) crossed polarization from simulations using a combination of NEP and TNEP models (red lines) as well as experiment (gray lines)[85].The spectra shown in (a,c) have been predicted for the nominal alignments used in the experimental measurements.The corresponding polarizations are indicated by the dashed horizontal lines in (b,d).

TABLE I .
RMSEs (in e • bohr) and RRMSEs (unitless) for µ for the validation sets using NEP rank-1 tensor models.For liquid water, the dipole moment is given per water molecule.

TABLE II .
RMSEs (in bohr3) and RRMSEs (unitless) for α (molecules) and χ/ρ (liquid water) for the validation sets using TNEP rank-2 tensor models.For liquid water, the χ/ρ is given per water molecule.Fe 2 O 3 (Sect.S5).The TNEP model for α-Fe 2 O 3 yields a very high R 2 score for the validation set close to one.