Ewald summation in periodic systems

In 3D periodic solids the Coulomb series \(\sum_{g,A,B} Z_A Z_B / |R_{AB} + g|\) is conditionally convergent: different truncation shapes give different answers. The standard fix is Ewald summation, which splits each \(1/r\) term into a short-range part (evaluated in real space, where it converges exponentially) plus a long-range part (evaluated in reciprocal space, where the Gaussian damping makes the sum converge fast too). vibe-qc’s full 3D Ewald infrastructure — nuclear repulsion, nuclear attraction \(V(g)\), and the electron-repulsion Coulomb \(J(\mu\nu)\) — is shipped through a single Coulomb-method dispatcher; this page tells you how to use it.

The user-facing recommendation

For any 3D-periodic SCF, set the Coulomb method to EWALD_3D and call the dispatcher:

import vibeqc as vq

opts = vq.PeriodicSCFOptions()
opts.lattice_opts.coulomb_method = vq.CoulombMethod.EWALD_3D

# Multi-k:
result = vq.run_rhf_periodic_scf(system, basis, kmesh, opts,
                                 omega=0.5, spacing_bohr=0.3)
# Γ-only:
result = vq.run_rhf_periodic_gamma_scf(system, basis, opts,
                                       omega=0.5, spacing_bohr=0.3)

The dispatcher routes on opts.lattice_opts.coulomb_method to the appropriate backend:

CoulombMethod

Routed backend (multi-k)

Routed backend (Γ-only)

Status

DIRECT_TRUNCATED

run_rhf_periodic

run_rhf_periodic_gamma

Default. Half-summed direct truncation — exact in the molecular-limit regime (large unit cell, cutoff excludes all non-zero cells). Conditionally convergent in 3D bulk.

EWALD_3D

run_rhf_periodic_multi_k_ewald3d

run_rhf_periodic_gamma_ewald3d

Full Gaussian-charge Ewald — unconditionally convergent. Use for any quantitative 3D bulk SCF. Requires orthorhombic cells (FFT-based long-range solver).

SLAB_EWALD_2D

—

—

Not yet wired into SCF dispatch (raises NotImplementedError).

NEUTRALIZED_1D

—

—

Not yet wired into SCF dispatch.

Nuclear repulsion routes through the same enum: when EWALD_3D is selected, nuclear_repulsion_per_cell automatically uses Ewald summation, so the electronic and nuclear contributions stay self-consistent.

Same dispatcher for KS and UKS

The pattern above is the HF dispatcher. The KS (Kohn-Sham DFT) and UKS (open-shell DFT) families ship the same dispatcher shape and the same CoulombMethod enum:

Driver family

Dispatcher (multi-k)

Dispatcher (Γ-only)

Phase

Closed-shell HF

run_rhf_periodic_scf

run_rhf_periodic_gamma_scf

12e-c-4

Closed-shell KS

run_rks_periodic_scf

run_rks_periodic_gamma_scf

15c-1, 15c-2

Open-shell UKS

run_uks_periodic_scf

run_uks_periodic_gamma_scf

15c-3

Each accepts the same keyword arguments (omega, grid_shape, origin, spacing_bohr, linear_dep_threshold) and routes on options.lattice_opts.coulomb_method identically. Pick the options class that matches your method (PeriodicSCFOptions / PeriodicRHFOptions for HF, PeriodicKSOptions for RKS/UKS), set the Coulomb method, and call the dispatcher.

The Ewald-specific keyword arguments to the dispatcher (omega, grid_shape, origin, spacing_bohr, linear_dep_threshold) are forwarded to the EWALD_3D backend and ignored on the DIRECT_TRUNCATED path. Defaults: omega=0.5 bohr⁻¹ (the short/long-range splitting parameter), spacing_bohr=0.3 (reciprocal-space FFT grid spacing). Most users never need to touch these — the only knob that matters in practice is the CoulombMethod choice itself.

A complete worked example

import vibeqc as vq

# Conventional cubic LiH rocksalt — 8 atoms per cell, lattice 4.084 Å.
a_ang = 4.084
a = a_ang / 0.529177210903
import numpy as np
lat = a * np.eye(3)

# Off-FFT-grid shift: atoms exactly at the corner of the box inflate
# the long-range residual ~100x; a small (0.05 bohr) shift fixes it.
shift = 0.05
unit_cell = []
for z in (3, 1):                                       # 3 = Li, 1 = H
    base = 0.5 if z == 1 else 0.0
    for fx, fy, fz in [(0, 0, 0), (0, 0.5, 0.5),
                        (0.5, 0, 0.5), (0.5, 0.5, 0)]:
        x = (fx + base) * a + shift
        y = (fy + base) * a + shift
        z_ = (fz + base) * a + shift
        unit_cell.append(vq.Atom(z, [x, y, z_]))

system = vq.PeriodicSystem(dim=3, lattice=lat, unit_cell=unit_cell)
basis = vq.BasisSet(system.unit_cell_molecule(), "sto-3g")

# Γ-only SCF with Ewald
opts = vq.PeriodicSCFOptions()
opts.lattice_opts.coulomb_method = vq.CoulombMethod.EWALD_3D
opts.lattice_opts.cutoff_bohr = 12.0          # short-range AO cutoff
opts.lattice_opts.nuclear_cutoff_bohr = 25.0   # nuclear-Ewald real cutoff

result = vq.run_rhf_periodic_gamma_scf(
    system, basis, opts,
    omega=0.5,
    spacing_bohr=0.4,
)

print(f"Converged:    {result.converged}")
print(f"Total energy: {result.energy:.8f} Ha per cell")
print(f"SCF iters:    {result.n_iter}")

For multi-k, build a :class:vibeqc.KPoints (or, for back-compat, a raw BlochKMesh) and call run_rhf_periodic_scf(system, basis, kpoints, opts). The result mirrors the Γ-only structure plus per-k orbital coefficients. KPoints exposes Monkhorst–Pack, Γ-centered, IBZ-reduced (via spglib), HPKOT band paths (via seekpath), explicit user lists, and density-based auto-meshes (KPPRA / kspacing / VASP Auto) — see user_guide/k_points.md for the full reference.

What “ω-invariance” means and why we test it

Ewald has one free parameter — the splitting \(\omega\) that decides how much of \(1/r\) is sent to real-space vs reciprocal-space. The total energy must be independent of \(\omega\) to whatever the real- and reciprocal-space cutoffs allow. This is the primary correctness test of any Ewald implementation. vibe-qc’s bulk-crystal benchmark suite verifies it explicitly across H₂ / LiH / MgO / Ne fixtures: scanning \(\omega \in \{0.3, 0.4, 0.5, 0.6, 0.7\}\) recovers the same total energy to within \(10^{-8}\) Ha when the cutoffs are adequate. If your run is \(\omega\)-sensitive at the \(10^{-4}\) level, the cutoffs are too tight — increase cutoff_bohr and nuclear_cutoff_bohr and re-test.

Choosing cutoffs and the FFT grid

The EWALD_3D path has three convergence knobs:

  • opts.lattice_opts.cutoff_bohr — real-space cutoff for AO pair / ERI integrals. Default 15 bohr; can usually drop to 10–12 for valence-only bases like sto-3g, pob-DZVP.

  • opts.lattice_opts.nuclear_cutoff_bohr — real-space cutoff for the nuclear-Ewald real-space term. Default 25 bohr; auto-derived from omega if you pass omega explicitly.

  • spacing_bohr= (kwarg on the dispatcher) — reciprocal-space FFT grid spacing. Default 0.3 bohr — gives a \(\approx 64^3\) grid for a 20-bohr cubic cell, which converges the long-range J piece to \(\sim 10^{-8}\) Ha. Tighten to 0.2 for benchmark-quality runs.

The auto-derived \(\alpha\) (Ewald splitting parameter) is chosen so that erfc\((\alpha \cdot R_{\text{cut}}) \approx 10^{-12}\). That makes both the real- and reciprocal-space sums converge fast on the supplied cutoffs without you having to tune \(\alpha\) by hand.

Low-level Ewald primitives

The dispatcher is the right entry point for SCF work. For quasi-classical lattice-electrostatics calculations (Madelung constants, point-charge potentials, charge-only models), use the primitives directly:

from vibeqc import EwaldOptions, ewald_point_charge_energy
import numpy as np

# NaCl-like rocksalt with ±1 charges
a = 5.6 / 0.529177210903                       # 5.6 Å in bohr
lattice = 0.5 * a * np.array([[0, 1, 1],
                              [1, 0, 1],
                              [1, 1, 0]], dtype=float).T
positions = np.column_stack([[0, 0, 0],
                             [0.5 * a, 0, 0]])
charges = np.array([+1.0, -1.0])

opts = EwaldOptions()
opts.real_cutoff_bohr = 30.0          # auto-derive α and recip cutoff
energy = ewald_point_charge_energy(lattice, positions, charges, opts)
print(f"NaCl Madelung / pair: {energy:.8f} Ha   (expected -0.330275485)")

Other low-level entry points:

  • vq.ewald_point_charge_potential — Madelung potential at given probe points.

  • vq.ewald_nuclear_potential / vq.ewald_nuclear_repulsion — nuclei-only versions.

  • vq.compute_nuclear_lattice_ewald — nuclear-attraction matrix elements via Ewald.

  • vq.build_j_ewald_3d — direct J(μν) Coulomb build at a single density (used inside the SCF dispatcher).

These all operate on the same EwaldOptions / LatticeSumOptions machinery; the dispatcher just wires them together at SCF time.

Theory

Splitting the Coulomb sum

For a 3D-periodic distribution of point charges \(\{Z_A\}\) at positions \(\{R_A\}\) in a unit cell with translations \(\{g\}\), the energy

\[ E_{\text{nn}} = \tfrac{1}{2} \sum_{A, B} {\sum_{g}}' \frac{Z_A Z_B}{|R_A - R_B + g|} \]

(primed sum excludes the \(A = B, g = 0\) term) is conditionally convergent. Ewald’s 1921 trick is to split each \(1/r\) as

\[ \frac{1}{r} = \frac{\operatorname{erfc}(\alpha r)}{r} + \frac{\operatorname{erf}(\alpha r)}{r} \]

The first term decays exponentially in real space; the second is smooth enough to Fourier-transform and sum in reciprocal space. The total energy under the tin-foil + jellium convention becomes

\[ E_{\text{nn}} = E_{\text{real}} + E_{\text{recip}} + E_{\text{self}} + E_{\text{bg}} \]

with

  • \(E_{\text{real}} = \tfrac{1}{2} \sum_{A,B} {\sum_{g}}' Z_A Z_B \cdot \operatorname{erfc}(\alpha \, r) / r\) — exponentially convergent real-space sum.

  • \(E_{\text{recip}} = \tfrac{2\pi}{V} \sum_{G \neq 0} |S(G)|^2 \, e^{-|G|^2/(4\alpha^2)} / |G|^2\) — k-space sum with Gaussian damping. \(S(G) = \sum_A Z_A \, e^{iG \cdot R_A}\) is the nuclear structure factor.

  • \(E_{\text{self}} = -(\alpha/\sqrt{\pi}) \sum_A Z_A^2\) — removes the spurious self-interaction from each nucleus’s Gaussian-smeared image.

  • \(E_{\text{bg}} = -(\pi/(2 \alpha^2 V)) (\sum_A Z_A)^2\) — background correction for a non-neutral cell (the “jellium” prescription).

From point charges to AO densities

For an AO basis the analogue construction extends the same idea to Gaussian charge distributions. Each AO pair \(\chi_\mu \chi_\nu^g\) is itself a sum of Gaussians; their \(1/r\)-mediated interaction with the rest of the lattice splits the same way:

  • Near field (AO pairs within a short-range cutoff): exact analytic Gaussian-Gaussian Coulomb integrals, no multipole expansion.

  • Far field: the AO pair is approximated by its multipole expansion (monopole + dipole + quadrupole) and the multipoles are Ewald-summed in reciprocal space via FFT.

This is the Saunders–Dovesi multipolar splitting used in CRYSTAL. The FFT step is the reason EWALD_3D requires an orthorhombic cell — vibe-qc’s FFT-based long-range Poisson solver doesn’t yet support triclinic lattices. (A triclinic-FFT update is on the roadmap; it’s a Cauchy-Born transform on the existing solver, not new physics.)

Nuclear-attraction \(V(g)\)

The same erfc / erf splitting applies to the nuclei-electron attraction matrix elements

\[ V_{\mu\nu}(g) = \sum_{A, h} -Z_A \, \langle \chi_\mu | 1 / r_{A,h} | \chi_\nu^g \rangle \]

The short-range piece is computed by compute_nuclear_erfc_lattice (real-space, exponentially convergent in \(h\)). The long-range piece adds the smooth erf\((\omega r)/r\) kernel, which integrates exactly over Gaussian AO pairs in reciprocal space. The dispatcher wires both together inside the SCF Fock build.

References

  • Ewald, P. P. “Die Berechnung optischer und elektrostatischer Gitterpotentiale,” Ann. Phys. 64, 253 (1921). The original Ewald paper.

  • Saunders, V. R.; Freyria-Fava, C.; Dovesi, R.; Salasco, L.; Roetti, C. “On the electrostatic potential in crystalline systems where the charge density is expanded in Gaussian functions,” Mol. Phys. 77, 629 (1992). Multipolar splitting used in CRYSTAL — the basis of vibe-qc’s far-field treatment.

  • Toukmaji, A. Y.; Board Jr., J. A. “Ewald summation techniques in perspective: a survey,” Comput. Phys. Commun. 95, 73 (1996). Modern Ewald review covering the tin-foil convention, \(\alpha\) optimization, and the various flavors used in MD codes.

  • Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids, 2nd ed., Oxford University Press (2017). Clear textbook treatment with the \(\alpha\) auto-tuning formula vibe-qc uses.

  • Pisani, C.; Dovesi, R.; Roetti, C. Hartree-Fock Ab Initio Treatment of Crystalline Systems. Lecture Notes in Chemistry vol. 48, Springer (1988). Foundational text for periodic LCAO HF with the Ewald/multipolar treatment of the lattice Coulomb sums CRYSTAL is built on.