BIPOLE handover - current state after integration sprint (2026-05-20)

This handover documents the post-dense-k-sign-off integration work that added RKS/UKS DFT drivers, UHF feature parity (ODA/MOM/level-shift), periodic_runner dispatch, and the multipole far-field J builder. The older handover docs/handover_bipole_v0_8_2026_05_18.md covers the dense-k sign-off history.

Current headline

All four BIPOLE drivers (RHF, UHF, RKS, UKS) are now integrated into the periodic_runner dispatch under jk_method="bipole":

  • run_pbc_bipole_rhf — closed-shell RHF (dense-k parity signed off)

  • run_pbc_bipole_uhf — open-shell UHF (ODA/MOM/level-shift parity to RHF)

  • run_pbc_bipole_rks — closed-shell RKS DFT (libxc V_xc + hybrids)

  • run_pbc_bipole_uks — open-shell UKS DFT (spin-polarised V_xc + hybrids)

All share the CRYSTAL-gauge Ewald J-split F²e build.

The multipole far-field J builder (bipole_fock_multipole.py) exists as a standalone module with diagnostic hooks in the RHF driver. It is not yet wired for computational acceleration.

Remaining work:

  1. Wire multipole far-pair into Fock build for acceleration,

  2. Complete compute_ext_el_spheropole higher-l terms (~17%→100%),

  3. BIPOLE gradient.

Commits since the older handover

The useful BIPOLE-specific commits at the top of main are:

  • f3d3dcc Implement CRYSTAL-gauge periodic BIPOLE SCF

    • Analytic 3D Ewald V_ne by shifted AO-pair Fourier transforms.

    • Native build_jk_2e_real_space component build.

    • Multi-k Ewald-J path with shifted-nu AO-pair FT, k-space rho_hat(K), real-space J_LR(g) blocks, and electron-electron background.

    • CRYSTAL ENECYCLE diagnostics and MgO dense-k parity tooling.

  • b8db49f Harden dense-k BIPOLE parity scripts

    • Diamond and Si dense-k parity scripts updated to the passing cutoff-14 setup.

    • Large k-mesh overlap logging condensed to one k-mesh summary unless a non-ok severity is detected.

  • 5ada267 Default BIPOLE RHF to CRYSTAL gauge

    • run_pbc_bipole_rhf(..., use_ewald_j_split=None) now defaults to the CRYSTAL-gauge Ewald-J path for 3D.

    • The shared CRYSTAL parity runner auto-disables k symmetry for Ewald-J multi-k runs, avoiding under-sampled IBZ input.

    • The original MgO SHRINK 8 parity script now targets the passing dense-k setup instead of the stale direct-only branch.

  • 83a5da8 Guard BIPOLE ODA against gauge mismatch

    • Temporary safety guard to prevent ODA from using a direct-only trial Fock on the Ewald-J path.

  • 8c1d3ee Support ODA on BIPOLE Ewald-J path

    • Replaced the guard with real support: the ODA trial Fock uses the same Ewald-J composition as the main Fock.

    • Tracks whether the current density is exactly represented by the current C(k) orbitals. Damped or ODA-mixed densities fall back to real-space density blocks for J_LR density transforms instead of reusing stale orbitals.

  • a947934 Expand IBZ meshes for BIPOLE Ewald-J

    • IBZ-reduced Monkhorst-Pack meshes carrying ir_mapping metadata are expanded internally to the full uniform mesh for the current Ewald-J density transform.

  • current working slice (after a947934, not yet part of the older dense-k sign-off):

    • python/vibeqc/pbc_bipole_uhf.py adds run_pbc_bipole_uhf, a spin-unrestricted CRYSTAL-gauge Ewald-J BIPOLE driver.

    • tests/test_pbc_bipole_uhf.py covers public export, closed-shell UHF -> RHF equivalence at gamma and full 2x2x2 k-mesh, and an open-shell doublet smoke test.

Dense-k CRYSTAL14 sign-off

All numbers below are RHF/STO-3G, SHRINK 8 8, full k-mesh (use_symmetry=False), cutoff 14 bohr, analytic Ewald V_ne, Ewald-J F^2e, DIIS starting at iter 2.

Demo

vibe-qc E_total (Ha/cell)

CRYSTAL14 reference

Delta

MgO

-271.2177748509

-271.2181437248

+0.369 mHa

diamond

-74.8771393842

-74.8769943962

-0.145 mHa

silicon

-571.3214715798

-571.3208122082

-0.659 mHa

Component diagnostics also look sane:

  • Diamond two-cycle ENECYCLE diagnostic:

    • CYC 0 total delta about -0.012 mHa.

    • CYC 1 total delta about -0.332 mHa.

  • Silicon two-cycle ENECYCLE diagnostic:

    • CYC 0 total delta about +0.233 mHa.

    • CYC 1 total delta about -1.053 mHa.

The MgO dense-k parity script is examples/regression/crystal_parity/parity_mgo_primitive_shrink88.py. There are also explicit Ewald-J scripts for MgO SHRINK 2/8 and gamma, plus diamond and silicon gamma/dense-k scripts.

Verification already run

Recent focused test slice after ODA support:

.venv/bin/python -m pytest \
  tests/test_pbc_bipole_ewald_split_integration.py \
  tests/test_pbc_bipole_multik_ewald_split.py \
  tests/test_runner_crystal.py \
  -m 'not slow' -q

Result: 52 passed, 2 deselected in 519.86s.

Additional focused UHF slice:

.venv/bin/python -m pytest tests/test_pbc_bipole_uhf.py -q

Result: 4 passed in 0.19s.

Post-UHF focused BIPOLE slice:

.venv/bin/python -m pytest \
  tests/test_pbc_bipole_uhf.py \
  tests/test_pbc_bipole_ewald_split_integration.py \
  tests/test_pbc_bipole_multik_ewald_split.py \
  -m 'not slow' -q

Result: 25 passed, 2 deselected in 394.38s.

Manual NiO smoke after the UHF scaffold used build_nio_sto3g(), a gamma mesh, cutoff 5 bohr, nuclear cutoff 10 bohr, HCORE guess, max_iter=1, and ewald_precision=1e-6. Result: E_total = -1559.9365798653657 Ha, S^2 = 2.00000000006 for the triplet target (S(S+1)=2). This only proves the path runs; it is not a CRYSTAL parity result.

CRYSTAL14 NiO STO-3G reference is now generated at examples/regression/crystal_parity/crystal_demos/nio_sto3g.out:

  • SHRINK 8 8, UHF, SPINLOCK 2 15, LEVSHIFT 3 1, FMIXING 30.

  • Final energy: -1563.9384656162 Ha in 24 cycles.

  • Final ENECYCLE components:

    • EXT EL-POLE -1135.0520786544 Ha

    • EXT EL-SPHEROPOLE +9.7016419049 Ha

    • BIELET ZONE E-E +1493.8680799096 Ha

    • TOTAL E-E +368.5176431601 Ha

    • TOTAL E-N + N-E -3218.2714523824 Ha

    • TOTAL N-N -278.0698772643 Ha

    • KINETIC ENERGY +1563.8852208704 Ha

  • GCALCO reports RMAX 58.88485 BOHR and 6999 direct lattice vectors. This is the important lesson: NiO is not a small-cutoff Ewald-J scaffold case. It needs the native SDR multipole far-pair branch for practical parity.

After this observation, run_pbc_bipole_uhf was corrected to mirror CRYSTAL’s UHF SAD convention for even-electron cells: start CYC0 from a spin-averaged total SAD (D_alpha = D_beta = D_SAD/2) and impose the requested spin state through alpha/beta occupations after the first diagonalisation. This matches CRYSTAL’s NiO SUMMED SPIN DENSITY 0.0 at CYC0.

The earlier post-rebase focused slice before ODA support was 50 passed, 2 deselected; after adding default-gauge and ODA coverage that slice became 52 tests. The new tests assert:

  • 3D BIPOLE defaults to the CRYSTAL-gauge Ewald-J path.

  • ODA on Ewald-J runs and its trial Fock exposes the same J/K component decomposition.

Syntax/hygiene checks run on the touched files:

.venv/bin/python -m py_compile \
  python/vibeqc/pbc_bipole.py \
  tests/test_pbc_bipole_ewald_split_integration.py
git diff --check

Important implementation details

Shared Ewald state

The driver intentionally uses one alpha for:

  • analytic V_ne reciprocal/background,

  • E_nn, and

  • J_LR.

This mirrors CRYSTAL’s COMMON/VRSMAD/ behavior. Do not reintroduce independent alpha choices in separate components.

Fock build helpers

python/vibeqc/pbc_bipole.py now has an internal _PBCBipoleFockBuild bundle and an inner _build_fock_for_density helper inside run_pbc_bipole_rhf. Both the main SCF Fock and the ODA trial Fock route through this helper. This is deliberate: it prevents future gauge drift between “normal” and “auxiliary” Fock builds.

python/vibeqc/pbc_bipole_uhf.py has the analogous _PBCBipoleUHFFockBuild. The UHF two-electron operator is:

  • J_total = J_SR(D_alpha + D_beta) + J_LR(D_alpha + D_beta) + V_bg*S

  • F2e_alpha = J_total - K(D_alpha)

  • F2e_beta = J_total - K(D_beta)

The UHF two-electron energy is evaluated as 0.5 * (tr[D_alpha F2e_alpha] + tr[D_beta F2e_beta]).

Density provenance for J_LR

The reciprocal long-range J path has two valid density routes:

  • If the density is exactly represented by the current orbitals, compute rho_hat(K) from k-space D(k).

  • If the density is local SAD, fixed-damped, or ODA-mixed, compute rho_hat(K) from the real-space density blocks.

The density_from_c_per_k flag in pbc_bipole.py tracks this. Be careful when changing damping/ODA/DIIS behavior: using stale C(k) for an ODA-mixed density is a composition bug.

Multi-k symmetry handling

The Ewald-J path needs full-mesh sampling for the Bloch-summed AO-pair FT. Symmetry-reduced IBZ meshes have non-uniform weights and would under-sample that transform if consumed directly.

The driver now accepts IBZ-reduced Monkhorst-Pack meshes when ir_mapping metadata is present and expands them internally to the corresponding full uniform mesh. This is a correctness/usability bridge, not a performance win. True IBZ operator/orbit expansion can come later. Explicit non-uniform custom k-meshes without orbit metadata are still rejected on the Ewald-J path.

Known limitations

  • 3D RHF has dense-k parity for MgO / diamond / Si STO-3G. 3D UHF now has a working Ewald-J scaffold, but NiO parity is still unsigned.

  • dim=1 / dim=2 BIPOLE parity remains blocked.

  • Metallic smearing in BIPOLE remains blocked; Be bulk remains blocked.

  • Inline CRYSTAL basis blocks remain blocked; several extended-basis demos are still stubs.

  • True IBZ acceleration is not implemented; IBZ meshes are expanded to full meshes internally on the Ewald-J path.

  • The native Saunders-Dovesi-Roetti far-pair multipole dispatch is not wired into the C++ Fock kernel. The current Ewald-J path is a CRYSTAL-gauge parity scaffold with good dense-k energies but not the final production-performance BIPOLE implementation.

  • EXT EL-SPHEROPOLE is not a standalone kernel. Its contribution is effectively captured by the current reciprocal-space construction for the dense-k sign-off cases, but a line-for-line CRYSTAL decomposition would still need additional work.

Useful commands

Full dense-k MgO parity:

.venv/bin/python -u examples/regression/crystal_parity/parity_mgo_primitive_shrink88.py

Diamond and silicon:

.venv/bin/python -u examples/regression/crystal_parity/parity_diamond_sto3g.py
.venv/bin/python -u examples/regression/crystal_parity/parity_sibulk_sto3g.py

Two-cycle ENECYCLE diagnostics against CRYSTAL:

.venv/bin/python -u examples/regression/crystal_parity/diagnose_bipole_enecycle.py \
  diamond --shrink 8 --cutoff 14 --maxcycle 2 --cycles 0,1 \
  --ewald-precision 1e-6 --crystal-exe ~/bin/crystal

.venv/bin/python -u examples/regression/crystal_parity/diagnose_bipole_enecycle.py \
  si --shrink 8 --cutoff 14 --maxcycle 2 --cycles 0,1 \
  --ewald-precision 1e-6 --crystal-exe ~/bin/crystal

Focused test slice:

.venv/bin/python -m pytest \
  tests/test_pbc_bipole_ewald_split_integration.py \
  tests/test_pbc_bipole_multik_ewald_split.py \
  tests/test_runner_crystal.py \
  -m 'not slow' -q