Robert Benda

Eric Cancès

Virginie Ehrlacher

Benjamin Stamm

Multi-center decomposition of molecular densities: a mathematical perspective

The aim of this paper is to analyze from a mathematical perspective some existing schemes to partition a molecular density into several atomic contributions, with a specific focus on Iterative Stockholder Atom (ISA) methods. We provide a unified mathematical framework to describe the latter family of methods and propose a new scheme, named L-ISA (for linear approximation of ISA). We prove several important mathematical properties of the ISA and L-ISA minimization problems and show that the so-called ISA algorithms can be viewed as alternating minimization schemes, which in turn enables us to obtain new convergence results for these numerical methods. Specific mathematical properties of the ISA decomposition for diatomic systems are also presented. We also review the basis-space oriented Distributed Multipole Analysis method, the mathematical formulation of which is also clarified. Different schemes are numerically compared on different molecules and we discuss the advantages and drawbacks of each approach.

Introduction

A successful study of intermolecular interactions is closely interwoven with the quality of the representation of the electronic density and subsequently of the electrostatic potential (ESP) around molecules. When a continuous representation of the density, such as in Density Functional Theory (DFT) [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF][START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF] or Wave-Function Theory (WFT) [START_REF] Jensen | Introduction to computational chemistry[END_REF], cannot be used because of a too high computational burden, the density has to be summarized by a finite number of parameters thanks to localization schemes. Classical force fields, such as GROMOS [START_REF] Walter | The GROMOS biomolecular simulation program package[END_REF], CHARMM [START_REF] Mackerell | All-atom empirical potential for molecular modeling and dynamics studies of proteins[END_REF] or AMBER [START_REF] Pearlman | AMBER, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules[END_REF][START_REF] Cornell | A second generation force field for the simulation of proteins, nucleic acids, and organic molecules[END_REF], among others, use such a representation, with (fixed) so-called partial charges assigned to all atomic sites, to represent the electronic density in a simplified way and account for electrostatic interactions with neighboring molecules. Molecular dynamics (MD) simulations based on these force fields are commonly used to gain insight into the structural or energetic properties of biologic or organic compounds. The physical relevance of MD studies heavily relies on their quality, and thus on the specific localization scheme used to compute them. The purpose of the present paper is precisely to review several localization schemes from a mathematical perspective, providing a sound mathematical basis to the latter. Let us first review the existing localization schemes that have been used to derive atomic partial charges or atomic multipole moments in the context of force field parameterization.

First, a popular method to derive partial charges is to fit the ESP generated by the set of atomic point charges to the exact quantum ESP (derived from DFT or WFT calculations) at a given geometry [START_REF] Singh | An approach to computing electrostatic charges for molecules[END_REF], possibly imposing additional constraints to lower the conformational dependence of the obtained optimal partial charges, such as in the RESP method [START_REF] Christopher | A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model[END_REF]. Adding constraints to reproduce the total dipole moments at different conformations is also possible [START_REF] Reynolds | Atomic charges for variable molecular conformations[END_REF]. So is to fit the ESP generated by the atomic charges to the approximate ESP generated by a series of distributed multipole moments [START_REF] Stone | Distributed multipole analysis, or how to describe a molecular charge distribution[END_REF][START_REF] Stone | Distributed multipole analysis[END_REF], as implemented e.g. in the Mulfit code [START_REF] György | Charges derived from distributed multipole series[END_REF]. In both cases, the partial charges are chosen as those leading to the smallest error on the ESP for a suitable norm. This family of ESP-derived charges methods amounts to find indirectly the partial charges from the quantum calculation, using the electrostatic potential as a proxy -and possibly the molecular dipole or quadrupole moments -as relevant target. However, they are associated to several drawbacks, such as conformational-dependent charges -although strategies to mitigate this undesired feature have been designed [START_REF] Reynolds | Atomic charges for variable molecular conformations[END_REF] -and lack of transferability. Chipot et al. have for instance shown that ESP derived partial charges are not transferable from one to another similar molecule [START_REF] Chipot | Transferable net atomic charges from a distributed multipole analysis for the description of electrostatic properties: a case study of saturated hydrocarbons[END_REF]. In other words, the partial charge assigned to a carbon atom in a given specific local chemical environment of a given molecule cannot be used as partial charge for a carbon atom in a similar local chemical environment of another molecule. It is also possible to derive atomic partial charges directly from a quantum calculation (namely, using the coefficients of the converged density matrix) using e.g. Mulliken [START_REF] Mulliken | Electronic population analysis on LCAO-MO molecular wave functions. i[END_REF] or Löwdin [START_REF] Löwdin | Approximate formulas for many-center integrals in the theory of molecules and crystals[END_REF] population analysis, although these methods are not commonly used for partial charges parameterizationor, if they are, charges are readjusted afterwards by fitting to the ESP, to minimize the deviation between the reference QM ESP and the atomic charges-generated ESP. Indeed, Mulliken and Löwdin schemes adopt a basis space approach, which results in significant basis-set dependency of partial charges and lack of convergence with increasing basis set size.

Real-space approaches, using a Voronoi-like partitioning of the space, such as the Bader Atom In Molecules [START_REF] Richard | Atoms in Molecules. A Quantum Theory[END_REF], the Becke space-decomposition (fuzzy Voronoi cells) approach [START_REF] Becke | A multicenter numerical integration scheme for polyatomic molecules[END_REF], or stockholder partitionning schemes originating from the Hirshfeld method [START_REF] Hirshfeld | XVII. Spatial partitioning of charge density[END_REF] (see below), have also been used to derive partial charges, and do not suffer as much from basis set dependencebut may suffer other drawbacks [START_REF] Heidar-Zadeh | Information-theoretic approaches to atomsin-molecules: Hirshfeld family of partitioning schemes[END_REF].

Finally, thanks to the recent availability of large databases of molecules already parameterized for classical force fields, modern machine-learning inspired methods have also been used to estimate partial charges, e.g. formulating the charge assignment problem as a multiple-choice knapsack problem [START_REF] Engler | Multiple-Choice Knapsack for Assigning Partial Atomic Charges in Drug-Like Molecules[END_REF]. Yet, although these methods allow to avoid explicit charge calculations, they may end up with atomic charges suffering from the same drawbacks as those of the database on which they have been fitted.

To correct for the too approximate nature of the ESP represented by atomic partial charges only [START_REF] Chipot | Transferable net atomic charges from a distributed multipole analysis for the description of electrostatic properties: a case study of saturated hydrocarbons[END_REF] (among other drawbacks), higher order moments of the charge distribution have been introduced within localization schemes, and consequently in force fields used for MD simulations. At their early times, multipole moments of the electronic density were introduced as a generalization of Mulliken population analysis, to analyse and interpret molecular wave-functions with a larger number of local descriptors [START_REF] Stone | Distributed multipole analysis, or how to describe a molecular charge distribution[END_REF][START_REF] Stone | Distributed multipole analysis[END_REF][START_REF] Vigné | The exact multicenter multipolar part of a molecular charge distribution and its simplified representations[END_REF]. Since the advent of modern polarizable force fields, such as AMOEBA [START_REF] Ren | Consistent treatment of inter-and intramolecular polarization in molecular mechanics calculations[END_REF][START_REF] Ren | Polarizable atomic multipole water model for molecular mechanics simulation[END_REF][START_REF] Pengyu Ren | Polarizable atomic multipole-based molecular mechanics for organic molecules[END_REF][START_REF] Wu | Automation of AMOEBA polarizable force field parameterization for small molecules[END_REF], these local descriptors have also been used successfully to improve the description of biomolecular and organic molecules, and of their mutual interactions, in particular thanks to an improved treatment of the water solvent. Local multipole moments -which are used e.g. up to quadrupole moments in AMOEBA force field -allow to capture the local anisotropy of the charge density (contrary to atomic charges alone) and are thus a key ingredient of polarizable force fields.

In a similar manner as for partial charges, local multipole moments can be derived either indirectly (by fitting the multipoles-derived ESP to the quantum ESP) or directly (e.g. using the density matrix corresponding to a certain basis expansion, or the values of the density on a grid) from the result of a quantum calculation. To the former rationale belongs the fitting method implemented e.g. in the Mulfit code [START_REF] György | Toward improved force fields. 2. Effective distributed multipoles[END_REF][START_REF] Winn | Toward improved force fields. 1. multipole-derived atomic charges[END_REF], which allows to derive atomic charges [START_REF] Winn | Toward improved force fields. 1. multipole-derived atomic charges[END_REF] or multipoles [START_REF] György | Toward improved force fields. 2. Effective distributed multipoles[END_REF] (up to user-specified orders, for every atom) by fitting on the ESP generated by higher order local multipoles -typically obtained previously from distributed multipole analysis (DMA) [START_REF] Stone | Distributed multipole analysis, or how to describe a molecular charge distribution[END_REF][START_REF] Stone | Distributed multipole analysis[END_REF]. To the latter family of methods belong Stone's DMA [START_REF] Stone | Distributed multipole analysis, or how to describe a molecular charge distribution[END_REF][START_REF] Stone | Distributed multipole analysis[END_REF], performing a redistribution of the density matrix coefficients and elementary local multipole moments to final expansion sites (e.g. all the atoms), and real-space partitioning schemes in the continuity of Bader concepts of atoms in molecules [START_REF] Bader | Virial field relationship for molecular charge distributions and the spatial partitioning of molecular properties[END_REF]. The family of Hirshfeld-like partitionning schemes [START_REF] Hirshfeld | Bonded-atom fragments for describing molecular charge densities[END_REF][START_REF] En Maslen | Atomic charges and electron density partitioning[END_REF][START_REF] Ayers | Atoms in molecules, an axiomatic approach. i. maximum transferability[END_REF][START_REF] Robert G Parr | What is an atom in a molecule?[END_REF], namely its iterative version Hirshfeld-I [START_REF] Bultinck | Uniqueness and basis set dependence of iterative hirshfeld charges[END_REF][START_REF] Bultinck | Critical analysis and extension of the hirshfeld atoms in molecules[END_REF], the iterative stockholder approach (ISA) [START_REF] Timothy | Redefining the atom: atomic charge densities produced by an iterative stockholder approach[END_REF][START_REF] Lillestolen | Atomic charge densities generated using an iterative stockholder procedure[END_REF] and its variants GISA [START_REF] Verstraelen | The conformational sensitivity of iterative stockholder partitioning schemes[END_REF], MB-ISA [START_REF] Verstraelen | Minimal basis iterative stockholder: atoms in molecules for force-field development[END_REF], and more recent basis-space implementation (BS-ISA) [START_REF] Misquitta | Distributed multipoles from a robust basis-space implementation of the iterated stockholder atoms procedure[END_REF], belong to this family of real-space oriented methods. This family of methods yields multipole moments which have been claimed to be rather transferable [START_REF] Heidar-Zadeh | Information-theoretic approaches to atomsin-molecules: Hirshfeld family of partitioning schemes[END_REF]. Note that in principle, real-space partitioning schemes end up with local multipole moments located at atomic sites only (to enable force computations in force field applications), although this is not a requirement. The DMA method can also naturally allocate contributions to atoms and to other non-atomic sites, such as bond centers.

In this article, we review and clarify several existing schemes to partition a molecular density into several (e.g. atomic) contributions or to derive directly distributed multipole moments. In Section 2, we first focus on the real-space oriented methods of the Hirshfeld / Iterative Stockholder Atom (ISA) family. We provide a unified mathematical framework to describe the existing ISA methods, through the formulation [START_REF] Goldfarb | A numerically stable dual method for solving strictly convex quadratic programs[END_REF]. This formulation is based on information theory, a link that was first introduced by Nalewajski and Parr (see [START_REF] Roman | Information theory, atoms in molecules, and molecular similarity[END_REF]) and then subsequently used in the works of Ayers, Bultinck, Heidar-Zadeh and co-workers (see [START_REF] Bultinck | Uniqueness and basis set dependence of iterative hirshfeld charges[END_REF][START_REF] Bultinck | Critical analysis and extension of the hirshfeld atoms in molecules[END_REF][START_REF] Heidar-Zadeh | Information-theoretic approaches to atomsin-molecules: Hirshfeld family of partitioning schemes[END_REF] and references therein), but is new to our knowledge. It is also worthwhile to stress the mathematical properties established in [START_REF] Ayers | Atoms in molecules, an axiomatic approach. i. maximum transferability[END_REF]. Driven by formulation [START_REF] Goldfarb | A numerically stable dual method for solving strictly convex quadratic programs[END_REF], we propose a new AIM scheme, named L-ISA (linear approximation of ISA), and presented in Section 2.5.5. We also formulate the Minimal-Basis ISA (MB-ISA) method as a minimization of the Kullback-Leibler entropy between atomic-shell densities and the corresponding pro-atomic-shell densities, under equality constraints involving zero and first order moments of the atomic densities (see section 2.7). We then focus in Section 3 on the mathematical analysis of the ISA and L-ISA methods. We first prove the existence and uniqueness of the solutions to the ISA and L-ISA optimization problems (Theorem 6). While uniqueness easily follows from the strict convexity properties of the Kullback-Leibler entropy, establishing the existence of a minimizer is more difficult due to possible loss of compactness. We then prove that the ISA and L-ISA algorithms can be interpreted as alternating minimization methods, that the associated entropy is a Lyapunov function of the algorithm, and that the L 2 norm between two successive iterates converges to 0. We also prove that, in the case of the L-ISA method, the L-ISA algorithm converges toward the unique minimizer of the L-ISA optimization problem. All these results are collected in Theorem 7. We finally prove specific mathematical properties of ISA decomposition for diatomic systems in Propositions 8 and 9. In Section 4, we review a popular basis-space oriented method, the Distributed Multipole Analysis (DMA), whose mathematical formulation is also clarified. Finally, we present in Section 5 some numerical results (with a focus on diatomic systems) of the different ISA schemes, namely GISA, L-ISA (section 2.5.5), MB-ISA and (historical) ISA (section 2.5.4), which have been implemented numerically. The DMA method, which has also been implemented in a more modern and modular form, enabling for more redistribution strategies to be explored (section 4), is also compared to the ISA method in terms of convergence of local multipole moments with increasing basis set size, for some test systems.

AIM methods based on objective functionals

We denote by

X := f ∈ L 1 (R 3 ) ∩ L ∞ (R 3 ) lim |r|→∞ f (r) = 0, ˆR3 |r||f (r)| dr < ∞
the Banach space of bounded, integrable, real-valued functions on the physical space R 3 with finite first moments, and vanishing at infinity. In the sequel, we consider a given nonnegative density ρ ∈ X such that

N := ˆR3 ρ(r) dr > 0,
and a collection R = (R a ) 1≤a≤M ∈ (R 3 ) M of M points in R 3 . Ground and excited state electronic densities of molecules and clusters are known to be continuous, positive, functions on R 3 decaying exponentially fast at infinity (see e.g. [START_REF] Fournais | The electron density is smooth away from the nuclei[END_REF]Theorem 1.3]), and are therefore elements of X.

Set of admissible AIM decompositions

The goal is to decompose ρ as a sum of translated nonnegative densities ρ opt a ∈ X, i.e.

ρ(r) = M a=1 ρ opt a (r -R a ), (1)
each ρ opt a being localized around the origin and optimized in some sense. In quantum chemistry, the R a 's are most often the positions of the M nuclei of a molecular system containing N ∈ N * electrons, and ρ an approximation of its ground-state (or k-th excited-state) electronic density obtained by a given electronic structure calculation method in a given basis set. We will then denote by z a ∈ N * the charge of nucleus a for 1 ≤ a ≤ M . However, in some applications, it can be useful to choose expansion centers away from the nuclei, for instance on a chemical bond, or at the center of mass of a functional group. We will therefore call the points (R a ) 1≤a≤M the expansion centers (rather than the atomic positions); for convenience, we will however use the standard terminology atoms-in-molecules (AIM) to refer to the family of functions (ρ opt a ) 1≤a≤M . For later purposes, we introduce the atomic charge linear map N :

X M → R M defined by ∀ρ := (ρ a ) 1≤a≤M ∈ X M , [N (ρ)] a = ˆR3 ρ a (r) dr, ∀1 ≤ a ≤ M. (2)
Remark 1. We use the decomposition (1) instead of the (equivalent) usual decomposition ρ(r) = M a=1 ρ opt a (r) because the former is more convenient for the analysis of the dependence of the AIM decomposition on the atomic positions. For instance, it is expected that the ρ opt a (centered at the origin) in (1) converge to the ground (or excited) state density of the isolated atom or ion in the dissociation limit.

The ρ opt

a 's can be used as such, or transformed into a collection of descriptors, typically their charges, and their first (dipolar) and second-order moments:

q a := ˆR3 ρ opt a (r) dr, p a := ˆR3 rρ opt a (r) dr, Q a := ˆR3 r ⊗ r ρ opt a (r) dr.
We denote by X + := {f ∈ X | f ≥ 0 a.e.} the convex cone of bounded integrable nonnegative densities vanishing at infinity and introduce

K ρ,R := ρ = (ρ a ) 1≤a≤M ∈ X M + M a=1 ρ a ( • -R a ) = ρ , (3)
the set of admissible AIM decompositions of the density ρ with respect to the expansion cen-

ters R = (R a ) 1≤a≤M ∈ (R 3 ) M .
The set K ρ,R is non-empty and convex. It has also some nice topological properties which will allow us to prove the existence of an optimal decomposition for the usual AIM decomposition methods (see the mathematical results in Section 3 and their proofs in Section 6).

Main ingredients of an AIM decomposition method

Let us denote by X r ⊂ X (respectively X r + ⊂ X + ) the subset of functions of X (respectively X + ) that are radially symmetric.

An objective-functional-based AIM decomposition method consists in seeking an optimal decomposition ρ opt := (ρ opt a ) 1≤a≤M ∈ K ρ,R , for a specific criterion. It is characterized by two main ingredients:

(i) a collection of sets of (radially symmetric) proatom densities (K 0 z ) z∈N * , so that for each z ∈ N * , K 0 z ⊂ X r + contains reasonable approximations of the ground state density of an isolated atom with nuclear charge z. The set of proatom densities is then used to build a set of admissible proatoms-in-molecule (or promolecule) densities

K 0 := K 0 z 1 × • • • × K 0 z M = ρ 0 := (ρ 0 a ) 1≤a≤M , ρ 0 a ∈ K 0 za ⊂ (X r + ) M (4)
(ii) an objective functional E : X M + × (X r + ) M → R ∪ {+∞}, usually in the form of a relative entropy.

The set K 0 z may, or not, depend explicitly on z. It may have the structure of a vector space but does not need to in general. For all 1 ≤ a ≤ M , it is expected that, ρ opt a should be close in some sense to one of the elements of K 0 za , but let us emphasize that, in general, ρ opt a will not be a radially symmetric function, hence will not belong to K 0 za . The quantity E(ρ, ρ 0 ) measures the discrepancy between a trial AIM decomposition ρ := (ρ a ) 1≤a≤M ∈ X M + and a trial promolecule density ρ 0 := (ρ 0 a ) 1≤a≤M ∈ K 0 . Most objectivefunctional-based AIM decomposition methods (including ISA, Hirschfeld, Hirschfeld-I, L-ISA, NL-ISA) consist in finding ρ opt , ρ 0,opt ∈ D ρ,R solution to the minimization problem

ρ opt , ρ 0,opt ∈ argmin (ρ,ρ 0 )∈D ρ,R E(ρ, ρ 0 ) (5)
where D ρ,R is a well-chosen subset of K ρ,R × K 0 .

In Sections 2.3 and 2.4, we detail various possible choices of sets of proatom densities K 0 z and objective functionals which are encountered in the literature. The aim of Section 2.5 is to present in more details the methods based on the Kullback-Leibler divergence, including Hirschfeld, Hirschfeld-I, ISA, L-ISA and NL-ISA. The GISA and MB-ISA methods cannot be expressed under the form (5); their mathematical structures are presented in Sections 2.6 and 2.7 respectively.

Proatom densities

As anticipated above, we now discuss several choices for the proatom densities K 0 z .

1 In the original Hirschfeld method [START_REF] Hirshfeld | Bonded-atom fragments for describing molecular charge densities[END_REF], the set K 0 z is a singleton and only contains the ground state density of the neutral atom K 0 z,Hirschfeld = ρ 0 z,z , where ρ 0 z,n denotes (an approximation of) the radially symmetric ground state density of the atomic system consisting of a single nucleus of charge z and n electrons 1 , being understood that ρ 0 z,n = ρ 0 z,N max z for n greater than N max z , the maximum number of electrons that a nucleus of charge z can bind. However, choosing for ρ 0 a electronic densities of neutral atoms does not seem appropriate for strongly polarized molecules where significant charge transfers occur. This and other drawbacks are detailed in [START_REF] Bultinck | Critical analysis and extension of the hirshfeld atoms in molecules[END_REF], such as the fact that the Hirshfeld method is limited to neutral molecules, and that Hirshfeld atomic charges are on average too small to accurately describe molecular polarization. A desirable feature of an AIM method is its ability to detect charge transfer in an automatic way, without requiring a priori knowledge from the user on the chemical system under consideration.

2 In the so-called Hirschfeld-I method [START_REF] Bultinck | Critical analysis and extension of the hirshfeld atoms in molecules[END_REF] (where I stands for iterative), this is achieved to some point by the following choice of proatom densities which contains not only the neutral groundstate atomic density, but also the ground-state densities of the ionized forms of the atom, as well as those of the mixed states corresponding to fractional (in the sense of non-integer) numbers of electrons:

K 0 z,Hirschfeld-I = ρ 0 z,n , n ∈ R + with ∀n ∈ R + , ρ 0 z,n (r) := ( n -n) ρ 0 z, n (r) + (n -n ) ρ 0 z, n (r)
, where for x ∈ R, x and x respectively denote the largest integer smaller than x, and the smallest integer larger than x.

3 In the iterative stockholder approach (ISA) [START_REF] Timothy | Redefining the atom: atomic charge densities produced by an iterative stockholder approach[END_REF][START_REF] Lillestolen | Atomic charge densities generated using an iterative stockholder procedure[END_REF], this set is independent of the chemical element and is chosen as the whole convex cone of bounded integrable nonnegative radially symmetric functions vanishing at infinity:

K 0 ISA = X r + .
4 In the finite-dimensional linear approximations of ISA (L-ISA), the set K 0 z is chosen to be a non-empty closed convex subset of X r + . More precisely,

K 0 z,L-ISA = ρ 0 (r) = mz k=1 c k g z,k (r), c k ∈ R + , (6)
where the g z,k 's, 1 ≤ k ≤ m z are z-dependent given linearly-independent, positive, L 1 -normalized functions of K 0 ISA = X r + .

5 In the finite-dimensional nonlinear approximations of ISA (NL-ISA),

K 0 z = K 0 z,NL-ISA = ρ 0 (r) = mz k=1 c k g z,k,α k (r), c k ∈ R + , α k > 0 , (7)
where the g z,k,α k 's are L 1 -normalized positive functions of K 0 ISA = X r + depending on a parameter α k in a non-affine manner and which needs to be optimized. For the sake of simplicity, here, we assume that the parameters α k are positive real numbers, but α k could be a vector subjected to equality and/or inequality constraints.

6 In GISA [START_REF] Verstraelen | The conformational sensitivity of iterative stockholder partitioning schemes[END_REF], the set K 0 z is defined as the finite-dimensional convex cone

K 0 z,GISA = ρ 0 z (r) = mz k=1 c k ζ α z,k (r), c k ∈ R + where ζ α (r) = α π 3/2 e -α|r| 2 , (8)
and where m z ∈ N * and α z,k > 0 are fixed empirical parameters.

Objective functionals

As mentioned above, the aim of the quantity E(ρ, ρ 0 ) is to measure the discrepancy (in some sense) between a trial AIM decomposition ρ := (ρ a ) 1≤a≤M ∈ X M + and a trial promolecule density ρ 0 := (ρ 0 a ) 1≤a≤M ∈ K 0 . Most of the popular information-theory-based AIM methods (Hirschfeld, Hirschfeld-I, ISA, L-ISA, NL-ISA -see below for further details) make use of the relative entropy

E(ρ, ρ 0 ) := S(ρ|ρ 0 ), (9)
where

S(ρ|ρ 0 ) := M a=1 ˆR3 ρ a (r) log ρ a (r) ρ 0 a (r) dr = M a=1 s KL (ρ a |ρ 0 a ), (10)
constructed from the Kullback-Leibler divergence commonly used in information theory:

∀f, g ∈ X + , s KL (f |g) := ˆR3 f (r) log f (r) g(r) dr. (11)
In order to make the above definition for L 1 functions consistent with the general definition of the Kullback-Leibler divergence for bounded positive measures, the following conventions must be used: 0 log 0 0 = 0 and for all p > 0, 0 log 0 p = 0, p log p 0 = +∞.

Note however that alternative methods have been proposed [START_REF] Heidar-Zadeh | Deriving the hirshfeld partitioning using distance metrics[END_REF], where the objective functional E(ρ, ρ 0 ) is chosen as

E(ρ, ρ 0 ) = M a=1 s(ρ a |ρ 0 a ),
where s(•|•) is either the Hellinger distance

s Hel (f |g) := ˆR3 f (r) -g(r) 2 dr = ˆR3 f (r) g(r) f (r) -1 2 , (13)
or a more general divergence of the form

s φ (f |g) := ˆR3 f (r)φ g(r) f (r) , (14)
where φ : R + → R ∩ {+∞} is a strictly convex function, smooth on R * + and such that φ(1) = 0. The GISA method [START_REF] Verstraelen | The conformational sensitivity of iterative stockholder partitioning schemes[END_REF] makes use of both the relative entropy S(ρ|ρ 0 ) and the L 2 -distance

d L 2 (ρ, ρ 0 ) := M a=1 ρ a -ρ 0 a 2 L 2 1/2
, in a more complex fashion as will be summarized below in Section 2.6.

Lastly, the Minimal Basis Iterative Stockholder algorithm (MB-ISA) [START_REF] Verstraelen | Minimal basis iterative stockholder: atoms in molecules for force-field development[END_REF] is in fact not based on a direct atomic decomposition on ρ. It rather provides an atomic-shell decomposition of ρ on proatomic-shell densities, from which AIM densities can be derived. As will be seen in Section 2.7, it nevertheless nicely fits in the unified mathematical framework described in this article.

Kullback-Leibler entropy methods

In this section, we more specifically focus on the AIM decomposition methods the objective functional of which is based on the Kullback-Leibler divergence S(•|•). Let us first recall the well-known formula (16) below, a rigorous proof of which is provided in Section 6.

Lemma 2. Let ρ 0 := (ρ 0 a ) 1≤a≤M ∈ (X r + ) M and ρ ∈ X + be such that

s KL (ρ|ρ 0 ) < +∞ where ρ 0 (r) := M a=1 ρ 0 a (r -R a ). (15)
Then, the constrained minimization problem

inf ρ∈K ρ,R S(ρ|ρ 0 ) admits an unique minimizer G ρ,R (ρ 0 ) := ([G ρ,R (ρ 0 )] a ) 1≤a≤M , given by ∀1 ≤ a ≤ M, [G ρ,R (ρ 0 )] a (r) = ρ 0 a (r) ρ 0 (r+Ra) ρ(r + R a ) if ρ 0 a (r) > 0, 0 if ρ 0 a (r) = 0, (16)
and it holds that min

ρ∈K ρ,R S(ρ|ρ 0 ) = s KL (ρ|ρ 0 ). (17)
Observe that if ρ 0 a (r) > 0, then ρ 0 (r + R a ) ≥ ρ 0 a (r) > 0, while if ρ 0 a (r) = 0, then we can have either ρ 0 (r + R a ) > ρ 0 a (r), in which case ρ 0 a (r) ρ 0 (r+Ra) ρ(r + R a ) is well-defined and equal to 0, or ρ 0 (r + R a ) = ρ 0 a (r) = 0, in which case ρ 0 a (r) ρ 0 (r+Ra) ρ(r + R a ) is not well-defined.

Remark 3. In the case when E(ρ, ρ 0 ) = S(ρ|ρ 0 ) and D ρ,R = K ρ,R ×K 0 , we infer from Lemma 2 that the minimization problem in (5) is equivalent to

inf ρ 0 :=(ρ 0 a ) 1≤a≤M ∈K 0 s KL ρ M a=1 ρ 0 a (r -R a ) . ( 18
)
Although this problem seems natural, it is not the one solved in most AIM methods, probably because it is computationally challenging in general. Only the basic Hirschfeld method fits into this framework, as will be seen in the next section, but in this case, (18) is trivial since K 0 Hirschfeld = {ρ 0 Hirschfeld } contains only one element.

The Hirschfeld method

The Hirshfeld method consists in solving

inf ρ∈K ρ,R S(ρ|ρ 0 Hirschfeld ), (19)
where ρ 0 Hirschfeld = (ρ 0 za,za ) 1≤a≤M is given. As already mentioned in Remark 3, this problem fits into the general framework (5) using the objective functional E defined in [START_REF] Benda | Dma[END_REF] and D ρ,R = K ρ,R × K 0 Hirschfeld . The atomic densities ρ 0 za,za are positive everywhere and decay as exponentials at infinity; in addition, exact molecular densities and their usual approximations decay (at least) exponentially fast at infinity. The assumptions in Lemma 2 are thus satisfied, which implies that the basic Hirschfeld problem [START_REF] Frisch | Gaussian 16 Revision C.01[END_REF] is well-posed and has an explicit solution:

ρ opt,Hirschfeld a (r) = [G ρ,R (ρ 0 Hirschfeld )] a (r) = ρ 0 za,za (r) N b=1 ρ 0 z b ,z b (r -R b + R a ) ρ(R a + r). (20)

AIM iterative methods based on the Kullback-Leibler divergence

The ISA, Hirschfeld-I, L-ISA and NL-ISA methods can be formulated using the following unified formalism: the optimal AIM decomposition is obtained by solving

inf (ρ,ρ 0 )∈C ρ,R S(ρ|ρ 0 ) (21)
where

C ρ,R = (ρ, ρ 0 ) ∈ K ρ,R × K 0 | N (ρ) = N (ρ 0 )
with N being the atomic charge map defined in [START_REF] Andersen | CVXOPT: A Python package for convex optimization[END_REF]. Note that the variation that determines the different methods (ISA, Hirschfeld-I, L-ISA, NL-ISA) lies in the definition of K 0 defined from the

K 0 z 's given in Section 2.3, i.e. K 0 z = K 0 z,Hirschfeld-I , K 0 z = K 0 ISA , K 0 z = K 0 z,L-ISA or K 0 z = K 0 z,NL-ISA .
These methods hence fit into the general framework (5) using the objective functional E defined in [START_REF] Benda | Dma[END_REF] and D ρ,R = C ρ,R with corresponding K 0 . The solution to [START_REF] Goldfarb | A numerically stable dual method for solving strictly convex quadratic programs[END_REF] is not explicit and must be computed numerically, by an iterative algorithm. The constraints N (ρ) = N (ρ 0 ) request that for all 1 ≤ a ≤ M , ρ a and ρ 0 a have the same charge, i.e. ´R3 ρ a (r) dr = ´R3 ρ 0 a (r) dr, which is not the case for the optimal AIM obtained with the basic (non-iterative) Hirschfeld method.

Remark 4. For this approach to make sense, the least one can ask is that there exists (ρ, ρ 0 ) ∈ C ρ,R such that S(ρ|ρ 0 ) < +∞.

(

) 22
This condition guarantees that the infimum in (21) is finite. Note that it is not satisfied, even for physical ground-state electronic densities ρ, if e.g. the K 0 z 's only contain compactly supported functions. For [START_REF] Heidar-Zadeh | Deriving the hirshfeld partitioning using distance metrics[END_REF] to be satisfied, the sets K 0 z must contain functions which do not decay extremely fast at infinity. This is always the case for ISA for which the sets K 0 ISA = X r + ' contains the function e -|r| , so that taking ρ a (r) = M -1 ρ(r + R a ) and ρ 0 a (r) = N 8πM e -|r| , we have

S(ρ|ρ 0 ) = ˆR3 ρ log ρ - N M log M -N log N 8πM + 1 M M a=1 ˆR3 |r|ρ(r + R a ) dr,
which is a finite quantity for all ρ ∈ X + . In the Hirschfeld-I, L-ISA, or NL-ISA settings, the sets K 0 za must contain functions decaying asymptotically not faster than e -α|r| (for some α > 0) for [START_REF] Heidar-Zadeh | Deriving the hirshfeld partitioning using distance metrics[END_REF] to be satisfied for all ρ ∈ X + . The requirements on the sets K 0 za can be weakened by imposing more conditions on ρ. For instance, if we have ´R3 |r| 2 ρ(r) dr < ∞, then condition [START_REF] Heidar-Zadeh | Deriving the hirshfeld partitioning using distance metrics[END_REF] is satisfied as soon as the sets K 0 za contain functions decaying asymptotically not faster than a Gaussian. Most of the iterative numerical schemes for solving [START_REF] Goldfarb | A numerically stable dual method for solving strictly convex quadratic programs[END_REF] that have been proposed in the literature amounts to minimizing S(ρ|ρ 0 ) alternatively with respect to each of the variables ρ 0 and ρ, under suitable constraints. They read as follows.

Generic iterative AIM algorithm:

• Initialization: Choose ρ 0,(0) ∈ K 0 such that s KL (ρ|ρ 0,(0) ) < +∞.

• Iteration m ≥ 0:

Step 1: set

ρ (m) = G ρ,R (ρ 0,(m-1) ) = arg min ρ∈K ρ,R S(ρ|ρ 0,(m-1) ), (23)
Step 2: find ρ 0,(m) ∈ arg min

ρ 0 ∈K 0 , N (ρ 0 )=N (ρ (m) ) S(ρ (m) |ρ 0 ). (24)
Recall that the solution to [START_REF] Heidar-Zadeh | Information-theoretic approaches to atomsin-molecules: Hirshfeld family of partitioning schemes[END_REF] in Step 1 is explicit, since the expression of the function G ρ,R (ρ 0 ) is given by [START_REF] György | Toward improved force fields. 2. Effective distributed multipoles[END_REF].

Due to the particular nature of S(•|•), the solution to Step 2 can then be obtained by solving M independent and local problems of the form

inf ρ 0 a ∈K 0 za , ´R3 ρ 0 a = ´R3 ρa s KL (ρ a |ρ 0 a ) (25)
for

ρ a = ρ (m-1) a .
At the noticeable exception of the ISA method, solving the optimization problem in Step 2 for the various AIM methods presented above require evaluating the integrals

N (m) a := [N (ρ (m) )] a = ˆR3 ρ (m) a (r) dr,
with ρ (m) = G ρ,R (ρ 0,(m-1) ).

Step 2 therefore requires an efficient numerical quadrature scheme to evaluate integrals of the form

Q a,ρ,R (ρ 0 ) := [N (G ρ,R (ρ 0 ))] a = ˆR3 ρ 0 a (r) N b=1 ρ 0 b (r -R b + R a ) ρ(r + R a ) dr, (26)
(with the convention that the integrand is equal to zero when the denominator vanishes), for ρ 0 = (ρ 0 a ) 1≤a≤M in the set K 0 of promolecule densities. The elements of the sets K 0 GISA , K 0 Hirschfeld-I or K 0 MB-ISA of promolecule densities constructed from the GISA, Hirschfeld-I, and MB-ISA proatom densities respectively, are smooth, fast-decaying, radially symmetric functions. The integrals of the form [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF] can therefore be efficiently computed using a one-dimensional quadrature scheme for the radial part and a Lebedev quadrature scheme for the angular part.

We successively present in details in the following sections how problems of the form [START_REF] Hirshfeld | XVII. Spatial partitioning of charge density[END_REF] are solved for four different methods: the Hirschfeld-I method, the general ISA method, the finitedimensional linear approximations of the ISA (L-ISA) method, the finite-dimensional nonlinear approximations of the ISA (NL-ISA) method.

Note that the abstract versions L-ISA and NL-ISA haven not been reported in the literature, although they are similar, but not identical, to the GISA and MB-ISA methods respectively. For L-ISA however, we can establish rigorous results (see Section 3) such as existence and uniqueness of the solution as well as global convergence of the AIM iteration scheme.

Hirschfeld-I iterations

For the Hirschfeld-I method [START_REF] Bultinck | Critical analysis and extension of the hirshfeld atoms in molecules[END_REF], solving (25) is trivial since the set {ρ 0 a ∈ K 0 za,Hirschfeld-I | ´R3 ρ 0 a = N a } contains a single element, namely ρ 0 za,Na . At the continuous level (i.e., without any discretization of the proatom densities), the Hirschfeld-I iterative procedure corresponding to the generic iterative AIM algorithm for

K 0 = K 0 Hirschfeld-I = K 0 z 1 ,Hirschfeld-I × • • • × K 0 z M ,Hirschfeld-I
can therefore be formulated more explicitly as follows:

Hirschfeld-I algorithm:

• Initialization: choose an initial guess of charges

N (0) = (N (0) a ) 1≤a≤M ∈ (R + ) M so that N a=1 N (0) a = N . • Iteration m ≥ 0: set ρ (m) a (r) := ρ 0 za,N (m) a (r) M b=1 ρ 0 z b ,N (m) b (r -R b + R a ) ρ(R a + r), (27)
where

ρ 0 za,N (m) a = N (m) a -N (m) a ρ 0 za, N (m) a + N (m) a -N (m) a ρ 0 za, N (m) a
; compute

N (m+1) a := ˆR3 ρ (m) a (r) dr. (28)
The atomic charges N (m) a at iteration m are fractional (i.e. non-integer) in general. The global convergence property of this iterative scheme was observed numerically in [START_REF] Bultinck | Critical analysis and extension of the hirshfeld atoms in molecules[END_REF] on a benchmark of 168 molecules. A convergence proof was proposed in [START_REF] Bultinck | Uniqueness and basis set dependence of iterative hirshfeld charges[END_REF].

From a numerical point of view, the Hirschfeld-I method can be seen as a fixed-point procedure of the form N (m+1) = F Hirschfeld-I (N (m) ), where F Hirschfeld-I : R M + → R M + is the function defined by ( 27)- [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF]. Using N = (N a ) 1≤a≤M ∈ R M + as main variables, it is not necessary to store in memory the functions ρ (m) a : they can be evaluated using ( 27) from ρ, N (m-1) , and the ρ 0 z,n 's. The implementation of the Hirschfeld-I method thus requires:

1. the pre-computation and storage of the radially symmetric functions ρ 0 z,n , for all chemical elements contained in the molecular system of interest, and integer values n ≤ N max z . Note that these functions can be pre-computed once and for all with high accuracy and then be used for any molecular system; 2. repeated evaluations of integrals [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF], where the ρ 0 a 's are of the form ρ 0 za,Na .

ISA iterations

The ISA method was originally introduced in [START_REF] Timothy | Redefining the atom: atomic charge densities produced by an iterative stockholder approach[END_REF][START_REF] Lillestolen | Atomic charge densities generated using an iterative stockholder procedure[END_REF] as the iterative scheme described below, which was then interpreted in [START_REF] Heidar-Zadeh | Information-theoretic approaches to atomsin-molecules: Hirshfeld family of partitioning schemes[END_REF] as a specific instance of the generic iterative AIM algorithm for

K 0 = K 0 ISA := ρ 0 = (ρ 0 a ) 1≤a≤M ∈ (X r + ) M . For K 0 za = K 0 ISA = X r + and
Step 2, the solution to problem [START_REF] Hirshfeld | XVII. Spatial partitioning of charge density[END_REF], when it exists, has a particularly simple form. Indeed, denoting by

ρ a s (r) : = S 2 ρ a (|r|σ) dσ (29)
= 1 4π ˆπ 0 sin θ ˆ2π 0 ρ a (r sin θ cos φ e x + r sin θ sin φ e y + r cos θ e z ) dφ dθ, with ffl

S 2 := 1 4π ´S2
, the spherical average of ρ a , we have the following result. Lemma 5. Let ρ a ∈ X + . Then, the constrained optimization problem

inf ρ 0 a ∈K 0 ISA , ´R3 ρ 0 a = ´R3 ρa s KL (ρ a |ρ 0 a ) (30)
is well-posed and its solution is ρ a s .

It is convenient to formulate the ISA method in terms of the radial functions

w a (r) := S 2 ρ a (rσ) dσ, (31)
which belong to the closed convex cone

Y + := {w ∈ L 1 (R + ; (r 2 + r 3 ) dr) ∩ L ∞ (R + ) | w ≥ 0 a.e., lim r→+∞ w(r) = 0} of the Banach space Y := L 1 (R + ; (r 2 + r 3 ) dr) ∩ L ∞ (R + )
, and are such that

∀r > 0, ∀σ ∈ S 2 , ρ a s (rσ) = w a (r).
The functions w a are the main variables used in the practical implementation of the ISA method. At the continuous level, the ISA algorithm can thus be formulated as follows:

ISA algorithm:

• Initialization: choose an initial guess for w (0) := (w

(0) a ) 1≤a≤M ∈ Y M + such that M a=1 w a (|r -R a |) > 0 in R 3 .
• Iteration m ≥ 1:

Step 1: set

ρ (m) a (r) := w (m-1) a (|r|) M b=1 w (m-1) b (|r -R b + R a |) ρ(R a + r),
Step 2: compute

w (m) a (r) := S 2 ρ (m) a (rσ) dσ. (32)
Introducing the linear map X M ρ → ρ s ∈ X M where for all ρ = (ρ a ) 1≤a≤M ∈ X M , ρ s := ( ρ a s ) 1≤a≤M with ρ a s defined by [START_REF] Vyacheslav | A quadrature formula for the sphere of the 131st algebraic order of accuracy[END_REF], the ISA optimization problem can be written as min

ρ∈K ρ,R S(ρ| ρ s ) (ISA optimization problem). (33)
A practical implementation of ISA can then be obtained by 1. discretizing the radially symmetric functions w (m-1) a on a suitable one-dimensional grid;

2. computing the spherical averages [START_REF] Liu | AMOEBA+ classical potential for modeling molecular interactions[END_REF] for each a and each grid point of the radial grid, using e.g. Lebedev integration method.

More details about the practical aspects of the implementation will be given in Section 5.

Linear approximation of the ISA (L-ISA) method

Instead of working with the three-dimensional, but radially symmetric, functions given in K 0 z = K 0 z,L-ISA , it is more convenient to work with their generating function defined on R + . For this purpose, let us denote by g za,k : R + → R + the function such that g za,k (r) = g za,k (|r|) for all r ∈ R 3 , 1 ≤ a ≤ M and 1 ≤ k ≤ m za . Let us also denote by

V := Span { g za,k , 1 ≤ a ≤ M, 1 ≤ k ≤ m za } , V + := Span + { g za,k , 1 ≤ a ≤ M, 1 ≤ k ≤ m za }
(where the notation Span + denotes the set of linear combinations with non-negative coefficients of the g za,k 's). Then, it holds that V ⊂ Y , V + ⊂ Y + and the L-ISA algorithm can be rewritten as follows:

L-ISA algorithm:

• Initialization: choose an initial guess for v (0) := (v

(0) a ) 1≤a≤M ∈ V M + such that M b=1 v b (|r - R b |) > 0 in R 3 .
• Iteration m ≥ 1:

Step 1: set

ρ (m) a (r) := v (m-1) a (|r|) M b=1 v (m-1) b (|r -R b + R a |) ρ(R a + r), (34)
Step 2: compute

v (m) a ∈ argmin v a ∈ V + , 4π ´R+ v a (r)r 2 dr = N (m) a ˆR3 ρ (m) a (r) log ρ (m) a (r) v a (|r|) dr, (35)
with

N (m) a := ˆR3 ρ (m) a (r) dr.
In this case, the subproblems (25) (and equivalently subproblems [START_REF] Mackerell | All-atom empirical potential for molecular modeling and dynamics studies of proteins[END_REF] with

ρ a = ρ (m) a ) to be solved in Step 2 are of the form inf v a ∈ V + , 4π ´R+ v a (r)r 2 dr = N a ˆR3 ρ a (r) log ρ a (r) v a (|r|) dr, (36)
with N a := ´R3 ρ a . The existence of minimizers to this problem is discussed in Lemma 11. Problem (36) can be rewritten equivalently, since the functions g za,k are L 1 -normalized, as

inf ca=(c a,1 ,••• ,ca,m za )∈R mz a + mz a k=1 c a,k = N a ˆR3 ρ a (r) log       ρ a (r) mz a k=1 c a,k g za,k (r)       dr. (37)
Note that the minimization problem (37) differs from the original GISA method [START_REF] Verstraelen | The conformational sensitivity of iterative stockholder partitioning schemes[END_REF] where an

L 2 norm distance between v a = mz a k=1 c a,k g za,k and the spherical average of ρ a (namely ρ (m) a , at iteration m) is minimized, though under the same constraint mz a k=1 c a,k = N a (see Section 2.6).
Using the function w a defined in [START_REF] Lillestolen | Atomic charge densities generated using an iterative stockholder procedure[END_REF], we see that

ˆR3 ρ a (r) log       ρ a (r) mz a k=1 c a,k g za,k (r)       dr = ˆR3 ρ a (r) log (ρ a (r)) dr -4π ˆ+∞ 0 r 2 w a (r) log mz a k=1 c a,k g za,k (r) dr.
It follows that the minimizers of (37) can be obtained by solving the problem

inf ca=(c a,1 ,••• ,ca,m za )∈R mz a + mz a k=1 c a,k = N a F a (c a ), (38)
where

F a (c a ) := - ˆ+∞ 0 r 2 w a (r) log mz a k=1 c a,k g za,k (r) dr. (39)
For functions g za,k decaying extremely fast at infinity, it may happen that F a (c a ) = +∞ for w a 's corresponding to physically admissible ρ a 's. On the other hand, if the g za,k decay at most exponentially fast (resp. not faster than a Gaussian function) and the first-order moments (resp. the second-order moments) of the function ρ to be decomposed is finite, then the function F a is strictly convex and continuous on the simplex

c a = (c a,1 , • • • , c a,mz a ) ∈ R mz a + mz a k=1 c a,k = N a
for all N a > 0, and therefore has a unique minimizer. In addition, it is smooth, so that the minimizer can be computed efficiently for small values of m za by standard numerical optimization algorithms. The components of the gradient and Hessian of F a are given respectively for all

1 ≤ i, j ≤ m za by ∂F a ∂c a,i (c a ) = - ˆ+∞ 0 r 2 w a (r) g za,i (r) mz a k=1 c a,k g za,k (r) dr, ∂ 2 F a ∂c a,i ∂c a,j (c a ) = ˆ+∞ 0 r 2 w a (r) g za,i (r) g za,j (r) mz a k=1 c a,k g za,k (r)
2 dr.

Nonlinear approximation of the ISA (NL-ISA)

Proceeding as in the previous section, the subproblems ( 25) to be solved in Step 2 are now of the form inf

ca∈R mz a + , αa∈(0,+∞) mz a mz a k=1 c a,k = N a F a (c a , α a ), (40)
denoting

c a = (c a,1 , • • • , c a,mz a ), α a = (α a,1 , • • • , α a,mz a )
and where

F a (c a , α a ) := - ˆ+∞ 0 r 2 w a (r) log mz a k=1 c a,k g za,k,α a,k (r) dr.
Still for the sake of simplicity, we assume that the α a,k 's are positive real numbers, but more complicated settings can be considered as well.

Not much can be said about the optimization [START_REF] Noel | cclib: A library for packageindependent computational chemistry algorithms[END_REF] without additional assumptions on the functions g za,k,α a,k . It may have no solution for some admissible w a 's, and multiple solutions for other admissible w a 's.

GISA iterations

The aim of this section is to describe in details the original GISA method. The GISA iterations can be formalized as follows.

GISA algorithm:

• Initialization: Choose ρ 0,(0) ∈ K 0 GISA .

• Iteration m ≥ 1:

Step 1: Find ρ (m) ∈ K ρ,R solution to

ρ (m) = arg min ρ∈K ρ,R S(ρ|ρ 0,(m-1) ). (41)
Step 2: Find ρ 0,(m) ∈ K 0 GISA solution to ρ 0,(m) = arg min

ρ 0 ∈K 0 GISA | N (ρ 0 )=N (ρ (m) ) ρ 0 -ρ (m) 2 L 2 . ( 42
)
The solution to ( 41) is explicit, while (42) can be split into the M independent minimization problems ρ 0,(m) a = arg min

ρ 0 a ∈K 0 za,GISA | ´R3 ρ 0 a = ´R3 ρ (m) a ρ (0) a -ρ (m) a 2 L 2 . ( 43
)
It is convenient to reformulate the GISA iterations using as main variables the coefficients

c = (c a ) 1≤a≤M ∈ R p + with c a = (c a,k ) 1≤k≤mz a ∈ R mz a +
and p = M a=1 m za (see Eq. ( 8)). The iterations ( 41)-( 42) can be rewritten as

ρ (m) a (r) = mz a k=1 c (m-1) a,k ζ α a,k (r) M b=1 mz b k=1 c (m-1) b,k ζ α b,k (r + R a -R b ) ρ(r + R a ), (44)
and

c (m) a = arg min ca∈R mz a + , mz a k=1 c a,k = ˆR3 ρ (m) a mz a k=1 c a,k ζ α a,k -ρ (m) a 2 L 2 = arg min ca∈R mz a + , mz a k=1 c a,k = ˆR3 ρ (m) a 1 2 c T a S za c a -c T a b (m) a , (45) (46)
with

[S za ] k = 2 ˆR3 ζ α a,k (r) ζ α a, (r) dr = 2 π √ π (α a,k α a, ) 3 2 √ α a,k + α a, 3 , b (m) a,k = ˆR3 ζ α a,k (r) ρ (m) a (r) dr.
Again, the functions ρ and the function ρ. In practice, m za 's are small so that (45) can be solved very easily once these integrals have been computed, using standard routines for quadratic programming problems [START_REF] Goldfarb | Dual and primal-dual methods for solving strictly convex quadratic programs[END_REF][START_REF] Goldfarb | A numerically stable dual method for solving strictly convex quadratic programs[END_REF].

The existence of a fixed point to the GISA iterations can be proved using Brouwer fixed-point theorem. On the other hand, the convergence of the GISA iterations for any initial guess is not guaranteed.

MB-ISA iterations

The MB-ISA method [START_REF] Verstraelen | Minimal basis iterative stockholder: atoms in molecules for force-field development[END_REF] was originally defined through a Lagrangian that does not seem to be canonically associated to a well-defined constrained optimization problem. It is more satisfactory from a mathematical point of view, to reformulate it as follows:

1. the density ρ is decomposed as a sum of atomic-shell densities. The set of admissible atomic-shell-in-molecule (ASIM) decompositions of the density ρ is

K ρ,R,z = ρ = (ρ a,k ) 1≤a≤M, 1≤k≤mz a ∈ X M + M a=1 mz a k=1 ρ a,k (• -R a ) = ρ .
The integer m za is the number of electronic shells in atom a; it is a function of the atomic number z a (hence the subscript z = (z 1 , . . . , z M ) ∈ (N * ) M collecting the atomic numbers of the atoms contained in the molecular system of interest);

2. the set of proatomic-shell densities is defined as

K 0 MB-ISA = ρ 0 (r) = c α 3 8π e -α|r| , c ≥ 0, α > 0 ,
(note that this set is independent of the chemical element and the shell) and the set of admissible pro-ASIM densities as

K 0 z = ρ 0 = (ρ 0 a,k ) 1≤a≤M, 1≤k≤mz a , ρ 0 a,k ∈ K 0 MB-ISA ;
3. the atomic-shell relative entropy is defined as

S sh (ρ|ρ 0 ) := M a=1 mz a k=1 s KL (ρ a,k |ρ 0 a,k );
4. the optimal MB-ISA decomposition is finally obtained by solving the constrained minimization problem inf

(ρ,ρ 0 )∈C ρ,R,z S sh (ρ|ρ 0 ), (47)
where

C ρ,R,z = (ρ, ρ 0 ) ∈ K ρ,R,z × K 0 z | M 0 (ρ) = M 0 (ρ 0 ), M 1 (ρ) = M 1 (ρ 0 ) . Denoting by K = M a=1
m za and

X n := ρ ∈ X ˆR3 |r| n |ρ(r)| dr < ∞ , the linear maps M n : X K n → R K are defined by ∀ρ = (ρ a,k ) 1≤a≤M, 1≤k≤mz a ∈ X K n , [M n (ρ)] a,k = ˆR3 |r| n ρ a,k (r) dr.
Note that X 0 = X and M 0 = N . As for any c ∈ R, ˆR3 c α 3 8π e -α|r| dr = c and

ˆR3 |r|c α 3 8π e -α|r| dr = 3c α ,
the function

G 0 K 0 z (ρ) := argmin S sh (ρ|ρ 0 ), ρ 0 ∈ K 0 z , M 0 (ρ 0 ) = M 0 (ρ), M 1 (ρ 0 ) = M 1 (ρ) , has an explicit expression: for all ρ = (ρ a,k ) 1≤a≤M, 1≤k≤mz a ∈ K ρ,R,z , [G 0 K 0 z (ρ)] a,k (r) = c(ρ a,k ) α(ρ a,k ) 3 8π e -α(ρ a,k )|r| , with c(ρ a,k ) = ˆR3 ρ a,k (r) dr and α(ρ a,k ) = 3 c(ρ a,k ) ˆR3 |r|ρ a,k (r) dr -1
.

MB-ISA can be reformulated in an easily implementable iterative algorithm in which the main variable is y := (c a,k , α a,k ) 1≤a≤M, 1≤k≤mz a ∈ R K + : MB-ISA algorithm:

• Initialization: choose (c (0) a,k ) 1≤a≤M,1≤k≤mz a and (α (0) a,k ) 1≤a≤M,1≤k≤mz a -see section 5.1 for details. • Iteration m ≥ 1: set c (m) a,k = ˆR3 c (m-1) a,k α (m-1) a,k 3 e -α (m-1) a,k |r| M b=1 K b =1 c (m-1) b, α (m-1) b, 3 e -α (m-1) b, |r+Ra-R b | ρ(R a + r)dr, (48) α (m)
a,k = 3 c (m) a,k        ˆR3 c (m-1) a,k α (m-1) a,k 3 e -α (m-1) a,k |r| M b=1 K b =1 c (m-1) b, α (m-1) b, 3 e -α (m-1) b, |r+Ra-R b | |r|ρ(R a + r)dr        -1 , (49)
which only requires repeated evaluations of integrals of the form (26) (with ρ(R a +r) replaced by |r|ρ(R a + r) in the evaluation of ( 49)).

Mathematical analysis of AIM decomposition methods

In this section, we make a specific focus on the AIM decomposition methods falling into the scope of the unified formulation [START_REF] Goldfarb | A numerically stable dual method for solving strictly convex quadratic programs[END_REF]

: find ρ opt , ρ 0,opt ∈ K ρ,R × K 0 solution to inf ρ, ρ 0 ∈ K ρ,R × K 0 N (ρ) = N (ρ 0 ) S ρ|ρ 0 , (50)
where K 0 is a non-empty subset of (X r + ) M . Note that the only freedom to change the method is the set K 0 and the remaining mathematical structure remains identical.

For the ISA and L-ISA method, the set K 0 is convex and closed for the topology of X M . Consequently, ( 50) is a strictly convex problem, so that its solution, if it exists, is unique. On the other hand, the existence of a minimizer is not obvious a priori due to possible loss of compactness; this issue is addressed in Theorem 6.

For the Hischfeld-I method, the existence of a minimizer follows from the very simple structure of K 0 Hirschfeld-I by a simple compactness argument, while the uniqueness of the AIM decomposition is not guaranteed due to the non-convexity of K 0 Hirschfeld-I . For NL-ISA, neither existence, nor uniqueness, can be proved in the absence of further assumptions on ρ and the functions g za,k,α a,k .

In the rest of the section, we focus on ISA and L-ISA. Let us emphasize that although the sets K 0 L-ISA and K 0 GISA have the same mathematical structure, the arguments below cannot be applied to GISA since the latter approach does not fit into the framework [START_REF] Shi | Polarizable atomic multipole-based AMOEBA force field for proteins[END_REF]. The existence and uniqueness of a minimizer to the corresponding optimization problem is established in Theorem 6. The aim of Theorem 7 is to study the convergence properties of ISA and L-ISA. We prove that the objective functional is non-increasing along the iterations and that the L 2distance between to successive iterates tends to zero as the iteration number tends to infinity. In the L-ISA case, we prove the convergence of the L-ISA algorithm towards the unique minimizer of the L-ISA optimization problem. Some additional mathematical properties of the minimizer of the ISA optimization problem in the diatomic case are collected in Propositions 8 and 9. All the proofs are postponed until Section 6. Theorem 6. Let ρ ∈ X + . Then, 1. For K 0 = K 0 ISA , problem (50) admits a unique minimizer ρ opt , ρ 0,opt ∈ K ρ,R × K 0 ISA .

2. The same holds for K 0 = K 0 L-ISA under the following additional assumption

∀1 ≤ a ≤ M, ρ( • + R a ) log min 1≤k≤mz a g za,k (•) ∈ L 1 (R 3 ). ( 51
)
3. In both cases (ISA and L-ISA), ρ opt := ρ opt a 1≤a≤M and ρ 0,opt := ρ 0,opt a 1≤a≤M satisfy the Hirschfeld relation

∀1 ≤ a ≤ M, ρ opt a (r) = ρ 0,opt a (r) M b=1 ρ 0,opt b (r -R b + R a ) ρ(r + R a ), (52)
with the convention that ρ opt a (r) = 0 whenever ρ 0,opt a (r) = 0.

4. In the ISA case, we have in addition

ρ 0,opt a = ρ opt a s . (53)
At first sight, it is not clear at all that the fixed-point iterations ( 23)-( 24) for solving problem (50) should converge. Indeed, the constraint N (ρ) = N (ρ 0 ) is not taken into account in the minimization subproblem [START_REF] Hirshfeld | Bonded-atom fragments for describing molecular charge densities[END_REF], and it is therefore not clear a priori that the algorithm ( 23)-( 24) possesses a Lyapunov functional (i.e. a functional which decreases at each iteration, a very useful tool in convergence proofs). It turns out that the function S itself is in fact a Lyapunov functional: S ρ (m+1) |ρ 0,(m+1) ≤ S ρ (m) |ρ 0,(m) for all m. This is a consequence of the special properties of the Kullback-Leibler divergence, which gives rise to the following result.

Theorem 7. Let ρ ∈ X + be such that ρ > 0 almost everywhere, and ρ 0,(0) = (ρ 0,(0) a

) 1≤a≤M ∈ K 0 (with K 0 = K 0 ISA or K 0 = K 0 L-ISA
) be an initial guess such that ρ 0,(0) a > 0 almost everywhere for all 1 ≤ a ≤ M , and s KL (ρ|ρ 0,(0) ) < +∞, where ρ 0,(m) = M a=1 ρ 0,(m) a ( • -R a ). In the L-ISA case, we assume in addition that (51) is satisfied and all the g za,k 's are positive almost everywhere. Then 1. For both ISA and L-ISA, the AIM iterations (23)-( 24) are well-defined and S is a Lyapunov functional of the algorithm. More precisely, we have for all m ≥ 1, (ρ (m) , ρ 0,(m) ) ∈ C ρ,R and the estimates

S(ρ (m-1) |ρ 0,(m-1) ) -S(ρ (m) |ρ 0,(m) ) ≥ 1 2 ρ L ∞ M a=1 ρ (m) a -ρ (m-1) a 2 L 2 ≥ 0, (54)
and

s KL (ρ|ρ 0,(m-1) ) -s KL (ρ|ρ 0,(m) ) ≥ 1 2 ρ L ∞ M a=1 ρ (m+1) a -ρ (m) a 2 L 2 ≥ 0. ( 55
)
As a consequence, it holds that

ρ (m) -ρ (m-1) -→ m→+∞ 0 strongly in L 2 (R 3 ) M . (56)
2. For L-ISA, the sequence (ρ 0,(m) , ρ (m) ) m∈N converges to the unique minimizer (ρ opt , ρ 0,opt ) of ( 50) in some sense. More precisely, (ρ 0,(m) ) m∈N converges strongly to ρ 0,opt in L p (R 3 ) for any 1 ≤ p ≤ +∞, and (ρ (m) ) m∈N converges to ρ opt for the weak-* topologies of M b (R 3 ) M and L ∞ (R 3 ) M , for the weak topology of L p (R 3 ) M for any 1 < p < ∞, and for the strong topology of L ∞ loc (R 3 ). Here, M b (R 3 ) denotes the Banach space of the bounded (signed) Radon measures on R 3 endowed with the total variation norm.

We finally establish additional mathematical properties of the ISA AIM decomposition ρ opt = (ρ opt 1 , ρ opt 2 ) of diatomic systems (K 0 = K 0 ISA and M = 2). We first show that one of the atomic densities ρ opt a is identically equal to zero if and only the total density ρ is a radial function centered on the other atom. Proposition 8. Let us assume that M = 2 and ρ ∈ X + is continuous and positive on R 3 . Let ρ opt , ρ 0,opt be the unique minimizer of problem [START_REF] Löwdin | Approximate formulas for many-center integrals in the theory of molecules and crystals[END_REF] for

K 0 = K 0 ISA . Then, ρ opt 1 = 0 if and only if ρ(• + R 2 ) = ρ opt 2 = ρ opt 2 s .
We also have the following result, which proves stronger regularity results on ρ under appropriate assumptions on the regularity of ρ. Proposition 9. Assume that M = 2 and ρ ∈ X + ∩ C 0,1 (R 3 ). Let (ρ opt , ρ 0,opt ) be the unique minimizer of problem (50) for K 0 = K 0 ISA . For a = 1, 2, let w a ∈ Y + be the function defined by ∀r > 0, ∀σ ∈ S 2 , w a (r) := ρ opt a s (rσ).

Then, 1. for almost all r ≥ 0, w a (r) ∈ {0, w + a (r)}, where the functions w + a : R + → R + are bounded by M ρ := max R 3 ρ, continuous on (0, +∞), and Lipschitz in the neighborhood of any point of (0, +∞) at which they are positive; 2. if in addition, the density ρ is C k away from the centers R a for a = 1, 2 and the functions w a are essentially bounded away from zero on any compact subset of [0, +∞), then the functions w a are Lipschitz on [0, +∞) and

C k on (0, R) ∪ (R, +∞) with R := |R a -R b |.

Distributed Multipole Analysis

Other schemes allow to represent the electronic density as a sum of atom-centered contributions.

Mulliken population analysis [START_REF] Mulliken | Electronic population analysis on LCAO-MO molecular wave functions. i[END_REF], and Distributed Multipole Analysis (DMA) [START_REF] Stone | Distributed multipole analysis, or how to describe a molecular charge distribution[END_REF][START_REF] Stone | Distributed multipole analysis[END_REF] -which can be viewed as a generalization of Mulliken scheme to higher-order multipole moments -belong to the latter. These schemes are fully basis-set oriented, in the sense that they use as a starting point a decomposition of the density of the form

ρ(r) = M k=1 f k (r -a k ) (57)
where A := {a k } 1≤k≤M is a set of points of the physical space R 3 , and f k : R 3 → R 3 functions for which the multipoles

(Q k m ) = K m ˆR3 R m (r)f k (r)dr, (58)
with R m (r) = |r| Y m r |r| (real solid harmonics using L 2 (S 2 )-orthonormal real spherical harmonics), are easily computable. The coefficients K m are real normalization coefficients (for instance,

K 0 0 = √ 4π and K m 1 = √ 4π
3 ). The purpose of DMA methods is to redistribute the multipoles Q k m carried by the points a k to multipoles Q j m carried by some chosen points S j , called DMA expansion sites. The set S = {S j } 1≤j≤J can coincide with the set of atomic positions, but can include other points. For instance, Mulliken population analysis defines overlap populations (that can be attributed to the bonds and allocated to bond centers), and the DMA can also include charge contributions from e.g. bond centers and from many other points.

The rationale of the distributed multipole analysis is thus the following:

1. first, the multipole moments (Q k m ) are computed. The latter will be coined multipole moments at natural centers in the following. In the case of an electronic density computed in a basis of Gaussian Type Orbitals (GTOs) -as used in most quantum chemistry applications -, there is an analytical expression for these elementary moments, allowing for a very fast computation; 2. then, these multipoles are redistributed to the DMA expansion sites S j , giving rise to the redistributed multipoles Q j m carried by the sites S j . These two aspects are discussed in the following two sections.

Multipole moments with respect to natural centers

We recall here the procedure for rewriting as [START_REF] Verstraelen | The conformational sensitivity of iterative stockholder partitioning schemes[END_REF] an electronic density ρ issued from a quantum chemistry calculation in a Gaussian type orbital (GTO) basis set, and computing the corresponding multipoles Q k m . First, ρ is decomposed as

ρ(r) = N b µ,ν=1 P µν χ µ (r -R aµ )χ ν (r -R aν ), (59)
where (χ µ ) 1≤µ≤N b is a basis of primitive Gaussian polynomial functions, a µ the label of the nucleus χ µ is attached to, and (P µν ) 1≤µ,ν≤N b the density matrix in this basis. The matrix (P µν ) 1≤µ,ν≤N b is easily obtained from the density matrix in the contracted GTO basis set used for the computation of ρ by a suitable linear transform. Each χ µ is of the form

χ µ (r) = N µ R mµ µ (r)e -ζµ|r| 2 , (60)
with µ and m µ its angular momentum quantum numbers, ζ µ the Gaussian exponent controlling its decay rate, and N µ a normalization constant. Since

e -ζµ|r-Ra µ | 2 e -ζν |r-Ra ν | 2 = K µν e -(ζµ+ζµ)|r-Rµν | 2 , ( 61
)
with

K µν = e - ζµζµ ζµ+ζν |Ra µ -Ra ν | 2 , and R µν = ζµRa µ +ζν Ra ν ζµ+ζν
, it holds that

χ µ (r -R aµ )χ j (r -R aν ) = K µν N µ N ν R mµ µ (r -R aµ ) R mν ν (r -R aν ) e -(ζµ+ζµ)|r-Rµν | 2 , (62)
so that the electronic density given by ( 59) can be written as [START_REF] Verstraelen | The conformational sensitivity of iterative stockholder partitioning schemes[END_REF], where the functions f k : R 3 → R are Gaussian polynomials of the form (62), and A := {a k } 1≤k≤M a set of points in the physical space, (strictly) containing the positions of the nuclei.

The key quantities for the calculation of the distributed multipole moments are thus the multipole moments with respect to natural expansion centers R µν :

Q µν m = P µν K m ˆR3 R m (r)χ µ (r -R µ + R µν )χ ν (r -R ν + R µν )dr. (63)
Given the expression (62), the series of multipole moments {Q µν m } ∈N,m=-,..., is finite and can be easily evaluated numerically, by expanding the product of the two solid harmonics in [START_REF] Wu | Automation of AMOEBA polarizable force field parameterization for small molecules[END_REF] as a weighted sum of homogeneous polynomials, multiplying it with the solid harmonic R m (r) (homogeneous polynomial) and finally using the series of tabulated integrals:

I q = ˆ+∞ 0 u q e -u 2 du = Γ q+1 2 2 , (64)
for all needed q integer values, where Γ is the Euler Gamma function.

Redistribution formula

Once multipole moments at natural centers are computed, a lower number of expansion sites (than the total number of overlap centers i.e. |A|) may be desired to perform the distributed multipole expansion, as the number of overlap centers grows quadratically with the number of atoms or basis size. For instance, only nuclei positions may be retained as final centers (e.g.

to use these local multipole moments in polarizable force fields such as AMOEBA [START_REF] Shi | Polarizable atomic multipole-based AMOEBA force field for proteins[END_REF]), or all nuclei and bond centers [START_REF] Stone | Distributed multipole analysis, or how to describe a molecular charge distribution[END_REF][START_REF] Vigné | The exact multicenter multipolar part of a molecular charge distribution and its simplified representations[END_REF]. The multipole moments at natural centers, which are easy to compute (see above), have thus to be redistributed to the chosen final expansion sites. We recall that the set of final expansion sites is a choice taken by the user.

This requires two ingredients. First, the introduction of redistribution weights in order to distribute the multipoles from natural sites not belonging to the set of final sites. Second, a transformation of a multipole corresponding to one expansion center to a multipole corresponding to a different expansion center which is a linear transformation that is known in the Fast Multipole Method (FMM) framework as M 2M -operator [START_REF] Fortin | Algorithmique hiérarchique parallèle haute performance pour les problèmes à N-corps[END_REF] (used here in the context of realvalued multipolar coefficients).

We first focus on the latter. We briefly present the resulting formula, and refer to [START_REF] Benda | Modélisation multi-échelle de nano-capteurs à base de polymères conjugués pour la qualité de l'eau[END_REF] for a detailed derivation:

Q µν→j m = Q µν→j ,-m C + (-1) m Q µν→j m C √ 2 if m > 0, (65)
Q µν→j m = (-1) |m| Q µν→j ,-m C - Q µν→j m C √ 2 if m < 0, (66)
Q µν→j ,0 = Q µν→j ,0 C , (67)
where the complex M 2M -operator is given by

Q µν→j m C = =0 m =- c m m |R µν -S j | -Y m-m - R µν -S j |R µν -S j | Q µν ,m C , (68)
where Y m denote the complex (L 2 (S 2 )-orthonormal) spherical harmonics and c m ,m = +m +m -m -m a normalization factor. Complex arithmetic operations are used since it simplifies the presentation of the formulae. Finally, the complex multipole moments are simply given by

Q µν m C = P µν K m ˆR3 |r| Y m r |r| χ µ (r -R µ + R µν )χ ν (r -R ν + R µν )dr. (69)
The multipole moments are redistributed to the final expansion sites S j using a formula of the form

Q j m = (µ,ν) C µν→j Q µν→j m R , (70)
where the C µν→j 's are user-specified redistribution weights, representing the part of the multipolar distribution centered at R µν allocated to the expansion site S j . It holds

C µν→j ≥ 0, p j=1 C µν→j = 1, C µν→j = 1 if S j = R µν .
In particular, if R µν is one of the final expansion sites, then Q µν m is fully allocated to this site. As a consequence, Q µν→j m needs to be computed only if R µν is not one of the final expansion sites. It is up to the user to choose both the final expansion sites and the values of the weight coefficients C µν→j . A complete flexibility is allowed in our DMA code [START_REF] Benda | Dma[END_REF]. For instance, following Stone's neighbor-takes-it-all strategy implemented in GDMA code [START_REF] Anthony | Distributed multipole analysis: Stability for large basis sets[END_REF] would amount to take:

C µν→j = 1 if S j is the final expansion site nearest to R µν , 0 otherwise, (71)
or 1 q for each of the q final sites equally close to R µν . This strategy means that multipole moments at natural centers are only redistributed locally, at the nearest neighbor final site. An obvious drawback of this approach is that the redistributed moments Q j m are clearly not differentiable, nor even continuous, with respect to the nuclear coordinates. We have also implemented, and provide as a possible user-choice, the more balanced rule of Vigne-Maeder et al. [START_REF] Vigné | The exact multicenter multipolar part of a molecular charge distribution and its simplified representations[END_REF], which amounts to choose the following weight coefficients:

C µν→j = 1 |S j -Rµν | k 1 |S k -Rµν | , ( 72
)
where the sum in the denominator runs over the p chosen final DMA sites. This time, multipole moments at natural centers that are overlooked in the set of DMA expansion sites are redistributed to all final DMA sites {S j } j , with larger contributions to nearer sites.

Multipolar expansion of the electrostatic potential

Once the local distributed multipole moments are computed, the electrostatic potential generated by the charge density distribution (that the latter summarize) can be evaluated by a simple sum, instead of a three-dimensional integral needed to evaluate the exact, quantum potential:

V (r) = ˆR3 ρ(r ) |r -r | dr .
The multipolar expansion of the total potential with respect to the expansion centers R µν then writes:

V M (r) = N b µ,ν=1   +∞ =0 4π (2 + 1) m=- Q µν m Y m r-Rµν |r-Rµν | K m |r -R µν | +1   , (73)
where the sum over is actually finite and terminates for = µ + ν , the sum of the degrees of the two real solid harmonics involved in the primitive GTOs χ µ and χ ν respectively. The multipolar expansion V M is an excellent approximation of V far enough from the molecule, and a very poor approximation in the region of space where the electronic density is large. In the intermediate region, it is valid up to a penetration term [START_REF] Liu | AMOEBA+ classical potential for modeling molecular interactions[END_REF][START_REF] Rackers | An optimized charge penetration model for use with the AMOEBA force field[END_REF].

After redistribution of the multipole moments to the DMA final expansion centers, the potential generated by the DMA distributed local multipole moments writes:

V DMA M (r) = J j=1   +∞ =0 4π (2l + 1) m=- Q j m Y m r-S j |r-S j | K m |r -S j | +1   ( 74
)
where the total multipole moments ( Q j m ) 1≤j≤J with respect to final expansion centers S j include the contributions of multipole moments naturally centered at S j (if any), as well as of the reallocated (e.g. neighboring) multipoles -see Eq. (70). The series of multipole moments

Q j m l,m
with origin S j is in general infinite, and must be truncated in practice. For this reason, V DMA M is always an approximation of V M as soon as one of the natural centers R µν is not kept as a final DMA expansion site, even away from the numerical support of the density ρ.

Let us finally mention that the new version of DMA [START_REF] Anthony | Distributed multipole analysis: Stability for large basis sets[END_REF] can be thought as a mix between this DMA basis-set oriented method and real-space oriented (partitionning) methods, as it associates a redistribution of the density matrix contributions for some (the least diffuse) pairwise products of GTOs with a real-space integration for the contributions of the others (the most diffuse) pairwise products of GTOs, thus allowing to stabilize local multipole moments in the limit of a very large Atomic Orbitals basis for the molecular density ρ.

Numerical results

We have implemented the four methods of the ISA family described in Section 2, namely ISA, GISA, L-ISA and MB-ISA, in Python language with a specific focus on accurate quadrature schemes for diatomic molecules. This serves as a proof-of-concept to highlight the theoretical results obtained in this article and can by far not be considered as an exhaustive computational study.

The Distributed Multipole Analysis (DMA) method (Section 4), has also been implemented in Python language in a new, more modular package, allowing for user-specified redistribution strategies (see Eq.( 70) and e.g. the two possible choices of coefficients ( 71) and ( 72)). Our code has been validated by comparing the results of local multipole moments (up to hexadecapoles) obtained on test molecules (H 2 O, NH 3 , CH 4 , benzene, and ions such as ClO -) with those obtained by the GDMA code [START_REF] Anthony | Distributed multipole analysis: Stability for large basis sets[END_REF], using Stone neighbor-takes-it-all redistribution strategy (equation 71). Several basis sets were tested for each molecule. We obtained identical distributed (atomic) local multipole moments up to hexadecapoles, up to numerical accuracy. All these codes are built upon the functionalities of cclib library [START_REF] Noel | cclib: A library for packageindependent computational chemistry algorithms[END_REF], which allows to parse in a generic way the output of most quantum chemistry packages (14 different packages are supported by this library), without any specific parsing required by the user. This new feature enables an easy comparison of DMA or ISA partitionning results among the most popular quantum chemistry codes, providing more interoperability. The parsing by the existing Stone's GDMA code [START_REF] Anthony | Distributed multipole analysis: Stability for large basis sets[END_REF] for DMA was indeed only implemented for Gaussian [START_REF] Frisch | Gaussian 16 Revision C.01[END_REF] or Psi4 [START_REF] Turney Justin | Psi4: an open-source ab initio electronic structure program[END_REF] outputs, and similarly for the existing implementations of the ISA methods in Horton code [58].

Implementational details

Solver. The ISA AIM decomposition has been computed using two different methods: (i) the fixed-point algorithm described in Section 2.5.4, and (ii) an alternative scheme based on the Newton method applied to a suitable formulation of the Euler-Lagrange equations of the problem (see [START_REF] Benda | Modélisation multi-échelle de nano-capteurs à base de polymères conjugués pour la qualité de l'eau[END_REF] for further details). Both methods yielded the same results (up to numerical accuracy) in our test cases.

Discretization and numerical quadrature. The AIM methods rely either on a real-space discretization, or on a basis-space representation of the functions spanning the linear or nonlinear approximation space.

As already detailed in Section 2, the GISA, L-ISA and MB-ISA methods rely on basis-set representations of pro-atomic densities, through linear combination of normalized Gaussians (for GISA and L-ISA) or Slater functions (for MB-ISA). For GISA and L-ISA, only the coefficients of the linear combination are optimized. For MB-ISA, the coefficients of the expansion and the exponents involved in the Slater basis functions are both optimized simultaneously. Typically, a small number of shells (e.g. 4 for H and 6 for C, O and N [START_REF] Verstraelen | The conformational sensitivity of iterative stockholder partitioning schemes[END_REF]) are assigned to each atom in GISA and MB-ISA [START_REF] Verstraelen | Minimal basis iterative stockholder: atoms in molecules for force-field development[END_REF] methods.

On the other hand, the original ISA algorithm introduced by Lillestolen et al. [START_REF] Timothy | Redefining the atom: atomic charge densities produced by an iterative stockholder approach[END_REF][START_REF] Lillestolen | Atomic charge densities generated using an iterative stockholder procedure[END_REF]] (Section 2.5.4) relies on a real space radial and angular discretization by point-evaluations.

In the family of ISA-methods, several integrals arise both when computing the total mass N a of the atomic densities (imposed as a constraint in the GISA, L-ISA and MB-ISA schemes), the first-order moments of the atomic densities (imposed as a constraint in the MB-ISA scheme), or when computing spherical averages in the original ISA method. This requires a numerical quadrature which is given by a tensor product of a radial and spherical grid. The radial grid can be e.g. an 1D integration grid based on N r Gauss quadrature points extending up to some maximal R max , while the angular grid can be e.g. a Lebedev grid [START_REF] Vyacheslav | A quadrature formula for the sphere of the 131st algebraic order of accuracy[END_REF] made of N s points on the unit sphere S 2 (with their associated weights). For each atom a, the discretization points are defined as follows:

∀a ∈ [[1, M ]], ∀i ∈ [[1, N r ]], ∀j ∈ [[1, N s ]], r a i,j := R a + r a i σ j , (75)
denoting [[1, M ]] = {1, 2, . . . , M } and where R a is the position of atom a, (r a i ) i∈[ [START_REF] Adamo | Toward reliable density functional methods without adjustable parameters: The PBE0 model[END_REF]Nr]] are the (atom-dependent) radial grid points (which can be e.g. equally-spaced or logarithmic), and (σ j ) j∈[ [START_REF] Adamo | Toward reliable density functional methods without adjustable parameters: The PBE0 model[END_REF]Ns]] are the Lebedev grid points grids used for the quadrature on the unit sphere. The molecular density ρ is precomputed at these points r a i,j , yielding tabulated values ρ a i,j = ρ(r a i,j ) used throughout the algorithm repeatedly.

In the special case of symmetric systems, such as for some diatomic molecules (associated to an electronic density invariant by any rotation around the axis of the molecule), the spherical average on the unit sphere S 2 can be replaced by a simple one-dimensional integration on the θ variable using N s points, e.g. using a Gauss-Legendre quadrature by aligning the atoms on the z-axis. All integrals on the radial variable r are computed likewise.

Interpolation. In the ISA-method, where (pro-)atomic densities are known only on a real space grid, we use a piecewise linear (sometimes referred to P 1 ) interpolation between the nearest neighbors. We use such an interpolation for instance to evaluate the pro-atomic density ρ 0 a of atom a at a radial discretization value r b i associated to a near atom b. Indeed, only the values ρ 0 a (r a i ) for i ∈ [1, N r ]], at the radial grid point r a i , are known. On the contrary, for the GISA, L-ISA and MB-ISA schemes, the expressions of the pro-atomic densities are (by construction) analytical due to the expansion on Gaussian or Slater basis functions and no interpolation is thus needed for these methods.

Optimization. For the resolution of the GISA problem, in particular for the minimization of the quadratic functional under linear (equality and inequality) constraints [START_REF] Ren | Polarizable atomic multipole water model for molecular mechanics simulation[END_REF], we use standard routines for quadratic programming problems [START_REF] Goldfarb | Dual and primal-dual methods for solving strictly convex quadratic programs[END_REF][START_REF] Goldfarb | A numerically stable dual method for solving strictly convex quadratic programs[END_REF] implemented in the QUADPROG Python function. For solving the convex minimization [START_REF] Mulliken | Electronic population analysis on LCAO-MO molecular wave functions. i[END_REF] arising in L-ISA, we use the CVXOPT Python package [START_REF] Andersen | CVXOPT: A Python package for convex optimization[END_REF].

Parameterization. For GISA and L-ISA methods, the atomic shell exponents

(α a,k ), with a ∈ [[1, M ]], k ∈ [[1, m za ]],
are fixed once and for all. The solutions of the GISA and L-ISA schemes thus a priori depend on the values of these exponents. For C, H, N and O atoms, we use the values of the exponents provided by Verstraelen et al. for the GISA method [START_REF] Verstraelen | The conformational sensitivity of iterative stockholder partitioning schemes[END_REF], obtained by fitting to the densities of isolated (neutral or ionized) atoms, so that the GISA Gaussian basis represents accurately enough all the pro-atomic densities associated to the quantum calculations on both isolated neutral atoms and ions. For other atoms, for which no values of the exponents are available, we apply the empirical formula provided by Verstraelen et al. in the MB-ISA paper [START_REF] Verstraelen | Minimal basis iterative stockholder: atoms in molecules for force-field development[END_REF]:

∀k ∈ [[1, m za ]], α a,k = 2z 1-k-1 mz a -1 a a 0 , ( 76
)
where a 0 is the value of 1 Bohr in Angström and m za the number of shells (derived from the atomic number z a ). For MB-ISA, the exponents

(α a,k ), with a ∈ [[1, M ]], k ∈ [[1, m za ]],
of the Slater functions used to expand the pro-atomic densities are optimized instead of attributing fixed values. As initial guess for these exponents, we simply take the values of the GISA (or L-ISA) exponents detailed above. These initial values of the exponents appear numerically not to matter (as regards the final fixed-point solution) in the MB-ISA scheme, which optimizes both the weight and the Slater exponent of every pro-atomic shell (see Section 2.7). Note that some choices of exponents lead to an ill-conditionned local overlap matrix S za , which can prove problematic in particular in the context of the GISA method (see Section 2.6), where a linear system involving this matrix has to be solved.

For the expansion coefficients

(c a,k ), with a ∈ [[1, M ]], k ∈ [[1, m za ]]
, which are optimized in the GISA, L-ISA and MB-ISA methods, several initial guesses that satisfy the property that the initial pro-atomic mass N

a,k = z a m za for example. Other less-balanced initial guesses, e.g. putting all the initial weight on a specific Gaussian (or Slater) basis function:

∀k ∈ [[1, m za ]], c (0)
a,k = δ k,k 0 z a , with k 0 ∈ [[1, m za ]
] and where δ i,j denotes the Kronecker delta, are also possible.

Uniqueness and dependency on initial condition. For the GISA method, there is no guarantee of the uniqueness of the solution for the expansion coefficients, and local minima are possible (i.e. a dependence of the computed solution on the initial guess for the expansion coefficients). An example found numerically for a very simple test density is provided in the Appendix. On the contrary, for L-ISA, the optimal expansion coefficients c opt a,k does not depend on the initial guess c (0) a,k due to the strict convexity. For MB-ISA, there is no guarantee of the uniqueness of the fixed-point solution, so that the numerically computed solution may depend on the initial guess. However, numerically, the MB-ISA solution appears quite robust to the guess (both in terms of expansion coefficients c a,k and of exponents α a,k ), and no local minimum has been unveiled in this study on diatomic systems.

Validation. Our implementation of the different existing ISA partitionning schemes has been validated on some test diatomic molecules, for several basis sets, by comparison to the Horton code [58], specifically for the ISA and MB-ISA scheme. Both implementations (ours and Horton's) yielded the same local multipole moments (up to quadrupoles) up to numerical accuracy on the tested molecules. For GISA and L-ISA, our implementation was validated based on custom test atomic densities made of a sum of one or several normalized Gaussians, centered at both atomic positions (diatomic molecules only).

Comparison of the different ISA schemes

In this section, we compare typical pro-atomic densities profiles and the associated atomic moments (up to second order) for diatomic molecules or ions, at the equilibrium inter-atomic distance for Hartree-Fock densities obtained with GAMESS [START_REF] Schmidt | General atomic and molecular electronic structure system[END_REF]. The GISA and L-ISA profiles look very similar but lead to slightly different atomic multipole moments (see Table 1). The GISA, L-ISA and MB-ISA pro-atomic densities are (visually) very close to the ISA reference near the first maximum of the curve -corresponding to the core electrons (1s) of the carbon atom. Near the second maximum -corresponding to the 2s electrons shell -the GISA and L-ISA profiles are also very close to the ISA reference profile, while the MB-ISA profile features a less pronounced local maximum. The GISA and L-ISA profiles then differ quite notably from the MB-ISA and ISA reference profile in the asymptotic regime.

CO molecule.

This proximity of GISA and L-ISA to the ISA reference profile seems rather due to the choice of optimized exponents α a,k for C and O atoms in GISA and L-ISA from Ref. [START_REF] Verstraelen | The conformational sensitivity of iterative stockholder partitioning schemes[END_REF]. For the N 2 molecule (also studied, results not reported), the results are indeed qualitatively similar to those obtained for CO. In the case of optimized exponents α a,k , the GISA and L-ISA Gaussian expansions of the pro-atomic densities variational space perform well to capture correctly the short-range features of ρ 0 a . However, this (visual) proximity of GISA and (quantitative) proximity of L-ISA to the reference ISA profile vanishes if the (fixed) exponents α a,k used in GISA and L-ISA are not optimized for one of the atoms in the molecule. Moreover, in the case of non-optimized exponents, only L-ISA recovers qualitatively the reference ISA profile, but not GISA (see below the case of Cl atom in ClO -ion, Figure 2, compared to Figure 3 The local multipole moments associated to the different ISA profiles of Figure 1 are reported Table 1. The results of the L-ISA scheme are overall the closest to the ISA reference, followed by MB-ISA and GISA, which appears less satisfactory (and for which the atomic monopole differs by more than 30 % from the reference ISA monopole).

ClO -ion. Figures 2 and3 display the r → log 4πr 2 ρ 0 a (r) profiles of GISA, L-ISA, MB-ISA and ISA solutions for the chlorine and oxygen atoms respectively, using a molecular density ρ computed at the RHF/6-31G(d) level of theory. Optimized exponents α a,k for Cl atom have not been reported in the literature to our knowledge. We first tried the empirical rule (76) but finally assigned six shells (in GISA and L-ISA) to the chlorine atom with the same (optimized) exponents as for the oxygen atom, as it yielded slightly better (although not satisfactory, see below) results for the GISA profiles and multipole moments. The solution of the MB-ISA method was found not to depend on this initial choice of exponents for the Cl atom (empirical rule (76) or same exponents as the O atom).

Contrary to the case of neutral molecules (CO, N 2 ) with optimized exponents for both atoms, the pro-atomic densities profiles of GISA, L-ISA and MB-ISA quantitatively differ from the reference ISA profile, with some common features (e.g. three local minima for the L-ISA and ISA reference schemes, vs. only two for MB-ISA). The GISA profile ends up with an nonphysical, too large (positive) charge (1.94 a.u.) on Cl atom (-2.84 a.u. on O atom) corresponding to a too quick (Figure 2) or (resp.) too slow (Figure 3) decrease in the Cl (resp. O) atomic density profile, despite a good short-range accordance with the ISA reference profile in the case of the O atom (using its optimized Gaussian exponents [START_REF] Verstraelen | The conformational sensitivity of iterative stockholder partitioning schemes[END_REF]).The other GISA atomic multipole moments are also either much too large or largely differing (more than 100 %) from the ISA reference multipoles. The MB-ISA profile is associated to local multipole moments in fair accordance with the ISA reference multipoles (see Table 2) despite differing qualitatively quite notably from the ISA reference profile, in particular in the asymptotic regime. The latter is a feature rather due to the expansion on Slater functions, leading to a decrease at a smaller rate.

Quite notably, unlike the GISA scheme, the L-ISA method is also associated to local multipole moments in fair accordance with the ISA reference multipoles (Table 2), similarly to the MB-ISA method, although the latter optimizes on the (Slater) exponents while the L-ISA uses fixed (non-optimized) Gaussian exponents. The L-ISA scheme, which has been shown here to be mathematically more grounded than GISA, seems to be more robust with respect to the Gaussian exponent choices and yields more physical (and closer to the ISA reference profile) results, although this assumption would have to be confirmed on larger test sets. and N s = 250 Gauss-Legendre angular points were used for both the quadratures and the ISA discretization.

Dissociation of diatomic molecules

In this section, we apply the ISA multi-center decomposition schemes implemented in this work to the study of the dissociation of LiH. We consider ground-state densities obtained by Full CI calculations in a small basis sets (computed with PySCF [START_REF] Sun | PySCF: the Python-based simulations of chemistry framework[END_REF][START_REF] Sun | Recent developments in the PySCF program package[END_REF]), for different inter-atomic distances, and extract the corresponding one-body density matrices. We then compute the corresponding ISA atomic charges and dipoles. The results are reported in Figure 4. As expected, the ISA atomic charge tends to zero in the dissociation limit, while the Li-H bond is very polarized in the covalent regime (with atomic charges of order 0.80 a.u. at the equilibrium distance 1.70 Å). The ISA Li and H atomic dipoles also tend to zero in the dissociation limit, but in a non-monotonic manner, with a minimum (resp. maximum) for H (resp. Li) atomic dipole at about 2.2 Å (resp. 3.0 Å).

Basis set convergence of DMA and ISA

In this section, we compare some typical features of the DMA and ISA methods, such as the convergence of local multipole moments with increasing basis sets. Figure 5 displays the evolution of the local atomic charges on H and O atoms of H 2 O computed with DMA (using both Stone's and Vigné-Maeder's redistribution strategies) and ISA (computed with Hipart code, now known as Horton [58]), from a molecular density computed at the DFT/PBE0/(aug)-cc-pVXZ (X ∈ {D, T, Q, 5}) level of theory [START_REF] Adamo | Toward reliable density functional methods without adjustable parameters: The PBE0 model[END_REF]. ISA atomic charges stabilize up to a few % for basis sets larger than aug-cc-pVDZ, the DMA atomic charges do not stabilize with increasing basis set size (e.g. four-fold increase of atomic charge on the oxygen atom). This was already known for DMA atomic charges computed with Stone's redistribution strategy and attributed to the presence of diffuse functions in the basis set, and had lead to the hybrid real-space and basis-space GDMA method [START_REF] Anthony | Distributed multipole analysis: Stability for large basis sets[END_REF]. We see here that the more balanced Vigné-Maeder redistribution strategy [START_REF] Vigné | The exact multicenter multipolar part of a molecular charge distribution and its simplified representations[END_REF] leads to a qualitatively similar non-converging behavior. A similar behavior is observed for atomic dipoles (Figure 6), with ISA dipoles stabilizing up to 10 % beyond aug-cc-pVDZ basis.

Proofs

Let us first elaborate on the assumptions on ρ and the ρ 0 a 's made in the lemmata, propositions and theorems to be proved. As already mentioned previously, electronic densities of molecules and clusters are continuous, positive, functions on R 3 , decaying exponentially fast at infinity. We can therefore assume without loss of generality that ρ ∈ X + , and this implies in particular that ρ log ρ ∈ L 1 (R 3 ). For the same reasons, it makes sense to restrict ourselves to proatom densities ρ 0 a belonging to X + . In addition, if the ρ a 's are nonnegative measurable functions satisfying M a=1 ρ a (• -R a ) = ρ for some ρ ∈ X + , then each ρ a satisfies 0 ≤ ρ a ≤ ρ(• + R a ), and therefore belongs to X + . This justifies a posteriori the condition ρ ∈ X M + in the definition (3) of K ρ,R .

The following technical remark is important as it allows one to justify several algebraic manipulations in the proofs below.

Remark 10. Let ρ ∈ X + . For ρ = (ρ a ) 1≤a≤M ∈ K ρ,R , we have ρ a ∈ X + , and therefore ρ a log ρ a ∈ L 1 (R 3 ). As a consequence, we have for all ρ = (ρ a ) 1≤a≤M ∈ K ρ,R and

ρ 0 = (ρ 0 a ) 1≤a≤M ∈ (X r + ) M , S(ρ|ρ 0 ) = M a=1 ˆR3 ρ a log ρ a ρ 0 a = M a=1 ˆR3
ρ a log ρ a -log ρ 0 a (in agreement with ( 12))

= M a=1 ˆR3   ρa log ρ a ∈L 1 (R 3 ) -(ρ a log(ρ 0 a )) + ∈L 1 (R 3 )∩L ∞ (R 3 ) + (ρ a log(ρ 0 a )) - ≥0    , (77)
where x + = max(x, 0) and x -= -min(x, 0) denote the positive and negative parts of x ∈ R (x + , x -≥ 0 and x = x + -x -). In addition, the function (ρ a log(ρ 0 a )) + has compact support since ρ 0 a vanishes at infinity. It follows that S(ρ|ρ 0 ) < +∞ if and only if ρ a log(ρ 0 a ) ∈ L 1 (R 3 ) for all 1 ≤ a ≤ M , and that we then have (using Fubini theorem)

S(ρ|ρ 0 ) = M a=1 ˆR3 ρ a log ρ a - M a=1 ˆR3 ρ a log ρ 0 a = M a=1 ˆR3 ρ a log ρ a - M a=1 ˆR3
ρ a s log ρ 0 a . (78)

Proof of Lemma 2

Let us define

ρ * a (r) := ρ 0 a (r) ρ 0 (r+Ra) ρ(r + R a ) if ρ 0 a (r) > 0, 0 if ρ 0 a (r) = 0 (79)
(recall that ρ 0 (r) := M a=1 ρ 0 a (r -R a )), and ρ * := (ρ * a ) 1≤a≤M as an Ansatz for the minimizer of [START_REF] Fortin | Algorithmique hiérarchique parallèle haute performance pour les problèmes à N-corps[END_REF]. Let

B 0 a := {r ∈ R 3 | ρ 0 a (r) = 0}, B 0 := {r ∈ R 3 | ρ 0 (r) = 0}, B = {r ∈ R 3 | ρ(r) = 0}.
It follows from the hypothesis (15) (i.e. since s KL (ρ|ρ 0 ) < +∞) and the conventions (12) that B 0 ⊂ B. In addition, we have B 0 = M a=1 (R a + B 0 a ), so that if r ∈ B 0 , then for all 1 ≤ a ≤ M , ρ 0 a (r -R a ) = 0, and therefore, using (79), ρ * a (r -R a ) = 0. Thus, on B 0 , there holds

M a=1 ρ * a (• -R a ) = 0 = ρ as B 0 ⊂ B. Besides, we have ∀1 ≤ a ≤ M, ∀r ∈ R 3 \ B 0 , ρ * a (r -R a ) = ρ 0 a (r -R a ) ρ 0 (r) ρ(r).
Thus, there holds that M a=1 ρ * a ( • -R a ) = ρ on R 3 \ B 0 as well, and therefore on the entire space R 3 . Since the ρ * a 's are obviously nonnegative, bounded (by ρ(r + R a )), integrable, with finite first moments, and vanishing at infinity, it holds that ρ * ∈ K ρ,R . Now, we have

S(ρ * |ρ 0 ) = M a=1 ˆR3 ρ * a log ρ * a ρ 0 a = M a=1 ˆR3 ρ * a ( • -R a ) log ρ * a ( • -R a ) ρ 0 a ( • -R a ) = M a=1 ˆR3 ρ * a ( • -R a ) log ρ ρ 0 = ˆR3 ρ log ρ ρ 0 = s KL (ρ|ρ 0 ) < +∞.
All these equalities can be rigorously justified by using arguments similar as in Remark 10 and the infimum in ( 17) therefore has a finite value. For all ρ = (ρ a ) 1≤a≤M ∈ K ρ,R such that S(ρ|ρ 0 ) < +∞, we have in addition, using the fact that ρ = M a=1 ρ a ( • -R a ),

S(ρ * |ρ 0 ) = ˆR3 ρ log ρ ρ 0 = M a=1 ˆR3 ρ a ( • -R a ) log ρ ρ 0 = M a=1 ˆR3 ρ a log ρ( • + R a ) ρ 0 ( • + R a ) .
Therefore, using again arguments as in Remark 10 to justify each equality, we get

S(ρ|ρ 0 ) -S(ρ * |ρ 0 ) = M a=1 ˆR3 ρ a log ρ a ρ 0 a - M a=1 ˆR3 ρ a log ρ( • + R a ) ρ 0 ( • + R a ) = M a=1 ˆR3 ρ a log ρ a ρ 0 ( • + R a ) ρ 0 a ρ( • + R a ) = M a=1 ˆR3 ρ a log ρ a ρ * a = S(ρ|ρ * ) ≥ 0,
for all ρ ∈ K ρ,R , which proves that ρ * is a minimizer of [START_REF] Fortin | Algorithmique hiérarchique parallèle haute performance pour les problèmes à N-corps[END_REF]. As this problem is strictly convex, it is the unique one.

Proof of Lemma 5

Obviously, ρ a s ∈ K 0 ISA = X r + and ´R3 ρ a s = ´R3 ρ a , so that ρ a s is in the minimization set of [START_REF] Timothy | Redefining the atom: atomic charge densities produced by an iterative stockholder approach[END_REF]. Reasoning as in Remark 10, we obtain that

s KL (ρ a | ρ a s ) = ˆR3 ρ a log(ρ a ) - ˆR3 ρ a s log( ρ a s ) < +∞,
and that, for all ρ 0 a ∈ K 0 ISA such that s KL (ρ a |ρ 0 a ) < +∞, it holds

s KL (ρ a |ρ 0 a ) = ˆR3 ρ a log(ρ a ) - ˆR3 ρ a s log(ρ 0 a ).
Therefore, for all ρ 0 a ∈ K 0 ISA such that s KL (ρ a |ρ 0 a ) < +∞, we have

s KL (ρ a |ρ 0 a ) = s KL (ρ a | ρ a s ) + ˆR3 ρ a s log( ρ a s ) - ˆR3 ρ a s log(ρ 0 a ) = s KL (ρ a | ρ a s ) + s KL ( ρ a s |ρ 0 a ) ≥ s KL (ρ a | ρ a s
). This implies that ρ a s is a minimizer of (30) and the infimum is finite. The set K 0 ISA is a nonempty closed convex subset of the vector space X, and the map J 0 := K 0 ISA ρ 0 a → s KL (ρ a |ρ 0 a ) ∈ R ∪ {+∞} is convex, and strictly convex on the convex set on which it takes finite values. The uniqueness of the solution to (30) is a direct consequence of this strict convexity property.

L-ISA counterpart of Lemma 5

For λ := (λ a ) 1≤a≤M ∈ R M , µ := (µ a ) 1≤a≤M ∈ R M , and ρ := (ρ a ) 1≤a≤M ∈ X M , we denote by λ µ ∈ R M and λ ρ ∈ X M the Kronecker products defined by

[λ µ] a = λ a µ a and [λ ρ] a = λ a ρ a . Lemma 11. Let ρ = (ρ a ) 1≤a∈M ∈ X M
+ such that for all 1 ≤ a ≤ M , there exists ρ 0 a ∈ K 0 z,L-ISA such that ´R3 ρ 0 a = ´R3 ρ a and s KL (ρ a |ρ 0 a ) < +∞. Then, the problem inf

ρ 0 ∈K 0 L-ISA , N (ρ 0 )=N (ρ) S(ρ|ρ 0 ) (80)
has a unique minimizer G 0

K 0 L-ISA (ρ). For all λ ∈ (R + ) M , it holds G 0 K 0 L-ISA (λ ρ) = λ G 0 K 0 L-ISA (ρ). (81)
Let us point out that problem (80) can be decomposed as M independent problems of the form inf

ρ 0 a ∈K 0 za,L-ISA , ´R3 ρ 0 a = ´R3 ρa s KL (ρ a |ρ 0 a )
for all 1 ≤ a ≤ M as in [START_REF] Hirshfeld | XVII. Spatial partitioning of charge density[END_REF].

Proof. The set {ρ 0 ∈ K 0 L-ISA | N (ρ 0 ) = N (ρ)} is non-empty and convex, and is included in a finite-dimensional vector subspace of X. The existence of a solution (80) follows from the compactness of K 0 ISA (for the strong topology of X M ) and the strong lower-semicontinuity of the functional J 0 introduced in the previous section. The second assertion of the lemma is a consequence of the fact that (i) for all λ ∈ R + M , ρ ∈ X M , and ρ 0 ∈ K 0 L-ISA , λ ρ 0 ∈ K 0 L-ISA , N (λ ρ 0 ) = N (λ ρ) and (ii) s KL (λf |λg) = λs KL (f |g) for all λ ∈ R + and f, g ∈ X + .

Topological properties of the set of admissible AIM decompositions

The following lemma provides some information about the topology of K ρ,R which will be useful for our analysis. Recall that M b (R 3 ) is the Banach space of the bounded (signed) Radon measures on R 3 endowed with the total variation norm and that M b (R 3 ) is the dual of the Banach space C 0 (R 3 ) of the real-valued continuous functions going to zero at infinity, and that a sequence

(µ n ) n∈N of elements of M b (R 3 ) weakly-* converge to some µ ∈ M b (R 3 ) if and only if ∀f ∈ C 0 (R 3 ), ˆR3 f dµ n -→ n→+∞ ˆR3 f dµ. Lemma 12. Let ρ ∈ X + and R = (R 1 , . . . , R M ) ∈ R 3M .
The set K ρ,R of admissible AIM decompositions defined by (3) is nonempty and convex. It is also bounded and closed for the norm topology of the Banach space X M , and sequentially compact in the following sense: from any sequence (ρ n ) n∈N = ((ρ n a ) 1≤a≤M ) n∈N of elements of K ρ,R , one can extract a subsequence (ρ n k ) k∈N converging toward an element ρ * of K ρ,R for the weak-* topologies of M b (R 3 ) M and L ∞ (R 3 ) M , and for the weak topology of L p (R 3 ) M for any 1 < p < ∞, and such that

∀1 ≤ a ≤ M, ˆR3 ρ n k a (r) dr -→ k→∞ ˆR3 ρ * a (r) dr, ( 82
) ˆR3 |r|ρ n k a (r) dr -→ k→∞ ˆR3 |r|ρ * a (r) dr. (83)
Proof. Let (ρ n ) n∈N denote a sequence of elements of K ρ,R so that for all n ∈ N, ρ n := (ρ n a ) n∈N .

Since for all n ∈ N and 1 ≤ a ≤ M , ρ n a ≥ 0 and ρ ≥ 0, the fact that

M a=1 ρ n a ( • -R a ) = ρ implies that ρ n a L 1 (R 3 ) ≤ ρ L 1 (R 3 ) and ρ n a L ∞ (R 3 ) ≤ ρ L ∞ (R 3 )
. Thus, up to the extraction of a subsequence, still denoted by (ρ n a ) n∈N for the sake of simplicity, there exists a bounded Radon measure ρ * a such that

ρ n a n→∞ ρ * a weakly-* in M b (R 3 ).
In particular, ρ * a ≥ 0 and

N a=1 ρ * a (• -R a ) = ρ. Since (1 + | • |)ρ ∈ L 1 (R 3 ), ρ * a is necessarily

Proof of Theorem 6

Using the arguments in Remark 4 for ISA, and assumption [START_REF] Singh | An approach to computing electrostatic charges for molecules[END_REF] for L-ISA, we see that the infimum in ( 50) is finite. Let ρ n , ρ 0,n n∈N be a minimizing sequence for [START_REF] Shi | Polarizable atomic multipole-based AMOEBA force field for proteins[END_REF]. By construction, the sequence ρ n , G 0 K 0 (ρ n ) n∈N , where

G 0 K 0 (ρ n ) := arg min ρ 0 ∈K 0 , N ( ρ 0 )=N (ρ n ) S(ρ n | ρ 0 )
(we use the notation ρ 0 to avoid a conflict of notation with the first term of the sequence (ρ n ) n∈N ), is uniquely defined for ISA and L-ISA by virtue of Lemmata 5 and 11, is also a minimizing sequence to [START_REF] Shi | Polarizable atomic multipole-based AMOEBA force field for proteins[END_REF]. We can therefore assume without loss of generality that ρ 0,n = G 0 K 0 (ρ n ).

We infer from Lemma 12 that there exists ρ * ∈ K ρ,R such that, up to extraction of a subsequence, (ρ n ) n∈N converges to ρ * for the weak-* topologies of M b (R 3 ) M and L ∞ (R 3 ) M , and for the weak topology of L p (R 3 ) M for any 1 < p < ∞, and

∀1 ≤ a ≤ M, N n a := ˆR3 ρ n a (r) dr -→ k→∞ ˆR3 ρ * a (r) dr =: N * a .
1 We want to show that G 0 K 0 (ρ n ) G 0 K 0 (ρ * ) for the weak-* topology of M b (R 3 ) M . i Let us first deal with ISA. We have ρ 0,n a := (G 0 K 0 ISA (ρ n )) a = ρ n a s (spherical average). We therefore have for all 1 ≤ a ≤ M and all n ∈ N, 0 ≤ ρ 0,n a ≤ ρ( • +R a ) s ∈ X + so that ρ 0,n a ∈ X + . Reasoning as in the proof of Lemma 12, we obtain the existence of a function ρ 0, * a ∈ X + such that ρ 0,n a n→+∞ ρ 0, * a , weakly-* in M b (R 3 ) and L ∞ (R 3 ), and weakly in L p (R 3 ),

for all 1 < p < ∞. The function G 0 K 0 ISA
is linear and strongly continuous on e.g. L 2 (R 3 ) M , hence weakly continuous on this space. We therefore have ρ 0,n = G 0

K 0 ISA (ρ n ) G 0 K 0 ISA (ρ * )
weakly in L 2 (R 3 ) M . Since (ρ 0,n ) n∈N also converges to ρ 0, * weakly in L 2 (R 3 ) M , we obtain that ρ 0, * = G 0

K 0 ISA (ρ * ). Thus, G 0 K 0 ISA (ρ n ) G 0 K 0 ISA (ρ * ) for the weak-* topology of M b (R 3 ) M .
ii Let us now turn to L-ISA. The above argument for ISA heavily relies on the fact that ρ → G 0

K 0 ISA (ρ) is a linear map, which is not the case of the map ρ → G 0 K 0 L-ISA (ρ).
We denote by K

0 za := ρ 0 a ∈ K 0 za,L-ISA | ´R3 ρ 0 a = 1 and F n a : K 0 za ρ 0 a → -´R3 ρ n a log ρ 0 a . Since K 0
za is a closed finite-dimensional simplex of the vector space (X r ) M , all the relevant topologies on (X r ) M (the strong and weak-* topologies of M b (R 3 ) M and L ∞ (R 3 ) M , the strong and weak topologies of L p (R 3 ) for 1 < p < ∞) are equivalent on K 0 za . It is then easily seen using the concavity of the logarithm that the functions F n a are convex and continuous, and strictly convex if ρ n a = 0. The same property holds for F * a . Let us prove that ρ 0,n a -→ n→∞ ρ 0, * a (for any of these topologies)

⇒ F n a (ρ 0,n a ) -→ n→∞ F * a (ρ 0, * a ). (84)
First, since for all n ∈ N, ρ 0,n a := mz a k=1 c n a,k g za,k where the g za,k are positive L 1 -normalized functions of X r + and mz a k=1 c n a,k = 1, it holds that g - a ≤ ρ 0,n a ≤ g + a so that log g - a ≤ log ρ 0,n a ≤ log g + a , where g - a := min 1≤k≤mz a g za,k and g + a := max 1≤k≤mz a g za,k . Thus, we obtain that

ρ n a | log ρ 0,n a | ≤ ρ(• + R a ) max | log g - a |, | log g + a | . Assumption (51) implies that ρ(• + R a )| log g - a | ∈ L 1 (R 3
). In addition, since g za,k ∈ X r + for all 1 ≤ k ≤ m za and are thus bounded and go to 0 at infinity, it holds that | log g + a (r)| ≤ | log g - a (r)| for sufficiently large values of |r|, which implies that ρ(• + R a )| log g + a | ∈ L 1 (R 3 ). Therefore, there exists a nonnegative function h a ∈ L 1 (R 3 ) such that for all n ∈ N, ρ n a | log ρ 0,n a | ≤ h a . Let > 0 and R > 0 such that ´R3 \B R h ≤ . Then,

|F n a (ρ 0,n a ) -F * a (ρ 0, * a )| ≤ 2 + ˆBR ρ n a log ρ 0,n a - ˆBR ρ * a log ρ 0, * a . Now, denoting by B R := {r ∈ R 3 | |r| < R}, we have that (ρ n a | B R ) n∈N weakly converges to ρ * a | B R in L 2 (B R )
, while (ρ 0,n a ) n∈N converges strongly to ρ 0, * a since all the relevant topologies are equivalent on K 0 za . The latter property also implies that (log ρ 0,n a ) n∈N converges strongly to log ρ 0, * a as we have g - a ≤ ρ 0,n a ≤ g + a with g + a bounded and g - a bounded away from zero on B R (this function is continuous and positive on the compact B R ). Hence (84). This properties implies in particular that the sequence of functionals (F n a ) n∈N Γ-converges to the functional F * a , and therefore that the minimizers of F n a converge to the minimizers of F * a , these minimizers being unique when ρ n a and ρ * a are not identically equal to zero. It follows from Lemma 11 that

[G 0 K 0 L-ISA (ρ n )] a = ρ 0,n a = N n a arg min K 0 za F n a if N n a = 0, 0 if N n a = 0, converges (in any relevant topologies) to ρ 0, * a = N * a arg min K 0 za F * a if N * a = 0, 0 if N * a = 0, = [G 0 K 0 L-ISA (ρ * )] a .
Hence the desired result.

2 Let us now prove that ρ * , G 0 K 0 (ρ * ) is a minimizer to problem [START_REF] Shi | Polarizable atomic multipole-based AMOEBA force field for proteins[END_REF]. To simplify the notation, let us set ρ 0, * = ρ 0, * a 1≤a≤M

:= G 0 K 0 (ρ * ). It is sufficient to prove that for all 1 ≤ a ≤ M , s KL ρ * a |ρ 0, * a ≤ lim inf n→+∞ s KL ρ n a |ρ 0,
→ R, ˆR3 f dρ n a -→ n→+∞ ˆR3 f dρ * a .
Besides, using similar arguments as above, it holds that the sequence ρ 0,n a n∈N with ρ 0,n a := . By the joint lower semicontinuity of the Kullback-Leibler divergence for the weak convergence of probability measures, it then holds that ˆR3 ρ * a log(ρ * a /ρ 0, * a ) ≤ lim inf n→+∞ ˆR3 ρ n a log(ρ n a /ρ 0,n a ).

ii Let us now assume that N * a = 0, which implies that ρ * a = ρ 0, * a = 0. Since N n a -→ n→+∞ N * a = 0, the sequences (ρ n a ) n∈N and (ρ 0,n a ) n∈N strongly converge to 0 in L 1 (R 3 ). Thus, up to the extraction of subsequences, both sequences converge almost everywhere to 0. Using the fact that • |x log x| ≤ |y log y| for all y ≤ 1/e and 0 ≤ x ≤ y;

• |x log x| ≤ e for all 0 ≤ x ≤ e,

• |x log x| ≤ |y log y| for all y ≥ e and 0 ≤ x ≤ y; and the bounds 0 ≤ ρ n a ≤ ρ(• + R a ) for all n ∈ N, we obtain that In the case when K 0 = K 0 L-ISA , since for all n ∈ N, ρ 0,n a := mz a k=1 c n a,k g za,k where the g za,k are positive L 1 -normalized functions of X r + , it holds that

|ρ n a log ρ n a | ≤ eχ 1 e ≤ρ(•+Ra)≤e + |ρ(• + R a ) log ρ(• + R a )| 1 -χ 1 e ≤ρ(•+Ra)≤e . In the case when K 0 = K 0 ISA , it holds that ρ 0,n a = ρ n a s for all n ∈ N, so that 0 ≤ ρ 0,n a ≤ ρ(• + R a ) s . Hence, |ρ 0,n a log(ρ 0,n a )| ≤ eχ 1 e ≤ ρ(•+Ra) s≤e + | ρ(• + R a ) s log ρ(• + R a ) s | 1 -χ 1 e ≤ ρ(•+Ra) s≤e , where χ E denotes the characteristic function of E. Note also that ρ(• + R a ) s log ρ(• + R a ) s ∈ L 1 (R 3
N n a g - a ≤ ρ 0,n a ≤ N n a g + a so that log g - a + log N n a ≤ log ρ 0,n a ≤ log g + a + log N n a . As a consequence, | log ρ 0,n a | ≤ | log N n a | + max | log g - a |, | log g + a | , and
ρ n a | log ρ 0,n a | ≤ ρ n a | log N n a | + ρ n a max | log g - a |, | log g + a | . On the one hand, ˆR3 ρ n a | log N n a | = N n a | log N n a | -→ n→+∞ 0.
On the other hand,

ρ n a max (| log g - a |, | log g + a |) ≤ ρ(• + R a ) max (| log g - a |, | log g + a |) ∈ L 1 (R 3 ). Thus, the Lebesgue dominated convergence theorem yields that ˆR3 ρ n a max | log g - a |, | log g + a | -→ n→+∞ 0. Finally, ´R3 ρ n a | log ρ 0,n a | -→ n→+∞ 0, and ˆR3 ρ n a log ρ n a ρ 0,n a -→ n→+∞ 0 = ˆR3 ρ * a log ρ * a ρ 0, * a .
We have thus proved that in any case,

ˆR3 ρ * a log ρ * a ρ * a s ≤ lim inf n→+∞ ˆR3 ρ n a log ρ n a ρ n a s
for all 1 ≤ a ≤ M , which implies that ρ * , ρ 0, * is a minimizer to (50).

3 The uniqueness comes from the strict convexity of the functional

K ρ,R × K 0 ρ, ρ 0 → S(ρ|ρ 0 ) ∈ R + ∪ {+∞}
on the convex subset of K ρ,R × K 0 on which it is finite.

Lastly, ( 52) is a direct consequence of the fact that ρ * := argmin ρ∈K ρ,R S ρ|ρ 0, * and Lemma 2.

Proof of Theorem 7

It is easy to check that the assumptions ρ > 0 and s KL (ρ|ρ 0,(0) ) < +∞ ensure that for both ISA and L-ISA, the first iteration is well defined. Assume now by induction that the first m iterations are well-defined for some m ∈ N * and that s KL (ρ|ρ 0,(m-1) ) < +∞ for all 1 ≤ m ≤ m.

For all 1 ≤ m ≤ m and 1 ≤ a ≤ M , we denote by

N (m) a := ˆR3 ρ (m) a = ˆR3 ρ 0,(m) a and N (m) = (N (m) a ) 1≤a≤M ∈ R M + .
From the assumptions of Theorem 7, it can be easily proven by recursion that N (m) a

> 0 for all 1 ≤ m ≤ m and all 1 ≤ a ≤ M .

Using Lemmata 2 and 11, we obtain that for all 1 ≤ m ≤ m, ρ (m) and ρ 0,(m) are respectively minimizers of ρ (m) ∈ argmin ρ∈K ρ,R S(ρ|ρ 0,(m-1) )

and ρ 0,(m) ∈ argmin

ρ 0 ∈K 0 , N (ρ 0 )=N (m) S(ρ (m) |ρ 0 ). (86)
In order to prove (54), we first show that S(ρ (m) |ρ 0,(m-1) ) ≥ S(ρ (m) |ρ 0,(m) ).

It follows from Lemmata 2 and 11 that µ (m) ρ 0,(m) , with µ (m) :=

N (m-1) a N (m) a 1≤a≤M
, is the unique minimizer of inf

ρ 0 ∈K 0 , N (ρ 0 )=N (m-1) S(ρ (m) |ρ 0 ). (88)
We therefore have

S(ρ (m) |ρ 0,(m-1) ) ≥ S(ρ (m) |µ (m) ρ 0,(m) ) = M a=1 ˆR3 ρ (m) a log    ρ (m) a N (m-1) a N (m) a ρ 0,(m) a    = M a=1 ˆR3 ρ (m) a log ρ (m) a ρ 0,(m) a + ˆR3 ρ (m) a log N (m) a N (m-1) a = M a=1 ˆR3 ρ (m) a log ρ (m) a ρ 0,(m) a + N (m) a log N (m) a N (m-1) a , = S(ρ (m) |ρ 0,(m) ) + M a=1 N (m) a log N (m) a N (m-1) a . It thus remains to show that M a=1 N (m) a log N (m) a N (m-1) a ≥ 0.
We now use the fact that -log

x ≥ 1 -x (since log x ≤ x -1) to get M a=1 N (m) a log N (m) a N (m-1) a = M a=1 -N (m) a log N (m-1) a N (m) a ≥ M a=1 N (m) a 1 - N (m-1) a N (m) a = M a=1 (N (m) a -N (m-1) a ) = 0, since M a=1 N (m) a = M a=1 ˆR3 ρ (m) a = ˆR3 ρ = M a=1 ˆR3 ρ (m-1) a = M a=1 N (m-1) a .
Hence (87). Let us now show that

S(ρ (m-1) |ρ 0,(m-1) ) ≥ S(ρ (m) |ρ 0,(m-1) ) + 1 2 ρ L ∞ M a=1 ρ (m) a -ρ (m-1) a 2 L 2 . (89)
Note that, by definition of ρ (m) , S(ρ (m-1) |ρ 0,(m-1) ) ≥ S(ρ (m) |ρ 0,(m-1) ). Using a second-order Taylor expansion formula with integral remainder, it holds that for all x, x , y > 0,

x log x y ≥ x log x y + (x -x) log x y + 1 + 1 max(x, x ) (x -x ) 2 2 .
As a consequence, since for all 1 ≤ a ≤ M ,

1 ρ (m) a L ∞ ≥ 1 ρ L ∞ , we obtain that S(ρ (m-1) |ρ 0,(m-1) ) ≥ S(ρ (m) |ρ 0,(m-1) ) + M a=1 ˆR3 (ρ (m-1) a -ρ (m) a ) log ρ (m) a ρ 0,(m-1) a + 1 + M a=1 1 2 ρ L ∞ ˆR3 (ρ (m-1) a -ρ (m) a ) 2 .
Using the fact that the expression

ρ (m) a (r) ρ 0,(m-1) a (r) = ρ(r + R a ) M b=1 ρ 0,(m-1) b (r -R b + R a )
is independent of a and the equality

M a=1 ˆR3 ρ (m-1) a = M a=1 ˆR3 ρ (m) a = ˆR3 ρ,
we thus obtain that S(ρ (m-1) |ρ 0,(m-1) ) ≥ S(ρ (m) |ρ 0,(m-1) )

+ M a=1 ˆR3 (ρ (m-1) a (r -R a ) -ρ (m) a (r -R a )) log ρ(r) M b=1 ρ 0,(m-1) b (r -R b ) dr + 1 2 ρ L ∞ M a=1 ρ (m-1) a -ρ (m) a 2 L 2 , = S(ρ (m) |ρ 0,(m-1) ) + ˆR3 (ρ(r) -ρ(r)) log ρ(r) M b=1 ρ 0,(m-1) b (r -R b ) dr + 1 2 ρ L ∞ M a=1 ρ (m-1) a -ρ (m) a 2 L 2 , = S(ρ (m) |ρ 0,(m-1) ) + 1 2 ρ L ∞ M a=1 ρ (m-1) a -ρ (m) a 2 L 2 .
Hence (89). Collecting this result together with (87), we obtain [START_REF] Stone | Distributed multipole analysis[END_REF]. It also follows from ( 17), (85), and (87) that s KL (ρ, ρ 0,(m-1) ) = S(ρ (m) |ρ 0,(m-1) ) ≥ S(ρ (m) |ρ 0,(m) ) ≥ inf ρ∈K ρ,R S(ρ|ρ 0,(m) ) = s KL (ρ, ρ 0,(m) ).

In particular s KL (ρ, ρ 0,(m) ) < +∞ for all 1 ≤ m ≤ m, so that the iteration m + 1 is well defined. By recursion, all the iterations are well defined and the above results are valid for all m ≥ 1. The bound ( 55) is obtained by a simple combination of the above inequalities.

Another consequence of ( 54) is that S(ρ (m) |ρ 0,(m) ) m≥0 is a non-increasing sequence of nonnegative real numbers, hence converges to a nonnegative real number as m tends to infinity. The bound (54) also implies that Let us finally prove the last assertion of Theorem 7. We infer from Lemma 12 that there exist ρ * ∈ K ρ,R and a subsequence (ρ (m l ) ) l∈N of (ρ (m) ) m∈N converging to ρ * for the weak-* topologies of M b (R 3 ) M and L ∞ (R 3 ) M , and for the weak topology of L p (R 3 ) M for any 1 < p < ∞, and such that ∀1 ≤ a ≤ M, N (m l ) a :=

ˆR3

ρ (m l ) a (r) dr -→ l→∞ ˆR3 ρ * a (r) dr =: N * a .

In addition, we can prove that there exists ρ 0, * ∈ K 0 , such that, up to extraction of a subsequence, (ρ 0,(m l ) ) l∈N converges to ρ 0, * for the weak-* topologies of M b (R 3 ) M and L ∞ (R 3 ) M , and for the weak topology of L p (R 3 ) M for any 1 < p < ∞. Using similar arguments as in the proof of Theorem 6, we obtain that ρ 0, * = G 0 K 0 L-ISA (ρ * ).

Moreover, since K 0 L-ISA ⊂ Span {g za,k , 1 ≤ a ≤ M, 1 ≤ k ≤ m za }, which is a finite-dimensional vector space, the convergence of the sequence (ρ 0,(m l ) ) l∈N to ρ 0, * holds strongly in any L p (R 3 ) for any 1 ≤ p ≤ +∞.

Using the fact that This relationship, together with the fact that ρ 0, * = G 0 K 0 L-ISA (ρ * ) implies that (ρ * , ρ 0, * ) is the unique minimizer of (50). Hence the desired result.

ρ (m l +1)

Proof of Proposition 8

If ρ(• + R 2 ) = ρ(• + R 2 ) s , it is clear that ρ opt 1 = 0 and ρ opt 2 = ρ opt s . Conversely, let us assume that ρ opt 1 = 0 (and therefore ρ opt 2 = ρ(• + R 2 ) > 0 on R 3 , which implies that ρ opt 2 s = ρ(• + R 2 ) s > 0 on R 3 ), and prove that necessarily ρ opt 2 = ρ(• + R 2 ) s . For all ρ ∈ K ρ,R , we denote by J (ρ) := S(ρ| ρ s ). Let us consider radially symmetric functions h 1 ∈ X r and h 2 ∈ X r such that ρ opt 2 s + h 2 > 0 a.e. for > 0 small enough, and consider the perturbed AIM densities To the first order in , we have, uniformly in L ∞ loc (R 3 ),

ρ 1 = h 1 ρ(• + R 2 ) s (• + R 1 ) ρ(• + R 1 ) + o( ), ρ 2 = ρ(• + R 2 ) - h 1 (• -R 1 + R 2 ) ρ(• + R 2 ) s ρ(• + R 2 ) + o( ),
and thus

J (ρ ) -J (ρ opt ) = - ˆR3 ρ(• + R 2 ) ρ(• + R 2 ) s log ρ(• + R 2 ) ρ(• + R 2 ) s h 1 (• -R 1 + R 2 ) + o( ).
Let us denote by ρ s,2 (r) := ρ(• + R 2 ) s (r -R 2 ) .

It then holds that Choosing a particular function h 1 such that h 1 (rσ) = g 1 (r) with g 1 > 0 on A and h 1 = 0 on [0, +∞) \ A would imply that ˆR3 ρ(• + R 1 ) ρ s,2 (• + R 1 ) log ρ(• + R 1 ) ρ s,2 (• + R 1 )

J (ρ ) -J (ρ opt ) = - ˆR3 ρ(• + R 1 ) ρ s,2 (• + R 1 ) log ρ(• + R 1 ) ρ s,2 (• + R 1 ) h 1 + o( ) = - ˆR3 ρ(• + R 1 ) ρ s,2
= ˆR3 ρ(• + R 2 ) ρ(• + R 2 ) s log ρ(• + R 2 ) ρ(• + R 2 ) s = ˆR3 ρ(• + R 2 ) ρ(• + R 2 ) s log ρ(• + R 2 ) ρ(• + R 2 ) s s ≥ ˆR3 ρ(• + R 2 ) s ρ(• + R 2 ) s log ρ(• + R 2 ) s ρ(• + R 2 ) s = 0,
s h 1 = 4π ˆ+∞ 0 w(r)g 1 (r)r 2 dr > 0.

Letting go to 0, we reach a contradiction since (ρ opt , ρ opt s ) is a minimizer of (50). Hence, it necessarily holds that

ˆR3 ρ(• + R 1 ) ρ s,2 (• + R 1 ) log ρ(• + R 1 ) ρ s,2 (• + R 1 ) s = ˆR3 ρ(• + R 2 ) ρ(• + R 2 ) s log ρ(• + R 2 ) ρ(• + R 2 ) s s = 0. Since ρ(•+R 2 ) ρ(•+R 2 ) s log ρ(•+R 2 )
ρ(•+R 2 ) s s ≥ 0 almost everywhere from Jensen's inequality, we obtain

ρ(• + R 2 ) ρ(• + R 2 ) s log ρ(• + R 2 ) ρ(• + R 2 ) s s = 0.
This equality implies that

ρ(• + R 2 ) log ρ(• + R 2 ) s = ρ(• + R 2 ) s log ρ(• + R 2 ) s a.e.
Using again Jensen's inequality and the strict convexity of the function x → x log x, we finally obtain ρ(• + R 2 ) = ρ(• + R 2 ) s . Hence the desired result.

Proof of Proposition 9

It is clear from their definition that the functions w a are nonnegative and bounded by M ρ . Taking spherical averages in the optimality condition [START_REF] Stone | Distributed multipole analysis, or how to describe a molecular charge distribution[END_REF], and using (53), we see that for all 1 ≤ a ≤ 2, I a (r) = Note that the function S 2 σ → ρ(R a + rσ) is continuous on S 2 , nonnegative, and bounded by M ρ . Consequently, F a is well-defined on (0, +∞) × (0, +∞), and for each r > 0, the function (0, +∞) y → F a (r, y) ∈ R is decreasing and converges to -1 at infinity, so that for each r ∈ (0, +∞), there exists at most one y ∈ R + such that F a (r, y) = 0. We denote by w + a (r) the unique solution in R + to F a (r, w + a (r)) = 0 if it exists, and we set w + a (r) = 0 otherwise. We therefore have w a (r) ∈ {0, w + a (r)} for almost all r ∈ R + . We have for all r, y, z ∈ (0, +∞), dσ.

Using the fact that ρ is continuous and positive everywhere on R 3 , hence bounded away from zero on any compact subset of R 3 , we get that for all (r, y, z) ∈ (0, +∞) × [0, M ρ ] The function w 1 is therefore defined implicitly on (0, +∞) by F 1 (r 1 , w 1 (r 1 )) = 0 or w 1 (r 1 ) = 0, where the function F 1 : (0, +∞) × R + → [0, +∞] is given by

F 1 (r 1 , y) = ˆR+r 1 |R-r 1 | h 1 (r 1 , r 2 ) y + w 2 (r 2 ) dr 2 -1.
Note that the function h 1 is continuous on D 1 and that we have for all (r 1 , r 2 ) ∈ D 1 ,

|h 1 (r 1 , r 2 )| ≤ M ρ r 2 2r 1 R .
By dominated convergence, F 1 is therefore locally bounded and continuous on (0, +∞)×(0, +∞).

Recall that 0 ≤ w a (r a ) ≤ M ρ for a.a. r a ∈ (0, +∞). We therefore have for all r 1 , r 1 , y, z ∈ (0, +∞), F 1 (r 1 , y) -F 1 (r 1 , z) ≤ F 1 (r 1 , y) -F 1 (r 1 , y) + F 1 (r 1 , y) -F 1 (r 1 , z) Assume that we can extract from (r 1,n ) n∈N a subsequence (r 1,n k ) k∈N with positive limit, i.e., such that w + 1 (r 1,n k ) → y * > 0 as k tends to infinity. Passing to the limit by the dominated convergence theorem, we get

  do not need to be stored in memory. The integrals ´R3 ρ (m) a and the entries of the vectors b (m) a can be computed using quadrature formulas from the vector c (m-1) a

  k equals the atomic number z a , are possible, such as∀k ∈ [[1, m za ]], c

Figure 1

 1 represents the numerically computed r → log 4πr 2 ρ 0 a (r) profiles solutions of GISA, L-ISA, MB-ISA and ISA for the carbon atom in the CO molecule, using a molecular density ρ computed at the RHF/aug-cc-pVDZ level of theory. Qualitatively similar profiles are obtained for the oxygen atom. Note that the quantity 4πr 2 ρ 0 a (r) represents the contribution of the atomic density ρ a (whose spherical average equals ρ 0 a ) to the atomic charge N a = ´+∞ 0 4πr 2 ρ 0 a (r)dr at a radial distance r.

Figure 1 :

 1 Figure 1: Profiles r -→ log 4πr 2 ρ 0 a (r) for the carbon atom in CO (1.131 Å interatomic distance, RHF/aug-cc-pVDZ level of theory), obtained by GISA, L-ISA, and MB-ISA and ISA methods. N r = 320 Gauss-Legendre radial (up to R max =8 Å) and N s = 150 Gauss-Legendre angular points were used for both the quadratures and the ISA discretization. The black vertical dashed line indicates the position of the oxygen atom.

Figure 2 :

 2 Figure 2: Profiles r -→ log 4πr 2 ρ 0 a (r) for the chlorine atom in ClO -ion (1.733 Å interatomic distance, RHF/6-31G(d) level of theory), obtained by GISA, L-ISA, and MB-ISA and ISA methods. N r = 500 Gauss-Legendre radial (up to R max =8 Å) and N s = 250 Gauss-Legendre angular points were used for both the quadratures and the ISA discretization. The black vertical dashed line indicates the position of the oxygen atom.

Figure 3 :

 3 Figure 3: Profiles r -→ log 4πr 2 ρ 0 a (r) for the oxygen atom of ClO -ion (1.733 Å interatomic distance, RHF/6-31G(d) level of theory), obtained by GISA, L-ISA, and MB-ISA and ISA methods. N r = 500 Gauss-Legendre radial (up to R max =8 Å) and N s = 250 Gauss-Legendre angular points were used for both the quadratures and the ISA discretization. The black vertical dashed line indicates the position of the chlorine atom.

Figure 4 :

 4 Figure 4: Quantities of interest for LiH as functions of the interatomic distance, with a density obtained at the FCI/6-31G level of theory. Top-left : ground state energy (Ha). Top-right : ISA atomic charge on Li. Bottom: ISA atomic dipoles along the dimer axis on Li (left panel) and H (right panel).

Figure 5 :

 5 Figure 5: Local atomic charges on hydrogen and oxygen atoms of the water molecule with increasing basis set size, computed from DMA (Stone's and Vigné-Maeder's redistribution strategies) and ISA. The number of basis functions refers to the whole molecule (two hydrogen atoms and oxygen atoms).

Figure 6 :

 6 Figure 6: Local atomic dipoles on hydrogen and oxygen atoms of the water molecule with increasing basis set size, computed from DMA (Stone's and Vigné-Maeder's redistribution strategies) and ISA. The number of basis functions refers to the whole molecule (two hydrogen atoms and oxygen atoms).

  absolutely continuous with respect to the Lebesgue measure and hence(1 + | • |)ρ * a ∈ L 1 (R 3). Besides, (ρ n a ) n∈N being bounded in every L p (R 3 ) for 1 < p ≤ ∞, it holds that, up to the extraction of a subsequence,ρ n a n→+∞ ρ * a weakly-* in L ∞ (R 3) and weakly in L p (R 3 ), for all 1 < p < ∞.Thus ρ * = (ρ * a ) 1≤a≤M ∈ K ρ,R , and the sequential compactness properties listed in the lemma hold true. Moreover, let us denote by N n a := ´R3 ρ n a and by N * a := ´R3 ρ * a . It holds that N * a ≤ lim inf n→+∞ N n a . The fact that for all n ∈ N * , M a=1 N n a = ´R3 ρ = M a=1 N * a implies (82). The convergence (83) is obtained in the same way.

  the sense of probability measures to ρ 0, * a := ρ 0, * a N * a

  ) by the convexity of the function R + x → x log x. Thus, since ˆR3 ρ

2 L 2 <

 22 +∞, which, in turn, implies[START_REF] Sun | Recent developments in the PySCF program package[END_REF].

  m l ) b (r -R b + R a ) ρ(r + R a )together with[START_REF] Sun | Recent developments in the PySCF program package[END_REF] yields thatρ * a (r) = ρ 0, * a (r) b=1 ρ 0, * b (r -R b + R a ) ρ(r + R a ).

  b )(• -R b + R a ) ρ(• + R a ).

0 .

 0 the last inequality being a consequence of Jensen's inequality. Let us reason by contradiction and assume thatˆR3 ρ(• + R 1 ) ρ s,2 (• + R 1 ) log ρ(• + R 1 ) ρ s,2 (• + R 1 ) s >There would then exist a bounded Borel subset A of [0, +∞) with positive Lebesgue measure such that w(r) > 0 for almost all r ∈ A where w ∈ Y + is the function such that ∀r > 0, ∀σ ∈ S 2 , w(r):= ρ(• + R 1 ) ρ s,2 (• + R 1 ) log ρ(• + R 1 ) ρ s,2 (• + R 1 )s (rσ).

S 2 ρ 2 ρ

 22 (R a + rσ) w a (r) + w b (|R a + rσ -R b |) dσ = 1, or w a (r) = 0, where b = 2 if a = 1 and b = 1 if a = 2.For given and fixed w b , the function w a is therefore defined implicitly on (0, +∞) by for almost all r ∈ (0, +∞), F a (r, w a (r)) = 0 or w a (r) = 0,where the function F a : (0, +∞) × R + → [0, +∞] is given byF a (r, y) = S (R a + rσ) y + w b (|R a + rσ -R b |) dσ -1.

2 ρ

 2 |F a (r, y) -F a (r, z)| = S 2 ρ(R a + rσ) 1 y + w b (|R a + rσ -R b |) -1 z + w b (|R a + rσ -R b |) dσ = K a (r, y, z) |y -z|,(90)withK a (r, y, z) := S (R a + rσ) (y + w b (|R a + rσ -R b |))(z + w b (|R a + rσ -R b |))

2 ,For r 1 , r 2 ∈ 1 + R 2 - 2 2 and r 2 2 + R 2 - 1 ,4π ˆπ 0 dθ ˆ2π 0 dφ ρ 1 1 (r 1 , θ 1 (r 1 , r 2 ), φ) w 1 (r 1 ) + w 2 (r 2 ) r 2 r 1 R = ˆR+r 1 |R-r 1 | h 1 (r 1 , r 2 ) w 1 (r 1 ) + w 2 (r 2 ) dr 2 ,r 2 4πr 1 R ˆ2π 0 ρ 1

 2121222210111112112111121122201 K a (r, y, z) ≥ m a , positive, decreasing function of r. R + such thatr 1 + r 2 ≥ R = |R a -R b |, we denote by θ 1 (r 1 , r 2 ) and θ 2 (r 1 , r 2 ), the unique real numbers in [0, π] such that r 2 2r 1 R cos(θ 1 (r 1 , r 2 )) = r 2r 2 R cos(θ 2 (r 1 , r 2 )) = r 2and by ρ a (r, θ, φ) the representation of the total density ρ in spherical coordinates relative to R a , with the z-axis aligned with e := R b -Ra R , i.e. ρ a (r, θ, φ) := ρ(R a + r(sin θ cos φ e x + sin θ sin φ e y + cos θ e)),where (e x , e y , e) forms an orthonormal basis of R 3 . Using the change of variables r 2 = (r 2 1 + R 2 -2r 1 R cos θ) 1/2 , for which r 1 R sin θ dθ = r 2 dr 2 , we haveI 1 (r 1 ) = 1 (r 1 , θ, φ) w 1 (r 1 ) + w 2 ((r 2 1 + R 2 -2r 1 R cos θ)where the function h 1 : D 1 → R + is defined byD 1 := {(r 1 , r 2 ), r 1 ∈ (0, +∞), r 2 ∈ [|R-r 1 |, R+r 1 ]}, h 1 (r 1 , r 2 ) := (r 1 , θ 1 (r 1 , r 2 ), φ) dφ.

≤

  2M ρ min(y, z) |r 1 -r 1 | + M ρ yz |y -z|, which proves that F 1 is locally Lipschitz on the open set (0, +∞) × (0, +∞). It follows from (90) and the implicit function theorem for Lipschitz functions (see Theorem 13 in Appendix B) that if w + 1 (r 1 *) > 0 for some r 1 * > 0, then w + 1 is Lipschitz in the neighborhood of r 1 * . Let r 1 * ∈ (0, +∞) be such that w + 1 (r 1 * ) = 0, and (r 1,n ) ∈ (0, +∞) N a sequence converging to r 1 * . We have for all n ∈ NF 1 (r 1,n , w + 1 (r 1,n )) = ˆR+r 1,n |R-r 1,n | h 1 (r 1,n , r 2 ) w + 1 (r 1,n ) + w 2 (r 2 )dr 2 -1 = 0.

F 1 (

 1 r 1 * , y * ) = ˆR+r 1 * |R-r 1 * | h 1 (r 1 * , r 2 ) y * + w 2 (r 2 )dr 2 -1 = 0.

Table 1 :

 1 for atom O in the same ion).

	Multipole moment (atomic units)	GISA	L-ISA	MB-ISA	ISA
	Charge q 1 = q C	0.129	0.182	0.225	0.183
	Dipole d 1 z	0.217	0.298	0.325	0.307
	Dipole d 2 z	-0.015	0.015	0.079	0.008
	Q 1 xx = Q 1 yy	4.109	3.903	3.835	3.887
	Q 2 xx = Q 2 yy	3.476	3.682	3.751	3.696
	Q 1 zz	5.100	4.823	4.789	4.786
	Q 2 zz	3.636	3.806	3.917	3.812

Local charges, dipoles and second-order moments carried by C (atom 1) and O (atom 2) in CO (1.131 Å interatomic distance, RHF/aug-cc-pVDZ level of theory) computed by GISA, L-ISA, MB-ISA and ISA methods. 320 Gauss-Legendre radial (up to R max =8 Å) and 150 Gauss-Legendre angular points were used for both the quadratures and for the ISA discretization.

Table 2 :

 2 Local charges, dipoles and second-order moments carried by Cl (atom 1) and O (atom 2) in ClO -(1.733 Å interatomic distance, RHF/6-31G(d) level of theory), obtained by GISA, L-ISA, and MB-ISA and ISA methods. N r = 500 Gauss-Legendre radial (up to R max =8 Å)

	Multipole moment (atomic units)	GISA	L-ISA	MB-ISA	ISA
	Charge q 1 = q Cl	1.941	-0.368	-0.328	-0.342
	Dipole d 1 z	-1.024	-0.117	-0.198	-0.232
	Dipole d 2 z	-6.569	0.086	0.033	0.115
	Q 1 xx = Q 1 yy	4.455	10.460	10.333	10.353
	Q 2 xx = Q 2 yy	10.581	4.576	4.703	4.682
	Q 1 zz	4.176	9.787	9.567	9.364
	Q 2 zz	28.460	4.026	4.159	3.983

  (• + R 1 ) log ρ(• + R 1 ) ρ s,2 (• + R 1 )

	First, we have ˆR3 ρ(• + R 1 ) ρ s,2 (• + R 1 )	log	ρ(• + R 1 ) ρ s,2 (• + R 1 )	s	= = =	ˆR3 ˆR3 ˆR3	ρ(• + R 1 ) ρ s,2 (• + R 1 ) ρ ρ ρ s,2 log ρ s,2 log ρ ρ s,2 ρ s,2 log ρ	ρ(• + R 1 ) ρ s,2 (• + R 1 ) s,2

s h 1 + o( )

If the ground state is degenerate and the radial symmetry broken, ρ 0 z,n is chosen equal to the radially symmetric mixed-state ground-state density obtained by averaging the pure-state ground-state densities (with respect to the Haar measure of the rotation group SO(3)).

Ackowledgements R.B. thanks Michael Herbst, Emmanuel Giner and Laurent Vidal for useful discussions and technical help. The authors are also very grateful to Antoine Levitt and Alston Misquitta for stimulating discussions. This publication is part of a project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Programme (Grant Agreement n • 810367).

Hence, y * = w + 1 (r 1 * ) = 0. We reach a contradiction. This proves that w + 1 is continuous on (0, +∞). The same holds for w + 2 by symmetry. This completes the proof of the first statement. Let us now establish the second one. If the w a 's are bounded away from zero on every compact subset of [0, +∞), then w a = w + a almost everywhere on [0, +∞), which, in view of the first statement, implies that the w a 's are Lipschitz on (0, +∞) and that for all r 1 , r 2 > 0,

where

By a simple continuity argument, we deduce from the fact that the w j 's are bounded away from zero on [0, 1] that w j is Lipschitz on the whole interval [0, +∞). If ρ is C 1 away from the centers R j , then the functions F j : [0, +∞) × (0, +∞) → R are C 1 on ((0, R) ∪ (R, +∞)) × (0, +∞). By the implicit function theorem for C 1 functions, the functions w j are C 1 on (0, R) ∪ (R, +∞) (we use (90) and its analogue for w 2 to show that the assumptions of the implicit function theorem are satisfied). By a simple bootstrap argument, if ρ is C k away from the centers R j , then the functions w j are C k on (0, R) ∪ (R, +∞).

Appendix A: Non-uniqueness of GISA fixed-points

Since GISA is covered by the framework [START_REF] Shi | Polarizable atomic multipole-based AMOEBA force field for proteins[END_REF], uniqueness of a minimizer, as stated in Theorem 6 for ISA and L-ISA, can not be guaranteed. We have therefore investigated the uniqueness of minimizers and have observed numerically a sensitivity of the GISA solution (obtained by the fixed-point iterations described in Section 2.6) to the initial guess that we report in this Appendix. More specifically, for a (non-dimensional) test density sum of two normalized Gaussian functions:

131, using six shells on each atom in the GISA pro-atomic densities variational space, and Gaussian exponents α 1,k ∈ {0.01, 0.1, 1, 2, 5, 10} , α 2,k ∈ {0.05, 0.5, 2, 4, 10, 50} , we found two local minima ρ 0, * a and ρ 0, * * a , whose profiles r → log 4πr 2 ρ 0 a (r) are represented in Figure 7. The expected solution, corresponding to one unit charge on each atom was obtained e.g. with the initial guess:

= (0, 0, 0, 1, 0, 0), c (0), * 2,k k=1..6

= (0, 0, 0, 0, 1, 0),

(i.e. with weights initially on more compact Gaussians than the expected solutions), while a more balanced initial guess:

leads to a fixed-point associated to slightly different pro-atomic densities (see Figure 7, right panel) and slightly distorted atomic charges and dipoles, see Table 3. In particular, F (x, φ(x)) = 0 for all x ∈ U .