Benjamin Stamm

Louis Lagardère

Étienne Polack

Yvon Maday

Jean-Philip Piquemal

A coherent derivation of the Ewald summation for arbitrary orders of multipoles: The self-terms

In this work, we provide the mathematical elements we think essential for a proper understanding of the calculus of the electrostatic energy of point-multipoles of arbitrary order under periodic boundary conditions. The emphasis is put on the expressions of the so-called self parts of the Ewald summation where different expressions can be found in literature. Indeed, such expressions are of prime importance in the context of new generation polarizable force field where the self field appears in the polarization equations. We provide a general framework, where the idea of the Ewald splitting is applied to the electric potential and subsequently, all other quantities such as the electric field, the energy and the forces are derived consistently thereof. Mathematical well-posedness is shown for all these contributions for any order of multipolar distribution.

INTRODUCTION

The computation of physical quantities involving the Coulomb potential is a challenging issue due to the slow decay of the interacting kernel as the inverse of the distance. This long-range potential often prevents the use of simple techniques like cutoffs methods that only take into account short-range interactions. This problem has been addressed with the use of hierarchical methods (of order O(N ) or O(N log N ) complexity) that approximate the long-range interactions and Fourier (of order O(N log N )) methods that compute part of the Coulomb interaction in the dual space by considering the physical system under periodic boundary conditions.

For molecular dynamics simulations of biological systems, the most widely used method is a Fourier method, the particle-mesh Ewald 1,2 -or shortly pme. This method is based on the Ewald summation 3 , which gives a well-posed definition for the energy of the system. This is indeed not granted at all, since the energy is not well defined due to the conditional convergence of the involved series of the infinite periodic system if the (neutral) unit cell has a non-zero dipolar moment. In this case, different orders of summation provide different energies.

Background on the Ewald summation. The mathematical derivation of the Ewald energy summation for point charges in three dimensions was carried out by Redlack and Grindlay 4 , de Leeuw, Perram, and Smith 5 . With respect to the focus of this paper involving multipoles of any order, Weenk and Harwig 6 and Smith 7 gave expressions for the energy using Ewald summation for density of charges expressed as a sum of multipoles up to quadrupoles. Those expressions have been used, for example, in the works of of Nymand and Linse 8 , Toukmaji et al. 9 , Wang and Skeel 10 for dipoles and by Aguado and Madden 11 for quadrupoles.

However some expressions in the paper by Smith 7 are justified using physical insight, and only the Ewald energies and forces are given. We think this is the reason why some other authors use other (inconsistent) expressions. For example Nymand and Linse 8 give an expression for the electric field that is different from the one by Toukmaji et al. 9 . This difference was then discussed by Laino and Hutter 12 and corrected in Stenhammar, Trulsson, and Linse 13 .

Moreover all the terms (potential, field, energy, forces) for the Ewald summation are in our knowledge never presented all together in one place consistently, and the derivation is seldom explained. For example, Wang and Skeel 10 , Aguado and Madden 11 don't give expression for the field and Toukmaji et al. 9 give no expression for the potential. Stenhammar, Trulsson, and Linse 13 builds an exception, however, the proposed self-energy differs for quadrupolar distributions as the work by Stenhammar, Trulsson, and Linse 13 does not include a quadruple-quadrupol interaction whereas Aguado and Madden 11 does. The latter is however with a different formula than what we propose later in this work. This may be explained by a missing double factorial in Aguado and Madden 11 and Nymand and Linse 8 that was pointed out by Laino and Hutter 12 . As only the net expressions are provided, it is difficult to trace back this difference. Recent developments have been made for efficient PME calculations using spherical harmonic point multipoles in Giese et al. 14 and Simmonett et al. 15 , where in particular the former also provides expressions for energies, potentials, and forces using arbitrary order point multipoles.

Contribution. This paper should be seen as an extension of the work of Smith 7 . Although not fully rigorous and lead by physical intuition, his reasoning for the expression of the self-energy can be proven with the use of some mathematical arguments, which can then be used to find the self-terms of any multipolar distribution. While we do not introduce a new theory, model or mathematical expressions, we introduce here a coherent mathematical framework to derive the self-terms of multipolar distributions of any order for the electric potential and field as well as the associated energy and forces and confirm the results proposed by Smith 7 . Further, we present proofs of the well-posedness of the self parts to the energy, electric potential and field for multipolar distributions of any order.

Our derivation is different from what has been proposed in the past, and emphasizes that the Ewald splitting should first be done on the potential or the field -and not directly on the energy. We derive the self-potential and self-field from scratch using Ewald splitting and deduce from those expressions the results for the self-energy and self-forces.

The purpose of the present article is to provide a coherent mathematically driven derivation of all self-terms, which, in consequence, provides a base for methodology developments of force-fields. We present in the appendix of a complete and precise derivation of all self-terms such that differences in expressions as highlighted above can be traced back. In particular, and this is our main motivation, a correct derivation of the self-field is indispensable for polarizable force-fields. Indeed, to solve the polarization equation, the total field, and thus the self-field, is required to compute the polarization field 16 . In practice, such terms are well implemented in production codes like Tinker and Amber. However, other codes exist and omitting these terms would result in highly different properties. Indeed, as shown in Fig- Outline. First, in Section I we introduce the notations that we use and review general results about the Ewald summation. In Section II, Section III and Section IV we give, respectively, a derivation for the potential, the field and the energy using Ewald summation.

Finally, in Section V we give explicit expressions of the self-terms and provide the proof that justifies the existence of the self-terms for any multipolar distribution.

I. EWALD SUMMATION FOR MULTIPOLES

In this article, we consider a system composed of a discrete distribution of N point multipoles in R 3 under periodic boundary conditions. The system consists in an electronically neutral primitive triclinic cell U with charges in form of multipoles located at r i ∈ U for i ∈ {1, . . . , N }. The set of positions r i is represented by the global vector r {N } := (r 1 , . . . , r N ).

The unit cell U is then duplicated in all directions and the system derived from r {N } is therefore composed of an infinity of charges.

The unit cell U is spanned by the three vectors (a 1 , a 2 , a 3 ) which is called the basis of U .

We then introduce the lattice-indices n and m of the form

n = 1≤γ≤3 n γ a γ and m = 1≤γ≤3 m γ a * γ , (1)
where n γ , m γ ∈ N and (a * 1 , a * 2 , a * 3 ) is the dual basis of (a 1 , a 2 , a 3 ): that is a * γ • a γ = δ γγ (the Kronecker symbol). We will also denote by V the volume of the primitive cell U and by U * the dual of the primitive cell.

Then, one can informally introduce "the" electrostatic interaction energy of r {N } up to

2 p -poles, p ∈ N as E(r {N } ) := 1 2 n 1≤i,j≤N L i L j 1 |r ijn | , (2)
where r ijn := r ir j + n, the sign on the sum means that for i = j when n = 0 the interaction is not counted (this avoids self-interaction of a point multipole with itself) and the multipolar operator L i is defined as

L i := 0≤k≤p M k i • D k i . (3)
Here,

M k i is a k-dimensional array of dimension 3 k describing moment of the point 2 k -pole, D k
i is the matrix of k-order partial derivatives with respect to the variable r i and • is the point-wise product which writes

(A • B) α = A α B α = A α 1 ,...,α k B α 1 ,...,α k , for two arbitrary k-dimensional arrays A, B ∈ R 3k and where α = (α 1 , . . . , α k ), α i ∈ {1, 2, 3}, is a k-dimensional multi-index.
For instance, k = 0 represents a point-charge of charge M 0 i at r i and k = 1 a dipole where D 1

i is equivalent to the usual ∇ notation with respect to r i and M 1 i denotes the dipolar moment for each location r i . Next, k = 2 represents a quadrupole, D 2 i denotes the Hessian matrix and M 2 i is a 3 × 3 matrix that incorporates the quadrupolar moments. We will see that the energy in (2) is actually not well defined: As in the case of single point charges, it can be shown, by a Taylor expansion with respect to n, that the series in equation ( 2) is what is called conditionally convergent. That implies that the result of the energy E(r {N } ) depends on the order of summation and is thus not uniquely defined.

The electrostatic energy can equivalently be stated in the following form

E(r {N } ) = 1 2 1≤i≤N L i φ i (r i ) with φ i (r i ) := n 1≤j≤N L j 1 |r ijn | , (4)
so that φ i (r i ) denotes the potential at r i which is generated by all multipoles different than the one located at r i . In consequence, equation ( 2) represents indeed the interaction energy between every multipole i in the unit cell with the potential created by all other multipoles (indexed by j and n) of the infinite lattice.

Let us make a subtle comment. While r i is the fixed position of the i-th multipole, the multipole operator L i involves derivatives which requires to consider the potential φ i in a local neighborhood of r i . We denote therefore by r the variable belonging to a local neighborhood of r i and write

L i φ i (r i ) = M k i • D k r φ i (r) r=r i , (5)
since we have to consider the potential φ i (r) and its derivatives ultimately evaluated at r = r i .

As anticipated above, equations ( 2) and ( 4) are not well-defined and hence the need to use a definition of an expression for the energy that is well defined. One possible remedy is the introduction of the Ewald energy to give a unique meaning of this expression by

E Ewald (r {N } ) := 1 2 n 1≤i,j≤N L i L j erfc(α 1/2 |r ijn |) |r ijn | + 1 2πV m =0 exp(-π 2 m 2 /α) m 2 S(m)S(-m) + E self (r {N } ), ( 6
)
where α is a positive real number,

S(m) := 1≤j≤N F(L j )(m) exp(2πim • r j ) (7)
is the structure factor and F is the discrete Fourier transform of the operator L j . For example for a point-multipoles up to order p = 2 (quadrupoles), F reads

F(L j )(m) = M 0 j + 2πiM 1 j • m -(2π) 2 M 2 j • M,
with M γγ = m γ m γ . The third term in ( 6) is commonly referred to as the self-energy. In the realm of polarizable force-fields, the commonly used definition of the self-energy is the one from Smith 7 .

The fundamental property of the Ewald energy is that it is independent on the order of summation due of the absolute convergence of the involved sums.

It can be shown 5,[START_REF] Darden | Extensions of the Ewald method for Coulomb interactions in crystals[END_REF] that the interaction energy (2) of the system is related to the Ewald energy through the relation

E(r {N } ) = E Ewald + J(D, M), (8)
where the surface term J(D, M) depends on the dipolar moment D = 1≤i≤N M 0 i r i and the sum of dipoles M = 1≤i≤N M 1 i of the primitive cell U . Only this term is responsible for the order of summation in equation ( 2), it reflects the macroscopic shape of the system (see the upcoming Remark 1 for a discussion on the notion of macroscopic shape). The order of summation of the conditionally convergent series is therefore a factor to choose in order to specify the exact value of the interaction energy E(r {N } ) and is often supposed to be spherical (by shells of n such that |n| is increasing).

By supposing that the macroscopic system is surrounded by a continuum dielectric with some dielectric permittivity ε, the interaction of the microscopic system with the continuum can be taken into account and explicitly dealt with for spherical summation orders. Further, in the limiting case of a perfect conductor ε = ∞ as surrounding environment (and still with spherical summation order), it can be proven that the surface term vanishes [START_REF] Darden | Extensions of the Ewald method for Coulomb interactions in crystals[END_REF][START_REF] De Leeuw | Proc. R. Soc. Lond. A[END_REF] . This model is called the tinfoil model. In consequence, this implies that the energy of the system is in this case the Ewald energy.

In this paper, we do not longer comment on the convergence issues, which will be subject of a forthcoming paper, and concentrate on the proper definition of the self-energy E self (r {N } ), which requires some subtle development if general multipoles are considered that go beyond the results for point charges.

More precisely, there are two aspects that we address in the work. First, we investigate a mathematically clean derivation of the self-potential (and thus of the energy thanks to equation ( 4)) and self-field when general multipoles are considered and not only pointcharges. We then deduce thereof the expression of the self-energy. Second, we present the proofs which demonstrates that these quantities are mathematically well-defined.

II. DERIVATION OF THE POTENTIAL

First, we revisit the derivation of the Ewald summation for the potential generated by the multipoles. The conditionally convergent series in (4) defining the potential φ j is given a precise meaning by considering the limit

lim k→∞ n∈Ω(P,k) 1≤i≤N L i 1 |r i -r j + n| ,
for some domain P in R 3 containing the origin that represents the macroscopic shape of the system (see Remark 1) and where

Ω(P, k) := n = 1≤γ≤3 n γ a γ (n γ ) 1≤γ≤3 ∈ Z 3 , n k ∈ P . (9)
At the base of the derivation of the potential is the splitting

1 |r| = erfc(α 1/2 |r|) |r| + 1 π m U * exp -π 2 |v + m| 2 /α |v + m| 2 exp(-2iπ(v + m) • r) d 3 v, (10)
for any positive α and which can be deduced [START_REF] Darden | Extensions of the Ewald method for Coulomb interactions in crystals[END_REF] from the integral expression of the gamma function at the point 1 2 for all r but at the origin. Using the present splitting and following the arguments presented in Darden 18 (Sections 3.5.2.3.2 and 3.5.2.3.1), one can define

ζ k (r) = n∈Ω(P,k) erfc(α 1/2 |r + n|) |r + n| + 1 πV m =0 exp(-π 2 |m| 2 /α) |m| 2 exp(-2iπm • r) - π αV + H k (r), (11)
such that n∈Ω(P,k)

1 |r + n| = ζ k (r) + o(1),
as k → ∞ and which consists of potential at r that is generated by unit point charges located at the vertices of the lattice indexed by n such that n ∈ Ω(P, k). For sake of completeness, we outline this step in Appendix A where we also give the definition of H k (r) in Eqn. (A2).

Based on ζ k , we now introduce the function

Φ k (r) := 1≤j≤N L j ζ k (r -r j ), (12)
defined everywhere but at the location of the point multipoles located in U . The function r → Φ k (r) represents the potential at r ∈ U generated by all the multipoles and their images contained in periodic lattice cells indexed by n such that n ∈ Ω(P, k).

In consequence, the limit as r tends to any point multipole r j is not finite. Note that this has been handled above with the sign after the sum since only the interaction energy is considered. Instead, if one considers the potential at position r generated by all multipoles except the multipole located at r i , then one has to substract the contribution for n = 0 in equation ( 11) for ζ k (rr i ) to get

φ i k (r) = Φ k (r) -L i 1 |r -r i | , (13)
with finite limit at r i given by

φ i k (r i ) = lim r→r i φ i k (r) = lim r→r i Φ k (r) -L i 1 |r -r i | . ( 14
)
The function r → φ i k (r) denotes the potential at an arbitrary position r generated by all multipoles contained in periodic lattice cells indexed by n such that n ∈ Ω(P, k) except multipole i in the unit cell (n = 0).

Hence, using the splitting introduced in equation ( 11) combined with ( 12) and ( 13), it follows that

φ i k (r) = n∈Ω(P,k) 1≤j≤N L j erfc(α 1/2 |r jn |) |r jn | + 1 πV m =0 exp(-π 2 |m| 2 /α) |m| 2 exp(-2πim • r)S(m) -L i erf(α 1/2 |r -r i |) |r -r i | + 1≤j≤N L j H k (r -r j ), (15)
with r jn := |rr j + n| and where S(m) was defined in equation (7).

We have therefore introduced the splitting (precisely defined in the listing below)

φ i k (r) = φ i 0,k (r) + φ i self (r) + φ i k,surf (r), (16)
where each individual term is defined and discussed in the following.

The absolutely converging part of the potential: Due to their quick convergence in k → ∞ the first two terms in equation (15), denoted by φ i 0,k in equation ( 16), do not depend on the order of summation. The first term is called the direct potential and the second the reciprocal potential.

The self-potential: The third term in equation ( 15) is what we call the self potential 16) and does not depend on the other nuclear position r i , i = j, and k, and is non-constant in r around r i . This term being independent on the other nuclear positions can by no means model the interaction potential, hence the name.

φ i self (r) in (
From the derivation it becomes clear that in the limit r → r i , φ i self (r i ) is the quantity to be subtracted from the reciprocal potential in order that the potential at r = r i is the potential created by all other multipoles except multipole i. Note that the contribution in the direct space has already been taken into account in equation (15) since the sign appears on the first sum.

We will provide in Section V explicit values of this terms in limit r → r i for arbitrary multipolar distributions.

The surface-potential: The fourth term (15), denoted by φ i k,surf in equation ( 16), is the surface potential which will be well-defined only if the sum converges as k tends to infinity. It is intimately linked with the order of summation and is related to subtle questions. We want to focus in this article to the self-terms and are therefore assuming convergence in k here.

Remark 1. It is not very intuitive to understand what is meant by the macroscopic shape of the system and its environment, and how this is mathematically accounted for. From the microscopic viewpoint, the sequence of shapes (Ω(P, k)) k∈N should be seen as the scaling of one macroscopic shape P , i.e. Ω(P, k) = kP . Then, the sequence Ω(P, k) covers larger and larger parts of the microscopic space R 3 as k increases. We would like to advocate also the viewpoint of introducing a change of variables from the microscopic variable n to the macroscopic variable n = n/k that can be used to rewrite sums of the form

n∈Ω(P,k)∩Z 3 f (n) as n∈P ∩(Z 3 /k) f (n), with f (n) := f (k n). ( 17
)
This means that the microscopic space contracts more and more within the macroscopic shape P , see Figure 2 for an illustration. The role of the macroscopic shape P becomes visible and the exterior of P is then the surrounding environment to P .

h fixed kP, k → ∞ (a) Microscopic space P fixed h = 1 k → 0 k → ∞ macroscopic environment (b) Macroscopic space
FIG. 2: The limiting process k → ∞ observed in the microscopic and macroscopic space In the following, we introduce

φ i 0 (r) = lim k→∞ φ i 0,k (r), φ i surf (r) = lim k→∞ φ i k,surf (r), (18)
where we have assumed that the second term converges as we want to study the self-terms.

Thus

φ i (r) = φ i 0 (r) + φ i self (r) + φ i surf (r), (19)
denotes the potential at position r generated by all multipoles except multipole i in the unit cell. Recall that the term φ i surf (r) depends on the order of summation represented by a particular shape P , whereas the other terms φ i 0,k (r) and φ i self (r) do not.

III. DERIVATION OF THE FIELD

The derivation we have given for the potential gives a straightforward one for the field.

It is based on the splitting of φ i k (r) developed in the previous section and uses the fact that the electric field is minus the gradient of the electric potential. Indeed, taking the derivative D r (thus with respect to r) in equation ( 15) yields

E i k (r) = -D r φ i k (r) = - n∈Ω(P,k) 1≤j≤N D r L j erfc(α 1/2 |r jn |) |r jn | - 1 πV m =0 exp(-π 2 m 2 /α) m 2 S(m)D r exp(-2πim • r) + D r L i erf(α 1/2 |r -r i |) |r -r i | - 1≤j≤N D r L j H k (r -r j ). (20)
Therefore, E i k (r) denotes the electric field at a general position r generated by all multipoles in a cell belonging to Ω(P, k) except multipole i in the unit cell.

In consequence, we define each term individually as for the potential:

E i 0,k (r) = - n∈Ω(P,k) 1≤j≤N D r L j erfc(α 1/2 |r jn |) |r jn | - 1 πV m =0 exp(-π 2 m 2 /α) m 2 S(m)D r exp(-2πim • r), (21a)
E i self (r) = D r L i erf(α 1/2 |r -r i |) |r -r i | , (21b)
E i surf,k (r) = - 1≤j≤N D r L j H k (r -r j ). (21c)
In particular, we defined the self electric field E i self (r) as the third term in (20), which will be shown in Section V to be well-defined, in particular at r i , and give explicit expressions.

Evaluating E i k (r) at r = r i then yields E i k (r i ) = E i 0,k (r i ) + E i self (r i ) + E i k,surf (r i ). (22)
Using classical results from convergence of series, we obtain that as soon as the surfacepotential converges (in the limit k → ∞) and the surface-field converges uniformly in r in a neighborhood of r i , the gradient of the limit of the surface-potential is exactly the surface-field:

E i surf (r) = -lim k→∞ 1≤j≤N D r L j H k (r -r j ) = -D r lim k→∞ 1≤j≤N L j H k (r -r j ) = -D r φ i surf (r) . (23)
Again, this is a subtle question related to the convergence in k that will be addressed in an upcoming work. The focus of this article is shed on the self-terms.

IV. DERIVATION OF THE ENERGY

Recalling equation ( 5) combined with the splitting (19) of the potential into different parts, we define the following energy contributions

E 0 (r {N } ) = 1 2 1≤i≤N 0≤k≤p M k i • D k r φ i 0 (r) r=r i , (24a)
E self (r {N } ) = 1 2 1≤i≤N 0≤k≤p M k i • D k r φ i self (r) r=r i , (24b)
E surf (r {N } ) = 1 2 1≤i≤N 0≤k≤p M k i • D k r φ i surf (r) r=r i . (24c)
Note that the self-potential is non-constant in a neighborhood of r i so that the higher multipolar moments, i.e. the derivatives, act on the self-potential φ i self . Further, notice that E 0 (r {N } ) can be written as

E 0 (r {N } ) = 1 2 1≤i≤N L i φ i 0 (r i ) = 1 2 n 1≤i,j≤N L i L j erfc(α 1/2 |r ijn |) |r ijn | + 1 2πV m =0 exp(-π 2 m 2 /α) m 2 S(m)S(-m), (25)
and in consequence, we write

E(r {N } ) = E 0 (r {N } ) + E self (r {N } ) E Ewald (r {N } ) + E surf (r {N } ) J(D,M) = E Ewald (r {N } ) + J(D, M). (26)
Note that we confirm with this derivation equation (8) and that the Ewald energy and the self-energy do not depend on the order of summation whereas the surface energy does.

The corresponding force-terms then naturally result from differentiating the different energies with respect to the nuclear coordinates. In particular, as we will see further below, the self-energy is independent on any nuclear coordinate and the self-energy therefore doesn't induce any force term. However, the correct term of the self-field is mandatory in the context of polarizable force-fields. 16 V. WELL-POSEDNESS OF THE SELF-TERMS

In this part, we outline the proofs that the self-potential, field and (24b) are well-defined in the limit r → r j and in consequence also the self-energy. As done by Smith 7 , we introduce recursively the functions B n for any n ∈ N and all r ∈ R + \ {0} by

B 0 (r) := - erf(α 1/2 r) r B n (r) := 1 r 2 (2n -1)B n-1 (r) + (2α) n √ απ exp -αr 2 . ( 27
)
Then, the following result holds.

Theorem 1. For any n ∈ N, there holds that

lim r→0 B n (r) = B n (0) = - α n+1/2 √ π 2 n+1 2n + 1 , (28)
and dB n dr (r) = -rB n+1 (r).

The proof of Theorem 1 is presented in Appendix B.

In order to give explicit formulae for the self-terms, we first note that from (29) follows

D 1 i B n (|r -r i |) = (r -r i )B n+1 (|r -r i |). (30)
Since we have derived the values of B n (0) in (28), we can give explicit formulae for the self-potential, the self-field and the self-energy in consequence.

For sake of a simple presentation, we consider a multipolar charge distribution up to quadrupoles in the following. Intrinsically, the quadrupolar moments M 2 i are symmetric matrices with zero trace.

The self-potential: Therefore, the i-th self-potential φ self at an arbitrary point r in a neighborhood of r i writes as

φ i self (r) = -L i erf(α 1/2 |r -r i |) |r -r i | = L i B 0 (|r -r i |) (31a) = M 0 i B 0 (|r -r i |) + M 1 i • D 1 i B 0 (|r -r i |) + M 2 i • D 2 i B 0 (|r -r i |) (31b) = M 0 i B 0 (|r -r i |) + B 1 (|r -r i |) M 1 i • (r -r i ) + B 2 (|r -r i |) M 2 i • (r -r i ) (r -r i ) . (31c)
Then, the evaluation of the i-th self-potential at r = r i is given by

φ i self (r i ) = lim r→r i φ i self (r) = M 0 i B 0 (0) = -M 0 i 2 α π , (32)
which does no longer depend on r i , only depends on the charge M 0 i and is well-defined. This formula is of course valid for any kind of multipolar distribution and not restricted to orders to up to quadrupoles only. Note that the self-potential is not constant in a neighborhood of r i in this derivation.

The self-field: We want to stress that in contrast to what is presented in Nymand and Linse 8 , there is indeed a non-zero self-contribution to the electric field. The i-th part of the self-field at r = r i is defined by

E i self (r) = -lim r→r i D r φ i self (r) = -M 1 i B 1 (0) = M 1 i α π 4α 3 , (33)
which only depends on the dipole moment at site i and is also valid for any kind of multipolar distribution.

The self-energy: Finally, the self energy as defined above writes as

E self (r {N } ) = 1 2 1≤i≤N 0≤k≤2 M k j • D k r φ i self (r) r=r i , (34)
= 1 2 1≤i≤N M 0 i • M 0 i B 0 (0) + M 1 i • M 1 i B 1 (0) + 2 M 2 i • M 2 i B 2 (0) (35) = - α π 1≤i≤N M 0 i • M 0 i + 2α 3 M 1 i • M 1 i + 8α 2 5 M 2 i • M 2 i . ( 36
)
We recognize the current practice that for relative energies and forces, the correct term of the self-energy is not needed since a constant misfit cancels out in energy differences.

However, for sake of having a complete theory based on a rigorous development, we think that it is important to state the self-energy as well.

VI. CONCLUSION

In that paper, we proposed a new mathematically clean and coherent derivation of the Ewald summation for a system consisting of N -body electrostatic interaction with multipolar charges of any order. The existing results in the literature differ between different authors and no common development of all quantities can be found. The essential differences lie in the self-term expressions. We presented a clean derivation and confirm the expressions proposed by Smith 7 for which we proved well-posedness. Our model is derived from a clean application of the Ewald splitting to the electric potential and the subsequent quantities such as the electric field, the energy and the forces are derived thereof. A complete derivation of all these quantities is mandatory in the context of next generation polararizable force-fields where in particular the self-field is required and needs to be consistent with the theory.

Overall, the new model which is mathematically sound maintains the use of the tinfoil model and provides simpler expressions for the self-energy that are closer to the original idea of Ewald to work on the potential and not on the energy.

APPENDIX A

Essentially, for sake of a complete presentation we present here the derivation of (11) in a compact way following the arguments presented in Darden 18 , see also [START_REF] Smith | [END_REF]21 . As briefly mentioned in Section II, we start with the following splitting

1 |r| = erfc(α 1/2 |r|) |r| + 1 π m U * exp -π 2 |v + m| 2 /α |v + m| 2 exp(-2πi(v + m) • r) d 3 v, (
1 |r + n| = n∈Ω(P,k) erfc(α 1/2 |r + n|) |r + n| + n∈Ω(P,k) m U * h m,r (v) exp(-2πiv • n) d 3 v with h m,r (v) = exp -π 2 |v + m| 2 /α π|v + m| 2 exp(-2πi(v + m) • r),
and where we used that

exp(-2πi(v + m) • n) = exp(-2πiv • n) since m • n ∈ N.
Case m = 0: We first recognize that

1 V h m,r (n) = U * h m,r (v) exp(-2πiv • n) d 3 v
where h m,r denotes the Fourier coefficient of h m,r and we recall that

V = 1 |U * | . Then, there holds that n h m,r (n) = h m,r (0) = exp -π 2 |m| 2 /α π|m| 2 exp(-2πim • r)
and thus

n∈Ω(P,k) h m,r (n) = h m,r (
- π αV + 1 V H k (r) + o(k),
as k → ∞. Note that we do not shed emphasis on the different arguments that guarantee existence of the different limits but put rather emphasis on the compact development to derive (11).

APPENDIX B: MATHEMATICAL PROOFS

Before we really tackle the proof of Theorem 1, we first prove some auxiliary results.

ure 1 ,

 1 FIG. 1: Computational experiment performed with the Tinker-HP 17 software and the AMOEBA force field. Removing the self-field terms in the computation of the polarization energy gives rise to strong differences in the oxygen-oxygen radial distribution function compared to the correct Tinker-HP initial implementation. Simulation settings: 1ns NVT simulation at 300K, 4000 water molecules within a 49.323 3 Angstrom square box.

π|v| 2 1 - π 2 α |v| 2 += 1 -

 121 0) + o(1), as k → ∞.Case m = 0: As visible from above, this development does not hold for m = 0 and is more subtle. We have thath m,r (v) = exp -π 2 |v| 2 /α π|v| 2 exp(-2πiv • r)and the combination of two Taylor expansions yieldsh m,r (v) = 1 O(|v| 4 ) 1 -2πiv • r -2π 2 |v • r| 2 + O(|v| 3 ) 2πiv • r -2π 2 |v • r| 2 π|v| 2 -π α + O(|v|).This motivates the definitionH k (r) = n∈Ω(P,k) U * 1 -2πiv • r -2π 2 |v • r| 2 |v| 2 exp(-2πiv • n) d 3 v, k) U * exp(-2πiv • n) d 3 v = -π αV + o(1),as k → ∞. Further there holdsn∈Ω(P,k) U * O(|v|) exp(-2πiv • n) d 3 v = o(1),so that combining all terms yieldsζ k (r) = n∈Ω(P,k) erfc(α 1/2 |r + n|) |r + n| + 1 πV m =0 exp -π 2 |m| 2 /α |m| 2 exp(-2πim • r)

Lemma 1 . 2 =Lemma 2 .r 2 (

 1222 For any integer n, the function B n is explicitly given byB n (r) = exp(-αr 2 ) nk) -1)!! -(2n -1)!! erf(α 1/2 r) r 2n+1 ,(B1)where (2n -1)!! := (2n -1) × • • • × 3 × 1 with the convention that for any non-positive integer k, k!! = 1 and that a sum from 0 to -1 is zero.Proof. The proof follows by induction. For n = 0 we see that the proposition holds by inspection. Now, let us assume that (B1) holds for a given n. Inserting (B1) into the definition of B n+1 in (27) impliesr 2 B n+1 (r) = (2n + 1)B n (r) + (2α) n+1 √ απ exp -αrThe functions B n can be rewritten for all positive r asB n (rα +n-k 2 n-k (2n -1)!! !(2(nk) -1)!! r 2( -k-1) -n) . (B2)Proof. The result is obtained by inserting the expression of the error and exponential functions as a power series, i.e., ) r 2 +1(2 + 1) ! , in equation (B1).

  A1) see Eq. (3.5.2.16) in Darden[START_REF] Darden | Extensions of the Ewald method for Coulomb interactions in crystals[END_REF] . Then, one can write

n∈Ω(P,k)

ACKNOWLEDGEMENTS

This work was made possible thanks to the French state funds managed by the CalSim-Lab LABEX and the ANR within the Investissements d'Avenir program (reference ANR-11-IDEX-0004-02) and through support of the Direction Générale de l'Armement (DGA) Maîtrise NRBC of the French Ministry of Defense. Benjamin Stamm acknowledges the funding from the German Academic Exchange Service (DAAD) from funds of the "Bundesministeriums für Bildung und Forschung" (BMBF) for the project Aa-Par-T (project-id 57317909). Yvon Maday and Etienne Polack acknowledge the funding from the PICS-CNRS (Project N • 230509) and the PHC PROCOPE 2017 (Project N • 37855ZK).

We now need to prepare some result that is used in a later proof. Lemma 3. For any n ∈ N 0 , there holds n+1 k=1 (-2) k (nk + 1)!(2k -1)!! + 2 (2n + 1)n! = 0.

Proof. We denote by Γ the usual gamma-function. Introduce as well the double-factorial for even numbers (2n)!! = 2n × (2n -2) × . . . × 2 and observe that the following identities hold

In consequence, there holds

where B(•, •) denotes the beta-function. Next, we observe that

by exploiting the change of variable s = 1t and the fact that t(1t) = s(1s). Further using the identity 4t(1t) = 1 -(2t -1) 2 and another change of variable s = 2t -1 yields

Now, we use (B6) in combination with the the binomial coefficient theorem as follows

Then replacing the double factorial (2k)!! = 2 k k! and shifting the index k by one yields the desired result.

The first claim of Theorem 1 is formulated in the following lemma.

Lemma 4. For any n ∈ N, there holds that

Proof. It remains only to study the limit as r tends to zero in (B2) of the terms wheren andk -1 are nonpositive, as positive powers of r will converge to zero.

In order to have a clear picture, we first reorder the sums over and k of the first term in (B2) by introducing the following change of indices = qs + 1 and k = ns as follows

where S( , k) denotes the summand of the first term in (B2). The second term in (B2) is modified by the change of indices = qs + 1 and k = ns resulting in the following expression

As outlined in the beginning of the proof, we focus on the non-negative powers of r, thus for non-negative qn. The coefficient for such a non-negative power qn of r is given by q+1 s=1 (-2) s (qs + 1)!(2s -1)!! + 2 (2q + 1)q! which vanishes by Lemma 3. This proves well-posedness of the limit r → 0 and the limit is given by the coefficient of the zero-th power in r, i.e. for q = n. Note that min(n, q + 1) = n = q -1 and we apply once again Lemma 3 to obtain the desired limit.

The second claim of Theorem 1 is formulated in the following lemma. Lemma 5. For any n ∈ N 0 , there holds that

Proof. This proof follows by induction. For n = 0, the claim can easily be proven by