BIPOLE: Ewald-J-split periodic HF and DFT

Feature status

The BIPOLE drivers (run_pbc_bipole_rhf, _uhf, _rks, _uks) form vibe-qc’s BIPOLE periodic workstream. All four methods support multi-k. The Saunders 1992 quartet-level bipolar far-field automatically activates for large cells (cutoff > 8 bohr) when the Ewald exchange split is active. Production gradients and structure optimization (atomic + cell relaxation) use the finite-difference BIPOLE force path by default. Analytic gradients are available as a gated research-preview path on the maintained regression surface.

For routine calculations, prefer:

• run_rhf_periodic_gamma_gdf for Γ-only hybrid / pure DFT with native Gaussian density fitting.

• run_rhf_periodic_multi_k_ewald3d for multi-k with 3D Ewald for the long-range Coulomb piece.

What BIPOLE is

CRYSTAL splits the periodic Fock build into two halves and keeps them in separate gauges:

• One-electron Coulomb (V_ne) and the nuclear-nuclear energy (E_nn), 3D Ewald (point-charge ↔ Gaussian-pair doesn’t decay on the nucleus side, so direct truncation diverges on charged-nucleus crystals like MgO). Both share a single Ewald α (a single shared Ewald state).

• Two-electron Coulomb + exchange (F^{2e} = J + K), CRYSTAL’s native code uses direct-space BIPOLE screening and multipole-far-pair replacement. The current vibe-qc parity path implements the same electrostatic gauge as J_SR(ω) + J_LR(ω) + V_bg·S - 1/2 K_full, where ω is the same Ewald α used by V_ne / E_nn and V_bg = -π N_e/(α²V). This replaces the older broken run_rhf_periodic_multi_k_ewald3d composition, whose long-range J did not share CRYSTAL’s gauge.

What today’s implementation does

vibeqc.run_pbc_bipole_rhf ships:

• The full real-space one-electron pipeline at Ewald gauge: \(S(g)\), \(T(g)\), \(V_{ne}(g)\) at opts.lattice_opts.cutoff_bohr, Bloch-summed to per-k \(S(k)\) and \(H_{core}(k)\), canonical orthogonalisation \(X(k)\). The singular erfc part of \(V_{ne}\) is analytic via libint, and the smooth reciprocal/background part is analytic by shifted AO-pair Fourier transforms. A tightened unpruned Lebedev grid is still available as an explicit diagnostic fallback via v_ne_grid_options.

• A shared-gauge Ewald-J two-electron build by default for 3D systems (use_ewald_j_split=None, or True explicitly): short-range J from direct erfc-screened ERIs, long-range J from reciprocal AO-pair Fourier transforms, the matching electron background potential, and full direct-space exchange.

• The standard SCF inner loop: Bloch-sum, direct real-space energy evaluation (Σ_g D(g)M(g)), optional _MultiKPulayDIIS, optional LEVSHIFT, diagonalisation, optional MOM occupied-subspace reorder, and density rebuild. ODA mixing builds its trial Fock with the same Ewald-J composition; for ODA-mixed densities that no longer have an exact C(k) representation, the long-range-J density transform falls back to the real-space density blocks.

• The matching nuclear-nuclear Ewald sum so \(E_{nn}\) shares the α used by \(V_{ne}\).

vibeqc.run_pbc_bipole_uhf uses the same one-electron and Ewald-J machinery with unrestricted spin densities:

• \(D_{tot}(g) = D_\alpha(g) + D_\beta(g)\) drives the Hartree operator.

• \(F^{2e}_\alpha(g) = J_{tot}(g) - K[D_\alpha](g)\) and \(F^{2e}_\beta(g) = J_{tot}(g) - K[D_\beta](g)\).

• The UHF energy is evaluated as \(E_{2e} = \frac{1}{2}\mathrm{tr}[D_\alpha F^{2e}_\alpha] + \frac{1}{2}\mathrm{tr}[D_\beta F^{2e}_\beta]\) in the same real-space lattice contraction convention as RHF.

• Even-electron SAD guesses are spin-averaged for CRYSTAL CYC0 parity: \(D_\alpha = D_\beta = D_{SAD}/2\). The requested multiplicity enters through the alpha/beta occupation counts after diagonalisation.

What is still gated:

• A certified/default Saunders-Dovesi-Roetti 1992 multipole-far-pair branch, at large distances CRYSTAL uses truncated multipole expansions rather than exact direct ERIs. vibe-qc has an opt-in experimental branch (use_multipole_far_field=True), but its interaction-tensor normalization and end-to-end accuracy are still being certified. The exact Ewald-J split captures the same long-range gauge and remains the production path.

• True IBZ orbit expansion for the multi-k long-range-J density transform. Until that lands, symmetry-reduced Monkhorst-Pack inputs are accepted by internally expanding them back to the full mesh for Ewald-J correctness. This is a usability bridge, not a performance speedup.

The driver is therefore a parity harness plus an increasingly useful SCF implementation. The first-cycle CRYSTAL comparison is now a high-signal regression test; final-energy sign-off still belongs in the CRYSTAL parity scripts.

Why ship it at all?

Two reasons:

1. CRYSTAL parity surface. The 15 demo geometries committed alongside the driver (commit c1e00ff, tests/demos/

   • LiH / NaCl / MgO / Al2O3 / C-diamond / Si-diamond / ZnO / TiO2 / SiO2 / Ne-fcc) and the matching CRYSTAL14 reference inputs let the periodic-dev chat track multipole-branch progress against ground-truth CRYSTAL energies. Each demo has a parity_<system>.py runner that produces a side-by-side comparison.

2. Multipole infrastructure. The Phase 1-4b commits (d1cc402 … 7de6431) added the shell-pair Cartesian multipole moments, the multipole-multipole interaction tensor, the IDIPC geometric dispatch per quartet, and the cell-level multipole moments from density. These are the prerequisites for certifying the SDR multipole branch and for the bipole_lattice_self_energy analytic dipole-dipole self-image correction. They are tested in isolation against analytical results; the exact Ewald-J split remains the production SCF route while the opt-in multipole branch is certified.

BIPOLE phase status

Phase	Status	What
1	landed	Shell-pair Cartesian multipole moments (C++ via libint)
2	landed	Multipole-multipole Coulomb interaction tensor
3	landed	IDIPC geometric dispatch per quartet
4a	landed	Cell-level multipole moments from density
4b	landed	Neutral-cell spherical-sample dipole surface diagnostic
4b proper	gated	CRYSTAL EXT EL-POLE shell-multipole/penetration decomposition; the former complete-AO-density reciprocal helper fails closed
5 (Γ-only)	landed (2026-05-18)	Γ-only end-to-end SCF with use_ewald_j_split=True
5 + DIIS	landed (2026-05-18 late)	DIIS-compatible via D-consistent energy/error formulation
5 (multi-k)	landed (dense-k sign-off, 2026-05-20)	Per-k F_J^LR(k), MgO/diamond/Si STO-3G within 1 mHa of CRYSTAL14
5K (RKS)	landed (2026-05-20)	BIPOLE RKS driver with libxc V_xc + hybrid support
5U	landed (2026-05-20)	UHF BIPOLE with ODA/MOM/level-shift parity to RHF
5UK	landed (2026-05-20)	BIPOLE UKS driver with spin-polarised V_xc
5b	landed (2026-05-20)	Multipole far-field J builder (bipole_fock_multipole.py); opt-in research-preview diagnostics
6a	gated	Default-on multipole far-field certification; exact Ewald-J remains production
6b	landed	compute_ext_el_spheropole higher-l terms
7	roadmap	CRYSTAL14 numerical-parity sign-off on the 15 demo set

Phase 1-5 plus the exact Ewald-J split give the production BIPOLE energy route; it does not call the gated EXT EL-POLE helper. The remaining Phase 4b-proper and Phase 6-7 work implements CRYSTAL’s term decomposition, certifies the optional multipole far-field branch, and widens the external parity surface.

Dense-k STO-3G sign-off results at cutoff 14 bohr:

Demo	vibe-qc E_total (Ha/cell)	Δ vs CRYSTAL14
MgO SHRINK 8	-271.2177748509	+0.369 mHa
diamond SHRINK 8	-74.8771393842	-0.145 mHa
silicon SHRINK 8	-571.3214715798	-0.659 mHa

Quick start: SCF with the Ewald J-split

The Ewald-J split is automatic for 3D parity work. Full multi-k meshes remain fastest to reason about, but IBZ-reduced Monkhorst-Pack inputs are accepted and expanded internally to the full mesh until true IBZ orbit expansion is wired. Pass use_ewald_j_split=False only when intentionally reproducing the legacy direct-only diagnostic branch.

import vibeqc as vq
from vibeqc import InitialGuess, monkhorst_pack
from vibeqc._vibeqc_core import PeriodicRHFOptions
from vibeqc.pbc_bipole import run_pbc_bipole_rhf

system, basis = build_mgo_sto3g()    # 15 demo geometries in crystal_demos/

# Γ-only:
kmesh = monkhorst_pack(system, [1, 1, 1])
# Multi-k:
# kmesh = monkhorst_pack(system, [2, 2, 2], use_symmetry=False)
# kmesh = monkhorst_pack(system, [2, 2, 2], use_symmetry=True)
# The second form is expanded internally to the full mesh for Ewald-J.

opts = PeriodicRHFOptions()
opts.lattice_opts.cutoff_bohr = 14.0
opts.lattice_opts.nuclear_cutoff_bohr = 14.0
opts.initial_guess = InitialGuess.SAD
opts.use_diis = True
opts.diis_start_iter = 2
opts.damping = 0.0
opts.max_iter = 15

result = run_pbc_bipole_rhf(
    system, basis, kmesh, opts,
    # use_ewald_j_split=None defaults to the CRYSTAL-equivalent
    # F^{2e} build for 3D systems.
    ewald_precision=1e-8,
)
print(f"E_total = {result.energy:+.6f} Ha (CRYSTAL: -271.218)")

For open-shell systems, use the UHF scaffold:

from vibeqc.pbc_bipole_uhf import run_pbc_bipole_uhf

result = run_pbc_bipole_uhf(
    system, basis, kmesh, opts,
    ewald_precision=1e-6,
)

Constraints (2026-05-20):

• True IBZ acceleration is not implemented. The k-space ρ̂(K) formula does not yet expand IBZ symmetry orbits directly. If an IBZ-reduced Monkhorst-Pack mesh carries ir_mapping metadata, the driver expands it to the corresponding full mesh internally. Explicit non-uniform custom k-meshes without that metadata are still rejected on the Ewald-J path.

• Ewald-J requires dim = 3. With the auto default, dim < 3 runs stay on the legacy direct-only diagnostic branch; passing use_ewald_j_split=True for 2D / 1D raises.

Restricted exact bielectronic zone (exact_zone_bohr)

On diffuse bases the S(k) fold gate legitimately demands a large cutoff_bohr (LiH/STO-3G needs ~37 bohr for a 1e-4 fold drift), but the fold requirement is one-electron and cheap; historically the expensive exact erfc J/K traversal grew with it. exact_zone_bohr (all four RHF/RKS/UHF/UKS drivers, opt-in; screened hybrids and ROHF/ROKS fail closed) bounds the exact-ERI output zone independently: operator cells beyond the zone carry the reciprocal J_LR channel plus the neutralising background alone – the Ewald far-field model of the density. This is the bielectronic-zone / monoelectronic-zone partition of the CRYSTAL Coulomb machinery (Dovesi, Pisani, Roetti, Saunders, Phys. Rev. B 28, 5781 (1983), Eq. (24a); Pisani, Dovesi, Roetti, Lecture Notes in Chemistry 48 (1988), Sec. II.4b/II.4d) expressed in the Ewald-split gauge.

result = run_pbc_bipole_rhf(
    system, basis, kmesh, opts,      # opts.lattice_opts.cutoff_bohr = 20.0
    exact_zone_bohr=12.0,            # exact erfc quartets only to 12 bohr
    ewald_precision=1e-8,
)

The omitted short-range tail is a measured, per-case quantity – choose the zone with a small ladder exactly like the SR pad (sr_image_extent_bohr). Measured fixed-density anchors (pad 20, 2026-08-06): MgO/STO-3G tail 2.7e-3 Ha at zone 10 and 3.2e-4 Ha at zone 12; LiH/STO-3G (most-diffuse exponent 0.048) keeps mHa-scale tail out to ~14 bohr and needs a correspondingly wider zone. Guard rails: requires the corrected exchange split and the padded SR path, must sit strictly below cutoff_bohr, is rejected with SYM3b symmetry reduction / pair-resolved domains, and the analytic-gradient preview fails closed on it (FD forces remain production). The resolved zone is recorded on the result (result.exact_zone_bohr) and in the run log.

Where to look in the code

python/vibeqc/pbc_bipole.py
    run_pbc_bipole_rhf — the closed-shell RHF driver.

python/vibeqc/pbc_bipole_uhf.py
    run_pbc_bipole_uhf — the spin-unrestricted UHF scaffold.

python/vibeqc/pbc_bipole_rks.py
    run_pbc_bipole_rks — the closed-shell RKS (DFT) driver.

python/vibeqc/pbc_bipole_uks.py
    run_pbc_bipole_uks — the open-shell UKS (spin-DFT) driver.

python/vibeqc/bipole_multipole.py
    Shell-pair multipole moments (Phase 1) + interaction tensor
    (Phase 2).

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

python/vibeqc/bipole_cell_moments.py
    Cell-level multipole moments from a density block (Phase 4a).

python/vibeqc/bipole_ext_el_pole.py
    AO-density Fourier helpers, a reciprocal electronic Ewald diagnostic,
    the spheropole, and a fail-closed EXT EL-POLE compatibility entry point.

python/vibeqc/bipole_fock_multipole.py
    Experimental opt-in multipole far-field J builder (Phase 5b) —
    cell-level multipole-multipole interactions for far cell pairs.
    The exact Ewald-J split is the production path.

python/vibeqc/bipole_fock_ewald.py
    Ewald J-split F²e build — short-range erfc + long-range
    reciprocal + background (used by all four drivers).

python/vibeqc/periodic_jk_method.py
    PeriodicJKMethod.BIPOLE enum — selectable via
    run_periodic_job(..., jk_method="bipole").

examples/periodic/input-bipole-*.py
    Example inputs: MgO RHF, MgO RKS PBE, Li UHF, and
    run_periodic_job with jk_method="bipole".

Forces and analytic-gradient preview

Production BIPOLE atomic forces use compute_bipole_gradient_fd. It finite-differences the same total energy reported by the SCF driver, so it is the default path used by geometry optimization and NEB for all four BIPOLE methods.

from vibeqc.bipole_gradient import (
    compute_bipole_gradient_fd,
    compute_bipole_gradient_rhf,
    compute_bipole_gradient_uhf,
    compute_bipole_gradient_rks,
    compute_bipole_gradient_uks,
)

# Production force path:
grad = compute_bipole_gradient_fd(system, "sto-3g", kmesh, opts, method="RHF")
print(f"max|grad| = {np.max(np.abs(grad)):.4e} Ha/bohr")

The analytic drivers are still a research-preview surface. The corrected (Ewald-exchange-split) gauge, the BIPOLE default, is the FD-validated analytic route: RHF/UHF/RKS/UKS at Γ and multi-k, including finite-temperature/fractional (Mermin free-energy) occupations, all pinned against FD. Meta-GGA τ-Pulay is landed (the term validated to 1e-9 against an independent reimplementation; SCAN/r2SCAN are well-behaved, TPSS/M06-L carry a known SCF-eigenvalue residual). In the legacy (use_exchange_ewald_split=False) gauge the maintained preview covers RHF/UHF Γ (general crystals) and maintained RHF/UHF multi-k regressions, Gamma-local zero-smearing RKS/UKS (LDA/GGA XC Pulay, moving-grid, and KS-response terms), and RKS/UKS multi-k (diagonal-Z + corrected W + J^LR + XC Pulay), which warns that the full multi-k KS coupled-perturbed response is not included. Legacy-gauge RKS/UKS calls with finite-temperature/fractional occupations raise NotImplementedError. The sole remaining gated analytic case is the legacy multi-k KS-CPHF: that path uses a diagonal-Z approximation and warns rather than solving the full coupled-perturbed response, and it is deferred (the corrected gauge already covers multi-k KS variationally). Use the FD force path for production forces in all cases.

compute_bipole_gradient_fd costs about 6N SCFs for N atoms. It fails fast if any displaced SCF point does not converge, rather than differentiating a failed iterate; use require_converged=False only for diagnostics.

Structure optimization

Atomic positions and lattice parameters can be optimized using :mod:vibeqc.bipole_optimize:

from vibeqc.bipole_optimize import (
    relax_atoms, relax_cell, relax_cell_gradient, relax_full,
)

# Atomic positions only (FD gradients by default, L-BFGS-B)
result = relax_atoms(system, "sto-3g", kmesh, method="RHF",
                     max_iter=30, conv_tol_grad=1e-4)

# Lattice only — energy FD (Nelder-Mead, no gradient needed)
result = relax_cell(system, "sto-3g", kmesh, method="RHF")

# Lattice only — FD strain gradient (L-BFGS-B)
result = relax_cell_gradient(system, "sto-3g", kmesh, method="RHF")

# Full cell + atom relaxation (alternating)
result = relax_full(system, "sto-3g", kmesh, method="RHF",
                    max_outer=5)

# One-shot via high-level API
from vibeqc.periodic_runner import run_periodic_job
result = run_periodic_job(system, basis, method="RHF", jk_method="bipole",
                          optimize=True, optimize_cell=True)

Atomic relaxation uses fractional coordinates for stability with general lattices. Gradient-based cell optimization finite-differences the SCF energy with respect to the 6 Voigt strain components, so the cell step follows the same objective that is being minimized. The force-virial helper is available only as a diagnostic and is not the periodic stress:

from vibeqc.bipole_gradient import compute_stress_tensor
virial = compute_stress_tensor(system, gradient)  # diagnostic 3×3, Ha/bohr³

Dimensionality support

Dim	Coulomb	Gradient	Optimization	Notes
3D	Ewald J-split (default) or DIRECT_TRUNCATED	FD production; analytic research preview	Atoms + cell	Production-ready energy/FD forces
2D (surfaces)	DIRECT_TRUNCATED only	Analytic RHF	Atoms only	Requires vacuum padding; Ewald J-split raises
1D (wires)	DIRECT_TRUNCATED only	Analytic RHF	Atoms only	Requires vacuum padding; Ewald J-split raises

For 2D/1D, the BIPOLE driver automatically falls back to DIRECT_TRUNCATED when use_ewald_j_split=None (default). This is adequate for molecular-limit cells with sufficient vacuum padding (~10+ bohr). For tight surfaces/wires, the truncated Coulomb sum may diverge, use GDF for those cases.

Orbital and population sidecars

BIPOLE runs (RHF, RKS, UHF, UKS, restricted-open, Γ-only and multi-k alike) write both sidecars. What each one means differs, and the difference matters:

Sidecar	Option	What BIPOLE writes
{out}.population.{txt,json}	write_population_file	A crystal population. Mulliken contracts the real-space lattice density and overlap blocks in the same convention as the BIPOLE SCF energy; Löwdin and Mayer contract over the full SCF k-mesh. No Γ restriction. Dipole stays unsupported: the bulk position operator needs a Berry-phase polarization.
{out}.molden	write_molden_file	The Γ-block orbitals only, truncated to the home cell. Requires the k-mesh to contain Γ.

Prefer QVF and vibe-view over Molden for crystals

Molden is a molecular format with no lattice vectors, no k-points, and no complex coefficients, so a periodic .molden is necessarily a narrow projection: the Γ block, evaluated over home-cell basis functions with no image sum, and therefore clipped wherever an orbital straddles a cell face. Orbitals at k ≠ 0 are irreducibly complex and cannot be written at all — vibe-qc refuses rather than emitting a file a viewer would misread.

Set output_qvf=True and open the archive with vibe-view. QVF carries the lattice, the k-point metadata, the complex Bloch wavefunction, and torus-periodic grids that wrap across cell faces. See QVF and vibe-view and Molden and periodic systems.

All three options are tri-state. None (the default) emits the sidecar when the route can produce it truthfully; True is a guarantee and raises before the SCF starts if the route cannot; False opts out. A shifted Monkhorst-Pack mesh has no Γ point, so write_molden_file=True raises there; pass None to let the run proceed without the orbital file.

References

• Saunders, Dovesi, Roetti 1992, V. R. Saunders, R. Dovesi, C. Roetti, M. Causà, N. M. Harrison, R. Orlando, C. M. Zicovich-Wilson, “CRYSTAL92 User’s Manual”, University of Torino (1992). The original BIELET / BIPOLE description of the multipole-far-pair Coulomb construction.

• CRYSTAL14, R. Dovesi, V. R. Saunders, C. Roetti, R. Orlando, C. M. Zicovich-Wilson, F. Pascale, B. Civalleri, K. Doll, N. M. Harrison, I. J. Bush, P. D’Arco, M. Llunell, M. Causà, Y. Noël, “CRYSTAL14: A program for the ab initio investigation of crystalline solids,” Int. J. Quantum Chem. 114, 1287 (2014). doi:10.1002/qua.24658. The modern reference for the EXT EL-POLE / EXT EL-SPHEROPOLE conventions vibe-qc is targeting.

• CRYSTAL23, A. Erba et al., “CRYSTAL23: a program for computational solid state physics and chemistry,” J. Chem. Theory Comput. 19, 6891 (2023). doi:10.1021/acs.jctc.2c00958.

BIPOLE runs through run_periodic_job(..., jk_method="bipole") now route the methodology references through the citation database (user_guide/citations) via routes.methods["bipole"]. Saunders 1992 and Dovesi 2014 fire automatically on every BIPOLE run.

See also

• Periodic-SCF methods (comparative tour): side-by-side comparison of BIPOLE vs GDF vs GPW/GAPW, plus the surface-reactions workflow walkthrough.

• Multi-k periodic SCF, the production multi-k RHF/RKS path today.

• Density fitting, the production Γ-only hybrid-DFT path via native GDF.

• Ewald summation, the long-range Coulomb machinery used by the legacy run_*_periodic_*_ewald3d drivers.

Quartet-level bipolar far-field (Saunders 1992)

For large cells (cutoff > 8 bohr), a quartet-level bipolar far-field automatically replaces exact four-centre ERIs for distant shell quartets. This is governed by a geometric penetration criterion (Saunders 1992, Sec. 5.3) applied per quartet.

Enabling / disabling

The far-field is auto-enabled by default when both conditions hold:

• lattice_sum_options.cutoff_bohr > 8.0

• The Ewald exchange split is active (the default for 3D systems)

Set use_multipole_far_field explicitly to override:

result = run_periodic_job(system, basis, method="RHF", jk_method="bipole",
    use_multipole_far_field=True,    # force-enable
    multipole_l_max=4,               # hexadecapole (default)
)

# Or disable:
result = run_periodic_job(system, basis, method="RHF", jk_method="bipole",
    use_multipole_far_field=False,   # exact four-centre ERIs only
)

Multipole order

multipole_l_max controls the highest spherical harmonic order (default 4):

L_max	Order name	Spherical components
0	Monopole only	1
1	Dipole	4
2	Quadrupole	9
3	Octupole	16
4	Hexadecapole	25

L=5 and L=6 are available via isotropic adjoined-Gaussian extension (extend_lattice_moments, even-order moments only).

Performance

The far-field infrastructure is built once per geometry and reused across SCF iterations:

• C++ OpenMP contractor: 14.9× faster than Python for 222k quartets

• Pre-computed Fock kernel: interaction tensors cached per unique (R, L, mu)

• Symmetry-reduced dispatch: 2.5× fewer unique quartets for MgO

• C++ moment-derivative buffer: dM/dC precomputed per shell pair for gradient

References

• Saunders, Freyria-Fava, Dovesi, Salasco & Roetti, Mol. Phys. 77, 629 (1992), doi:10.1080/00268979200102671

• Pisani & Dovesi, Int. J. Quantum Chem. 17, 501 (1980), Sec. 4