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 multipole far-field branch is opt-in, experimental, and off by default; the exact Ewald-J split remains the production energy path. 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, but finite differences remain the production recommendation.

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:

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.
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	Analytic dipole-dipole self-image energy
4b proper	landed (2026-05-18)	Ewald reciprocal-sum form of EXT EL-POLE, full multipole content via \(\hat{\rho}(K)\)
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. The remaining Phase 6-7 work 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.

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
    Ewald reciprocal-sum EXT EL-POLE (Phase 4b proper).

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.

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.