Design: self-consistent long-range gamma for GFN2-xTB-SECCM

Status: implemented (2026-08-18) as the opt-in ewald_gamma=True option, shipped behind the default-off flag in v0.15.137.dev. The construction below supersedes the earlier proposal in this file’s git history: the original erf/erfc Ewald split was numerically evaluated and rejected (see “The rejected split”, below), and the tblite get_amat_3d convention was adopted as the convergent reference.

Problem

The GFN2-xTB-SECCM electrostatics kernel is assembled in two pieces:

  1. the Wigner-Seitz folded Klopman-Ohno shell gamma (short range, partition-of-unity weights), and

  2. the opt-in Madelung/Ewald tail of the atomic Mulliken fluctuations (SMADEL subtraction of the WS-cell bare Coulomb), shipped in v0.15.136 as madelung=True.

The two pieces are inconsistent at the WS boundary: inside the cell the kernel is the damped Klopman-Ohno form gamma(R) = 1/sqrt(R^2 + eta^2) (eta = shell-pair hardness), outside it is the bare Coulomb 1/R. The error is small in the energy but large in the charge response, because the damping affects the derivative of the kernel. The observed failures share one signature: a strongly negative charge-response eigenvalue of the WS-truncated supercell map:

  • fcc Cu 2x2x2 at compressed volumes and 4x4x4 at all volumes: the over-polarized shell-transfer basin dominates the map (IID 130);

  • corundum with the embedding: the ionic mode runs away;

  • the historical aligned Mg/O-plane stress fixture with 4+ layers: period-2/period-3 limit cycles, with and without the embedding. This fixture is an artificial two-coordinate P4/mmm stack, not B1 rocksalt(100); physical B1 study geometries are maintained separately.

No mixer (simple, charge-DIIS, bounded Aitken) and no deposit convention (diagonal vs S-weighted) removes the instability: the map itself carries it. The xtb reference converges every one of these cells with the same GFN2 parameters, so the defect is in the kernel assembly, not in the method.

Why a naive lattice sum of the KO kernel does not exist

The first candidate fix was an Ewald split of the Klopman-Ohno kernel,

Gamma(si, sj) = sum_n gamma_eta(|R(si,sj) + n T|)

summed over all supercell images. Numerical evaluation (studies/seccm-bulk3d/probe_kernel_split.py, probe_kernel_element.py) shows this object is IR-ill-defined for multi-element cells:

  • gamma_eta(R) - 1/R decays as -eta^2/(2 R^3), and a 3-D lattice sum of 1/R^3 diverges logarithmically. With pair-dependent eta = 2/(h_a + h_b), the charge-weighted leading image term sum_ab q_a q_b eta_ab^2 does not cancel, so both the per-element matrix and the neutral-cell energy diverge. Measured on MgO 2x2x2: 1/2 dq . Gamma . dq drifts -0.03 -> -1.75 -> -3.24 Ha as the real-space cutoff grows 25 -> 40 -> 60 bohr.

  • The erfc/erf split proposed in the original version of this file double-counts the erf(alpha r)/r piece (once in the real-space bracket, once in the reciprocal series), shifting the kernel by the conditionally convergent sum of (2 erf - 1)/r over the images (about -47 mHa on the MgO 2x2x2 energy; probe_kernel_split.py).

  • A pure-reciprocal representation with the exact Fourier transform F[gamma_eta](G) = 4 pi eta K1(G eta)/G converges (the finite eta is its own short-range regulator), but dropping G=0 subtracts a pair-dependent divergent constant, so the result is a different, equally arbitrary regularization of the same divergence.

The physical reason: the KO kernel is an approximation to the exact density-based Coulomb integral, whose gamma - 1/R decays exponentially. The KO power-law tail R^-3 makes the periodic model IR-divergent; every implementation must choose a regularization. The xtb binary is cutoff-stable to 1e-9 across its 20-28 bohr cutoffs, so tblite has a convergent convention - and it is the reference we match.

Implemented construction (tblite get_amat_3d)

Reading tblite/coulomb/charge/effective.f90 (the module behind the xtb binary’s GFN2 electrostatics) gives the convergent convention:

Gamma(si,sj) = Ewald_Coul(d)
             + sum_WS w * [gamma_eta(R_img) - 1/R_img]
             + delta_same_atom * gamma_eta(0)

with the standard full Ewald Coulomb kernel (no background term; neutral cells only) and gamma_eta(0) = 1/eta_ab the molecular hardness. The Wigner-Seitz weights truncate the KO-minus-Coulomb correction inside the WS cell (finite by construction, the CCM contract), while the Coulomb tail rides the converged Ewald series. The Ewald part is dimension-specific:

  • 3-D: `sum_n erfc(alpha r_n)/r_n

    • (4 pi / V) sum_{G!=0} e^{-G^2/4 alpha^2}/G^2 cos(G.d)

    • delta(d=0) 2 alpha/sqrt(pi)with the volume-adaptivealpha = sqrt(pi)/cbrt(V)(theindo::_madkonst_3d` convention);

  • 2-D: the Parry/Heyes Ewald (the validated indo::_madkonst_2d arithmetic): the erfc real-space sum, the reciprocal series `(pi/A) sum_K cos(K.d)/K [e^{K z} erfc(alpha z + K/2 alpha)

    • e^{-K z} erfc(-alpha z + K/2 alpha)], the K=0 term -(2 pi/A)[z erf(alpha z) + e^{-alpha^2 z^2}/(alpha sqrt(pi))], and the self term, with alpha = 0.85 sqrt(pi/A)`;

  • 1-D: not implemented yet (the wire KO kernel is the next increment).

Properties:

  • the eta-to-zero off-diagonal limit is the bare-Coulomb Ewald kernel (the validated Madelung machinery) - the regression anchor;

  • single-element cells use the same finite WS regularization. They do not generally define a full KO lattice sum: GFN2 hardness is shell-dependent, and the KO-minus-Coulomb correction retains its non-absolutely-convergent power-law tail;

  • multi-element cells: the WS weights are the regularization, exactly as in tblite - the xtb parity reference;

  • H0, S, and the pair repulsion keep their WS-folded assembly (they are short-ranged and that is the CCM contract), while the charge functional (the gamma channel) uses this kernel, exactly the role the MSINDO CCM reserves for its Madelung embedding. The AES channel keeps its WS-weighted assembly for this milestone.

Validation (2026-08-18)

  • Kernel-vs-formula: Python reimplementations of both constructions reproduce the C++ engine’s cyclic_scc_energy at the converged shell charges, in 3-D (MgO 2x2x2) and in 2-D (the two-layer synthetic aligned-plane stress cell) (tests/test_ccm_semiempirical.py::test_ewald_gamma_*).

  • MgO 2x2x2: converges (no limit cycles), Mulliken q_rms 0.40-0.44 vs xtb 0.39-0.41 across a = 4.0-4.5 A; the energy tracks xtb’s curve shape with a constant offset of about -88 mHa/atom (the WS-folded H0/AES/repulsion channels, not the kernel). xtb’s own Gamma-only 2x2x2 curve is monotone on this range, so the SECCM curve shape is the parity signal, not a minimum.

  • Corundum hex cell: converges in ~180 iterations (was runaway), q_rms 0.369 vs xtb 0.41.

  • 2-D: the two-layer synthetic aligned-plane stress cell converges through the Parry/Heyes kernel (115 iterations, gap 0.96 eV, q_rms 0.356) and matches the validated atomic-embedding (madelung=True) construction on the same cell to 2e-3 Ha (energy) / 2e-2 e (charges). This agreement cross-validates the two tail constructions but is not a rocksalt surface result. The four-layer and six-layer versions cycle instead of converging under every tested linear mixer. The cycle mechanism is now pinned (studies/seccm-bulk3d/probe_slab_jacobian.py): a Python reconstruction of the SCC map (validated bit-identical against the native iteration, 5e-15 on the 2-layer trace) shows the linear map at the neutral point is contractive for both cells (J min eigenvalue -0.99 / -1.34), while trial iterations near the four-layer branch have a strongly repelling response. The measured T=0 trajectory develops a near-degenerate frontier and a period-2/3 cycle. These measurements show numerical nonconvergence of the tested iteration, not nonexistence of an exact fixed point.

  • Newton SCC mixer (scc_mixer="newton", 2026-08-19): a line-searched frozen-lagged-multipole chord/quasi-Newton step is wired into the GFN2-SECCM loop (factored charge map, finite-difference charge response with lagged multipoles held fixed, constrained Newton solve, trust-region cap, and residual-reduction line search). On the production T=0 map, every tested solver (the chord step, warm starts, Levenberg-Marquardt, and a GAM3 homotopy) plateaus at a roughly 0.24 e residual when the near-degenerate Aufbau frontier switches branches. A small Fermi temperature smooths that measured discontinuity and scc_mixer="newton" plus electronic_temperature=0.002 converges the four-layer synthetic cell in 6 iterations (energy -9.11722295 Ha/prim cell, q_rms 0.429). The two-layer cell converges under the same mixer at T=0 to the simple-mixing state (energy difference 3.8e-10, 115 to 8 iterations). None of these numerical probes proves that an exact T=0 fixed point does or does not exist.

  • fcc Cu 2x2x2: the new map reaches a charge-density-wave basin (atomic q_rms 2.7 e) whose signature the shipped shell-rms gate did not cover; the gate now also fails atomic q_rms > 2.0 e closed, and the Cu runs fail closed. Reaching the sane basin with this kernel is the open follow-up (IID 130).

  • 1-D cells fail closed (the wire KO kernel is not implemented yet).

  • Molecular limit: bit-identical with the option off (unchanged).

Escalation questions for the maintainer

  1. Option surface: ewald_gamma ships as a separate flag next to madelung; should the two merge into one graded option later?

  2. Scope: keep the AES WS-weighted for this milestone (done), or fold its tail into the same treatment?

  3. The released v0.15.136 wording: the changelog describes the embedding as the Madelung tail; the new kernel replaces that mechanism. Acceptable as a dev-cycle change, or gate behind the option default until a release note can be written?