Lesson 5.2 · 18 min

Ewald summation

Ewald summation computes the full periodic Coulomb energy exactly. It splits 1/r into a short-ranged part, summed directly with a cutoff, and a smooth part, summed in Fourier space with a few wave vectors. A self-energy term completes it.

The Gaussian screening idea behind the erfc/erf split
The real-space, reciprocal-space and self terms, and their forces
How to choose α, the cutoff and k_max; why the total does not depend on α

The idea

Lesson 5.1 showed that the Coulomb sum over a periodic system converges only conditionally. Ewald (1921) found a way to rewrite it as two sums that both converge fast and absolutely.

The trick. Surround every point charge with a Gaussian cloud of the opposite charge. The point charge plus its cloud is screened: from outside it looks neutral, so its potential, erfc(αr)/r, dies off quickly. Sum these screened interactions directly in real space, with a cutoff and neighbour lists, as for Lennard-Jones.

Then add back Gaussian clouds of the original sign. Their potential, erf(αr)/r, is smooth everywhere: no singularity at r = 0. A smooth periodic function needs only a few Fourier components, so this part is summed in reciprocal space over wave vectors k = 2πm/L, using the structure factor S(k) = Σ q_j e^ik·r_j.

Finally, each point charge has been made to interact with its own Gaussian in the Fourier part. That artificial self-energy, −(α/√π) Σ q_i², is subtracted. (With the usual "tin-foil" boundary condition there is no extra surface term.)

The parameter α decides how the work is divided. A large α makes the real-space part short-ranged (small cutoff) but needs more k vectors; a small α does the opposite. The total does not depend on α: a powerful check. In practice r_c is set (about 1 nm), α is chosen so that erfc(αr_c) ≈ 10⁻⁵, and k_max is raised until the reciprocal error matches. The cost scales as N^3/2 at best; PME (Lesson 5.3) brings it to N log N.

Key idea. E = E_real (screened, short-ranged) + E_recip (smooth, Fourier) + E_self. The split parameter α moves work between the two sums but does not change the answer.

The algorithm

Energy and forces for a neutral cubic box of side L (units with 1/(4πε₀) = 1):

# real space (with cutoff rc and neighbour lists)
for pairs i < j with r = |r_ij (minimum image)| < rc:
    E_real += q_i q_j erfc(α r) / r
    F_ij    = q_i q_j [erfc(α r)/r + 2α/√π exp(-α² r²)] r_ij / r²

# reciprocal space
for k = 2π m / L, m ≠ 0, |m_x|,|m_y|,|m_z| ≤ k_max:
    A(k) = exp(-k² / 4α²) / k²
    S(k) = Σ_j q_j exp(i k·r_j)
    E_recip += (2π/V) A(k) |S(k)|²
    F_i     += (4π q_i/V) A(k) k Im[S*(k) exp(i k·r_i)]

E_self = -(α/√π) Σ_i q_i²
E = E_real + E_recip + E_self

Use k and −k symmetry to halve the reciprocal work, and update S(k) once per step for all atoms.
Excluded pairs (bonded 1–2, 1–3) are inside the reciprocal sum anyway, so their erf(αr)/r part must be subtracted in real space.
Check: the energy must be independent of α, and the forces must match a numerical derivative (both in the code tab).
GROMACS sets α from rcoulomb and ewald-rtol (default 10⁻⁵); LAMMPS from kspace_style ewald 1e-5. Plain Ewald is fine for a few hundred atoms; use PME/PPPM beyond that.