CRYSTAL-style BIPOLE periodic HF and DFT¶

Feature status

The BIPOLE drivers (run_pbc_bipole_rhf, _uhf, _rks, _uks) form the CRYSTAL-gauge 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 CRYSTAL-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, CRYSTAL-style 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. RHF/UHF Gamma and maintained RHF/UHF multi-k regressions are pinned against FD; Gamma-local zero-smearing RKS/UKS include the current LDA/GGA XC Pulay, moving-grid, and KS-response terms. Broader KS functional/cell coverage and meta-GGA tau-Pulay remain gated. RKS/UKS analytic calls with finite-temperature/fractional occupations or multi-k KS meshes raise NotImplementedError and should use the FD force path.

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.