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 periodic Coulomb treatment is being built up in three sub-phases; today’s main branch has the first two shipped.

Dispatch — `CoulombMethod`¶

Every LatticeSumOptions has a coulomb_method field that selects how 1/r lattice sums are evaluated:

from vibeqc import LatticeSumOptions, CoulombMethod

opts = LatticeSumOptions()
opts.coulomb_method = CoulombMethod.DIRECT_TRUNCATED   # default
# or:
opts.coulomb_method = CoulombMethod.EWALD_3D           # quantitative 3D

Enum value	Meaning	Status
`DIRECT_TRUNCATED`	Half-summed direct truncation. Default. Exact in the molecular-limit regime (large unit cell, cutoff excludes all non-zero cells). Conditionally convergent in 3D bulk.	Works for all sums (molecular-limit correctness)
`EWALD_3D`	Ewald summation under the tin-foil + jellium convention. Unconditionally convergent.	Nuclear repulsion: ✅ shipped (Phase 12e-a) Nuclear attraction V(g): k-space part pending (12e-c) ERI Coulomb J: pending (12e-c)
`SLAB_EWALD_2D`	2D slab Ewald.	Not yet implemented
`NEUTRALIZED_1D`	1D line with neutralising background.	Not yet implemented

Phase 12e-a: nuclear-nuclear Ewald (shipped)¶

The Madelung 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) becomes unconditionally convergent under Ewald summation with the tin-foil boundary convention:

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

where:

\(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).

The Ewald parameter \(\alpha\) can be set explicitly or auto-derived from a real-space cutoff so that erfc\((\alpha R_{\text{cut}}) \approx\) tolerance.

Usage (low-level)¶

For an arbitrary 3D lattice of point charges:

import numpy as np
from vibeqc import EwaldOptions, ewald_point_charge_energy

# 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)")

Usage (through `nuclear_repulsion_per_cell`)¶

For a real periodic SCF, route through LatticeSumOptions.coulomb_method:

import numpy as np
from vibeqc import (
    Atom, PeriodicSystem, LatticeSumOptions, CoulombMethod,
    nuclear_repulsion_per_cell,
)

a = 5.6 / 0.529177210903
lattice = 0.5 * a * np.array([[0, 1, 1],
                              [1, 0, 1],
                              [1, 1, 0]], dtype=float).T
sysp = PeriodicSystem(
    dim=3, lattice=lattice,
    unit_cell=[Atom(11, [0, 0, 0]),
               Atom(17, [0.5 * a, 0, 0])],
)

opts = LatticeSumOptions()
opts.coulomb_method = CoulombMethod.EWALD_3D
opts.nuclear_cutoff_bohr = 30.0
print(nuclear_repulsion_per_cell(sysp, opts))

Choosing cutoffs and \(\alpha\)¶

If you pass only real_cutoff_bohr, \(\alpha\) and the reciprocal cutoff are chosen automatically to achieve tolerance (default 1e-12) in both sums.
If you pass an explicit alpha, use a larger real cutoff or the real-space sum may not converge.
Final energy is independent of \(\alpha\) to the specified tolerance — this is the primary correctness test of any Ewald implementation. vibe-qc’s test suite verifies it explicitly for NaCl, CsCl, ZnS, and simple-cubic jellium.

Phase 12e-b: erfc-screened nuclear attraction (building block)¶

compute_nuclear_erfc_lattice(basis, system, omega, options) returns the short-range component of \(V(g)\):

\[ V^{\text{erfc}}_{\mu\nu}(g; \omega) = \sum_{A, h} -Z_A \, \langle \chi_\mu | \operatorname{erfc}(\omega \, r_{A,h}) / r_{A,h} | \chi_\nu^g \rangle \]

At \(\omega \to 0\) this reduces to the ordinary 1/r nuclear attraction (compute_nuclear_lattice); at finite \(\omega > 0\) it is exponentially convergent as a real-space lattice sum. This is the plumbing 12e-c will combine with a reciprocal-space long-range piece to provide the full \(\text{EWALD\_3D}\) treatment of \(V(g)\).

from vibeqc import compute_nuclear_erfc_lattice, LatticeSumOptions

opts = LatticeSumOptions()
opts.cutoff_bohr = 15.0
opts.nuclear_cutoff_bohr = 25.0
V_erfc = compute_nuclear_erfc_lattice(basis, system, omega=0.3, opts)

Today this function is primarily useful for tests and for users exploring Ewald internals — the SCF drivers do not call it yet. Until 12e-c lands, the default path (direct-truncated V) remains correct in the molecular limit and accurate-enough-for-qualitative in genuine 3D bulk; quantitative 3D results arrive with 12e-c.

Phase 12e-c: full Gaussian-charge Ewald (pending)¶

The long-range parts of both the nuclear attraction \(V(g)\) and the electron-repulsion integrals \((\mu\nu | \lambda\sigma)\) are smooth Gaussian-smeared potentials; their matrix elements fold as reciprocal- space sums involving the Fourier transform of Gaussian AO pairs. The planned implementation follows the standard Saunders–Dovesi multipolar splitting used in CRYSTAL:

Near field (nuclei / AO pairs within a short-range cutoff of each other): compute exact integrals, no multipole expansion.
Far field: expand the charge distribution in a multipole series (monopole + dipole + quadrupole typically sufficient) and Ewald-sum the multipoles in reciprocal space. The long-range tail converges exponentially in \(|G|\).

This phase also adds FFTW3 as a build-time dependency to accelerate the reciprocal-space grid sums. Until it ships, 3D bulk calculations with EWALD_3D are limited to the nuclear-repulsion piece — sufficient to validate Madelung constants and correctly compute total energies in the molecular-limit regime, but not yet to give quantitative bulk SCF energies.

References¶

P. P. Ewald, Ann. Phys. 64, 253 (1921). Original Ewald paper.
V. R. Saunders, C. Freyria-Fava, R. Dovesi, L. Salasco, C. Roetti, Mol. Phys. 77, 629 (1992). Multipolar splitting used in CRYSTAL.
A. Y. Toukmaji, J. A. Board Jr., Comp. Phys. Commun. 95, 73 (1996). Modern Ewald review.
M. P. Allen, D. J. Tildesley, Computer Simulation of Liquids, 2nd ed., Oxford, 2017. Clear textbook treatment with the \(\alpha\) auto-tuning formula vibe-qc uses.