Electrostatic pair-potentials based on the Poisson equation

Electrostatic pair-potentials within molecular simulations are often based on empirical data, cancellation of derivatives or moments, or statistical distributions of image-particles. In this work we start with the fundamental Poisson equation and show that no truncated Coulomb pair-potential, unsurprisingly, can solve the Poisson equation. For any such pair-potential the Poisson equation gives two incompatible constraints, yet we find a single unique expression which, pending two physically connected smoothness parameters, can obey either one of these. This expression has a general form which covers several recently published pair-potentials. For sufficiently large degree of smoothness we find that the solution implies a Gaussian distribution of the charge, a feature which is frequently assumed in pair-potential theory. We end up by recommending a single pair-potential based both on theoretical arguments and empirical evaluations of non-thermal lattice- and thermal water-systems. The same derivations have also been made for the screened Poisson equation, i.e. for Yukawa potentials, with a similar solution.


Introduction
Electrostatics is an important part of molecular interactions and is often included in simulations of such systems. However, the long-ranged nature of the Coulomb potential often makes it impossible to include the exact electrostatic contribution between a particle and its (approximately) infinite surrounding. For a replicated finite simulation-cell, the Ewald summation method [1] and varieties of it [2] are perhaps the most widely used schemes to account for electrostatic interactions. These procedures are nonetheless computationally expensive which is one reason truncated pair-potentials are often used. A truncated pair-potential has the property that it gives zero potential outside of some bounded region, most often a sphere with radius R c . This effectively makes interactions finite-ranged, which in turn reduces the order of the computational complexity. Many such pairpotentials have been developed based on effective interactions [3], electrostatic moment cancellation [4][5][6], and statistical distribution of anisotropic regions [7].
In this study we aim to find a pair-potential which most accurately obey the Poisson equation, which is the most fundamental classical description of electrostatics as of yet. Unsurprisingly we find that no truncated Coulomb pair-potential can solve the Poisson equation exactly, yet we propose approximate solutions which each obeys different aspects of it. After the derivations of the mentioned pair-potentials, we discuss their respective limitations based on both thermal and non-thermal results, and end up by recommending a single pair-potential.

Theory
The fundamental equation in classical electrostatics is the Poisson equation which relates the electric potential ψ at a given location r to the total volume charge-density ρ as e y r -  = · ( where ∇is the gradient operator and ε 0 the vacuum permittivity. For a single point-particle with charge ez, where e is the elementary charge and z the valency, this equation simplifies to where ∇ 2 is the Laplace operator, the delta function d = ¥ ( ) r for = ( ) r 0, 0, 0 , zero otherwise, and . The exact solution to equation (2) is the Coulomb potential . This potential is long-ranged and can thus many times not be used explicitly since the number of interactions is essentially infinite, i.e. unmanageable for any computer. In this work we therefore assume that an approximate solution to the Poisson equation exists and can be described by where ( ) q is a short-ranged function (i.e. zero outside of some bounded region) and q=r/R c . By using this approximate solution to the Poisson equation, equation (2) transforms to where the superscript of ( ) q index the order of derivative with respect to q. By assuming  = ( ) 0 1, which in appendix A is justified for any short-ranged effective potential, equation (5) simplifies to The short-ranged function is now assumed as where Q(q) is twice differentiable (or class  2 ), and θ(q) is the Heaviside step function defined in equation (7).
q thus has the property that it is zero for q>1 (or equivalently r>R c ) which gives it a finite range. For such a short-ranged function where the two far-right terms are, for general Q(q), non-zero at q=1. Yet, by constraining = using the integers dä{0, 1}, designated as the cut-off constraint, these contributions disappear and equation (6) transforms to which we denote as the second derivative constraint. The solutions given by the cut-off constraint are incompatible with the solutions given by the second derivative constraint, cf the divergence theorem. That is, equation (5) has no solution using the assumption of a short-ranged effective potential. The physical implication of this is that the flux through the surface of the bounded cut-off region does not match the charge-density inside of the same. Therefore, one has to choose which physical property (i.e. constraint) to match and which to disregard. There are multiple functions which would suit as approximate solutions to the given problem, yet here we restrict Q(q) to a Kth degree polynomial as in equation (11), where a k are polynomial coefficients. Note that since  = ( ) 0 1 we get a 0 =1.
The cut-off constraint By requiring equation (9) using dä{0, 1} to hold, and thus neglecting the requirement of equation (10), the solution will include the factor - is an integer determining the smoothness of Q(q) at the cut-off. Though not obeying equation (10) strictly, it can be violated in different ways, especially at the singularity of y ( ) r but also at the discontinuity of θ(1−q). By iterative procedure it is straight-forward to acknowledge that the (C−2)th derivative of the left-hand-side of equation ( and equation (9) holds using c, dä{2, 3, K, C}. Note that the former implies a c =0, and that C determines the smoothness of Q(q) at the origin. In order to get the left-hand-side of equation (10) as smooth as possible for qä[0, 1], C and thus implicitly also D must be maximized. By choosing the lowest possible degree of Q(q) given these requirements, the unique solution present itself in equation (12).
We acknowledge that equation (12) is indeed the unique lowest degree solution based on the following: the factor of (1−q) D+1 which requires D+1 restrained a k coefficients, the expanded right-hand-side of equation (12) gives a c =0 for cä{2, 3, K, C} and thus C−1 restrained a k coefficients, and finally that a 0 is restrained to unity, making the total sum of D+C+1 restrained coefficients which corresponds exactly to the number of coefficients in equation (12). By using equation (12) we furthermore note that which by using the variable substitutions =m C 1 and n=D+C−2 transforms to The right-hand-side of this equation has the shape of a binomial distribution and for large n, i.e. large C or D, this binomial distribution can be well-approximated by a normal distribution ( ) q with mean value μ=m/n and variance σ 2 , see equation (15). The factor (n + 1) of the right-hand-side is used to force the integral over the interval q ä [0, 1] to be unity for any m and n.
By integrating ( ) q twice we get an analytic expression for Q(q), see equation (16) where A and B are constants, and erfc is the complementary error-function.
In order to increase the smoothness of equation (10) at the entire interval qä[0, 1], an increase of either C or D is insufficient (since equation (10) is class  2 ), instead they both need to be equal and thus we get μ=1/2. The resulting ( ) Q q is presented in equation (17) where we for sufficiently large C and = ( ) S 0 1get  / B 1 2, and by using equation (9) We have now seen that Q(q) can be described among others by the complementary error-function. The complementary error-function has been utilized in many previous summation schemes without solid theoretical justification [4][5][6][8][9][10][11][12][13], whereas we here connect this (approximate) Gaussian distribution of the charge directly to requirements of the Poisson equation. However, whereas seemingly all such schemes includes erfc(η q) in their short-ranged functions, where η is some arbitrary parameter, we argue that it should rather be the shifted

The second derivative constraint
By instead obeying equation (10), then equation (9) using dä{0, 1} must be violated. Solely neglecting equation (9) using d=0 gives Q(q)=1, and solely neglecting equation (9) using d=1 gives Q(q)=1−q. Note that equation (12) using C=−D gives the former, and using C=1 and D=0 gives the latter. Therefore equation (12) provides all possible lowest order degree polynomial approximate solutions to the Poisson equation, irrespective of which constraint we choose to obey.

Comments
We will now show that in the limit of infinite cancellation of derivatives (i.e.  ¥ C or  ¥ D ), the incompatible cut-off and second derivative constraints gives qualitatively equivalent solutions. In this limit we get σ → 0 and thus  d m - q. This means that equation (12) (i.e. a solution to equation (9)) gives Q (2) (q)∼δ(μ−q)/μ. By using  = ( ) 0 1and equation (9) using either d=0 or d=1, we then get . This expression is for all intent and purposes qualitatively equivalent to which is as solution to equation (10). The solutions using different incompatible constraints to the Poisson equation have now been explored and we found that they can all be described by equation (12). Based on this equation and different smoothness parameters C and D, we therefore move on to present results in two systems: energies of a non-thermal lattice, and densities and radial distribution functions in thermal liquid water. Finally, we want to mention that the same derivation also has been made for the screened Poisson equation, i.e. for Yukawa potentials, see appendix B. There we found that equation (12)

Methods
The non-thermal system was composed of a CsCl lattice, based on a primitive cubic cell, with the Madelung constant 2.04 [14]. For the thermal setup we used the SPC/E water model [15], 2000 molecules and the production runs spanned 10ns. Molecular dynamic simulations was utilized in the isobaric-isothermal ensemble at pressure 1 bar and temperature 298.15 K, using OpenMM 7.0.1 [16]. The pressure was kept constant using a Monte Carlo barostat [17] and we used a Langevin integrator [18], with a friction coefficient of 1.0ps and a time step of 2.0fs, to keep constant temperature. Before production runs, the systems were energy minimized, then equilibrated during 1 ns. The reference Ewald simulations utilized a fractional force error tolerance of 5×10 −4 , and we used R c =0.96 nm for all potentials.

Madelung energies
In figure 1 we present the Madelung energy for a CsCl lattice, and approximated values using the derived pairpotentials (including the self-energy retrieved in appendix A). When observing the D-lines, we note that their oscillations increases in amplitude as C grows, and for C2 they do so around the Madelung energy. Increasing D results in lower oscillations and the curves seems all but shifted horizontally. By observing the extreme curves D=1 and D=5 we note that the optimal lines generally lies somewhere in-between, or more specifically at C;D. The grey lines, based on equation (17), gives the most accurate results using low C2, and for higher C the results are indistinguishable from the corresponding more erroneous equation (12) results. We have excluded results using Q(q)=1 or Q(q)=1−q since they show no sign of converging within the presented range. For the interested reader we refer elsewhere for such results on a NaCl lattice [4]. We have made the same analysis on such a NaCl lattice with qualitatively equal results as for the CsCl lattice.

Simulations
In table 1 we present the densities for the simulated water-systems. We note a somewhat linear relationship between C and D giving accurate results, C;D+{0, 1}, which becomes more lax as D increases. Like in the previous section, the Q(q)=1 and Q(q)=1−q pair-potentials gives fairly erroneous results why we exclude them from further analysis in this section. For results using the latter on the same system as used here we refer elsewhere [19].
In figure 2 we present the logarithm of the ratio between the pair-potential and Ewald radial distribution functions between oxygen atoms. An increase of D generally seems to increase the accuracy at the cut-off whereas decreasing the accuracy at small distances. The reverse hold for C where an increase yield more accurate results at small distances but larger errors at the cut-off. Other ranges are fairly accurate as long as both these regions have low errors. We note that results using ³ D 3 resembles Ewald at the cut-off for some C2, and that C2 is able to reproduce Ewald results at small distances for some D2. For low values of C, C;D+1 gives accurate results on the entire interval whereas for higher values of C the choice of D is, again, more lax. Equation (17) gives accurate results using C=4, for lower C the results primarily diverge from Ewald at the cutoff whereas higher C gives worse results at small distances.

Discussion
Starting with the cut-off constraint, we then moved on to increase the smoothness of Q (2) (q) at q=0 and q=1, which equation (12) does to an arbitrary degree. However, the aim should also be to decrease the absolute Q (2) (q) function-values at the entire interval as to decrease the (maximum) absolute errors. The maximum value of  (17) and Cä [3,4,5]. For lower values of C the radial distribution functions are off by 50% over the entire interval. ( ) q (and thus Q (2) (q) for reasonable large C or D) is inversely proportional to σ. Hence, by maximizing σ, one may decrease the maximal absolute error of equation (10), yet by doing so fewer higher order derivatives will be (approximately) zero at q=0 and q=1. By changing σ, one may thus increase the smoothness at the expense of the maximal error, or vice verse. Similar arguments holds for low C and D yet not in terms of the normal distribution which we used for simplicity to make our case. Since μ only shift this maximum error peak, we find that no smoothness parameter can improve every feature, i.e. if you win some, you lose some.
In the previous section we found that C;D gives the most accurate results in the investigated non-thermal systems and C;D+{0, 1} gives the most accurate results in the tested thermal systems. Thus optimal parameters for ordered systems seems to imply optimal parameters in unordered systems, yet the reverse does not necessarily hold. For reasonably large values of C (or D) C;D implies that equation (12) can be wellapproximated by equation (16) which is a relatively computationally inexpensive function and shows great resemblance with other short-ranged functions suited for summation of long-ranged electrostatics in thermal systems [4][5][6][8][9][10][11][12][13]. In fact, some of these are special cases of equation (12) with specific values of C or D. For example [9][10][11] using C=D=1, C=D=2, and C=D+1=4, which thus support our conclusions on the optimal parameters. Since C and D are directly linked to the smoothness at q=0 and q=1 respectively, and CD seems to hold for all systems, we note that the impact of the singularity of y ( ) r in equation (10) is generally greater (or equal) to the impact of the irregularity of the Heaviside step function at the cut-off.
We now address the impact of the cut-off value R c . For any regular finite Q (2) (q) we note that the left-handside of equation (10) can be made arbitrary small given an increase of R c . Thus for large R c the parameterization can almost be neglected whereas for small R c it becomes increasingly important. The feature of diminishing the maximum value of the left-hand-side of equation (10) as R c increases, can also be found by increasing σ, using the normal distribution analog, however as we discussed earlier the width of the distribution will then also widen, leading to a lower degree of smoothness. Thus no parameter can match R c in its ability to exclusively decrease the maximal absolute error of equation (10) with its general shape fixed.
For isotropic setups using periodic boundary conditions, R c is argued to be no larger than a fourth of the shortest simulation-cell side-length in order to avoid inter-cell correlations [20]. Yet it still needs to be larger than any local anisotopic region which thus restrict the lower bound of the simulation-cell size. In this study we found that for highly ordered systems like the CsCl and NaCl lattices with lattice constant a, no more than R c =2a is needed to get accurate energies. This makes sense as the structure has then repeated once, which effectively provides all needed information about the infinite lattice setup, if the pattern is periodic that is, which indeed is the case for a lattice.

Conclusion
We have presented truncated pair-potentials for electrostatic interactions which as closely as possible solves the Poisson equation. Starting from the Poisson equation, we find two incompatible constraints on Q(q) (∼the short-ranged function): cancellation of the zeroth and first order derivative at the cut-off, and zero-valued second order derivative at the entire interval. The results strongly suggests the cut-off constraint to be more important than the second derivative constraint. By using the cut-off constraint and increasing the smoothness of the second order derivative (by means of derivative cancellation), we find a general expression for a pairpotential which surprisingly can meet either of the constraints pending two physically connected smoothness parameters. In the limit of infinite cancellation we, by using the cut-off constraint, qualitatively regain one of the pair-potentials derived by using the second derivative constrain, which thus connects the two incompatible constraints. Since the constraints are equivalent in the limit of infinite cancellation, and the second derivative constraint generally yield pair-potentials with inaccurate results, we conclude that an optimal choice in the procedure of constructing an accurate physical truncated pair-potential must use the cut-off constraint and finite cancellation. Furthermore, we find that for sufficiently large smoothness, the short-ranged function can be well-approximated by a Gaussian distribution of the charge, a feature which has previously been widely used without theoretical justification. We find that such a distribution does indeed produce accurate thermal results for various low variances (high C=D), yet non-thermal results are accurate for higher variance (low C=D2). We end up by noting that many other pair-potentials are special cases of the derived equation (12), however we recommend it using C=D=3, i.e. Equation (20), as the optimal short-ranged function for both non-thermal and thermal systems.
The first term in the sum is a short-ranged interaction and the second term is a longranged interaction. In this work we have assumed a short-ranged effective pair-potential, which gives the energy Note that the energy is neglected by solely using  E compared to E. However,  E is not necessarily zero which is why we now will investigate its value.
The assumption of an effective short-ranged potential is exact when the sum of all long-ranged interactions cancel, i.e.
i j i j ij ij Self is denoted as the self-energy. Ideally, the long-ranged function ( ) q should be constructed such that each separate such contribution is zero for  | | r 0. This includes every pair-interaction, including charges at distance = | | r 0 ij , which in practice is only the case for i=j due to the Pauli principle. The terms in equation (A6) diverges for any  ¹ ( ) 0 0(or equivalently  ¹ ( ) 0 1) which is why we constraint ( ) 0 to be zero from here on. By assuming  q = -( ) ( ) ( ) q Q q q 1 and Q(q) as in equation (11), we thus get a 0 =1 and where by using equation (12) we note that = -+ ( ) a C D C 1 . With E Self  known, the total interaction-energy is then Here the superscript of (˜) q index the order of derivative with respect toq. This equation has the same form as equation (6) (save the quadratic term which is non-zero for κ>0) and thus the same approximate solutions applies, i.e. equation (12) where q is replaced byq.