Generalized Moment Correction for Long-Ranged Electrostatics

TL;DR

This paper introduces a rigorous, parameter-free method for accurately approximating long-range electrostatics in molecular simulations by canceling multipole moments, offering an efficient alternative to traditional Ewald summation.

Contribution

A new formally exact method that cancels infinitely many multipole moments without damping parameters, improving long-range electrostatics modeling in molecular dynamics.

Findings

Method accurately reproduces radial distribution functions.

Achieves comparable accuracy to Ewald summation.

Applicable to dipole-dipole interactions.

Abstract

Describing long-ranged electrostatics using short-ranged pair potentials is appealing since the computational complexity scales linearly with the number of particles. The foundation of this approach 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 is applicable also for dipole-dipole interactions.

Figures

Tables

Table 1. Table 1 : Electrostatic schemes and their short-ranged functions, 𝒮 ​ ( q ) 𝒮 𝑞 \mathcal{S}(q) . η = α ​ R c 𝜂 𝛼 subscript 𝑅 𝑐 \eta=\alpha R_{c} where α 𝛼 \alpha is a damping-parameter. SP means “shifted potential” and the suffix denotes the number of potential derivatives that are zero at the cutoff.

label	$𝒮\left(q\right)$	ref.
Coulomb	1
Reaction-field	$1+q^3\left(\frac{{\epsilon }_{RF}-1}{2{\epsilon }_{RF}+1}\right)$	2
Ewald, real space	$\mathrm{erfc}\left(\eta q\right)$	4
SP0 (Wolf)	$\mathrm{erfc}\left(\eta q\right)-q\mathrm{erfc}\left(\eta \right)$	9
SP1	${\left(1-q\right)}^2$	12
SP3	$\left(1+2.25q+3q^2+2.5q^3\right){\left(1-q\right)}^4$	13
$q$-potential	${\prod }_{n=1}^P\left(1-q^n\right)$	this work

Table 2. Table 2 : Density ρ 𝜌 \rho [kg/m 3 ], relative dielectric constant ε r subscript 𝜀 𝑟 \varepsilon_{r} [unitless], Kirkwood factor G K subscript 𝐺 𝐾 G_{K} [unitless], standard deviation of total energy σ E subscript 𝜎 𝐸 \sigma_{E} [kJ/mol], and diffusion coefficient D 𝐷 D [m 2 /s/10 -9 ], for the different potentials applied on a bulk water-system and experimental reference 40 , 41 . The standard deviations for the density, dielectric constant, and diffusion coefficient are: σ ρ ∼ 2 similar-to subscript 𝜎 𝜌 2 \sigma_{\rho}\sim 2 kg/m 3 , σ ε r ∼ 1 similar-to subscript 𝜎 subscript 𝜀 𝑟 1 \sigma_{\varepsilon_{r}}\sim 1 (over last half of simulation), and σ D ∼ 0.1 similar-to subscript 𝜎 𝐷 0.1 \sigma_{D}\sim 0.1 . The exceptions are all for P = 1 𝑃 1 P=1 which gives higher values.

Potential	$R_c=0.96$ nm					$R_c=1.28$ nm					$R_c=1.60$ nm
	$\rho$	${\epsilon }_r$	$G_K$	${\sigma }_E$	$D$	$\rho$	${\epsilon }_r$	$G_K$	${\sigma }_E$	$D$	$\rho$	${\epsilon }_r$	$G_K$	${\sigma }_E$	$D$
$q\left(P=1\right)$	1137	—	4.8	—	—	1106	—	3.2	—	—	1073	—	3.0	—	—
$q\left(P=2\right)$	1005	78	3.0	24	2.2	1002	75	3.0	15	2.2	1000	69	2.9	11	2.7
$q\left(P=3\right)$	1002	67	2.9	10	2.3	1001	70	3.0	8	2.5	1000	76	3.0	3	2.4
$q\left(P=4\right)$	1002	67	2.9	9	2.4	1000	67	2.9	6	2.5	1000	69	2.9	10	2.4
$q\left(P=5\right)$	1003	71	3.0	8	2.2	1001	67	2.9	8	2.6	1000	72	2.9	12	2.3
$q\left(P=6\right)$	1004	76	3.2	5	2.4	1001	71	3.0	4	2.4	1000	68	2.9	3	2.7
$q\left(P=7\right)$	1005	73	3.0	8	2.5	1001	68	2.9	4	2.4	1000	69	2.9	3	2.5
$q\left(P=8\right)$	1006	72	3.0	9	2.5	1001	70	2.9	6	2.6	1000	71	3.0	6	2.5
$q\left(P=\infty \right)$	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

Table 3. Table 3 : Fitting parameters used to obtain activity coefficients as a function of salt molality using Equation 24 .

Salt	$b_1$	$b_2$	$b_3$	$b_4$
NaCl	1.4369	0.0054	0.0495	0.0092
NaI	1.4681	0.1361	0.0344	-0.0102

Table 4. Table 4 : Lennard-Jones parameters for Na + , Cl + , and I - : charge z 𝑧 z , diameter d 𝑑 d , and interaction-strength ϵ italic-ϵ \epsilon . To the right the number of water-molecules in each simulation is presented for the different molal salt concentrations (mol -1 kg -1 ).

Salt	Ion	$z$	$d$/nm	$ϵ$/kJ$\cdot$mol-1		${𝒩}_{H_{2}O}$
					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

Equations

$$u_{ij}=\frac{e^2{z}_i{z}_j}{4\pi {\epsilon }_0{\epsilon }_{r}r}\mathcal{S}(q).$$

$$\mathcal{S}(q)=\prod _{n=1}^P(1-q^n).$$

$$E_{\mathrm{Tot}}=\frac{e^2}{4\pi {\epsilon }_0{\epsilon }_r}\frac{1}{2}\sum _{\mathbf{n}\in {\mathbb{Z}}^3}^{\prime }\sum _{i=1}^{\mathcal{N}}\sum _{j=1}^{\mathcal{N}}{z}_i{\mathbf{T}}_0\left({\mathbf{r}}_{ij}+\mathbf{n}\circ \mathbf{L}\right)z_j.$$

$${\mathbf{T}}_0(\mathbf{r})=\frac{1}{|\mathbf{r}|}.$$

$$E_{\mathrm{Tot}}=\frac{e^2}{4\pi {\epsilon }_0{\epsilon }_r}\frac{1}{2}\sum _{i=1}^{\mathcal{N}}\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^{\mathcal{N}}{z}_i{\mathbf{T}}_0({\mathbf{r}}_{ij})z_j.$$

$$V_q(\mathbf{r},z)=V(\mathbf{r},z)+\sum _{p=1}^{P}V({\mathbf{r}}_p,z_p).$$

$$\left[\begin{array}{c}1\\ r\\ ⋮\\ r^M\end{array}\right]z+\left[\begin{array}{cccc}1& 1& \cdots & 1\\ r_1& r_2& \cdots & r_P\\ ⋮& ⋮& \ddots & ⋮\\ r_1^M& r_2^M& \cdots & r_P^M\end{array}\right]\left[\begin{array}{c}{z}_1\\ z_2\\ ⋮\\ z_P\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ ⋮\\ 0\end{array}\right].$$

$$z_p=z\cdot {\left[\genfrac{}{}{0ex}{}{P}{p}\right]}_q(-1{)}^p{q}^{p(p-1)/2}$$

$${\mathbf{T}}_0^q(\mathbf{r})={\mathbf{T}}_0(\mathbf{r})\left(1+\sum _{p=1}^P\frac{{\left[\genfrac{}{}{0ex}{}{P}{p}\right]}_q(-1{)}^p{q}^{p(p-1)/2}}{q^{-p}}\right).$$

$${\mathbf{T}}_0^q(\mathbf{r})={\mathbf{T}}_0(\mathbf{r})\prod _{n=1}^P(1-q^n)={\mathbf{T}}_0(\mathbf{r})(q;q{)}_P$$

$${\mathbf{T}}_0^q(\mathbf{r})={\mathbf{T}}_0(\mathbf{r})(1-q{)}^P\prod _{n=1}^P\sum _{k=0}^{n-1}{q}^k$$

$${\mathbf{T}}_0^q(\mathbf{r})={\mathbf{T}}_0(\mathbf{r})(q^s;q^s{)}_P$$

$$E_{\mathrm{Self}}=-\frac{e^2}{4\pi {\epsilon }_0{\epsilon }_r{R}_c}\sum _{i=1}^{\mathcal{N}}{z}_i^2.$$

$$\left[\begin{array}{c}1\\ r\\ ⋮\\ r^M\end{array}\right]z+\left[\begin{array}{cccc}1& 1& \cdots & 1\\ r_1& r_2& \cdots & r_P\\ ⋮& ⋮& \ddots & ⋮\\ r_1^M& r_2^M& \cdots & r_P^M\end{array}\right]\left[\begin{array}{c}{z}_1\\ z_2\\ ⋮\\ z_P\end{array}\right]=\left[\begin{array}{c}{\Psi }^{(0)}\\ {\Psi }^{(1)}\\ ⋮\\ {\Psi }^{(M)}\end{array}\right]$$

$$D=\underset{t\to \infty }{\lim }\left(\frac{1}{6}\frac{\text{d}⟨|\mathbf{r}(t)-\mathbf{r}(0){|}^2⟩}{\text{d}t}\right)$$

$$G_{AB}=4\pi {\int }_0^{\infty }[g_{AB}^{\mu VT}(r)-1]r^2\text{d}r$$

$$a_c^{\prime }={\frac{\partial \ln a_c}{\partial \ln {\rho }_c}|}_{\mathcal{P}T}=\frac{1}{1+{\rho }_c(G_{cc}-G_{wc})}.$$

$$G_K=⟨\sum _{i=1}^{\mathcal{N}}\sum _{j=1}^{\mathcal{N}}\cos ({\hat{\mu }}_i\cdot {\hat{\mu }}_j)⟩$$

$${\hat{G}}_{AB}(R)=4\pi {\int }_0^R[g_{AB}^{\mathcal{NP}T}(r)-1]r^2\text{d}r.$$

$$C_{AB}^{(1)}(r)=\frac{N_B\left(1-\frac{V(r)}{V_{cell}}\right)}{N_B\left(1-\frac{V(r)}{V_{cell}}\right)-\Delta N_{AB}(r)-{\delta }_{AB}}.$$

$$G_{AB}=\frac{1}{V}{\int }_V{\int }_V[g_{AB}^{\mu VT}(r_{12})-1]\text{d}{\mathbf{r}}_1\text{d}{\mathbf{r}}_2$$

$$C_{AB}^{(2)}(r)=1-{\left(\frac{r}{R}\right)}^3.$$

$$G_{AB}\approx {\hat{G}}_{AB}^{\ast }(R)=4\pi {\int }_0^R[g_{AB}^{\mathcal{N}pT}(r)C_{AB}^{(1)}(r)-1]r^2{C}_{AB}^{(2)}(r)\text{d}r$$

$$\ln \gamma (m)=-\frac{1.18\sqrt{m}}{1+b_1\sqrt{m}}-\ln (1-b_{2}m)+b_{3}m+b_4{m}^2.$$

Full text

Lund University] Department of Chemistry, Division of Physical Chemistry, Lund University, Sweden \alsoaffiliation[Lund University] Department of Chemistry, Division of Theoretical Chemistry, Lund University, Sweden

Lund University] Department of Chemistry, Division of Theoretical Chemistry, Lund University, Sweden

Lund University] Department of Chemistry, Division of Theoretical Chemistry, Lund University, Sweden \alsoaffiliation[LINXS] LINXS – Lund Institute of Advanced Neutron and X-ray Science, Scheelevägen 19, SE-223 70 Lund, Sweden

Table of Contents Graphic

1 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, $R_{c}$, is applied. The potential is zero beyond the cutoff, and hence the computational complexity scales as $\mathcal{O}(\mathcal{N})$ where $\mathcal{N}$ is the number of charged particles in the system. However, to correctly parameterize the reaction field method, the surrounding dielectric constant, $\varepsilon_{RF}$, needs to be estimated and artifacts may arise due to discontinuities at the cutoff 3.

Lattice approaches, including the Ewald method4, 5, provides an accurate electrostatic description within a well-known parameter space6, 7 but are inherently limited to periodic systems. The computational complexity of Ewald is $\mathcal{O}(\mathcal{N}^{2})$ but reduces to $\mathcal{O}(\mathcal{N}^{3/2})$ 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 $\mathcal{O}(\mathcal{N}\log\mathcal{N})$ scaling.

Alternative cutoff-based methods with $\mathcal{O}(\mathcal{N})$ scaling have been shown accurate for isotropic systems and include the charge neutralizing Wolf method9. This and many derivates thereof have been thoroughly tested in the litterature 10, 11, 12, 13. Evolutions from the original Wolf formalism are primarily based on cancellation of the derivative of the potential at the cutoff, which is convenient since this is a prerequisite for Molecular Dynamics simulations. Such cancellation of the derivative may be equivalent to canceling an introduced dipolar artefact14.

By expanding from the Wolf method, we here show that cancellation of every electrostatic moment – monopole, dipole, quadrupole, octupole, etc. – 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, cancels equally many higher derivatives of the potential at the cutoff. From a physical view-point this is reasonable as to to avoid truncation errors; a conjecture that has previously been brought forward11. A similar approach of higher order moment-cancellation has recently been presented15, 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 cancel. 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 interfaces16 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.

2 Theory

The electrostatic pair-interaction energy $u_{ij}$ between charges $z_{i}$ and $z_{j}$ can for all here examined models be written as a function of separation $r$,

[TABLE]

Here $e$ is the elementary charge, $\varepsilon_{0}$ is the vacuum permittivity, and $\varepsilon_{r}$ the relative permittivity of the dispersing medium. The short-ranged function, $\mathcal{S}(q)$, is defined in Table 1 using $q=r/R_{c}$ and $\mathcal{S}(q)\equiv 0$ for $q>1$ where the latter ensures the potential to be finite-ranged. Applying $\mathcal{S}(q)\neq 1$ 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 the “$q$-potential”, cancel an arbitrary number of electrostatic moments and derivatives of the potential at the cutoff, and is defined by the short-ranged function

[TABLE]

In the next section we formally derive the theory.

2.1 Generalized Moment Cancellation

The total interaction energy between $\mathcal{N}$ charged particles using periodic boundary conditions (PBC) is

[TABLE]

Here the prime indicates that $i\neq j$ when $\bm{n}=\bm{0}$, $\circ$ denotes the Hadamard product, $\bm{r}_{ij}$ is the distance-vector between point-charges $i$ and $j$, the size of the cuboid cell is described by its side-lengths $\bm{L}=(L_{x},L_{y},L_{z})$, and the zeroth order interaction-tensor $\bm{T}_{0}(\bm{r})$ is

[TABLE]

By assuming $|\bm{L}|/2<\min(\bm{L})$ 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 moments17, 18 as in Fig. 1. This conversion is valid as the center of any reciprocal point-multipole particle, positioned at $\bm{n}\circ\bm{L}$, is located further away from the origin of the centered cell ($\min_{\bm{n}}(|\bm{n}\circ\bm{L}|)\geq\min(\bm{L})$) than any point-charge particle in the centered cell itself ($\max(|{\bf r}_{i}|)=|\bm{L}|/2$ where ${\bf r}_{i}$ 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, i.e. no long-ranged interactions. For such a system the interaction-energy would simply be

[TABLE]

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 box. Based on the above description the procedure can be regarded as a reaction field approach where the induced field cancels moment locally. Even if elements in higher order tensors does not fully cancel and therefore are not exactly zero, Eq. 5 might be a highly accurate approximation when the short-ranged interactions (described by higher order tensors) from reciprocal cells only gives minor contributions. Consequently, even imperfect conductors with finite $P$ can give accurate results for the long-ranged electrostatics. In fact, for imperfect conductors like 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, i.e. $P=3$ which we also see in Section 4. 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.

2.2 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 $V_{q}$, which is a function of the distance-vector ${\bf r}$ and the charge $z$,

[TABLE]

Here $V$ represents the Coulomb potential, $z$ the charge of the original particle, $z_{p}$ the charge of image particle $p$, and ${\bf r}_{p}=c_{p}{\bf r}$ where $c_{p}$ 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 as

[TABLE]

Here $r$ is any component of ${\bf r}$, and $r_{p}$ is the corresponding component of ${\bf r}_{p}$. The meaning of row $m$ in Eq. 7 is that the $m-1$ order moment from the original particle (first term; $r^{m-1}z$) together with the sum of the same order moments of the image particles (second term; $\sum_{p=1}^{P}r_{p}^{m-1}z_{p}$) should equal zero (right-hand-side). Assuming that the image particle positions $r_{p}$ are unique and non-zero, and that the matrix is square ($M=P-1$), we ensure that a solution exists.19 By then choosing the positions $r_{p}$ we can extract the charges, $z_{p}$. The positions of the image particles may be arbitrary chosen, with the exceptions mentioned earlier, however in this work we use $c_{p}=q^{-p}$. That is, the position of image particle $p+1$ ($r_{p+1}$) is the mirror image of the position of image particle $p-1$ ($r_{p-1}$) in the position of image particle $p$ ($r_{p}$), i.e. $r_{p}^{2}=q^{-2p}r^{2}=q^{-(p-1)}rq^{-(p+1)}r=r_{p-1}r_{p+1}$. This choice is closely related to the method of image charges20, 21. Together with Eq. 7, the above assumptions yield the solution for the image charges as

[TABLE]

where ${P\brack p}_{q}$ is the q-analogue of the binomial coefficients22, 23. For details see the supplemental information (SI), Section S1. By using Eq. 8 it is further possible to present an expression for the aggregated potential $V_{q}$, and the modified interaction-tensor as

[TABLE]

The expression within the parenthesis is composed of the contributions from the original particle (the $1$) and all image particles (the sum). Rearranging and simplifying23 this expression gives

[TABLE]

where $(a;q)_{P}$ is the q-Pochhammer Symbol. In Fig. 2 we show this multiplicative modification to the original interaction-tensor for different $P$. Since the multiplicative term is a $q$-analogue of the Pochhammer symbol we address the proposed potential as the q-potential.

By re-writing Eq. 10 as

[TABLE]

and specifically note the $(1-q)^{P}$-term, we acknowledge that the number of cancelled higher order moments, $P-1$, equals the number of higher order derivatives with respect to $|\bm{r}|$ (or equivalently $q$) that are zero at the cutoff. Finally, using the already described theory, but with image particle positions $c_{p}=q^{-sp}$ we retrieve the short-ranged function

[TABLE]

where $s\in\{1,2,...\}$, $s-1$ higher order derivatives are cancelled also at $q=0$. This is straight-forward to acknowledge as the lowest power of $q$ in $(q^{s};q^{s})_{P}$ (save $q^{0}=1$) is $q^{s}$. Consequently, the $(s-1)$th derivative is the highest order derivative to vanish at $q=0$ (given that $P>0$). The expression in Eq. 12 can now be tuned in how many derivatives to cancel at $q=0$ (using $s$) and at $q=1$ (using $P$). Due to the physical interpretation of the image charge positions20, 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, yet 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 the SI Section S2 the detailed derivation is presented which culminates into

[TABLE]

2.3 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 previsously noted that to cancel every moment within a cutoff region gives the right-hand side of Eq. 7 as the zero-vector. If, however, there would be a non-zero moment within this region, one gets Eq. 14 where $\Psi^{(p)}$ is the $p$th higher order moment in the region projected onto the $r$-component of ${\bf r}$. The solution to this equation is more comprehensive than the solution to Eq. 7, yet solvable, at least given the mentioned restrictions on $r_{p}$ and $M$. Here we have not explicitly expanded the given framework for such cases, yet we highlight the possibility to develop a $q$-potential for non-zero moment environments.

[TABLE]

2.4 Choosing the cutoff

It is desirable to use a minimal cutoff to speed up 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. Yet, most systems do display a certain local anisotropy, or equivalently, non-vanishing higher order moments. Accordingly, the cutoff needs to enclose this whole space. To ensure no inter-cell correlations between anisotropic regions, it should not exceed one fourth of the shortest cell-length24. 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.

3 Methods

3.1 Molecular Dynamics Simulation

All simulations were performed in the isobaric-isothermal ensemble at pressure $\mathcal{P}=$1 bar and temperature $T=$298.15 K, using OpenMM 7 25. The pressure was kept constant using a Monte Carlo barostat26 and we used a Langevin integrator 27, with a friction coefficient of 1.0 ps and a time step of 2.0 fs, to keep constant temperature. Before production runs, the systems were energy minimized, then equilibrated. The reference Ewald summation/PME simulations used a fractional force error tolerance of $5\cdot 10^{-4}$.

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

The diffusion coefficient defined in Eq. 15 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 micro-canonical 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 micro-canonical production simulation also provided the base for the analysis of the energy.

[TABLE]

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

3.2 Kirkwood-Buff Theory

The Kirkwood-Buff (KB) theory provides a general way to obtain macroscopic properties of a solution from its microscopic properties 29. The central property is the Kirkwood-Buff integral (KBI) between components $A$ and $B$, defined as

[TABLE]

where $g_{AB}^{\mu VT}(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 procedure32 in which two corrections33, 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

[TABLE]

Here $G_{cc}$ and $G_{wc}$ are the corresponding truncated, corrected KBIs obtained from simulations, $\hat{G}^{*}_{cc}(R)$ and $\hat{G}^{*}_{wc}(R)$ respectively (see Appendix A). Subscripts $c$ and $w$ denote electrolyte and water respectively, $a_{c}=\gamma_{c}\rho_{c}$ is the electrolyte activity, $\gamma_{c}$ is the molar mean activity coefficient of the electrolyte, and $\rho_{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 $g_{Oc}(r)$. Experimental activity derivatives for the simulated electrolytes were calculated from previously reported activity coefficients36, using the fitting procedure described in Appendix A. For the simulations, we used one force field37 for the Na+ and I- ions whilst another38 was used for Cl-. The parameters of these force fields are presented in Appendix A.

4 Results

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

4.1 Bulk water-system

A representation of the radial distribution functions between water oxygen atoms are presented in Fig. 3, 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 Ewald. However, by using $q(P=\infty)$ and $R_{c}\leq 1.28$ nm we diverge from the Ewald-like results retrieved by intermediate $P$-values.

The SP1-approach cancels one derivative at the cutoff and yields results more closely resembling Ewald than the similar $q(P=2)$-potential. Yet, the $q(P=3)$-potential yields results even closer to Ewald. The SP3-approach, which cancels three derivatives at the cutoff, produces among the most Ewald-like results akin the $q(P=5)$- and $q(P=7)$-potentials.

The density, dielectric constant (see Section S3 in the SI), Kirkwood factor $G_{K}$ (see Eq. 18), standard deviation of the total energy, and diffusion coefficient, for the different potentials are presented in Table 2. 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 $\sim 2$, nearly all potentials (for $R_{c}\geq 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\neq 1$ are consistent within the range $3.0\pm 0.1$, except $q(P=6)$ which for $R_{c}=0.96$ nm gives $G_{K}=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 yields particularly low values, even compared to Ewald. Yet, a seemingly odd behavior is the increase in $\sigma_{E}$ for $P=4$ and $P=5$ while expanding the cutoff from $R_{c}=1.28$ nm to $R_{c}=1.60$ nm. A similar trend is also observed for the Ewald summation result. The diffusion coefficients are all similar to results of previous studies42, 43 with $D=2.5\pm 0.3$ for all pair potentials with $R_{c}\geq 1.28$ nm.

Finally as stated in Section 2.4 we note that the cutoff should obey $R_{c}<\min(L_{x},L_{y},L_{z})/4$ and therefore, since the box-length is $\sim 3.9$ nm, artifacts may be present in systems using excessive cutoffs. However, in the SI Section S4 we present simulation results based on larger systems and see no effects not previously pointed out.

[TABLE]

4.2 Salt solutions

Fig. 5 shows activity derivatives for NaCl and NaI solutions and the $q(P=3)$-potential generally gives the most accurate results as compared to PME. Yet, the presented error-bars for all methods do overlap which indicates comparable results. We have performed simulations using the $q(P=1)$-potential with poor outcome (akin for the RDF), wherefore we conclude 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.

Fig. 6 shows RDFs between water and salt for the different potentials. Similarly to the oxygen-oxygen RDFs (Fig. 3), 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=\infty$ again yields results showing increased deviation from PME.

5 Discussion

For the simulated water-system using $R_{c}=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 cutoff45. There is thus a limit on how many cancellations are appropriate. Though 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 (Fig. 4) is more sensitive to electrostatic interactions than the radial distribution function (Fig. 3). For example, in previous studies using the Wolf formalism it has been shown that suboptimal combinations of cutoff and damping parameter can lead to unphysical results, where water dipoles align towards 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 Fig. 4) at $R_{c}=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 $R_{c}=0.96$ nm. Ergo, this suggests $R_{c}=1.28$ nm to be the lowest cutoff to give Ewald-like results on every distance. Further, the discrepancies in Fig. 3 and Fig. 4 seems to imply deviations from Ewald in densities, which is most clearly seen for the results using $R_{c}=0.96$ nm.

From a theoretical perspective we expect all long-ranged interactions to vanish by using $P=\infty$, i.e. 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 $P$ too large results in over-damping 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 $R_{c}\to\infty$ the $q$-potential collapses to the unmodified exact Coulomb potential which is independent of $P$. Conversely, for $R_{c}\to 0$ no $P$-values give a proper description of the electrostatics and the optimal interval for $P$ varies with the cutoff. Optimal values are about $P=5$ for $R_{c}=0.96$ nm; $P\in\{4,5,6\}$ for $R_{c}=1.28$ nm; and $P\in\{3,4,...\}$ for $R_{c}=1.6$ nm. Note that even the most precise results using $R_{c}=0.96$ nm are still inferior to a mediocre parameterization at $R_{c}=1.6$ nm. The optimal interval seemingly widen with increasing cutoff, and as discussed in Section 2.1 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 $\mathcal{N}\cdot R_{c}^{3}$, where the prefactor is implementation dependent. Evaluation of the specifically used short-ranged function $\mathcal{S}(q)$ is part of this prefactor. The complementary error-function, ${\rm erfc}$, used in Ewald is commonly approximated47 using seven multiplications and one ${\rm 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)_{P}$ includes $2(P-1)$ multiplications which for any relevant choice of $P\lesssim 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 evaluated6 to $R_{c}\approx 1.6$ nm. We find that a lower cutoff can be used for the $q$-potential. Therefore, for all valid $P$ values and $R_{c}\leq 1.6$ nm we find that the $q$-potential is computationally less costly than both Ewald summation and PME, and scales with $P\cdot\mathcal{N}\cdot R_{c}^{3}$.

6 Conclusion

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 $\mathcal{O}(\mathcal{N})$. 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, i.e. number of cancelled moments. We find that $R_{c}=1.28$ nm ($\sim 4$ molecular diameters) and $P=\{4,5,6\}$ gives 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.

\acknowledgement

We thank LUNARC in Lund for computational resources, and for financial support: Södra Research foundation; the Swedish Research Council; and the Swedish Foundation for Strategic Research. Finally, we thank the referees for valuable comments.

\suppinfo

The Supporting Information is available free of charge on the ACS Publications website at DOI:XYZ and includes: derivation of image charges, derivation of the self-energy, derivation of the dielectric constant expression, and simulation results of larger systems.

The scripts used to produce the data for this work can be found at https://github.com/mlund/SI-qpotential. An implementation of the derived potential as well as other truncation based schemes are also available in the CoulombGalore C++ package48.

Appendix A Kirkwood-Buff Theory

For the simulated $\mathcal{N}\mathcal{P}T$ ensemble, we considered the following running integral

[TABLE]

where the hat indicates that RDFs in the $\mathcal{N}\mathcal{P}T$ ensemble is used in place of the ones in the grand canonical ensemble, and $R$ is the distance of truncation. The excess coordination number of component $B$ around component $A$, i.e. the excess number of particles $B$ around particles $A$ compared to bulk composition, is obtained as $\Delta N_{AB}=\rho_{B}G_{AB}$, where $\rho_{B}$ is the average number density of particles $B$ in the system.

Due to the finite size of simulated systems, some complications arise when applying KB theory. Firstly, the calculated RDFs usually does not converge within the simulation cell. This is most easily understood considering the excess coordination number $\Delta N_{AB}$. In a finite system containing non-ideal particles, the excess coordination numbers between the particles are non-zero due to interactions. Since the system is finite, an excess or deficit of for instance particles $B$ around particles $A$ will cause a corresponding deficit or excess of particles $B$ in the bulk. Hence the bulk composition, $\rho_{AB}^{(bulk)}$, will be different from $\rho_{B}$ which is used for normalization according to $g_{AB}(r)=\rho_{AB}(r)/\rho_{B}$. To make the RDFs asymptotically approach unity, we used the following correction factor 31, 33 which scales the RDF so that it is normalized with $\rho_{AB}^{(bulk)}$ rather than $\rho_{B}$

[TABLE]

Here $N_{B}$ is the total number of particles $B$, $V(r)$ is the volume of a sphere with radius $r$, $V_{cell}$ is the volume of the simulation cell, and $\delta_{AB}$ is the Kronecker delta. In the expression, the numerator is the number of particles $B$ that would occupy the volume beyond $r$ if the composition would be uniform, while the denominator approximates the true number of particles $B$ in the same volume.

[TABLE]

Secondly, for closed systems one cannot reduce the double integral in Eq. 21 to the single integral in Eq. 16 since the integration domain of $r$ is no longer independent of $\mathbf{r_{1}}$. To account for this, a method has been suggested34 where the exact KBI is found by applying a geometry dependent factor to the single integral and evaluating the expression in the thermodynamic limit, $1/R\rightarrow 0$. Alternatively 34, it has been suggested that an extrapolated expression of the exact KBI can be obtained by applying the following factor directly in the single integral, which is the method we used,

[TABLE]

To obtain approximate KBIs accounting for the system being both closed and of finite size, we simultaneously apply the two corrections32 according to

[TABLE]

where the star indicates that correction factors are applied. For the activity derivatives, the integral was then evaluated at a value of $R$ where convergence is obtained, corresponding to the distance from the center particle at which bulk composition is reached.

The experimental activity coefficients36 were used to find the values of the parameters $b_{j}$ in the following fitting function44, relating the activity coefficient to the concentration on the molal scale,

[TABLE]

The obtained values of $b_{j}$ are reported in Table LABEL:tab:S1.

Since the experimental activity coefficients used in this study are on the molal scale, we converted the concentration scale to be in the unit of number density in order to calculate the activity derivatives as $a_{c}^{\prime}=\left(\partial\ln a_{c}/\partial\ln\rho_{c}\right)_{\mathcal{P}T}$. For the conversion, we used experimental density data49, 50.