BIPOLE handover — v0.8.x complete, v0.9.0+ scoped (2026-05-18)

This handover documents the state of the BIPOLE periodic SCF re-implementation after the 2026-05-17/18 push that delivered the Γ-only Ewald-J-split production path. It also captures the half-baked items that need follow-up chats to close out before the original v0.8.0 BIPOLE deliverable (sub-mHa CRYSTAL parity on the canonical demo set) is fully landed.

What’s complete on origin/main

Driver + tests

  • vibeqc.pbc_bipole.run_pbc_bipole_rhf — production driver with use_ewald_j_split=True + DIIS supported. CRYSTAL-gauge- equivalent F^{2e} build: J^SR(direct erfc-screened) + J^LR(reciprocal-sum analytic) − ½K(direct), single shared α across V_ne / E_nn / J_LR (CRYSTAL’s COMMON/VRSMAD/ pattern).

  • vibeqc.bipole_fock_ewald.build_fock_2e_ewald_j_split_gamma — standalone Γ-only Fock build for direct use.

  • Phase 1-4b infrastructure:

    • cpp/src/multipole_moments_lattice.cpp (Phase 1 — libint emultipole{1,2,3} shell-pair multipoles).

    • python/vibeqc/bipole_multipole.py (Phase 2 — multipole- multipole interaction tensor).

    • python/vibeqc/bipole_dispatch.py (Phase 3 — IDIPC geometric dispatch).

    • python/vibeqc/bipole_cell_moments.py (Phase 4a — cell-level moments from density).

    • python/vibeqc/bipole_lattice_self_energy.py (Phase 4b — dipole-dipole depolarization placeholder).

    • python/vibeqc/bipole_ext_el_pole.py (Phase 4b proper — Ewald reciprocal-sum framework).

  • 73 BIPOLE-branch tests passing + 1 intentional skip.

  • 15 CRYSTAL eval-version demo geometries via examples/regression/crystal_parity/crystal_demos/builders.py.

  • 6 working parity scripts: 3 Γ-only Ewald-J-split (MgO, diamond, Si — all PASS) + 3 legacy ones (mgo_primitive_shrink88, etc).

  • 12 stubbed parity scripts for blocked demos (UHF / 2D / 1D / inline-basis / SG-expansion).

  • User docs updated: docs/user_guide/bipole.md includes quick-start + empirical parity table.

  • 3 memory files documenting the work + findings + roadmap.

Empirical CRYSTAL parity (final)

All 3 RUNNABLE STO-3G demos PASS via run_pbc_bipole_rhf(use_ewald_j_split=True, use_diis=True):

Demo

iters

E_total

Δ vs CRYSTAL SHRINK 8 8

Si

8

-570.876 Ha

+444 mHa

MgO

9

-269.337 Ha

+1881 mHa

diamond

12

-72.348 Ha

+2529 mHa

All converged at dE < 1e-7. Si has the smallest Δ (largest lattice → smoothest k-mesh landscape); diamond the largest (covalent + small lattice + steep band dispersion).

What’s HALF-BAKED — v0.9.0+ followups

#1: Multi-k F_J^LR(k) proper [HIGHEST IMPACT]

Problem. The Γ-vs-converged-k-mesh shift (+0.4 to +2.5 Ha above) is the physics gap between vibe-qc and CRYSTAL. Closes when multi-k is wired properly.

What was tried + reverted. A k-independent J^LR(Γ) shortcut (added at g=0 only) was tried — empirically OVER-BINDS at iter 2+ on MgO SHRINK 2 2 (-287 Ha at iter 3 vs CRYSTAL’s -271 Ha). Reverted in commit 6cc5eab. Two compounding errors in the shortcut:

  1. F_J^LR(k) is k-dependent for non-trivial bra-pair-overlap density blocks. The k-independent treatment is only exact at iter 1 SAD (where D(g≠0)=0 → ρ̂(K) is k-independent).

  2. ρ̂(K) for D(g≠0) needs proper shifted-ν AO-pair FT, not the bra-pair-at-home phase approximation in compute_cell_density_fourier. For Gaussian basis, shifting ν by R_g changes the Gaussian-product center P_g non-linearly, so FT_μν(K; g) ≠ FT_μν(K; 0) · exp(-iK·R_g).

What’s needed. Two closed-form derivations + implementations:

  1. Per-k F_J^LR(k) via Bloch-summed bra-pair FT: F_J^LR(k)_μν = (4π/V) Σ_K kernel(K) · ρ̂_full(K) · FT*_μν^{bloch}(K; k) where FT_μν^{bloch}(K; k) = Σ_g exp(ik·R_g) · FT_μν(K; g). The Gaussian-primitive Bloch-summed FT has a Jacobi-theta-like structure that becomes a Dirac-comb-like selection rule in reciprocal space — non-trivial closed form.

  2. Correct ρ̂(K) for D(g≠0) via shifted-ν AO-pair FT: ρ̂_full(K) = Σ_{μν, g} P(g)_μν · FT_μν(K; g) where FT_μν(K; g) = K_pq^g · (π/γ)^{3/2} · exp(-K²/4γ) · exp(-iK·P_g) with P_g and K_pq^g both depending on R_g.

References. Stone “Theory of Intermolecular Forces” §3. Saunders / Dovesi / Roetti 1992. PySCF’s pyscf.pbc.df.aft for the AO Fourier transform machinery.

Effort estimate. ~1 day of math + implementation + validation against CRYSTAL on MgO SHRINK 2 2 → SHRINK 8 8.

Files to touch.

  • python/vibeqc/_aopair_ft.py (extend with Bloch-summed FT).

  • python/vibeqc/bipole_ext_el_pole.py (update compute_cell_density_fourier to use shifted-ν FT).

  • python/vibeqc/bipole_fock_ewald.py (extend compute_J_long_range_gamma to a per-k variant).

  • python/vibeqc/pbc_bipole.py (remove the n_k != 1 guard; wire the per-k F_J^LR build).

  • Add tests/test_pbc_bipole_multik_ewald_split.py.

  • Add examples/regression/crystal_parity/parity_mgo_shrink88_ewald_jsplit.py.

#2: EXT EL-SPHEROPOLE-equivalent [DEFERRED TO bipole-multik — 2026-05-18]

Problem. CRYSTAL prints ::: EXT EL-SPHEROPOLE = +4.119 Ha on MgO at CYC 0 (+3.964 Ha converged), separate from EXT EL-POLE. This is the K=0 spheropole-coupling term of the Ewald reciprocal- space sum.

Investigation (2026-05-18, this entry). CRYSTAL source review at ~/gitlab/crystal/src/libmondep.f::CJAT33 (3D) + libxy.f::CJAT3 shows the per-primitive-pair s-s integrand is::

SHIFT_ab  =  CFAC1(LOFA) · (3/γ_p + (FACA² + FACB²)·R²)
                        · γ_p^(-1/2) · TOTCHR / (2V)

with TOTCHR = N_electrons / 6 (ints2.f:6004), γ_p the per- primitive-pair Gaussian exponent sum (α_a + α_b), and the integrand being the bond-symmetrised second moment ⟨(r−A)² + (r−B)²⟩ divided by the s-s overlap.

The full memory entry is reference_bipole_ext_spheropole_2026-05-18.md.

Why naive (constant · Tr(Q)) doesn’t match. The γ_p^(-1/2) weighting is per-primitive-pair — different contracted pairs carry different prefactors — so the formula does NOT factorise into (const) · Tr(M_total). Plus higher-l shells get extra bond-derivative terms. Implementing it exactly requires either a custom C++ primitive-pair integral kernel or a proper K→0 Taylor expansion of the Bloch-summed Ewald reciprocal sum.

Decision: deferred to bipole-multik. Reasons:

  1. Dominant residual is multi-k, not spheropole. The Δ vs CRYSTAL on MgO Γ-only is +1.88 Ha; the spheropole at converged density is +3.96 Ha. Adding it naively pushes Δ to +5.85 Ha (worse).

  2. The bipole-multik per-k F_J^LR(k) work should naturally subsume it. When the K=0 limit of |ρ̂(K+k)|²·exp(-K²/4α²)/K² is taken correctly (Taylor expansion of ρ̂ around K=0 for a neutral cell), the dipole-layer and quadrupole/spheropole contributions emerge without a separate analytic term.

  3. Implementation cost is high. A standalone EXT EL-SPHEROPOLE kernel would need a custom primitive-pair integral with γ_p^(-1/2) weighting — not a libint emultipole2.

Sanity check for the bipole-multik chat. Once per-k F_J^LR(k) is wired, isolate the K=0 contribution and verify it matches CRYSTAL’s printed EXT EL-SPHEROPOLE value within a few mHa on MgO/diamond/Si. If the post-fix Γ-vs-SHRINK 8 8 residual is still >1 mHa, the spheropole pathway may need its own analytic term after all — but that’s an empirical decision, not a structural prerequisite.

Effort estimate. Subsumed by bipole-multik (~1 day). No separate implementation in this chat.

#3: BIPOLE C++ multipole-far-pair branch [PERFORMANCE]

Problem. Phase 5 ships as Python orchestration on top of the existing C++ build_fock_2e_real_space direct kernel. For production use on larger systems / extended basis sets, the direct-space ERI cost dominates.

What’s needed. Promote the Phase 1-4 multipole infrastructure to a C++ build_fock_2e_bipole kernel that:

  • Per-quartet IDIPC dispatch (Phase 3) inside the C++ loop.

  • For near quartets: call libint exact ERI (existing path).

  • For far quartets: compute via multipole expansion (Phase 1 moments × Phase 2 interaction tensor).

  • Phase 4b EXT EL-POLE correction applied as a post-pass.

Effort estimate. ~2-4 weeks of C++ work + testing.

Files. New cpp/src/build_fock_2e_bipole.cpp + binding + extensive testing.

#4: UHF + extended-basis + dim<3 unblocks [PARALLEL TRACKS]

These would activate the 12 currently-stubbed parity scripts:

  • UHF BIPOLE driver (run_pbc_bipole_uhf) — unblocks parity_nio_sto3g.py. Mirror of periodic_uhf_multi_k_ewald.py but on the BIPOLE / Ewald-J- split path. ~1 day of work.

  • Inline-basis parser for CRYSTAL .d12 atom-specific basis blocks — unblocks 6 extended-basis variants (mgo_bulk, diamond_extended, sibulk_extended, etc.). Likely belongs on basissetdev. ~1-2 days.

  • Metallic SCF (Fermi-Dirac smearing + chemical potential bisection) — unblocks parity_be_sto3g.py. ~1 day.

  • dim = 1 / dim = 2 BIPOLE path — unblocks parity_graphite_sto3g.py + parity_sn_polym_sto3g.py. ~2 days (needs low-dim Ewald handling).

  • Ghost-atom Z=108 + asymmetric slabs — unblocks parity_mgo001co.py. ~1 day.

  • spglib-driven SG asymmetric-unit expansion — unblocks parity_urea_321G.py. ~4 hours.

Files to consult

  • python/vibeqc/pbc_bipole.py — current driver, contains the Γ-only short-circuit + multi-k NotImplementedError guard.

  • python/vibeqc/bipole_fock_ewald.py — Phase 5 Ewald J-split build (currently Γ-only).

  • python/vibeqc/bipole_ext_el_pole.py — Phase 4b reciprocal-sum framework, including compute_cell_density_fourier (the bra-pair-at-home phase approximation that needs fixing).

  • python/vibeqc/_aopair_ft.py — existing AO-pair FT (home shells only). Needs Bloch-summed extension for multi-k.

  • examples/regression/crystal_parity/parity_{mgo,diamond,sibulk}_gamma_ewald_jsplit.py — Γ-only validated parity scripts (templates for multi-k variants).

  • ~/.claude/projects/-Users-mpei-gitlab-vibeqc-bipole/memory/ — three relevant memory files:

    • reference_bipole_multipole_required_2026-05-17.md — structural finding that BIPOLE multipole branch is required.

    • reference_bipole_phase5_ewald_jsplit_works_2026-05-18.md — Phase 5 success.

    • reference_ewald_composition_bug_2026-05-17.md — the bug that motivated the BIPOLE rewrite.