For a unit cell of N point charges: q₁, q₂, … , q_N located at r₁, r₂, … , r_N. The Coulomb energy of this unit cell in the straightforward cutoff method is written as

where r_cut is the cutoff distance, the Heaviside step function is θ (x) = 0 for x < 0 and θ (x) = 1 for x ≥ 0, and L(i) is the intramoleclar neighbor list of the i-th charge. L(i) usually includes the i-th charge itself and any other charge that has already interacted with the i-th charge via the intramolecular interaction. The interaction between the i-th charge and any of its intramolecular neighbors is therefore excluded in Eq.  (1). The distance between the two charges r_ij ≡ | r_i − r_j| is usually taken as the minimum distance among their images subject to the periodic boundary condition (PBC). This direct consequence of the PBC is also called the minimum image convention.^[42] In the above equation and the following equations, we have omitted the prefactor 1/(4π ɛ ₀) for simplicity.

Alternatively, one may introduce the group-based cutoff (GB-cutoff) method such that any charge always interacts with particles that are overall neutral. Suppose there are M neutral groups in the unit cell, each of which includes N_m charges for m = 1, 2, … , M, we denote the charge and the coordinate of the j-th charge in the m-th group as q_j(m) and r_j(m) respectively. The electroneutrality condition for each group is then written as . Using this notation, the Coulomb energy in the group-based truncation method can be written as

where is the group distance which can be defined as either the distance between any specified pair of charges or the distance between the two center of mass positions. Note that it is in general nontrivial to divide the total charges of the system into M non-overlapping neutral groups for an arbitrary system of interest. Although it is straightforward to divide the N charges of a system of small molecules (e.g., pure water or a water-acetonitrile mixture) into M neutral groups (e.g., each molecule is a neutral group for water-acetonitrile mixture), it is difficult to do such a similar division for systems of large molecules and ions. Definitions that differ from the above Eq.  (2) might be used and exclusion of the Coulomb interaction from the intramolecular list might have to be considered too as in Eq.  (1). Nonetheless, the truncation type methods are very efficient and its implementation is simple.

One disadvantage of the above truncation methods is that the potential is not continuous at the boundary of the cutoff distance, which can be seen from the use of the Heaviside function in the Coulomb energy. For the cutoff approach, it has been intuitively believed that increasing the cutoff length and using the group-based cutoff strategy might make the simulation result more accurate and reliable. But in 2005 Yonetani^[43] reported that by setting the cut-off length as large as 1.8  nm, a bulk water simulation of TIP3P model water surprisingly produces a spurious layer structure. These results have indicated that for the uniform and neutral system, the truncation method unexpectedly generates significant artifacts. Soon after that, van der Spoel and van Maaren^[44] simulated a series of water systems and calculated the electrostatic interaction via different methods to give a reasonable explanation for the spurious structure found by Yonetani.^[43] These works concluded that only the appropriate simulation parameters can avoid the artificial phase transition. When the layer structure is produced by the use of the truncation method, the density increases, the potential energy reduces and water diffuses anisotropically. From the dipole– dipole correlation, they found that there are a large number of water molecules oriented similarly around a given one, and after some width of such a cluster, a larger number of water molecules orient oppositely. The underlying reason is that each water molecule and its surrounding molecules form a cluster in the vacuum followed by the decrease of the free energy through minimizing its net polarization. Simultaneously the periodic boundary condition strengthens the formation of the artificial layer structure.

Although the truncation methods are efficient and easy to implement, the associated defects seem inevitable at least for some specific systems relevant to biological simulations. The above interesting finding of the artifact of the truncation method is definitely not intuitive. Therefore, it indicates that properties of specific systems could potentially be very sensitive to the accurate treatment of the electrostatics. As water is the most important system relevant to biological functionality and activity, the above finding essentially suggests that very careful checking needs to be done when the truncation methods are used for simulations of biological systems.

The basic idea of the Ewald sum is to divide the Coulomb interaction into the short-ranged part and the long-ranged part using a screening factor α

where erf(x) and erfc(x) ≡ 1 − erf(x) are the error function and the complementary error function respectively. The corresponding short-ranged part (the complementary error function part) of the Coulomb lattice sum is then summed in real space as usual but the long-ranged part (the error function part) is summed in reciprocal space. By defining an infinite boundary term associated with the behavior that the system asymptotically approaches the infinite, it has been recently shown that the Ewald sum for electrostatics in a system that is periodic in one (Ewald1D), two (Ewald2D) or three (Ewald3D) dimensions can be written as a sum of pairwise interactions.^[37] For a system of point charges in bulk, the usual Ewald3D sum for N charges in the rectangular unit cell is thus written as^[37]

where the real space sum, the reciprocal space sum and the infinite boundary term for the system of 3D periodicity are,

and

respectively. The vector n stands for (n_xL_x, n_yL_y, n_zL_z) in which n_x, n_y, and n_z are integers and L_x, L_y, and L_z are the length of the unit cell in the three directions respectively. V = L_xL_yL_z is the volume of the unit cell. The corresponding reciprocal vector k is defined as k ≡ 2π (k_x/L_x, k_y/L_y, k_z/L_z), k = | k| , and k_x, k_y, k_z are all integers. In fact, the pairwise form and the infinite boundary term are not unique for 3D periodicity and the similar definitions can be made for the 3D system with reduced 2D (Ewald2D) or 1D (Ewald1D) periodicity.^[37] The clear definition of the pairwise potential might be important for various understandings of the Ewald sum. For example, the purpose of a biological simulation is often to provide molecular level insight for the effect of a particular selected group of particles on another selected group of particles. However, the presence of the explicit images make it difficult to exactly split the total electrostatic energy into pairs. The usual way of analyzing the real space part of the Ewald sums depends on the choice of the screen factor α . In contrast, the unique definition of v(r) in bulk immediately solves this problem because v(r) does not depend on α any more. Similarly, v^2D(r) and v^1D(r) for systems in reduced periodic dimensions^[37] are useful for the analysis of the trajectory when the corresponding Ewald2D or Ewald1D method are used to simulate interfacial biological systems.

Although the pairwise potential v(r) is conceptually useful, the pairwise form by itself does not directly produce an efficient algorithm. In practice, the widely used particle mesh Ewald (PME) algorithms^{[45– 47]} interpolate the usual reciprocal (Fourier) space term part of the Ewald3D sum on a grid and then utilizes the computationally (or more rigorously Nlog₂N) fast Fourier transform (FFT) technique. This technique is in spirit close to the much earlier particle-particle– particle mesh method.^{[47, 48]}

In the popular smoothed particle mesh Ewald algorithm, ^[46] a powerful B-spline interpolation scheme is used for the structure factor and FFT is implemented for all convolution operations.

The Ewald3D sum consisting of the real space sum and the reciprocal space sum (v^IB = 0) is usually regarded as the standard method to validate other algorithms. One particular truncation method that is closely related to the Ewald method is the Wolf truncation method^[49] in which the real space sum of the Ewald sum is used and the reciprocal space sum is ignored. For biological systems in a crowded environment under the condition of constant volume, the success of the Wolf truncation method or the real space part of the Ewald sum might be understood because the force associated with the long-ranged component in the reciprocal space part essentially cancels each other. However, the virtues and defects of the above discussed truncation methods or any other spherical truncations of Coulomb interaction for biological simulations remains a subject of ongoing debate.^[50] Nonetheless, the real space part of the Ewald sum provides an alternative efficient truncation method that is practically useful at least for bulk simulations of liquid systems.

Although the Ewald sum without the infinite boundary term has been regarded as the standard method in computer simulations of liquid systems in general, a number of artifacts associated with the explicit presence of images were discussed in the literature. We list three types of the potential artifacts in the following. Because the Coulomb lattice sum is a conditionally convergence series, the Ewald sum might in principle have an extra conditional term. This term was called the surface or dipole term^[42] but ambiguity exists because the dipole of the unit cell is not well defined. In most applications of the Ewald-type methods including the PME methods, the conditional term is usually ignored and the resulting Ewald sum with the real and reciprocal sums only is called the Ewald sum with the tin-foil condition boundary corresponding to the physical situation that the macroscopic system is surrounded by a conducting metal.^{[42, 51]} When the boundary condition is written as the pairwise infinite boundary term of Eq.  (7), the alternative formulation in principle removes the ambiguity such that the computation for a given infinite boundary term is not problematic any more. However, it is still unknown which infinite boundary term or which path of k → 0 should be chosen for certain purposes in computer simulations of liquid systems.^[37]

Besides, Smith and Pettitt have reported that there exists an energy barrier when a dipole rotates in the simulation box with its center fixed.^[52] In fact, this broken symmetry of rotational invariance can be clearly seen from the fact that the pairwise potential v(r) depends on the vector r instead of the scalar r. This effect might contribute little to the property of the system because the electrostatic interaction is screened by a large dielectric constant for most bulk liquids (e.g., the dielectric constant of bulk water ɛ ₀ ≃ 80).

For the charged system, the classical Ewald summation should be added by a correction term as a result of adding a uniform neutralizing charged background plasma. Many researchers^{[53– 56]} have pointed out the periodicity-induced artifacts when dealing with the charged system with the Ewald summation method, and given their corrections. However Bogusz et al.^[57] pointed out that such a uniform neutralizing background plasma brings forth serious artifacts to both the system energy and pressure. They also instituted a net-charge correction term to correct these problems. What is more, Hü nenberger and McCammon^[55] performed the detailed investigation on the Ewald artifacts in computer simulations of ionic solvation and ion-ion interactions. Based on continuum electrostatics and two examples both comparing the non-periodic boundary condition and periodic boundary condition, they summarized three factors which may enlarge the periodicity-induced artifacts. No matter what, the finite size effect generated by the periodic boundary condition is widely believed to be corrected by self-energy, which was proposed^[58] and examined on the analysis of free energy of ionic hydration. Ekimoto et al.^[56] recently presented a similar correction scheme for the free energy of charged protein in an explicit solvent.

The above three aspects of the Ewald artifacts are ongoing topics in this field. We believe that the problem of removing the artifacts associated with the Ewald sum has not been settled in principle although the Ewald type method is by far the standard method for the bulk system.

LMF theory maps the structural and thermodynamical properties of a complex system with long-ranged intermolecular interaction in a general external field into those of a simple system with only short-ranged intermolecular interaction but under a reconstructed effective external field. When the basic Coulomb interaction 1/r is separated into a short-ranged component v₀ and a long-ranged component v₁: 1/r = v₀(r) + v₁(r). The LMF theory suggests that the reconstructed external electrostatic potential for a system under the general external electrostatic potential can be written as

where is the (ensemble-averaged) singlet charge density associated with the reconstructed external field. This equation can be solved self-consistently^[32] and the efficient method to produce the solution has been developed.^[33]

When the external field has a certain symmetry in the direction of β (e.g. planar symmetry r^β = z or spherical symmetry r^β = r), one can alternatively introduce the mean field approximation for the instantaneous configuration:^[35]

where 〈 〉 _sp is the symmetry-preserving operator defined as the normalized integration over the degrees of freedom in the directions with no broken symmetry (other than β ). R̄ denotes the collective coordinates of all charges and ρ ^q(r, R̄ ) is the instantaneous charge density. An expansion of using an orthogonal basis always exists

It has been argued that when the long-ranged component is slowly varying, only a small number of terms in the above series are needed to achieve high accuracy.^[35] Therefore, the reconstructed field reaches the virtue of single-particle potential as it is in the LMF theory. Both the LMF and SPMF theories yield algorithms facilitating efficient calculations of electrostatic forces and potentials. SPMF theory differs from LMF in that the reconstructed external field is dynamical and it can be regarded as a special case of dynamical LMF theory when the external field has certain symmetry.

For the system of bulk solutions with translational and rotational invariance symmetry, the symmetry-preserving operator yields a constant and therefore the long-ranged components essentially cancel each other, which is equivalent to the strong coupling approximation in LMF theory^[34] or the Wolf truncation method. Obviously, both the SPMF and LMF theories are more physically based with a clear explanation of the cancellation effect in bulk. This suggests that analysis of the defects of the truncation methods might be achieved in the framework of the LMF theory via considering the failure of the strong coupling approximation.

The Ewald separation suggests the use of v₀(r) = erfc(α r)/r for the short-ranged component and correspondingly v₁(r) = erf(α r)/r. Although the general framework of the LMF theory can be applied to any other separations as well, the Ewald choice is a very convenient choice in practice. Because erf(α r)/r essentially arises from a normalized Gaussian distribution with width 1/α , the reconstructed external field in the LMF theory exactly satisfies Poisson’ s equation but with a Gaussian-smoothed charge density given by the convolution of the point charge distribution with the normalized Gaussian distribution. This property allows many convenient and accurate applications of the LMF and SPMF theories as well. Moreover, corrections to thermodynamics in strong coupling approximation of the LMF theory can be done analytically when the Ewald choice is used.^[61] This previous finding would suggest that perhaps the Ewald separation is the optimal choice for the practical application of the mean-field-type methods to accurately treat electrostatic interactions. However, a solid justification in principle has yet to be developed.