BIPOLE handover — EXT EL-SPHEROPOLE partial implementation (2026-05-20)¶

This handover covers the EXT EL-SPHEROPOLE work item from docs/handover_bipole_v0_8_2026_05_18.md §2. A partial implementation landed; this doc tells a follow-up chat exactly what works, what doesn’t, and how to finish it.

TL;DR¶

compute_ext_el_spheropole exists in python/vibeqc/bipole_ext_el_pole.py (commit ab5fef8).
It captures only the s-s-equivalent term of CRYSTAL’s per-cycle EXT EL-SPHEROPOLE. On MgO/STO-3G SAD it returns +0.71 Ha vs CRYSTAL’s +4.119 Ha (~17%).
It is NOT wired into run_pbc_bipole_rhf — adding an 80%-deficit correction to E_total would worsen parity.
The remaining work is the higher-l (l_a + l_b ≥ 1) Cartesian-moment-derivative terms of CRYSTAL’s CJAT33.

What CRYSTAL’s EXT EL-SPHEROPOLE is¶

CRYSTAL’s per-cycle ENECYCLE block breaks the electron-electron energy into three pieces:

TOTAL E-E = BIELET ZONE E-E + EXT EL-POLE + EXT EL-SPHEROPOLE
MgO/STO-3G CYC 0:  46.211 = +570.696 + (-528.604) + (+4.119)

EXT EL-SPHEROPOLE is the K=0 spheropole-coupling term of the Ewald reciprocal-space sum — the finite K→0 limit of |ρ̂(K)|²·exp(-K²/4α²)/K² for a neutral cell, evaluated by CRYSTAL as a real-space per-shell-pair integral.

Sealed CRYSTAL references (in examples/regression/crystal_parity/crystal_demos/):

System	CYC 0 (SAD)	converged
MgO STO-3G	+4.119189 Ha	+3.964 Ha
diamond STO-3G	+3.315182 Ha	—

CRYSTAL source — the exact formula¶

The accumulator is SHIF (historical name; not level-shift). Traced through:

~/gitlab/crystal/src/libmondep.f — CJAT33 (3D, line 691), MONMO3/TOTENY energy printers (SHIF accumulation at :1568, :1591; FORMAT 202 print at :1825).
~/gitlab/crystal/src/libxy.f — CJAT3 (line 824), the alternative 3D path with the same SHIFT integrand.
~/gitlab/crystal/src/libx5_com.f — MADEL2 Ewald setup (POISS = π/V, RPI = 1/2π, TOTCHR).
~/gitlab/crystal/src/ints2.f:6004 — TOTCHR = INF(9)/6 where INF(9) = total nuclear charge (= N_electrons, neutral).
~/gitlab/crystal/src/libxa.f::S1S — confirms CRYSTAL’s C1 contraction coefficients are the raw Pople d^std values (NO primitive normalisation absorbed).

s-s primitive pair (CJAT33 case 33000)¶

SHIFT_ab = d_a^std · d_b^std · γ_p^{-3/2} · exp(-FACX·R²)
         · (3 + γ_p·(FACA² + FACB²)·R²) · TOTCHR/(2V)

with γ_p = α_a + α_b, FACA = α_a/γ_p, FACB = α_b/γ_p, FACX = α_a·α_b/γ_p, R the shell-pair separation, TOTCHR = N_electrons/6, V the cell volume.

The bracket (3/γ_p + (FACA²+FACB²)·R²) is the bond-symmetrised second moment ⟨a|(r−A)² + (r−B)²|b⟩ ÷ ⟨a|b⟩.

higher-l primitive pairs (cases 33001, 33002) — NOT implemented¶

For l_a + l_b = 1 (libxy.f:792-804):

SHIFT = CFAC1(LOFA+IXXX)·F100 + CFAC1(LOFA-MAXP)·F010
      + CFAC1(LOFA-1)·F001 + CFAC1(LOFA)·E100

For l_a + l_b = 2 (libxy.f:806-820): adds a Laplacian-like BILBO·VSQR term where BILBO = CFAC1(LOFA+ITXXX) + CFAC1(LOFA-ITYYY) + CFAC1(LOFA-2).

F100/F010/F001 = R_{x,y,z}·(scaled VSQR) and the CFAC1(LOFA±offset) are higher Cartesian-moment-derivative coefficients of the Gaussian product (DFAC4 output layout). These bond-derivative terms are what the current vibe-qc implementation is missing.

What landed in vibe-qc (commit `ab5fef8`)¶

python/vibeqc/bipole_ext_el_pole.py:

_libint_primitive_norm(l, alpha) — libint’s spherical-Gaussian primitive norm N = (2α/π)^(3/4)·(4α)^(l/2)/√((2l-1)!!).
_split_basis_into_primitives(basis, system) — splits the contracted basis into single-primitive child shells, returns (prim_basis, prim_alpha, prim_origins, prim_N_lib, contraction_map).
compute_ext_el_spheropole(P_real, basis, system, lat_opts) — the energy correction. Algorithm:
1. Split into primitive basis.
2. compute_multipole_moments_lattice(..., L_max=2) → per- primitive Cartesian moments (S, D, Q).
3. Per primitive-AO pair, build the bond-symmetrised second moment `M_ab^g = 2(Q_xx+Q_yy+Q_zz) − 2(A+B+R_g)·D
  - (|A|²+|B+R_g|²)·S`.
4. Weight by γ_ab / (N_a·N_b) (the N undoes libint’s primitive norm baked into the emultipole2 output, recovering CRYSTAL’s d^std·d^std contraction).
5. Map primitive → contracted AO via contraction_map, trace with P_real.
6. Multiply by TOTCHR / (2V·π^(3/2)).

This is exact for s-s shell pairs and a known-incomplete approximation for everything else.

Why it’s energy-only (not in the Fock)¶

CRYSTAL adds SHIFT to the printed energy but not to the Fock matrix — libmondep.f:1592 adds only FRODO + BILBO (= KE + V_ne + EPOLE) to FG_IRR; SHIF is a pure energy accumulator. vibe-qc should follow the same convention: a post-SCF E_total += compute_ext_el_spheropole(...) correction with no ∂/∂D term in F^2e.

Empirical state (measured at `ab5fef8`)¶

examples/regression/crystal_parity/probe_spheropole_sad.py runs the function on the SAD density and compares to CRYSTAL CYC 0:

System	vibe-qc	CRYSTAL CYC 0	ratio
MgO STO-3G	+0.706833 Ha	+4.119189 Ha	0.17
diamond STO-3G	+0.646975 Ha	+3.315182 Ha	0.20

The deficit is dominated by the Mg/O/C 2p shells (l=1), which need the case-33001/33002 terms.

Probe-iteration history (debugging trail)¶

Three probe versions, each adding a correction term:

Version	Per-pair weight	MgO result
v1	`1` (contracted only)	+0.118 Ha
v2	`γ_ab`	+0.190 Ha
v3 (landed)	`γ_ab / (N_a·N_b)`	+0.707 Ha

Each moves toward CRYSTAL but the s-s-only structural limit caps it at ~17%.

Tests¶

tests/test_bipole_ext_el_pole.py:

test_ext_el_spheropole_runs_cleanly — PASS (sanity bounds on MgO SAD: positive, < CRYSTAL upper bound, > 0.3 Ha).
test_ext_el_spheropole_dim_check — PASS (raises for 1D/2D).
test_ext_el_spheropole_matches_crystal_mgo — @pytest.mark.skip (target ±1 mHa vs +4.119 Ha; un-skip once higher-l lands).

How to finish it — two options¶

Option A: C++ CJAT33-mirror kernel (most faithful)¶

Port CRYSTAL’s CJAT33 Cartesian-moment-derivative accumulator directly. This needs the DFAC4 CFAC1 layout (Cartesian-monomial indexing of the Gaussian product) and the LOFA±offset index arithmetic. New cpp/src/ kernel + binding. Highest fidelity, ~1-2 weeks.

Option B: Obara-Saika bond-symmetrised second moment (recommended)¶

Derive ⟨a|(r−A)² + (r−B)²|b⟩ for general (l_a, l_b) via the Obara-Saika recursion, then weight per primitive pair by γ_ab/(N_a·N_b). The current code already has the s-s case via libint emultipole2; the missing piece is that libint’s contracted emultipole2 cannot expose the per-primitive γ_ab weight for l ≥ 1 correctly because the angular factors mix. A clean path: compute the primitive-level ⟨a|…|b⟩ analytically (the integrand is a degree-(l_a+l_b+2) Cartesian polynomial × Gaussian product — closed form), skipping libint entirely. ~2-3 days.

Validation target¶

probe_spheropole_sad.py should report Δ < 1 mHa vs the sealed CRYSTAL CYC 0 values for MgO and diamond. Then un-skip test_ext_el_spheropole_matches_crystal_mgo and wire the term into run_pbc_bipole_rhf as a post-SCF E_total correction.

Coordination¶

The bipole-multik per-k F_J^LR(k) work (commit 602c315, handover §1) closed the dominant Γ-vs-k-mesh shift. The K=0 spheropole partially overlaps with the per-k Ewald K=0 limit — before wiring the spheropole into E_total, verify it is not double-counted against the multi-k path’s K=0 handling.
CRYSTAL-gauge BIPOLE SCF was substantially reworked since ab5fef8 (f3d3dcc Implement CRYSTAL-gauge periodic BIPOLE SCF and follow-ups). compute_ext_el_spheropole only depends on compute_multipole_moments_lattice + BasisSet splitting — stable APIs — so it should still run, but re-measure the probe against current main before trusting the numbers above.

Files¶

python/vibeqc/bipole_ext_el_pole.py — implementation.
tests/test_bipole_ext_el_pole.py — 2 passing + 1 skipped test.
examples/regression/crystal_parity/probe_spheropole_sad.py — reproducible CRYSTAL-CYC-0 probe.
docs/handover_bipole_v0_8_2026_05_18.md §2 — original problem statement + the deferral-then-implement decision trail.
~/.claude/projects/.../memory/reference_bipole_ext_spheropole_2026-05-18.md — line-by-line CRYSTAL source references.