GDF periodic-SCF handover — compcell + AFT state (2026-05-20)

⚠️ TARGET MOVED — v0.9.0 → v0.10.0 (Mike, 2026-05-20). This file keeps its _v0_9_ filename for link stability, but the GDF chain’s release target is now v0.10.0, at the earliest. Rationale: GDF still does not reach PySCF parity on real systems (LiH primitive FCC +22 Ha, MgO 8-atom conv cell ~+1500 Ha — see TL;DR below). It has now been deferred twice (v0.8.0 → v0.9.0 → v0.10.0). There is no schedule pressure on this track — the next GDF chat should treat correctness as the only gate and take the time it needs. Do not let a future release deadline force a half-validated GDF into a tag; the v0.8.0 prep-doc reconciliation (overstated GDF maturity) is the cautionary precedent. The Ewald-3D periodic stack remains vibe-qc’s shipping periodic surface in the meantime; the BIPOLE route is the nearer-term periodic-correlated path. When GDF does reach parity, the release chat picks it up for whatever minor is current then.

This handover documents the state of the GDF route of vibe-qc’s periodic SCF re-implementation (the parallel route to BIPOLE; see docs/handover_bipole_v0_8_2026_05_18.md for that one). It is written so another chat can pick the work up cold.

GDF is the Gaussian-density-fitting periodic SCF stack whose goal is to reproduce PySCF’s KRHF + .density_fit() energies to µHa on real systems. Release target: v0.10.0 at the earliest (deferred from v0.9.0 on 2026-05-20; originally slated for v0.8.0).

TL;DR status

  • H2/STO-3G (Γ-only and multi-k): sub-mHa parity vs PySCF. Done.

  • LiH primitive FCC (2-atom, multi-k kmesh=2,2,2, STO-3G): SCF converges, lands +22 Ha above PySCF at full numerical convergence. Bug localized (see below) — NOT yet at parity.

  • MgO 8-atom conv cell (Γ-only): ~+1500 Ha off PySCF. Separate, larger-signature bug. Not yet diagnosed.

The compcell construction, exxdiv=’ewald’ Madelung shift, multi-k Bloch assembly, AFT correction (now at any q), and pair-FT are all implemented and on origin/main. The remaining gap is a convention bug in the bare 3c lattice sum or chg single-AO FT for L≥1 shells.

What’s complete on origin/main

Drivers

  • vibeqc.run_pbc_gdf_rhf (python/vibeqc/pbc_gdf.py) — Γ-only RHF compcell driver. PBCMethod / PBCExxDiv enums, PBCGDFResult dataclass. Hardening: rejects dim≠3, kmesh≠(1,1,1), charged cells; energy-sanity post-condition.

  • vibeqc.run_krhf_periodic_gdf (python/vibeqc/periodic_k_gdf.py) — multi-k GDF driver. Kwargs: use_compcell, compcell_eta, apply_aft_correction, aft_ft_convention, aft_precision, rcut_strategy, rcut_precision, exxdiv.

Compcell GDF pipeline (python/vibeqc/aux_basis.py)

Sun 2017 compensated-charge construction — fixes the bare-aux divergence on diffuse JKfit auxiliaries:

  1. make_modrho_aux_basis — modrho-renormalised aux.

  2. make_compensating_basis — chg compensating Gaussians (η knob).

  3. make_fused_basis + fuse_transform_matrix — fused basis + A.

  4. build_lpq_compcell — Γ-only compcell Lpq.

  5. build_lpq_bloch_compcell — multi-k compcell Lpq(k) via cell-resolved blocks + Bloch phases.

AFT long-range correction (PySCF _CCGDFBuilder)

  • _compcell_aft_correction — 2c j2c_p subtraction.

  • _compcell_aft_correction_3c — 3c add_ft_j3c subtraction.

  • Both now accept an optional q parameter (commit landed 2026-05-20). For q=0 they keep the historical real-valued Γ-only path bit-exactly; for q≠0 they use the q-shifted kernel 4π/|G+q|² · F̂(G+q) and return complex. This lifted the Γ-only restriction on AFT for multi-k.

exxdiv=’ewald’ (python/vibeqc/madelung.py)

  • madelung_constant_for_cell, apply_exxdiv_ewald_to_K, exxdiv_ewald_energy_shift (closed-shell -¼·ξ·tr(D·S·D·S)).

  • Multi-k K-shift: K(k) += ξ·S(k)·D(k)·S(k) per k.

rcut auto-tuning (python/vibeqc/lattice_screening.py)

  • RcutStrategy enum (FLAT, PYSCF_AUTO), estimate_rcut_pyscf_per_shell, make_lattice_opts.

  • PySCF-style per-shell rcut = √(max(log_prec/(-α), 0.1)).

Tests

tests/test_pbc_gdf_compcell.py — 25 tests, all green. Covers compcell metric plateau, AFT shape/symmetry, q=None≡q=0 parity, q≠0 complex + Hermitian-under-q→-q, H2 sub-mHa (Γ + multi-k), rcut strategies, driver hardening.

The three convention bugs found + fixed this session (2026-05-20)

1. AFT only worked at q=0 → now works at any q

build_lpq_bloch_compcell used to raise NotImplementedError for apply_aft_correction=True at k≠0. Implemented the q-shifted FT kernel. Result on LiH primitive multi-k: bare compcell diverged at +10,582 Ha (not converged); AFT-q converges in 9 iter. Commit: aux_basis: AFT compcell correction at q != 0.

2. Pair-FT (+iG)^t sign bug → fixed to (-iG)^t

_aopair_ft.cartesian_gaussian_product_ft used (+iG_x)^t in the McMurchie-Davidson Hermite-chain FT. The standard identity is F[H_t(r-P;γ)](G) = (-iG_x)^t · … because (∂/∂P_x)^t exp(-iG·P) = (-iG_x)^t exp(-iG·P).

For same-center pairs only t=L_total survives; the odd-L sign flip was being patched by a (-1)^L per-AO factor in _per_ao_pair_libint_to_libcint_scale. For off-center pairs both even-t and odd-t E-coefficient terms contribute and the blanket per-AO sign over-corrected the even-t parts.

Two-part fix: (-iG)^t in cartesian_gaussian_product_ft + dropped the (-1)^L sign in _per_ao_pair_libint_to_libcint_scale (kept only the 1/√(4π/(2L+1)) magnitude — Y_lm normalization).

Verification: examples/debug/gdf_pair_ft_pyscf_diff.py — all 16 LiH STO-3G shell-pair pair-FT blocks now match PySCF’s ft_ao.ft_aopair to 1e-15 (was 1.58 rel-err on the off-center Li 2p × H 1s block). Commit: aopair_ft: fix (+iG)^t → (-iG)^t in pair-FT.

3. LiH residue localized to bare 3c / chg single-AO FT

After fixes 1+2, LiH primitive multi-k still had a residue. A joint precision sweep (3c-rcut × AFT-G-mesh) revealed it is non-monotonic — tightening convergence makes it WORSE:

mode

Δ vs PySCF

iters

bare (no AFT)

+10,582 Ha

50 (✗)

AFT / flat 30-bohr cutoff

+109 Ha

9

AFT / auto rc=1e-8 ε=1e-10

+3.6 Ha

21

AFT / auto rc=1e-10 ε=1e-12

+16.4 Ha

8

AFT / auto rc=1e-12 ε=1e-14

+22.4 Ha

8

The rc=1e-8 “good” number was a coincidental cancellation from under-converged image sums (note the 21-iter struggle). At full convergence the true residue is +22 Ha. Per CLAUDE.md §7 this is a bug, not a convergence-aid problem.

H2 doesn’t exercise this (L=0/L=0 only), which is why H2 sub-mHa parity is preserved.

4. (2026-05-20 cont.) Element-wise AFT diff — FT convention pinned

examples/debug/gdf_aft_pyscf_diff_lih.py (new) compares the compcell-AFT intermediates against a PySCF _CCGDFBuilder subprocess element-wise. Two distinct bugs isolated:

Bug 4a — _aft_fourier_transform_libcint_convention mis-scaled. Probing F̂_fused at G = b0, b1, b0+b1 vs PySCF’s ft_ao:

  • rsgdf_aux_fourier_transform == PySCF ft_ao bit-exactly for all L (ratio +1.00000 at L=0,1,2,3), once the L=1 (py,pz,px) → libcint (px,py,pz) m-permutation is applied. This is the TRUE PySCF / libcint FT convention.

  • _aft_fourier_transform_libcint_convention is mis-named: it does NOT match PySCF’s ft_ao. It is off by exactly √((2L+1)/(4π)) per shell — L=0: 0.28209; L=1: 0.48860; L=2: 0.63078; L=3: 0.74635. It equals rsgdf / √(4π/(2L+1)).

  • Consequence: _compcell_aft_correction with the default ft_convention="libcint" builds j2c_p from the mis-scaled FT → j2c_p comes out ~10× too small (block-norm ratio ≈ 0.10 vs PySCF).

Bug 4b — bare M_fused integral-output normalization mismatch. The bare 2c metric compute_2c_eri_lattice on the fused (modrho_aux ∪ chg) basis does NOT match PySCF’s int2c2e on the same basis (ratio: LiH L=0 diag 0.748, L=1 diag 1.0; system- dependent overall).

IMPORTANT — the fused BASIS is correct. The shell-by-shell parity check (gdf_aft_pyscf_diff_lih.py) confirms: all exponents match PySCF bit-for-bit (0 mismatches / 30 shells), and — decisive — rsgdf_aux_fourier_transform == PySCF ft_ao bit-exactly for all L proves the basis FUNCTIONS χ(r) are identical (the FT uniquely determines the function). The per-shell coefficient numbers differ, but that is only the libint-vs-libcint internal contraction- normalization representation — two encodings of the same χ(r).

So Bug 4b is NOT a basis-construction bug.

Bug 4b — the bare-M_fused mismatch is a lattice-sum cutoff artifact, NOT a compute_2c_eri_lattice bug. This corrects the overstated “structural cross-L bug” framing of commit b43e5b3.

The decisive isolation (gdf_aft_pyscf_diff_lih.py, the MOLECULAR-vs-LATTICE per-(region,L) grids):

  • MOLECULAR 2c metric — vibe-qc compute_2c_eri vs PySCF fused_cell.intor('int2c2e') (single cell, no lattice sum): every (region, L_P, L_Q) block matches at ratio 1.00000. The basic libint 2c call is bit-perfect, all L.

  • LATTICE 2c metric — vibe-qc compute_2c_eri_lattice vs PySCF pbc_intor('int2c2e'): the weird non-factorizable grid ((0,2)=0.04, (1,2)=0.17, (2,3)=2.4, …).

The key realization: the bare fused 2c metric is a conditionally-convergent lattice sum. The modrho_aux block carries net charge — Σ_g (P_0|Q_g) for the aux block does NOT converge; only the charge-neutral (modrho_aux − chg) combination (i.e. the compensated metric, or the AFT-corrected metric) converges. So comparing the raw bare M_fused at vibe-qc’s cutoff_bohr=30 vs PySCF’s own cell.rcut compares two different partial sums of a (near-)divergent series — and different L-blocks converge at different rates, which is exactly the non-factorizable grid seen.

Each individual image term (P_0|Q_g) vibe-qc computes IS correct (it is a molecular integral at a shifted centre, and compute_2c_eri

  • shift_shells are bit-perfect). compute_2c_eri_lattice is not buggy. The bare-metric comparison was simply the wrong thing to compare.

The quantity that must match PySCF is the compensated, AFT- corrected metric A · (M_fused − j2c_p) · Aᵀ, which is the convergent one. With Bug 4a fixed (correct FT scale) the AFT subtraction makes the (cutoff-dependent) bare metric into a (cutoff-independent) compensated metric.

Why these can’t be fixed independently. Bug 4a and Bug 4b currently partially cancel. The H2 “sub-mHa” result uses ft_convention="libint" (the correct rsgdf routine) at the one η=1.0 where the bare-M_fused error happens to cancel. Fixing Bug 4a alone (making “libcint” == “libint” == correct) makes the multi-k default diverge on LiH, because it exposes Bug 4b. The two must be fixed together with full SCF verification at every step.

Next steps (open tasks)

Task #28 — fix Bug 4a; compare the convergent compensated metric

compute_2c_eri_lattice is NOT buggy — the bare-metric mismatch is a conditionally-convergent-lattice-sum cutoff artifact (see Bug 4b above; the MOLECULAR 2c metric matches PySCF bit-exactly for every L-block, ratio 1.0). The actual remaining work:

  1. Fix Bug 4a. _aft_fourier_transform_libcint_convention is off by √((2L+1)/(4π)) per L vs PySCF ft_ao. The proven- correct routine is rsgdf_aux_fourier_transform (== PySCF ft_ao bit-exactly, all L). Route _compcell_aft_correction’s ft_convention="libcint" path through rsgdf_aux_fourier_transform (or add the missing √(4π/(2L+1)) per-shell scale to _aft_fourier_transform_libcint_convention).

  2. Compare the CONVERGENT quantity. Extend gdf_aft_pyscf_diff_lih.py to compare the compensated AFT- corrected metric A·(M_fused − j2c_p)·Aᵀ vs PySCF’s _CCGDFBuilder.get_2c2e() output — that is the cutoff- independent target. The bare M_fused per-L grid is NOT a meaningful target (it is a near-divergent sum; the modrho_aux block carries net charge, only the modrho−chg compensated combination converges).

  3. Verify cutoff-independence of the compensated metric as cutoff_bohr grows (15 → 30 → 45). If it drifts, the AFT mesh (aft_precision) and the lattice cutoff are not matched.

Verify with the H2 sub-mHa tests AND the LiH primitive multi-k bench at every step. Target: LiH primitive multi-k < 1 mHa vs PySCF.

Task #27 — 3c-tensor element-wise PySCF diff (after #28)

Analogous element-wise diff for the bare 3c tensor T[P, μ, ν] (compute_3c_eri_lattice vs PySCF int3c2e_sph). Lower priority than #28 — fix the 2c metric first, then re-measure the 3c residue.

Task #25 — MgO conv-cell ~1500 Ha residue

MgO 8-atom conv cell Γ-only is ~+1500 Ha off PySCF (η-sweep flat across 0.3–2.0). Different signature from LiH (def2-svp-jkfit on Mg/O has higher-L aux shells). May share a root cause with #27 or be a distinct higher-L issue. Diagnose after #27.

Reproduction / environment

The worktree’s editable install points at a _vibeqc_core.so built from the main checkout’s third_party/, so you must set the dyld fallback path:

export DYLD_FALLBACK_LIBRARY_PATH=\
<main-checkout>/third_party/libint/install/lib:\
<main-checkout>/third_party/libxc/install/lib:\
<main-checkout>/third_party/spglib/install/lib:\
<main-checkout>/third_party/fftw/install/lib:\
<main-checkout>/third_party/libecpint/install/lib

PySCF subprocess runner (CLAUDE.md §10 — PySCF is never imported inside python/vibeqc/; parity scripts spawn it):

export VIBEQC_PYSCF_PYTHON=<a-venv-with-pyscf>/bin/python

Key debug scripts under examples/debug/:

  • gdf_pair_ft_pyscf_diff.py — pair-FT element-wise diff (the diagnostic that cracked bug 2). Now bit-parity.

  • gdf_lih_primitive_multik_parity.py — LiH primitive multi-k bench with the joint precision sweep.

  • gdf_aft_calibration.py — L>0 single-AO FT convention.

  • gdf_bare_aux_divergence.py — bare-aux divergence + compcell plateau demo.

Run the compcell test suite:

cd <worktree> && python -m pytest tests/test_pbc_gdf_compcell.py -q

Discipline notes (CLAUDE.md)

  • §7 — periodic SCF landing at impossible energies is a bug; don’t paper over with damping / level-shift / empirical thresholds. The +22 Ha LiH residue is being chased as a bug.

  • §10 — never import pyscf inside python/vibeqc/. All parity validation is out-of-process subprocess runners.

  • §2 — rebase-then-push; never force-push main.