Using Gaussian Process Regression to Simulate the Vibrational Raman Spectra of Molecular Crystals

TL;DR

This paper demonstrates that Gaussian process regression can efficiently predict vibrational Raman spectra of molecular crystals, reducing computational costs while maintaining high accuracy, and enabling transferability across different polymorphs.

Contribution

The study introduces a symmetry-adapted GPR framework for predicting vibrational spectra, showing its effectiveness and transferability in molecular crystal simulations.

Findings

Accurate spectra with less than 1000 training points.

Symmetry-adapted GPR outperforms other GPR variants.

Model trained on one polymorph predicts spectra of another.

Abstract

Vibrational properties of molecular crystals are constantly used as structural fingerprints, in order to identify both the chemical nature and the structural arrangement of molecules. The simulation of these properties is typically very costly, especially when dealing with response properties of materials to e.g. electric fields, which require a good description of the perturbed electronic density. In this work, we use Gaussian process regression (GPR) to predict the static polarizability and dielectric susceptibility of molecules and molecular crystals. We combine this framework with ab initio molecular dynamics to predict their anharmonic vibrational Raman spectra. We stress the importance of data representation, symmetry, and locality, by comparing the performance of different flavors of GPR. In particular, we show the advantages of using a recently developed symmetry-adapted version of GPR. As an examplary application, we choose Paracetamol as an isolated molecule and in different crystal forms. We obtain accurate vibrational Raman spectra in all cases with less than 1000 training points, and obtain improvements when using a GPR trained on the molecular monomer as a baseline for the crystal GPR models. Finally, we show that our methodology is transferable across polymorphic forms: we can train the model on data for one structure, and still be able to accurately predict the spectrum for a second polymorph. This procedure provides an independent route to access electronic structure properties when performing force-evaluations on empirical force-fields or machine-learned potential energy surfaces.

Equations

$$I(\omega )$$

$$I_{\text{iso}}(\omega )$$

$$I_{\text{aniso}}(\omega )$$

$$\mathbf{\alpha }$$

$$\bar{\alpha }$$

$${\alpha }_{\gamma \delta }^{\text{ML}}(\mathcal{A})={\bar{\alpha }}_{\gamma \delta }^{\text{ai}}+\sum _{j=1}^N{w}_j^{\gamma \delta }k(\mathcal{A},{\mathcal{A}}_j),$$

$$k(\mathcal{A},{\mathcal{A}}_j)=e^{-\frac{|\mathbf{u}(\mathcal{A})-\mathbf{u}({\mathcal{A}}_j){|}^2}{2{\sigma }^2}},$$

$${\mathbf{w}}^{\gamma \delta }={\left(\mathbb{K}+\eta 1\right)}^{-1}\cdot \Delta {\mathbf{\alpha }}_{\gamma \delta },$$

$${\rho }_s(\mathbf{r})=\sum _{i\in s}\exp \left(-\frac{|\mathbf{r}-{\mathbf{r}}_i{|}^2}{2{\gamma }_s^2}\right),$$

$${\mathbf{k}}^{(2)}(\mathcal{A},\mathcal{B})=\int \mathrm{d}\hat{R}{\mathbf{D}}^{(2)}(\hat{R})\kappa (\mathcal{A},\hat{R}\mathcal{B}),$$

$${\mathbf{\alpha }}^{(2)}(\mathcal{A})$$

$${ϵ}_{\gamma \delta }[\sigma ,\eta ]=100\times \sqrt{\frac{\frac{1}{N}{\sum }_j^N{\left({\alpha }_{\gamma \delta }^{\text{ML}}({\mathcal{A}}_j)-{\alpha }_{\gamma \delta }^{\text{ai}}({\mathcal{A}}_j)\right)}^2}{\frac{1}{N-1}{\sum }_j^N({\alpha }_{\gamma \delta }^{\text{ai}}({\mathcal{A}}_j)-{\bar{\alpha }}_{\gamma \delta }^{\text{ai}}{)}^2}},$$

$${\nu }^2=\frac{1}{N_{\text{val}}}\sum _{j=1}^{N_{\text{val}}}\frac{||{\mathbf{\alpha }}_{\text{pred}}(j)-{\mathbf{\alpha }}_{\text{ref}}(j)|{|}^2}{{\sigma }^2(j)}$$

$${\mathbf{\alpha }}_k^{\prime }=\bar{\mathbf{\alpha }}+\nu ({\mathbf{\alpha }}_k-\bar{\mathbf{\alpha }})$$

$${\alpha }_{\gamma \delta }^{\text{ML}}(\mathcal{A})={\bar{\alpha }}_{\gamma \delta }^{\text{ai,crys}}+{\alpha }_{\gamma \delta }^{\Sigma }(\mathcal{A})-{\bar{\alpha }}_{\gamma \delta }^{\Sigma }+\sum _{j=1}^N{w}_{j}k(\mathcal{A},{\mathcal{A}}_j),$$

$$G_{im}^{\text{ref}}=R_{im}{G}_{im},$$

$${\mathbf{\alpha }}_{i,m}^{\text{mol}}=R_{im}^{-1}{\mathbf{\alpha }}_{i,m}^{\text{ref,mol}}{R}_{im}.$$

$${\mathbf{\alpha }}^{\Sigma }=\sum _{i=1}^{N_{\text{mol}}}{\mathbf{\alpha }}_i^{\text{mol}},$$

$$\Delta\alpha_{\gamma\delta}{\color[rgb]{0,0,0}{{(\mathcal{A})}}}=\alpha^{\textrm{\scriptsize crys}}_{\gamma\delta}{\color[rgb]{0,0,0}{{(\mathcal{A})}}}-\bar{\alpha}^{\textrm{\scriptsize crys}}_{\gamma\delta}-(\alpha^{\Sigma}_{\gamma\delta}{\color[rgb]{0,0,0}{{(\mathcal{A})}}}-\bar{\alpha}^{\Sigma}_{\gamma\delta})\ ,$$

$${\alpha }_{\gamma \delta }^{\text{ML}}(\mathcal{A})={\bar{\alpha }}_{\gamma \delta }^{\text{ai,crys}}+{\alpha }_{\gamma \delta }^{\Sigma }(\mathcal{A})-{\bar{\alpha }}_{\gamma \delta }^{\Sigma }+\sum _{j=1}^N{w}_{j}k(\mathcal{A},{\mathcal{A}}_j).$$

Full text

††: New J. Phys.

\ioptwocol

1 Introduction

Machine-learning (ML) models are becoming increasingly popular in the field of atomistic simulations, providing a way to obtain data-driven physical insights [1, 2, 3] and reduce the cost of simulations [4, 5]. Most efforts have been concentrated into predicting total energies and forces from atomic coordinates [6, 7, 8, 9, 5, 10, 11, 12, 13], which are most often the largest cost in a first-principles simulation. More recently, machine-learning models have been also applied to the prediction of response properties of molecules [14, 15, 16, 17, 18]. When dealing with the response of a material to an applied field, the cost of a first-principles calculation is often larger than that of force evaluation. This is thus an area where one can also take advantage of supervised learning techniques in order to reduce the cost ab initio simulations that make use of such response properties.

Vibrational Raman spectra are a good example of properties that requires the knowledge of the response of the system to electric field perturbations. The Raman signal is very useful to monitor phase transitions, as well as for the identification of global and local structural patterns [19, 20, 21]. Any technique used to calculate this property requires the calculation of several instances of the polarizability tensor (in molecules) or the dielectric susceptibility (in crystals). Previously, some of the present authors have shown that anharmonic vibrational Raman spectra calculated through a time-correlation formalism can be a powerful tool to identify structural fingerprints in molecular crystals [22, 23]. Within this formalism, it is necessary to calculate ab initio molecular dynamics trajectories and compute the response quantities for subsequent atomic configuration, employing, for instance, density-functional perturbation theory (DFPT) [24, 25, 26, 27, 22]. These calculations are computationally demanding, not only because of the tens of thousands of force evaluations that need to be performed to provide sufficient statistical sampling, but also because each DFPT calculation is typically four times more expensive than a force evaluation [23]222This number is clearly system- and settings-dependent, and simply represents the estimate reported in Ref.[23].. Furthermore, while there are several empirical potentials available that can be used to simulate the dynamics of molecular crystals [28], empirical models of the polarizability tensors are rare and poorly transferable [29].

Here we investigate frameworks to obtain accurate predictions of the dielectric response properties for a number of consecutive molecular dynamics (MD) configurations, that are necessary to converge simulated vibrational Raman spectra of molecules and molecular crystals. We compare different flavors of Gaussian Process Regression (GPR), which is a method that has already been proven to be efficient in predicting dielectric response properties [15, 17, 14]. In particular, we compare standard GPR schemes with symmetry-adapted GPR (SA-GPR)[30], which is advantageous when describing tensorial quantities. For the former, we exploit the internal structural rigidity of the system in order to model each individual component of the polarizability tensor, therefore remapping the tensor learning problem onto many separate scalar regression tasks. For the latter, we employ a symmetry-adapted representation of the system in order to learn the irreducible spherical components of the tensor in a covariant fashion. As a trial system we consider the Paracetamol molecule and form I and form II of the Paracetamol crystal, represented in Fig. 1. As we demonstrate below, SA-GPR comes out as the methodology with the best performance. When predicting the Raman spectra of crystals, both methods can benefit from using a GPR trained on the monomer as a baseline. We also find that our model is transferable between different polymorphic forms. Given that empirical and machine-learned potential energy surfaces are becoming more accurate for molecular crystals, the methodology proposed here can be combined in a straightforward manner to such potentials, giving access to the electronic polarization and polarizability of crystals.

In the following, Section 2 introduces the general machinery describing GPR and SA-GPR, as well as the different representations used to map structures to the inputs of the ML scheme. In Section 3, we present applications on the Paracetamol molecule and monoclinic form I crystal, showing for the latter how our models can be refined by the inclusion of molecular polarizability tensors. Finally, we illustrate in Section 4 the transferability of the SA-GPR method, by predicting the Raman spectrum of the orthorhombic form II crystal without prior knowledge of the corresponding values of the dielectric response.

2 Theory

2.1 Vibrational Raman Spectra

The central quantity needed for simulating vibrational Raman spectra is the static polarizability tensor $\bm{\alpha}$. For simplicity, we will refer to the polarizability tensor for all the rest of the paper, but one should keep in mind that for solids the quantity of interest is rather the electric susceptibility tensor of the system.

As discussed in Ref. [23] and others [31, 32], the vibrational Raman spectrum can be calculated using several approximations, the simplest of which is the harmonic approximation. The framework developed in this paper can be applied to the harmonic case, as shown in the supplemental information, Fig. S8, but we will here focus on the more challenging task of applying it to the linear-response time-correlation formalism, which fully takes into account the anharmonicity of the potential energy surface. In this formalism, vibrational Raman intensities can be obtained from the Fourier transform of the static polarizability autocorrelation function [33] at thermodynamic equilibrium. In particular, the so-called powder spectrum intensity is given by a combination of the anisotropic and isotropic contributions as

[TABLE]

where $n$ is the number of atoms in the system, the brackets $\langle\cdot\rangle$ denote an ensemble average and Tr represents the trace. $\bar{\alpha}$ and $\bm{\tilde{\alpha}}$ are, respectively, the isotropic and anisotropic parts of the polarizability tensor, defined as follows,

[TABLE]

In this paper, we do not address the problem of obtaining forces, which are necessary to calculate a full Raman spectrum from either molecular dynamics trajectories or in the harmonic approximation. Instead, we focus on predicting only the Raman intensities and combine them with pre-computed ab initio trajectories and displacements.

2.2 Component-wise Gaussian Process Regression

Gaussian process regression (GPR) is a well-established method based on a kernel function that measures the similarity between structures. It is formally equivalent to kernel-ridge regression (KRR), to which it differs mostly by the fact that GPR frames the construction of the model in a Bayesian language. In this sense, the kernel represents the prior distribution for the statistical correlations of the property we aim to predict. In the usual supervised-learning framework, a dataset of structures (i.e., atomic coordinates) and associated polarizabilities is used to train the model. Once the training is complete, the model is tested on a different set of configurations for which the polarizability is in principle not known. If it is possible to align the system to a reference structure, a straightforward way to apply this procedure to tensorial quantities is to learn each component of the polarizability tensor separately. In particular, a GPR prediction of each individual polarizability component $\alpha_{\gamma\delta}$ reads

[TABLE]

where $N$ is the number of configurations included in the training set, $\gamma$ and $\delta$ represent Cartesian coordinates, $w_{j}^{\gamma\delta}$ are the regression weights that need to be determined from the training data for each component, and $k$ is the kernel that couples the target system $\mathcal{A}$, with the training structure $\mathcal{A}_{j}$. The quantity $\bar{\alpha}^{\textrm{\scriptsize ai}}_{\gamma\delta}$ is the average (over the training set) of the $\gamma\delta$ polarizability component computed from an ab initio method, which effectively allows the training to focus on the fluctuation of the property with respect to a known baseline value.

The kernel entering GPR is based on a Gaussian similarity measure between structures $\mathcal{A}$, given by

[TABLE]

with $\sigma$ being a hyperparameter that controls the magnitude of the correlation between training points. We optimize the value of $\sigma$ - as well as that of other hyperparameters entering the definition of the features or the kernel - by cross-validation (CV), even though in a GPR framework hyperparameters are often determined by likelihood maximization (see SI, section 1 for further details). $\bm{u}$ is a vector that has the role of mapping the atomic coordinates of the structure $\mathcal{A}$ to a given representation of dimension $M$.

The regression weights $\mathbf{w}^{\gamma\delta}$ are obtained by minimizing a loss function regularized by an $L2$-norm [34] over the training set. This procedure leads to the following expression

[TABLE]

where $\mathbb{1}$ is the identity matrix, $\mathbb{K}$ is the $N\times N$ kernel matrix associated with the reference structures, such that $\mathbb{K}_{ij}=k(\mathcal{A}_{i},\mathcal{A}_{j}),i,j=1,...,N$, and $\eta$ is the regularization parameter which controls to which extent the fitted data can deviate from the training points. The quantity $\Delta\bm{\alpha}_{\gamma\delta}$ represents the vector containing all $N$ entries in the training set of \Delta\alpha_{\gamma\delta}^{j}{\color[rgb]{0,0,0}{{(\mathcal{A}_{j})}}}=\alpha^{\textrm{\scriptsize ai}}_{\gamma\delta}(\mathcal{A}_{j})-\bar{\alpha}^{\textrm{\scriptsize ai}}_{\gamma\delta}.

The efficiency of any GPR model strongly depends on the quality of the representation that is encoded in $\bm{u}(\mathcal{A})$ (see Eq. 3). For a representation to be efficient, it should contain the least possible number of elements to express the identity of a given structure, while avoiding redundant information. Among the vast choice of representations one could conceive, we will employ two that are very similar in spirit. First, we will consider a representation of each structure in terms of its atomic density (AD), explicitly evaluated on a grid. This AD is combined with the GPR framework, and with an alignment procedure that makes it possible to learn tensor components independently. Second, we will consider a symmetry-adapted Gaussian process regression (SA-GPR) scheme [30], whose application is summarized in the next section. The framework is based on $\lambda$-SOAP kernels, that are built by covariant integration of the atom density [35], and therefore automatically incorporate molecular symmetries.

Atomic-density grids consist of a conceptually simple representation. The construction introduces a Cartesian reference frame centered on the system under study, and defines a 3-dimensional grid around it. For each grid-point $\mathbf{r}$, we calculate the atomic density distribution, defined as

[TABLE]

where $s$ identifies a specific atom type, $\mathbf{r}_{i}$ are the nuclear positions and $\gamma_{s}$ is an adjustable parameter. The feature vector is given by $\bm{u}(\mathcal{A})=(\rho_{s}(\mathbf{r})),s=1,...,N_{s}$, where $N_{s}$ is the number of different atomic species. For the applications using the standard GPR scheme shown in this paper, we have chosen the same $\gamma_{s}=0.5$ Å for all species. This choice was based on the fact that we observed very minute changes when employing different values of $\gamma_{s}$ for different species, but this observation is likely to be system- and method- dependent. We have used a grid of evenly-spaced points spanning the maximal extension of the system, although we note that more refined methods could be utilized to define physically-motivated grid points, based on the possible directions of vibrations of atomic species [36]. Such a representation is illustrated in Fig. 2, which shows a 2-dimensional cut of the density distribution associated with a paracetamol molecule. Note that one could use a combination of different descriptors instead of a single one, whose elements would be concatenated to form the feature vector $\bm{u}$.

As mentioned above, the components of the polarizability tensor are not invariant to rotations in cartesian space, and therefore a Cartesian space representation of this quantity, like the atomic densities, requires an alignment to a reference structure. To do so, we have used the Kabsch algorithm considering only heavy atoms [37, 38]. We note that this alignment procedure is not applicable to all systems, in particular to very flexible molecules. For simple relatively rigid molecules, it is known to work well, as has been shown previously for the case of hyperpolarizabilities of water molecules [15]. We will show below that even for a more complex and relatively flexible molecule like paracetamol, this representation still yields accurate predictions.

2.3 Symmetry-adapted Gaussian Process Regression with $\lambda$-SOAP Kernels

SA-GPR [30] represents a generalization of the GPR formalism, where the tensorial nature of the target property, together with its covariant 3D transformations, are naturally incorporated within the regression algorithm [12]. This technique removes the need of often arbitrary alignment procedures of the systems (like the one described above) in order to predict tensorial quantities of any rank. As such, the model focuses only on the portion of variability of the tensor connected with an internal structural distortion of the molecular geometry, greatly improving the regression performances. Furthermore, by effectively learning an atom-centered model of the polarizability, a SA-GPR scheme makes it possible to build models that are transferable between different molecules [17], although we didn’t exploit this possibility in the present work.

A simplification of the learning problem can be obtained if the target property is first decomposed in its irreducible spherical tensor components. The static polarizability (or the static susceptibility), in particular, being a symmetric rank-2 tensor, can be formally decomposed into an isotropic contribution that transforms as a spherical harmonic of angular momentum $\lambda=0$, and an anisotropic contribution that transforms as a spherical harmonic of angular momentum $\lambda=2$. The former contribution, $\alpha^{(0)}_{0}\propto\bar{\alpha}$, being directly proportional to the trace of the tensor (which is rotationally and translationally invariant), can be learned in the usual manner by a standard GPR. The tensorial nature of the polarizability is contained in the $\lambda=2$ term, $\bm{\alpha}^{(2)}=(\alpha^{(2)}_{-2},\alpha^{(2)}_{-1},\alpha^{(2)}_{0},\alpha^{(2)}_{1},\alpha^{(2)}_{2})$, which is related to the anisotropic part of the Cartesian $\bm{\alpha}$-tensor $\bm{\tilde{\alpha}}$ of Eq. 2 by a linear transformation. Within SA-GPR, the prediction of this contribution is carried out by making use of a tensorial kernel function $\bm{k}^{(2)}(\mathcal{A},\mathcal{B})$, which is a matrix of size 5$\times$5, that describes at the same time the structural similarity between molecular configurations and the tensorial geometric relationship of order $\lambda=2$ between these configurations, as detailed in Ref. [30]. In general terms, such a kernel can be thought of as a generalization of a scalar kernel function such as the one introduced in Eq. 4, which makes use of the following covariant integration,

[TABLE]

where $\hat{R}$ represents the rotation operator and $\bm{D}^{(2)}(\hat{R})$ the associated Wigner-D matrix which has the role of expressing the rotation of $\lambda=2$ spherical harmonics. In this definition, the scalar kernel $\kappa(\mathcal{A},\mathcal{B})$ only needs to be invariant under rigid translations and rotations of the laboratory reference frame with respect to which both $\mathcal{A}$ and $\mathcal{B}$ are defined. In the present work, $\kappa(\mathcal{A},\mathcal{B})$ is given by a superposition of atom-centered Gaussian densities, equivalent to the one used in the popular Smooth Overlap of Atomic Position (SOAP) kernel [39]. As such, the kernel of Eq. (7) represents the tensorial generalization of SOAP, usually called $\lambda$-SOAP, which recovers the scalar case of $\lambda=0$ as a special limit [16].

Formally, a covariant $\lambda=2$ prediction performed with a symmetry-adapted kernel function of the same order reads

[TABLE]

where $i$ and $j$ run over the $N$ reference configurations used to train the regression model and $\eta$ is the regularization parameter. The set of tensorial regression weights $\{\bm{w}^{(2)}_{i}\}$ are determined by inverting the $(5\times N)^{2}$ kernel matrix $\mathbb{K}^{(2)}$ associated with the training structures, and projecting it on the vector of reference tensors $\{\bm{\alpha}^{(2)}_{j}\}$. In doing so, the statistical average of $\{\bm{\alpha}^{(2)}_{j}\}$ over the training set is assumed to be zero by symmetry. As such, no baseline of the anisotropic part of the polarizability is adopted when doing the regression.

2.4 Errors and uncertainty estimations

In order to gauge the accuracy of the machine-learned polarizability components, we calculate the root mean square error (RMSE) normalized by the standard deviation (STD) of the set we want to evaluate the error on,

[TABLE]

where $j$ represents the $j$th configuration, $N$ is the amount of points in the data set we consider, and the bar here denotes the average of the polarizability component over the dataset of interest.

Whenever the reference properties of the testing data are not available, one cannot make use of Eq. (9) to evaluate the error of the predicted polarizabilities. In these circumstances, one typically needs to estimate the error made on the predicted properties by making use of some a priori probabilistic criteria. In the particular case of GPR, the expected error associated with a testing structure $\mathcal{A}$, can be computed as $\varepsilon(\mathcal{A})=k(\mathcal{A},\mathcal{A})-\sum_{IJ}k(\mathcal{A},I)k^{-1}(I,J)k(J,\mathcal{A})$, where $I$ and $J$ run over the training structures. As detailed in Ref. [40], this strategy is however not very practical because of its computational expense, so that other kind of methods such as bootstrapping or subsampling can rather be used to estimate the prediction errors. In addition, in the particular case of the present work one would like to propagate the error that occurs in the prediction of $\bm{\alpha}$ to the Raman spectrum. While this error propagation would be difficult to carry out on top of the GPR intrinsic covariance, it comes naturally from a bootstrapping or subsampling scheme.

In this work, $N_{\textrm{\scriptsize RS}}$ subselections of the training dataset are considered to generate an ensemble of predictions for the polarizability. From these, $N_{\textrm{\scriptsize RS}}$ Raman spectra are computed by Fourier transforming the time series of each model in the ensemble. Finally, the average and the standard deviation of the predicted spectra over the $m$ subselections give the final Raman spectrum prediction and the propagated estimated error respectively. The downside of this approach is that this model works under the assumption that the training data corresponds to independent identically distributed samples. This is of course not true in general, so that one needs to correct the model to take into account for the underlying correlations. Following Ref. [40], a maximum likelihood recipe can be adopted to linearly scale the variance of the predictions by a constant factor $\nu^{2}$. The calibration of this scaling factor is carried out by computing the actual prediction errors of the polarizabilites over a suitably selected validation set $N_{\textrm{\scriptsize val}}$, for which the reference polarizabilities are known, and then considering

[TABLE]

where $\sigma^{2}(j)$ are the variances of the predicted polarizabilities. Once the value of $\nu$ has been determined, each polarizability prediction of a given training model $k$ can be updated by considering

[TABLE]

where $\bar{\bm{\alpha}}$ is the predicted polarizability averaged over the $N_{\textrm{\scriptsize RS}}$ models. This scaling procedure guarantees that the variance of the models is consistent with the outcome of the likelihood maximization. By computing the Raman spectrum for each scaled model $k$, the propagated uncertainty estimation associated with the spectra will automatically take into account the calibration of the variance.

2.5 Simulation and training details

To perform the ab initio calculations, we used the FHI-aims [41] program package with light settings for all atomic species. We obtained aiMD trajectories using the PBE functional with many-body dispersion (MBD) corrections [42, 43], and employing a time step of 0.5 fs. The polarizabily/susceptibility was instead computed every 1 fs. Most of the data used here was already available from references [22] and [23]. For each system, we had 20 picoseconds of simulation in the NVT ensemble at 300 K. From this trajectory, a few thousands of configurations were selected to define the training and test set of our model. A full trajectory of 15 picoseconds in the NVE ensemble was instead considered as our validation set to assess the quality of our predictions, by comparing the predicted Raman spectra to the ones obtained with ab initio polarizabilities. We chose to train our model on the NVT ensemble and predict on the NVE ensemble. This choice is justified by the fact that - particularly for this relatively short trajectory - canonical sampling should be more ergodic than microcanonical sampling, and therefore sample a larger portion of configurational space.

3 Results

In the following, we apply the previously described methods to the calculation of the Raman spectra of Paracetamol. Each time we mention GPR, it is used in combination with the AD representation. Similarly, each time SA-GPR is mentioned, it is associated to $\lambda$-SOAP kernels.

3.1 Paracetamol molecule

We first consider the case of the isolated Paracetamol monomer (Fig. 1 (a)). We constructed the training dataset by selecting 2000 structures with farthest point sampling (FPS) using the scalar SOAP metric from the full NVT trajectory. For GPR, the three-dimensional density field was constructed within a box of $6\times 4\times 2.5$ Å, where the molecule had its longest axis along the $x$ direction and the equilibrium geometry lied on the $xy$ plane, and the grid spacing was $dr=0.5$ Å. The GPR was computed using $\sigma=10$ and $\eta=10^{-3}$. Details about the optimization of the hyperparameters are given in the supplemental information. Regression performances are reported in Fig. 3 (a), where the error $\epsilon$, given by Eq. 9, computed on 500 randomly selected test structures (that are excluded from training) out of the total of 2000 is shown as a function of the number of training monomers. 333Learning performance when using only unprocessed atomic coordinates are also shown in the SI, Fig. S2.

Figure 3(a) shows that the learning capability of the model does not reach saturation within the training set sizes explored. The learning of all the polarizability components follow a similar slope, but they are predicted with different accuracy because of the strong anisotropy of the dielectric response of the Paracetamol molecule. Because of the $\pi$-conjugation of the system in the molecular plane, the system is much more polarizable along the $x$-axis rather than along other directions in this particular alignment, making it harder for the learning algorithm to capture the corresponding variations across the dataset. The $\alpha_{xx}$ component presents the largest error, going from about 40% with 300 training points to 17% with 1500 training points. The best learning performance is instead obtained for the $\alpha_{yz}$ component, where the prediction error can be brought down to about 6%.

In order to validate the model, we show in Fig. 4 the correlation between predicted and computed polarizability components on the validation set composed by the full independent NVE trajectory, for a representative training size of 900 molecular configurations, including the $\epsilon$ for each component. Although we observe a worsening of the predictions when compared to the previous case, where we trained and predicted on the same ensemble (and trajectory), the predicted polarizabilities are still well correlated to the reference values. The remaining question is how these errors translate to the actual prediction of the vibrational Raman spectrum.

In Fig 5, we show a machine-learned Raman spectrum averaged over 16 subselections of the training set of 900 configurations each, along with its standard deviation (see Sec. 2.4 for a detailed discussion about this procedure), and compare it to the one calculated from fully ab initio data. We find that the estimated variance has to be scaled by a factor of $\nu^{2}=2.0$. Despite a relatively small amount of training points, the agreement with the reference spectrum is excellent in the entire frequency range. As shown in Fig. 3, increasing the number of training points in the model would decrease errors even further. From Figs. 3 and 5, it is clear that the error we make on the polarizability components does not translate directly into an error of similar magnitude on the spectrum. This is a consequence of the fact that the Raman intensities depend on the derivatives of the polarizability components with respect to atomic coordinates, and not on their absolute value. This simple procedure is able to reproduce almost perfectly a reference Raman spectrum with fewer than 1000 training points on a desktop computer in just a few minutes.

In order to assess whether the use of a model that incorporates symmetries can be advantageous even for a relatively rigid molecular system, we then contrast the GPR and the SA-GPR model. To this end, $\lambda$-SOAP kernels were constructed using a Gaussian width of 0.3 Å and an environment cutoff of 4.0 Å. Details about the SOAP parameters optimization can be found in the SI. The corresponding learning curve is shown in Fig. 3(b). The improvement over a standard GPR scheme is systematic at any training set size, underlining the importance of automatically incorporating the $\mathcal{O}$(3)-covariance of the tensor at the scale of individual atomic environments. Similarly to the case of GPR, the accuracy of predictions differs between tensor components. With 300 training points, the $\alpha_{xx}$ component presents the largest error (about 14%), that reduces to only 6% with 1500 training structures. The best learning performance is again obtained for the $\alpha_{yz}$ component, for which the prediction error remarkably can be reduced to less than 2%. As shown in Fig. 5, SA-GPR reproduces very accurately also the ab initio Raman spectrum with only 300 training points. Obviously, as reflected by the learning curves, increasing the amount of points reduces this error even further, as exemplified in Fig. S4 of the SI, where the machine-learned and ab initio spectra are virtually indistinguishable.

3.2 Paracetamol crystal

We now turn our attention to the first crystalline form of Paracetamol, containing four individual Paracetamol molecules per unit cell, as shown in Fig. 1 (b).

3.2.1 Direct approach

Since we are now dealing with a periodic system, we first build a supercell corresponding to the appropriate crystal structure. Then, the AD representation is constructed following the same procedure discussed in the previous section. A value of $\sigma=40$ has been selected to build the Gaussian kernel, and the regularization parameter has been set to $\eta=10^{-4}$ by CV optimization. The three-dimensional density field was evaluated within a box of $12\times 14\times 20$ Å3, using a grid spacing of $d\mathbf{r}=0.75$ Å444No noticeable improvement was observed when using a finer grid with $d\mathbf{r}=0.5$ Å.. For SA-GPR, the $\lambda$-SOAP kernels were constructed using the same parameters as before.

The training set is built by considering a random selection of 2500 configurations extracted from a NVT trajectory. A full NVE trajectory is once again considered to test the quality of the predicted Raman spectrum. Learning curves for both regression models are shown in Fig. 6 (solid lines). When comparing the two methods, we always use the same configurations in both cases.

We observe that both the learning capability of GPR and SA-GPR do not reach saturation when increasing the number of training data, going from 81% (respectively 73% with SA-GPR) of error with 25 training points to 17% (respectively 11%) with 2000 of them; Again, making use of a kernel that is built on a symmetry-adapted comparison between local environments brings a substantial improvement. Overall, however, for the same amount of training points, the errors are much (typically between two and three times) larger than for the monomer case.

3.2.2 Incorporating molecular polarizability

To improve our results, we refined our models by using the predictions of the non-interacting monomers included in the molecular crystal, as we explain in the following.

Suppose we have a molecular crystal made of $N_{\textrm{\scriptsize mol}}$ molecular units (in the case of Paracetamol I, $N_{\textrm{\scriptsize mol}}=4$, while $N_{\textrm{\scriptsize mol}}=8$ for Paracetamol II). Since we have already learnt the polarizability tensors of the individual molecules, we can grasp most of the variability of the polarizability tensor of the crystal by summing up the predictions for individual monomers. Equation 3 is modified according to

[TABLE]

where $\bm{\alpha}^{\Sigma}$ denotes the sum of the molecular polarizability tensors and $\bar{\alpha}_{\gamma\delta}^{\textrm{\scriptsize ai,crys}}$ is the average of the ab initio polarizability tensors of the full crystal over the training set. Specific details about the procedure are explained in A. An analogous expression to Eq. 12 is obtained for SA-GPR.

Figure  6 shows the advantage of using the molecular baseline for the regression models (dashed lines). The improvement is most noticeable for models based on few training points. Molecular baselining leads to a large decrease in error of about 25% for both GPR and SA-GPR. Upon increasing the number of training points, the difference diminishes, but an improvement remains visible. It is worth noticing that, from 250 training structures on, the direct application of SA-GPR (without any molecular baseline) performs better than the GPR scheme with the baseline. Note that the prediction accuracy can be improved even further if one scales the molecular polarizability tensors so that their average matches that of the full crystal, as illustrated in the SI, Fig. S8.

Figure 7 shows the effect of the baselining procedure on the predicted SA-GPR Raman spectrum, when one increases the amount of training points.

Several observations can be made. First, just like for the monomer, high frequencies require more training points to be reproduced, be it with a direct prediction or by including molecular polarizability tensors. Second, including $\bm{\alpha}^{\Sigma}$ greatly enhances the spectrum intensity accuracy when few training points are used. This is especially true at high frequencies, where the improved model already gives the right structure of the peaks, albeit not with the right intensity, while the direct learning does not show any peak in this region. Overall, the predicted spectrum is extremely well reproduced when employing enough training points.

Figure  8 shows the predicted Raman spectrum and the corresponding estimated error. In this case, different learning models have been first defined by considering 16 random subselections, each of them made of 80% of the training dataset. Then, for each of these learning models, the polarizabilities of the full NVE trajectory have been predicted and the associated Raman spectra have been computed. We then estimated the standard deviation of these predictions according to the procedure detailed in Sec. 2.4. We find that in this case the estimated variance has to be increased by roughly an order of magnitude, i.e., $\nu^{2}=10.9$. One can observe that the excellent agreement between the reference and predicted spectrum at low frequencies is consistent with a negligible estimated error, while larger discrepancies and error bars can be observed in the high-frequency domain.

Results this far show that, even for a relatively rigid molecular system, incorporating symmetries and learning all the components of the polarizability tensor in a covariant fashion can improve the accuracy and the efficiency of these ML schemes. As anticipated, an additional advantage of SA-GPR is that it is based on local atomic environments. This results in a greater transferability, as we will discuss in the following section.

4 Extrapolation on other polymorphic forms

Within the $\lambda$-SOAP formalism, the polarizability of the system is effectively decomposed in local atom-centered contributions that are summed in order to obtain the final predicted value of $\bm{\alpha}$, making it possible to model the susceptibility of a molecule or a crystal through the definition of effective atomic polarizabilities. This implies that the information is learned at the local scale and, as such, can be transferred across systems that share a similar chemical nature. In the case of paracetamol polymorphs, one can think of predicting the polarizability of the the form II crystal (Fig. 1 (c)) with the model trained on form I only. Since different polymorphic forms are mainly distinguished by the different intermolecular interactions, major difficulties in this extrapolation procedure are expected to be associated with the low-frequency (intermolecular) modes of the molecular crystal. To put this idea to the test, we used the SA-GPR trained on 2000 structures of form I to make predictions for a NVE trajectory of form II. We used an ensemble of 16 subsampled models to estimate uncertainty. As expected, we observed that the rather large error in the prediction of the polarizability tensor is mostly associated to a small offset in the time series of some of the polarizability components (detailed in the SI, Fig. S10), that however does not have a substantial impact on the Raman spectrum.

As shown in Fig. 9, the general lineshape is excellent, and all the main features of the ab initio spectrum are reproduced, even though few discrepancies in terms of intensity can be observed, and the error in the intensities is overall higher than for the direct prediction of the first polymorph. We underline the difference in behaviour in contrast to Fig. 8: now, high frequencies are better described and errors are more pronounced at low frequencies. This suggests that the model can reproduce accurately changes in the polarizability associated with intra-molecular vibrations, but is less accurate in predicting low-frequency components that are specific to the molecular packing of form II, which is not represented in the training set for form I. We also observe that the discrepancy between predictions and ab initio spectrum is reflected accurately in the estimated uncertainty, that can therefore be used as a reliable measure of the accuracy of the model also in an extrapolative regime.

5 Conclusions

In this work, we proposed GPR models to predict vibrational Raman spectra, based on learning polarizability and susceptibility tensors obtained from density-functional perturbation theory. As an example, we applied our methodology to predict anharmonic Raman spectra of the Paracetamol molecule and two polymorphs of the Paracetamol crystal. The methodology also works for simpler harmonic Raman intensities. The use of an ensemble of models to estimate the uncertainty in the polarizability tensors allows us to propagate the error estimation from the ML prediction of the polarizability tensors to the vibrational spectra, by generating an ensemble of spectra out of which it is simple to compute frequency-dependent confidence intervals.

We showed that for the molecule a standard GPR scheme that takes as input a nuclear density representation on a 3D grid works extremely well and enables one to reproduce Raman spectra almost perfectly with a low number of training points. For the crystal, such a scheme, albeit possible, is more difficult to apply for several reasons: The difficulty to compare crystal structures with different unit cell sizes, the redundancy of information contained in a fixed grid-based representation, and the increase of the number of grid points with system size. We have shown that a more effective solution is to use a symmetry-adapted GPR scheme, used here in combination with the $\lambda$-SOAP representation [30]. Such a scheme yields accurate predictions for molecules and crystals, due to its capability to better capture the local structural information in a covariant fashion. Moreover, since $\lambda$-SOAP is a local representation, it is easy to treat larger systems sizes and even transfer models to other polymorphs. We have shown this transferability by successfully predicting the Raman spectrum of Paracetamol II with a model trained only on Paracetamol form I. This suggests the possibility of predicting Raman spectra of any polymorphic form, as long as a model trained for one of them is available. In addition, for all models presented, we observe a considerable improvement when using previously-trained GPR models for the molecular units as a baseline for the crystal prediction, thus reducing the amount of more costly condensed-phase calculations that must be performed to train the bulk model. In a similar manner, it is also straightforward to extend this framework to other ensembles (e.g. NPT) or path-integral molecular dynamics simulations, which include the quantum nature of the nuclei.

The models we presented regard the electronic electric-field response properties, and can be extended to dipoles and higher-order responses. They can thus be seamlessly combined with empirical potentials or other machine-learned potentials that give access to forces. This presents an alternative route to including the training of such quantities directly into these potentials [44], which can present a higher level of complexity. Finally, we remark that even though we applied our framework to polarizabilities and Raman spectroscopy, applying it to any other kind of spectroscopy, like infrared, sum-frequency generation, etc., would be straightforward, as long as reference electronic-structure data is available.

6 Acknowledgments

The authors acknowledge funding from the Max Planck-EPFL center. MC acknowledges support from the European Research Council (Horizon 2020 Grant Agreement No. 677013-HBMAP).

Appendix A Molecular baselining with GPR

When building a model to predict the polarizability tensor of the crystal from its individual molecular components within the GPR and AD representation scheme, one should consider that each molecule is oriented in a specific direction that differs from the one we trained the molecule on. In order for the regression model to recognize the orientation at hand, we first need to find the relationship between the orientation of the molecules in the crystal and the reference one used in the molecular GPR procedure described in section 2.2. The scheme is depicted in Fig. 10.

Each geometry $G_{im}$ corresponding to the $i$th molecule in the crystal of the $m$th structure of the set is thus rotated by a rotation matrix $R_{im}$ as

[TABLE]

where $i=1\cdots N_{\textrm{\scriptsize mol}}$, with $N_{\textrm{\scriptsize mol}}$ the number of molecular blocks per unit cell. Finally, once a molecular polarizability \bm{\alpha}^{\textrm{\scriptsize ref{\color[rgb]{0,0,0}{{,mol}}}}}_{i} is predicted, we rotate the tensor back to its original orientation inside the crystal, i.e.,

[TABLE]

Having defined the sum of molecular polarizabilities $\bm{\alpha}^{\Sigma}$ as

[TABLE]

we consider the regression target

[TABLE]

where, once again, the bar denotes the average over the training set. Finally, we modify accordingly equation 3, which becomes (see also Eq. 12)

[TABLE]

Note that the weights in the previous equation are different than the ones in Eq. 3, as this time the regression target $\Delta\alpha$ contains the molecular polarizability tensors.

An analogous expression is obtained for SA-GPR, but the rotation and alignments previously described do not need to be carried out explicitly since the rotational covariance of the tensor is built in the structure of the method.