Non-Ewald methods for evaluating the electrostatic interactions of charge systems: similarity and difference

In molecular simulations, it is essential to properly calculate the electrostatic interactions of particles in the physical system of interest. Here we consider a method called the non-Ewald method, which does not rely on the standard Ewald method with periodic boundary conditions, but instead relies on the cutoff-based techniques. We focus on the physicochemical and mathematical conceptual aspects of the method in order to gain a deeper understanding of the simulation methodology. In particular, we take into account the reaction field (RF) method, the isotropic periodic sum (IPS) method, and the zero-multipole summation method (ZMM). These cutoff-based methods are based on different physical ideas and are completely distinguishable in their underlying concepts. The RF and IPS methods are “additive” methods that incorporate information outside the cutoff region, via dielectric medium and isotropic boundary condition, respectively. In contrast, the ZMM is a “subtraction” method that tries to remove the artificial effects, generated near the boundary, from the cutoff sphere. Nonetheless, we find physical and/or mathematical similarities between these methods. In particular, the modified RF method can be derived by the principle of neutralization utilized in the ZMM, and we also found a direct relationship between IPS and ZMM.


Introduction
Molecular dynamics (MD) and Monte-Carlo (MC) calculations enable investigation of macroscopic systems in terms of microscopic viewpoint with fine resolutions of time and space (Allen and Tildesley 1987;Hoover 1991;Schlick 2010). They are now widely used for studying the nature of biomolecular systems, including proteins, lipids, nucleotides, and surrounding solvent and ions, to reveal their structures, chemical properties, and biological functions. Electrostatic interactions among the particles of the system much affect these properties and so are critical in these molecular simulations (Nakamura 1996;Karttunen et al. 2008;Ren et al. 2012). The electrostatic interactions  Graduate School of Information Science, University of Hyogo, 7-1-28 Minatojima, Minamimachi, Chuo-ku, Kobe, Hyogo, 650-0047, Japan 2 Institute for Protein Research, Osaka University, 3-2 Yamadaoka, Suita, Osaka 565-0871, Japan are originally defined by the Coulombic function that is long ranged, so that its sufficient approximation is nontrivial. In addition, calculation of the electrostatic interactions is the most time-consuming part in the simulation. Thus, considerable attention should be paid for calculating the electrostatic interactions.
For calculating electrostatic interactions, the lattice sum (LS) method (Hünenberger 1999;Sagui and Darden 1999) is often used in many studies. The periodic boundary condition (PBC) is assumed and the interactions from the infinitely expanded lattice is considered. The Ewald method (Ewald 1921;de Leeuw et al. 1980) is one of the representatives of LS methods. It decomposes the energy into a short range part working around the MD cell and a long range part expanded on the lattice, for which the former is treated by cutoff and the latter is treated as Fourier expansion. A detailed review on the LS method is given by Cisneros et al. (2014), and the recent development has been given in literature (Tornberg 2015;Wang et al. 2016a;Shamshirgar and Tornberg 2017;Hammonds and Heyes 2022).
However, the LS method takes into account the interaction of particles infinitely far from each other during a non-zero unit timestep t of MD, which means that information beyond the speed of light is being exchanged. In addition, the concept of an infinitely extended lattice itself does not actually exist in reality. Thus, it is necessary to remind that the LS method is only an approximation for mimicking the nature.
In this paper, we focus on the non-Ewald method for calculating electrostatic interactions, which is also called the real-space method or the cutoff-based (CB) method . The most simple CB method is the straight truncation method that calculates the electrostatic energy of N charged particle system (the non-SI unit) by where r c is a cutoff length. In the CB method, in general, the interactions are computed by a finite pairwise sum, where U ij is a pair potential function considered in a physically reasonable manner, based on the Coulombic function q i q j /r, and the second term should be a certain constant. Now, what is the artifact of the CB method? The often observed very large magnitude of artifact in the straight truncation method is due to the discontinuities of the pair potential and pair force functions, as discussed by . Then, what is the artifact of the CB method endowed with continuities for those quantities? It is often thought that the artifact is caused by the truncation of far distant information (Allen and Tildesley 2017), but to what extent is this argument correct?
In fact, we know that the effects exerted from afar are not instantaneous by the cutoff. However, as mentioned above, there is no need for such instantaneous coverage of distant information. In effect, if we take in information within a space, e.g., r c = 10Å = 10 −9 m during an MD timestep t = 10 −15 s, we are taking in information at a speed of 10 6 m/s, or 1/300 of the speed of light. The effect of the distance over r c is not instantaneous due to the cutoff, but it is not blocked forever, but comes with a time delay through the medium of particles that exist in between. For example, the information of a particle at n times the cutoff distance can be delivered by the mediation of (n − 1) particles at the shortest distance, so that the information at 10 times the cutoff distance r c , which is about 100Å away, can be transmitted at the fastest speed of about 10 t = 10 −14 s. This transmission speed is thus the same as before, 10 6 m/s, and the transmission would be robust via the surroundings of each particle in a crowded particle system, such as a biomolecular system and condensed matter system. We would like to consider whether this speed and transfer method of information in CB methods really produce truncation artifacts.
The theme of this paper is to discuss the sophisticated CB method based on a suitable pairwise finite sum. However, since we have already reviewed the work presented earlier , this paper mainly reports what has been clarified since then. First, one of the attractions of the CB method is that it is often much faster than LS. Another reason for expectation is that the artificial regularity, PBC, is not necessarily needed. In this paper, we will review in detail the theoretical aspects such as the latter issue rather than the practical aspects of the former issue, i.e., the computational speed. More specifically, the discussion will focus on the reaction field (RF) method (Barker and Watts 1973;Hünenberger and van Gunsteren 1998;Schulz et al. 2009), isotropic periodic sum (IPS) method (Wu and Brooks 2005;2008;2009;, and zero multipole summation method (ZMM) Fukuda et al. 2011Fukuda et al. , 2014, on the basis of the following reasons.
The RF method has been expected as a method to compensate for the shortcomings of the straight truncation method. It has a long history (Onsager 1936) and is still one of the representative CB methods. This is the primary reason why we focus on this method in this paper. The characteristic of this method is that the far state is considered as a dielectric media and the interaction from the dielectric is taken into account, wherein the physical and chemical background is well established (Fröhlich 1958). For this reason and the computational efficiency, the RF method has been used in many studies as an alternative to the LS method. In fact, it has been implemented in GROMACS (Berendsen et al. 1995) and applied to many biomolecular systems (Schreiber and Steinhauser 1992;Norberg and Nilsson 2000;Mulakayala et al. 2013;Soares et al. 2016).
The artifact caused by the group-based cutoff in the RF method was reported in detail (Ni and Baumketner 2011;Wang et al. 2016b). Refer to our previous work  for the artifacts that had been pointed out before that. In this paper, we reconsider the RF method by adding a new theoretical point of view. For example, the relationship between the RF and ZMM methods has already been discussed , but here we consider a deeper relationship than that in our previous article.
The IPS method, developed by Wu and Brooks (2005) and implemented in CHARMM (Brooks et al. 1983), incorporates distant information with a different idea from the RF and LS methods. Specifically, assuming homogeneous nature, local states are copied non-regularly to represent distant states. The advantage of this method is that it can avoid introducing artificial translational symmetry, which is required in the LS method. In addition, unlike the RF method, the amount of the external information, such as a permittivity parameter, is not included in the final energy expression. Because the IPS method is a promising method that has been updating its theory and description with various ideas, it is discussed in this paper. The applications of the IPS method are found in literature (Klauda et al. 2007;Venable et al. 2009;Takahashi et al. 2010;Nakamura et al. 2013;Ojeda-May and Pu 2014b;Nozawa et al. 2015aNozawa et al. 2015b. ZMM is a method that we and our collaborators have developed Mashimo et al. 2013;Fukuda et al. 2011;Kasahara et al. 2016;Sakuraba and Fukuda 2018). Unlike the LS, RF, and IPS methods described above, the ZMM does not actively consider the incorporation of distant information and argues that the artifact comes from the ad hoc computation of all interactions within the cutoff sphere, rather than from truncating the distant effects. ZMM gives an idea of which of the interactions in the cutoff sphere should be excluded. To achieve this, it is necessary to switch the viewpoint from the details of the energy function to the domain of definition of the function. The basic concept is that local electrostatic neutrality should always be taken into account for calculating the interactions. ZMM is implemented in myPresto (Mashimo et al. 2013), omega-gene (Kasahara et al. 2016), and GROMACS (Sakuraba and Fukuda 2018). Its applications to biomolecular systems are found in literature Kasahara et al. 2014;Nishikawa et al. 2014;Iida et al. 2016;Kamiya et al. 2016;Nishigami et al. 2016;Bekker et al. 2017;Kasahara et al. 2017;Kasahara et al. 2018b;Hayami et al. 2019;Higo et al. 2020;Hayami et al. 2021;Ikebe et al. 2021;Higo et al. 2022).
Although the ideas on which the above RF and IPS are based are different, they have in common that they are "addition methods" that try to incorporate influences outside the cutoff sphere. On the contrary, ZMM is a "subtraction method" that tries to remove the artificial influence near the boundary from the cutoff sphere. In order to make such a contrast, we focus on the ZMM in this paper and discuss necessary principle and mathematics for this comparison. Since "addition methods" and "subtraction methods" are based on completely different ideas, they may seem to give completely different results (Allen and Tildesley 2017). However, they are in fact closely connected, as shown in this paper by revisiting the physical concepts and mathematical structure of the individual methods.
In "Physical concepts in non-Ewald methods," we review the RF method, PBC, the IPS method, and the ZMM. We derive the RF (Appendix B) and IPS methods (Appendix C), on the basis of the original work, by adding our original interpretations. This helps to keep the consistencies of the descriptions, capture the essence, and simplifies the understanding. We give a brief derivation of the ZMM in Appendix D. In "Relationships among non-Ewald methods," we then discuss similarities among the non-Ewald methods, utilizing mathematical common feature. Within this review, we synthesize the following concepts: (1) derivation of the modified RF (MRF) by neutralization principle. The notion of "cavity" conforms with the notion of the neutralized subset, which is taken into account in the ZMM. (2) A direct relationship between the ZMM and IPSp, which is one of the extended IPS methods discussed in "Isotropic periodicity principle." Although one may think the Wolf type method and the IPS method are totally different, the ZMM can bridge over these methods. (3) Approximate relationship between the generalized RF (GRF) and ZMM. In "Conclusive discussion," we briefly summarize the results and give discussions on them.

Physical concepts in non-Ewald methods
We consider electrostatic interactions among N atoms (in this paper, atom means a classical point particle) with charges q 1 , . . . , q N and coordinates r 1 , ..., r N , which we will simply denote by x ≡ (r 1 , ..., r N ) ∈ R 3N . For ease, denotes all the atoms in the cutoff sphere with the radius of a given cutoff length r c centered at the position of atom i, and denotes those atoms except i. Here r ij ≡ r ij ≡ r i − r j is the distance between atoms i and j . Thus, Eq. 2 can be simply rewritten as Note that the second term represents a constant that is independent of charge configurations in an ideal sense but sometimes allowed to be dependent if it could be approximated as a constant.

Reaction field principle
As stated in the "Introduction," the RF method (Onsager 1936;Barker and Watts 1973;Watts 1974;Neumann 1983) has a long story and has been studied with various applications (Alper and Levy 1989;Schreiber and Steinhauser 1992;Gil-Villegas et al. 1997;Heinz et al. 2001;Baumketner and Shea 2005;van der Spoel and van Maaren 2006;Míguez et al. 2010;Reif and Oostenbrink 2014;Pechlaner et al. 2022). The RF method takes into account the interactions between each molecule involved within a "cavity" and the environment outside the cavity, as well as the interactions between molecule pairs inside the cavity. Here, "cavity" of each molecule a is regarded as a set of molecules included in the cutoff sphere with radius r c (or often referred as the cutoff sphere itself) in the conventional idea. Note that, while it resembles the conventional idea, we will provide below an alternative sophisticated definition about this idea. In any case, the region outside the cavity is assumed to be a dielectric continuum with a dielectric constant RF . The continuum polarizes via reacting with the molecules inside the cavity, and the polarization generates an electric field (reaction field) which is proportional to the total dipole in C a , the cavity of molecule a (Fröhlich 1958;Allen and Tildesley 1987). Here μ b ≡ i∈Mol b q i r i is a dipole of molecule b in the cavity. The dipole in the cavity reacts with this field E a and the energy − (μ a |E a ) is produced. Thus, the total energy of M molecules considered in a system becomes Here the suffix "Out" means that this energy represents the interactions that comes from the environment outside each cavity. Several derivations of the RF method are known. In convention, interactions between molecules inside the cavity are treated as a simple truncation within the cutoff sphere with radius r c , where the system consists of N atoms, and the total energy of this system is estimated as where is a pair potential function of the RF method.
One issue is the discontinuity of the energy function (9), which comes from V RF (r c ) = 0 for RF = 0. As discussed in our previous work , this causes serious problems in molecular simulation, and for a remedy, we have heuristically proposed a modified RF (MRF) summation method, which is defined by the energy formula, where V RF is the same as in Eq. 10. As well as obtaining the continuity, the self-energy term arises in the last term of Eq. (11).
In this paper, we give justifications of Eq. 11 via two routes. One of the routes is a practical one, where E RF (x) itself is used and an approximation described in Appendix A is applied to it. This means subtracting the boundary term or adopting the charge neutrality condition for E RF (x), and we immediately have Eq. 11 by substituting V RF into V (see Appendix A, (81)). Another route is more systematic, which is detailed in Appendix B. It is based on redefining the idea of "cavity" as stated above (see Fig. 1), for which, simply speaking, the cavity can be treated as a charge neutralized subset, as will be introduced in "Neutralization principle." Following this idea, starting from a quantity that corresponds (but not equivalent) to Eq. 8 and from expression (7), it leads us to Eq. 11 in a simple and straightforward manner, which is not through (9). It thus gives an alternative proof of energy formula of the RF energy, which is not exactly the original one but its analogue.
Even though the RF or MRF is supposed to incorporate the external information outside each cutoff sphere, the method is regarded as a CB method because it is written in a form of Eq. 5 using the internal information inside the cutoff sphere along with the data RF . An issue is that the only information outside the cutoff region, RF , is often unknown a priori. It is also noted that the force is not continuous unless RF → ∞. Following extensions for the RF method (b) New interpretation of cavity: cavity is the set of all atoms of solid ellipsoidal molecules, and it is common for atoms i and i . Open ellipsoidal molecules, including c, d, e, f , and g, do not affect both i and i . A molecule e, which is supposed to be totally charged, is excluded from the cavity; otherwise, it will violate the charge neutrality condition for the cavity (condition (b') in the text). A molecule f , which is supposed to be charge neutral, is also excluded from the cavity; otherwise, the partial charge will violate the charge neutrality condition for the cavity (Perram and Smith 1987;Tironi et al. 1995;Hünenberger and van Gunsteren 1998;Heinz and Hünenberger 2005;Lin et al. 2009), a simple approach to remedy this force discontinuity by shifting or switching for function V RF is recently investigated by Kubincová et al. (2020).

Periodicity principle
An alternative to the RF method for taking into account the effect of the environment surrounding the target system involves the use of LS methods which impose PBC. We first review here concepts associated with PBC in order to better introduce another periodicity principle known as isotropic periodicity.
In a typical PBC, the system is defined by an atom set distributed in a cubic box, referred as a "basic cell" (or MD cell), surrounded with its copies that are arranged along the three orthogonal direction under the space-filling and disjoint manner in R 3 . The energy of a charged atom i can be defined as follows: where r ij n = r i − r j + Ln with n ∈ Z 3 − {0} indicates three triplet integers except n = (0, 0, 0), with L being the cubic box length, and ε ij (r ij ) (i, j ∈ N ) is a suitably defined pair potential function value such that the series in Eqs. (12a) and (12b), i.e., n∈Z 3 −{0} ε ij (r ij n ), converge for is utilized for the Coulombic interaction in the Ewald method real part (pure Coulomb potential is conditionally convergent). The total energy is obtained by 1 may ordinarily be interpreted as the interaction energy between i and all other atoms, which include every atom except i in the basic cell and include every atom in every image cell (or called PB cell) (see Fig. 2a), an alternative interpretation can be given as follows. First, since r ij n (j ∈ N i ) can be interpreted as the distance between atom j in the basic cell and atom i in the image cell n, E PBC ij can be viewed as the interaction energy between j and all "i"s that are composed from the original i (viz, i in the basic cell) and image "i"s (viz., original i's copies positioned in image cells). Thus , which appears in the first term of the RHS of Eq. 12b, indicates the interaction energy from original i and all image "i"s to the all atoms (indexed by j ) in the basic cell except i. Second, the so-called the selfinteraction n∈Z 3 −{0} ε ii (r iin ), which appears in the second term of Eq. 12b, indicates the interaction energy from all image "i" s to the original i. Thus, E PBC i (x) in Eq. 12b can be interpreted as the interaction energy from original i and all image "i" to the all atoms in the basic cell (see Fig. 2b). In other words, E PBC i (x) is the energy felt by the system from all "i" (where the "system" is composed from all atoms in the basic cell). This is also true for Eq. 12c. However, the way the components are divided differs  (c) correspond Eqs. 12b and 12c, respectively: E PBC i (x) is the energy felt by the system from all "i"s (i in the basic cell and all image "i"s). Equation 12b represents the energy felt by "j "s in the system plus the energy felt by i in the system, while Eq. 12c represents the energy influenced from i in the system and the energy influenced from "i"s outside the system between Eqs. 12b and 12c. Equation 12b decomposes the energy into the one felt by the system other than i and the one felt by i in the system, while Eq. 12c decomposes the energy into the one influenced from i in the system, j ∈N i ε ij (r ij ), and the one influenced from "i"s outside the system, Fig. 2c). These alternative viewpoints are one of the key points to understand the isotropic periodicity, as demonstrated below.
The PBC can be generalized with using a space filling shape, instead of the cube, and the copies are distributed under the space-filling, disjoint, and periodic manner. In any type of PBC, a certain symmetry arises where the basic cell is copied. Thus, the physical property in the basic cell should be governed by this symmetry. This is not physically sound in systems such as liquids (Kurtovié et al. 1993;Wu and Brooks 2005).

Isotropic periodicity principle
To solve the problem that the PBC is not physically sound in systems such as liquids, the IPS method does not feature the PBC, but uses the isotropic boundary condition (IBC). The method utilizes, instead of the basic cell, a sphere around each atom i with a radius r c , for which the sphere shape should be harmonic with the isotropy and the sphere is called the "local region" of i. The method then distributes copies of the local region, where the copies are indexed by m ∈ N, under an isotropic manner in which the space filling and disjoint property are not necessarily met (see Fig. 3).
. Shown is only for m = 1, 2. Dotted circles represent copies of the local region (called "image regions"), while, for calculating the interactions, the continuum approximation 2 , · · · are actually used instead, which can be viewed as continued images of i , i , · · · , respectively (see Appendix C (ii)). The "spherical interaction" (or "random interaction") is interactions to atom j from spheres S Each copied region is called an image region. The energy of atom i is defined as follows: Here, we restrict our discussion in the case of but the IPS method can be formulated for other interactions; R i and R i are given in Eqs. 3 and 4, respectively. Note that there is a case that a factor 1/2 is attached for Eqs. 13b and 13c (e.g., Brooks 2005, 2019) for formulating theories, while not attached (Ojeda-May and Pu 2014a) for calculating Madelung energies, wherein we follow the latter definition, which is also parallel to Eqs. 12b and 12c in the PBC. The index m ∈ N labels the radial direction centered at r i , with letting R m = 2mr c , which captures the radius of concentric spheres. The image regions are placed inside each m-indexed spherical shell These shells decompose the total space such that 0 be the local region. On the other hand, the spherical directions are captured as an integration of continuously distributed individual copy. ε ij includes all the effects such as this integration, as described below and detailed in Appendix C. m∈N ε ii (r ii ; m) in the second term of Eq. 13b represents self interaction.
The common feature between Eqs. 12 and 13 is clear, as interpreted by the correspondence: Namely, (i) the all atoms in the basic cell, N , in the PBC corresponds to all the atoms in the local region, R i in the IBC. Accordingly, N i corresponds to R i ; (ii) while the interaction energy between atom j in the basic cell and all the copies of atom i in the PBC is n∈Z 3 −{0} ε ij (r ij n ), it corresponds to the interaction energy between atom j in the local region and all the copies of atom i in the IBC, m∈N ε ij (r ij ; m); (iii) the self interaction in the PBC, which is obtained by replacing ij →ii in the left side of the third equation in (16), corresponds to the self interaction in the IBC, which is obtained by replacing ij →ii in the right side of the last equation in Eq. 16. On the other hand, the difference between Eqs. 12 and 13 is seen in: (1) For the PBC, the fundamental region (basic cell) and the atoms (N ) are common for all i, but for the IBC, the fundamental region (local region) and the atoms (R i ) depend on each i; (2) the IBC assumes the special dependence of the energy, ε ij (r ij ; m), which is on a separated contribution from distance r ij and index m.
As discussed in Appendix C, ε ij is defined through several procedures, and we have where ψ is the digamma function, and γ is the Euler-Mascheroni constant. Using the notation we have the total energy of the IPS method (see Appendix C), Here the right first term is a pairwise summation form with a cutoff procedure and the second term is called the selfenergy term. The last term is called the "boundary energy," whose approximation method is discussed in Appendix C (vi). Note that, for ease in the computational handling, several approximations of have been proposed to avoid the use of ψ Brooks 2005, 2008), such as IPS (r) being a polynomial (Wu and Brooks 2008).

IPSp method
The method introduced in "Isotropic periodicity principle" is derived based on a nonpolar homogeneous approximation, so that it is labelled as IPSn method. The IPSp method derives a polar electrostatic 3D-IPS potential (Wu and Brooks 2009). The pair potential function in the IPSp method is defined by adding a certain "screening" function Screen ij (r ij ) into the original IPSn pair potential function (18). In other word, the IPSp method is obtained by replac- and we have instead of E IPS ij (r ij ) (see Eqs. 13a and 18). In the IPSp method, Screen ij is defined by Screen and the coefficients b m s are determined by constraints: Here we have denoted as E A physical interpretation for Eqs. 21 and 24 is not so straightforward. Original work is based on considering that Screen (r) represents a screening effect by oppositecharge distribution around each charge in polar systems and that there are interactions with dipoles, which yield charge-dipole interactions, resulting in m = 2 condition in Eq. 24, and dipole-dipole interactions (m = 3), as well as charge-charge interactions (m = 1).
Three remarks are made.
(1) Note that Screen is specific to electrostatic interactions, whereas IPS ij can be constructed within the same procedure for any interactions such as the LJ type, instead of the electrostatic interaction.
(2) It should be stressed that Screen (r) is treated to work only as the axial interactions (see Appendix C (iii)), not as the spherical interactions. This is a characteristic choice, because a non-Coulombic function could contribute to spherical interactions (while pure Coulombic interactions cannot contribute to spherical interactions, as stated in Appendix C (v)). From this choice, a symmetry with respect to the line i − j arises, and the force between j and all "i"s in the line becomes zero when j is at the local region boundary (i.e., r ij = r c ). (3) Another conceivable idea based on a physical viewpoint regarding Eq. 24 is that the screening effect can be measured not only for a charge that is really placed at the site j but also for a multipole that is imaginary placed at the site j . Considering that the chargecharge interaction is described by E IPSp (r), an interaction potential and a force between a charge q i and an mthmoment multipole (m = 0, . . . , L − 1) are proportional to D m E IPSp (r) and D m+1 E IPSp (r), respectively. Here we have assumed that the multipoles have nonvanishing components only for the i − j direction (line multipole). This idea contributes to a possible supplementary explanation to justify Eq. 24.
The energy functions for the IPSp method are obtained by the replacement of IPS ij into IPSp ij . Thus, in the IPSp method, where we have used IPS

Neutralization principle
Each of the RF, IPS, and LS methods is an "addition method" that tries to incorporate the effects beyond the truncation region (which corresponds to the cavity, local region, and basic cell, respectively), whereas the ZMM discussed in this subsection is a "subtraction method" that tries to remove artificial effects from the truncation region (which corresponds to a cutoff sphere), where the artificial effects are generated near the boundary of the cutoff sphere. As discussed in "Introduction," we do not need to consider interactions from so far away that the rate of the information transfer exceeds the speed of light. In the most simple cutoff, i.e., the straight truncation, each charge i feels the interactions from every charge j only in the cutoff sphere with the radius of a given cutoff length r c : Then, is a cutoff length of 10Å too short? Given that the potential function value decays with the distance between atoms, is it more accurate to increase the cutoff length even a little? In fact, it has been shown that this is not the case; Yonetani (2006) shows that the longer the r c , the worse the results can be. What, then, is the essential problem with the CB method?
The reason is that the usual CB method takes in just all interactions within the cutoff length, which is pure distance clipping of interactions. A specific problem is that high-energy states, e.g., configurations gathering only positive charges, may occur just inside the cutoff region, which is artificially prepared. Although no such highenergy states can be directly observed with a slight change in the region surrounding such a state, such states are virtually present when applying the straight truncation ( Fig. 4). High-energy configurations (or also extremely lowenergy configurations) should be excluded, considering the canonical distribution at the room temperature, the most encountered situation in equilibrium molecular simulation. However, because the straight truncation incorporates the interactions from such a state, this can give false information to the central atom i. Moreover, immediate decay of such a state, as expected in a sufficiently wide region under the canonical distribution if it were generated, does not necessarily occur inside the clipping region.
The basic policy of ZMM (Fukuda et al. 2011; corresponds to discarding the virtual high-energy states arising just inside the artificial cutoff boundary. In other words, it is a method of "subtraction," i.e., not counting interactions from such states, and discarding them. Since these virtual high energy states derive from the violation of the local electrostatic neutrality, ZMM takes into account electrostatically neutralized configurations inside the cutoff sphere. The neutrality is considered locally and is acceptable even if it is approximately attained. In addition, the neutralization should not only be for the charges but also dipole and higher multipole if needed. That is, (nearly) neutralized configuration, which is called a neutralized subset, is utilized in the ZMM (Fukuda et al. 2011;: where M (l) i is the neutralized subset in the sense of lth-(multi)pole with respect to charge i (see Fig. 4). However, seeking the neutralized subset M (l) i may be time consuming in practice. Thus, the calculation of M (l) i is not required, but instead, just its existence is assumed. Our strategy is that we return to a simple cutoff sum form but instead deform the Coulomb pair potential function into a certain function U (l) in order to approximate the neutralized sum (28) such that To justify this idea and specifically formulate U (l) , it relies on the following three conditions for characterizing M (l) i : The first two conditions are just the ones explained above. The last one indicates that the charges ( = i) excluded from the neutralized subset M (l) i in the cutoff sphere of i are located close to the cutoff surface. In other words, the neutrality is broken near the cutoff surface that is the boundary of the ad hoc chosen perfect sphere shape (which by itself does not necessarily ensure the neutrality). Using these three conditions along with a consistency condition (which is just the action-reaction condition), Eq. 29 can be mathematically deduced .
To briefly review it, we first fix some notations: ε(r) = 1/r is the Coulombic kernel, and is an excess subset, which is a subset of charges whose interactions to charge i should be subtracted. That is, the total energy E(x) is represented as wherein the second line means that interactions from the excess subset. 1 2 i∈N j ∈J (l) i q i q j ε(r ij ) which is called the excess energy, is subtracted from the interactions from all charges in the cutoff sphere, 1 2 i∈N j ∈R i q i q j ε(r ij ). As reviewed in Appendix D, by using a suitable polynomial the excess energy can be approximated by a pairwise form for this polynomial, allowing the subtraction to be accomplished using a pairwise form, This formula accords with the RHS of Eq. 29 with U (l) (r) = ε(r) − ε l (r).
The above procedure utilizes the bare (undamped) Coulombic function ε(r) = 1/r in the starting point of Eq. 27. In general, the same procedure can be applied to a damped function V (r), by suitably treating the remaining part 1 r −V (r). Specifically, the energy function of the ZMM can be represented by Here, in the first term, is the damped Coulombic function, and the remaining function 1 r − V (r) contributes the last term of Eq. 35, viz., − α √ π i∈N q 2 i . In the pairwise term, is the polynomial whose degree depends on the multipole moment l, and the coefficients a (l) m are uniquely determined and depend on α, l, and r c . Note that α = 0 reduces E (l) ZM (x) into the RHS of Eq. 34 , and in particular for all l and m.
The argument that "subtraction methods have artifacts and addition methods do not" is incorrect. The first reason is that these criticisms are based on the results of older literature (e.g., Koehl 2006;Izvekov et al. 2008) which mostly refer to artifacts in the straight truncation or groupbased cutoff (Baumketner 2009) and do not refer to the artifact of the CB method in general, nor to the artifact of the subtraction method in the CB method. In fact, for example, ZMM, a subtraction CB method, shows that there is no unignorable artifact in polar fluids such as water systems Wang et al. 2016b), DNA in explicit solvent as a highly charged system (Arakawa et al. 2013), and a ferroelectric crystal as a dipolar system (Wang et al. 2016b). Rather, the method also helped to clarify the artifact in the PBC (Kasahara et al. 2018b). The second reason is rather conspicuous, that is, the subtraction method, ZMM, can be equal to an addition method, IPS method, as will be shown in the next section. Finally, it is also noted that the subtraction method idea allows us to obtain the continuity for energy functions (see Appendix A), but there is no known case where this can be achieved by use of an addition method idea.

Other methods
One of other promising electrostatic interaction calculation methods is the fast multipole method (Greengard and Rokhlin 1987), which uses the multipole expansion in a suitably divided cell, utilizing the feature that the pair potential function decays as the atomic length is larger. Its recent development is provided in literature ( A quantitatively different approach for calculating electrostatic interactions is to utilize non-Euclidean geometry, such as 3-dimensional sphere S 3 (Caillol 1999;Ramírez EV and Elvingson 2022), instead of R 3 . Recently, Clifford torus is also considered and applied to calculate Madelung energy of crystal systems (Tavernier et al. 2021).

Mathematically unified form
There is a common mathematical structure in the form of the energy function among the non-Ewald method stated in the previous sections. As stated in literature , the energy function follows the form: Here, the first term is the pairwise sum using a basic function u that defines a pair potential function which is at least continuous. The second term − 1 2 u(r c ) i∈N q 2 i is the self-energy term, which can be viewed as the i = j term arising inside the first term. The last term comes from the choice of the kernel as the Ewald real function (see Eq. 36) and works as a correction of the energy on this choice.
CB methods taking this form are summarized in Fig. 5, and the relationships among these methods will be discussed in the following subsections.

Relationship between the RF method and the ZMM
The energy function of the ZMM can be rewritten (see Appendix D) as where u (l) is the basic function represented by Thus, it exactly takes the form of Eq. 39 with u ≡ u (l) . The ZMM with l = 1 and l = 0 are the zero dipole (ZD) method (Fukuda et al. 2011; and the Wolf method (Wolf 1992;Wolf et al. 1999), respectively. The ZMM with l = 0 and α = 0 is also considered to be the same as the scheme by Harrison (2006) (see ). ZMM with α = 0 is indicated as non-damped ZMM (Nd-ZMM) in Fig. 5 the scheme modifying the Wolf method to use it in MD with a smooth force function using a switching length r 1 . The energy of the modified RF (MRF), Eq. 11, also takes the form of Eq. 39 with u ≡ V RF and α = 0. As pointed out before , we can see the basic function correspondence between V RF of the RF and u (1) of the ZMM (ZD), V RF | RF →∞ = u (1) α=0 . Note that an interpretation of the ZD method relating to the RF principle is discussed by Elvira and MacDowell (2014), where the dipole in the medium causes an electric field, which defines the force and its integration leads to the a potential function u (1) (r ij ) − u (1) (r c ).

Relationship between the IPSp method and the ZMM
We show an equivalence between the pair potential function of the IPSp method with a polynomial approximation and the pair potential function of the ZMM with α = 0 and l = 3. As shown in Appendix D, the pair function of the ZMM, u (l) for l = 3, is given by as well as Screen (r) (see Eq. 23); (ii) M = 3 is chosen (see Eq. 8 in Wu and Brooks (2009)). Corresponding to this polynomial approximation, replace E IPSp (r) by and ( Comparing Eqs. 43 and 46, we see the equality of the basic functions between the ZMM with l = 3, α = 0 and the IPSp method, i.e., for all r, and we also observe the equality of the pair potential functions of both methods (see the first term of RHS of Eq. 26; see also Eq. 19 in Wu and Brooks (2009)): Now we consider the reason of the correspondence, (47). In the ZMM with α = 0, to obtain the lth-order approximation of the excess energy, we require the following condition (see Eq. 33 and Appendix D): m } m=0,1,...,3 , and thus uniquely solved . While for the IPSp, the method to determine polynomial coefficients via Eq. 24 is rewritten by the strategies (i)-(iii) such that where we have put and defined d 0 by Considering the fact that the pair potential function of the Eqs. 37 and 38 and Appendix D,Eq. (116)) and that of the IPSp method is given by Eq. 52, they are the same owing to Eq. 56 . Thus, we have obtained the proof of Eq. 48a. Equation 47 follows from this and from the fact that E IPSp (r c ) = u (3) | α=0 (r c ), which is derived from Eqs. 54 ,56,38,and 115. To summarize, between the IPSp and the ZMM with α = 0 and l = 3, the equations determining the polynomials that define the common pair potential function form are mathematically the same. Recall that, however, the physical motivations for the polynomial determination are different between the two methods: the ZMM with α = 0 and l = 3 motivates to have the third order approximation of the excess energy error, while the IPSp motivates to realize the IBC and adapt to polar systems. It should be noted that the equivalence, (47), is owing to strategies (i)-(iii) stated above. If we do not choose these strategies, the exact equivalence is not expected to be obtained. In fact, (i) is not chosen by Wu and Brooks (2005), for which IPS (r) is not approximated as a pure polynomial (see Eq. 58 of Wu and Brooks (2005), which will also be shown in below Eq. 61, for the approximated form of ε(r) + IPS (r), and see also Eq. 38 of the literature for other potentials) (ii) is not chosen in (Wu and Brooks 2008), for which M = 6 is chosen (see Eq. 12 in Wu and Brooks (2008)). In general, to improve the accuracy of the approximation to the original pair function in numerical fitting, it will generally be advantageous to increase the number of coefficients involved in the approximate function (e.g., by increasing a polynomial degree, L or M). In fact, in IPSn, the deviation from the original function is measured by RMSD (Wu and Brooks 2005). In the IPSp, in contrast, the degrees of freedom of the fitting is lost by ignoring the higher order terms of the polynomial.
Regarding the energy function of the ZMM with α = 0 and l = 3 (see Eq. 41), and the approximate energy function of the IPSp (cf. (26)) is Their pairwise summation terms and the self terms are the same as seen above. The difference between these two energy formulas is in the last term of Eq. 58, but these two will be the same if we use an approximation of subtracting the boundary term or adopting the charge neutrality, as detailed in Appendix A, noting that E IPSp (r c ) = 35/16r c = 0 owing to Eq. 46. Namely, if we denote an energy of a modified IPSp method by E MIPSp (x) that is defined, using the above approximation, by then it is the same as the ZMM energy with α = 0 and l = 3. It is reported in the study of a model liquid-crystal system (Nozawa et al. 2015b) that IPSn and IPSp methods indicate similar properties with that of the particle mesh Ewald method (Essmann et al. 1995) in a range of temperatures except for away from equilibria.

Other relationships
To compare pair potential functions in different methods, u − u(r c ), as employed in Eq. 39, it is often convenient to use, instead of u : where s = r r c . This is because the original functions u often contain the parameter r c so that r c dependence should also be studied for their comparison, while functions v, which are scaled but completely capture the original functions u, are all parametrized by s ∈ (0, 1) in a unit scale and may not contain r c . In fact, the function v in the following methods, including the IPS, RF, and ZMM with α = 0, does not contain r c .
(1) IPSn is similar to the RF with low RF -.
As indicated by Wu and Brooks (2005), the basic function of the IPS method, E IPS , approximated such that where both scaled functions are independent of r c . Figure 6 indicates (2) GRF and ZMM with α = 0 and l = 4 are similar-.
The generalized RF (GRF) method by  is derived from an attempt to eliminate the asymmetry in the RF method by assuming that the interacting distant atoms are also surrounded by a uniform background charge, and it has the following pair potential function: Note that u GRF is based on V RF ((10)) with RF → ∞ and that r c in Eq. 64 corresponds to twice the cutoff length for V RF , as used in e.g., bulk water (Teplukhin 2018), ice-water interface , ion solvation (Soper and Weckstrom 2006), and Amyloid β peptide (Gnanakaran et al. 2006). Here, u GRF (r) = 1 r + 4 r 2 r 3 c − 3 r 3 r 4 c + 2 5 r 5 r 6 c so that its scaled function is v GRF (s) = 1 s + 4s 2 − 3s 3 + 2 5 s 5 .
The pair potential function of the IPSM3 (IPSMm with m = 3) with a polynomial approximation and that of the ZMM with α = 0 and l = 4 is also shown to be equal. Here, the IPSMm method (Wu et al. 2016) defines "background multipole IPS potential"ˆ BGMm and it is subtracted fromˆ IPS (r ij ) to represent a pair potential functionˆ IPSMm , i.e., which is analogous to Eq. 21 with replacement of Screen into −ˆ BGMm but different in that it should be zero at r ij = r c . An approximated pair potential function of the IPSM3, E IPSM3 , becomes (see order 3, Table II  This is an approximation ofˆ IPSM3 (r) obtained in a similar manner stated for the IPSp method in "Relationship between the IPSp method and the  (60)), which completely characterize the pair potential functions and become irrelevant to the cutoff length, are shown for the IPS method, the MRF method with RF = 5 (very similar to IPS), that with RF → ∞, the GRF method, ZMM with α = 0 and l = 4, ZMM with α = 0 and l = 3, and the general shifting method (GSM) with l = 1, 2, and 3. Panel (a) indicates for almost full scale 0.1 ≤ s ≤ 1 and (b) for 0.8 ≤ s ≤ 1 ZMM." Namely, E IPSM3 (r) is constructed from a polynomial approximation with the type of Eq. 52 within the replacement of 3 by 4 in the sum and from considering the boundary condition corresponding to Eq. 55 (replace 3 → 4). Thus, from the discussion in "Relationship between the IPSp method and the ZMM," E IPSM3 (r) should be equal to the pair potential function of ZMM with α = 0 and l = 4 (see Eq. 119 in Appendix D). However, note also that this correspondence is due to the strategy of the polynomial approximation with approximation with the specific choice of the degree in the IPSMm, as stated in "Relationship between the IPSp method and the ZMM." (4) PA method: averaging principle-.
The pre-averaging (PA) method (Yakub and Ronchi 2003) connects the Ewald and non-Ewald methods (Fig. 5). The energy formula comes from averaging the quantities over the spherical angular coordinates expressed in the Ewald method, wherein a more mathematical discussion on the derivation appears very recently (Demyanov and Levashov 2022). The total energy can be represented by Eq. 39 with α = 0 and u ≡ V PA , where with r m being the radius of the volume-equivalent sphere of the cubic basic cell. It is shown ) that V PA = u (1) α=0, r c →r m , leading to the connection between the PA method and ZMM (Fig. 5).
Each of the above methods, including the RF, IPS, and ZMM, is built on a solid physicochemical background, and their relationships have been discussed so far. On the other hand, there may be a standpoint that the physical interpretation of the methods is not always necessary for their approximate computation (e.g., AI-like tools are useful for practical use). As a specific example, we can mention the general shifting method (GSM), which is a practical method, but in fact there are many application studies using a particular GSM, the shifted-force method, which has been developed and intensively studied by Fennell and Gezelter (2006). We will discuss the relevance of these methods to the ZMM.
As seen from Eq. 35, the pair potential function of the ZMM takes the form of V minus a function V l whose Taylor expansion around r c is the same as that of V up to the lth order (as shown by Eq. 107 in the α = 0 non-damped case; see  for a general case). A similar approach when focusing only on this point is taken in the GSM (Levitt et al. 1995;McCann and Acevedo 2013). In that method, the energy is defined as the sum of the pair potential function 1/r − g l (r), where g l (r) is the Taylor expansion of 1/r around r c up to the lth order. This method is possible to attenuate the straight-truncation method artifact often appearing in the vicinity of r c , such as that in the radial distribution function (Kale and Herzfeld 2011). Then, if we ask whether the pair potential function 1/r − ε l (r) of ZMM with α = 0 and the pair potential function 1/r − g l (r) of GSM are similar, it is true that they are similar at very near r = r c . However, they can be different, as r is shorter. Figure 6b shows the detailed view in a range 0.8 ≤ s ≤ 1, which represents the behavior relatively near r c , and that the curves, including ZMM with α = 0 and the GSM with l = 1, 2, are different. An exception is for ZMM with α = 0 and l = 4 and the GSM with l = 3, wherein they are similar in 0.8 ≤ s ≤ 1. However, they differ for shorter s, as clearly seen in Fig. 6a. It is a reflection of the variety of C l functions that, unlike analytic functions, they can have the same behavior in the neighborhood of one point (r c in this case), but globally they can be completely different.
What is characteristic of ZMM is the establishment of Eq. 111 (see Appendix D), which is deduced from the physical viewpoint represented in Eq. 30b. One of the differences between the GSM and the ZMM is that the polynomial that establishes (111) is chosen in the ZMM, instead of the simple Taylor expansion polynomial of 1/r. In addition, because of the presence of Eq. 111, the second term of Eq. 34 (i.e., the self term) arises in the ZMM. Such a term is necessary to obtain accurate estimates of the Madelung energy and is important for energetic properties, but it is not derived from the GSM alone (even if the method in Appendix A is used) and is often ambiguously treated. This is another difference between the two methods. Numerical comparison (Elvira and MacDowell 2014) on the accuracy between the ZD method (which is ZMM with l = 1) and the Fennell-Gezelter shifted-force method (which is GSM with l = 1) indicates that the former method is comparable or superior than the latter method in crystal systems. For more specific features of the shifting method, the reader can refer to Brooks et al. (1985), Steinbach and Brooks (1994), and Lamichhane et al. (2014).

Conclusive discussion
We reviewed non-Ewald methods, emphasizing the RF, IPS, and ZMM, with regard to their physical background concepts in relation to their derivations. We also discussed their physical and/or mathematical relationships and found new results. Our synthesis has shed light on a technical point of view that has not been discussed clearly so far.
It is not appropriate to criticize "subtraction methods" such as the Wolf method and ZMM by claiming that they have artifacts because they ignore distant contributions, and to praise "addition method" such as RF and IPS methods as if they are saviors because they are supposed to incorporate distant effects in order to avoid such artifacts. In fact, the present result that the ZMM and IPS method can be consistent is the antithesis of such a claim.
To take a more constructive view, we would like to explore why these "addition method" (RF and IPS) and "subtraction methods" (ZMM) are closely related or indeed why they are equivalent for certain parameter values. First, the distant influence brought about by the addition method is considered to enter the interior of the cutoff sphere through its boundary. On the other hand, the subtraction method, ZMM, manipulates the information on the boundary to consider the neutralization condition. This manipulation may actually be a reflection of the distant influence. It may be a future theme to investigate whether this can be quantitatively expressed.
On the other hand, the LS method tries to incorporate the distant influence under the PBC. However, the PBC is only a model, where the PBC is close to the reality in some systems but not in others (see (Kasahara et al. 2018a) and the references therein). Before we become obsessed with the PBC and try to reproduce the energy under the PBC as precisely as possible, we should not forget the points made in relation to the non-Ewald methods discussed in this review.

Appendix A: Continuity of energy function
Consider an atomic energy function E i defined by a cutoff sum (see Eq. 4) of a continuous pairwise function V ij = q i q j V such that For energetic studies and for keeping energy conservation, e.g., in NEV simulations, the continuity of the energy function is required. We will demonstrate that this requirement can be met by subtracting the boundary term (which is the last term of Eq. 72; see below for detail) or adopting the charge neutrality condition. Specifically, we state that they are equivalent, i.e.,

Subtracting the boundary term =
Adopting the charge neutrality condition.
In the following, we assume V ij (r c ) = 0, because V ij (r c ) = 0 makes E i continuous with respect to x so that the requirement is already satisfied. Obviously, the continuity of E i indicates the continuity of the total energy function To access the problem, just using a simple equivalence, we observe where the relation (4) has been used in the second line. The continuity of E i and so a continuous time development t → E i (x(t)) seem to be guaranteed by the last line of Eq. 72, but they do not hold. This is because the quantity j ∈R i q j in the final term, often called the "boundary energy" term, is a fluctuating quantity, since R i depends on the configuration x (viz., R i = R i (x) by definition (3)), so that the sum over it, j ∈R i (x) q j , is not ensured to be constant in time.
To proceed ahead, one may consider using a uniform/homogeneous approximation (Wu and Brooks 2005) in the boundary energy term such that for any i, where V is the volume of the basic cell and V 0 is the volume of the cutoff sphere. Another possible approximation may be for any i, which implies a subtraction of the boundary term in the original energy function, viz., The two approximations (73) and (74) are equivalent (Ojeda-May and Pu 2014a) for a large system with V 0 V or for a charge neutral system with According to either approximation, we get where the constant is −V (r c )q 2 i + V (r c ) V 0 V q i j ∈N q j or −V (r c )q 2 i , and hence we reach the continuity of E i . However, considering the neutralization principle ("Neutralization principle"), we see that Eq. 74 does not often hold and large deviations from it are often observed, as Wolf et al. (1992) deeply investigate this issue in ionic systems under condition (76). Therefore, the accuracies in the above approximations would become problematic. Namely, these are approximations that transfer the discontinuity from the original cutoff sum j ∈R i V ij (r ij ) into the boundary energy term and ignore the latter as a constant. In other words, approximation (74), or (73) for a large/neutral system, allows large intrinsic errors against the originally defined energy E i (x), in exchange for the gain of the continuity of E i .
Nevertheless, investigation (Ojeda-May and Pu 2014a) of the Madelung energy of an ionic crystal (in the IPS method) suggests that the accuracy of E i (x) becomes better if we take the approximation (74), which should have a large error against the originally defined energy E i (x). A possible explanation for this contradiction is that an error that is intrinsically involved in E i (x) can be extracted by removing the boundary energy. In fact, this explanation can be justified by the charge neutralization principle (Wolf et al. 1999), as follows.
The neutralization principle in the lowest order, discussed in "Neutralization principle" (which is recommended to be read to continue reading below), suggests the replacement of N i into the neutralized subset M (l) i with l = 0 for calculating the interactions, indicating that E i (x) should be replaced by We see that (30a) and Eq. (31)) Equation 79 indicates the subtraction approximation (75), so that the statement (71) is proven. On the total energy, we have Here, the final term of the second line is still not necessarily zero even for a neutral system, but it can be treated as zero in the approximation (75), or in adopting the neutrality condition, to get Finally note that adopting the charge neutrality condition, or the approximation with subtracting the boundary term, makes the energy function continuous, but it does not ensure the force function's continuity. The force function's continuity is ensured by, e.g., the ZM scheme with l ≥ 1 (Fukuda 2013).

Appendix B: Derivation of the reaction field method
We give a simple and systematic derivation of the RF energy formula, Eq. 11, using Eq. 7. For this purpose, we reinterpret the "cavity" using the notion of the charge neutralization, which turns out to be an idea that is really in harmony with this.

A. Cavity
Here, the "cavity" of molecule a, C a , is formally treated as a set of atoms, but not as a set of molecules (for the latter, a resembled notation C a is used in "Reaction field principle"), and the notation in Eqs. 3 and 4 is followed ("Neutralization principle" should be read before reading below). In conventional simulations, cavity is often considered as atom-wise to be the all atoms contained within the sphere of radius r c centered at atom i that is contained in a molecule a (see Fig. 1a).
The new interpretation of the cavity is given molecularwise and also in a manner that enables atom-based interaction treatment, as follows. Cavity C a is a subset of atoms contained within the sphere of radius r c centered at an atom i in a molecule a, so that it can be denoted as M i , viz., This subset M i needs not necessarily contain all the atoms of R i but is assumed to satisfy the following conditions: (b) The total charge of these atoms is zero, j ∈M i q j = 0, and (c) Atoms in R i but excluded from M i are placed near the boundary of R i .
Note that the conditions (a), (b), and (c) imply that conditions (30a), (30b), and (30c), respectively, hold for l = 0 when we view M i as However, for M i to be defined as a subset associated to each atom i and for the viewpoint of Eq. 82 to be appropriate, it is necessary that M i = C a is always valid, even when any other atom i ( = i) ∈Mol a is chosen. We will see this is justified in the following subsection. Before doing so, we specify the definition of interactions. The interaction on individual atom i in a molecule a is considered on the basis of this cavity C a (rather than the basis of R i as before), coming from inside C a (by its whole atoms) and from outside C a (by the electric field generated from a dielectric continuum). That is, we consider the interactions acting on atom i ∈Mol a as (i) an energy E Int coming from inside the cavity C a and (ii) an energy E Out coming from outside C a .

B. On (i)
By the definition, the interaction of type (i) acting on an atom i ( = i) in a molecule a is only from inside C a , wherein C a does not necessarily coincide with R i that is the total atoms in the cutoff sphere centered at i (Fig. 1b). Specifically, we assume that a molecule a is sufficiently small such that (a ) C a is a subset of R i (c ) Atoms in R i but excluded from C a are placed near the boundary of R i . Note also that it holds from above (a) and (b) that (b ) j ∈C a q j = 0. These conditions (a'), (b), and (c') also imply that conditions (30a), (30b), and (30c), respectively, hold for l = 0, when C a is treated as M i . That is, we reach a consistent characterization of C a such that for any molecule a. Note that it is should be C a ⊂ ∩ i∈Mol a R i , but not necessarily C a = ∩ i∈Mol a R i ; Fig. 1b illustrates the case that the equality does not hold. Hence, M i is defined with respect to every atom i, and conditions (30a)-(30c) with l = 0 hold for any i, where e.g., condition (30b) is denoted by We see that the concept of the neutralized subset is compatible with the concept of the cavity of molecule and will also see below that it fits well with the RF principle. Now, the energy E Int , coming from inside C a , is thus represented as

C. On (ii)
The energy E Out , which is coming from outside C a , is described by an electric field (reaction field) E a defined by Eq. 6, which is now represented, using the constant γ ≡ 2( RF −1) where any i ∈Mol a can be used in the most right equation (see Eq. (83)). Thus, where the last line is obtained as follows: To get Eq. 88c, we have used Eq. 31 and the consistency condition ) for l = 0; and to get Eq. 88d, we have used Eq. 84.

D. Total energy
Last, combining Eqs. 85 and 87, we have where V RF is defined by Eq. 10. Thus, a similar procedure to get Eq. 79 leads us to Eq. 11 such that

Appendix C: Derivation of the IPS method
We derive Eqs. 17-19, based on the original work (Wu and Brooks 2005) and our additional interpretation.
(i) First, the number of copies, n(m), of the local region distributed in a shell D (i) m should be defined (see Fig. 3). The number density is set to be irrelevant to m so that n(m)/V m = n(0)/V 0 = 1/V 0 , where V m is the volume of D (i) m , leading to the result n(m) = 24m 2 + 2.
By this setting, the isotropy along the radial direction could be met.
(ii) To attain the isotropy along the spherical direction, several strategies can be considered. A complete isotropy along the spherical direction cannot be directly attained by an arrangement, such as proposed in (i), of a finite number of spherical regions. Thus, a continuum approximation is defined: each copy of i in D (i) m is extended uniformly along the spherical direction. That is, the charge q i is supposed to exist not as a point charge, but as a uniform surface density ρ (m) i spread on a sphere S (i) m ≡ {r ∈ R 3 | r − r i = R m } (see Fig. 3). The potential energy, generated by this density, acting on the charge q j is where ε(r) = 1/r, and s ij m is the distance between r j and a surface area dσ on S (i) m . After this integration, we get which we here refer to the spherical interaction (rather than the "random interaction" (Wu and Brooks 2009)). This discussion can be generalized for any function ε ij (r) and for d-dimensional space R d , and this generalized ij (m) depends on r ij , in principle. Thus, the result, (93), is a special case for the Coulombic function in the three-dimensional space, known as the elementary fact that Coulombic/gravitational potential, which is produced by charge/mass uniformly distributed on the shell, acting on any point charge/mass inside the shell is constant. This simplicity, however, yields a difficulty to construct a method responsible for the IPS principle. That is, although one would like to straightforwardly use these good ideas by applying a definition ε ij (r ij ; m) = n(m) ij (m) to Eq. 13a, one finds a divergence in the series m∈N ε ij (r ij ; m) = q i q j m∈N (24m 2 + 2)/2mr c . Certain remedies are thus needed, as stated below (iii)-(v). (iii) The IPS method requires that there should be interactions that are symmetric with respect to the "ji axis," which is the axis including both r i and r j (i = j ). To meet this requirement, we put (at least two) image regions into each D (i) m along the j -i axis. So the image regions arrange, stick close each other, along the j -i axis (see Fig. 3). These image regions yield the energy term m∈N ε ij (R m − r ij ) + ε ij (R m + r ij ) , which is called the "axial interaction." Thus, a candidate of ε ij (r ij ; m) becomes (n(m)−2) ij (m)+ε ij (R m −r ij )+ε ij (R m +r ij ), (95) which recovers the dependence on r ij . (iv) Introduce a parameter ξ to ensure that the force at the distance of r c becomes zero, which is required for a stable MD simulation. Although we expect the axial interaction defined by Eq. 95 to yield the results that the force at the distance of r c is zero, it is not automatically ensured. To solve this problem, recall the strategy that we may duplicately put copies of the local region. Namely, we can put ξ copies along the j -i axis for the two opposite side, so that Eq. 95 will be changed into • ε ij (r ij ; m) ≡ (n(m) − 2ξ) ij (m) + ξε ij (R m − r ij ) +ξε ij (R m + r ij ).
The specific process to determine the value of ξ ensuring the force continuity is seen in (v), and ξ will be 1 for the three-dimensional Coulombic interaction. Note that ξ will not be an integer for a general interaction (e.g., LJ), so that its interpretation as the number of copies should also be generalized. (v) Regardless of the introduction of ξ , the sum of Eq. 96 with respect to m ∈ N does not converge due to the presence of the right first term. To remedy this problem, a certain regularization is introduced. Since ε ij (r ij ; m) should be a potential energy (acting on j from all copies of i), the energy zero can be freely chosen. We define then the first term of Eq. 96 vanishes, and the sum is convergent: where ψ is the digamma function and γ is the -Mascheroni constant. Equation 13a is thus E IPS ij (r ij ) = q i q j E IPS (r ij ) = ε ij (r ij ) + IPS ij (r ij ) = q i q j r ij − q i q j ξ 2r c ψ 1 − r ij 2r c + ψ 1 + r ij 2r c + 2γ .
Note that this result implies that all the spherical interactions (93) do not contribute to the IPS potential and that the only axial interactions (94) define the IPS potential. Using the properties of the digamma function, we observe that −D E IPS (r c ) = 0 is attained only if ξ = 1.
Thus, we have Eq. 18 along with Eq. 17. (vi) Continuity of the energy function is required in general, as detailed in Appendix A. The condition of the energy continuity in the IPS method, E IPS (r c ) = 0, is attained only if ξ = 1/(1 − 2 ln 2) (due to ψ (3/2) = ψ (1/2) + 2 = −γ − 2 ln 2 + 2), which is not equal to the condition of the force continuity, Eq. 101. Thus, if we use Eq. 101, an alternative treatment for the energy continuity is required. For this, going back to Eq. 13b (see also Eq. (17a)), we have using IPS ii (0) = 0 for the Coulomb function (owing to ψ (1) = −γ ). Thus, applying (72) in the scheme described in Appendix A, the atomic energy is and the last term may be approximated as via Eq. 73 or approximated just as zero via Eq. 74. The total energy of the IPS method is described by Eq. 19, i.e., where the last term may be approximated similarly above. In particular, subtracting the boundary term,