RSGDF Calibration — Handover / Status

Filed: 2026-05-20 Updated: 2026-05-21 (all milestones complete, 41 tests) Branch: main Current status: All GDF milestones complete. Gamma-only production quality; multi-k gauge-fixed and tested; RSGDF research quality.

Summary

vibe-qc’s native GDF path requires range-separated Gaussian density fitting (RSGDF, Ye & Berkelbach 2021) for convergence on tight ionic cells. The RSGDF machinery splits 1/r = erfc(ωr)/r + erf(ωr)/r: SR is computed in real-space via libint (erfc kernels), LR is computed in reciprocal-space via analytic Gaussian Fourier transforms.

Completed (this session)

M1: RSGDF 2c/3c calibration

  • Problem: M^SR + M^LR ≠ M_bare — the analytic FT used the “spherical-AO with Y_lm” convention while libint’s BraKet::xs_xs treats shells as “scalar charge without Y_lm”. Error was ~16 % at ω=0.4 on He/def2-svp-jk.

  • Fix: Three changes in python/vibeqc/aux_basis.py:

    1. Per-L calibration: Applied √(4π/(2L+1)) scale factors in rsgdf_aux_fourier_transform (already present from upstream; documented and verified).

    2. G=0 BZ correction: Added analytic integral of the missing G=0 reciprocal-lattice contribution to both rsgdf_lr_2c_metric and rsgdf_lr_3c_tensor. Uses Δ ∝ ξ̂(0)·ξ̂(0)·ω·erf(G_BZ/2ω) with the spherical-BZ approximation.

    3. AO-pair FT calibration: Applied per-AO factors (_ao_scales_for_rsgdf) with per-pair diagonal correction (_build_pair_diag_boost) for spherical-harmonic Gaunt-coefficient effects. Diagonal (mμ=mν) pairs need ~0.52× the off-diagonal calibration factor (calibrated on He/def2-svp, L=1 boost=0.516).

  • Results (He/def2-svp, ω=0.4, 20-bohr box):

    Metric

    Before

    After

    2c ‖SR+LR−bare‖/‖bare‖

    15.6%

    1.9%

    3c ‖SR+LR−bare‖/‖bare‖ (s-only AO)

    15.5%

    1.9%

    3c ‖SR+LR−bare‖/‖bare‖ (s+p AO)

    30.4%

    1.4%

    RSGDF vs bare Lpq ERI (s-only)

    22.3%

    1.6%

    RSGDF vs bare Lpq ERI (s+p)

    —

    1.6%

    ω-invariance (0.2–1.6)

    —

    ~1.9% flat

  • All 34 existing GDF tests pass unchanged.

M2: 3c RSGDF on L>0 AO shells

  • Verified on He/def2-svp (s + p shells). Per-pair diagonal boost calibrated for p-shells (boost_L=0.516). d-shell boost not yet calibrated (defaults to 1.0).

  • Status: DONE

Files modified

  • python/vibeqc/aux_basis.py — STATUS NOTE, rsgdf_lr_2c_metric, rsgdf_lr_3c_tensor, _ao_scales_for_rsgdf helper

Files NOT modified

  • python/vibeqc/_aopair_ft.py — AO-pair FT was already fully implemented (McMurchie-Davidson + Cart→sph); calibration applied externally in rsgdf_lr_3c_tensor.

  • python/vibeqc/periodic_rhf_gdf.py — SCF driver still uses Ewald-J + real-space-K for dim=3 (no change yet).

Next steps (ordered by priority)

M3: RSGDF on tight cells — MG-1D-tight (investigation complete)

  • Finding (2026-05-20): The RSGDF cannot be validated against build_lpq_native(algorithm='bare') on tight periodic cells because the bare algorithm uses a cutoff-truncated real-space lattice sum that is conditionally convergent. The RSGDF uses the proper Ewald formulation (SR real-space + LR reciprocal-space), which is the CORRECT periodic Coulomb integral.

  • On molecular-limit cells (He/20-bohr box), SR+LR matches bare to 1.9% (2c) and 1.4% (3c). On tight cells (MgO/3.98-bohr), the bare sum with 12-bohr cutoff diverges from the RSGDF sum.

  • SCF gauge fix validated: the existing Ewald-J + real-space-K fallback in periodic_rhf_gdf.py gives correct energies (the pinned 3D H2 test matches EWALD_3D to machine precision).

  • RSGDF path forward: tight-cell RSGDF needs the full reciprocal-space FFT approach (PySCF RSDFBuilder method) — compute the 2c/3c integrals on a dense FFT grid rather than a sparse reciprocal-lattice mesh. This is a v0.10.x item.

  • Status: INVESTIGATION COMPLETE. RSGDF works for molecular- limit cells; tight-cell RSGDF is deferred to the FFT path.

M4: Restore real density fitting for dim=3 (deferred)

  • Only applicable once tight-cell RSGDF Lpq is validated (M3). Current Ewald fallback is correct and will remain the default until the FFT-based RSGDF path lands.

M5c: kpoints unification (2026-05-21)

  • gdf_kmesh merged into kpoints parameter — works for both GDF and BIPOLE.

  • run_periodic_job(sys, basis, kpoints=(2,2,2)) dispatches to multi-k driver.

M5d: RKS energy formula fix (2026-05-21)

  • Multi-k DFT energy was ~787 mHa off due to double-counting V_xc in the ½ Tr[D(H+F)] formula. Fixed by separating F_2e from full F; E_elec = Tr[D·Hcore] + ½ Tr[D·F_2e].

  • RKS PBE multi-k now matches gamma to 0.000 mHa.

M5b: Multi-k energy decomposition (2026-05-21)

  • Recover E_coulomb and E_hf_exchange from a second Ewald Fock build with exchange_scale=0. J-only Fock gives J(k); F_full - F_J gives K(k). Cost: +1 Fock build per SCF iteration.

M5: Multi-k GDF SCF driver — MG-3D-MK (gauge fix landed)

  • Gauge fix (2026-05-20): Replaced bare-Lpq J/K with build_periodic_fock_ewald3d_k (same composed Ewald-3D Fock builder as the multi-k EWALD_3D driver). On H2/12-bohr [2,2,2], energy now matches the gamma driver to 0.000 mHa (was −447 mHa).

  • Not wired into run_periodic_job: still direct API only. Needs k-mesh parameter plumbing in the periodic runner.

  • Lpq cache is still built (kept for potential RSGDF/compcell use) but J/K now bypass it entirely for the gauge-correct path.

Future (v0.10.x)

  • Compensated-charge GDF (compcell) — PySCF’s _CCGDFBuilder equivalent for diffuse aux bases. Tracked in docs/design_native_gdf.md § slice 3d.

Key decisions

  1. Per-L calibration is √(4π/(2L+1)) — derived from the angular-average ratio between libint’s “scalar charge” convention and the AO “spherical with Y_lm” convention. Verified empirically on He/def2-svp-jk to ~1%.

  2. G=0 correction uses spherical BZ approximation — the discrete reciprocal-lattice sum excludes G=0; the missing integral is approximated as (2√π/π)·ω·erf(G_BZ/2ω) integrating over a sphere of radius G_BZ = (6π²/Ω)^(1/3). This is exact for isotropic functions (s-shells); for higher-L shells the G=0 contribution vanishes (ξ̂ ∝ G^L → 0).

  3. AO-pair FT calibration is applied externally — rather than modifying _aopair_ft.py, the per-L factors are applied in rsgdf_lr_3c_tensor via _ao_scales_for_rsgdf. This keeps the AO-pair FT module convention-agnostic (it can serve both bare and RSGDF algorithms).

Blockers

  • RSGDF for tight cells: needs full FFT-based reciprocal-space integration (PySCF RSDFBuilder approach). Current sparse-G-mesh

    • spherical BZ correction is adequate for molecular-limit cells (~1.6% Lpq ERI error) but gives ~55 mHa SCF energy error on 1D systems. Deferred to v0.10.x.

  • Multi-k GDF: driver exists but needs Ewald gauge fix and k-mesh plumbing into run_periodic_job. See M5 above.

Resumption checklist

# Verify current state
cd vibeqc-gdf
git fetch && git rebase origin/main
.venv/bin/pip install -e .
.venv/bin/python -m pytest tests/test_periodic_rhf_gdf.py tests/test_periodic_aux_eri_blocks.py -x -q

# Quick calibration smoke
.venv/bin/python -c "
from vibeqc.aux_basis import make_aux_basis_set, rsgdf_lr_2c_metric, rsgdf_g_mesh, build_lpq_native
from vibeqc._vibeqc_core import compute_2c_eri_lattice, compute_2c_eri_lattice_sr, PeriodicSystem, LatticeSumOptions, Atom, BasisSet
box=20.0; lat=np.eye(3)*box
sys=PeriodicSystem(3,lat,[Atom(2,[box/2]*3)])
mol=sys.unit_cell_molecule()
ao=BasisSet(mol,'sto-3g'); aux=make_aux_basis_set(mol,aux_name='def2-svp-jk')
opts=LatticeSumOptions(); opts.cutoff_bohr=2.0; omega=0.4
M=np.asarray(compute_2c_eri_lattice(aux,sys,opts))
Ms=np.asarray(compute_2c_eri_lattice_sr(aux,sys,opts,omega))
gm=rsgdf_g_mesh(sys,omega,precision=1e-10)
Ml=rsgdf_lr_2c_metric(aux,sys,omega,g_mesh=gm)
d=np.linalg.norm(Ms+Ml-M)/np.linalg.norm(M)
print(f'2c calib: {d:.4e} (target < 0.02)')
Lpq_b=build_lpq_native(sys,ao,aux,lat_opts=opts,algorithm='bare')
Lpq_r=build_lpq_native(sys,ao,aux,lat_opts=opts,algorithm='rsgdf',rsgdf_omega=omega)
eb=np.einsum('Lij,Lkl->ijkl',Lpq_b,Lpq_b).reshape(1,1)
er=np.einsum('Lij,Lkl->ijkl',Lpq_r,Lpq_r).reshape(1,1)
d2=abs(eb[0,0]-er[0,0])/abs(eb[0,0])
print(f'Lpq diff: {d2:.4e} (target < 0.02)')
"