BIPOLE Implementation Status

Last updated: 2026-09-28

Summary

Finite-mesh exchange accuracy

The Ewald/Madelung exchange correction does not eliminate finite-mesh error. For three-dimensional insulating systems, uniform refinement along all three mesh axes leaves an asymptotic error of order 1/N_k; see Xing and Lin, Phys. Rev. X 14, 011059 (2024), Introduction, items 1-2. Agreement with a reference using the same Ewald convention therefore does not establish agreement with a different exchange summation convention at the same mesh.

The H2/STO-3G diagnostic in issue #703 illustrates this limitation. Its orthorhombic cell is 6 x 7 x 8 angstrom, with fractional positions (0.25, 0.25, 0.25) and (0.25, 0.25, 0.3425) and conversion 0.5291772083 angstrom/bohr. Restricted HF BIPOLE at fixed 15-bohr pair support, requested 20-bohr nuclear support (resolved Ewald real-space radius 15 bohr) and 1e-12 Ewald/SR precision gives the following converged energies with full Gamma-centered meshes. SCF energy and gradient tolerances are 1e-11 Ha and 1e-8, respectively:

Mesh

Energy (Ha)

Difference from original CRYSTAL23 reference (microHa)

1 x 2 x 4

-1.1166643767403097

+92.5442

2 x 4 x 8

-1.1167422580671185

+14.6629

4 x 8 x 16

-1.1167550671549666

+1.8538

The original reference remains -1.1167569209709 Ha and the acceptance tolerance remains 1e-6 Ha: even the 512-point row fails that comparison. The last refinement reduces the offset by about eight, consistent with residual inverse-volume exchange error. Explicit matched basis primitives, tighter CRYSTAL integral thresholds, higher multipole order and disabling its bipolar integral approximation do not remove the discrepancy. A direct PySCF Ewald-exchange calculation at 64 points agrees with tight BIPOLE within 4.1e-12 Ha. These are bounded numerical controls, not general CRYSTAL23 accuracy or performance acceptance. Preserve mesh and zero-mode identities in comparisons; an extrapolated energy must not replace a measured row.

Issue #705 develops a Wigner-Seitz truncated exchange alternative following Sundararaman and Arias, Phys. Rev. B 87, 165122 (2013), Appendix A1, Eq. (A4). The internal _ws_exchange kernel implements the minimum-image FFT quadrature for orthogonal BvK supercells, including rotated cells. Skew cells, off-mesh reciprocal vectors, Nyquist violations and FFT allocations beyond the declared limit are refused. Its finite zero mode and long-range remainder subtract the existing erfc kernel consistently; an additional Madelung correction must not be applied.

The direct Python RHF/UHF/RKS/UKS drivers expose an experimental WS selection for full-range HF and global hybrids. Pure DFT needs no exchange truncation; screened exchange and chi routes refuse WS. Existing Ewald exchange remains the default. Select WS with explicit resolution controls:

from vibeqc._ws_exchange import WSExchangeOptions
from vibeqc.pbc_bipole import run_pbc_bipole_rhf

result = run_pbc_bipole_rhf(
    system, basis, kmesh, options,
    exchange_exxdiv="ws",
    ws_exchange=WSExchangeOptions(reciprocal_cutoff=6.0),
)

The cutoff is in inverse bohr and applies only to exchange. The example value is not a general convergence recommendation. Increase it independently of the FFT grid_scale and quadrature precision; the old Gaussian-only Ewald cutoff was insufficient for this kernel. The default maximum FFT grid is 8,000,000 points. Results record the resolved WS grid, BvK mesh, lattice, quadrature parameter, exchange cutoff and finite zero mode in result.resolved_domain.ws_exchange, included in the settings digest. result.exchange_exxdiv identifies the WS convention. The original J reciprocal cutoff and Ewald settings remain separate. The shared executed_exchange_q0 helper reads this executed BIPOLE field and reports wigner-seitz, even if the caller supplies an Ewald fallback. Comparison checks reject mixing this label with Ewald or strict-zero results. Matching the label alone does not establish matching kernel resolution, physical support or numerical convergence.

The public Python run_periodic_job accepts the same exchange_exxdiv="ws" and ws_exchange=WSExchangeOptions(...) pair with jk_method="bipole" for RHF/UHF and full-range global-hybrid RKS/UKS. Preflight checks the orthogonal 3D cell, full uniform mesh, FFT allocation bound and exchange cutoff without allocating the FFT. The cutoff must remain strictly below every FFT Nyquist face, min_i(pi * N_i / L_i), with BvK lengths L_i and FFT sizes N_i. Increase grid_scale or refine precision if the requested cutoff reaches this bound; admission alone does not establish numerical convergence. The .system manifest records requested resolution controls, the executed convention and the resolved kernel domain. Fixed-cell atomic geometry optimization (relax_atoms, run_periodic_job(optimize=True)) runs on the finite-difference force path, which differentiates the WS energy; OptimizeResult records the executed convention (exchange_convention / ws_exchange_identity). On H2/STO-3G these forces agree with the isolated molecule to 5e-5 Ha/bohr in a 12 bohr box and 2e-6 in a 16 bohr box, where the Ewald convention is still 1.9e-3 and 7.6e-4 away. Periodic run_neb accepts the same exchange_exxdiv="ws" / ws_exchange pair; each image energy and force equals the WS single point and finite-difference gradient at that geometry. Pure or screened KS exchange, variable-cell optimization, Hessians, TDDFT/TDA and READ restarts explicitly refuse this selection. Default Ewald calculations are unchanged. The molecular response kernel in the periodic TDDFT/TDA wrapper does not differentiate the WS SCF operator; both public preflight and direct response calls reject this combination before computing excitations.

Executable QVF inputs accept the same controls as a JSON record in job.spec.options. For example, pass the following options to write_pending_qvf(system, ..., method="RHF", basis="sto-3g", kpoints=[2, 1, 1]) and execute the container with vibeqc run:

options = {
    "jk_method": "bipole",
    "exchange_exxdiv": "ws",
    "ws_exchange": {"reciprocal_cutoff": 6.0, "precision": 1e-12,
                    "grid_scale": 2, "max_grid_points": 8_000_000},
}

The cutoff remains an example requiring convergence checks. It is required in the record; no Ewald-derived default is substituted. Unknown fields, nonfinite numbers, boolean integers and invalid resolutions are refused during typed reconstruction. This input support does not qualify READ restarts or turn a standalone result archive into a complete replay request. Re-executing an archive with WS result provenance now requires explicit periodic k-mesh counts and the BIPOLE/WS controls in its saved request. An incomplete standalone result is refused before archive mutation or computation, even with force=True. Original declarative WS inputs retain their replay path. This admission check prevents fallback to Ewald/Gamma; it does not certify that every physical support or convergence setting matches an earlier calculation.

Both the direct RHF driver and the public Python runner gave -1.1167569182733417 Ha on the matched 64-point H2 case, including exact-Fock convergence confirmation. The public run completed its output lifecycle: the manifest records the WS request, resolved kernel domain, matching settings digest and terminal convergence. Its difference from the original CRYSTAL23 reference -1.1167569209709 Ha is about 2.70 nHa, within the unchanged 1 microhartree tolerance. Agreement on this example does not certify convergence to nanohartree precision or general numerical or performance parity. Complete artifact/restart convention handling and broader qualification remain open. Both analytic and finite-difference derivative entry points explicitly refuse this experimental selection pending validation.

The separate Gamma H2/STO-3G comparison does not pass the same 1 microhartree target. In the 6 by 7 by 8 angstrom cell, with H fractional coordinates (0.25, 0.25, 0.25) and (0.25, 0.25, 0.3425), WS RHF gave -1.1167800790881228 Ha versus CRYSTAL23 -1.5449409530363 Ha. Refined-grid WS PBE0 RKS gave -1.154349912800185 Ha versus CRYSTAL23 XXXLGRID -1.2613898336484 Ha. Both programs converged, and independent PBE0 grid changes were below 1e-7 Ha, but the cross-code offsets remain approximately 0.428161 and 0.107040 Ha. Closed-shell UKS agrees with the native RKS result; that is not independent magnetic validation. The Gamma-reference CRYSTAL executable also reproduced the original 64-point input energy exactly at printed precision, despite its different binary checksum. These Gamma rows remain excluded from CRYSTAL23 parity acceptance. No energy shift or tolerance relaxation is applied, and the original 64-point reference and completed 0.18 release scope are unchanged.

The subsequent fixed-Gamma reference support ladder (#743) changes this diagnosis: SCF convergence did not establish a converged exchange series. At TOLINTEG 8/8/8/8/16, 10/10/10/10/20 and 12/12/12/12/24, CRYSTAL HF gives -1.4787026888051, -1.5449409530363 and -1.8605198026943 Ha; PBE0 XXXLGRID gives -1.2448302739884, -1.2613898336484 and -1.3402845352735 Ha. The spans are about 0.382 and 0.0955 Ha, far above the original 1 microhartree target. On a 2x2x2 mesh at the middle threshold, HF is -1.1167569259362 Ha and PBE0 is -1.1543271318230 Ha. These different-mesh controls diagnose sampling; they do not replace the Gamma references or qualify a matched Gamma comparison.

The mechanism is consistent with CRYSTAL23 manual Section 18.4, Eqs. (18.19)-(18.20): long-range exchange depends on the squared density matrix divided by image distance, and pseudo-overlap thresholds truncate its real-space series. Sundararaman and Arias, Section II A, Eqs. (10)-(12), DOI 10.1103/PhysRevB.87.165122, show that finite k sampling makes pair densities periodic on the BvK supercell; WS removes the replica interaction by truncating the potential. Thus matching the mesh and SCF stopping criterion does not match those finite-domain operators. The measured threshold sensitivity supports this diagnosis, without proving that every contribution has been isolated.

Independent native Gamma controls retained converged HF/PBE0 energies under split-parameter changes from 0.4 to 1.0 inverse bohr (changes below 1.4e-9 Ha). Increasing the WS reciprocal cutoff from 6 to 8 inverse bohr and grid scale from 2 to 3 changed the tested energies by less than 9.1e-8 Ha. This supports the native WS implementation on this example, not CRYSTAL Gamma parity. No native energy correction is justified by the unstable reference. Comparison records can now carry a measured support-energy span; a known unstable reference is refused rather than promoted by accidental energy agreement. A converged, explicitly matched finite-domain reference remains necessary before accepting this Gamma row.

That matched-operator check is now available for the bounded WS Gamma H2 case, using independent PySCF 2.14.0 FFTDF with exxdiv='vcut_ws'. The PySCF periodic SCF documentation identifies this as WS truncation, supported with FFTDF. Explicit STO-3G primitives, cell, positions, neutral closed shell and Gamma sampling match the native calculation above. The reference uses precision 1e-12, energy and gradient stopping criteria 1e-11 and 1e-8, and unpruned periodic Becke quadrature for PBE0.

Method

Independent WS reference (Ha)

Native refined WS result (Ha)

RHF

-1.1167801730605524

-1.1167801672837765

PBE0 RKS

-1.1543499192166702

-1.1543499348490986

Both differences are below the original 1e-6 Ha tolerance. Increasing the reference AO FFT mesh from 81x95x109 to 101x119x135 changes either energy by less than 9e-12 Ha. PBE0 quadrature refinement from 99 radial/1454 angular to 150 radial/2030 angular points changes the reference by less than 1.2e-10 Ha. Separately increasing only the reference WS quadrature grid from 120x140x160 to 180x210x240 changes HF/PBE0 by 3.77e-9/9.42e-10 Ha. This last control instruments only PySCF’s WS grid-size expression; the unmodified reference results above remain the pins. The original and instrumented kernel source, native-library hashes, inputs and outputs are retained with #743. No reference-program implementation is imported by the production code.

test_ws_gamma_h2_matches_independent_fft_reference exercises the native production drivers at WS cutoff 6 and grid scale 2 against these independent pins. The existing native support/resolution ladder and this independent reference qualify only this closed-shell Gamma WS example. No native energy defect was demonstrated. CRYSTAL’s original Gamma comparison remains failed and excluded because its exchange support is unstable; the independent WS acceptance is not CRYSTAL Gamma parity, general magnetic-state acceptance, or a change to completed 0.18.

The public-runner reference used the recorded existing native runtime with the tested Python source overlay. It is scoped numerical and manifest evidence, not qualification of a rebuilt native binary from current main. New .system manifests include run.bipole_resolved_domain_json, the canonical JSON string used for run.bipole_domain_settings_sha256. Hash its UTF-8 bytes to verify the digest, then parse the JSON to retain null values and nested option pairs exactly. Historical digest semantics are unchanged. The readable bipole_resolved_domain TOML table remains available, but its string representation of a nullable omitted-shell budget is not a lossless substitute. Older manifests without the JSON field still need their separately retained canonical record. This settings record alone does not supply a complete restart request or certify numerical accuracy.

The internal exchange-cache builder accepts this kernel independently of J. It uses a separate diagonal kernel, replaces nonzero-transfer kernels, and provides the matching finite zero-mode coefficient through the cache. Contractions require the declared full, uniformly weighted regular mesh; common shifts are allowed, while duplicate or reduced meshes are refused. The mesh-wide contraction retains the most recent ordered k-pair grouping for an exact point list and lattice. Later density or weight updates reuse that grouping; in-place changes to point coordinates or ordering invalidate it. The mapping retains O(Nk squared) incidence metadata separately from the existing byte-capped Fourier panels. Reciprocal canonicalization, physical screening and the order of exchange contractions are unchanged. Gamma and multi-k AO contraction tests verify energy-direction derivatives, batched/per-target agreement and preservation of the original J kernel. The direct drivers construct this cache on a separate exchange reciprocal envelope, leaving Coulomb J tables unchanged.

Restricted and unrestricted shared Fock builders obtain the zero mode from an explicit WS exchange cache at Gamma and multi-k. At Gamma, the exchange contraction uses the cache’s separate diagonal kernel rather than J’s kernel; without a WS cache the existing Ewald path is unchanged. The reported Madelung finite-size component is unset for WS, since that term does not belong to its convention. The total exchange energy still includes the finite WS zero mode and the short-range subtraction. Gamma and two-point closed-shell RHF/UHF consistency tests exercise this boundary; they do not constitute independent Gamma reference acceptance.

Production route

BIPOLE is vibe-qc’s own Ewald-J-split periodic SCF path. All four methods support Gamma and multi-k SCF and the exact finite-difference atomic-force path. Fixed-cell atomic relaxation is available through run_periodic_job(optimize=True). Variable-cell and coupled atom/cell optimization fail closed pending a certified strain convention and terminal same-geometry force/stress convergence check.

The exact Ewald-J-split path is the production route. Both multipole far-field implementations are unavailable from the SCF drivers. Passing use_multipole_far_field=True now fails before setup because the dormant quartet implementation does not preserve the exact contraction’s periodic translation domain. Independent measurement corroborates the retirement: with the quartet path exercised against exact ERIs, LiH/STO-3G shows 6.9e-2 Ha (14 bohr), 6.3e-2 Ha (18 bohr), and 5.9e-2 Ha (24 bohr) error that barely decays with separation, consistent with a domain mismatch rather than a truncation artefact.

RHF, UHF, RKS and UKS BIPOLE drivers now resolve physical AO-product separations for the one- and two-electron operators by default. Their cutoff_bohr controls pair separation; eri_interaction_cutoff_bohr controls the distance between product midpoints (zero uses cutoff_bohr). An explicit or automatically resolved BIPOLE SR extent supplies the latter distance. Exchange may therefore need output cells beyond the one-electron support. Physical mode uses the full direct Fock build: the legacy symmetry projector does not cover that wider exchange support. Explicit projection requests, periodic COSX and the diagnostic grid Hartree backend refuse this combination. The explicit historical_domain=True replay option restores the index-ball domain; the public runner exposes it as bipole_historical_domain=True. The native LatticeSumOptions default remains unchanged for other consumers. See the domain contract for the migration and its outstanding numerical acceptance matrix. The default change does not establish cutoff convergence or performance acceptance. For nonzero momentum-transfer exchange channels, physical support now uses the production Fourier kernel one supplied cell at a time, with the same AO-product mask. This removes the dense all-cell temporary and the repeated Python Fourier recursion for supported native shells. Cached channel tensors remain resident; this does not establish the full-node performance target. Finite-domain agreement alone does not certify infinite-lattice convergence; pair support, interaction reach and reciprocal resolution must all converge. The standalone continuation record distinguishes completed local checks from remaining acceptance work.

LiH and Si PBE cold starts at physical 20-bohr support agree with reference-density warm starts within 1e-8 Ha; 20 fixed-density grid cases and the default-screening comparisons are recorded there. Those earlier opt-in measurements do not certify the new production default or a native fine-grid SCF limit.

⚠ Absolute energies on ionic crystals, FIXED at Γ (all four drivers); multi-k corrected gauge available opt-in (2026-06-11). The root-caused defects (spheropole gauge double-count; Γ-locality projection dropping real cross-cell 1e/J terms; a formally divergent full-Coulomb direct-K series with the Bloch density) are fixed by the Ewald exchange split (option (b) redesign): K = K_SR(erfc ω) + K_LR(reciprocal, K≠0) + (ξ_M − π/(Vω²))·S·D·S, full-Bloch density, no spheropole term. Validated out-of-process vs PySCF GDF: exchange element-wise to 8e-4 / E_K to 0.3 mHa at the converged MgO density (ω-invariant); H₂ 12-bohr box converges to PySCF to ~1 µHa (RHF), ~0.05 mHa (RKS/SVWN + PBE, the projection removal also fixes the XC grid density), ~0.08 mHa (PBE0, the α_HF-scaled hybrid split path), +0.07 mHa (UHF triplet) and ~1 mHa (UKS). All four drivers default to the corrected gauge at BOTH Γ and multi-k (Phase-5 flip, 2026-06-13). The multi-k q = k−k′ ≠ 0 LR-exchange channels (option (b) Phase 3, 2026-06-11/12), per-(k,k′) momentum-transfer channels on q-shifted reciprocal meshes + the BvK-supercell Madelung correction (α_HF-scaled for hybrids, per spin for UHF/UKS), are validated on the H₂ box [2,1,1] vs out-of-process PySCF (RHF +0.12 mHa at cutoff 7, −0.002 at 12; RKS SVWN/PBE0 +0.04/+0.06; UHF +0.03; UKS k-shift 0.013) and to +0.0001 mHa/cell against the explicitly-doubled supercell at Γ (the supercell-unfolding identity, the rigorous internal proof). The multi-k corrected gauge needs a Monkhorst-Pack mesh. Physical-domain BIPOLE refuses an ad-hoc multi-k list that would fall back to the unsupported legacy gauge. Explicit legacy-gauge replay requires historical_domain=True together with use_exchange_ewald_split=False. Outstanding heavy confirmation: the MgO [2,2,2] c12 RHF+RKS parity run (examples/regression/bipole_parity/mgo_multik_parity.py) is staged for the queue, belt-and-suspenders on top of the supercell identity, fleet-blocked at the 2026-06-13 flip. Two practical caveats for ionic Γ work: (1) the corrected gauge contracts full Bloch folds, diffuse AO tails (e.g. STO-3G Mg 3sp) need fold-converged cutoffs (MgO/STO-3G: ≥12-16 bohr; the driver fails before SCF above 1e-2 S(Γ)-fold drift and reports intermediate truncation until the quantitative 1e-4 target is met); (2) minimal-basis ionic Γ SCF has multiple genuine solutions (PySCF itself lands on different totals from different starts on MgO/STO-3G Γ), check the basin against a reference density for absolute claims. Forensics + tables: the measurements above; examples/regression/bipole_parity/mgo_exchange_split_probe.py.

Production forces use finite differences of a converged driver energy. compute_bipole_gradient_fd differentiates the actual driver energy, so it follows whichever exchange gauge the SCF used (corrected or legacy); bipole_optimize, run_periodic_job(optimize=True), and periodic NEB default to it. Asymmetric legacy-gauge HF is not a supported converged-SCF route. The BeH2 and BeH legacy-gauge regression rows failed the exact terminal stationarity check even after their overlap folds were converged; those capability claims are retired (#66). Increasing the cutoff does not make a nonstationary density valid for a finite-difference gradient. The analytic BIPOLE gradient began with the LEGACY Gamma-local gauge. Its maintained scope is narrower than the corrected-gauge paths below. Corrected-gauge analytic-gradient prototypes landed 2026-06-17 for all four methods (RHF/UHF/RKS/UKS) at both Γ and multi-k (G3 row below). They remain a gated preview; the legacy-gauge analytic preview remains available and covers:

  • RHF/UHF Gamma and multi-k, maintained symmetric controls; the asymmetric legacy-gauge HF rows no longer certify this preview.

  • RKS/UKS Gamma-local, maintained LDA asymmetric cells with XC Pulay, moving-grid correction, and KS Bloch-CPHF.

  • RKS/UKS multi-k, diagonal-Z + corrected-W (the v_bg/spheropole W convention) + J^LR + XC Pulay; warns. KS CPHF gated to Gamma.

  • Finite-temperature/fractional-occupation analytic gradients are now landed in the corrected gauge (Mermin free-energy form; FD-validated RKS Gamma ~2e-8 / multi-k ~4e-9). The legacy gauge above still rejects them. The padded M5 Gamma composition also fails closed for fractional occupations until its finite-domain J/K map has a variational density-transpose completion.

  • Meta-GGA tau-Pulay is landed (un-gated 2026-06-20; validated to 1e-9 against an independent reimplementation; see G3).

  • The production M5 padded erfc domain is now differentiated for corrected-gauge RHF, UHF, and integer-occupation RKS/UKS Gamma with radial or pair-resolved SYM3b output. Pure-RKS radial parity on the post-exact-FT production SR+LR surface is 1.14e-7 Ha/bohr; pair-resolved UHF is 4.1e-6 Ha/bohr and pair-resolved LDA/PBE0 RKS pass below 1e-5 Ha/bohr. The radial finite-domain density adjoint closes the fully occupied triplet-H2 UHF discriminator from 2.906e-3 to below 1e-6 Ha/bohr and brings padded LDA/PBE0 UKS below 1e-5 Ha/bohr. The pair-mask transpose is also complete for UHF: its off-boundary saturated-spin discriminator improves from 2.31e-3 to 9.85e-8 Ha/bohr. Pair-resolved LDA/PBE0 UKS pass at 1.11e-7 Ha/bohr. Corrected-gauge integer-occupation RHF also supports the padded radial domain at multi-k: the finite-domain density adjoint and full real-linear k-coupled response pass the asymmetric BeH2 [2,1,1] production-energy finite difference at the fold-reliable 13-bohr cutoff. The RHF response now solves the weighted transpose equation with consistent complex rotations; the former negative-curvature expected failure is an ordinary passing regression at its unchanged 5e-5-Ha/bohr tolerance. See the continuation evidence for the operator diagnostic and validation scope. Corrected-gauge integer-occupation UHF also supports padded radial multi-k: preserving the full BvK support of D_alpha + D_beta plus the spin-coupled response closes asymmetric BeH [2,1,1] parity to 6.83e-8 Ha/bohr. Padded multi-k KS, pair-resolved, and fractional routes remain fail-closed. Fixed-density J/K and method-level full-SCF finite differences pin the supported slices; FD remains the production force path.

  • The legacy multi-k KS-CPHF remains a ratified deferral (2026-06-20), legacy-gauge-only, since the corrected-gauge multi-k KS path is variational and needs no CPHF.

  • Historical-domain space-group symmetry: attached symmetry auto-enables both the bit-exact SYM2c one-electron S/T reduction and, on the corrected Ewald exchange split, the SYM3b pair-resolved Fock reduction. The M5 Fock default is safe because every erfc SR build now uses a precision-derived padded ket-image ball with charge-pair Schwarz screening that includes angular and contraction tails. The full pair-resolved output, density, and internal domains are group-invariant: reduce equals enforce end to end, and J-block orbit asymmetry is 6.3e-15 on MgO, 5.4e-15 on diamond, and 4.4e-15 on Si at cutoff 6. Pass use_fock_symmetry_reduce=False to disable the automatic reduction; explicit enforcement-only diagnostics use use_fock_symmetry=True, use_fock_symmetry_reduce=False. Physical-domain AUTO uses the full direct build; explicit Fock symmetry requests are refused until that projection covers the wider exchange support.

  • The M5 erfc SR default is sr_image_precision=1e-6. It derives an absolute radius as cutoff_bohr + bipole_sr_image_extent(...) and enables LatticeSumOptions.sr_range_screening. Pass sr_image_precision=None to disable automatic image padding while retaining the selected physical or historical pair policy. On a historical erfc SR route this also suppresses automatic SYM3b reduction unless it is explicitly requested. An explicit sr_image_extent_bohr remains the absolute-radius oracle and takes precedence. HSE screened exchange uses the same padded internal domain; the longer-ranged of its erfc K kernel and the Ewald erfc J kernel sets the precision-derived radius. Under physical pair support that radius is also the exchange pair-interaction reach, so a Gamma-only run includes the screened exchange of its periodic images out to the resolved precision (1.7e-4 Ha on H2 in a 20-bohr box); the multi-k Ewald drivers reach the same image shells only once cutoff_bohr passes them. Periodic run_neb and relax_atoms expose and forward the same sr_image_precision=1e-6 production default. Their None setting explicitly disables automatic padding. run_neb retains physical pair support; relax_atoms also accepts historical_domain=True for explicit historical replay. All production force paths default to the exact FD gradient compute_bipole_gradient_fd.

Completed

Milestone

Description

M1-M12

Foundation, spheropole diagnosis, multipole infrastructure

M13

Multi-k via kpoints kwarg

M14

RKS/UKS/UHF promoted to user API

M15

Multipole near/far split fix

G2

C++ spheropole kernel (all cases through d-d)

G3

FD gradient (all 4 methods, exact) + analytic research preview

M16

2D/1D gradient support documented

M17

optimize=True in run_periodic_job

M18

Force-virial diagnostic; production cell optimization fails closed

User API

from vibeqc.periodic_runner import run_periodic_job

# SCF + atomic relaxation
result = run_periodic_job(system, basis, method="RHF", jk_method="bipole",
                          optimize=True, optimize_max_iter=30)

# Standalone force virial diagnostic (not the periodic stress)
from vibeqc.bipole_gradient import compute_stress_tensor
virial = compute_stress_tensor(system, gradient)

Remaining Gaps

#

Task

Priority

Status

G1

Cell-level multipole far-pair

HIGH

RETIRED / DRIVER-UNREACHABLE. The G1 builder is retained only as a research artifact. No direct or public driver can activate it; explicit multipole requests fail before SCF.

G2

Spheropole higher-l terms

HIGH

LANDED, +4.112 Ha vs CRYSTAL +4.119 Ha

G3

Production force/stress coverage

MEDIUM

Finite-difference atomic forces differentiate the supported exact driver energy. Variable-cell optimization and the standalone far-field gradient prototypes are not certified and fail closed.

G4

Space-group symmetry for one- and two-electron integrals

MEDIUM

Historical-domain SYM2c/SYM3b support is retained. Physical-domain BIPOLE uses the full direct Fock build; explicit legacy Fock projection requests fail closed because the projector does not cover the wider exchange support.

G5

Exact-zone / fold-range decoupling (BIPOLE-EXACT-ZONE-UNBOUNDED)

HIGH

Increment 1 LANDED 2026-08-06 (exact_zone_bohr, now wired through RHF, RKS, UHF and UKS, opt-in): the exact erfc output zone is bounded independently of the operator/fold cutoff; far operator cells carry J_LR + background (the Ewald far-field model – the Dovesi 1983 Eq. (24a) bielectronic/monoelectronic partition in the split gauge; NOT the retired G1 cell-level branch, which stays retired). Measured fixed-density tails: MgO/STO-3G 2.7e-3 Ha @ zone 10, 3.2e-4 @ 12; LiH-class diffuse bases keep mHa tails to ~14 bohr – per-case ladder knob. Auto SYM3b reduction yields to a zone request; explicit reduction and analytic gradients fail closed. Unrestricted SCF/provenance and spin-schedule forwarding regressions cover the existing UHF/UKS wiring. Remaining: automatic fold-cutoff selection coupled to the zone, then the quartet-level screened bipolar far field (pair-resolved, S92 Eq. (92) Boys-function radial ladder) as increment 2 – the scaling fix proper. Any future quartet far-field increment must first preserve the common translation sum in Pisani-Dovesi-Roetti 1988 Eq. II.4.10; this exact-zone result is not evidence for that route.

G6

Quartet-level bipolar far-field

P0

FAIL-CLOSED 2026-08-13. The driver gate was dormant, and enabling it directly would be unsafe: the exact Fock contraction spans three lattice translations and contracts density at their relative translation, while the dispatch, skip mask, and add-back span only two. Additional blockers include the penetration rule, screened kernel, total-order ceiling, Coulomb/exchange composition, symmetry reconstruction, and gradients. All four drivers default to False; explicit True raises before SCF. Independent exact-vs-far-field measurement corroborates the retirement: LiH/STO-3G shows 6.9e-2 Ha (14 bohr), 6.3e-2 Ha (18 bohr), 5.9e-2 Ha (24 bohr), barely decaying with separation – a domain mismatch, not truncation. The restricted Fock builder’s dormant G6 branch was also made internally consistent (skip-masked near-K for both the Fock and the exchange energy) so any future re-wiring starts from a coherent reference.

G7

Penetration dispatch and contractor prototypes

HIGH

NOT ROUTE-CERTIFIED. Isolated Python/C++ parity tests cover only internal representations. They do not establish a literature basis for the implemented classifier, do not validate the published penetration criterion or the exact/add-back periodic domain. Order zero is currently overloaded as an exact sentinel, the classifier omits the second pair width, and screened caches and symmetry reconstruction need independent equation-level oracles before reuse.