Madelung constants via Ewald summation

The classical Madelung sum for an ionic crystal converges only conditionally in 3D — different truncation shapes give different answers. vibe-qc’s 3D Ewald implementation computes the infinite sum unconditionally.

NaCl rocksalt

import numpy as np
from vibeqc import EwaldOptions, ewald_point_charge_energy

# NaCl primitive cell: 2 ions (Na⁺, Cl⁻) on a face-centered-cubic
# lattice with conventional constant a = 5.64 Å.
A_TO_BOHR = 1.0 / 0.529177210903
a = 5.64 * A_TO_BOHR

# Primitive FCC lattice vectors (columns)
lattice = 0.5 * a * np.array([[0, 1, 1],
                              [1, 0, 1],
                              [1, 1, 0]], dtype=float).T

# One ion of each sign in the primitive cell
positions = np.column_stack([[0, 0, 0],               # Na⁺ at origin
                             [0.5 * a, 0, 0]])         # Cl⁻ at (a/2, 0, 0)
charges = np.array([+1.0, -1.0])

opts = EwaldOptions()
opts.real_cutoff_bohr = 30.0    # α and k-space cutoff auto-derived
energy = ewald_point_charge_energy(lattice, positions, charges, opts)

# Extract the Madelung constant: E = −M / r_nn, where r_nn = a/2 is the
# nearest Na-Cl distance.
r_nn = 0.5 * a
M = -energy * r_nn
print(f"NaCl Madelung constant: {M:.10f}")
print(f"Literature value:        1.7475645946")

Output:

NaCl Madelung constant: 1.7475645946
Literature value:        1.7475645946

Matches to every printed digit.

Why α doesn’t matter

Ewald’s power lies in splitting 1/r into a short-range plus a long-range piece so that both converge fast. The width of the split — the Ewald parameter α — can be tuned to balance work between the two sums, but the total energy must be independent of α. Let’s verify:

for alpha in (0.2, 0.3, 0.5, 1.0):
    opts = EwaldOptions()
    opts.alpha = alpha
    opts.real_cutoff_bohr = 40.0
    opts.recip_cutoff_bohr_inv = 12.0
    E = ewald_point_charge_energy(lattice, positions, charges, opts)
    print(f"α = {alpha:.2f}:  E = {E:.12f}")

All four lines agree to 1e-9 — the hallmark of a correct Ewald implementation.

Dispatching through the SCF driver

For an actual periodic HF/DFT calculation, you don’t need to call ewald_point_charge_energy yourself. Set the Coulomb method on the LatticeSumOptions and the driver handles the rest:

from vibeqc import (
    Atom, PeriodicSystem,
    LatticeSumOptions, CoulombMethod,
    PeriodicSCFOptions,
    # ... run_rhf_periodic etc.
)

opts = PeriodicSCFOptions()
opts.lattice_opts.coulomb_method = CoulombMethod.EWALD_3D
opts.lattice_opts.nuclear_cutoff_bohr = 30.0
# pass opts to run_rhf_periodic(...) as usual

As of the current main, only the nuclear–nuclear lattice sum is fully Ewald-summed; the nuclear-attraction V(g) and the electron repulsion J still use direct truncation. That means bulk 3D total energies are quantitatively correct for the nuclear piece but still cutoff-limited overall. Quantitative bulk SCF lands with Phase 12e-c. See user_guide/ewald.md for the state of the Ewald rollout.