Generalized Moment Correction for Long-Ranged Electrostatics.

Describing long-ranged electrostatics using short-ranged pair potentials is appealing because the computational complexity scales linearly with the number of particles. The foundation of the approach presented here is to mimic the long-ranged medium response by cancelling electric multipoles within a small cutoff sphere. We propose a rigorous and formally exact new method that cancels up to infinitely many multipole moments and is free of operational damping parameters often required in existing theories. Using molecular dynamics simulations of water with and without added salt, we discuss radial distribution functions, Kirkwood-Buff integrals, dielectrics, diffusion coefficients, and angular correlations in relation to existing electrostatic models. We find that the proposed method is an efficient and accurate alternative for handling long-ranged electrostatics as compared to Ewald summation schemes. The methodology and proposed parameterization are applicable also for dipole-dipole interactions.

Introduction

Accurate and efficient schemes to calculate long-ranged electrostatic interactions are important to expand time and space in atomistic simulations. One of the earliest schemes is the reaction field method[1,2] where a spherical cutoff, Rc, is applied. The potential is zero beyond the cutoff, and hence, the computational complexity scales as , where is the number of charged particles in the system. However, to correctly parameterize the reaction field method, the surrounding dielectric constant, εRF, needs to be estimated, and artifacts may arise because of discontinuities at the cutoff.[3] Lattice approaches, including the Ewald method,[4,5] provide an accurate electrostatic description within a well-known parameter space[6,7] but are inherently limited to periodic systems. The computational complexity of Ewald is but reduces to with optimal parameters. This growing cost with increasing system size makes Ewald methods demanding for large systems, albeit derived versions such as Particle Mesh Ewald[8] (PME) has improved scaling.

Alternative cutoff-based methods with scaling have been shown accurate for isotropic systems and include the charge neutralizing Wolf method.[9] This and many derivates thereof have been thoroughly tested in the litterature.[10−13] Evolutions from the original Wolf formalism are primarily based on cancellation of the derivative of the potential at the cutoff, which is convenient because this is a prerequisite for molecular dynamics simulations. Such cancellation of the derivative may be equivalent to canceling an introduced dipolar artifact.[14]

By expanding from the Wolf method, we show here that cancellation of every electrostatic moment—monopole, dipole, quadrupole, octupole, and so forth—is equivalent to taking all long-ranged interactions into account in a conducting medium. We present a potential which in addition to cancelling arbitrarily (including infinitely) many moments P, cancels P − 1 derivatives of the potential at the cutoff. From a physical viewpoint, this is reasonable as to avoid truncation errors, a conjecture that has previously been brought forward.[11] A similar approach of higher order moment cancellation has recently been presented,[15] albeit an operational damping parameter is introduced in addition to the cutoff distance. The approach presented in this work requires (i) a cutoff distance and (ii) choosing how many electric moments, P, to be cancelled. The integer P thus connects directly to physical properties of the examined system, while avoiding the operational damping parameter in similar formalisms. Short-ranged potentials, as the one presented in this work, commonly assume isotropy beyond the given cutoff region and should be cautiously used in systems where this is not a valid premise, for example, at interfaces[16] or in ferroelectric systems. In the next section, we derive and explore the introduced pair potential, whereafter we test it in molecular dynamics simulations of aqueous systems.

Theory

The electrostatic pair-interaction energy u between charges z and z can be written as a function of separation r for all models examined hereHere e is the elementary charge, ε0 is the vacuum permittivity, and εr is the relative permittivity of the dispersing medium. The short-ranged function, , is defined in Table using q = r/Rc and for q > 1, where the latter ensures the potential to be finite-ranged. Applying results in a modified Coulomb potential which still captures the general 1/r distance dependence of the original Coulomb potential. The new method presented here, which we denote as the “q-potential”, cancels an arbitrary number of electrostatic moments and derivatives of the potential at the cutoff and is defined by the short-ranged function

Table 1

Electrostatic Schemes and Their Short-Ranged Functions, . η = αRc Where α is a Damping Parametera

label		refs
Coulomb	1
reaction field		(2)
Ewald, real space	erfc(ηq)	(4)
SP0 (Wolf)	erfc(ηq) – q erfc(η)	(9)
SP1	(1 – q)2	(12)
SP3	(1 + 2.25q + 3q2 + 2.5q3) (1 – q)4	(13)
q-potential	∏n=1P(1–qn)	this work

SP means “shifted potential”, and the suffix denotes the number of potential derivatives that are zero at the cutoff.

In the next subsection, we formally derive the theory.

Generalized Moment Cancellation

The total interaction energy between charged particles using periodic boundary conditions (PBCs) isHere, the prime indicates that i ≠ j when = 0, ◦ denotes the Hadamard product, is the distance vector between point charges i and j, the size of the cuboid cell is described by its side lengths = (L, L, L), and the zeroth-order interaction tensor 0() isBy assuming ||/2 < min(), it is possible to convert all point-charge particles in each replicated cell into a single-centered point-multipole particle represented by infinitely many higher order moments,[17,18] as in Figure . This conversion is valid as the center of any reciprocal point-multipole particle, positioned at ○ , is located further away from the origin of the centered cell (min(| ○ |) ≥ min()) than any point-charge particle in the centered cell itself (max(|r|) = ||/2, where r is the distance vector from point-charge i to the origin). If every moment of every such point-multipole particle is a zero tensor, then there would not be any interactions between reciprocal cells, that is, no long-ranged interactions. For such a system, the interaction energy would simply be

Figure 1

Single point-multipole particle (left) describing a point-charge distribution (right).

Ergo, if an effective potential cancels all the electrostatic moments of the cell, then the total interaction energy would be given solely from the particles in the centered cell. Based on the above description, the procedure can be regarded as a reaction field approach where the induced field cancels the moments locally. Even if elements in higher order tensors do not fully cancel and therefore are not exactly zero, eq might be a highly accurate approximation when the short-ranged interactions (described by higher order tensors) from reciprocal cells only give minor contributions. Consequently, even imperfect conductors with finite P can give accurate results for the long-ranged electrostatics. In fact, for imperfect conductors such as water, an infinite P might cancel too many interactions and thus give erroneous results. For higher order interactions than dipole–dipole and ion–quadrupole interactions, the long-ranged interactions are absolutely convergent. Therefore, the highest order moment that in theory needs to be cancelled, for sufficiently large systems, is the quadrupole moment, that is, P = 3, which we also see in Section . The same arguments, as explained in this section, hold for interactions between any regions with zero or small tensor moments, and the cancellation approach is therefore not limited to PBC systems.

Derivation of the q-Potential

We begin the derivation by assuming that an electrostatic potential based on moment cancellation can be described by the potential from the original particle and P image particles. The aggregated potential from the original and all image particles is denoted as V, which is a function of the distance vector r and the charge z,Here, V represents the Coulomb potential, z is the charge of the original particle, z is the charge of the image particle p, and r = cr, where c is a proportionality factor to be discussed. The requirement that up to M higher order moments of all image particles should cancel the original particle moments can be formulated asHere, r is any component of r and r is the corresponding component of r. The meaning of row m in eq is that the m – 1 order moment from the original particle (first term; rz) together with the sum of the same order moments of the image particles (second term; ∑rz) should equal zero (right-hand side). Assuming that the image particle positions r are unique and non-zero and that the matrix is a square (M = P – 1), we ensure that a solution exists.[19] By then choosing the positions r, we can extract the charges, z. The positions of the image particles may be arbitrarily chosen, with the exceptions mentioned earlier; however, in this work, we use c = q–, that is, the position of image particle p + 1 (r) is the mirror image of the position of image particle p – 1 (r) in the position of image particle p (r), that is, r2 = q–2r2 = q–(rq–(r = rr. This choice is closely related to the method of image charges.[20,21] Together with eq , the above mentioned assumptions yield the solution for the image charges aswhere is the q-analogue of the binomial coefficients.[22,23] For details, see the Supporting Information, Section S1. Using eq , it is further possible to present an expression for the aggregated potential V and the modified interaction tensor asThe expression within the parenthesis is composed of the contributions from the original particle (the 1) and all image particles (the sum). Rearranging and simplifying[23] this expression giveswhere (a;q) is the q-Pochhammer symbol. In Figure , we show this multiplicative modification to the original interaction tensor for different P values. Because the multiplicative term is a q-analogue of the Pochhammer symbol, we address the proposed potential as the q-potential.

Figure 2

Short-ranged function for different numbers of cancelled moments, and (inset) the same scaled with the Jacobian 4πq2 as to assess its volumetric influence.

By rewriting eq asand specifically noting the (1 – q) term, we acknowledge that the number of cancelled higher order moments, P – 1, equals the number of higher order derivatives with respect to || (or equivalently q) that are zero at the cutoff. Finally, using the already described theory, but with image particle positions c = q–, we retrieve the short-ranged functionwhere s ∈ {1, 2, ...} and s – 1 higher order derivatives are cancelled also at q = 0, since the lowest power of q in (q;q) (save q0 = 1) is q. Consequently, the (s – 1)th derivative is the highest order derivative to vanish at q = 0 (given that P > 0). The expression in eq can now be tuned in how many derivatives to cancel at q = 0 (using s) and at q = 1 (using P). Because of the physical interpretation of the image charge positions,[20,21] we, nonetheless in this work, use s = 1.

A pair interaction entails two particles, and above, we cancelled moments of one of them, specifically the one at distance r from the origin (where the other particle is located). The charge of the origin-located particle has thus not been cancelled; however, we now address this issue in a similar manner as the previously explained moment cancellation procedure. The resulting energy is defined as the self-energy, and in Section S2, the detailed derivation is presented, which culminates into

Potential for Systems with Moments

The no net-moment assumption is fair for many liquids, but generally not so at interfaces or in systems with dipolar or ferroelectric properties. For such systems, we suggest a generalization of the presented methodology. We previously noted that cancelling every moment within a cutoff region gives the right-hand side of eq as the zero vector. If, however, there would be a non-zero moment within this region, one gets eq , where Ψ( is the pth higher order moment in the region projected onto the r-component of r. The solution to this equation is more comprehensive than the solution to eq , yet solvable, at least given the mentioned restrictions on r and M. Here, we have not explicitly expanded the given framework for such cases; however, we highlight the possibility to develop a q-potential for non-zero moment environments.

Choosing the Cutoff

It is desirable to use a minimal cutoff to speed up the calculation time. The physical interpretation of such a cutoff which still accurately describes the system is the smallest region which exhibits small enough fluctuations in the aggregated moments to be corrected for by induced image particles. It is formally sound to use a larger cutoff as larger regions maintain this quality. However, most systems do display a certain local anisotropy or equivalently nonvanishing higher order moments. Accordingly, the cutoff needs to enclose this whole space. To ensure no intercell correlations between anisotropic regions, it should not exceed one-fourth of the shortest cell length.[24] Finally, the no net-moment approach does not impose an isotropic cutoff and any closed space shape can be used. The choice of cutoff shape, we speculate, may be particularly important in strongly anisotropic regions such as interfaces.

Methods

Molecular Dynamics Simulation

All simulations were performed in the isobaric–isothermal ensemble at pressure bar and temperature T = 298.15 K, using OpenMM 7.[25] The pressure was kept constant using a Monte Carlo barostat,[26] and we used a Langevin integrator,[27] with a friction coefficient of 1.0 ps and a time step of 2.0 fs, to maintain a constant temperature. Before production runs, the systems were energy-minimized and then equilibrated. The reference Ewald summation/PME simulations used a fractional force error tolerance of 5 × 10–4.

For simulations of only SPC/E water molecules, Ewald summation was used as a reference, and the production runs spanned 10 ns. Three different atom-based (real space) cutoff distances were used for all electrostatic schemes: 0.96, 1.28, and 1.60 nm, which roughly correspond to 3, 4, and 5 oxygen (∼water) diameters. The Lennard–Jones interactions used the same cutoff and a long-ranged correction term to compensate for the neglected contributions.[28] Systems containing SPC/E water molecules and NaCl/NaI ions were simulated to retrieve activity derivatives according to Kirkwood–Buff (KB) theory, as described in the next subsection. The simulations were carried out using Rc = 1.28 nm in a cubic box with water molecules depending on the type and concentration of the electrolyte presented in Appendix A. The production runs spanned 40 ns, and PME was used as a reference.

The diffusion coefficient defined in eq was calculated by least-square fitting of the mean-square displacement (MSD) of the particles throughout the simulation to a linear curve, where the slope is the sought after value. The MSD was sampled in the microcanonical ensemble, which continued upon the last state of the above-explained isobaric–isothermal simulations with a 1 ns equilibration and 2 ns production run. The microcanonical production simulation also provided the base for the analysis of the energy.

The standard deviations were based on samples from ten consecutive, equally sized intervals that spanned the whole simulation.

KB Theory

The KB theory provides a general way to obtain macroscopic properties of a solution from its microscopic properties.[29] The central property is the KB integral (KBI) between components A and B, defined aswhere gμ(r) is the radial distribution function (RDF) in the grand canonical ensemble (i.e,. an open system) and r is the distance between the components. To obtain approximate KBIs in a closed system, the integral is typically truncated at a distance R after which the RDFs converge, and the grand-canonical RDF is replaced by the RDF computed in the closed system.[30,31] To obtain accurate approximate KBIs according to this procedure, corrections are needed to compensate for the introduced errors. We followed a procedure[32] in which two corrections[33,34] are applied simultaneously to the KBIs. These correction factors are further explained in Appendix A. The derivative of the electrolyte activity at constant pressure and temperature is obtained according to[35]Here, Gcc and Gwc are the corresponding truncated, corrected KBIs obtained from simulations, Ĝcc*(R) and Ĝwc*(R), respectively (see Appendix A). Subscripts c and w denote the electrolyte and water, respectively, ac = γcρc is the electrolyte activity, γc is the molar mean activity coefficient of the electrolyte, and ρc is the number density of the electrolyte. For the water electrolyte KBI, we chose to use the RDF between the water oxygen and the electrolyte, in the following denoted as gOc(r). Experimental activity derivatives for the simulated electrolytes were calculated from previously reported activity coefficients,[36] using the fitting procedure described in Appendix A. For the simulations, we used one force field[37] for the Na+ and I– ions, while another[38] was used for Cl–. The parameters of these force fields are presented in Appendix A.

Results

To evaluate the developed q-potential, we investigate a bulk water system and an aqueous salt solution by analyzing, among others, radial distributions functions, angular correlations, and KBIs. Although the Ewald/PME methods assume a replicated environment and may therefore not necessarily represent a true isotropic system,[24,39] we choose this as a reference system because of its widespread use in molecular simulations.

Bulk Water System

A representation of the radial distribution functions between water oxygen atoms is presented in Figure , where the q(P = 1) potential stands out with a clear peak at the respective cutoff distance. By initially increasing the order P from 2, we get results closer resembling those of Ewald. However, using q(P = ∞) and Rc ≤ 1.28 nm, we diverge from the Ewald-like results retrieved by intermediate P-values. The dipole-dipole correlation for water depicted in Figure shows a similar trend.

Figure 3

Deviation from the bulk presented as gOO(r) – 1 scaled by the Jacobian 4πr2, where gOO(r) is the oxygen–oxygen RDF. For clarity, we have faded the strongly oscillating q(P = 1) lines. The SP1 approach cancels one derivative at the cutoff and yields results more closely resembling Ewald than the similar q(P = 2) potential. However, the q(P = 3) potential yields results even closer to those of Ewald. The SP3 approach, which cancels three derivatives at the cutoff, gives the most Ewald-like results akin the q(P = 5) and q(P = 7) potentials.

Figure 4

Dipole–dipole correlation as a function of separation using normalized dipole μ̂ and (inset) the difference to Ewald where the ordinate spans ±5 × 10–3.

The density, dielectric constant (see Section S3 in the Supporting Information), Kirkwood factor GK (see eq ), standard deviation of the total energy, and diffusion coefficient for the different potentials are presented in Table . Again, we note that q(P = 1) gives distinctly different results than the others and will therefore not be commented from here on. The densities for all other q-potentials are generally slightly above a thousand. However, with a standard deviation of ∼2, nearly all potentials (for Rc ≥ 1.28 nm) are statistically indistinguishable. Although the dielectric constant is seen to be irregular in P, the results using q(P > 1), SP1, and SP3 are in reasonable agreement with Ewald and experimental values. For the Kirkwood factor, all potentials with P ≠ 1 are consistent within the range 3.0 ± 0.1, except q(P = 6), for which Rc = 0.96 nm gives GK = 3.2. When comparing the standard deviation of the energies of the pair potentials to the, for PBC formally exact, Ewald reference, we observe low values for all approaches. Especially, SP1 and SP3 yield particularly low values, even when compared to Ewald. However, a seemingly odd behavior is the increase in σE for P = 4 and P = 5 while expanding the cutoff from Rc = 1.28 nm to Rc = 1.60 nm. A similar trend is also observed for the Ewald summation result. The diffusion coefficients are all similar to results of previous studies[42,43] with D = 2.5 ± 0.3 for all pair potentials with Rc ≥ 1.28 nm.

Table 2

Density ρ [kg/m3], Relative Dielectric Constant εr, Kirkwood Factor GK [Unitless], Standard Deviation of Total Energy σE [kJ/mol], and Diffusion Coefficient D [m2/s/10–9] for the Different Potentials Applied on a Bulk Water System and Experimental Refs (40) and (41)a

	Rc = 0.96 nm					Rc = 1.28 nm					Rc = 1.60 nm
potential	ρ	εr	GK	σE	D	ρ	εr	GK	σE	D	ρ	εr	GK	σE	D
q(P = 1)	1137		4.8			1106		3.2			1073		3.0
q(P = 2)	1005	78	3.0	24	2.2	1002	75	3.0	15	2.2	1000	69	2.9	11	2.7
q(P = 3)	1002	67	2.9	10	2.3	1001	70	3.0	8	2.5	1000	76	3.0	3	2.4
q(P = 4)	1002	67	2.9	9	2.4	1000	67	2.9	6	2.5	1000	69	2.9	10	2.4
q(P = 5)	1003	71	3.0	8	2.2	1001	67	2.9	8	2.6	1000	72	2.9	12	2.3
q(P = 6)	1004	76	3.2	5	2.4	1001	71	3.0	4	2.4	1000	68	2.9	3	2.7
q(P = 7)	1005	73	3.0	8	2.5	1001	68	2.9	4	2.4	1000	69	2.9	3	2.5
q(P = 8)	1006	72	3.0	9	2.5	1001	70	2.9	6	2.6	1000	71	3.0	6	2.5
q(P = ∞)	1008	73	3.1	9	2.4	1001	66	2.9	4	2.7	1000	67	2.9	4	2.5
SP1	993	67	2.9	1	2.9	996	69	2.9	1	2.8	996	68	2.9	1	2.8
SP3	998	73	3.0	0	2.4	998	66	2.9	0	2.5	998	71	3.0	0	2.5
Ewald	998	67	2.9	4	2.5	998	73	3.0	2	2.6	998	69	3.0	6	2.9
Exp	997	79			2.3	997	79			2.3	997	79			2.3

The standard deviations for the density, dielectric constant, and diffusion coefficient are σρ ≈ 2 kg/m3, σε ≈ 1 (over last half of simulation), and σD ≈ 0.1, respectively. P = 1 is an exception and gives higher values.

Finally, as stated in Section , we note that the cutoff should obey Rc < min(L, L, L)/4, and therefore, because the box length is ∼3.9 nm, artifacts may be present in systems using excessive cutoffs. However, in Section S4, we present simulation results based on larger systems and see no effects not previously pointed out.

Salt Solutions

Figure shows activity derivatives for NaCl and NaI solutions, and the q(P = 3) potential generally gives the most accurate results as compared to PME. However, the presented error bars for all methods do overlap, which indicates comparable results. We have performed simulations using the q(P = 1) potential with a poor outcome (akin for the RDF), wherefore we conclude that cancellation larger than one is needed to capture the nature of the system. For m = 3 mol kg–1, all potentials—including PME—yield systematically lower activity derivatives compared to experimental data.

Figure 5

Activity derivatives for NaCl (left) and NaI (right) using different pair potentials, with Rc = 1.28 nm. Solid curves represent experimental data.[36,44] For visual clarity, the PME results are plotted on top of the experimental curves (in black), whereas the results from the pair potentials are shifted in steps of 0.3 and plotted on top of the correspondingly shifted experimental curves.

Figure shows RDFs between water and salt for the different potentials. Similarly to the oxygen–oxygen RDFs (Figure ), the lowest order of P, here P = 2, yields results deviating the most compared to PME. Increasing the order up to P = 7 improves the agreement, whereas P = ∞ again yields results showing increased deviation from PME.

Figure 6

Water oxygen–electrolyte correlation plotted as gOc(r) – 1 scaled with the Jacobian 4πr2, where gOc(r) is the water oxygen–electrolyte RDF, for different pair potentials including PME. The results are for NaCl solutions with a salt concentration of 2.0 mol kg–1, with Rc = 1.28 nm. The inset shows a magnification to emphasize the discrepancies between pair potentials at long distances.

Discussion

For the simulated water system using Rc = 1.28 nm, we have found P = {4, 5, 6} to closely match Ewald and PME results. This agrees with results for the dipolar q-potential in Stockmayer systems using a similar cutoff.[45] There is thus a limit on how many cancellations are appropriate. Although the discrepancies between the pair potentials and Ewald are oscillating, the peaks consistently occur around the cutoff distance. Arguably, the dipole–dipole correlation as a function of separation (Figure ) is more sensitive to electrostatic interactions than the radial distribution function (Figure ). For example, in previous studies using the Wolf formalism, it has been shown that suboptimal combinations of cutoff and damping parameters can lead to unphysical results, where water dipoles align toward each other.[3,46] We do not see such tendencies with the q-potential, which for all conditions show the same trends as Ewald. The results of the dipole–dipole correlation difference to Ewald (inset of Figure ) at Rc = 1.28 nm display distinguishable oscillating trends. However, the amplitudes are of the same magnitude as the noisy segments at other distances, in contrast with the results using Rc = 0.96 nm. Ergo, this suggests Rc = 1.28 nm to be the lowest cutoff to give Ewald-like results on every distance. Furthermore, the discrepancies in Figures and 4 seem to imply deviations from Ewald in densities, which is most clearly seen for the results using Rc = 0.96 nm.

From a theoretical perspective, we expect all long-ranged interactions to vanish using P = ∞, that is, cancelling every moment locally. From the point-of-view of long-ranged interactions, this would solve the problem of summation if indeed the surrounding medium did induce perfect moment cancellation. In this study, we observe that even water with its rather high relative permittivity does not act as a perfect conductor. An infinite P suppresses all local fluctuations and gives misrepresented short-ranged interactions. Thus, choosing too large P results in overdamping of the system compared to the reference, which suggests that moment cancellation gives physical results only up to a point. This is nonetheless only the case for the two smaller cutoffs, while for the largest one, there is no upper limit for the optimal P. As Rc → ∞, the q-potential collapses to the unmodified exact Coulomb potential which is independent of P. Conversely, for Rc → 0, no P-values give a proper description of the electrostatics, and the optimal interval for P therefore varies with the cutoff. Optimal values are about P = 5 for Rc = 0.96 nm, P ∈ {4, 5, 6} for Rc = 1.28 nm, and P ∈ {3, 4, ...} for Rc = 1.6 nm. However, note that even the most precise results using Rc = 0.96 nm are still inferior to a mediocre parameterization at Rc = 1.6 nm. The optimal interval seemingly widens with increasing cutoff, and as discussed in Section , we find the expected lower value P = 3 for a sufficiently large cutoff.

The reason why all pair potentials reproduce PME activity derivatives may either be similarities in their RDFs or that the differences cancel while integrating to get the KBIs. For the q-potential, even the lowest order depicted, P = 2, with RDFs deviating significantly more from PME relative to higher orders yields accurate activity derivatives. This implies that the oscillating differences observed between RDFs using PME and the other pair potentials cancel when integrated and do not necessarily pose an obstacle for retrieving accurate activity derivatives.

As for the computational complexity, the Ewald summation real space pair potential is proportional to , where the prefactor is implementation dependent. Evaluation of the specifically used short-ranged function is part of this prefactor. The complementary error function, erfc, used in Ewald is commonly approximated[47] using seven multiplications and one exp call where the latter requires at least ten times more clock cycles than a single multiplication. The corresponding exact evaluation of the q-potential short-ranged function (q;q) includes 2(P – 1) multiplications, which for any relevant choice of P ≲ 5 is significantly faster than the real space part of Ewald and PME. In addition to the real space evaluation, Ewald and PME include a reciprocal space. The optimal real space cutoff for Ewald summation in the simulated bulk system can roughly be evaluated[6] to Rc ≈ 1.6 nm. We find that a lower cutoff can be used for the q-potential. Therefore, for all valid P values and Rc ≤ 1.6 nm, we find that the q-potential is computationally less costly than both Ewald summation and PME and scales with .

Conclusions

We have presented a new theory for summation of long-ranged electrostatic interactions which relies on electrostatic moment and derivative cancellation and has computational complexity of . The new method is a generalization of the Wolf method[9] with the advantages that it (i) allows for cancellation of any number of moments using an integer; (ii) is free of operational damping parameters needed in similar algorithms; (iii) is mathematically rigorous and has a simple form; and (iv) reproduces results from Ewald and PME summation techniques for water and electrolyte solutions.

Results for solution structures, dielectric properties, diffusion coefficients, and activity derivatives of aqueous electrolyte systems agree well with Ewald and PME for a range of P values, that is, the number of cancelled moments. We find that Rc = 1.28 nm (∼ 4 molecular diameters) and P = {4, 5, 6} give appropriate results, where the interval size widens as the cutoff increases. The q-potential method is therefore a valid alternative to Ewald and PME at a lower computational cost. The methodology has been expanded to dipolar interactions and is applicable also for explicit polarization simulations.

Table 3

Fitting Parameters Used to Plot Experimental Activity Coefficients as a Function of Salt Molality Using Eq

salt	b1	b2	b3	b4
NaCl	1.4369	0.0054	0.0495	0.0092
NaI	1.4681	0.1361	0.0344	–0.0102

Table 4

Lennard–Jones Parameters for Na+, Cl+, and I–: Charge z, Diameter d, and Interaction Strength ϵ.a

salt	ion	z	d/nm	ϵ/kJ mol–1	1.0	2.0	3.0
NaCl	Na+	+1.0	0.255	0.280	5963	5252	4672
	Cl–	–1.0	0.440	0.418
NaI	Na+	+1.0	0.255	0.280	5796	5026	4402
	I–	–1.0	0.491	0.158

To the right, the number of water molecules in each simulation is presented for the different molal salt concentrations (mol kg–1).