Introduction

Drop a handful of charged atoms into a periodic box. Each atom feels a Coulomb force from every other atom in the box and from every one of their infinitely-many periodic images. Write that down as a sum and you get something alarming:

E=12∑i≠jqiqjrijE = \frac{1}{2}\sum_{i \neq j} \frac{q_i q_j}{r_{ij}}

The trouble is that the sum barely converges — the Coulomb potential decays as 1/r1/r, just barely fast enough. Adding terms in the wrong order gives a completely different answer. Your computer would take longer than the age of the universe to get a reliable result from the raw series.

In 1921, the crystallographer Paul Peter Ewald found a way out. He did not try to sum the series differently; he split it into two new series, each of which converges much faster than the original. Real-space terms handle the nearby atoms; reciprocal-space terms handle the far-away ones. Together they reproduce the original sum exactly, yet each half is easy to compute.

That split-and-conquer idea became the heartbeat of every classical molecular-dynamics and Monte Carlo simulation of ionic systems — from salt crystals to proteins in water.

Split the Series

The demo below shows a 1-D version of Ewald's idea. A point charge sits at the origin, and we sum 1/r1/r over its periodic images out to NN neighbors on each side — a sum that converges very slowly.

Ewald's trick adds and subtracts a Gaussian screening charge around each image. This splits the sum into a real-space part (short-range, fast-decaying) and a reciprocal-space part (handled by a cosine series). Use the slider to tune the Gaussian width α\alpha and watch both halves converge to the same total.

<p class="hint">We split the slow series for pi^2/6 = 1/1 + 1/4 + 1/9 + ... into two fast halves. The real-space part (teal) uses a Gaussian screen; the reciprocal-space part (orange) uses its Fourier dual. Together they hit the exact value (dashed line) in just a few terms.</p>
<div class="controls">
  <label for="alpha-slider">Gaussian width alpha = <strong id="alpha-val">0.5</strong></label>
  <input id="alpha-slider" type="range" min="0.2" max="2.0" step="0.1" value="0.5">
  <label for="n-slider">Terms N = <strong id="n-val">20</strong></label>
  <input id="n-slider" type="range" min="3" max="60" step="1" value="20">
</div>
<canvas id="chart" width="560" height="240"></canvas>
<div class="legend">
  <span class="dot raw"></span> Raw sum (slow)
  <span class="dot real"></span> Real-space part
  <span class="dot recip"></span> Reciprocal-space part
  <span class="dot total"></span> Ewald total
</div>
<div class="status" id="status"></div>
<div class="btns">
  <button id="btn-reset" type="button" class="ghost">Reset</button>
</div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; }
.hint { font-size: .9rem; color: #444; margin: 0 0 .8rem; line-height: 1.5; }
.controls { display: flex; flex-direction: column; gap: .4rem; margin-bottom: .8rem; font-size: .88rem; }
.controls label { display: flex; align-items: center; gap: .6rem; }
.controls input[type=range] { flex: 1; max-width: 320px; }
canvas { display: block; width: 100%; max-width: 560px; border: 1px solid #d0d8e0; border-radius: 8px; background: #f8fafc; }
.legend { display: flex; flex-wrap: wrap; gap: .5rem 1rem; font-size: .82rem; margin: .5rem 0; align-items: center; }
.dot { display: inline-block; width: 12px; height: 12px; border-radius: 50%; }
.dot.raw   { background: #aaa; }
.dot.real  { background: #2a9d8f; }
.dot.recip { background: #e76f51; }
.dot.total { background: #1d3557; }
.status { font-size: .95rem; font-weight: 600; min-height: 1.4em; margin: .3rem 0; color: #1d3557; }
.btns { display: flex; gap: .5rem; flex-wrap: wrap; margin-top: .4rem; }
button { font: 600 14px system-ui, sans-serif; padding: .45rem .9rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }

const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Ewald splitting demo: pi^2/6 = sum(1/k^2) split into two fast-converging halves.
// Real part: sum exp(-alpha^2*k^2)/k^2. Reciprocal part: sqrt(pi)*alpha - alpha^2/2 + Fourier correction.
// Real-space part: Gaussian screen kills large-k terms fast.
// Reciprocal-space part: Poisson-summation dual, converges super-fast via erfc.

const canvas = document.getElementById('chart');
const ctx = canvas.getContext('2d');
const W = canvas.width, H = canvas.height;
const PAD = { top: 16, right: 16, bottom: 32, left: 52 };

// Abramowitz & Stegun 7.1.26 erfc approximation, max error ~1.5e-7.
function erfc(x) {
  if (x < 0) return 2 - erfc(-x);
  const t = 1 / (1 + 0.3275911 * x);
  const p = t * (0.254829592 + t * (-0.284496736 + t * (1.421413741 + t * (-1.453152027 + t * 1.061405429))));
  return p * Math.exp(-x * x);
}

// The exact target value pi^2/6 ~ 1.6449 (Basel problem).
const TARGET = Math.PI * Math.PI / 6;

// Raw partial sum: cumulative sum of 1/k^2 (converges as O(1/N)).
function rawPartial(n) {
  let s = 0;
  for (let k = 1; k <= n; k++) s += 1 / (k * k);
  return s;
}

// Real-space part: sum exp(-alpha^2*k^2)/k^2 (Gaussian kill, very fast).
function realPartial(n, a) {
  let s = 0;
  for (let k = 1; k <= n; k++) s += Math.exp(-a * a * k * k) / (k * k);
  return s;
}

// Reciprocal-space part: pi^2/6 - real_inf, expressed as Poisson dual series.
function recipSum(mMax, a) {
  let s = Math.sqrt(Math.PI) * a - a * a / 2;
  for (let m = 1; m <= mMax; m++) {
    const c = Math.PI * m / a;
    s += 2 * Math.pow(Math.PI, 1.5) * (Math.exp(-c * c) / c - Math.sqrt(Math.PI) * erfc(c));
  }
  return s;
}

function computeSeries(nMax, a) {
  const raw = [], real = [], recip = [], total = [];
  for (let n = 1; n <= nMax; n++) {
    raw.push(rawPartial(n));
    const r = realPartial(n, a);
    const k = recipSum(n, a);
    real.push(r);
    recip.push(k);
    total.push(r + k);
  }
  return { raw, real, recip, total };
}

function drawChart(data) {
  ctx.clearRect(0, 0, W, H);
  const all = [...data.raw, ...data.real, ...data.recip, ...data.total, TARGET];
  const mn = Math.min(...all) - 0.05, mx = Math.max(...all) + 0.05;
  const n = data.raw.length;
  const sx = (i) => PAD.left + (i / (n - 1)) * (W - PAD.left - PAD.right);
  const sy = (v) => PAD.top + (1 - (v - mn) / (mx - mn)) * (H - PAD.top - PAD.bottom);

  // Draw horizontal grid lines with y-axis labels.
  ctx.strokeStyle = '#dde3ea';
  ctx.lineWidth = 1;
  for (let t = 0; t <= 4; t++) {
    const yv = mn + t * (mx - mn) / 4;
    const y = sy(yv);
    ctx.beginPath(); ctx.moveTo(PAD.left, y); ctx.lineTo(W - PAD.right, y); ctx.stroke();
    ctx.fillStyle = '#777'; ctx.font = '10px system-ui';
    ctx.fillText(yv.toFixed(2), 2, y + 4);
  }

  // Dashed reference line at pi^2/6.
  ctx.strokeStyle = '#1d3557';
  ctx.lineWidth = 1;
  ctx.setLineDash([4, 3]);
  ctx.beginPath();
  ctx.moveTo(PAD.left, sy(TARGET));
  ctx.lineTo(W - PAD.right, sy(TARGET));
  ctx.stroke();
  ctx.setLineDash([]);

  // Draw all four partial-sum curves.
  const series = [
    { data: data.raw,   color: '#aaa',    lw: 1.5 },
    { data: data.real,  color: '#2a9d8f', lw: 2 },
    { data: data.recip, color: '#e76f51', lw: 2 },
    { data: data.total, color: '#1d3557', lw: 2.5 }
  ];
  for (const s of series) {
    ctx.strokeStyle = s.color; ctx.lineWidth = s.lw;
    ctx.beginPath();
    s.data.forEach((v, i) => { i === 0 ? ctx.moveTo(sx(i), sy(v)) : ctx.lineTo(sx(i), sy(v)); });
    ctx.stroke();
  }

  // X-axis tick labels (number of terms N).
  ctx.fillStyle = '#555'; ctx.font = '11px system-ui';
  const step = Math.max(1, Math.round(n / 6));
  for (let i = 0; i < n; i += step) {
    ctx.fillText(i + 1, sx(i) - 4, H - PAD.bottom + 14);
  }
}

function update() {
  const alpha = parseFloat(document.getElementById('alpha-slider').value);
  const nMax  = parseInt(document.getElementById('n-slider').value, 10);
  document.getElementById('alpha-val').textContent = alpha.toFixed(1);
  document.getElementById('n-val').textContent = nMax;
  const data = computeSeries(nMax, alpha);
  drawChart(data);
  const finalTotal = data.total[data.total.length - 1];
  const err = Math.abs(finalTotal - TARGET);
  document.getElementById('status').textContent =
    i18nFmt('Ewald total = %val% (exact pi^2/6 = %ref%, error = %err%)', { val: finalTotal.toFixed(5), ref: TARGET.toFixed(5), err: err.toExponential(2) });
}

document.getElementById('alpha-slider').addEventListener('input', update);
document.getElementById('n-slider').addEventListener('input', update);
document.getElementById('btn-reset').addEventListener('click', () => {
  document.getElementById('alpha-slider').value = '0.5';
  document.getElementById('n-slider').value = '20';
  update();
});

update();
Result
Console

Notice how the raw partial sum zigzags slowly toward the limit, while the split parts converge almost immediately. That is the essence of Ewald's insight: rearrange, not eliminate, the difficult terms.

The Real Complexity

Ewald's math is elegant, but the complexity story is where it gets practically important.

  • Naive direct sum: O(N2)O(N^2) per time step — every pair of atoms, every step. For thousands of atoms that is billions of operations per nanosecond of simulation.
  • Classic Ewald: splitting into real-space and reciprocal-space sums brings the cost down to O(N3/2)O(N^{3/2}). This is the published 1921 result, and it was already a huge win for hand calculations on crystals.
  • Particle-Mesh Ewald (PME): by projecting charges onto a grid and using the Fast Fourier Transform (FFT) for the reciprocal part, the cost drops to O(Nlog⁡N)O(N \log N). This is the method inside virtually every modern molecular-dynamics engine (AMBER, GROMACS, NAMD).
  • Fast Multipole Method (FMM): a different philosophy — hierarchically group distant charges — reaches O(N)O(N) in theory, though with a larger constant.

The interplay between these approaches connects to deep ideas in fast multiplication: when your computation has structure (periodicity, smoothness, locality), you can exploit it algorithmically rather than brute-forcing through every pair.

The parameter α\alpha that splits the Gaussian is a tuning knob: large α\alpha shifts work to reciprocal space (more Fourier modes needed); small α\alpha shifts it to real space (longer cutoff radius). Optimal α\alpha balances both halves and minimizes total cost.

Where It Matters

Ewald summation is not a curiosity — it is load-bearing infrastructure for computational science:

  • Molecular dynamics of biomolecules: simulating proteins, DNA, and lipid membranes in water requires accurate long-range electrostatics at every femtosecond time step. PME makes this tractable.
  • Drug design: free-energy calculations that predict how a drug molecule binds to its target depend on correct electrostatic energies. Errors in the long-range sum corrupt the binding affinity.
  • Materials science: predicting the stability of ionic crystals (salts, ceramics, battery electrolytes) needs the Madelung energy — the exact quantity Ewald's method computes.
  • Plasma physics and astrophysics: the same periodic-image idea applies whenever you model a representative chunk of a much larger system.
  • Lattice-QCD: even quantum chromodynamics simulations of quarks and gluons on a spacetime lattice use Ewald-style techniques for long-range color charges.

The computational cost reduction from O(N2)O(N^2) to O(Nlog⁡N)O(N \log N) is not a minor speedup — it is the difference between simulating a handful of atoms for a few steps and simulating millions of atoms for microseconds, the timescale where biology actually happens.

Conclusion

Ewald's insight is a masterclass in a principle that appears everywhere in algorithms: if a computation is hard in one domain, transform it to another where it is easy, then combine the results.

He did not make the infinite Coulomb sum simpler — he made it converge twice, once quickly in real space and once quickly in reciprocal space, so that the two fast sums together give the same answer as the single slow one. A century later, his splitting identity still drives every production molecular-dynamics code on the planet.

The same philosophy — split a hard problem, solve each piece in its natural domain, recombine — echoes through fast multiplication and the Fast Fourier Transform. Whenever you see an O(Nlog⁡N)O(N \log N) where brute force would be O(N2)O(N^2), there is probably an Ewald-style insight hiding underneath.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/ewald-summation/Content licensed under CC BY-NC 4.0.