Lesson 5.3 · 16 min

Particle-mesh Ewald

Particle-mesh Ewald makes Ewald summation fast enough for millions of atoms. Charges are spread onto a grid, the Fourier sum is done with fast Fourier transforms, and forces are interpolated back. The cost drops from N^1.5 (or worse) to N log N.

The four PME steps: spread, FFT, multiply by the influence function, transform back
B-spline interpolation, the grid spacing and the interpolation order
Why PME scales as N log N, and the standard parameter choices (and PPPM)

The idea

In plain Ewald (Lesson 5.2) the reciprocal sum needs the structure factor S(k) = Σ_j q_j e^ik·r_j for every wave vector: N atoms × N_k vectors. For large systems that is the bottleneck.

The PME idea (Darden, York and Pedersen 1993; smooth PME by Essmann et al. 1995): if the charges sat on a regular grid, S(k) for all k at once would be a single fast Fourier transform, costing K³ log K³. So:

Spread: assign each charge to the nearby grid points with smooth B-spline weights (order 4 means 4×4×4 points per atom).
FFT the charge grid.
Multiply by the Ewald influence function e^−k²/4α²/k², corrected by a factor |b(k)|² that undoes the smoothing of the B-splines. The energy is the sum of this product with |FFT|².
Inverse FFT to get the smooth potential on the grid, and interpolate the forces back to the atoms with the derivatives of the same B-splines.

The real-space part stays as in Ewald (cutoff, neighbour lists). Because the FFT is cheap, α can be large and the real-space cutoff short (about 1 nm), and the total cost scales as N log N.

Accuracy is controlled by the grid spacing (finer is better) and the interpolation order (higher is better). The defaults in most codes, 0.10–0.12 nm and order 4 with a relative tolerance of 10⁻⁵, give forces accurate to about 10⁻³ relative, which is plenty for MD. PPPM (Hockney and Eastwood, used by LAMMPS) is a closely related method with an optimised influence function.

Key idea. PME = Ewald with the Fourier sum done on a grid by FFTs. Spread → FFT → multiply → inverse FFT → interpolate. Finer grids and higher orders give higher accuracy; the cost grows only as N log N.

The algorithm

Smooth PME, reciprocal energy (forces follow from the same grids):

u_i = K * r_i / L                         # scaled fractional coordinates
Q = zeros(K, K, K)
for each atom i:                           # spread with B-splines of order n
    for a, b, c in 0..n-1:
        Q[floor(u_i) - (a, b, c)] += q_i M_n(w_x + a) M_n(w_y + b) M_n(w_z + c)
F = FFT(Q)
for each grid frequency m ≠ 0:
    G(m) = exp(-π² |m/L|² / α²) / |m/L|² * |b_x|² |b_y|² |b_z|²
E_recip = Σ_m G(m) |F(m)|² / (2π V)
φ = inverse FFT(G * F)                     # smooth potential on the grid
F_i = -q_i Σ (gradient of the B-spline weights) * φ   # interpolate back

Choose K as a product of small primes (2, 3, 5) for fast FFTs; codes do this automatically from the spacing.
Spreading and interpolation cost N × n³; the FFT costs K³ log K³ with K³ ∝ N. Both are linear (or N log N).
In parallel runs the FFT needs all-to-all communication; GROMACS can put PME on separate ranks or the GPU.
GROMACS: coulombtype = PME, fourierspacing = 0.12, pme-order = 4, ewald-rtol = 1e-5. LAMMPS: kspace_style pppm 1e-5 or pppm/tip4p.