Lesson 5.1 · 14 min

Coulomb's law in a periodic box

Coulomb interactions decay so slowly that a plain cutoff gives wrong answers, and in a periodic box every ion interacts with infinitely many images. The rock-salt Madelung constant shows the problem, and why neutral groups and damping help.

Why the 1/r Coulomb interaction cannot simply be cut off like Lennard-Jones
Conditional convergence: the result depends on the order of summation (Madelung constant)
Neutral groups, reaction field and damped-shifted (Wolf) sums as approximations, leading to Ewald

The idea

For Lennard-Jones we used a cutoff of about 2.5σ and a small tail correction (Lesson 2.5). That works because the number of neighbours at distance r grows as r², while the LJ energy falls as r⁻⁶: the far region contributes almost nothing. Coulomb energy falls only as 1/r. The far shells hold more and more charges, and their total contribution does not die away at all. The sum is finite only because positive and negative charges cancel.

Conditional convergence. The classic example is the Madelung constant of rock salt, M = 1.7476: the electrostatic energy of an ion in an infinite NaCl crystal. Add up the ions inside a growing sphere and the answer swings wildly (2.1, −4.2, 8.0, …) and never settles, because each new shell is not neutral. Add them in neutral blocks (expanding cubes with boundary ions shared, Evjen's method) and it converges within a few shells. For such sums, the answer depends on how you add the terms.

In a simulation with periodic boundaries (Lesson 2.4), each charge interacts with every other charge and with all their periodic images. A plain cutoff of 1–1.2 nm makes the included region randomly charged and produces known artefacts: wrong ion–ion structure, artificial ordering of water near the cutoff, and wrong dielectric behaviour.

Quick fixes keep a cutoff but soften its effect:

Charge groups: include or exclude whole neutral groups (old GROMACS, GROMOS).
Reaction field: treat everything beyond the cutoff as a continuum with dielectric constant ε_RF.
Damped, shifted sums (Wolf, DSF): damp 1/r with erfc(αr), shift it to zero at the cutoff, and neutralise the cutoff sphere. Cheap and often accurate for dense systems.

The proper solution, used in almost all modern simulations, is Ewald summation and its fast version PME (Lessons 5.2–5.3). It splits 1/r into a short-ranged part summed directly and a smooth part summed exactly in Fourier space.

Key idea. Never use a plain cutoff for charged or polar systems. The Coulomb sum over a periodic system is only conditionally convergent; use Ewald/PME, or at least a damped-shifted or reaction-field method.

The algorithm

The Madelung constant three ways (reference ion +1 at the origin, nearest-neighbour distance 1):

# 1. plain sphere: does not converge
M = -Σ_{r_j < R} q_j / r_j

# 2. neutral cubes (Evjen): converges fast
M = -Σ_{|x|,|y|,|z| ≤ L} w_j q_j / r_j
    w = 1 inside, 1/2 on faces, 1/4 on edges, 1/8 at corners

# 3. Wolf: damp, shift and neutralise
e_i = ½ Σ_{r_j < Rc} q_i q_j [erfc(α r_j)/r_j - erfc(α Rc)/Rc]
      - [erfc(α Rc)/(2 Rc) + α/√π] q_i²

The shift term makes each pair interaction zero at Rc. Together with the self term, the cutoff sphere behaves as if it were neutral.
α trades accuracy at short range against convergence at long range; α ≈ 0.2–0.3 Å⁻¹ with Rc ≈ 10–12 Å is common for ionic liquids and molten salts.
In GROMACS, coulombtype = Reaction-Field or PME; plain Cut-off is a reaction field with ε_RF = 1. LAMMPS: pair_style coul/wolf, coul/dsf, or coul/long with kspace_style pppm.