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

Abstract

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–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⁹ 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 [scr S, script letter 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 [scr S, script letter S] = 0. The choice of the short-ranged function is however delicate. Many variants have been developed and Table 1 gives a non-exhaustive list of such [scr S, script letter S](q). In this table we have selected three main groups of pair-potentials: reaction field methods; damping 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.

Table 1 Short-ranged functions, [scr S, script letter 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. η = αR_c, where α is the commonly used damping-parameter. For the top three reaction field type methods ε_RF is the relative permittivity of the surrounding medium, and for the bottom scheme P is the number of cancelled moments. Note that [scr S, script letter S](q) = 0 for q > 1, i.e. for r > R_c

Short-range function, [scr S, script letter S](q)	Ref.
Reaction field methods
	Reaction field^10
	7
	11
Damping based methods
erfc(qη)	Real-space Ewald^5
[η = 0 → 1]
erfc(qη) − q erfc(η)	Wolf^1
[η = 0 → 1 − q]
	12
[η = 0 → 1 + q − q^2]
	13
[η = 0 → (1 − q)^2]	14
	3
Damping free methods
(1 − q)^3	15
(1 − q)^4	16
(1 + q)(1 − q)^3	17
(1 + 2q + 2q^2)(1 − q)^4	4
	18
	q-Potential^9

---

In addition to Table 1 there are many other schemes: effective potentials,¹⁹ 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–28 Each scheme has (dis)advantages and the choice of the summation method is sensitive to the specific system. Often short-range pair-potentials 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.²⁹ 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.

Theory

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 isHere the prime indicates that i ≠ 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 ∈ [Doublestruck N]₀, is introduced aswhere r = |r|, which gives ion–ion interactions using l = 0, ion–dipole 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,⁹ 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⁻⁵.

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₀) or image-like dipoles (q₂). The former is found to be more appropriate.

By enforcing total moments described by zero-tensors,^1,3,9i.e. achieving a short-ranged potential, eqn (2) is simplified to a 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 ∈ [Doublestruck N] index the number of cancelled moments, is the q-analogue of the binomial coefficient, and (a;q)_P is the q-Pochhammer symbol.⁹ 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,⁹ 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 higher-order interactions are described by applying the gradient(s) to the same, that is ∇^kT^q₀(r) where k ∈ [Doublestruck N]₀. This approach thus cancels the (implicit) charges of the moments and their higher-order 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 gradient-operator on T^q₀(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 ≤ 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 ∇²T^q₀ to be zero at the cut-off distance (i.e. no higher-order cancellation, or p = 0), then we must use 0 + 2 ≤ P − 1 → P ≥ 3. This is equivalent to cancelling the implicit charge (P ≥ 1), dipole (P ≥ 2) and quadrupole (P ≥ 3) moment.

The second approach simply makes use of the dipole–dipole interaction-tensor T^q₂ directly, and thus 0 + 0 ≤ P − 1 → P ≥ 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⁸ with an upper limit of R_c ≤ min(L)/4. Thus, the long-range interactions determine P and the short-range ones together with the length-scale 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 μ not perpendicular to r, there is an angular alteration to its image-dipole μ′; see eqn (6) where the hats indicate normalized vectors.

μ′ ∝ μ − 3(r̂·μ)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 μ′ ∝ μ. 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₂(r)a(q) + I|r|⁻³b(q) (7)

where I is the identity matrix,for l = 0 and k = 2 (i.e. using image charges), whereas a(q) = (q³;q)_P and b(q) = 0 for l = 2 and k = 0 (i.e. using image-like dipoles).

In the original study⁹ ion–ion interactions (i.e. the q₀-potential) were investigated, and indeed the potential was found to be isotropic enough for water-systems. From a theoretical point-of-view the q₀ approach is more isotropic than that of q₂, due to the isotropy of the image charges and anisotropy of the image-like dipoles. However, this does not exclude the q₂-potential from being isotropic itself. Similarly it does not prove that the q₀ approach is isotropic enough for being an efficient pair-potential using dipole–dipole 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³² in the canonical ensemble were performed using the Faunus software,³³ with N = 3000, density ρ* = 0.924, temperature T* = 1.333, and squared dipole moment μ*² = 3.470, where stars indicate reduced units. The equilibration used N × 10⁵ 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 damping-parameter equal to π/R_c, and a conducting surrounding dielectric medium. The q₀-potential was tested from P = 3 and higher values, and the q₂-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₀ and P = 1 for q₂ are directly comparable in this regard.

Results

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₀ is closer to the reference than q₂, and that the largest discrepancy is between the smallest P-values for both q₀ and q₂. However, while for q₀ going from P = 3 to P = 4 gives more accurate results, for q₂ going from P = 1 to P = 2 gives (slightly) less accurate results. Yet, no q₂-potential rivals the results from any q₀-potential.

Fig. 2 The logarithm of the ratio between pair-potential and reference Ewald (black lines) radial distribution functions using q₀ (top) and q₂ (bottom). The insets show the (overlapping) radial distribution functions.

Dipole–dipole correlations

Fig. 3 shows dipole–dipole correlation differences to Ewald results, and it is clear that q₂ is incapable of capturing the system properties, even being negative at small distances. We also see that q₀(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₀(P = 5) (see ESI‡), and that increasing P beyond these values gives worse agreement.

Fig. 3 Dipole–dipole correlation differences to the reference result using q₀ (top) and q₂ (bottom). The insets show the dipole–dipole correlation 〈[small mu, Greek, circumflex](0)·[small mu, Greek, circumflex](r*)〉. The black lines are Ewald results.

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₀ and q₂ as shown so far and thus further verify the outcome of the analysis in the next section.

Discussion

The q₂-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₀-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₂-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₀-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₀- and q₂-approaches beyond the ones mentioned in the previous paragraph. First, the q₀ image particles can be interpreted as having a non-localized charge-distribution on the shell of a sphere, whereas the q₂ image particles are true point-particles (see Fig. 1). Thus the nature of the q₀-potential is more isotropic than that of the q₂-potential, a feature which seemingly affects the results. Second, the q₀ image particles are exact image particles, whereas the q₂ image particles are image-like as discussed earlier. Since we have found that q₀ provides more accurate results than q₂, which might be a consequence of the just mentioned facts, there may be a corresponding version of the q₂-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 reaction-field method¹⁰ using ε_RF = ∞. By comparing our results to those in a previous study using the described reaction-field method on an identical system as simulated here,²⁸ we note that such an approach is more accurate than q₂ yet less so than q₀.

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⁹). 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₀ using P = 4–5 gives most accurate results for the simulated bulk Stockmayer-system, like for the ionic q-potential for water-systems.⁹ From theory we might however expect P = ∞ 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 short-range 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₀ will render more accurate results than q₂. 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 damping-parameter.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank Swedish Research Foundation (grant number 2017-04372) for financial support and LUNARC in Lund for computational resources.

Notes and references

1. D. Wolf, P. Keblinski, S. R. Phillpot and J. Eggebrecht, J. Chem. Phys., 1999, 110, 8254–8282.
2. X. Wu and B. R. Brooks, J. Chem. Phys., 2005, 122, 044107.
3. I. Fukuda, Y. Yonezawa and H. Nakamura, J. Chem. Phys., 2011, 134, 164107.
4. B. Stenqvist, New J. Phys., 2019, 21, 063008.
5. P. P. Ewald, Ann. Phys., 1921, 369, 253–287.
6. A. Ladd, Mol. Phys., 1977, 33, 1039–1050.
7. I. G. Tironi, R. Sperb, P. E. Smith and W. F. van Gunsteren, J. Chem. Phys., 1995, 102, 5451–5459.
8. B. Stenqvist and M. Lund, EPL, 2018, 123, 10003.
9. B. Stenqvist, V. Aspelin and M. Lund, Generalized Moment Cancellation for Long-Range Electrostatics, 2019, arXiv:1904.10335.
10. J. Barker and R. Watts, Mol. Phys., 1973, 26, 789–792.
11. B. Stenqvist, Electric interactions: a study of cellulose, Department of Theoretical Chemistry, Lund University, Lund, 2016.
12. D. Zahn, B. Schilling and S. M. Kast, J. Phys. Chem. B, 2002, 106, 10725–10732.
13. C. J. Fennell and J. D. Gezelter, J. Chem. Phys., 2006, 124, 234104.
14. M. Levitt, M. Hirshberg, R. Sharon and V. Daggett, Comput. Phys. Commun., 1995, 91, 215–231.
15. S. Kale and J. Herzfeld, J. Chem. Theory Comput., 2011, 7, 3620–3624.
16. B. W. McCann and O. Acevedo, J. Chem. Theory Comput., 2013, 9, 944–950.
17. T. E. Markland and D. E. Manolopoulos, Chem. Phys. Lett., 2008, 464, 256–261.
18. G. S. Fanourgakis, T. E. Markland and D. E. Manolopoulos, J. Chem. Phys., 2009, 131, 094102.
19. S. Izvekov and G. A. Voth, J. Phys. Chem. B, 2005, 109, 2469–2473.
20. H. L. Friedman, Mol. Phys., 1975, 29, 1533–1543.
21. J. Lekner, Phys. A, 1991, 176, 485–498.
22. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
23. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593.
24. A. W. Appel, SIAM J. Sci. Stat. Comput., 1985, 6, 85–103.
25. L. Greengard and V. Rokhlin, J. Comput. Phys., 1987, 73, 325–348.
26. H. Kornfeld, Z. Phys. A: Hadrons Nucl., 1924, 22, 27–43.
27. M. Lamichhane, K. E. Newman and J. D. Gezelter, J. Chem. Phys., 2014, 141, 134110.
28. B. Stenqvist, M. Trulsson, A. I. Abrikosov and M. Lund, J. Chem. Phys., 2015, 143, 014109.
29. F. N. Mendoza, J. López-Lemus, G. A. Chapela and J. Alejandre, J. Chem. Phys., 2008, 129, 024706.
30. D. Wolf, Phys. Rev. Lett., 1992, 68, 3315–3318.
31. A. N. Tichonov and A. A. Samarskij, Equations of mathematical physics, Pergamon, New York, 1963.
32. W. H. Stockmayer, J. Chem. Phys., 1941, 9, 398–402.
33. B. Stenqvist, A. Thuresson, A. Kurut, R. Vácha and M. Lund, Mol. Simul., 2013, 39, 1233–1239.

---

Footnotes

† A C++ implementation of the derived potential as well as other truncation based schemes are available at https://doi.org/10.5281/zenodo.3522058.

‡ Electronic supplementary information (ESI) available: Detailing cutoff effects; derivation of self-energies and the dielectric constant. See DOI: 10.1039/c9cp03875b

---