Design — native vibe-qc GDF (Method 1 of the periodic JK roadmap)

Status (2026-05-09, feature/native-gdf-chain): in progress. Slice 1 (C++ periodic 2c/3c ERI kernels) and slice 2 (Python aux loader

  • Cholesky-fitted Lpq) have shipped on the feature/native-gdf-chain branch. Slice 3 (modrho compensation) is empirically confirmed mandatory — see “Slice 3 — modrho compensation is required” below. Slices 4–8 (V_ne via aux, driver rewrite, parity sweep, docs) are queued behind slice 3.

Goal: drop the PySCF dependency from the GDF spike. Build the Lpq tensor entirely from vibe-qc’s libint pipeline + an in-house auxiliary basis manager.

What Lpq is, recap

M_PQ      = (P|Q)_periodic           ← 2c periodic ERI on aux basis
(μν|P)_pp = Σ_T (μ_0 ν_0 | P_T)      ← 3c periodic ERI w/ ket image sum
Lpq[L, μ, ν] = Σ_P (μν|P)_pp · M^{-1/2}_{PL}      (Cholesky factor)

Then for closed-shell RHF (or KS):

ρ_L  = Σ_λσ Lpq[L, λ, σ] · D_λσ
J_μν = Σ_L Lpq[L, μ, ν] · ρ_L
K_μν = Σ_L (Lpq[L] @ D @ Lpq[L].T)_μν
F_2e = J − ½ K  (+ Madelung correction for exxdiv='ewald')

Pieces vibe-qc needs to add

1. C++ — 3c and 2c molecular ERI engines

Wraps libint2’s existing 3-center and 2-center engines.

cpp/include/vibeqc/aux_eri.hpp (new):

// 3c molecular: (μν|P) tensor (n_ao, n_ao, n_aux), no lattice sum.
Eigen::Tensor<double, 3> compute_3c_eri(
    const BasisSet& ao_basis,
    const BasisSet& aux_basis);

// 2c molecular: (P|Q) matrix (n_aux, n_aux), no lattice sum.
Eigen::MatrixXd compute_2c_eri(const BasisSet& aux_basis);

cpp/src/aux_eri.cpp:

  • Use libint2::Operator::coulomb with Engine::compute(s1, s2, s3) for 3c and Engine::compute(s1, s2) for 2c.

  • Mirror the molecular compute_eri already in cpp/src/integrals.cpp for the threading + buffer layout.

  • Output is ROW-major libint AO order; permutation to PySCF order (if needed) lives in Python.

2. C++ — periodic versions

cpp/include/vibeqc/aux_eri.hpp:

// Periodic 3c: Σ_T (μ_0 ν_0 | P_T) tensor (n_ao, n_ao, n_aux).
Eigen::Tensor<double, 3> compute_3c_eri_lattice(
    const BasisSet& ao_basis,
    const BasisSet& aux_basis,
    const PeriodicSystem& system,
    const LatticeSumOptions& opts);

// Periodic 2c: Σ_T (P_0 | Q_T) matrix (n_aux, n_aux).
Eigen::MatrixXd compute_2c_eri_lattice(
    const BasisSet& aux_basis,
    const PeriodicSystem& system,
    const LatticeSumOptions& opts);

Implementation: mirror cpp/src/lattice_integrals.cpp cell loops with shifted shells, accumulate. Schwarz screening on 3c can prune many quartet calls.

For the long-range part of the periodic 2c metric (where the 1/r sum diverges in real space), apply Ewald-style splitting: (P|Q)_periodic = (P|Q)_SR(ω) + (P|Q)_LR(ω) where SR is real-space erfc(ω)/r and LR is reciprocal-space erf(ω)/r. The LR part is a small G-mesh sum, NOT an FFT (the aux basis is small enough that direct G-summation works).

3. C++ — pybind11 bindings

Add to cpp/src/bindings.cpp:

m.def("compute_3c_eri",         &compute_3c_eri,
      "AO 3-center ERI tensor (μν|P).");
m.def("compute_2c_eri",         &compute_2c_eri,
      "AO 2-center ERI matrix M_PQ = (P|Q).");
m.def("compute_3c_eri_lattice", &compute_3c_eri_lattice,
      "Periodic 3-center ERI tensor with image sum on the ket.");
m.def("compute_2c_eri_lattice", &compute_2c_eri_lattice,
      "Periodic 2-center ERI metric matrix.");

4. Python — auxiliary basis manager

python/vibeqc/aux_basis.py (new):

def make_aux_basis_set(ao_basis: BasisSet, *, aux_name: Optional[str] = None,
                      drop_eta: float = 0.0) -> BasisSet:
    """Build an auxiliary BasisSet for an AO basis.

    Mirrors PySCF's pyscf.df.addons.predefined_auxbasis +
    pyscf.pbc.df.df.make_modrho_basis:
      - If aux_name is None, look up the standard JKfit aux for the
        AO basis (sto-3g → universal-jkfit, def2-tzvp → def2-jkfit, ...).
      - If aux_name is a string, load it from the libint .g94 path.
      - Drop primitives with exponents below drop_eta to avoid
        linear-dependent metric (matches PySCF make_modrho_basis).
      - Renormalize aux primitives so the monopole moment ∫ χ_P dr =
        √(1/(4π)). This convention simplifies the compensated-charge
        algebra (see PySCF make_modrho_basis docstring).
    """
    ...

Reference: pyscf/pbc/df/df.py:make_modrho_basis (lines 64–120).

5. Python — Lpq builder

In python/vibeqc/periodic_rhf_gdf.py, replace _lpq_tensor_from_pyscf with:

def _lpq_tensor_native(system, ao_basis, *, aux_name=None,
                      lat_opts=None, linear_dep_thr=1e-9):
    """Build Lpq(L, μ, ν) entirely with vibe-qc machinery.

    1. Build aux basis from ao_basis (or load aux_name).
    2. Compute periodic 2c metric: M[P, Q] = Σ_T (P_0 | Q_T)
       via compute_2c_eri_lattice (or with Ewald split if needed).
    3. Compute periodic 3c: T[μ, ν, P] = Σ_T (μ_0 ν_0 | P_T)
       via compute_3c_eri_lattice.
    4. Cholesky factor M = L L^T (Lehtola-style pivoted Cholesky for
       robustness against the diffuse-aux linear dependence).
    5. Solve L · X = T.reshape(naux, nao*nao) → Lpq.
    """
    aux_basis = make_aux_basis_set(ao_basis, aux_name=aux_name)
    M = compute_2c_eri_lattice(aux_basis, system, lat_opts)
    T = compute_3c_eri_lattice(ao_basis, aux_basis, system, lat_opts)

    # Pivoted Cholesky tolerates near-singular M (diffuse aux).
    from scipy.linalg import cho_factor, cho_solve
    L_chol, low = cho_factor(M, lower=True)
    Lpq = cho_solve((L_chol, low), T.reshape(M.shape[0], -1))
    return Lpq.reshape(M.shape[0], T.shape[0], T.shape[1])

6. Validation

Hard test: compare _lpq_tensor_native(cell, basis) to the existing _lpq_tensor_from_pyscf(cell) elementwise. Should match to ~1e-10 (Cholesky residual, identical aux basis, identical integrals).

Then drop the PySCF Lpq path. The full SCF chain becomes self-hosted:

S, T, V_ne   from vibe-qc compute_*_lattice (already)
Lpq          from native vibe-qc 3c/2c        (shipped for Γ)
cell blocks  from compute_*_eri_lattice_blocks (multi-k storage boundary)
J, K          from native einsum on Lpq        (already for Γ)
V_xc          from native libxc on Becke grids

PySCF and CRYSTAL are no longer in-process periodic backends. They are external reference programs whose outputs are parsed for parity checks. The remaining native multi-k GDF work is the k-dependent metric/tensor phase assembly, IBZ-weighted J/K contractions, and disk-backed storage.

Estimated effort

  • 3c/2c molecular bindings: 2–3 hours.

  • Periodic versions with Schwarz screening: 4–6 hours.

  • Pybind bindings + CMake: 30 min.

  • Aux basis manager + JKfit table loading: 2–3 hours.

  • Lpq builder + Cholesky: 1 hour.

  • Validation + debug: 4–8 hours (typical: aux normalization, linear- dep handling for diffuse bases, Schwarz threshold tuning).

Total: 1.5–2 focused days. Doable in one follow-up session.

Auxiliary basis library — what to ship

For sto-3g and minimal bases: PySCF auto-generates “even-tempered” aux when no explicit aux is registered. We can do the same: emit even-tempered Gaussians from the AO exponent range. Quick fix that unblocks the full sweep.

For the Bredow-group periodic bases (pob-tzvp, pob-tzvp-rev2): we already have auxiliary basis design as a paper-worthy item (see memory: project_pob_tzvp_aux_basis.md). Ship a prototype auto-aux for v0.7.x; the proper hand-tuned aux is a separate deliverable.

For def2-* bases: copy PySCF’s def2-jkfit / def2-jkfit-tight tables into vibe-qc’s basis_library/aux/.

Slice 3 — modrho compensation is required

The image-summed periodic 2c metric M_PQ = Σ_T (P_0 | Q_T) on a non-charge-compensated aux basis is mathematically divergent in 3D — the leading monopole-monopole term sums as Σ 1/|T| over a 3D lattice, which is the classical Madelung divergence. This was an unresolved question at design time; the empirical confirmation lives in examples/debug/scratch_lpq_native_smoke.py and examples/debug/scratch_lpq_native_energy.py. On MgO/sto-3g with def2-universal-jkfit (PySCF’s auto-pick):

cutoff (bohr)

‖M‖_F

min eig(M)

Cholesky

4

1.6e+3

+4.8e−6

✓

8

6.0e+3

−7.6e+1

✗

12

1.2e+4

−2.2e+1

✗

16

1.8e+4

−9.5e+0

✗

20

3.4e+4

−4.4e+0

✗

‖M‖_F grows monotonically with cutoff (no convergence) and the metric loses positive-definiteness past the Cholesky-OK cutoff (~4 bohr — barely larger than the 8 bohr cubic cell, so the SR aux contributions dominate while the LR divergence has not yet kicked in).

Iter-1 SCF energy with the SAD initial density:

Lpq source

E_iter1 (Ha)

comment

(A) PySCF Lpq baseline

−1084.108

converges to −1085.231

(B) Native cut=4 (Cholesky)

−270.229

off by ~800 Ha

(C) Native cut=20 (SVD pseudo)

−12,653.152

catastrophic

so the bare lattice sum is unusable for SCF. Modrho compensation is structurally required, not an optimisation.

Modrho — algorithmic note

PySCF’s pyscf.pbc.df.df.make_modrho_basis rescales each contracted aux shell so that its monopole moment is √(1/(4π)):

s_i  = Σ_p c[p, i] · gaussian_int(2L+2, α_p)
c[p, i] *= √(1/(4π)) / s_i

For l > 0 the renormalisation is conventional (the actual electric monopole vanishes by symmetry); for l = 0 it caps each S-shell’s monopole at the same finite value. Combined with PySCF’s separate charge-compensation on the AO-pair side (pyscf.pbc.df.df_builder._CCGDFBuilder.make_j3c, smooth Gaussians that cancel the AO pair’s monopole), the resulting (modified-aux | compensated-AO-pair) integral converges in real space.

For vibe-qc to reproduce this:

  1. Ship a way to construct a BasisSet with custom contraction coefficients. Two options:

    • (preferred) Add BasisSet(Molecule, vector<ShellInfo>) C++ constructor that builds libint2::Shells directly. Bypasses the .g94 round-trip; avoids file I/O races; clean API.

    • (workaround) Write a temp .g94 file under python/vibeqc/basis_library/basis/<unique_name>.g94 with denormalised coefficients, load via BasisSet(mol, name), clean up. Works without C++ changes but ugly. The denormalisation dance (libint applies primitive norm at load; PySCF’s modrho coeffs are post-normalised; the conversion factor must be exact) is error-prone.

  2. Implement modrho in vibeqc.aux_basis using vibe-qc’s own gaussian_int(2l+2, α) helper (port the 3-line PySCF function). Apply per-shell rescaling on the libint-internal coefficients.

  3. Implement compensated-AO-pair handling on the 3c side. The compensating smooth Gaussians become a separate compcell BasisSet; the periodic 3c is then computed on (aux | AO_pair − comp) and the missing (aux | comp) piece is added back analytically (single G-mesh sum, see PySCF’s _CCGDFBuilder.make_j3c ~600 LOC for the canonical algorithm).

Estimated additional effort: 1–2 days (modrho + compcell + a Python-side reciprocal-space sum for the compensating term + parity validation). The slice-1 C++ kernels are reusable as-is.

Resumption checklist for the next session

  1. Read this section + check the live periodic-SCF chats in vibeqc-queue for current state; the in-tree docs/handover_periodic_scf_followup.md brief was removed in the 2026-05-17 archive sweep (recover via git show HEAD~:docs/handover_periodic_scf_followup.md).

  2. Confirm slice-1 baseline still passes: .venv/bin/python examples/debug/scratch_aux_eri_lattice_smoke.py.

  3. Confirm slice-3b (modrho) reduces metric ill-conditioning by ~2 orders of magnitude but does NOT make the metric SPD past cutoff=4 bohr — this is the open question for slice 3c–d.

Update 2026-05-09 — slices 3a and 3b shipped

Slice 3a (BasisSet from explicit ShellInfo) — committed (0ed5854). Adds BasisSet(molecule, shells, name, coefficients_pre_normalized=True) C++ constructor + Python-mutable ShellInfo so modrho can run in pure Python without .g94 file round-trips. Round-trip on H2/sto-3g shows S and ERI agree exactly; scaling all coefficients by 2 produces overlap × 4, validating that the no-renorm path is honoured.

Slice 3b (modrho normalisation) — committed (5778e39). Direct port of PySCF make_modrho_basis:

s_i = Σ_p c[p, i] · gaussian_int(2L+2, α_p)
c[p, i] *= √(1/(4π)) / s_i

plus an optional drop_eta primitive cull. Wired into make_aux_basis_set via compensate=True (default).

Empirical effect on MgO/sto-3g/def2-universal-jkfit (image-summed 2c metric, increasing cutoff):

cutoff (bohr)

min eig(M) raw

min eig(M) modrho

min eig(M) modrho+drop=0.2

4

+4.8e−6

+3.5e−7

+5.3e−6

8

−7.6e+1

−7.2e−1

−3.6e−1

12

−2.2e+1

−1.0e−1

−2.3e−2

16

−9.5e+0

−2.6e−2

−4.5e−3

20

−4.4e+0

−8.6e−3

−4.7e−4

Modrho improves conditioning by ~2 orders of magnitude, but the metric remains non-SPD past cut=4 bohr. The 3D Madelung divergence is not cured by normalisation alone — this is the expected behaviour per the PySCF algorithm structure: modrho is a normalisation convention that simplifies downstream algebra; the actual convergence treatment lives in the AO-pair-side compensating-charge mechanism (slice 3d below).

Slice 3c–d — open algorithmic questions

Two paths are plausible from here; the literature sweep commissioned in parallel should resolve which is correct:

(i) Per-shell precision-truncated lattice cutoff

Hypothesis: PySCF’s pbc_intor("int2c2e") on the modrho’d auxcell returns a finite matrix (max |M| ≈ 20 on MgO/sto-3g) by using its cell.precision-controlled image cutoff via _estimate_rcut(es, l, c_max, precision). Each aux shell contributes to images only out to its own precision-justified radius; diffuse primitives have small radii so their long-range divergent contribution is naturally truncated.

If true, vibe-qc’s bare 2c lattice sum on modrho aux + per-shell rcut should reproduce PySCF’s pbc_intor("int2c2e") matrix to ~1e-10. Then the metric is “PSD by truncation” — physically accurate to the precision target, not to infinity.

Implementation cost: small. Add per-shell-pair rcut to the cell loop in compute_2c_eri_lattice / compute_3c_eri_lattice.

(ii) Full AO-pair compensating-charge mechanism

The textbook fix per Sun et al. 2017 § II.B: introduce a “compcell” of smooth Gaussians (one per atom × per l, with global smooth_eta). For each AO pair (μ, ν) at center A, compute its multipole moments; subtract a matching combination of compcell Gaussians; compute (aux | compensated AO pair) in real space (converges); add back analytical (aux | smooth-Gaussian) terms in reciprocal space.

Substantially more code:

  • Compcell BasisSet construction (make_chgcell in PySCF).

  • AO-pair multipole moments (libint emultipoleN — already used in compute_dipole; extend to higher l).

  • (aux | compcell) integrals — special case of the existing 3c lattice with compcell as the orbital pair.

  • Reciprocal-space (aux | smooth) sum — small G-mesh, analytic Fourier transforms of Gaussians.

Estimate: 2–3 days, plus parity validation.

Decision rule

Try (i) first — small change, may suffice for sub-µHa parity if our hypothesis about PySCF’s truncation is right. If (i) doesn’t reach sub-µHa, the lit sweep’s PySCF algorithm walkthrough should clarify exactly what (ii) needs to look like.

Resumption checklist (updated 2026-05-09)

  1. Read this section + the handover docs.

  2. Verify slices 1–3b still work:

    • examples/debug/scratch_aux_eri_lattice_smoke.py (slice 1 structural check).

    • Quick modrho check via make_aux_basis_set(mol, aux_name="def2-universal-jkfit", compensate=True) and a manual ‖M‖_F sweep — should reproduce the table above.

  3. Implement per-shell rcut (path (i) — start here).

  4. Compare vibe-qc’s 2c metric to PySCF’s auxcell.pbc_intor(int2c2e) elementwise. If matched to ~1e-10, ship as slice 3c.

  5. If not matched, defer to lit-sweep results for the proper compcell algorithm and implement path (ii) as slice 3d.

  6. Once Lpq parity is achieved, resume slices 4 (V_ne via aux), 5 (driver rewrite), 6 (_basis_g94 cleanup), 7 (full parity sweep), 8 (docs + PR).