Ewald summation

Background

Different from the short-ranged interactions, the electrostatic interaction describing by the coulomb law will not vanish beyond specific range for a dense system as:

(2)\[U_{\mathrm{ele}}(r+\mathrm{d}r) = \frac{q_sq_n}{4\pi\varepsilon_0 r}\times n_0 \pi r^2 dr = n_0 \frac{q_sq_n}{4\varepsilon_0} r \mathrm{d}r\]

According to Equation (2), the interaction energy between central particles and their neighboring particles within the range of \([r, r + \mathrm{d}r]\) increases linearly as the distance increases. Therefore, it is not feasible to simplify the calculation of electrostatic interaction by setting a cutoff radius and ignoring particles beyond it. Additionally, the presence of periodic boundary conditions (PBC) results in an infinite sum that cannot be easily calculated, further adding to the complexity of the electrostatic interaction calculation.

Fig. 6 Idea of ewald summation [4]

Ewald summation is a mathematical method utilized in computational chemistry and physics to precisely calculate the long-range electrostatic interactions between charged particles in a periodic system. The fundamental concept involves transforming the original delta-like charge distribution into the superposition of a screening distribution and a Gaussian-like distribution, as depicted above. The screening distribution is composed of the original distribution and a Gaussian with an opposite charge sign, similar to dipoles. Approximately, the energy is altered to:

(3)\[U_{\mathrm{ele}}(r+\mathrm{d}r) = \frac{q_sq_n}{4\pi\varepsilon_0 r^3}\times n_0 \pi r^2 dr = n_0 \frac{q_sq_n}{4\varepsilon_0r} \mathrm{d}r\]

As the distance increases, the energy now diminishes. Precisely, the energy of the screening distribution can be expressed as follows:

(4)\[\begin{split}\begin{equation} \begin{cases} \displaystyle \phi _{\text{G}}(r)=\frac{q}{4\pi \varepsilon _0r}\text{erf(}\alpha r\text{),erf(}x)=\frac{2}{\pi ^{\text{1/}2}}\int_0^x{\text{e}^{-x^2}}\text{d}x \\\\ \displaystyle \phi _{\delta}(r)=\frac{q}{4\pi \varepsilon _0r} \end{cases} \label{equation} \end{equation}\end{split}\]

(5)\[E_{\text{short}}=\sum_{i\ne j,\lVert \mathbf{r}_i-\mathbf{r}_i \rVert \le r_c}^{N_{\text{p}}}{\frac{q_iq_j}{4\pi \varepsilon _0 \lVert \mathbf{r}_i-\mathbf{r}_j \rVert }\left[ 1-\text{erf}\left(\alpha \lVert\mathbf{r}_i-\mathbf{r}_j\rVert\right) \right]}\]

Moreover, for the Gaussian-like distribution, we can substitute the summation in real space with the summation in reciprocal space, since the Gaussian distribution is periodic owing to the presence of PBC. First, considering the Poisson equation in reciprocal space:

(6)\[\lVert \mathbf{k} \rVert ^2\tilde{\phi}(\mathbf{k})=\frac{\tilde{\rho}(\mathbf{k})}{\varepsilon _0}\]

, the electric potential of Gaussian-like distribution in the reciprocal space expresses as:

(7)\[\tilde{\phi}_{\text{G}}(\mathbf{k})=\frac{1}{\varepsilon _0\lVert \mathbf{k} \rVert ^2}\int{\sum_{j=1}^{N_{\text{p}}}{q_j}}G(\mathbf{r}-\mathbf{r}_j\text{)e}^{-i\mathbf{kr}}\text{d}\mathbf{r}=\frac{1}{\varepsilon _0\lVert \mathbf{k} \rVert ^2}\sum_{j=1}^{N_{\text{p}}}{q_j}\exp \left( -\frac{\lVert \mathbf{k} \rVert ^2}{4\alpha ^2}-i\mathbf{kr} \right)\]

Next, we utilize the Poisson summation to convert the potential in reciprocal space into real space and give the potential energy:

(8)\[\begin{split}\begin{array}{rcl} \displaystyle\phi _{\text{G}}(\mathbf{r}) & = & \displaystyle \sum_{\mathbf{n}\in \mathbb{Z}^3}{\phi _{\text{G}}}'(\mathbf{r}+\mathbf{nL})=\frac{1}{V}\sum_{\mathbf{k}}{\tilde{\phi}_{\text{G}}}(\mathbf{k}\text{)e}^{i\mathbf{kr}}\\\\ & = & \displaystyle \frac{1}{\varepsilon _0V}\sum_{\mathbf{k}}{\sum_{j=1}^{N_{\text{p}}}{\frac{q_j}{\lVert \mathbf{k} \rVert ^2}}}\exp \left[ -\frac{\lVert \mathbf{k} \rVert ^2}{4\alpha ^2}-i\mathbf{k}(\mathbf{r}-\mathbf{r}_j) \right] \end{array}\end{split}\]

(9)\[E_{\text{G}}=\frac{1}{\varepsilon _0V}\sum_{\mathbf{k}}{\sum_{i,j}^{N_{\text{p}}}{\frac{q_iq_j}{\lVert \mathbf{k} \rVert ^2}}}\exp \left[ -\frac{\lVert \mathbf{k} \rVert ^2}{4\alpha ^2}-i\mathbf{k}(\mathbf{r}_i-\mathbf{r}_j) \right]\]

The summation over \(\mathbf{k}\) converges rapidly owing to the presence of \(\lVert \mathbf{k} \rVert ^2\) in the denominator and the exponential term. However, there are still two issues with Equation (9): (a) the value of the summation term when k = 0, and (b) the existence of self-interaction. The first question can be resolved when simulating an electroneutral system [5]. In contrast to disregarding all \(i=j\) pairs in Equation (5), we must manually eliminate the self-interaction term for the second question. The self-interaction term describes the interaction between the delta-like charge and the surrounding same charged Gaussian:

(10)\[E_{\text{self}}=\sum_{i=1}^{N_{\text{p}}}{\lim_{r\rightarrow 0}}q_i\phi _{\text{G}}(r)=\sum_{i=1}^{N_{\text{p}}}{\frac{\alpha q_{i}^{2}}{2\varepsilon _0\pi ^{\text{3/}2}}}\]

Then the total electric potential energy can be written as:

(11)\[U_{\text{ele}}=E_{\text{short}}+E_{\text{long}}-E_{\text{self}}\]

In practical applications, the Fast Fourier Transform (FFT) algorithm is utilized to hasten the calculation of long-range electrostatic interactions, and its implementation necessitates the discretization (griding) of space. Two common algorithms, namely the the Particle Mesh Ewald [6, 7] (PME) method and Particle-Particle Particle Mesh [8, 9] (PPPM or P3M) method, have been introduced to tackle this issue. In MDPy, the PME method is implemented to evaluate electrostatic interactions.

Particle Mesh Ewald

Fig. 7 Illustration of: (Right) B-spline interpolation, (Left) Electric potential.

Ewald coefficient

As shown in Equation (4), the Ewald coefficient \(\alpha\) determines the width of the Gaussian-like charge distribution, which in turn affects the electric potential expression of the screening distribution and the choice of the cutoff radius. In the PME method, the user specifies the cutoff radius \(r_c\) and error tolerance \(e\) as the input parameters to calculate the Ewald coefficient. The Ewald coefficient must satisfy the following equation:

(12)\[e = 1-\text{erf}\left(\alpha \lVert\mathbf{r}_i-\mathbf{r}_j\rVert\right)\]

In MDPy, the Ewald coefficient is calculated by solving Equation (12) via the Newton method.

Bspline interpolation

As FFT requires a discretization of space, the particles’ charge should be assigned to the grid points. B-spline interpolation is a widely used technique in numerical analysis and computer graphics to approximate complex functions with a smooth curve. It involves constructing a piecewise polynomial function that passes through a given set of control points, while ensuring a high degree of smoothness and continuity between adjacent curves.

B-splines are defined by a set of basis functions, which are polynomial functions of a fixed degree that are defined over a specific interval. By combining these basis functions with a set of control points, a smooth curve can be constructed that accurately approximates the original function.

Fig. 8 Example of charge spreading in 2D for 4th order B-spline; circle sizes signify absolute values of charge contributions [10].

MDPy employs a 4th order B-spline curve to map the charge to a cubic grid. The above diagram provides a 2D illustration of this process.

Tip

More detailed information can be found in reference [10].