On short-ranged pair-potentials for long-range electrostatics

In computer simulations, long-range electrostatic interactions are surprisingly well approximated by truncated, short-ranged pair potentials. Examples are reaction field methods; the Wolf method; and a number of schemes based on cancellation of electric multipole moments inside a cut-off region. These methods are based on the assumption that the polarization of the neglected surroundings can be inferred from a local charge distribution. Multipole moments themselves are only approximations to the true charge distribution, approximations which many times are needed to simplify calculations in complex systems. In this work we investigate a new, generalized pair-potential based on the idea of moment cancellation that covers interactions between electrostatic moments of any type. We find that moment cancellation in itself is insufficient to generate accurate results and a more restricted formalism is needed, in our case to cancel the virtual charges of the imposed moments. Thus, it is unfeasible to cancel higher-order moments with explicit higher-order moments such as dipoles and instead image charges are needed. The proposed pair-potential is general and straight forwardly implementable for any electrostatic moment – monopole, dipole, quadrupole, etc. – with a computational complexity scaling with the number of particles in the system.


Introduction
Accounting for long-ranged electrostatic interactions in computer simulations is an exquisite task [1][2][3][4] and although formally exact theories exist for repetitive structures subjected to periodic boundary conditions (PBC), 5,6 these are computationally expensive and may impose artificial symmetry from periodicity. Finite-ranged or truncated pair-potentials are fast alternatives to such lattice-sum models and can be more relevant (and often only valid) for isotropic systems. 7,8 Research of truncated pair-potentials covers many different approaches and we here highlight the concept of electric multipolar moment cancellation. That is, the effect of long ranged electrostatic interactions is approximated by cancelling one, two, or more electric multipolar moments in the cut-off sphere, 1,3 mimicking polarization of the surroundings. The newly developed q-potential 9 for ion-ion interactions is a generalization of this concept, and allows for cancellation of an arbitrary number of moments while being free from empirical damping parameters. In this work we expand the q-potential model to include any type of higher-order interactions, i.e. ion-dipole, dipole-dipole, etc. The methodology applied on dipole-dipole interactions is validated via fluid phase simulations and two distinct strategies for cancelling moments are critically investigated and compared with results from Ewald summation. Our discussion entails connections to existing methodologies and, hence, the next section is a brief overview of the field.

Overview
Developing and deriving approximate pair-potentials for long range electrostatic interactions is an active and long-studied research topic. Instead of rigorously obeying the Poisson equation for the system of interest, an approximate solution is used. For example, for ion-ion interactions, the Coulomb potential is simply scaled with a short-ranged function S(q), see eqn (1).
Here k c is the Coulomb constant, e is the elementary charge, z i and z j are the valence of particles i and j respectively, r is the inter-ionic distance, and q = r/R c , where R c is a spherical cut-off after which S = 0. The choice of the short-ranged function is however delicate. Many variants have been developed and Table 1  based methods; and damping free methods. The reaction field methods are based on a cavity occupied by explicit particles which in turn is surrounded by a dielectric continuum. The particles induce electrostatic moments in the implicit medium which in turn induce a field which interacts with the particles. The damping based methods are mainly evolved from the Ewald summation method, neglecting reciprocal space. To avoid a discontinuity at the cut-off distance, the potential and/or its derivative is shifted to zero. The damping free methods are somewhat differently derived: some are empirically fitted; some are derived as approximations to the Poisson equation; and some implicitly incorporate the damping parameter into a polynomial form. Like the damping based methods they usually also shift the potential and its higher order derivatives at the cut-off. In addition to Table 1 there are many other schemes: effective potentials, 19 image methods, 2,20 lattice-methods, 5,6,21-23 and the fast-multipole-method 24,25 among others. However, many times the charge distribution of a system is anisotropic and it is convenient to approximate it by electric point multipole moments instead of point charges. These moments each account for a portion of the charge distribution and do therefore inherently contain the (now implicit) charges. Several of the mentioned ion-ion approaches have also been expanded to, or do cover, multipole moment (i.e. higher-order) interactions. 6,10,[26][27][28] Each scheme has (dis)advantages and the choice of the summation method is sensitive to the specific system. Often short-range pairpotentials for electrostatics are derived based on the assumption of a homogeneous and isotropic distribution beyond the cut-off region. Therefore such approaches should be used with care in systems where this might not be true, for example at interfaces. 29 In the following we expand the theory of the q-potential to higher order interactions and connect the formalism to existing schemes presented in Table 1.

Multipolar energy in a periodic system
For a cuboidal unit cell with periodic boundaries, the total interaction energy between N particles with electric point moments of type v and w is Here the prime indicates that i a j when n = 0, r ij is the distance-vector between the centers of moments i and j, is the Hadamard product, and the size of the cuboid cell is described by side-lengths L = (L x ,L y ,L z ). The lth order interaction-tensor T l (r), using l A N 0 , is introduced as where r = |r|, which gives ion-ion interactions using l = 0, iondipole using l = 1, dipole-dipole using l = 2, and so forth. Note that the factor 1/2 in eqn (2) is included only for like type moments v and w.
Derivation of the multipolar q-potential Consider a local region of an isotropic system -a spherical cut-off sphere, for example -where the contained charge distribution is self-consistently polarized by its surroundings. We now describe both the local charge distribution and the surroundings as two multipoles. These will be oppositely polarized and, given that the local region is sufficiently large, perfectly cancel each other. This physical relation is the main idea behind the q-potential, 9 where image charges (or moments) of the local region are used to generate the opposing multipole of the surroundings (see Fig. 1). Perfect cancellation of the local and surrounding multipole moments is achieved by requiring that their sum results in zero-tensors, leading to an effectively short-ranged potential. Therefore, for a reasonably large cut-off region, the presented approach is valid for both PBC and non-PBC systems. A similar observation has been made in lattice systems 6,30 where the effective Coulomb interaction rapidly decays as r À5 . By enforcing total moments described by zero-tensors, 1,3,9 i.e. achieving a short-ranged potential, eqn (2) is simplified to a Table 1 Short-ranged functions, S(q), for various electrostatic schemes for ion-ion interactions where q = r/R c with r being the distance between the charges and R c the cut-off. Z = aR c , where a is the commonly used damping-parameter. For the top three reaction field type methods e RF is the relative permittivity of the surrounding medium, and for the bottom scheme P is the number of cancelled moments. Note that S(q) = 0 for q 4 1, i.e. for r 4 R c
Reaction field methods Damping based methods erfc(qZ) Damping free methods (1 À q) 3 15 (1 À q) 4 16 (1 + q)(1 À q) 3 17 View Article Online single cell n = 0, or for non-PBC systems, a local region. This has previously been done for ion-ion interactions (see especially eqn (9) in ref. 9) and generalizing this procedure, we obtain the modified interaction tensor for any type of electrostatic interactions, Here P A N index the number of cancelled moments, P p h i q is the q-analogue of the binomial coefficient, and (a;q) P is the q-Pochhammer symbol. 9 Since the multiplicative factor to the original interaction-tensor is a q-analogue, we index the modified interaction-tensor with a superscript q and from here on include the derived potential in the q-potential notation. Like for the Coulomb q-potential, we deduce the number of cancelled higher-order moments P À 1 to equal the number of derivatives of T q l (r) with respect to r to be zero at the cut-off. From now on, when we mention the derivative of T q l (r), it is with respect to r.

Moment cancellation schemes
We now detail two moment cancellation schemes -see Fig. 1 differing in how they cancel moments. The first approach, which has previously been validated for ion-ion interactions, 9 uses image charges (l = 0) while the second approach uses image-like dipoles (l = 2) to cancel higher-order moments. We later revisit the meaning of image-like dipoles.
The image charge approach is derived using the ion-ion interaction-tensor from the Coulomb q-potential, where higherorder interactions are described by applying the gradient(s) to the same, that is r k T q 0 (r) where k A N 0 . This approach thus cancels the (implicit) charges of the moments and their higherorder moments by using image charges. Note that there is an arbitrariness in image charge positioning while still achieving moment cancellation, yet the q-potential is based on a physically inspired scheme. 9,20,31 For each use of the gradientoperator on T q 0 (r), one loses the highest order derivative of the interaction-tensor to be zero at the cut-off. This is expressed in eqn (5), for any l, where the gradient-operator is applied k times, leaving p r P À 1 À k higher-order derivatives zero at the cut-off.
Therefore, if for example we require the dipole-dipole interaction-energy (l + k = 2) using r 2 T q 0 to be zero at the cutoff distance (i.e. no higher-order cancellation, or p = 0), then we must use 0 + 2 r P À 1 -P Z 3. This is equivalent to cancelling the implicit charge (P Z 1), dipole (P Z 2) and quadrupole (P Z 3) moment.
The second approach simply makes use of the dipole-dipole interaction-tensor T q 2 directly, and thus 0 + 0 r P À 1 -P Z 1 is enough for the interaction-energy to be zero at the cut-off distance. Here the explicit dipole moment is cancelled but not the implicit charges, nor necessarily the quadrupole moment.
We define q l (P) to be the q-potential, where the cancellation is based either on image charges (l = 0) or on image-like dipoles (l = 2), which cancel P À 1 higher-order moments where P is counted from what type of base moment is used in the cancellation procedure. When not using the index (P), like in q l , we refer to the general group of pair-potentials using l but for any P.
The presented q-potential approaches can be viewed as a generalization of global electroneutrality. Since long-range electrostatic interactions solely come in the form of ion-ion, ion-dipole, dipole-dipole, and ion-quadrupole interactions, we recognize that to accurately compensate for the neglected contributions outside of the cut-off region we generally need to cancel at least up to the quadrupole-moment. Higher order interactions are short-ranged and thus given a large enough cut-off region, these interactions will be explicitly accounted for. However, the cut-off region will also need to be large enough as to capture the local anisotropy of the system 8 with an upper limit of R c r min(L)/4. Thus, the long-range interactions determine P and the short-range ones together with the lengthscale of the local anisotropy determine R c . Now we address the use of image-like dipoles. For a charge at position r relative to the origin, its image charge in a conductive sphere with radius R c centered in the origin has a charge scaled by ÀR c /r as compared to the charge at r. Since potentials due to charges are isotropic, the charge and image charge have the same angular (in)dependence with regard to the origin. For dipoles however this is not true. Given a dipole moment l not perpendicular to r, there is an angular alteration to its image-dipole l 0 ; see eqn (6) where the hats indicate normalized vectors. l 0 p l À 3(rÁl)r (6) In this work we have used the exact image moment positions, yet fixed the angular part of the image moments, in the dipole example that is l 0 p l. Thus these dipoles are image-like. We later discuss how this may affect the results. Summarizing the above, the final interaction tensor for dipole-dipole interactions using the q-potential is given by T q (r) = T 2 (r)a(q) + I|r| À3 b(q) (7) Fig. 1 Illustration of how dipole moment cancellation is used to approximate the exact electrostatic energy using two different schemes. The pair-potential truncation neglects all interactions beyond the cut-off radius which is corrected for using either image charges (q 0 ) or imagelike dipoles (q 2 ). The former is found to be more appropriate.
In the original study 9 ion-ion interactions (i.e. the q 0 -potential) were investigated, and indeed the potential was found to be isotropic enough for water-systems. From a theoretical point-ofview the q 0 approach is more isotropic than that of q 2 , due to the isotropy of the image charges and anisotropy of the image-like dipoles. However, this does not exclude the q 2 -potential from being isotropic itself. Similarly it does not prove that the q 0 approach is isotropic enough for being an efficient pair-potential using dipoledipole interactions. That is why we in the next section investigated the presented potentials using numerical simulations.

Metropolis Monte Carlo simulations
Metropolis Monte Carlo simulations of dipolar Stockmayer fluids 32 in the canonical ensemble were performed using the Faunus software, 33 with N = 3000, density r* = 0.924, temperature T* = 1.333, and squared dipole moment m* 2 = 3.470, where stars indicate reduced units. The equilibration used N Â 10 5 steps each consisting of a combined translational and rotational Monte Carlo attempt. Production runs were ten times longer. The Ewald summation method was used as a reference, and comparisons will always be with regard to this scheme if not stated otherwise. The used cutoff was R c * = 4 for the pair-potentials while Ewald summations used a real-space cut-off of half the box-length, a reciprocal-space spherical integer cut-off of 9, a dampingparameter equal to p/R c , and a conducting surrounding dielectric medium. The q 0 -potential was tested from P = 3 and higher values, and the q 2 -potential from P = 1. These values were chosen since their respective interaction tensors are zero at the cut-off distance for such P, with no higher-order derivatives cancelled. Thus P = 3 for q 0 and P = 1 for q 2 are directly comparable in this regard.

Radial distribution functions
To validate the developed pair-potentials, we first compare radial distribution functions, g(r), with reference Ewald summation results. Fig. 2 shows the logarithm of the ratio between q l -potential and Ewald radial distribution functions. It is clear that q 0 is closer to the reference than q 2 , and that the largest discrepancy is between the smallest P-values for both q 0 and q 2 . However, while for q 0 going from P = 3 to P = 4 gives more accurate results, for q 2 going from P = 1 to P = 2 gives (slightly) less accurate results. Yet, no q 2 -potential rivals the results from any q 0 -potential. Fig. 3 shows dipole-dipole correlation differences to Ewald results, and it is clear that q 2 is incapable of capturing the system properties, even being negative at small distances. We also see that q 0 (P = 4), i.e. up until the first derivative of the interaction-tensor is zero at the cut-off, is closest to the reference akin to q 0 (P = 5) (see ESI ‡), and that increasing P beyond these values gives worse agreement.

Dipole-dipole correlations
In the ESI, ‡ we have included more results from the simulations (energies, mean squared dipole moments), and presented results for different cut-offs. Those data qualitatively show the same behaviour of q 0 and q 2 as shown so far and thus further verify the outcome of the analysis in the next section.

Discussion
The q 2 -potentials cancel moments by using image-like dipoles, yet the simulation results are far from the Ewald summation reference. By expanding the ion-ion interaction-tensor to dipolar systems, i.e. by using the q 0 -potentials which cancel moments by means of image charges, we get more reasonable outcomes. Though both approaches make use of moment cancellation there is a salient difference in how they do so, which is reflected in the simulation results. We conclude that moment cancellation in itself is insufficient whereas neutralizing the (implicit) charges of the moments and cancelling their higher-order moments provide accurate results. This contrasts other works which like the q 2 -potential uses explicit dipole-cancellation, but needs to introduce an arbitrary damping parameter to remedy poor results. 27,28 A single higher-order cancellation is sufficient in many regards to accurately retrieve valid results when using the q 0 -potential. This is convenient from a computational point-of-view since molecular dynamics equivalently requires zero force at the cutoff. A modest increase in the cancellation-order improves the results slightly. There are two main theoretical differences between the q 0 -and q 2 -approaches beyond the ones mentioned in the previous paragraph. First, the q 0 image particles can be interpreted as having a non-localized charge-distribution on the shell of a sphere, whereas the q 2 image particles are true point-particles (see Fig. 1). Thus the nature of the q 0 -potential is more isotropic than that of the q 2 -potential, a feature which seemingly affects the results. Second, the q 0 image particles are exact image particles, whereas the q 2 image particles are image-like as discussed earlier. Since we have found that q 0 provides more accurate results than q 2 , which might be a consequence of the just mentioned facts, there may be a corresponding version of the q 2 -potentials which utilizing explicit dipole cancellation and generates accurate results using non-localized image dipole distributions, or exact image particles. For example, by using a non-localized dipole moment at the cut-off sphere (P = 1) instead of image-like dipoles, we regain the reactionfield method 10 using e RF = N. By comparing our results to those in a previous study using the described reaction-field method on an identical system as simulated here, 28 we note that such an approach is more accurate than q 2 yet less so than q 0 .
Furthermore, we note that the currently proposed formalism is only one of formally infinitely many using the concept of electrostatic cancellation of moments. This since the position of the image moments can be altered while keeping the feature of moment-cancellation intact (see eqn (7) in the derivation of the q-potential 9 ). The positions chosen in this work, and for the Coulomb q-potential, are however physically based and in line with previous works on image moments. 9,20, 31 We have found that q 0 using P = 4-5 gives most accurate results for the simulated bulk Stockmayer-system, like for the ionic q-potential for water-systems. 9 From theory we might however expect P = N to give most accurate results since all moments then are cancelled and no explicit interactions with the surrounding are sustained. Nonetheless, most systems exhibit some form of local fluctuations and a too large value of P will therefore suppress this inherent quality. Higher order interactions than dipole-dipole or ion-quadrupole are shortrange and thus can be captured with a large enough cut-off region. Hence, infinite cancellation is not necessarily needed. Assuming that the fluctuations of the cut-off region will decay while increasing its size, a large P in combination with a large enough R c will still give accurate results. Thus, the order of the fluctuations in the system seems to determine the upper limit of P while the order of long-range interactions determines the lower limit, P = 3.
Finally, the presented formalism based on the T q l (r)-tensors is applicable, though not tested, for other electrostatic interactions. For example, the ion-quadrupole interaction makes use of the same interaction-tensor as in the dipole-dipole interaction. We therefore surmise our presented results to be generalizable to this case, i.e. q 0 will render more accurate results than q 2 . Though there is little need to use summation methods for even higher-order electrostatic interactions, the procedure is valid for similar cases outside the electrostatic framework in which arbitrary higher-order moments need cancellation.

Conclusion
We have expanded a formalism accounting for long-ranged Coulomb (ion-ion) interactions by means of moment cancellation to include higher-order electrostatic interactions. The expansion was done in two fundamentally different ways: cancellation of the implicit charges inherent to any higher-order electrostatic moment and their higher-order moments by means of image charges; and cancellation of the explicit moments and their higher-order moments by means of images of said moments. The results unambiguously show moment cancellation to be insufficient to generate accurate results, rather we find that cancellation of the (implicit) charges inherent to the used moments gives comparable results to the standard Ewald summation.
The expanded validated formalism is general and straight forwardly implementable for higher-order electrostatic interactions with a computational cost proportional to the number of particles in the system. The method utilizes a cut-off region large enough to represent the sample, and a cancelling parameter directly connected to the relevant higher-order moments in this region, and it therefore avoids any arbitrary dampingparameter.

Conflicts of interest
There are no conflicts to declare.