Electrostatics

ExTinyMD's electrostatics standard library implements Ewald summation and the image-charge method (ICM) for dielectrically confined slabs. Every method has two interfaces:

• a core API that takes a plan object plus plain arrays of positions and charges, and requires no ExTinyMD type at all;
• an MD adapter that wraps the same plan object as an AbstractInteraction so it can be added to an MDSys and driven by simulate!.

The two-layer API

The core API is a plan object (built once, holding cutoffs, k-vectors, and scratch buffers) queried with coulomb_energy, coulomb_force and coulomb_force!:

using ExTinyMD, StaticArrays

n, L = 100, (20.0, 20.0, 20.0)
poses   = [SVector(rand() * L[1], rand() * L[2], rand() * L[3]) for _ in 1:n]
charges = [isodd(i) ? 1.0 : -1.0 for i in 1:n]

inter = Ewald3D(n, L; α = 0.5, s = 4.0)   # r_c = s/α = 8.0, below min(L)/2 = 10.0
E = coulomb_energy(inter, poses, charges)
F = coulomb_force(inter, poses, charges)

Nothing above is an ExTinyMD Atom, Boundary, or MDSys — poses and charges are plain Vectors. This is deliberate: the core layer can be tested, benchmarked, and used from a plain script without any MD framework at all.

The same inter object also plugs into an ordinary MDSys. The adapter (src/interactions/electrostatics/adapter.jl) implements update_acceleration! and energy for it exactly like any other interaction — gathering positions and charges from sys/info in particle-slot order and calling the same coulomb_* functions above:

using ExTinyMD, StaticArrays, Random

n, L = 100, 20.0
boundary = CubicBoundary(L)
atoms = create_atoms([(n ÷ 2, Atom(type = 1, mass = 1.0, charge = 1.0)),
                      (n ÷ 2, Atom(type = 2, mass = 1.0, charge = -1.0))])
info = SimulationInfo(n, atoms, (0.0, L, 0.0, L, 0.0, L), boundary; min_r = 1.0, temp = 1.0)

inter  = Ewald3D(n, (L, L, L); α = 0.5, s = 4.0)
finder = CellList3D(info, inter.short.r_c, boundary, 1)

sys = MDSys(
    n_atoms = n, atoms = atoms, boundary = boundary,
    interactions = [(inter, finder)],
    loggers = [TemperatureLogger(100; output = false)],
    simulator = VerletProcess(dt = 0.001),
)

simulate!(sys.simulator, sys, info, 1000)

update_acceleration! is now called for you every step; no other code path needs to know that the interaction is electrostatic. Internally, the adapter gathers positions and charges from sys/info in particle-slot order with two small helpers:

EwaldInteraction(short, long, n_atoms)

coulomb_energy(interaction, poses, charges; neighbor_list = nothing) -> T

coulomb_force(interaction, poses, charges; neighbor_list = nothing) -> Vector{SVector{3,T}}

coulomb_force!(F, interaction, poses, charges; neighbor_list = nothing) -> F

Ewald3D(n_atoms, L; α, s, ϵ = 1.0, ϵ_inf = Inf)

julia> using StaticArrays

julia> inter = Ewald3D(2, (10.0, 10.0, 10.0); α = 1.0, s = 4.0);

julia> poses = [SVector(0.0, 0.0, 0.0), SVector(5.0, 0.0, 0.0)];

julia> round(coulomb_energy(inter, poses, [1.0, -1.0]); digits = 6)
-0.021815

Ewald2D(n_atoms, L; α, s, ϵ = 1.0)

Periodic3D()

PeriodicQ2D()

Method selection

N is the particle count and K the number of reciprocal-space vectors kept below the cutoff k_c. PME3D and ICMPME3D trade that K dependence for an FFT-driven O(N log N) cost — see Particle mesh below for what that trade actually buys, and costs, in practice.

Use Ewald2D only as an accuracy reference for small systems: its O(N²K) direct double sum over all particle pairs makes it unsuitable for production work. For large quasi-2D systems without dielectric walls, reach for QuasiEwald.jl or SoEwald2D.jl instead.

Particle mesh

PME3D and ICMPME3D live in a package extension so that FINUFFT — a binary artifact — never lands on users who only want the direct Ewald methods above. Load it explicitly:

using ExTinyMD, FINUFFT

Without that using FINUFFT, ExTinyMD still defines PME3D and ICMPME3D (so isdefined sees them and dispatch to the wrong method is not possible), but calling either raises immediately:

ERROR: PME3D requires FINUFFT. Run `using FINUFFT` (and add it to your project)
to load ExTinyMD's particle-mesh extension.

(ICMPME3D raises the same sentence with its own name substituted.)

The most important fact: it is the same sum, not an approximation

PME3D does not compute an approximation to Ewald3D's reciprocal-space sum that happens to converge to the same answer — it computes the identical sum. FINUFFT's type-1 transform naturally returns values on a rectangular k-grid, but PME3DLong zeroes the Green's function D_k everywhere outside the spherical shell 0 < |k| ≤ k_c, which is exactly the k-set Ewald3DLong sums over directly. The two solvers are therefore summing over the same finite set of k-vectors with the same per-vector weight; any difference between them is floating-point roundoff, not truncation. Measured on this implementation:

The practical consequence: choosing between Ewald3D and PME3D (or between ICMEwald3D and ICMPME3D) is a question of cost, not accuracy. K grows with the system, so Ewald3D's O(N·K) cost grows faster than PME3D's O(N log N) as N increases; where exactly that starts to matter depends on your machine, your FFT size, and how tight s needs to be, so measure it on your own system rather than trusting a quoted crossover. Ewald3D remains useful even once PME3D is faster: it is a good correctness reference to check a new PME3D parameter choice against, and it has no FFT grid to size.

Parameters mean what they mean for Ewald3D

PME3D takes the same α and s as Ewald3D, with the same meaning (see Parameter selection below) and the same constraint:

r_c = s / α  <  min(L) / 2

Nothing about the particle-mesh grid changes this — it is still CellListMap building the real-space neighbour list for EwaldShort, so the same ArgumentError: UNIT CELL CHECK FAILED fires if it is violated.

using ExTinyMD, FINUFFT, StaticArrays, Random

n, L = 30, (12.0, 12.0, 12.0)
poses   = [SVector(rand() * L[1], rand() * L[2], rand() * L[3]) for _ in 1:n]
charges = [isodd(i) ? 1.0 : -1.0 for i in 1:n]

inter = PME3D(n, L; α = 0.8, s = 3.5)   # r_c = 4.375 < min(L)/2 = 6.0
E = coulomb_energy(inter, poses, charges)
F = coulomb_force(inter, poses, charges)

ICMPME3D: N_image and N_pad behave as they do for ICMEwald3D

ICMPME3D is ICMEwald3D with PME3DLong swapped in for Ewald3DLong underneath the same image-charge reflection and ELC padding, so N_image and N_pad carry over unchanged: N_image controls how many reflections the image series keeps, and N_pad controls how much the z-period is inflated before ELC treats the slab as ordinary triply-periodic Ewald. Underpadding is still the dominant failure mode, not an artifact of the particle-mesh solver — at N_image = 3, going from N_pad = 1 to N_pad = 2 takes the disagreement against the exact ICMEwald2D route from 9.4e-3 down to 4.7e-8. (These two figures were measured on the ICMEwald3D route, in the ICM guidance table below, not on ICMPME3D itself — but the two particle-mesh and direct-sum routes agree with each other to 4.0e-14, so the padding behaviour transfers unchanged.)

using ExTinyMD, FINUFFT, StaticArrays, Random

n, L = 8, (5.0, 5.0, 10.0)
poses   = [SVector(rand() * L[1], rand() * L[2], 2.0 + 6.0 * rand()) for _ in 1:n]
charges = [isodd(i) ? 1.0 : -1.0 for i in 1:n]

# r_c = s/α = 2.353 < min(L[1], L[2])/2 = 2.5
inter = ICMPME3D(n, L; α = 1.7, s = 4.0, γ = (0.3, 0.3), N_image = 3, N_pad = 2)
E = coulomb_energy(inter, poses, charges)

PME3DLong(n_atoms, L; α, s, ϵ = 1.0, ϵ_inf = Inf)

PME3D(n_atoms, L; α, s, ϵ = 1.0, ϵ_inf = Inf)

ICMPME3D(n_atoms, L; α, s, γ, N_image, N_pad, ϵ = 1.0)

Parameter selection

Every method here is built from a splitting parameter α and a dimensionless accuracy parameter s, related to the real- and reciprocal-space cutoffs by

r_c = s / α          k_c = 2 α s

(see ewald_cutoffs). Larger s means smaller truncation error — roughly exp(-s²) — in both spaces at once; s = 4 gives an error near 1e-7 and is a reasonable default. α only trades work between real and reciprocal space: increasing α shrinks r_c (fewer real-space pairs) while growing k_c (more k-vectors), and vice versa. The total energy must not depend on it.

That last property is also the recommended self-check for a new parameter choice: since the split point between real and reciprocal space is arbitrary, if two different α (at fixed s) give different total energies, something is wrong with the parameters or the setup, not with the physics. Measured on this implementation: Ewald3D's total energy is independent of α to about 1e-8 at s = 4 (in practice closer to 1e-9 in the cases checked), and Ewald2D agrees to the same order.

total3d(α) = coulomb_energy(Ewald3D(n, L; α = α, s = 4.0), poses, charges)
total3d(0.45), total3d(0.5), total3d(0.55)   # r_c = 8.89, 8.0, 7.27 (all < min(L)/2 = 10); should agree to ~1e-8

ewald_cutoffs(s, α) -> (r_c, k_c)

All the methods here also assume overall charge neutrality: Ewald summation of a non-neutral system is only conditionally convergent and its energy picks up a box-volume-dependent offset. Check with check_neutrality before trusting an energy from an unfamiliar system.

check_neutrality(charges; atol) -> Σq

r_c is not free of the box

Because r_c = s/α, the pair (α, s) is not independent of the box: r_c must be strictly less than half the smallest periodic side, i.e.

s / α  <  min(L) / 2

For Ewald2D, ICMEwald2D, ICMEwald3D and ICMPME3D only L[1] and L[2] enter this bound, since the z axis is not periodic (for ICMEwald3D and ICMPME3D this applies to the unpadded slab size passed in, not the internal z-padded box used for ELC). Violating the bound does not silently degrade accuracy — CellListMap, which builds the real-space neighbour list, refuses to construct a unit cell where the cutoff exceeds half a side, and raises

ArgumentError: UNIT CELL CHECK FAILED: unit cell dimension ... must be greater than 2*cutoff

This single constraint accounted for every failing parameter set in the original design draft of this library, so it is worth checking first whenever a construction call raises anything unexpected.

The prefactor asymmetry

Reading the source, the long-range prefactors look inconsistent with the short-range one. They are not — the 4π from Gaussian units is folded differently into each piece:

• EwaldShort (real space, both 3D and the ICM real-space kernel) carries 1/(4πϵ) explicitly.
• Ewald2DLong carries 1/ϵ — the 4π is already absorbed into the in-plane reciprocal-space expressions it evaluates.
• Ewald3DLong is written as 1/(2Vϵ) — again with the 4π folded in, alongside the 1/2 from the k-space sum running over a full k ↔ -k-symmetric set.

The α-independence check above is exactly the test that pins this normalisation: if the prefactors of the short- and long-range parts of a method disagreed, the total energy would visibly drift with α, since more of the energy would shift into whichever piece carries the wrong constant.

EwaldShort(n_atoms, L; α, s, ϵ = 1.0, convention = Periodic3D())

E_s = 1/(4πϵ) [ Σ_{i<j, r<r_c} q_i q_j erfc(α r_ij)/r_ij − (α/√π) Σ_i q_i² ]

short_energy(short, poses, charges; neighbor_list = nothing) -> T

short_force!(F, short, poses, charges; neighbor_list = nothing)

Ewald3DLong(n_atoms, L; α, s, ϵ = 1.0, ϵ_inf = Inf)

E_l = 1/(2Vϵ) Σ_{0<|k|≤k_c} |ρ_k|² exp(−k²/4α²)/k²,   ρ_k = Σ_j q_j exp(i k·r_j)

Ewald2DLong(n_atoms, L; α, s, ϵ = 1.0)

E_k  = 1/ϵ Σ_i Σ_j q_i q_j cos(k·ρ_ij)
         [ e^{k z_ij} erfc(k/2α + α z_ij) + e^{−k z_ij} erfc(k/2α − α z_ij) ]
         / (8 L_x L_y k)

E_k0 = −1/ϵ Σ_i Σ_j q_i q_j [ e^{−(α z_ij)²}/(α√π) + z_ij erf(α z_ij) ]
         / (4 L_x L_y)

long_energy(long, poses, charges; n_target = long.n_atoms) -> T

long_energy(long::Ewald2DLong, poses, charges; n_target = long.n_atoms) -> T

long_force!(F, long, poses, charges; n_target = long.n_atoms)

long_force!(F, long::Ewald2DLong, poses, charges; n_target = long.n_atoms)

ICM guidance

ICM and its three constructors, ICMEwald2D, ICMEwald3D and ICMPME3D, model a slab confined between two dielectric walls by reflecting each real particle into an image series along z (icm_reflect!, internal).

• γ = (γ_up, γ_down) are the dielectric contrast ratios at the upper and lower walls, conventionally γ = (ϵ_mid - ϵ_out) / (ϵ_mid + ϵ_out) for each interface. γ = 0 is an index-matched (invisible) wall; γ → ±1 approaches a perfect conductor or insulator.

• N_image sets how many reflections deep the image series goes. Each increment adds 2n more image charges (n = number of real particles) and the reflection recurrence shows the series converging geometrically in γ_up * γ_down, so a handful of images is usually enough once |γ_up * γ_down| < 1.

• N_pad, used by ICMEwald3D and ICMPME3D, controls how much the z-period of the padded box is inflated (L_pad[3] = (2·N_pad + 1)·L[3]) before the slab is treated as ordinary triply-periodic Ewald (or particle-mesh Ewald) plus an ELC correction. N_pad must be large enough that periodic replicas of the whole image stack fail to interact — not merely the real slab. The image stack already spans (2·N_image + 1)·L[3] in z, so a N_pad sized only for the bare slab under-pads once N_image is not tiny.

  Measured on this implementation, comparing ICMEwald2D (exact) against ICMEwald3D (Ewald3D + ELC) for the same system:

  N_image	N_pad	disagreement
  3	1	9.4e-3 (too small — not an algorithm error)
  3	2	4.7e-8
  3	3	4.8e-8
  5	2	8.5e-4 (still too small)
  5	3	converges to the same order as N_image = 3, N_pad ≥ 2

  In short: start from N_pad = 2 and re-check against ICMEwald2D (or watch the energy stabilise as N_pad increases) whenever N_image grows past 3.

ICM(long, short, γ, N_image, n_atoms, L; elc = false, N_pad = 0)

ICMShort(n_atoms, L; α, s, ϵ = 1.0, N_image)

ICMEwald2D(n_atoms, L; α, s, γ, N_image, ϵ = 1.0)

ICMEwald3D(n_atoms, L; α, s, γ, N_image, N_pad, ϵ = 1.0)

Error estimates

The following closed-form error estimates are not part of the ExTinyMD source — no code path in this phase calls them — but are reproduced here as sizing guidance, ported from EwaldSummations.jl/src/error_estimate.jl. They are derived in Gao, Zhou, Gan & Liang, "Accurate Error Estimates and Optimal Parameter Selection in Ewald Summation for Dielectrically Confined Coulomb Systems", arXiv:2503.18126.

Lx, Ly are the in-plane box lengths, H the slab half-thickness (or wall separation, depending on convention — see the paper), ϵ the background permittivity, Cq a system-dependent charge-magnitude constant, and gu, gd the two wall contrast ratios (γ_up, γ_down above).

# Truncation error of the image-charge series, after keeping M reflections.
function icm_energy_error(Lx, Ly, H, M, ϵ, Cq, gu, gd)
    return (8 * π^2 / (Lx * Ly * ϵ)) *
           (Cq * (abs(gu * gd))^(floor((M + 1) / 2)) *
            exp(-(4 * π * H * floor((M + 1) / 2)) / max(Lx, Ly))) /
           (1 - abs(gu * gd) * exp(-(4 * π * H) / max(Lx, Ly)))
end

# Error from the ELC slab correction alone (periodic replicas along z),
# for a padded box of height Lz.
function elc_energy_error(Lx, Ly, Lz, H, ϵ, Cq)
    t = 4 * π^2 * Cq / (ϵ * Lx * Ly * (1 - exp(-(2 * π * Lz) / max(Lx, Ly))))
    return t * exp(-(2 * π * (Lz - H)) / max(Lx, Ly)) / (Lz - H)
end

# Combined error when ICM (M reflections) and ELC (padded height Lz) are both in play,
# i.e. the ICMEwald3D route.
function icm_elc_energy_error(Lx, Ly, Lz, H, M, ϵ, Cq, gu, gd)
    t = 4 * π^2 * Cq / (ϵ * Lx * Ly * (1 - exp(-(2 * π * Lz) / max(Lx, Ly))))
    e = t * exp(-(2 * π * (Lz - H)) / max(Lx, Ly)) / abs(Lz - H)
    for l in 1:M
        Cl = abs(gu^(ceil(l / 2)) * gd^(floor(l / 2))) +
             abs(gu^(floor(l / 2)) * gd^(ceil(l / 2)))
        e += t * Cl * exp(-(2 * π * (Lz - (l + 1) * H)) / max(Lx, Ly)) /
             abs(2 * (Lz - (l + 1) * H))
    end
    e += t * 4 * exp(-(2 * π * Lz) / max(Lx, Ly)) /
         (Lz * (1 - exp(-(4 * π * H) / max(Lx, Ly))))
    return e
end

Use these to pick M (N_image) and the padding depth analytically rather than by trial and error; the empirical ICM guidance above is a starting point, not a substitute, for a system whose accuracy requirement is known in advance.

Worked examples

Every figure below is generated when this page is built, so the numbers are whatever the current code actually produces.

A charged system under MD

A molten-salt-like box of equal numbers of +1 and -1 ions, with Lennard-Jones repulsion and Ewald3D electrostatics, integrated with velocity Verlet. Both interactions share one neighbour finder, whose cutoff must cover the larger of the two.

using ExTinyMD, StaticArrays, Random, Plots
Random.seed!(20260915)

n, L = 60, 14.0
boundary = Boundary((L, L, L), (1, 1, 1))

atoms = Atom{Float64}[]
for _ in 1:(n ÷ 2);       push!(atoms, Atom(type = 1, mass = 1.0, charge =  1.0)); end
for _ in (n ÷ 2 + 1):n;   push!(atoms, Atom(type = 2, mass = 1.0, charge = -1.0)); end

info = SimulationInfo(n, atoms, (0.0, L, 0.0, L, 0.0, L), boundary;
                      min_r = 1.0, temp = 1.0)
info.running_step = 1

coulomb = Ewald3D(n, (L, L, L); α = 0.8, s = 3.5)   # r_c = 4.375 < L/2 = 7
lj      = LennardJones(ϵ = 1.0, σ = 1.0, cutoff = 3.0)
finder  = CellList3D(info, max(coulomb.short.r_c, 3.0), boundary, 1)

sys = MDSys(n_atoms = n, atoms = atoms, boundary = boundary,
            interactions = [(lj, finder), (coulomb, finder)],
            loggers = [TemperatureLogger(1000; output = false)],
            simulator = VerletProcess(dt = 1e-4))

steps, E_coulomb = Int[], Float64[]
for block in 0:20
    block > 0 && simulate!(sys.simulator, sys, info, 25)
    push!(steps, block * 25)
    push!(E_coulomb, energy(coulomb, finder, sys, info))
end

plot(steps, E_coulomb; lw = 2, marker = :circle, markersize = 3, legend = false,
     xlabel = "step", ylabel = "electrostatic energy",
     title = "Ewald3D energy during an MD run")

This is the electrostatic component alone, so it is expected to vary as the ions move — it is not a conserved quantity on its own. What it should not do is jump discontinuously or diverge; if it does, the neighbour-finder cutoff is probably smaller than r_c.

The same system without ever building an MDSys:

poses   = [SVector(p.position[1], p.position[2], p.position[3])
           for p in info.particle_info]
charges = [atoms[p.id].charge for p in info.particle_info]

coulomb_energy(coulomb, poses, charges)

-0.04724164014403254

The Ewald split is arbitrary — the total is not

α decides how much work goes to real space and how much to reciprocal space. Move it and the two parts trade off dramatically, while their sum does not move at all. This is the single most useful check on a parameter choice.

Random.seed!(1)
n, L, s = 40, (12.0, 12.0, 12.0), 4.0
poses   = [SVector(rand() * L[1], rand() * L[2], rand() * L[3]) for _ in 1:n]
charges = [isodd(i) ? 1.0 : -1.0 for i in 1:n]

αs = 0.70:0.05:1.20        # r_c = s/α runs 5.71 down to 3.33, all < L/2 = 6
short_part = [short_energy(EwaldShort(n, L; α = α, s = s), poses, charges) for α in αs]
long_part  = [long_energy(Ewald3DLong(n, L; α = α, s = s), poses, charges) for α in αs]
total      = short_part .+ long_part

plot(αs, [short_part long_part total]; lw = 2,
     label = ["real space" "reciprocal space" "total"], legend = :right,
     xlabel = "α", ylabel = "energy", title = "The split moves; the total does not")

Each part swings by roughly a factor of two across this range. The total varies by:

(maximum(total) - minimum(total)) / abs(total[1])

3.3102721187776236e-7

If your own system does not behave this way, α, s and the box are inconsistent — check r_c < min(L)/2 first.

Choosing s

s sets accuracy, with error falling roughly like exp(-s²). Measured against a converged s = 7 reference:

Random.seed!(2)
n, L, α = 40, (12.0, 12.0, 12.0), 1.3
poses   = [SVector(rand() * L[1], rand() * L[2], rand() * L[3]) for _ in 1:n]
charges = [isodd(i) ? 1.0 : -1.0 for i in 1:n]

total_at(s) = short_energy(EwaldShort(n, L; α = α, s = s), poses, charges) +
              long_energy(Ewald3DLong(n, L; α = α, s = s), poses, charges)

reference = total_at(7.0)
svals = 2.0:0.5:5.5
err = [abs(total_at(s) - reference) / abs(reference) for s in svals]

plot(svals, err; lw = 2, marker = :circle, yscale = :log10, legend = false,
     xlabel = "s", ylabel = "relative error", title = "Accuracy versus s")

Thirteen orders of magnitude across a narrow range of s, which is why s between 3 and 5 covers most practical accuracy targets. Remember that raising s at fixed α raises r_c, so it can walk you into the box-size limit.

How deep the image series needs to go

For a dielectrically confined slab, N_image controls how many reflections are kept. The series is geometric, and a stronger contrast γ needs more of it.

Random.seed!(20260927)
n, L = 6, (5.0, 5.0, 4.0)
poses   = [SVector(rand() * L[1], rand() * L[2], 0.8 + 2.4 * rand()) for _ in 1:n]
charges = [isodd(i) ? 1.0 : -1.0 for i in 1:n]

p = plot(; yscale = :log10, xlabel = "N_image", ylabel = "|E(m+1) − E(m)|",
         title = "Convergence of the image series", legend = :topright)
for γv in (0.5, 0.7, 0.9)
    Es = [coulomb_energy(ICMEwald2D(n, L; α = 1.7, s = 4.0,
                                    γ = (γv, γv), N_image = m), poses, charges)
          for m in 1:7]
    plot!(p, 1:6, abs.(diff(Es)); lw = 2, marker = :circle, label = "γ = $γv")
end
p

Each extra reflection buys two to three orders of magnitude here, and the curves flatten once they reach floating-point noise — past that point a larger N_image costs time and buys nothing.