Andrea Grisafi
email: andrea.grisafi@ens.psl.eu

Alan M Lewis
email: alan.lewis@mpsd.mpg.de

Mariana Rossi

Michele Ceriotti

¶ †pasteur

Electronic-structure properties from atom-centered predictions of the electron density

come

Introduction

The electron density ρ plays a special role in electronic structure methods, besides being a well-defined physical observable that can be experimentally measured. [START_REF] Koritsanszky | Chemical Applications of X-ray Charge-Density Analysis[END_REF] In particular, within density-functional theory (DFT), ρ provides direct access to not just the ground-state electronic energy, but in principle to all of the ground-state properties of a system.

In the last few years, many efforts have been devoted to the integration of supervised learning techniques within state-of-the-art electronic-structure methods, aiming at accelerating the calculation of properties beyond ground-state total energies and forces. [START_REF] Hermann | Deep-neural-network solution of the electronic Schrödinger equation[END_REF][START_REF] Westermayr | Perspective on integrating machine learning into computational chemistry and materials science[END_REF][START_REF] Kirkpatrick | Pushing the frontiers of density functionals by solving the fractional electron problem[END_REF][START_REF] Kulik | Roadmap on Machine learning in electronic structure[END_REF] Different quantities have been considered as the target of machine learning (ML) models, including molecular multipoles [START_REF] Bereau | Transferable Atomic Multipole Machine Learning Models for Small Organic Molecules[END_REF][START_REF] Staacke | Kernel charge equilibration: efficient and accurate prediction of molecular dipole moments with a machine-learning enhanced electron density model[END_REF][START_REF] Veit | Predicting Molecular Dipole Moments by Combining Atomic Partial Charges and Atomic Dipoles[END_REF] and polarizabilities, 9 electronic density of states and bandgaps, [START_REF] Pilania | Multi-fidelity machine learning models for accurate bandgap predictions of solids[END_REF][START_REF] Ben Mahmoud | Learning the electronic density of states in condensed matter[END_REF][START_REF] Westermayr | Machine Learning for Electronically Excited States of Molecules[END_REF] kinetic density functionals, [START_REF] Snyder | Finding Density Functionals with Machine Learning[END_REF][START_REF] Meyer | Machine Learning Approaches toward Orbital-free Density Functional Theory: Simultaneous Training on the Kinetic Energy Density Functional and Its Functional Derivative[END_REF][START_REF] Ryczko | Toward Orbital-Free Density Functional Theory with Small Data Sets and Deep Learning[END_REF] exchange-correlation potentials, [START_REF] Zhou | Toward the Exact Exchange-Correlation Potential: A Three-Dimensional Convolutional Neural Network Construct[END_REF][START_REF] Schmidt | Machine Learning the Physical Nonlocal ExchangeâCorrelation Functional of Density-Functional Theory[END_REF][START_REF] Suzuki | Machine learning exchange-correlation potential in time-dependent density-functional theory[END_REF] single-particle wavefunctions [START_REF] Welborn | Transferability in Machine Learning for Electronic Structure via the Molecular Orbital Basis[END_REF][START_REF] Schütt | Unifying Wavefunctions[END_REF][START_REF] Gastegger | A deep neural network for molecular wave functions in quasi-atomic minimal basis representation[END_REF] and Hamiltonians. [START_REF] Schütt | Unifying Wavefunctions[END_REF][START_REF] Westermayr | Physically inspired deep learning of molecular excitations and photoemission spectra[END_REF][START_REF] Nigam | Equivariant representations for molecular Hamiltonians and N-center atomic-scale properties[END_REF] For problems where DFT is the method of choice, it is attractive to think that one could access all ground-state properties of a system at once by simply being able to predict its real-space electronic density. However, it is not always obvious whether one can actually benefit from computing derived electronic properties from the predicted densities. In fact, under many circumstances it appears more convenient to target the prediction of the property at hand directly, without going through intermediate electronic-structure ingredients. [START_REF] Grisafi | Transferable Machine-Learning Model of the Electron Density[END_REF] We have recently shown that using the electron density as a stepping stone to predict the electronic energy of a system can be beneficial when working in highly extrapolative regimes, which are specifically associated with a substantial increase of the system size. [START_REF] Lewis | Learning Electron Densities in the Condensed Phase[END_REF] These are the most sought-after regimes for ML models, since they are where ML can substantially decrease computational cost, thus allowing a quantum-mechanical description of systems that were unattainable just a few years ago. This good extrapolative behaviour is to be expected, because computing the electronic energy as a functional of ρ can effectively recover electrostatic and non-local effects, without the need to add an explicit long-range description to the model. For these reasons, designing efficient and accurate learning algorithms that specifically target the prediction of the electron density is becoming an increasingly attractive possibility both to provide useful insights on the electronic structure of the system and to complement the predictions of local ML potentials.

Various density-learning models currently available primarily differ by the nature of the difficult task of discretizing the scalar field ρ for numerical handling by a computer.

Algorithms can work with a finite number of real-space grid points, [START_REF] Alred | Machine Learning Electron Density in Sulfur Crosslinked Carbon Nanotubes[END_REF][START_REF] Chandrasekaran | Solving the Electronic Structure Problem with Machine Learning[END_REF][START_REF] Jørgensen | Neural Message Passing Network for Accurate Charge Density Prediction[END_REF][START_REF] Jørgensen | Graph neural networks for fast electron density estimation of molecules, liquids, and solids[END_REF] or with expanding the density over a suitable set of basis functions, [START_REF] Grisafi | Transferable Machine-Learning Model of the Electron Density[END_REF][30][START_REF] Bogojeski | Quantum chemical accuracy from density functional approximations via machine learning[END_REF][START_REF] Rackers | Cracking the Quantum Scaling Limit with Machine Learned Electron Densities[END_REF] for example. When it comes to the learning architectures, both neural networks and kernel-based approaches have been adopted over the years. While the former are more suited to handle a large number of training data, the latter can be used to interpret more transparently the relation between the learning model and the underlying physical problem, thus facilitating the prediction of complex targets such as the electron density. This paper addresses current shortcomings of kernel-based models of ρ that effectively limits its training in large datasets. We propose an extension of the symmetry-adapted learning of three-dimensional electron densities (SALTED) method, [START_REF] Lewis | Learning Electron Densities in the Condensed Phase[END_REF] where the density is discretized through a linear atom-centered expansion analogous to that of density-fitting (DF) approximations routinely used in quantum chemistry codes. [START_REF] Baerends | Self-consistent molecular Hartree-Fock-Slater calculations I. The computational procedure[END_REF][START_REF] Weigend | Accurate Coulomb-fitting basis sets for H to Rn[END_REF][START_REF] Golze | Local Fitting of the Kohn-Sham Density in a Gaussian and Plane Waves Scheme for Large-Scale Density Functional Theory Simulations[END_REF][START_REF] Ren | Resolution-of-identity approach to Hartree-Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions[END_REF] Within SALTED, a kernel-based parametrization of the learning problem is provided through a local representation of the atomic structure that mirrors the rotational symmetry of the atomic density components. The atom-centered nature of the density expansion is a critical aspect to ensure the transferability of the SALTED predictions across atomic environments that share similar local structural patterns. The non-orthogonal character of common DF basis sets, however, implies that particular care must be given to the spatial overlap between any pair of off-centered basis functions, which has the effect of coupling all the atomic density components together. In fact, the problem dimensionality depends not only on the number of atomic environments used to train the model but also on the basis set size associated with each atomic center. This downside of the method hinders the calculation of explicit solutions at an increasing number of training structures, effectively limiting the application of SALTED to relatively small datasets. Furthermore, as the number of training atomic environments can grow very quickly, suitable dimensionality reduction schemes must be adopted in order to limit the size of the fitting problem. In this respect, converging the prediction accuracy with respect to the number of sparse atomic inputs is known to be a particularly hard numerical problem that requires the implementation of suitable expedients. [START_REF] Foster | Stable and Efficient Gaussian Process Calculations[END_REF] These shortcomings become particularly severe whenever a large number of atomic density components are required to train the model, e.g., when targeting the prediction of the electron density of proteins [START_REF] Fabrizio | Learning (from) the Electron Density: Transferability, Conformational and Chemical Diversity[END_REF] or generic disordered, liquid systems.

In this paper, we overcome the aforementioned limitations in two ways: i) we make use of a featurization of the atomic structure expressed explicitly in the reproducing kernel Hilbert space (RKHS) of symmetry-adapted sparse kernel approximations, thus enabling a systematic convergence of the prediction accuracy with the number of sparse atomic environments, and ii) we implement an iterative gradient-based minimization of the SALTED loss function, which allows us to bypass the computational bottleneck associated with the inversion of prohibitively large regression matrices. Taking bulk liquid water as an example, we show that the combination of these two ingredients allow us to apply SALTED to arbitrarily complex systems, paving the way for the inexpensive calculation of electronic structure properties and accurate electronic energies.

Theory

Consider a generic density-fitting approximation of the electron density of a given structure A, expressed as a linear expansion over a set of atom-centered basis functions φ anλµ :

ρ DF A (r) = i,anλµ,U c anλµ (A i ) φ anλµ (r -R i + U) . ( 1
)
where r is the spatial coordinate, R i labels the position of atom i, U is the cell translation vector, assuming the system is periodic, and c anλµ (A i ) are the non-orthogonal expansion coefficients. The basis φ anλµ is defined as the atom-centered product of spherical harmonics Y µ λ , with λ and µ labelling the azimuthal and projection indices, and a set {n} of radial functions R aλ n which depend on the possible choices of atomic types a and λ values.

Similarly
to what is reported in previous works, [START_REF] Grisafi | Transferable Machine-Learning Model of the Electron Density[END_REF][START_REF] Lewis | Learning Electron Densities in the Condensed Phase[END_REF][START_REF] Fabrizio | Electron Density Learning of Non-Covalent Systems[END_REF] we then provide an atom-centered approximation to the expansion coefficients through the definition of a covariant similarity measure between atomic environments, which takes into account the transformation properties of spherical harmonics under rotations:

c anλµ (A i ) ≈ {M anλ j } |µ |≤λ b anλ jµ k anλ µµ (A i , M anλ j ). (2)
Here, b anλ jµ are the regression weights and k anλ µµ (A i , M Inserting Eq. (2) into Eq. ( 1), a ML approximation of the electron density ρ ML can be defined to linearly depend on the vector of regression weights b M , where we use the subscript M to indicate that the vector corresponds to the full set of reference environments. The vector b M can then be obtained from the minimization of a quadratic loss function that measures the integrated error over the N training structures between ρ ML and the reference DF definition of Eq. ( 1):

(b M ) = N A=1 u.c. dr ρ ML A (r; b M ) -ρ DF A (r) 2 + η b T M K M M b M , (3)
where in the second line we introduced a regularization term modulated by a hyperparameter η. The matrix K M M includes all the symmetry-adapted kernels k anλ µµ (M anλ i , M anλ j ) that couple the sparse set of atomic environments {M anλ j } and it is thus defined to be block diagonal in the possible basis function types {anλ}.

RKHS formulation

According to standard subset of regressors 41 (SoR) approaches, selecting the number of sparse atomic environments M anλ as a proper subset of the total number of training environments implies approximating the full kernel matrix of the SALTED problem as follows:

K N N ≈ K N M K -1 M M K T N M (4)
where we used the subscript N to indicate that the dimension of K N N is given by the sum of the number of DF basis functions associated with each atom of the training set. Most of the time, the matrix K M M is found to be ill-conditioned, yielding severe numerical instabilities when solving the regression problem. Commonly adopted solutions to circumvent this issue consist in introducing an ad hoc jitter term that provides a lower bound to the eigenvalues of the matrix to be inverted. [START_REF] Lewis | Learning Electron Densities in the Condensed Phase[END_REF] While effectively fixing the numerical instabilities, this stratagem carries the drawback of hindering the systematic convergence of the SoR approximation when taking the limit of M anλ towards the total number of training environments. To solve this problem, we first consider a truncation of the eigen-decomposition of the K M M matrix. Given that K M M is block-diagonal in the basis function type, for each triplet (anλ) we can write

K anλ N N ≈ K anλ N M   D anλ d v d M (v d M ) T λ d   (K anλ N M ) T , (5)
where λ d are the non-negligible eigenvalues of K anλ M M and v d M the corresponding eigenvectors, so that D anλ indicates the final truncated dimension. Note that each outer product

v d M (v d M ) T
conserves the transformation properties of K anλ M M , so that, individually, each eigenvector v d M has the same rotational symmetry as the corresponding basis functions (aλn). From the equation above, we can then define the reproducing kernel Hilbert space (RKHS) of the sparse kernel approximation as the space spanned by the following feature vectors:

Ψ anλ N D ≡ K anλ N M V anλ M D (Λ anλ DD ) -1 2 , (6)
which, consistently with Mercer's theorem, [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF] gives

K anλ N N ≈ K anλ N M (K anλ M M ) -1 (K anλ N M ) T ≈ Ψ anλ N D (Ψ anλ N D ) T . (7)
Finally, we can reformulate the density-learning problem as a linear regression task parametrized according to the feature vectors Ψ anλ N D . To do so, we can first invert Eq. ( 6) to approximate the K anλ N M matrices as

K anλ N M = Ψ anλ N D (Λ anλ DD ) 1 2 (V anλ M D ) T , (8)
where we made use of the unitary nature of the matrix of eigenvectors,

V anλ M D (V anλ M D ) T = 1 anλ M M .
Then, a sparse RKHS approximation of the density expansion coefficients can finally be obtained by rewriting Eq. ( 2) as follows

c anλµ (A i ) ≈ D anλ d banλ d ψ anλ d (A i ; anλµ) , (9)
where, for each basis function type (anλ), we defined a new vector of regression weights as

banλ D ≡ (Λ anλ DD ) 1 2 (V anλ M D ) T b anλ M . (10)
The so derived RKHS reformulation of SALTED allows us to effectively make use of a sparse symmetry-adapted kernel approximation of the density expansion coefficients that avoids the numerical instabilities associated with normal equations solutions. [START_REF] Foster | Stable and Efficient Gaussian Process Calculations[END_REF]

Iterative solution

To practically solve the SALTED-RKHS problem, it is convenient to define the global problem dimensionality D as the sum of all the truncated dimensions D anλ associated with each basis function type. In so doing, a single vector of regression weights bD constructed by stacking together each of the individual banλ D can be used to compactly write the SALTED-RKHS approximation of the density coefficients as

c anλµ (A i ) ≈ bT D Ψ D (A i ; anλµ) , (11)
where the global feature vector Ψ D (A i ; anλµ) is defined to be diagonal in the various basis function types. When inserted into Eq. ( 1), this allows us to rewrite the loss function of Eq. ( 3) as follows:

( bD ) =(Ψ N D bD -c DF N ) T S N N (Ψ N D bD -c DF N ) + η bT D bD . (12)
with c DF N the vector of training DF expansion coefficients. The matrix S N N , which is block diagonal in the training set, contains the spatial overlaps between the basis functions of each structure, φ|φ , and is responsible for coupling all the basis function types and therefore all the atomic density components together. Importantly, the sparsity pattern of Ψ N D differs from that of S N N , implying that the product

Ψ T N D S N N Ψ N D is, generally, a full D × D matrix.
A schematic representation of this is reported in Fig. 1.

An explicit minimization of the loss function of Eq. ( 12) would be hindered by the quickly growing dimensionality of the problem when considering an increasing number of sparse atomic environments, which would make the inversion of the resulting regression matrix prohibitively expensive when considering a large D and/or a large basis set size. For this reason, we implement a basic conjugate gradient (CG) algorithm to iteratively minimize the loss function of Eq. ( 12) directly. We rely on the analytical definition of the gradient of

∂ ( bD ) ∂ bD = Ψ T N D S N N (Ψ N D bD -c DF N ) + η bD , (13)
which is updated at each CG step from the knowledge of the Hessian matrix:

H DD = Ψ T N D S N N Ψ N D + η1 DD . ( 14
)
Note that as the dimension D can easily be > 10 5 when treating complex datasets, precomputing H DD would imply running into similar issues to those associated with solving the regression problem explicitly. For this reason, we estimate the gradient and curvature on the minimization surface at each CG step by directly computing the projection of Eq. ( 14) on the CG search direction. In this work, an efficient implementation of the latter quantity is obtained by parallelizing the sum over the training set, preloading, for each configuration, the associated dense block of the global overlap matrix S N N and RKHS descriptor, and exploiting the high sparsity pattern of Ψ N D (in our dataset about 5% of entries are non-zero)

to perform sparse matrix multiplications.

Orthogonal approach

As a further justification of the need to tackle the non-orthogonal regression problem, we consider an alternative path that could in principle be undertaken to bypass the computational difficulties discussed above. A suitable orthogonal transformation could be introduced in order to recast the SALTED problem into as many independent regression tasks as the number of the different basis function types {anλ}. In particular, one could formally rewrite the DF approximation of the electron density in terms of set of orthogonal basis functions, defined as follows

φk (r) = k S -1 2 k k φ k (r) , (15)
10 with S -1 2 the orthogonalization matrix. This transformation is known as Löwdin orthogonalization [START_REF] Löwdin | On the Non-Orthogonality Problem Connected with the Use of Atomic Wave Functions in the Theory of Molecules and Crystals[END_REF] and it is typically used in the context of single-particle theories of the electronic wavefunction to work in the space of atomic natural orbitals, [START_REF] Almlöf | In Atomic Natural Orbital (ANO) Basis Sets for Quantum Chemical Calculations[END_REF] as well in the localization procedure of maximally-localized Wannier functions. [START_REF] Marzari | Maximally Localized Generalized Wannier Functions for Composite Energy Bands[END_REF][START_REF] Marzari | An introduction to maximally-localized Wannier functions. Highlight of the Month[END_REF] Under this transformation, the orthogonal expansion coefficients can be obtained from the non-orthogonal DF coefficients as

ĉ = S 1 2 c . ( 16
)
Unlike hierarchical orthogonalization schemes, such as Gram-Schmidt, the Löwdin orthogonalization preserves both the locality of the basis functions about the atomic positions and the symmetry properties of spherical harmonics. As such, for each of the distinct triplets (anλ) that labels the new set of orthogonal basis functions, a ML approximation of Eq. ( 16) can be defined to follow the same sparse RKHS formulation already reported in Eq. ( 9).

As the basis functions are now orthogonal to each other, S N N = 1 N N in Eq. ( 12) and the regression problem decouples, leaving as many independent solutions as the number of orthogonal basis function types {anλ}. This means that one can afford to compute the regression weights for each (anλ) explicitly, i.e., baλn

D = ((Ψ aλn N D ) T Ψ aλn N D + η1 DD ) -1 (Ψ aλn N D ) T ĉaλn N . (17)
Note that because the learning target still remains the entire density field, the value of the regularization parameter η is defined to be the same for each independent solution.

In spite of the clear computational advantages, the presented orthogonal approach has two main downsides. First, since the basis set of each structure now depends on the specific atomic configuration that determines the Löwdin orthogonal transformation, the resulting regression task is effectively required to learn this transformation, making the method more data hungry. Moreover, reconstructing the density field from the predicted density coefficients necessarily requires to undo the orthogonal transformation as c = S -1 2 ĉ, adding a considerable computational burden at the prediction level. In fact, especially when targeting the prediction of systems with a large number of atoms, the calculation of S -1 2 can be particularly expensive, undermining the rationale for applying a ML model in the first place.

Accuracy metrics

As already discussed in Ref. 25, the mean spherical average of the density over the dataset, ρ0 (r), is used as a baseline field for the learning of ρ(r). This implies that very much like when doing regression of electronic energies, we can measure the root mean square error associated with the density predictions carried out on any given test set as a fraction of the standard deviation of the reference test densities about ρ0 (r), i.e.,

%RMSE =

A

dr |ρ ML A (r) -ρ DF A (r)| 2 A dr |ρ DF A (r) -ρ0 (r)| 2 × 100 (18)
This choice carries the advantage of providing a measure of the learning capability of the method, rather than an absolute estimation of the accuracy of the predicted density. In practice, Eq. ( 18) can be efficiently computed as follows:

%RMSE = A (c ML A -c DF A ) T S A (c ML A -c DF A ) A (c DF A -c0 ) T S A (c DF A -c0 ) × 100 ( 19
)
where c0 is the vector of spherical (λ = 0) coefficients associated with the expansion of the average density ρ0 (r). The overlap matrices S A of the test structures are set to the identity matrix when considering the prediction of orthogonal coefficients.

To provide an estimate of the absolute accuracy of the predicted density we use the mean absolute error (MAE) of the density as a fraction of the total electronic charge:

%MAE = 1 N e dr ρ ML (r) -ρ QM (r) × 100, (20)
where ρ QM is the reference Kohn-Sham density. This error is computed using FHI-aims's dense logarithmic real-space grids centred on each atom, ensuring an accurate integration of the density close to the nuclei. [START_REF] Blum | Ab Initio Molecular Simulations with Numeric Atom-Centered Orbitals[END_REF] 3 Results

Electron density of bulk liquid water

To test the enhanced SALTED framework, we begin by considering a dataset consisting of 1000 configurations of 32-molecule supercells of bulk liquid water. This example presents enough complexity to require training the model on a large number of atomic density components, hence also requiring a large value of M , highlighting the importance of the methodological advances discussed here. In the process, we also compare the performance of the method with the alternative orthogonal approach outlined in Sec. II-C. The reference density coefficients are computed at the DFT/LDA level using FHI-aims's tight settings by projecting the Kohn-Sham electron density on a set of auxiliary Slater-type atomic functions as implemented in FHI-aims. [START_REF] Ren | Resolution-of-identity approach to Hartree-Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions[END_REF] This results in a density field that is represented using up to λ = 8 spherical harmonics. The corresponding DF strategy and the assessment of the DF error is discussed in detail in Ref.; 25 in this case the error introduced to the total energy by DF is on the order of 1 µeV/atom.

Given the additional complexity of treating an extended condensed-phase system, we adjust the SALTED model according to the spatial resolution of the radial functions involved in the expansion of the density. For that, several λ-SOAP kernels are initially computed for various Gaussian widths σ and radial cutoffs r c , in such a way that a ratio of r c /σ = 10 is always conserved. In particular, we consider six different cutoffs which span r c = (2.0, 3.0, 4.0, 5.0, 6.0) Å and compute RKHS descriptors as in Eq. ( 6 measure the distance between pairs of atomic environments. [START_REF] Imbalzano | Automatic Selection of Atomic Fingerprints and Reference Configurations for Machine-Learning Potentials[END_REF] For each of the aforementioned models, we then perform orthogonal-learning predictions which are cross-validated for each type of atomic functions (anλ) on a two-fold 500/500 random partition of the dataset. Finally, the models that give minimal errors for each basis function type are selected to construct an optimal representation of the atomic structure to be used both in orthogonal and nonorthogonal approaches. Note that we will use throughout the same sparse environments selections for each (anλ). training structures, which carry a total of 19200 atomic density components. While the orthogonal approach yields optimal predictions associated with a RMSE of 7.0%, making use of the fully coupled non-orthogonal framework allows us to bring the error down to 3.3%. Note that this level of accuracy corresponds to a mean absolute error of 0.36% when measured as a fraction of the total electronic charge. In spite of the orthogonal approach being computationally more convenient to train, these results highlight the importance of adopting a fixed set of basis functions in order to solely focus regression effort on the density-based variations of the learning target.

Prediction of electronic-structure properties

Having generated the reference data at the DFT level, the predicted densities can be used to compute the electronic energy of the system. Within the Kohn-Sham DFT formalism, however, one can only have direct access to the electrostatic and exchange-correlation energy, while the kinetic energy must be expressed in terms of the Kohn-Sham molecular orbitals.

Therefore, we obtain the total energy by computing the Kohn-Sham potential as a functional of ρ ML (r) and solving the resulting eigenvalue problem. Doing so requires performing a single diagonalization step of the Kohn-Sham Hamiltonian, which can be done with almost linear-scaling cost thanks to dedicated libraries used in FHI-aims. [START_REF] Yu | ELSI -An open infrastructure for electronic structure solvers[END_REF] As shown in Fig. 2b, this procedure can in fact be used to yield extremely accurate electronic energies. In particular, for the maximum training set size adopted the predicted densities obtained through the nonorthogonal framework are associated with an error of 0.13meV/atom, greatly outperforming the accuracy that can be achieved through the orthogonal-learning framework. Interestingly, this level of accuracy is made possible by a large error cancellation between the kinetic energy T and the potential energy V of the system. As reported in the Supplementary Information (SI), computing separately the kinetic, electrostatic and exchange-correlation energy shows in fact that each of these individual terms are associated with errors that are from one to two orders of magnitude larger than those associated with the total electronic energy. The observed error cancellation is related to the fact that the first functional derivative of the total energy is zero when evaluated at the reference self-consistent density. In particular, when considering a functional expansion of the total energy E in increasing powers of the density error ∆ρ, the linear term in the error is expected to vanish at self-consistency, so that the leading order in the total energy error ∆E is in fact quadratic in ∆ρ. A formal derivation together with a numerical test of the quadratic relation ∆E ∝ ∆ρ 2 is reported in the SI.

In principle, being able to use the predicted density to solve the Khon-Sham eigenvalue problem allows us to access all kind of electronic-structure properties. In Fig. 3 we report an example of this by showing how the predicted HOMO-LUMO gaps of the 500 test water structures compared with the reference DFT values. Because of its non-local nature, the accurate prediction of this quantity is known to be a particularly hard problem. [START_REF] Schütt | Unifying Wavefunctions[END_REF][START_REF] Nigam | Equivariant representations for molecular Hamiltonians and N-center atomic-scale properties[END_REF] As shown in the Figure, in spite of having learned a local property using local structural information only, our density-based predictions correlate well with the reference DFT values, resulting in a low root mean square error of 84 meV (about 2% of the DFT band gap). We observe a small systematic overestimation of the gaps, arising primarily from a systematic underestimation of the HOMO energy. A 3D visualization of the HOMO-LUMO states of a given test water frame, shown in Fig. 4, demonstrates an excellent agreement between the reference and predicted Kohn-Sham orbitals.

The QM9 dataset

As an even more challenging test, we learn the electron density of the QM9 dataset, containing more than 130,000 different molecules containing H, C, N, O and F atoms. Having demonstrated the superior performance of the non-orthogonal learning framework against the computationally more convenient Löwdin approach, we will focus from here on exclusively on the former method. The vast chemical heterogeneity of the dataset requires dealing at the same time with a great number of sparse atomic environments and training structures, making it virtually impossible to use a naive implementation of SALTED where the regression exercise is carried out via the inversion of a kernel matrix. In testing the accuracy of our predictions we rely on the same partitioning of the dataset reported in Ref., [START_REF] Jørgensen | Graph neural networks for fast electron density estimation of molecules, liquids, and solids[END_REF] i.e., using the same random selection of 10,000 structures as a test set, while using the remaining configurations for training the model. In order to directly compare our results with what is In this case the reference density coefficients are computed at the DFT/PBE level, again using FHI-aims's tight settings. Although we observe superior learning performances when compared with the dataset of liquid water, the huge variance of the electron density in the QM9 dataset implies a lower absolute accuracy with respect to the previous example. In fact, as shown in the lower panel of the Figure, absolute density errors integrated on a real-space grid are distributed around a mean absolute error of 0.45%, which is about 25% larger than those observed for liquid water. When compared with the results of Ref. [START_REF] Jørgensen | Graph neural networks for fast electron density estimation of molecules, liquids, and solids[END_REF] for the very same QM9 test molecules, we find that our errors are about 1.5 times larger; this is still remarkable when considering that we used only 6% of the training set and that we predict the all-electron density, whereas

Ref. [START_REF] Jørgensen | Graph neural networks for fast electron density estimation of molecules, liquids, and solids[END_REF] uses almost the entire QM9 dataset and uses the pseudo-valence density as target.

In Figure 6 we also report learning curves and error histogram for the derived total energy predictions, showing a mean absolute error that is brought down to 1.57 kcal/mol, with 65% of our predictions falling within chemical accuracy. Although not reported in the Figure for visual purposes, we find that there are 196 molecules with an error >10 kcal/mol up to a maximum of ∼228 kcal/mol. These outliers are typically associated with specific atomic environments (such as 3-or 4-membered rings, 5-membered rings containing heteroatoms, or alkyne groups) which are not well represented in the training set, consistent with having used only a tiny fraction of the available data.

Learning curves for both densities and derived energies clearly show that there is still margin for improvement by increasing both M and N , so that our method is likely to surpass the performance of Ref. [START_REF] Jørgensen | Graph neural networks for fast electron density estimation of molecules, liquids, and solids[END_REF] and M and overcoming the need to use massive computational resources will be subject of future investigation.

Conclusions

The data-driven prediction of the electron density of a system is a challenging computational task which requires balancing the need for an optimal discretization of the scalar field with a regression model that is at the same time highly transferable and that can be easily integrated with electronic-structure programs. Adopting an atom-centered decomposition of the scalar field represents an optimal choice to simultaneously meet these requirements, but the implementation of an efficient machine-learning algorithm carries exceptional computational challenges. We have shown how to tackle these challenges through an iterative minimiza-tion scheme of the density loss function which makes use of an equivariant sparse kernel approximation of the expansion coefficients. In doing so, we bypassed the need to invert regression matrices, increasing the capabilities of the method by 1-2 orders of magnitude with respect to previous implementations. [START_REF] Lewis | Learning Electron Densities in the Condensed Phase[END_REF] It is worth mentioning that an equivariant neural network approach that works on a similar density-fitting basis has recently been proposed. [START_REF] Rackers | Cracking the Quantum Scaling Limit with Machine Learned Electron Densities[END_REF] In this method, the loss function is defined via the sum of square errors made on the individual expansion coefficients, rather than on the density itself. Although powerful, this approach may be limited in accuracy, as the coupling between pairs of non-orthogonal density components is only implicitly obtained through the network architecture, overshadowing the role of the correlations between the expansion coefficients in the final prediction of the scalar field. Conversely, we have shown that that within a kernel-based method the correlations between non-orthogonal coefficients are solely brought by the spatial overlap between the basis functions, consistently with the pure geometric nature of the coupling between off-centered density components. As an alternative approach, we have suggested that a decoupling of the learning problem could be obtained through a Löwdin orthogonalization strategy analogous to that routinely used for the atomic orbitals in quantum-chemistry calculations. However, because of the system-dependent nature of the orthogonalization procedure, this strategy is outperformed by the (fully coupled) non-orthogonal framework, especially when using the predicted density to compute derived electronic energies.

As long as a single diagonalization step of the Kohn-Sham Hamiltonian can be performed at low computational cost, [START_REF] Yu | ELSI -An open infrastructure for electronic structure solvers[END_REF] the possibility of predicting the electron density to access all kinds of ground-state electronic-structure properties is particularly attractive, as it sidesteps all intermediate learning exercises that might be carried out to predict quantities such as the density of states, band-gaps, wavefunctions and the Hamiltonian matrix itself. Moreover, the presented protocol could be used as an alternative to standard machine-learning potentials, where the total energy is learned and predicted directly from the atomic coordinates; while the computational cost of going through the density is necessarily higher, it allows one to more easily capture long-range effects that come from a more or less local integral of ρ and, importantly, incorporates the permanent electrostatics exactly. As already shown in Ref., [START_REF] Lewis | Learning Electron Densities in the Condensed Phase[END_REF] this is extremely convenient when predicting the properties of systems that are much larger than those used during the training phase.

In this respect, a relevant next step will be investigating the role of long-range descriptors in complementing the capabilities of the method to predict all kind of physical interactions.

For example, making use of far-field representations such as LODE [START_REF] Grisafi | Multi-scale approach for the prediction of atomic scale properties[END_REF] to learn the electron density would help when reproducing long-range phenomena such as induced polarization and charge-transfer effects, which could not be captured through a local learning of the electron density around the atoms. Finally, the accurate prediction of ρ could be used to feed existing orbital-free kinetic-energy functionals, which are typically limited by the quality of the density, as well as to correct specific classes of KS-DFT calculations which present pathological density-driven errors.

Figure 1 :

 1 Figure 1: Schematic representation of the sparsity pattern of the RKHS feature matrix Ψ N D (left) and basis-set overlap matrix S N N (right) entering the iterative minimization of the electron-density loss function. Ψ N D can be built from separate blocks, each associated with a separate type of atomic basis functions.

Figure 2 :

 2 Figure 2: Learning curves associated with the prediction of 500 randomly selected liquid water structures using models constructed with an increasing number of sparse environments M . (top) % RMSE of the predicted electron density as a function of the number of training structures. (bottom) Absolute RMSE in meV/atom of the electronic energies obtained by feeding the predicted electron densities into the density functional used to generate the reference data. Full and dashed lines correspond to non-orthogonal and orthogonal learning models, respectively.

Fig. 2a reports

  Fig.2areports learning curves for the electron-density prediction of 500 randomly selected liquid water frames, comparing regression exercises carried out at an increasing number of sparse atomic environments M . Both for orthogonal and non-orthogonal learning approaches,

Figure 3 :

 3 Figure 3: Correlation plot between the reference HOMO-LUMO gaps of 500 test water frames and those computed from the predicted electron density as obtained from a non-orthogonal learning exercise at N = 200 and M = 2000. The root mean square error in the predicted gaps is 84 meV.

Figure 4 :

 4 Figure 4: (top) Two pairs of HOMO isosurfaces at ±0.25 and ±0.05 a.u. associated with the reference DFT calculation (a) and those derived from the predicted density (b). (bottom) Two LUMO isosurfaces at 0.15 and 0.01 a.u. associated with the reference DFT calculation (c) and those derived from the predicted density (d).

Figure 5 :

 5 Figure 5: (top) % RMSE of the electron density prediction carried out on a total of 10,000 test molecules as a function of the number of training molecules at increasing values of M . (bottom) Histogram of density absolute errors integrated in real space as a fraction of the total electronic charge computed at M = 10000 and N = 8000.

Figure 5

 5 Figure 5 reports learning curves for the prediction of the electron density using a maximum of 8000 randomly selected training structures reaching a small error of 1.22% for M = 10000.

Figure 6 :

 6 Figure 6: (top) % MAE of the total energy prediction as a function of the number of training molecules at a increasing values of M . (bottom) Histogram of absolute energy deviations computed at M = 10000 and N = 8000. The vertical dashed line at 1 kcal/mol indicates the chemical accuracy threshold.

  This leap forward was enough to successfully tackle the regression of disordered condensed-phase systems such as liquid water and to achieve chemical accuracy on most of the test structures of the QM9 dataset with a handful of training structures. Further work should be undertaken to reduce the computational burden required to keep in memory all the overlap matrices and feature vectors needed during the training phase, thus allowing us to increase the training set size indefinitely. This could be done either by adopting better sampling algorithms that allow the selection of smaller active sets M , which are fine-tuned for each basis function type, or by implementing more efficient minimization strategies that compute the gradient by loading the training data on-the-fly at each iteration.

The authors thank Michael J. Willatt for an in-depth discussion of the orthogonal learning framework. We thank Peter B. Jørgensen and Arghya Bhowmik for sharing the QM9 test configurations used as benchmark. A.G. acknowledges funding from the Swiss National Science Foundation. A.M.L. is supported by the Alexander von Humboldt Foundation. M.C. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 101001890-FIAMMA). A.L. and M.R. acknowledge computational time from the Max Planck Computing and Data Facility (MPCDF).

Source code and data availability

An open-source implementation of the method, together with the dataset of liquid water, can be downloaded at https://github.com/andreagrisafi/SALTED. A worked example of the workflow outlined in this work using FHI-aims with SALTED is also provided in this repository.

The reference calculations for both datasets are available at 10.17172/NOMAD/2022.10.17-1.

Supporting Information

Supporting Information includes: i) learning curves for the indirect prediction of kinetic, electrostatic and exchange-correlation energies of liquid water, ii) an analytical and numerical analysis of the error made in the indirect prediction of the total electronic energy.


Christophe Leclercq
email: christophe.leclercq@chu-rennes.fr

Haran Burri

Heart failure is a clinical syndrome resulting from structural and/or functional abnormalities of the heart that result in elevated intracardiac pressures and/or inadequate cardiac output. In patients with heart failure with reduced (<0.40) left ventricular ejection fraction (HFrEF), impaired left ventricular filling due to atrioventricular (AV) dyssynchrony resulting from prolonged PR intervals may contribute to pump failure. The mechanisms which have been suggested are transmitral EA fusion leading to shortening of left ventricular filling time, diastolic mitral regurgitation, and in case of very prolonged PR intervals, simultaneous atrial and ventricular contraction. [START_REF] Salden | Atrioventricular dromotropathy: evidence for a distinctive entity in heart failure with prolonged PR interval?[END_REF] Computer simulations have suggested that AV timing plays a crucial role in cardiac resynchronization therapy (CRT) by impacting ventricular filling. [START_REF] Jones | Cardiac resynchronization therapy: mechanisms of action and scope for further improvement in cardiac function[END_REF] In the COMPANION trial, a PR > 200 ms was observed in 52% of the subjects and was associated with increased mortality and heart failure hospitalization in patients under medical therapy (HR 1.41, P=0.44). Cardiac resynchronization therapy (CRT) had a greater impact on reducing these endpoints in patients with a prolonged PR (HR 0.54 P<0.01) than in those with a normal PR (HR 0.71, P=0.019). A subgroup analysis of MADIT-CRT showed that patients with non-left bundle branch block QRS morphology and a PR>230 ms benefitted from CRT in terms of heart failure and mortality, whereas those with intervals <230 ms did not. [START_REF] Kutyifa | PR interval identifies clinical response in patients with non-left bundle branch block: a Multicenter Automatic Defibrillator Implantation Trial-Cardiac Resynchronization Therapy substudy[END_REF] Patients with a narrow QRS have been shown in the EchoCRT trial to have increased mortality when treated with CRT, 4 so biventricular pacing in patients with a prolonged PR interval in this population is unlikely to be of benefit. Right ventricular pacing (RVP) is known to have a detrimental effect on left ventricular ejection fraction (LVEF) resulting from dyssynchrony. [START_REF] Cicchitti | Heart failure due to right ventricular apical pacing: the importance of flow patterns[END_REF][START_REF] Fang | Deterioration of left ventricular systolic function in extended Pacing to Avoid Cardiac Enlargement (PACE) trial: the predictive value of early systolic dyssynchrony[END_REF] His bundle pacing (HBP) on the other hand provides a more physiological form of pacing and provides an opportunity to correct AV dyssynchrony without inducing ventricular dyssynchrony.

In the current issue of the journal, Whinnett et al. [START_REF] Whinnet | Effects of haemodynamically atrio-ventricular optimized His-pacing on heart failure symptoms and exercise capacity: The His Optimized Pacing Evaluated for Heart Failure (HOPE-HF) randomized, double-blind, cross-over trial[END_REF] present the results of the HOPE-HF study, which evaluated whether HBP improves the primary endpoint of peak oxygen uptake in patients with HFrEF, PR interval > 200ms and either QRS <140ms or right bundle branch block. A total of 167 patients were randomized using a cross-over design to 6-months of DDD pacing (with non-invasive hemodynamic optimization of AV intervals) and 6-months of backup RVP, with double-blinding of treatment allocation. There was no significant difference in the primary endpoint, nor in the secondary endpoints of LVEF, left ventricular diameters, BNP or quality of life as measured by EQ5D5L VAS. However, there was a slight improvement in quality of life as measured by the MLWHF questionnaire and 112/167 (67%) of patients preferred the HBP pacing period.

Although the results are disappointing, the authors may be congratulated for their endeavour and for having conducted the largest randomized controlled study published to date. A number of points may be raised regarding the findings. 1) This was not a study randomizing HBP versus no pacing. During the HBP-off phase, patients were programmed to "backup" RVP at 30 bpm, presumably mostly in the DDD mode as atrial pacing was on average 16±28 % and ventricular pacing 13±28 %. It has been shown that as little as 20% RVP may result in pacing-induced cardiomyopathy, and that pacing percentage as a continuous variable is significantly associated with outcome (i.e. that there is no minimum percentage of RVP which has been demonstrated to be harmless). [START_REF] Kiehl | Incidence and predictors of right ventricular pacing-induced cardiomyopathy in patients with complete atrioventricular block and preserved left ventricular systolic function[END_REF][START_REF] Glikson | ESC Guidelines on cardiac pacing and cardiac resynchronization therapy[END_REF] Therefore, the only positive findings of the study may have resulted from unnecessary RVP, which may have been avoided by turning off pacing (ODO/OVO mode) or minimized by programming to DDI/VVI 30 bpm.

2) the average PR interval in the study was only 249 ms with a mean optimized AV delay of 195 ms. Therefore, as acknowledged by the authors, the reduction in AV was modest, and treatment effect may have been greater in a cohort with longer PR intervals. Indeed, patients who would have been most likely to benefit are those with extremely long PR intervals in whom a "P on T" phenomenon may be observed at rest or during exercise, leading to a "pseudo-pacemaker syndrome". In these patients, atrial and ventricular contraction occur simultaneously, thereby reducing ventricular filling, increasing atrial pressure and resulting in a backflow of blood. Very long PR intervals may also result from algorithms designed to minimize ventricular pacing. [START_REF] Barold | Conventional and biventricular pacing in patients with first-degree atrioventricular block[END_REF] RVP has been shown to improve symptoms in these patients, [START_REF] Carroz | Pseudo-pacemaker syndrome in a young woman with first-degree atrio-ventricular block[END_REF] with class IIa level of evidence C recommendation for pacing in symptomatic first degree AV block (PR >300 ms) according to current guidelines. [START_REF] Glikson | ESC Guidelines on cardiac pacing and cardiac resynchronization therapy[END_REF] In these patients with HFrEF and a narrow QRS, HBP is indeed likely to be superior to RVP and even to CRT, but there was no sub-analysis in HOPE-HF of efficacy stratified according to baseline PR interval.

3) The authors did not evaluate transmitral flow by echocardiography. This is important as patients with HFrEF may have a restrictive pattern, with low-velocity A waves in whom atrial filling contributes little to ventricular filling and are thus less likely to benefit from AV coupling. EA fusion resulting from PR prolongation is often labelled as being detrimental, and that prolongation of ventricular filling time with EA separation resulting from ventricular pacing is beneficial. However it is important to remember that it is the EA velocity-timeintegral and not EA duration (i.e. ventricular filling time) which determines ventricular filling. Also, EA fusion not only depends on the PR interval but also on heart rate and is physiological during tachycardia. Furthermore, diastolic mitral regurgitation could also have been evaluated by echocardiography, although it is usually only minimal in degree and the correction of this factor is likely to play a minor role in improving hemodynamics. [START_REF] Barold | Conventional and biventricular pacing in patients with first-degree atrioventricular block[END_REF] In an elegant acute hemodynamic study on 8 pigs and 22 CRT patients, Salden et al. [START_REF] Salden | Pacing therapy for atrioventricular dromotropathy: a combined computational-experimental-clinical study[END_REF] showed that at AV interval programming impacted left ventricular filling and cardiac output by increasing forward flow over the mitral valve (area under the E-and A-waves) and reducing diastolic mitral regurgitation, although how these factors react to physiological conditions and impact clinical outcome over follow-up is unknown. ) Interatrial conduction delay may mitigate the effect of a long PR interval by delaying left atrial conduction and thereby re-establish left AV synchrony. [START_REF] Barold | Conventional and biventricular pacing in patients with first-degree atrioventricular block[END_REF] The HOPE-HF study did not present any data on P-wave duration, and this important factor was not evaluated.

Finally, hemodynamic impact of varying AV delays is not always straightforward. In the BRAVO study, performed by the same group as HOPE-HF, as many as 60% of CRT patients had little impact on EA separation by decrementing AV intervals. [START_REF] Jones | Applicability of the iterative technique for cardiac resynchronization therapy optimization: full-disclosure, 50-sequential-patient dataset of transmitral Doppler traces, with implications for future research design and guidelines[END_REF] On a positive note, HOPE-HF has shown that HBP was successful in 93% of patients with a single procedure and with reasonable thresholds at randomization and 12-months (1.33V ± 1.07 and 1.47V ± 1.21 V @ 1 ms), and did not adversely affect outcome in patients with HFrEF. Lead dislodgment or significant rise in threshold compromising participation in the study was observed in 5.7% of patients, which is in line with previous data reporting a 6% rate of HBP lead revision, [START_REF] Beer | His bundle pacing capture threshold stability during long-term follow-up and correlation with lead slack[END_REF] which is nevertheless high. Left bundle branch area pacing (LBBAP) is gaining adoption and has overtaken HBP according to a recent survey, [START_REF] Kircanski | Conduction system pacing in everyday clinical practice: EHRA physician survey[END_REF] mainly due to superior electrical parameters. It is possible that conduction system pacing (CSP) will ultimately replace RVP in many indications, and be used in lieu of BiVP in selected cases, e.g. in HFrEF patients with a narrow QRS requiring ventricular pacing.

To conclude, HOPE-HF failed to find a benefit of pacing in HFrEF patients with first-degree AV block, but whether CSP may benefit these patients with more prolonged PR intervals (e.g. > 300 ms), deserves further study.

4

 4

  18790844, ja, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/ejhf.2775 by Universite de Rennes, Wiley Online Library on [20/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  18790844, ja, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/ejhf.2775 by Universite de Rennes, Wiley Online Library on [20/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  18790844, ja, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/ejhf.2775 by Universite de Rennes, Wiley Online Library on [20/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

has received speaker honoraria and institutional research/fellowship grants from Abbott, Biotronik, Boston Scientific, Medtronic and Microport. C.L. has received speaker honoraria and institutional research/fellowship grants from Abbott, Biotronik, Boston Scientific, Medtronic and