# Learning the Electrostatic Response of the Electron Density through a Symmetry-Adapted Vector Field Model

Mariana Rossi, Kevin Rossi, Alan M Lewis, Mathieu Salanne, Andrea Grisafi

► To cite this version:

Mariana Rossi, Kevin Rossi, Alan M Lewis, Mathieu Salanne, Andrea Grisafi. Learning the Electrostatic Response of the Electron Density through a Symmetry-Adapted Vector Field Model. Journal of Physical Chemistry Letters, 2025, pp.2326-2332. 10.1021/acs.jpclett.5c00165. [https://doi.org/10.1021/acs.jpclett.5c00165]

HAL Id: hal-04976897

https://hal.science/hal-04976897v1

Submitted on 4 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Learning the Electrostatic Response of the Electron Density through
a Symmetry-Adapted Vector Field Model

Mariana Rossi,¹ Kevin Rossi,²,³ Alan M. Lewis,⁴ Mathieu Salanne,⁵,⁶ and Andrea Grisafi⁶,∗

¹Max Planck Institute for the Structure and Dynamics of Matter,
Luruper Chaussee 149, 22761 Hamburg, Germany

$^2$Department of Materials Science and Engineering,
Delft University of Technology, 2628 CD, Delft, The Netherlands

$^3$Climate Safety and Security Centre, TU Delft The Hague Campus,
Delft University of Technology, 2594 AC, The Hague, The Netherlands

⁴Department of Chemistry, University of York, Heslington, York, YO10 5DD, UK

$^5$Institut Universitaire de France (IUF), F-75231 Paris, France

$^6$Physicochimie des Électrolytes et Nanosystèmes Interfaciaux,
Sorbonne Université, CNRS, F-75005 Paris, France

A current challenge in atomistic machine learning is that of efficiently predicting the response of the electron density under electric fields. We address this challenge with symmetry-adapted kernel functions that are specifically derived to account for the rotational symmetry of a three-dimensional vector field. We demonstrate the equivariance of the method on a set of rotated water molecules and show its high efficiency with respect to number of training configurations and features for liquid water and naphthalene crystals. We conclude showcasing applications for relaxed configurations of gold nanoparticles, reproducing the scaling law of the electronic polarizability with size, up to systems with more than 2000 atoms. By deriving a natural extension to equivariant learning models of the electron density, our method provides an accurate and inexpensive strategy to predict the electrostatic response of molecules and materials.

The electrostatic response of the electron density is a
fundamental property of matter. It determines vibrational
Raman and sum-frequency cross-sections [1, 2], and the
optical dielectric constant of materials [3, 4]. Furthermore,
it defines the induced electronic polarization of metallic
systems [5, 6], of central interest for electrochemical energy
storage [7] and electrostatic catalysis [8–10]. This elec-
trostatic property is defined as the static linear response
function of the electron density under the application of
a homogeneous electric field along the three Cartesian
directions; as such, it can be formally represented as a
continuous three-dimensional vector field.

The calculation of electronic density responses via
density-functional theory (DFT) requires the evaluation of
the first-order perturbative term in the density expansion
around the ground-state, in absence of any external fields.
This response can be obtained through suitable finite-
difference calculations. Density-functional perturbation
theory (DFPT), however, provides a more elegant route to
compute this quantity, allowing the treatment of both pe-
riodic and aperiodic systems on an equal footing [11–13].
These calculations come with an added simulation cost
involving a larger prefactor with respect to ground-state
calculations, and a scaling with the number of electrons
that is the same as that of DFT; this aspect vastly lim-
its the system sizes that can be afforded, as well as the
number of calculations that can be performed.

Bypassing the need for large computational resources is
one of the central goals of machine-learning (ML) method-
ologies in atomistic simulations. In this context, ML mod-
els that are made equivariant under three-dimensional

spatial symmetries have proved to be very efficient in learn-
ing electronic densities of molecules [14–17], solids [18, 19],
liquids [20, 21], and surfaces [22–25]. The equivariant na-
ture of these methods drastically decreases the amount
of data necessary to train reliable models, and typically
enhances the overall accuracy by making predictions more
robust to numerical noise. Nevertheless, building a fully
equivariant framework for the prediction of vector-field
properties represents a current challenge of atomistic ma-
chine learning. This problem has been tangentially faced
in the description of spin degrees of freedom, as those
entering the construction of ML interatomic potentials for
magnetic materials [26, 27]. In that case, the atomistic
representation must be able to capture the symmetries of
the vector-field, but the learning target, i.e., the electronic
energy, is invariant under rotations.

In this work, we address the challenge of achieving
an equivariant atom-centered representation of the vec-
tor field and preserving the rotational symmetry in the
prediction target, which is here the density response to
applied electric fields. From a ML model perspective,
an advantage of directly targeting the density response
rather than its derived properties (e.g., polarizabilities)
is that the response is inherently local, making it more
likely to yield robust predictions for highly polarizable
systems. Previous works by some of the authors have
proposed methods to address this problem in the context
of symmetry-adapted Gaussian process regression [23, 28]
(SA-GPR), by extending the SALTED method for predict-
ing electron densities [18, 21]. These methods relied on a
modification of the descriptor of the local atomic environ-
ment (features), in order to encode the information about the direction (and the strength) of the applied electric field. Recently, a similar approach has been adopted in the context of neural-network architectures [25], which also makes use of field-induced embedded atom density (FI-EAD) features [29]. In Refs. 23 and 25, these models were successfully applied to predict the density response of metallic slabs along the direction perpendicular to the surface, thus representing, effectively, the prediction of a scalar field. In Ref. 28, the prediction target was the full vector field of the density response of molecules and insulating solids or liquids. However, the prediction of each Cartesian component was obtained through independent SA-GPR models. Learning and prediction were proven to be successful, but the rotational equivariance of the model was not formally obeyed.

As we will see, an equivariant formulation of the learning problem capable of consistently taking into account the vector-field symmetries presents significant advantages over previous methods. We will show the superior accuracy of our approach in predicting the density response of benchmark datasets that include isolated water molecules, as well as supercells of liquid water and molecular crystals. Additionally, we will demonstrate applications to metallic nanoparticles of increasing size that display a nonlocal polarization, enabling the inexpensive calculation of polarizability tensors at large length scales.

In line with electron-density learning methods based on a linear atomic-orbital representation [14, 18, 20], the response function of the electron density $n_e$ under an applied electric field $E_k$ can be written as a linear combination over a set of atom-centered basis functions $\phi_{n\lambda\mu}$. Assuming a non-periodic system for simplicity, we can express the continuous vector field at any given point $\mathbf{r}$ as

$$ \frac{\partial n_e(\mathbf{r})}{\partial E_k} = \sum_{in\lambda\mu} c_{in\lambda\mu}^k \phi_{n\lambda\mu}(\mathbf{r} - \mathbf{R}_i), \quad (1) $$

where $k$, $n$ and $\lambda\mu$ are the Cartesian, radial and spherical harmonics indexes, respectively, and $\mathbf{R}_i$ are atomic positions. Reference quantum-chemistry calculations of this quantity can be performed using a DFPT approach, as well as by finite differences of the electron density under a sufficiently small electric field along the three Cartesian directions. For periodic charge densities, these calculations need to be performed taking into account the modern theory of polarization [30]. The set of expansion coefficients, $c_{in\lambda\mu}^k$, can then be extracted by using standard density-fitting techniques. Our goal is to derive a kernel-based approximation of the coefficients $c_{in\lambda\mu}^k$, by extending well-established SA-GPR methods for the prediction of $n_e$ [14, 18, 21].

We start by noticing that the vectorial character of Eq. (1) can be equivalently expressed in terms of $|1k\rangle$ angular momentum states. This property allows us to represent the rotational symmetry of the expansion coefficients

within the tensor product space $|\lambda\mu, 1k\rangle \equiv |\lambda\mu\rangle \otimes |1k\rangle$. A similar scenario is encountered in the context of molecular Hamiltonians [31], where the matrix elements between atomic orbitals transform under rotations as the tensor product between pairs of angular momenta. In that case, one can take advantage of an orthogonalization procedure to recast the learning target into a set of irreducible angular components that can be regressed independently of each other. Because of the non-orthogonal nature of the basis functions used for the expansion of Eq. (1), however, this strategy cannot be effectively adopted for the present problem [21]. The rotational symmetry of the coefficients $c_{in\lambda\mu}^k$ must be handled altogether in the learning model.

From the previous discussion, the problem of providing an equivariant approximation of the $c_{in\lambda\mu}^k$ can be addressed by deriving kernel functions that express the structural similarity between pairs of atomic environments $i$ and $j$, together with the geometric correlation between spherical tensors of order $\lambda \otimes 1$. Following the rationale of Refs. 32, 33, this can be achieved with SO(3) integrals that include tensor products of Wigner $D$-matrices associated with the rotational symmetries we aim to enforce. After some algebraic manipulations, detailed in the Supporting Information (SI), we find that kernel functions adapted to $\lambda \otimes 1$ symmetries can be written in terms of an irreducible angular momentum decomposition. Specifically, we obtain the central result:

$$ K_{\mu k, \mu' k'}^{\lambda \otimes 1} = \sum_{l=|\lambda-1|}^{\lambda+1} \langle \lambda\mu, 1k|lm\rangle \langle \lambda\mu', 1k'|lm' \rangle K_{mm'}^l, \quad (2) $$

where $\langle \lambda\mu, 1k|lm\rangle$ are the Clebsch-Gordan coefficients used for the composition of angular momenta, such that $l = \{|\lambda - 1|, \lambda, \lambda + 1\}$ and $m = \mu + k$. Crucially, the irreducible components $K^l$ in Eq. (2) can be recognized as spherical tensor kernels of order $l$, like those originally derived in the prediction of molecular tensors [33], and later used for the prediction of electron densities [14]. This makes the calculation of Eq. (2) straightforward by accessing the relevant triplets of $K^l$ for each value of $\lambda$ that enters the representation of Eq. (1). Equivariance in $O(3)$ can be finally achieved by enforcing $\lambda \otimes 1$ inversion symmetry, which implies considering symmetric $K^l$ for $l = |\lambda \pm 1|$ and antisymmetric $K^l$ for $l = \lambda$.

Following the latest formulation of the SALTED method [21], the $\lambda \otimes 1$ kernels so derived can be used to perform equivariant predictions of the coefficients $c_{in\lambda\mu}^k$ by recasting the learning problem into a low dimensional reproducing kernel Hilbert space (RKHS). This is done by first selecting a sparse set of $M$ atomic environments that best represent the structural variability in the dataset, following the subset of regressor approximation commonly used in GPR methods [34]. Upon including this selection within the construction of $K^{\lambda \otimes 1}$, we find that the dimensionality of the learning problem is increased by
<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) Visualization of the vector field rotation corresponding to the electron density response of a water molecule to an applied homogeneous electric field. Vectors at a distance of 2 Å from the closest atomic position are shown. (b) Percentage root mean square errors (% RMSE) of the density response prediction on five water conformers that are only rotated with respect to one another. The equivariant method of this work correctly yields the same error for these structures, while the method from Ref. 28 (SALTER) displays larger errors which are different for each conformer. (c) Box plot of all errors on the validation set of 900 structures, including rotated and distorted structures. While the average errors are small for both methods, the equivariant method of this work presents errors that are one order of magnitude smaller.</figcaption></figure>

roughly a factor of 3 with respect to the case of the electron density, consistently with having incorporated into the model the vectorial symmetry of the density response. To limit the curse of dimensionality, we introduce in this work a Gaussian damping factor that decreases the value of $M$ for increasing values of angular momentum, i.e., $M_λ = M_0 e^{-0.05λ^2}$, thus considerably reducing the final RKHS size. A thorough discussion about the learning method reported in the SI.

We test our method on a dataset of 1000 isolated water molecules presenting 100 rigid configurations that are randomly rotated, plus 900 that have a distorted geometry with a dipole moment aligned along the z-axis. This dataset has proven to be a useful benchmark for testing various equivariant methods [31, 33], as it allows us to simultaneously assess both the learning capability of the model and the effective incorporation of the desired rotational symmetries. Reference calculations of the density response function were obtained from DFPT, using the FHI-aims electronic-structure package. [11, 35] An example of the computed vector field representing the density

response between two mutually rotated water molecules is depicted in Fig. 1-a.

We train our model on a dataset of $N = 100$ randomly selected configurations – thus leaving the remaining 900 for testing – and compare the results with those that can be obtained using the SALTER method from Ref. 28. As both methods are derived as an extension of SALTED, we are able to perform this comparison with completely equivalent parameters, most notably a rather small number of sparse environments ($M_0 = 100$). As a consequence of learning independent models for each Cartesian component of the density response, SALTER is not capable of automatically recognizing pairs of vector fields that are simply rotated with respect to each other.

This problem is well illustrated in Fig. 1-b by a comparison of the prediction errors reported for five water molecule configurations that are rigidly rotated with respect to one another. While SALTER displays inhomogeneous errors that fluctuate around 1% RMSE, our method correctly produces entirely equivalent predictions that are almost 100 times more accurate. This highlights the effectiveness of symmetry adaptation in capturing the otherwise complex variations of the vectorial field under a three-dimensional rotation of the atomic structure. For completeness, we report in Fig. 1-c a summary of the prediction errors of all the 900 test configurations, 90% of which are aligned along the $\hat{z}$ direction. We observe a 50-fold improvement in prediction accuracy, when comparing SALTER with the fully equivariant approach of this work. We expect that if most of the molecules were not aligned, we would observe a more pronounced difference between the accuracy of both methods.

<table><caption>TABLE I. % RMSE in the predicted electron density responses across a test set of 100 liquid water structures and 100 naphthalene crystal structures using both SALTER and the symmetry-adapted vector field method presented in this work. Results are shown for different training set sizes $N$, highlighting the data-efficiency obtained through symmetry adaptation against the more data-hungry SALTER method.</caption><thead><tr><th rowspan="2">Method \ N</th><th colspan="2">Liquid Water</th><th colspan="2">Naphthalene</th></tr><tr><th>40</th><th>200</th><th>40</th><th>200</th></tr></thead><tbody><tr><td>SALTER</td><td>13.22</td><td>11.64</td><td>4.72</td><td>4.57</td></tr><tr><td><strong>This work</strong></td><td><strong>9.26</strong></td><td><strong>8.87</strong></td><td><strong>2.86</strong></td><td><strong>2.78</strong></td></tr></tbody></table>

We continue by comparing the performance of the equivariant approach and SALTER for more realistic systems, namely the liquid water and naphthalene datasets described in Ref. 28. The liquid water dataset consists of cubic cells of side 9.67 Å, containing 32 water molecules; the naphthalene structures are 2 × 2 × 1 supercells of the P2$_1$/a naphthalene crystal, with 8 molecules per supercell. Each dataset contains 500 structures, of which 100 are reserved for validation. The % RMSE in the electron density responses predicted by each method are summarized in Table I. Once again, we see a significant improvement
<figure><img src="image_3.png" /><figcaption>FIG. 2. Comparison between the reference and predicted electron density response ($\partial_k n_e$), along a given Cartesian direction $k$, for a test nanoparticle made of 55 gold atoms. Response profiles along the three Cartesian directions are obtained by averaging the vectorial field components along the other two complementary Cartesian directions. Dashed black lines: DFT reference. Red and blue lines: ML predictions.</figcaption></figure>

upon moving to the fully equivariant method, which is clearly manifested by requiring a smaller amount of data to achieve a comparable prediction accuracy. In particular, we find that predictions that can be obtained by training our model on $N = 40$ configurations are more accurate than those associated with training SALTER on $N = 200$ configurations. This property is observed not only in the size of the training set, but also in the number of sparse atomic environments used to define the RKHS of the learning problem. In fact, while SALTER models were defined using a relatively large value of $M$, i.e., $M = 3000$ for water and $M = 2000$ for naphthalene for every value of $\lambda$, our symmetry-adapted method uses for both systems $M = 500$ environments when $\lambda = 0$, falling to just $M = 42$ environments when $\lambda = 4$.

We conclude our discussion by applying the method on a dataset of gold nanoparticles of various sizes, reaching prediction targets that are challenging to calculate with full ab initio theory. These systems are expected to display a large polarization effect that is related to a non-local charge rearrangement when approaching the metallic limit of a vanishing HOMO-LUMO gap [36, 37]. For this reason, we include long-distance equivariant (LODE) features [38] into the construction of each $K^l$, a choice that has already proven successful to predict the electronic polarization of metallic slabs [23]. To investigate different size regimes, we generated a dataset of 255 relaxed nanoparticles, ranging from 28 to 139 atoms. We consider magic-size archetypes (icosahedra, decahedra, and fcc-truncated octahedra), as well as geometries obtained by randomly eliminating atoms from the cluster surface, and a successive geometry relaxation. Details about reference DFT calculations are reported in the SI.

We train our vector-field model on 200 randomly selected configurations and retain the remaining 55 for validation. The learning problem is then recast into a low dimensional space of $M_0 = 200$ sparse atomic environments, following the Gaussian decay already described for $M_\lambda$ up to $\lambda = 7$ spherical harmonics. A collective 15% RMSE is measured over the entire test set, which we find to result in highly accurate predictions. As an example, we report in Figure 2 the predicted density response in real space for a test 55 atoms cuboctahedron. We observe an excellent agreement with DFT by comparing the density response along the 3 Cartesian axes, reproducing the expected symmetry property of the selected structure along the $xy$-plane.

As a relevant application of the model previously trained, we consider the calculation of the electronic polarizability of the Au nanoparticles. When considering isolated systems – or systems presenting a distribution of electronic charge that vanishes before the cell periodic boundaries – the polarizability tensor can in fact be directly computed as the first moment of the density-response function:

$$ \alpha_{kk'} = \int_{\mathbb{R}^3} dr \, r_k \frac{\partial n_e(r)}{\partial E_{k'}}. \qquad (3) $$

Notably, analytical formulas can be derived that allows us to inexpensively compute $\alpha_{kk'}$ from the predicted $\lambda = 0$ and $\lambda = 1$ density-response coefficients; an explicit derivation is reported in the SI. We remark, however, that all spherical harmonics components of the density response (up to $\lambda = 7$ in this example) must still be considered during the training phase to account for the spatial overlap between basis functions.

An exemplary important question is that of assessing the capability of our model to reproduce the scaling law of the isotropic component of the polarizability tensor, defined as $\alpha_0 = \frac{1}{3}(\alpha_{xx} + \alpha_{yy} + \alpha_{zz})$, with an increasing system size [39]. For this purpose, we perform reference DFT calculations for an additional set of 30 test relaxed configurations that include larger nanoparticles than those used for training the model, up to 201 atoms. These serve as benchmarks of our prediction accuracy in the extrapolative regime. Five unrelaxed cuboctahedron geometries of 309, 561, 923, 1415, and 2057 atoms are finally included to the prediction targets, in order explore an asymptotically large size regime. For these, no reference value is computed because of the large computational cost.

Figure 3 reports prediction results for the whole set of 90 test structures considered. We find that our model can be used to accurately predict the expected value of $\alpha_0$ both within the training and extrapolation size domains with a comparable level of accuracy, achieving a collective error of only 1.2% RMSE. In agreement with the classical electrodynamic response of ideal metallic spheres, $\alpha_0$ is expected to follow a cubic scaling law with the nanoparticle radius, which manifests as a linear increase of the polariz-
<figure><img src="image_4.png" /><figcaption>FIG. 3. Isotropic components ($α_0$) of the polarizability tensors of 90 gold nanoparticles of increasing size, as obtained from the corresponding electron density response coefficients. Empty circles: DFT reference. Blue dots: ML prediction based on SOAP descriptors. Red dots: ML prediction based on LODE descriptors. Gray line: linear fit of $α_0$ on the reference DFT values. The vertical dashed line separates the domain of sizes used for training the vector-field model ($n_{atoms} < 140$) from the large-scale domain of sizes used to validate the model under extrapolation conditions ($n_{atoms} > 140$). No DFT reference is reported for the 5 largest nanoparticles, which are chosen as unrelaxed cuboctahedron configurations made of 309, 561, 923, 1415, and 2057 atoms.</figcaption></figure>

ability with the number of gold atoms. At a fundamental level, this result derives from the localization of the density response at the nanoparticle surface when approaching the classical metallic limit [40]. A local model that relies on nearsighted SOAP descriptors cannot capture this collective surface effect and shows a sublinear scaling of $α_0$ with respect to the number of gold atoms, as shown in Figure 3. Instead, when including long-range information into the model through LODE structural features, the expected linear behaviour is recovered throughout the whole extrapolation set, which includes a tenfold increase in size with respect to the training geometries (see also Fig. S1 in the SI). To estimate the computational speedup obtained with our model, we parallelize the calculation on 192 AMD EPYC 9654 2.4 GHz CPU-cores. The observed prediction times are of 1.5 seconds for $n_{atoms} = 309$ and of 11.6 seconds for $n_{atoms} = 2057$. These timings represent a remarkable acceleration compared to DFPT or finite-field DFT calculations, estimated to be > $10^3$ using comparable computational settings.

Accurate results are also found in predicting the anisotropic part of the electronic polarizability. In order to obtain a single measure of this quantity, we compute the Frobenius norm of the traceless $α$-tensor. While most configurations present a relatively small anisotropic component, a few nanoparticles display a strong anisotropic response. We find that our method can effectively reprod

uce both limits, including, once again, those configurations that fall outside the regime of sizes used for training. A collective error of 11.9 % RMSE on the anisotropic part is obtained throughout the test set. A more in-depth discussion is reported in the SI (section VI).

Following benchmarks ranging from molecular and periodic systems with diverse chemistry, to metallic clusters of increasing size that exhibit nonlocal polarizability, we conclude that our equivariant vector-field approach offers an accurate and scalable framework for predicting the electrostatic linear response of the electron density in a cost-efficient manner. In particular, we confirm the vastly superior performance of this method with respect to similar approaches that do not incorporate equivariance with the learning target. A key feature of our method is the streamlined incorporation of the vector-field rotational symmetry by a single Clebsch-Gordan sum over state-of-the-art spherical tensor kernels. We recognize that an iterative version of this procedure has been similarly adopted in the context of increasing the structural body-order within the construction of kernel-based ML potentials [41]. Looking ahead, unifying these two approaches could enhance the capabilities of kernel methods in predicting generic electronic-structure properties within a shared mathematical framework. Moreover, training on highly heterogeneous datasets will require the implementation of suitable strategies aimed at reducing the size of
the learning problem [42].

When it comes to the calculation of derived proper-
ties, the local nature of the density response allows us to
perform highly transferable predictions of dipolar polar-
izability tensors via suitable integral procedures. While
these tensors can directly enter the calculation of standard
Raman signals, higher order polarizability moments could
also be similarly computed to enable the calculation of
surface-enhanced Raman spectra [43, 44]. Furthermore,
having access to the real-space distribution of the re-
sponse function could make it possible to treat spatially
inhomogeneous induced polarizations, such as those en-
countered in tip-enhanced Raman scattering [45]. Beyond
optical dielectric properties, the prediction of spatially
inhomogeneous responses could also be used for tuning
dielectric-consistent DFT functionals [46], relevant for
improving the performance of DFT applied to semicon-
ductor interfaces. Finally, we envisage the application of
our method to model the induced nonlocal polarization at
electrified metal/electrolyte interfaces [47], thus enhanc-
ing the capabilities of first-principles atomistic simulations
of electrochemical systems and electrocatalytic processes.

SUPPORTING INFORMATION

Please see the Supporting Information for a thorough
derivation of symmetry-adapted kernels, a discussion of
the vector-field learning model, details about the dataset
generation and ML parameters, the calculation of the
polarizability tensor, complementary results.

Training and test geometries, together with the DFT in-
puts used to generate the reference density-response data
are free to download at https://zenodo.org/records/
14678822. An open source implementation of the density-
response learning method is available at https://github.
com/andreagrisafi/SALTED, including example input
files for training the model on the water monomer dataset.

ACKNOWLEDGMENTS

This work was supported by the French National Re-
search Agency under the France 2030 program (Grant
ANR-22-PEBA-0002). This article is also based upon
work from COST Action CA22154 - Data-driven Ap-
plications towards the Engineering of functional Ma-
terials: an Open Network (DAEMON) supported by
COST (European Cooperation in Science and Technol-
ogy). This work used the Dutch national e-infrastructure
with the support of the SURF Cooperative using grant
No. EINF-8985. The authors acknowledge HPC resources
granted by GENCI, France (resources of CINES, Grant
No. A0170910463).

* andrea.grisafi@sorbonne-universite.fr

[1] N. Raimbault, V. Athavale, and M. Rossi, Phys. Rev. Mater. **3**, 053605 (2019).

[2] H. Shang, N. Raimbault, P. Rinke, M. Scheffler, M. Rossi, and C. Carbogno, New J. Phys. **20**, 073040 (2018).

[3] T. Sohier, M. Calandra, and F. Mauri, Physical Review B **91**, 165428 (2015).

[4] P. Cudazzo, I. V. Tokatly, and A. Rubio, Physical Review B **84**, 085406 (2011).

[5] T. Dufils, G. Jeanmairet, B. Rotenberg, M. Sprik, and M. Salanne, Phys. Rev. Lett. **123**, 195501 (2019).

[6] C. Zhang, T. Sayer, H. Hutter, and M. Sprik, Journal of Physics: Energy **2**, 032005 (2020).

[7] G. Jeanmairet, B. Rotenberg, and M. Salanne, Chemical Reviews **122**, 10860 (2022).

[8] F. Che, J. T. Gray, S. Ha, and J.-S. McEwen, ACS Catalysis **7**, 551 (2017).

[9] F. Che, J. T. Gray, S. Ha, N. Kruse, S. L. Scott, and J.-S. McEwen, ACS Catalysis **8**, 5153 (2018).

[10] C. Ke, Z. Lin, and S. Liu, ACS Catalysis **12**, 13542 (2022).

[11] H. Shang, N. Raimbault, P. Rinke, M. Scheffler, M. Rossi, and C. Carbogno, New Journal of Physics **20**, 073040 (2018).

[12] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. **73**, 515 (2001).

[13] X. Gonze, S. Rostami, and C. Tantardini, Phys. Rev. B **109**, 014317 (2024).

[14] A. Grisafi, A. Fabrizio, B. Meyer, D. M. Wilkins, C. Corminboeuf, and M. Ceriotti, ACS Cent. Sci. **5**, 57 (2019).

[15] A. Fabrizio, K. Briling, A. Grisafi, and C. Corminboeuf, CHIMIA International Journal for Chemistry **74**, 232 (2020).

[16] B. Cuevas-Zuviría and L. F. Pacios, Journal of Chemical Information and Modeling **61**, 2658 (2021).

[17] P. B. Jørgensen and A. Bhowmik, npj Computational Materials **8**, 183 (2022).

[18] A. M. Lewis, A. Grisafi, M. Ceriotti, and M. Rossi, Journal of Chemical Theory and Computation **17**, 7203 (2021).

[19] S. K. Achar, L. Bernasconi, and J. K. Johnson, Nanomaterials **13**, 1853 (2023).

[20] J. A. Rackers, L. Tecot, M. Geiger, and T. E. Smidt, Machine Learning: Science and Technology **4**, 015027 (2023).

[21] A. Grisafi, A. M. Lewis, M. Rossi, and M. Ceriotti, Journal of Chemical Theory and Computation **19**, 4451 (2023).

[22] T. Lv, Z. Zhong, Y. Liang, F. Li, J. Huang, and R. Zheng, Phys. Rev. B **108**, 235159 (2023).

[23] A. Grisafi, A. Bussy, M. Salanne, and R. Vuilleumier, Phys. Rev. Mater. **7**, 125403 (2023).

[24] A. Grisafi and M. Salanne, The Journal of Chemical Physics **161**, 024109 (2024).

[25] C. Feng, Y. Zhang, and B. Jiang, “Efficient sampling for machine learning electron density and its response in real space,” (2024), arXiv:2410.04977 [physics.chem-ph].

[26] M. Domina, M. Cobelli, and S. Sanvito, Physical Review B **105**, 214439 (2022).

[27] R. Drautz, Physical Review B **102**, 024104 (2020).
[28] A. M. Lewis, P. Lazzaroni, and M. Rossi, The Journal of Chemical Physics **159**, 014103 (2023).

[29] Y. Zhang, C. Hu, and B. Jiang, The Journal of Physical Chemistry Letters **10**, 4962 (2019).

[30] N. A. Spaldin, Journal of Solid State Chemistry **195**, 2 (2012), polar Inorganic Materials: Design Strategies and Functional Properties.

[31] J. Nigam, M. J. Willatt, and M. Ceriotti, The Journal of Chemical Physics **156**, 014115 (2022).

[32] A. Glielmo, P. Sollich, and A. De Vita, Phys. Rev. B **95**, 214302 (2017).

[33] A. Grisafi, D. M. Wilkins, G. Csányi, and M. Ceriotti, Phys. Rev. Lett. **120**, 036002 (2018).

[34] J. Quiñonero Candela and C. E. Rasmussen, J. Mach. Learn. Res. **6**, 1939–1959 (2005).

[35] V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter, and M. Scheffler, Comput. Phys. Commun. **180**, 2175 (2009).

[36] T. Ishida, T. Murayama, A. Taketoshi, and M. Haruta, Chemical Reviews **120**, 464 (2020).

[37] R. Sinha-Roy, P. García-González, and H.-C. Weissker, Nanoscale **12**, 4452 (2020).

[38] A. Grisafi and M. Ceriotti, J. Chem. Phys. **151**, 204105 (2019).

[39] A. Ojha, S. S. Bulusu, and A. Banerjee, Artificial Intelligence Chemistry **2**, 100080 (2024).

[40] D. R. Snider and R. S. Sorbello, Phys. Rev. B **28**, 5702 (1983).

[41] F. Bigi, S. N. Pozdnyakov, and M. Ceriotti, The Journal of Chemical Physics **161**, 044116 (2024).

[42] N. Lopanitsyna, G. Fraux, M. A. Springer, S. De, and M. Ceriotti, Phys. Rev. Mater. **7**, 045802 (2023).

[43] G. S. Kedziora and G. C. Schatz, Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy **55**, 625 (1999).

[44] X. X. Han, R. S. Rodriguez, C. L. Haynes, Y. Ozaki, and B. Zhao, Nature Reviews Methods Primers **1**, 87 (2022).

[45] Y. Litman, F. P. Bonafé, A. Akkoush, H. Appel, and M. Rossi, The Journal of Physical Chemistry Letters **14**, 6850 (2023), 2211.09475.

[46] J. Zhan, M. Govoni, and G. Galli, Journal of Chemical Theory and Computation **19**, 5851 (2023).

[47] L. Andersson, M. Sprik, J. Hutter, and C. Zhang, Angewandte Chemie International Edition **64**, e202413614 (2025).