# 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:

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.

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:

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

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:

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

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

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:

Then the total electric potential energy can be written as:

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#

### 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:

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.

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].