Periodic-SCF methods: BIPOLE, GDF, GPW / GAPW

This page is the comparative tour of the periodic Hartree-Fock and Kohn-Sham machinery in vibe-qc. It is meant to be the page you read first when you need to pick a route for a new crystalline (or periodic slab / wire) calculation. The per-method deep dives are in the specialised pages linked from each section.

Periodic boundary conditions are vibe-qc’s headline focus, and surface chemistry (slab single-points, adsorption energies, reaction-path search, transition-state characterisation) is the flagship workflow. Section 7 walks the surface-reactions stack end-to-end; sections 2 to 5 are the per-kernel theory plus implementation reference that workflow sits on top of.

Route status at a glance

Route	Status
GDF (jk_method="gdf", section 3)	Production: Gamma and multi-k RHF / RKS / UHF / UKS in 3-D, and the bounded closed-shell 2-D slab route with analytic gradients. Periodic ROHF on GDF is a preview. gdf_method="mdf" on compact dense-core cells is gated (section 3.4).
GPW (jk_method="gpw", section 4)	Production within its shipped envelope: Gamma RHF / UHF / RKS / UKS, multi-k pure-DFT RKS / ROKS / UKS, Gamma analytic gradients, DFT+U, D3-BJ. The .out route description carries an experimental tag while the GPW / GAPW family is completed (section 4.2 says what is emitted). GPW ROHF / ROKS are a preview.
GAPW (jk_method="gapw", section 4)	Experimental.
BIPOLE (jk_method="bipole", section 2)	Production energies and finite-difference forces (3-D, Gamma and multi-k RHF / UHF / RKS / UKS). Analytic gradients are a preview (a plain UserWarning naming the “maintained preview”); the quartet multipole far field is gated.
Periodic RIJCOSX (jk_method="rijcosx", section 4.8)	Experimental.
2-D slab Ewald (SLAB_EWALD_2D, what AUTO selects on dim=2)	Production for RHF / RKS / UKS; UHF is refused.
Periodic ROHF / ROKS (Ewald, GDF, GPW; section 3.8)	Preview.
AICCM (method="aiccm", every variant)	Experimental.

The label vocabulary is defined in Methods: theory and status.

Note

The cyclic cluster model (CCM): a real-space alternative. Beyond the k-point routes on this page, vibe-qc’s signature periodic approach is the cyclic cluster model, which treats a crystal as a finite cyclic cluster rather than sampling the Brillouin zone (Bredow, Thomas, Geudtner, Gerald, and Jug, Karl (2001). *Development of the cyclic cluster approach for ionic systems.* Journal of Computational Chemistry, **22**(1), 89–101. doi:[10.1002/1096-987X(20010115)22:1<89::AID-JCC9>3.0.CO;2-7](https://doi.org/10.1002/1096-987X(20010115)22:1<89::AID-JCC9>3.0.CO;2-7); at the Hartree-Fock level, Peintinger, Michael F. and Bredow, Thomas (2014). *The cyclic cluster model at Hartree-Fock level.* Journal of Computational Chemistry, **35**(11), 839–846. doi:[10.1002/jcc.23550](https://doi.org/10.1002/jcc.23550)). Two independent experimental ab-initio CCM lines ship in the current release: Γ-CCM / aiccm2026dev-a (union-and-weight/Wigner–Seitz integral weighting) and χ-CCM / aiccm2026dev-b (finite-translation-group characters; 3-D SCF and correlation only). They are distinct approaches; a matched exchange-q=0 convention is necessary for a comparison but does not define an approach delta. Both are selected through one front door, run_periodic_job(method="aiccm", variant=...), documented on the AICCM page. The semi-empirical (MSINDO) CCM is also available for bulk, surface, and adsorption calculations.

See also

• periodic_systems.md: the PeriodicSystem container, Born-von Karman PBC contract, and basic 1D / 2D / 3D setup.

• crystal_lattices.md: the 14 Bravais lattices, with visualisations and worked examples for the common cases (rocksalt, diamond, perovskite, HCP, wurtzite, rutile, corundum, alpha-quartz, graphene).

• k_points.md: Monkhorst-Pack mesh generation and Brillouin-zone sampling.

• bipole.md: full BIPOLE driver reference.

• cyclic_cluster_model.md: the cyclic cluster model, vibe-qc’s real-space (Born-von Karman) alternative to the k-point routes on this page.

• density_fitting.md: RIJ / RIJK / RIJCOSX for molecules, plus the GDF kernel reuse for periodic Gamma-only.

• multi_k_scf.md: the multi-k SCF surface, smearing for metals, and the multi-k GDF roadmap.

• ewald.md: the Ewald summation primitives that the one-electron and the long-range J builds share.

For a direct grid-based KS call to run_krhf_periodic_gdf or run_krks_periodic_gdf, supply PeriodicKSOptions() or omit options to use the periodic XC grid defaults. Supplying PeriodicRHFOptions with a functional raises TypeError before numerical setup: HF options do not carry the required KS grid settings. Ordinary HF calls continue to use PeriodicRHFOptions.

1. Dimensionality model and the periodic Coulomb problem

Every periodic calculation in vibe-qc lives in a single container, PeriodicSystem(dim, lattice, unit_cell):

import numpy as np
import vibeqc as vq
from vibeqc import Atom, PeriodicSystem

# 1D atomic wire (H-H along x). Columns 1 and 2 are vacuum. AUTO and the
# GDF drivers refuse dim=1 (section 3.4); run a wire as dim=3 with this
# transverse vacuum.
sys_1d = PeriodicSystem(
    dim=1,
    lattice=np.diag([4.0, 30.0, 30.0]),
    unit_cell=[Atom(1, [0, 0, 0]),
               Atom(1, [1.4, 0, 0])],
)

# 2D graphene-like slab. Prefer vq.slab_2d: it takes only the two
# in-plane vectors and synthesizes the (non-physical) third column for
# you. The SCF is invariant to that column; it is not a vacuum gap.
sys_2d = vq.slab_2d(
    [4.6487, 0.0000, 0.0],       # a1 (bohr)
    [2.3244, 4.0259, 0.0],       # a2 (bohr)
    [Atom(6, [0.0000, 0.0000, 0.0]),
     Atom(6, [2.3244, 1.3420, 0.0])],
)

# 3D crystalline bulk (silicon, FCC primitive cell, a = 10.26 bohr).
a_si = 10.26
sys_3d = PeriodicSystem(
    dim=3,
    lattice_vectors=[[0, a_si / 2, a_si / 2],
                     [a_si / 2, 0, a_si / 2],
                     [a_si / 2, a_si / 2, 0]],
    unit_cell=[Atom(14, [0, 0, 0]),
               Atom(14, [a_si / 4, a_si / 4, a_si / 4])],   # Cartesian bohr
)

dim selects how many of the three lattice columns are physical. The lattice matrix is always 3x3 and may be triclinic. Reciprocal vectors are generated automatically as \(2\pi A^{-T}\) so that \(a_i \cdot b_j = 2\pi \delta_{ij}\) holds for any cell metric. There are no cubic-only or orthorhombic-only fast paths in the public periodic API: methods accept the full lattice matrix and the test suite pins triclinic, monoclinic, orthorhombic, tetragonal, trigonal, hexagonal, and cubic metrics on every periodic route.

Why a naive lattice sum diverges

The fundamental difficulty is that the periodic two-electron operator,

\[ J_{\mu\nu}[D] = \sum_{\mathbf{T}}\sum_{\lambda\sigma} D_{\lambda\sigma}\, (\mu_{\mathbf{0}}\nu_{\mathbf{0}} | \lambda_{\mathbf{T}}\sigma_{\mathbf{T}}), \]

contains a tail that decays like \(1/|\mathbf{T}|\). In 3D the sum of \(1/|\mathbf{T}|\) over a Bravais lattice is conditionally convergent: it depends on the order of summation, and a brute real-space truncation can be off by tens or hundreds of Hartree on ionic crystals like MgO. The same holds for the electron-nucleus and nucleus-nucleus interactions, where it is just the classical Madelung problem. Periodic codes therefore always do one of three things:

1. Split the Coulomb interaction. Write \(1/r = \mathrm{erfc}(\omega r)/r + \mathrm{erf}(\omega r)/r\). The short-range piece is summed in real space; the long-range smooth piece is summed in reciprocal space. This is the Ewald construction (Ewald, Paul Peter (1921). *Die Berechnung optischer und elektrostatischer Gitterpotentiale.* Annalen der Physik, **369**(3), 253-287. doi:[10.1002/andp.19213690304](https://doi.org/10.1002/andp.19213690304)). BIPOLE (built on the Gaussian-density electrostatics of Saunders, V. R., Freyria-Fava, C., et al. (1992). *On the electrostatic potential in crystalline systems where the charge density is expanded in Gaussian functions.* Molecular Physics, **77**, 629–665. doi:[10.1080/00268979200102671](https://doi.org/10.1080/00268979200102671)) and the GDF Ewald-3D path do this; so does the existing FFT-Poisson driver.

2. Density-fit on a finite auxiliary basis. Project the AO pair density onto a charge-compensated auxiliary basis so that the divergent monopole contribution is cancelled at the source. The resulting fit Coulomb metric is finite. This is the GDF approach (Gaussian density fitting, Sun, Qiming, Berkelbach, Timothy C., et al. (2017). *Gaussian and plane-wave mixed density fitting for periodic systems.* Journal of Chemical Physics, **147**(16), 164119. doi:[10.1063/1.4998644](https://doi.org/10.1063/1.4998644)).

3. Carry the density on a uniform grid and solve the Poisson equation with FFT. Pin \(V(\mathbf{G} = 0) = 0\) for the neutral cell. The smooth-grid Hartree-J is then an \(\mathcal{O}(N_g \log N_g)\) operation, and a Gaussian augmentation around each nucleus restores all-electron accuracy. This is the GPW / GAPW approach (Lippert, Gerald, Hutter, Jürg, and Parrinello, Michele (1997). *A hybrid Gaussian and plane wave density functional scheme.* Molecular Physics, **92**(3), 477-487. doi:[10.1080/00268979709482119](https://doi.org/10.1080/00268979709482119), Lippert, Gerald, Hutter, Jürg, and Parrinello, Michele (1999). *The Gaussian and augmented-plane-wave density functional method for ab initio molecular dynamics simulations.* Theoretical Chemistry Accounts, **103**(2), 124-140. doi:[10.1007/s002140050523](https://doi.org/10.1007/s002140050523)).

vibe-qc ships an implementation of each of these three families. They coexist, and the active method (with all of its parameters) is recorded in the .system manifest so reproducibility never relies on the AUTO heuristic.

Where each method plugs in

The user-facing dispatch lives in PeriodicJKMethod:

import vibeqc as vq
from vibeqc import PeriodicJKMethod

# Three ways to be explicit about the Coulomb route:
result = vq.run_periodic_job(sysp, basis, method="RHF",
                             jk_method="bipole")
result = vq.run_periodic_job(sysp, basis, method="RHF",
                             jk_method="gdf")
result = vq.run_periodic_job(sysp, basis, method="RHF",
                             jk_method=PeriodicJKMethod.GDF)

# The ab initio cyclic cluster model (experimental): method="aiccm" plus a
# variant; the SCF reference is inferred (functional given -> KS, singlet -> RKS).
result = vq.run_periodic_job(sysp, basis, method="aiccm", variant="chi",
                             functional="pbe0",
                             aiccm_lattice_extension=(4, 4, 4),
                             aiccm_backend="rijcosx")

jk_method="auto" (the default) branches on the system’s dimensionality first:

• dim=2 (a slab): always SLAB_EWALD_2D, the vacuum-free 2D Ewald gauge (Parry, D. E. (1975). *The electrostatic potential in the surface region of an ionic crystal.* Surface Science, **49**, 433–440. doi:[10.1016/0039-6028(75)90362-3](https://doi.org/10.1016/0039-6028(75)90362-3), de Leeuw, S. W. and Perram, J. W. (1979). *Electrostatic lattice sums for semi-infinite lattices.* Molecular Physics, **37**(4), 1313–1322. doi:[10.1080/00268977900100951](https://doi.org/10.1080/00268977900100951)), for RHF / RKS / UKS. This AUTO choice remains the general route. Closed-shell RHF/RKS may explicitly select the bounded slab-GDF route; unsupported bulk builders still raise rather than returning an a3-dependent energy.

• dim=1: no route. AUTO raises NotImplementedError for every method, at Gamma and with a k mesh (no consistent 1-D Coulomb gauge; see section 3.4). Model a wire as a dim=3 cell with vacuum on the two transverse axes.

• dim=3: GDF for closed-shell RHF / RKS and BIPOLE for open-shell UHF / UKS.

The resolved method (never the literal "auto") is logged in the banner so a re-run does not depend on the AUTO heuristic version.

A method-by-method status table, kept in sync with python/vibeqc/periodic_jk_method.py:

For the experimental CCM selectors, Γ-CCM and χ-CCM[1] are distinct construction approaches: union-and-weight/Wigner-Seitz integral weighting and the finite-translation-group character construction, respectively. Comparisons hold the declared exchange-q=0 convention fixed as a comparison constraint; that convention does not identify the constructions.

jk_method= (AICCM rows: method="aiccm", variant=)	What it dispatches to	Status
"slab_ewald_2d"	Rigorous vacuum-free 2D (slab) Coulomb, the Parry / de Leeuw-Perram gauge. RHF / RKS / UKS at Gamma and multi-k. This is what AUTO selects for dim=2, so you rarely name it. dim=2 only; raises on bulk or polymer cells. Total energy is invariant to the synthesized a3; E_nuclear is k-mesh independent. No analytic gradients on this direct route (use explicit "gdf" for closed-shell slab forces), no smearing, no UHF yet. Pure and screened-exchange envelopes are supported; full-range multi-k RKS fails closed, while the legacy RHF/UKS full-range exchange paths warn and remain provisional. Aliases: "slab", "ewald_2d".	Production for pure/screened envelopes; full-range RHF/UKS provisional
"gdf"	Native Gaussian density fitting. Closed-shell RHF / RKS / full-range hybrids and maintained-preview ROHF are public at Gamma and multi-k; ROHF uses run_krohf_periodic_gdf and integer 2/1/0 occupations. On dim=2, explicit GDF uses the dedicated signed slab-truncated metric and rigorous Parry-gauge one-electron channels on a full Gamma-centered (n1,n2,1) tuple mesh. Symmetry-reduced (IBZ) meshes from KPoints.monkhorst_pack(..., symmetry=True) are accepted: the driver expands them to the full BZ before the HF exchange sum (correctness, not a speed win: every full-mesh point is still built and diagonalised); a mesh that cannot account for its own full-mesh parent (e.g. an explicit KPoints list without that provenance) still fails closed rather than being expanded against fabricated metadata. Slab analytic gradients ship (2026-07-30): compute_gradient=True / optimize=True on the closed-shell slab route. Slab open shell, ROKS/GDF, smearing, restarts, and DFT+U fail closed. In 3D, open-shell UHF / UKS remains available at Gamma and multi-k; a dim=1 cell raises NotImplementedError for every method, at Gamma and with a k mesh. AUTO refuses every dim=1 job and names the dim=3 vacuum-box model instead, selects slab_ewald_2d for every dim=2 job, and BIPOLE for 3D open shell. Full-grid external XC is experimental and zero-temperature only on the complete periodic AO-density route.	Production in 3D; ROHF maintained preview; verified closed-shell energy + analytic-gradient route in 2D; external XC experimental
method="aiccm", variant="chi" (legacy jk_method="aiccm2026dev-b", "chi", "chi-ccm" resolve with a DeprecationWarning)	χ-CCM, the finite-translation-group character construction evaluated on a Γ-centred character mesh. aiccm_lattice_extension=(N1,N2,N3) defines the real cyclic dimensions and derives the character net; aiccm_wigner_seitz_shells=s requests odd dimensions 2s+1. aiccm_backend= selects "four_center", "ri", or "rijcosx". All backends are 3-D only; every 1-D/2-D request fails closed. Pure full-grid external XC is accepted for zero-temperature RKS/UKS only with four_center; fitted backends and external-provider hybrids remain gated. See the derivation and the AICCM page.	Experimental, opt-in; 3-D only
method="aiccm", variant="real-gamma" (legacy jk_method="real-gamma", "real_gamma" warn)	The neutral Γ-CCM construction in its k-free real-Γ supercell representation, through the per-unit-cell adapter periodic/ccm/real_gamma_runner.py. RHF / RKS / UHF / UKS including screened hybrids; aiccm_lattice_extension (or the legacy Γ-centred kpoints alias, not both) is the BvK torus; fixed-cell optimize=True works for ordinary functionals. Zero-temperature pure full-grid external XC is accepted for RKS/UKS, but its gradients and optimization remain gated. No DFT+U, smearing, Hessian, or cell relaxation; an explicit damping, fock_mixing, fmixing_percent, density_mixer or level_shift fails closed because the direct-torus loop would drop it.	Experimental, opt-in; 3-D only
method="aiccm", variant="neutral-bloch"	The same neutral construction in the Bloch representation, through the producer run_ccm_*_gdf (the unit-cell GDF drivers on the torus mesh) and the adapter periodic/ccm/neutral_bloch_runner.py. RHF / RKS / UHF / UKS; the Gamma-centred mesh is the BvK torus; the producer refuses a vacuum-padded cluster and a positive converged neutral energy, and the arm forwards no mixing control. Full-grid external XC fails closed pending proof that its atom-block periodic partition is representation invariant. No gradients, DFT+U or restart. The former jk_method aliases "neutral-bloch", "bloch-control", "gdf-control", "aiccm-ri" (plain unit-cell GDF) are retired and fail closed.	Experimental, opt-in; 3-D only
method="aiccm", variant="four-center" (legacy jk_method="aiccm2026dev-a" warns)	The union-and-weight / Wigner-Seitz four-centre Γ-CCM lineage through the per-unit-cell adapter. RHF, RKS, UHF, and UKS run in 3-D; HF uses the scalable lattice-sum builder, while the front door always selects the symmetric aiccm2026dev-a weighting. Zero-temperature pure full-grid external providers are accepted for RKS/UKS. External-provider hybrids and ECP-bearing external-XC jobs remain gated, as do DFT+U, smearing, gradients and optimization, Hessians, and explicit mixing controls. The bare "gamma", "gamma-ccm", "gamma_ccm" spellings fail closed with the ruling-R1 message on both surfaces.	Experimental, opt-in; 3-D only
"bipole"	an exact Ewald J-split; multi-k RHF / UHF / RKS / UKS, plus maintained-preview ROHF and ROKS (Gamma + full Monkhorst-Pack meshes, energies only; ROHF and nondegenerate ROKS use integer 2/1/0 occupations, while an exactly degenerate one-electron ROKS shell uses a fixed equal-occupation ensemble; ROKS covers pure functionals and global hybrids) on the corrected-Ewald-exchange EWALD_3D engine; production FD forces plus fixed-cell atomic optimisation; variable-cell work fails closed and analytic gradients remain a gated preview. The public route is 3D-only and rejects dim=1 and dim=2; model a wire as a dim=3 cell with transverse vacuum, use "auto" / "slab_ewald_2d" for surfaces, or a physical dim=3 model when BIPOLE forces are required. Full-grid external XC is experimental, zero-temperature only, and requires the corrected gauge plus a complete Monkhorst-Pack mesh or valid expandable IBZ.	Production (3D exact route); external XC experimental
"direct"	4-centre real-space lattice sum with Schwarz screening. Valid only for vacuum-padded molecular-limit cells.	Diagnostic
"fft_poisson"	Retired (v0.13.0). The Γ-only EWALD_3D path returned a ~2 Ha-wrong energy on cells where periodic AO images overlap (its molecular-limit density convention D(g≠0)=0 breaks down, e.g. LiH rocksalt / STO-3G, whose Li 2s/2p are very diffuse). Requesting it now raises; use "gdf" (default) / "bipole" / "gpw". The internal Γ-only EWALD drivers remain for dilute periodic systems and fail closed on the image-overlap regime.	Retired
"gpw"	Gaussian-on-grid collocation + FFT-Poisson Hartree-J on a smooth real-space grid. The dispatcher runs a full SCF for Gamma RHF / maintained-preview ROHF / RKS / UHF / UKS and multi-k RKS, with forces, FD Hessian, +U for RHF/RKS/UHF/UKS at Gamma and RKS at multi-k, smearing, and D3(BJ). Multi-k open-shell +U is not wired. ROHF is 3D Gamma-only with integer occupations and fails closed for gradients, smearing, and the other unsupported envelopes. The standalone run_periodic_rohf_gpw entry point remains available. Full-grid external XC has a narrower contract: zero-temperature 3D RKS at Gamma or on a complete unreduced multi-k mesh. It bypasses pointwise FFT XC but retains FFT-Poisson Hartree-J; positive smearing, ECP-bearing jobs, and symmetry-reduced representative meshes fail closed. The .out route line carries an experimental tag and the grid infrastructure still emits GAPWExperimentalWarning.	Production within its shipped envelope; ROHF maintained preview; external XC experimental
"gapw"	GPW plus per-atom Gaussian augmentation on a fine radial grid for all-electron accuracy. The standard drivers cover Gamma RHF / RKS / UHF / UKS and multi-k RKS, with +U for RHF/RKS/UHF/UKS at Gamma and RKS at multi-k, smearing, and D3(BJ). Multi-k open-shell +U is not wired. At single Gamma, the method-aware one_centre="auto" policy selects the fit-free analytic-ERI augmentation for RHF/UHF only when the caller declares molecular_limit=True (gapw_molecular_limit=True on run_periodic_job); otherwise HF fails closed instead of guessing from geometry or falling back to the known-wrong block functional. RKS/UKS retain block. Analytic HF is validated only for vacuum-padded molecular-limit cells, not compact crystals or multi-k GAPW; use GDF or BIPOLE there. The high-level GAPW HF route rejects geometry/cell optimization and Hessians; VibeqcGAPW supplies fixed-cell forces, using full-SCF finite differences in analytic mode. HF phonons require explicit block mode, while DFT phonons default to block. Full-grid external XC is narrower still: zero-temperature 3D Gamma-only RKS, with the split smooth/hard/soft XC construction bypassed in favor of one complete AO-density evaluation. Positive smearing, multi-k, open-shell, and ECP-bearing external-XC jobs fail closed. The .out route line carries an experimental tag.	Experimental (analytic HF molecular-limit only); external XC experimental
"rijcosx"	RI-J plus seminumerical chain-of-spheres exchange K (Neese 2009, Neese, Frank, Wennmohs, Frank, et al. (2009). *Efficient, approximate and parallel Hartree-Fock and hybrid DFT calculations. A ‘chain-of-spheres’ algorithm for the Hartree-Fock exchange.* Chemical Physics, **356**(1–3), 98–109. doi:[10.1016/j.chemphys.2008.10.036](https://doi.org/10.1016/j.chemphys.2008.10.036)) on a periodic Becke grid. Gamma RHF uses the dedicated periodic RIJCOSX driver. True multi-k RHF / RKS / UHF / UKS route through the native GDF loop with k_exchange="cosx"; RKS/UKS hybrids exercise COSX exchange, while pure functionals have zero exact-exchange weight and therefore validate the RI-J plumbing. One-cell RKS/UHF/UKS requests fail closed rather than being mislabeled as RIJCOSX. Full-grid external XC is experimental and zero-temperature only on an unreduced complete multi-k mesh.	Implemented, experimental periodic COSX
"rsgdf"	Top-level route not implemented (enum-reserved): range-separated GDF as a standalone jk_method (Ye and Berkelbach, 2021, Ye, Hong-Zhou and Berkelbach, Timothy C. (2021). *Fast periodic Gaussian density fitting by range separation.* Journal of Chemical Physics, **154**(13), 131104. doi:[10.1063/5.0046617](https://doi.org/10.1063/5.0046617)). Distinct from the internal gdf_method="rsgdf" backend, which is the production default the GDF route uses for tight ionic cells (see § 3.2.1).	Roadmap
"cfmm"	Continuous fast multipole (White, Christopher A., Johnson, Benny G., et al. (1994). *The continuous fast multipole method.* Chemical Physics Letters, **230**(1–2), 8–16. doi:[10.1016/0009-2614(94)01128-1](https://doi.org/10.1016/0009-2614(94)01128-1)). Enum-reserved; not implemented.	Roadmap

Two periodic method families also expose dedicated entry points in addition to the public runner where noted:

• RSGAPW (range-separated GAPW) for range-separated hybrids (HSE06, ωB97X, CAM-B3LYP): run_periodic_rhf_rsgapw and run_periodic_rks_rsgapw in periodic_gapw_range_sep.py. It splits the Coulomb operator with the error function and carries the long-range piece on the plane-wave grid. See § 4.7.

• Periodic ROHF / ROKS (restricted open-shell HF and KS, Roothaan, Clemens C. J. (1960). *Self-Consistent Field Theory for Open Shells of Electronic Systems.* Reviews of Modern Physics, **32**(2), 179–185. doi:[10.1103/RevModPhys.32.179](https://doi.org/10.1103/RevModPhys.32.179), Guest, M. F. and Saunders, V. R. (1974). *On methods for converging open-shell Hartree-Fock wave-functions.* Molecular Physics, **28**(3), 819–828. doi:[10.1080/00268977400102171](https://doi.org/10.1080/00268977400102171)): run_rohf_periodic_gamma_ewald3d (Gamma) and run_rohf_periodic_multi_k_ewald3d (multi-k) in periodic_rohf_ewald.py / periodic_rohf_multi_k_ewald.py, run_roks_periodic_multi_k_ewald3d (Gamma + multi-k) in periodic_roks_multi_k_ewald.py, run_krohf_periodic_gdf (Gamma + multi-k, native GDF) in periodic_rohf_gdf.py, plus the standalone 3D Gamma GPW routes run_periodic_rohf_gpw / run_periodic_roks_gpw in periodic_gapw_open_shell.py. See § 3.8.

The method-aware one-centre promotion applies only to the standard Gamma run_periodic_rhf_gapw / run_periodic_uhf_gapw drivers and their public run_periodic_job(jk_method="gapw") dispatch. The dedicated RSGAPW and OT research entry points use separate experimental builders and do not inherit that default or its molecular-limit validation claim.

Screened hybrids (HSE06, Heyd, Jochen, Scuseria, Gustavo E., and Ernzerhof, Matthias (2003). *Hybrid functionals based on a screened Coulomb potential.* Journal of Chemical Physics, **118**(18), 8207–8215. doi:[10.1063/1.1564060](https://doi.org/10.1063/1.1564060), Krukau, Aliaksandr V., Vydrov, Oleg A., et al. (2006). *Influence of the exchange screening parameter on the performance of screened hybrid functionals.* Journal of Chemical Physics, **125**(22), 224106. doi:[10.1063/1.2404663](https://doi.org/10.1063/1.2404663)) per route. Every periodic driver resolves its functional through the shared CAM assembly (vibeqc.periodic_screened_exchange): K_HF = c_full·K_full + c_sr·K_erfc(ω_screen). Global hybrids (PBE0, Adamo, Carlo and Barone, Vincenzo (1999). *Toward reliable density functional methods without adjustable parameters: The PBE0 model.* Journal of Chemical Physics, **110**(13), 6158–6170. doi:[10.1063/1.478522](https://doi.org/10.1063/1.478522)) keep the full-range fraction unchanged. Screened hybrids (HSE06, c_full = 0) build the genuine erfc short-range exchange (a pure direct lattice sum with no Madelung/exxdiv seam and no bipolar far-field, the CRYSTAL treatment) on the routes that have that kernel: BIPOLE, the Ewald drivers (Γ + multi-k, RKS + UKS), and the slab route (slab_ewald_2d / AUTO on dim=2). Routes whose exact exchange is full-range only (GDF Γ + multi-k, GPW/GAPW, COSX) fail closed on any range-separated functional instead of silently running HSE06 as PBE0 (the pre-2026-07-12 behaviour, on every route). Long-range-corrected RSH (ωB97X, CAM-B3LYP, LC-ωPBE) fail closed everywhere periodically: their K_full arm needs an exxdiv treatment that is not validated. The experimental RSGAPW entry points above are unchanged by this and remain the GAPW-side option.

The full-range (c_full) arm on a multi-k mesh uses the corrected Ewald split. On a finite k mesh the SCF density is Born-von-Kármán-torus periodic, so D(g) does not decay with image distance and a bare 1/r image sum does not converge: it drifts without limit as the lattice cutoff grows. The multi-k EWALD_3D drivers (RHF, UHF, RKS, UKS) therefore build that arm as erfc(α r)/r in real space plus a reciprocal channel per momentum transfer q = k − k' plus the probe-charge q + G = 0 term: the exxdiv='ewald' convention (Sundararaman, R. and Arias, T. A. (2013). *Regularization of the Coulomb singularity in exact exchange by Wigner-Seitz truncated interactions: Towards chemical accuracy in nontrivial systems.* Physical Review B, **87**, 165122. doi:[10.1103/PhysRevB.87.165122](https://doi.org/10.1103/PhysRevB.87.165122); the integrable q = 0 singularity treatment goes back to Gygi, F. and Baldereschi, A. (1986). *Self-consistent Hartree-Fock and screened-exchange calculations in solids: Application to silicon.* Physical Review B, **34**, 4405–4408. doi:[10.1103/PhysRevB.34.4405](https://doi.org/10.1103/PhysRevB.34.4405)), shared with BIPOLE and implemented once in vibeqc.periodic_corrected_exchange. Consequences worth knowing:

• A hybrid or HF multi-k job no longer matches the Γ-only drivers bit-for-bit. Those keep their own truncated gauge, which is legitimate at Γ where the molecular-limit density does decay. Pure functionals build no full-range arm and are unaffected, as are screened hybrids (c_full = 0).

• One gauge at every mesh, (1,1,1) included (GitLab #23). The multi-k RHF and UHF drivers take the corrected split by default on a single Γ point too, so a k-convergence sequence run with the defaults stays in one gauge and matches PySCF exxdiv='ewald' at every step. Until #23 a single Γ point took the molecular-limit kernel instead, and stepping (1,1,1) to (2,2,2) on H₂/STO-3G in an 8-bohr cube mixed 4.4 mHa of real k convergence with a 21.7 mHa gauge switch. That kernel is still available as exchange_exxdiv='none', accepted only at a single Γ point (it zeroes every image block, so the sum has one term) and refused on any real mesh, where that sum does not converge. The two are different finite-size conventions with the same infinite-k limit (Carrier, Rohra & Görling, Phys. Rev. B 75, 205126 (2007), Sec. IV, Carrier, Pierre, Rohra, Stefan, and Görling, Andreas (2007). *General treatment of the singularities in Hartree-Fock and exact-exchange Kohn-Sham methods for solids.* Physical Review B, **75**, 205126. doi:[10.1103/PhysRevB.75.205126](https://doi.org/10.1103/PhysRevB.75.205126)): the molecular limit converges exponentially to the isolated molecule as a box grows, the corrected split as 1/L³ (0.37 mHa at L = 30 bohr). A Γ-supercell energy band-folds onto a multi-k energy only when both sides use the same gauge. See tests/test_periodic_exchange_gauge_band_folding.py.

• The Γ-only drivers are the molecule-in-a-box routes and keep the molecular limit: run_rhf_periodic_gamma_ewald3d and its UHF, ROHF, RKS and UKS siblings, run_rhf_periodic_gamma_scf, run_rks_periodic_gamma_scf, the legacy Γ GDF driver, and the ASE calculator at its default kpts=(1, 1, 1). An option that makes such a route borrow the multi-k driver keeps the gauge: smearing on run_rhf_periodic_gamma_scf passes exchange_exxdiv='none', and run_rks_periodic_gamma_scf refuses smearing with a hybrid (the multi-k KS driver has no molecular-limit kernel) and names run_rks_periodic_scf. The k-mesh dispatchers run_rhf_periodic_scf (which forwards exchange_exxdiv) and run_rks_periodic_scf use the multi-k drivers at every mesh; until #23 run_rks_periodic_scf sent a Γ mesh to the Γ-only driver unless smearing or a restart density was set, which moved a PBE0 energy by 18.9 mHa on H₂/STO-3G in a 6-bohr cube. Pure functionals give the same energy on both (to 1e-14 Ha there).

• Every result records the gauge it ran as result.exchange_q0 (BvK-ewald, molecular-limit, or bare-image-sum on the gauges the corrected split does not exist for yet, 2D slab and dim != 3). The analytic gradient (compute_gradient_periodic_rhf_multi_k) differentiates that gauge rather than inferring one from the mesh, and run_periodic_job writes it to .system [run].exchange_q0, the QVF provenance, and an exchange q=0 gauge: line in the .out.

• A convergence aid no longer changes the Γ GDF Hamiltonian (GitLab #400). A default-Γ bulk RHF/RKS job that the exact-Γ driver cannot evaluate (a level shift, Fock mixing for RHF, smearing, or a charged cell) runs the k-mesh engine at kpoints=(1, 1, 1), exactly what an explicit one-point mesh selects, and records BvK-ewald. Before, it fell back to the legacy Ewald-J driver, whose exchange carries no divergence correction: 75.65 mHa on H₂/STO-3G in a 6-bohr cube for any of those aids, 18.4 mHa on a charged HeH⁺ cell. symmetry_stabilize and symmetry_reduce_fock=True are refused on bulk GDF, since no GDF engine implements them (a k-mesh job used to drop them silently).

• Which engine ran, and why, is in the record. run_periodic_job decides the GDF engine once (vibeqc.periodic_gdf_route.resolve_gdf_route) and writes [run].gdf_engine (exact-gamma, k-mesh, legacy-gamma, slab-gdf), gdf_driver and gdf_route_reason to the .system manifest next to exchange_q0; the .out carries a gdf_engine line and a GDF engine: progress line. The rule for convergence aids: an explicitly requested aid the exact-Γ driver lacks moves the job to the k-mesh engine at (1, 1, 1), which implements it; an AUTO-resolved Fock mixing or level shift never moves the engine and is dropped (the strategy block says so) only when the chosen driver lacks it. Smearing always moves the engine, whatever its source. Calling an open-shell k-mesh driver directly with a level shift or Fock mixing raises instead of dropping the option.

• Until GitLab #399 the multi-k UHF driver kept the bare image sum on every mesh, so a closed-shell singlet on a compact cell did not reproduce RHF: H₂/STO-3G in a 6-bohr cube at (2,1,1) read -4.43 Ha against RHF’s -1.19 Ha, and kept falling as the lattice cutoff grew. It now takes the same gauge selection as RHF, and UHF(M=1) equals RHF on the same mesh and options.

• A symmetry-reduced (IBZ) Monkhorst-Pack mesh works, wedge-native: the SCF diagonalises only at the irreducible points and unfolds the density to the full Brillouin zone per iteration with the Pisani/Dovesi star transport (D(Vk) = W D(k) W^H plus time reversal, Dovesi, Roberto (1986). *On the role of symmetry in the ab initio Hartree-Fock linear-combination-of-atomic-orbitals treatment of periodic systems.* International Journal of Quantum Chemistry, **29**, 1755–1774. doi:[10.1002/qua.560290608](https://doi.org/10.1002/qua.560290608)). This covers 3D full-range and screened hybrids and 2D slab screened hybrids, including oblique lattices. The exchange quadrature is exact. K-dependent diagonalisation and Bloch work shrink with the wedge; 3D exchange can also use orbit-reduced real-space builds. Slabs retain the raw finite-cutoff exchange operator used by their full mesh, so their dominant real-space exchange build does not shrink by the star multiplicity; only the density is reconstructed. Caveat: at loose lattice cutoffs the truncated operators themselves break the point symmetry, so wedge and full-mesh answers can still differ by the truncation asymmetry; both coincide once the sums are fold-converged (the default auto_optimize_truncation regime). An explicit k list carries no star metadata: a full-range arm fails closed, while a screened-only arm treats its supplied weights as the user’s quadrature and makes no symmetry-unfolding claim.

• 2D slab and dim != 3 cells cannot use the corrected full-range split (the reciprocal q-channel cache is 3D-only). Full-range behavior is route-specific: RKS fails closed, legacy RHF/UKS paths warn and remain provisional, and long-range-corrected functionals are rejected by the periodic exchange resolver. HSE-type screened hybrids have c_full = 0, need no reciprocal q channel, and are supported by the direct erfc exchange route, including the 2D IBZ unfolding above.

The rest of this page is the side-by-side treatment of the three families that carry most day-to-day work: BIPOLE, GDF, and the GPW / GAPW plane-wave route, with RIJCOSX and RSGAPW noted alongside.

Validating a periodic route? Each JK route is its own implementation (never an external QC program at runtime). GDF and GPW/GAPW have same-family out-of-process references in PySCF and CP2K/GPAW. CRYSTAL is a component and convergence reference for BIPOLE, but the supported vibe-qc exact Ewald-J route does not execute CRYSTAL’s quartet replacement. The periodic JK routes & parity policy page is the map of which reference belongs to which route (and why cross-family comparison is a trap). It is the periodic backbone reference for the method paper.

2. BIPOLE: an Ewald J-split

2.1 Theory

BIPOLE is vibe-qc’s own direct-space Gaussian periodic SCF route. CRYSTAL famously splits the Fock build into two halves that live in the same electrostatic gauge:

• The one-electron Coulomb \(V_{ne}\) and the nuclear-nuclear energy \(E_{nn}\) are evaluated with 3D Ewald. The same Ewald parameter \(\omega\) is shared by both (a single shared Ewald state). Direct truncation diverges on charged-nucleus crystals because the point-charge / Gaussian-pair interaction on the nucleus side does not decay.

• The two-electron Fock matrix \(F^{2e} = J - \tfrac{1}{2}K\) is assembled in CRYSTAL by direct BIPOLE screening and multipole far-pair replacement (Dovesi, R., Pisani, C., et al. (1983). *Treatment of Coulomb interactions in Hartree-Fock calculations of periodic systems.* Physical Review B, **28**, 5781–5792. doi:[10.1103/PhysRevB.28.5781](https://doi.org/10.1103/PhysRevB.28.5781)). The supported vibe-qc route instead uses an exact Ewald-3D J-split in the same electrostatic gauge,

  \[ J(\omega) = J^{SR}(\omega) + J^{LR}(\omega) + V_{bg}\cdot S, \qquad V_{bg} = -\frac{\pi N_e}{\omega^2 V}, \]

  where \(\omega\) is the same Ewald parameter used by \(V_{ne}\) and \(E_{nn}\). \(J^{SR}\) uses erfc-screened real-space ERIs; \(J^{LR}\) is built from a small reciprocal-space sum over AO-pair Fourier transforms; the uniform background restores charge neutrality. Exchange \(K\) is assembled fully in direct space (Causà, M., Dovesi, R., et al. (1988). *Treatment of the exchange interactions in Hartree-Fock LCAO calculations of periodic systems.* The Journal of Physical Chemistry, **92**, 909–913. doi:[10.1021/j100315a010](https://doi.org/10.1021/j100315a010)) via shell-pair Schwarz screening (Häser, Marco and Ahlrichs, Reinhart (1989). *Improvements on the direct SCF method.* Journal of Computational Chemistry, **10**(1), 104–111. doi:[10.1002/jcc.540100111](https://doi.org/10.1002/jcc.540100111); the same machinery that the molecular four-index path uses).

The periodic quartet expansion is given in Chapter II.4c of Pisani, Dovesi, and Roetti (1988) (Pisani, Cesare, Dovesi, Roberto, and Roetti, Carla (1988). *Hartree-Fock Ab Initio Treatment of Crystalline Systems.* Lecture Notes in Chemistry, Vol. 48, Springer, **48**. doi:[10.1007/978-3-642-93385-1](https://doi.org/10.1007/978-3-642-93385-1)), where Eq. II.4.10 retains the common lattice-translation sum. Saunders et al. (1992) (Saunders, V. R., Freyria-Fava, C., et al. (1992). *On the electrostatic potential in crystalline systems where the charge density is expanded in Gaussian functions.* Molecular Physics, **77**, 629–665. doi:[10.1080/00268979200102671](https://doi.org/10.1080/00268979200102671)) develops the related Gaussian-density electrostatic potential, shell-penetration, and spheropole treatment; it does not derive vibe-qc’s quartet dispatcher. Dovesi et al. (2014) (Dovesi, R., Orlando, R., et al. (2014). *CRYSTAL14: A program for the ab initio investigation of crystalline solids.* International Journal of Quantum Chemistry, **114**, 1287–1317. doi:[10.1002/qua.24658](https://doi.org/10.1002/qua.24658)) and Erba et al. (2023) provide modern CRYSTAL context for the EXT EL-POLE / EXT EL-SPHEROPOLE conventions. The BIPOLE name in vibe-qc explicitly nods to that lineage.

2.2 What ships today

The BIPOLE driver is split across four method-flavour entry points in pbc_bipole.py and friends:

from vibeqc.pbc_bipole     import run_pbc_bipole_rhf
from vibeqc.pbc_bipole_uhf import run_pbc_bipole_uhf
from vibeqc.pbc_bipole_rks import run_pbc_bipole_rks
from vibeqc.pbc_bipole_uks import run_pbc_bipole_uks

All four support:

• 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)\) (Löwdin, Per-Olov (1970). *On the Nonorthogonality Problem.* Advances in Quantum Chemistry, **5**, 185–199. doi:[10.1016/S0065-3276(08)60339-1](https://doi.org/10.1016/S0065-3276(08)60339-1)). The historical enumeration keeps the translations with \(|g| \le\) cutoff_bohr; the physically right criterion is the pair separation \(|O_\mu - O_\nu - g|\) (Sharma and Beylkin, JCTC 17, 3916 (2021), Eqs. 17-18, Sharma, Sandeep and Beylkin, Gregory (2021). *Efficient Evaluation of Two-Center Gaussian Integrals in Periodic Systems.* Journal of Chemical Theory and Computation, **17**(7), 3916-3926. doi:[10.1021/acs.jctc.0c01195](https://doi.org/10.1021/acs.jctc.0c01195)), and the two differ on any cell with an intra-cell offset: LiH rocksalt / def2-SVP moves by 0.72 in an overlap element at the shipped cutoff when the same crystal is described with an atom one lattice vector away. Since GitLab #429 stage 2, opts.lattice_opts.pair_complete_1e = True enumerates the pair-complete term set for \(S\), \(T\), the direct \(V_{ne}\) and their gradient partners (invariant to round-off; the cell list grows, MgO at the default from 135 to 369 cells, and the padded cells that no pair reaches stay exactly zero). The option remains off by default pending cutoff convergence and performance review. vibeqc.pair_complete_lattice_cells returns the one-electron enclosure. The two-electron enclosure also includes the interaction reach between physical product midpoints: J restricts the AB/CD products, while K restricts AC/BD. Exchange output and density cells can extend beyond the one-electron support. Consumers align these lists by integer cell key and extend one-electron matrices with zeros where appropriate. eri_interaction_cutoff_bohr sets the midpoint reach; zero uses cutoff_bohr, and the BIPOLE SR extent overrides that reach when supplied. The Ewald-split nuclear family (analytical erfc, grid and FT \(V_{ne}\), SAP and their gradient partners) uses the same AO-pair filter in each complementary term. Nuclear source images are centered on each physical AO-product midpoint. The reciprocal cache and derivatives preserve the same pair support. Physical mode uses the full Fock build. The legacy symmetry projector, periodic COSX and diagnostic grid Hartree backend refuse this option because they do not implement its exchange/domain contract. Image relabelling preserves a fixed physical domain; convergence of the periodic operator additionally requires converging the pair cutoff, interaction reach and reciprocal resolution.

• The shared-gauge Ewald J-split two-electron build by default on 3D systems (use_ewald_j_split=None -> True): short-range \(J\) from direct erfc-screened ERIs, long-range \(J\) from reciprocal AO-pair Fourier transforms, matching electron background potential, full direct-space exchange.

• SCF inner loop with Bloch-sum -> direct real-space energy evaluation (\(\sum_g D(g) M(g)\)) -> optional multi-k Pulay DIIS (Pulay, Péter (1980). *Convergence acceleration of iterative sequences. The case of SCF iteration.* Chemical Physics Letters, **73**(2), 393–398. doi:[10.1016/0009-2614(80)80396-4](https://doi.org/10.1016/0009-2614(80)80396-4)) -> optional level-shift -> diagonalisation -> optional MOM (Gilbert, Andrew T. B., Besley, Nicholas A., and Gill, Peter M. W. (2008). *Self-Consistent Field Calculations of Excited States Using the Maximum Overlap Method (MOM).* Journal of Physical Chemistry A, **112**(50), 13164–13171. doi:[10.1021/jp801738f](https://doi.org/10.1021/jp801738f)) occupied-subspace reorder -> density rebuild.

• Multi-k acceleration in the dense-k regime. Symmetry-reduced Monkhorst-Pack (Monkhorst, Hendrik J. and Pack, James D. (1976). *Special points for Brillouin-zone integrations.* Physical Review B, **13**(12), 5188–5192. doi:[10.1103/PhysRevB.13.5188](https://doi.org/10.1103/PhysRevB.13.5188)) inputs are accepted and internally expanded back to the full mesh until true IBZ orbit expansion lands.

• Parameter-free bz_integration="gilat" occupations (Gilat, G. and Raubenheimer, L. J. (1966). *Accurate Numerical Method for Calculating Frequency-Distribution Functions in Solids.* Physical Review, **144**, 390-395. doi:[10.1103/PhysRev.144.390](https://doi.org/10.1103/PhysRev.144.390)) on the public RKS, UHF, and UKS routes. UHF/UKS integrate the fixed alpha and beta populations independently with single-spin degeneracy and zero entropy. This is a supported self-consistent path for gapped systems; metallic calculations use smearing for SCF and Gilat post-SCF because the sharp T = 0 occupation map is discontinuous at band crossings. RHF/ROHF/ROKS remain gated here.

• Periodic DFT quadrature uses an atom-image neighborhood centered on each integration point. becke_image_radius_bohr sets its minimum physical radius; it expands to twice the nearest-image distance where needed to keep vacuum regions covered. It does not truncate the atomic radial grid. Converge this radius together with radial and angular grid resolution, especially for diffuse densities and large vacuum regions. A radius of zero explicitly requests molecular partitioning. Moving an atom to an equivalent image preserves the physical partition; raw atomic weights and atom ownership retain the same layout.

• Production finite-difference atomic forces for all four flavours and fixed-cell atomic relaxation. Variable-cell and coupled atom/cell optimization fail closed; compute_stress_tensor is only a force-virial diagnostic. See bipole.md for the gradient and optimisation surface.

Dense-k STO-3G sign-off vs CRYSTAL14 at cutoff = 14 bohr:

Demo	vibe-qc \(E_{\text{total}}\) (Ha / cell)	\(\Delta\) vs CRYSTAL14
MgO SHRINK 8	\(-271.2177748509\)	\(+0.369\;\text{mHa}\)
Diamond SHRINK 8	\(-74.8771393842\)	\(-0.145\;\text{mHa}\)
Silicon SHRINK 8	\(-571.3214715798\)	\(-0.659\;\text{mHa}\)

These recorded sign-off numbers predate the repository split. The old tests/demos/ path is no longer available in this tree; current diagnostics are in examples/regression/bipole_parity/. Do not treat the historical figures as a fresh validation of the current build.

2.3 Advantages

• Multi-k production today, including multi-k hybrids. Multi-k SCF on UHF / UKS as well as RHF / RKS, and the route that carries exact-exchange K at \(\mathbf{k}\neq 0\) for hybrid functionals. (Open-shell GDF also runs multi-k UHF / UKS via run_kuhf_periodic_gdf / run_kuks_periodic_gdf; In 3D, BIPOLE remains the AUTO open-shell default and the multi-k hybrid route; 1D AUTO refuses, naming the dim=3 vacuum-box model, and 2D AUTO uses the slab route.)

• Global and screened hybrid DFT work on the supported envelope. Global hybrids and the validated HSE screened-exchange branch (Heyd, Jochen, Scuseria, Gustavo E., and Ernzerhof, Matthias (2003). *Hybrid functionals based on a screened Coulomb potential.* Journal of Chemical Physics, **118**(18), 8207–8215. doi:[10.1063/1.1564060](https://doi.org/10.1063/1.1564060)) share BIPOLE’s direct-K machinery. Long-range-corrected range-separated hybrids fail closed, as do VV10/nonlocal-correlation and double-hybrid functionals whose additional energy terms are absent from the BIPOLE driver.

• Production FD forces and fixed-cell optimisation. Atom optimisation runs in fractional coordinates. Variable-cell optimisation fails closed pending a certified stress and coupled convergence implementation.

• External convergence reference. Matched basis and k-mesh calculations provide useful component and convergence checks, but the supported exact Ewald-J route is not the same truncation algorithm as a bipolar-expansion calculation and no universal per-mHa cross-code guarantee is claimed.

• Diagnostic decomposition. Each iteration logs \(E_{kinetic}\), \(E_{nuclear\_attract}\), \(E_{J^{SR}}\), \(E_{J^{LR}}\), \(E_{exchange}\), \(E_{nn}\) separately so a parity mismatch can be localised to a single Coulomb contribution.

2.4 Disadvantages

• 3D only through the public BIPOLE route. AUTO refuses a dim=1 wire, and a dim=2 slab routes to SLAB_EWALD_2D. An explicit public jk_method="bipole" request rejects both dimensions before setup. The low-dimensional direct-only BIPOLE branch is a direct-API diagnostic.

• Cost scales with the real-space cutoff. Direct \(K\) at large cutoff is the dominant per-iteration cost. Corrected-gauge direct RHF/RKS/UHF/UKS runs therefore stop before the first Fock build when the cheap overlap-fold preflight is in the critically truncated drift > 1e-2 regime; increase bipole_cutoff_bohr until the drift is below 1e-4 for quantitative work. If that converged cutoff is too expensive, use a supported GDF or GPW route. The experimental multipole-far-pair branch is a retired research artifact, not a supported shortcut for this cost.

• No true IBZ acceleration yet. Symmetry-reduced Monkhorst- Pack inputs are accepted, but they are internally expanded to the full mesh before the Ewald-J build, so the cost is the same as a full-mesh run.

2.5 When to use BIPOLE

• You need an open-shell (UHF / UKS) periodic SCF.

• You need production FD forces or fixed-cell atom optimisation.

• You are doing CRYSTAL parity work.

• You need multi-k DFT today, including hybrid DFT, on dense meshes (2x2x2 through ~8x8x8) and modest basis sizes (STO-3G through pob-TZVP / def2-SVP class).

2.6 Quick start

import numpy as np
import vibeqc as vq
from vibeqc.pbc_bipole import run_pbc_bipole_rhf

ANG2BOHR = 1.0 / 0.529177210903
a = 4.21 * ANG2BOHR
lattice = (a / 2.0) * np.array(
    [[0.0, 1.0, 1.0],
     [1.0, 0.0, 1.0],
     [1.0, 1.0, 0.0]]
)
sysp = vq.PeriodicSystem(
    3, lattice,
    [vq.Atom(12, [0.0, 0.0, 0.0]),
     vq.Atom(8,  [a / 2.0, a / 2.0, a / 2.0])],
)
basis = vq.BasisSet(sysp.unit_cell_molecule(), "sto-3g")

kmesh = vq.monkhorst_pack(sysp, [2, 2, 2], use_symmetry=False)

opts = vq.PeriodicRHFOptions()
opts.lattice_opts.cutoff_bohr = 14.0
opts.lattice_opts.nuclear_cutoff_bohr = 14.0
opts.initial_guess = vq.InitialGuess.SAD
opts.use_diis = True
opts.max_iter = 30

result = run_pbc_bipole_rhf(
    sysp, basis, kmesh, opts,
    use_ewald_j_split=True,
    ewald_precision=1e-8,
)
print(f"E_total = {result.energy:+.6f} Ha")

End-to-end MgO RKS / PBE input lives at examples/periodic/input-bipole-mgo-rks-pbe.py, open-shell Li UHF at examples/periodic/input-bipole-li-uhf.py, and the unified-runner workflow at examples/periodic/input-bipole-workflow-mgo.py.

Periodic atomic grids

For routes using an atom-centered periodic grid, each quadrature point normalizes its Becke (Becke, A. D. (1988). *A multicenter numerical integration scheme for polyatomic molecules.* Journal of Chemical Physics, **88**(4), 2547–2553. doi:[10.1063/1.454033](https://doi.org/10.1063/1.454033)) or Stratmann (Stratmann, R. Eric, Scuseria, Gustavo E., and Frisch, Michael J. (1996). *Achieving linear scaling in exchange-correlation density functional quadratures.* Chemical Physics Letters, **257**(3-4), 213–223. doi:[10.1016/0009-2614(96)00600-8](https://doi.org/10.1016/0009-2614(96)00600-8)) weight over the physical atom images within max(becke_image_radius_bohr, 2*d_nearest). Here d_nearest is the distance from that point to its nearest atom image. This makes the image neighborhood independent of which lattice image labels a home atom, and retains coverage in vacuum for diffuse densities. The raw atomic points, weights and ownership keep the selected grid profile.

The default minimum radius remains 10 bohr. It is a truncation control, not a total-energy tolerance: converge it together with radial and angular resolution. For forces, compare more than one displacement step as well; a finite image neighborhood can introduce small boundary steps. A zero radius explicitly requests the molecular partition. These grid checks do not qualify an otherwise unsupported periodic method or gradient route.

2.7 Citations

• C. Pisani, R. Dovesi, C. Roetti, Hartree-Fock Ab Initio Treatment of Crystalline Systems, Lecture Notes in Chemistry 48, Springer, 1988. DOI 10.1007/978-3-642-93385-1. Chapter II.4c gives the periodic quartet bipolar expansion.

• V. R. Saunders, C. Freyria-Fava, R. Dovesi, L. Salasco, C. Roetti, “On the electrostatic potential in crystalline systems where the charge density is expanded in Gaussian functions,” Molecular Physics 77, 629 (1992). DOI 10.1080/00268979200102671. Source for the related shell-penetration and spheropole construction.

• 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. Causa, Y. Noel, “CRYSTAL14: A program for the ab initio investigation of crystalline solids”, Int. J. Quantum Chem. 114, 1287 (2014). DOI 10.1002/qua.24658.

• 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.

3. GDF: native Gaussian density fitting

3.1 Theory

Gaussian density fitting (Whitten 1973, Whitten, J. L. (1973). *Coulombic potential energy integrals and approximations.* Journal of Chemical Physics, **58**(10), 4496–4501. doi:[10.1063/1.1679012](https://doi.org/10.1063/1.1679012); Dunlap 1979, Dunlap, B. I., Connolly, J. W. D., and Sabin, J. R. (1979). *On some approximations in applications of Xα theory.* Journal of Chemical Physics, **71**(8), 3396–3402. doi:[10.1063/1.438728](https://doi.org/10.1063/1.438728); Eichkorn et al. 1995, Eichkorn, Karin, Treutler, Oliver, et al. (1995). *Auxiliary basis sets to approximate Coulomb potentials.* Chemical Physics Letters, **240**(4), 283–290. doi:[10.1016/0009-2614(95)00621-A](https://doi.org/10.1016/0009-2614(95)00621-A); Vahtras, Almlof, Feyereisen 1993) projects the AO pair density onto a finite auxiliary basis \(\{\chi_P\}\):

\[ \rho_{\mu\nu}(\mathbf{r}) = \phi_\mu(\mathbf{r})\phi_\nu(\mathbf{r}) \approx \sum_P d^{\mu\nu}_P\,\chi_P(\mathbf{r}). \]

The fit coefficients minimise the Coulomb error and require only the two-centre auxiliary metric \(M_{PQ} = (P|Q)\) and the three-centre integral \((\mu\nu|P)\). The factorisation

\[ L_{P,\mu\nu} = \sum_Q (\mu\nu|Q) (M^{-1/2})_{QP} \]

reduces the four-centre ERI to a single tensor, and the closed-shell Fock pieces collapse to

\[ J_{\mu\nu} = \sum_L L_{L,\mu\nu}\;\rho_L, \qquad \rho_L = \sum_{\lambda\sigma} L_{L,\lambda\sigma}\,D_{\lambda\sigma}, \]

\[ K_{\mu\nu} = \sum_L \left(L_L\,D\,L_L^\mathrm{T}\right)_{\mu\nu}. \]

This is the standard molecular RIJK kernel that vibe-qc’s DFJKBuilder implements; the periodic version requires three extra ingredients on top:

1. Periodic auxiliary metric \(M^{per}_{PQ} = \sum_T (P_0 | Q_T)\). Naively, the 3D lattice sum of \(1/r\) across a charge-neutral set of aux primitives is conditionally convergent. The fix is PySCF-style modrho charge compensation (Sun, Qiming, Berkelbach, Timothy C., et al. (2017). *Gaussian and plane-wave mixed density fitting for periodic systems.* Journal of Chemical Physics, **147**(16), 164119. doi:[10.1063/1.4998644](https://doi.org/10.1063/1.4998644)): each contracted aux shell is renormalised so its monopole equals \(\sqrt{1/(4\pi)}\), and the AO-pair side is charge-balanced by smooth compensating Gaussians. The compensated integral converges in real space; the missing \((\text{aux}|\text{comp})\) piece is added back analytically in reciprocal space.

2. Periodic three-centre tensor \(T^{per}_{\mu\nu,P} = \sum_T (\mu_0 \nu_0 | P_T)\). Implemented as a libint lattice sum with Schwarz screening (Häser, Marco and Ahlrichs, Reinhart (1989). *Improvements on the direct SCF method.* Journal of Computational Chemistry, **10**(1), 104–111. doi:[10.1002/jcc.540100111](https://doi.org/10.1002/jcc.540100111)), mirroring the molecular compute_3c_eri layout. For each k-point the Bloch-summed \(T^{per}_{\mu\nu,P}(\mathbf{k})\) is the same contraction.

3. Cholesky factor \(L^{per}_{P,\mu\nu} = \sum_Q (\mu_0\nu_0|Q)\,(M^{per})^{-1/2}_{QP}\), computed once per SCF.

vibe-qc’s native GDF stack lives in pbc_gdf.py, periodic_rhf_gdf.py, periodic_gdf_blocks.py, and periodic_k_gdf.py. The full design walkthrough, including the modrho compensation empirical record, is in design_native_gdf.md.

3.2 What ships today

Native GDF in vibe-qc covers:

• Gamma-only RHF, RKS, and hybrids in 3D (a dim=1 wire and a dim=2 slab raise) via run_rhf_periodic_gamma_gdf (and the runner-level run_periodic_job(..., jk_method="gdf")).

• Multi-k RHF / RKS in 3D (a dim=2 slab goes through the bounded slab route; a dim=1 wire on a multi-k mesh fails closed on the current build, see multi-k SCF) via run_krhf_periodic_gdf(..., kmesh, ...) and run_krks_periodic_gdf(..., kmesh, ...). A closed-shell dim=2 slab runs the dedicated slab route, and a dim=1 k-mesh fails closed (no consistent Coulomb gauge for dim < 3). A bulk (1, 1, 1) mesh runs this engine at Nk=1 on the default SR/LR fit. On gdf_method="compcell" or "mdf" it is handed to the exact-Gamma run_pbc_gdf_rhf when that driver implements every requested control, and otherwise also runs at Nk=1 (vibeqc.periodic_gdf_route); only dim < 3 reaches the legacy Gamma driver. Any other Monkhorst-Pack mesh (Monkhorst, Hendrik J. and Pack, James D. (1976). *Special points for Brillouin-zone integrations.* Physical Review B, **13**(12), 5188–5192. doi:[10.1103/PhysRevB.13.5188](https://doi.org/10.1103/PhysRevB.13.5188)) runs the full multi-k loop with k-dependent 3c / 2c blocks, Bloch-phase assembly, and multi-k DIIS (Pulay, Péter (1980). *Convergence acceleration of iterative sequences. The case of SCF iteration.* Chemical Physics Letters, **73**(2), 393–398. doi:[10.1016/0009-2614(80)80396-4](https://doi.org/10.1016/0009-2614(80)80396-4)). Every k-mesh run contracts its Hartree J from the cached fit of its gdf_method; an explicit use_compcell=False is overridden and logged. Pure-functional runs with use_compcell=False used to take an unfitted analytic-FT Ewald Hartree instead. Since 2026-09-29 no k-mesh job reached it (see the explicit gdf_method entry below), and it was removed from this driver (GitLab #415). The analytic-FT J remains the Hartree of the Ewald-3D drivers.

• AO-pair support of the classic fit. On a 3-D cell, fit_source="classic" sums every AO pair whose centres lie within lattice_opts.cutoff_bohr of each other, the term set of the native three-centre integrals (classic_support_policy="physical_pairs", the default since GitLab #466). Before, it summed the cells within cutoff_bohr of the origin for every shell pair. That ball drops the pairs that cross the face of a cell whose extent approaches the cutoff, so a Γ supercell and the k mesh it folds to built different operators: a 15.6-bohr H₈ chain cell at Γ (STO-3G, 3.9-bohr pair period, default 15-bohr cutoff) sat 9.79 mHa per H₂ below its H₂ cell on the \((4,1,1)\) mesh, and now agrees with it to \(3\times10^{-13}\) Ha. The H₂ cell itself moves by 11 nHa. classic_support_policy="historical_index_ball" replays the former operator, and the executed fit source records which one ran. The physical fit is unscreened, so a positive fit_screen_threshold needs the historical replay. The experimental AICCM/CCM and chi producers pin the historical ball, because their own classic folds use it, and the CCM fold refuses a driver cache built on any other support.

• Translation-resolved DF block accessors via compute_2c_eri_lattice_blocks and compute_3c_eri_lattice_blocks, plus the Bloch-sum helpers bloch_sum_2c_eri_blocks / bloch_sum_3c_eri_blocks and the fit helper build_lpq_bloch_native. Summing blocks reproduces the Gamma metric exactly.

• Native gauge sharing. \(V_{ne}\) and \(E_{nn}\) are forced to Ewald-3D on 3D systems regardless of the user’s J/K route choice so the total energy is gauge-consistent (the _gauge_lat_opts_for_v_ne_and_e_nuc helper in periodic_rhf_gdf.py).

• Fermi-Dirac smearing (Mermin, N. D. (1965). *Thermal Properties of the Inhomogeneous Electron Gas.* Physical Review, **137**(5A), A1441–A1443. doi:[10.1103/PhysRev.137.A1441](https://doi.org/10.1103/PhysRev.137.A1441)) at the Gamma point and on multi-k meshes via opts.smearing_temperature (\(k_B T\) in Hartree).

The driver picks one of two backends per call and records it on the PeriodicRHFGDFResult.backend field:

• "native-gamma-gdf": true density fitting via Lpq. Used for 3D molecular-limit cells whose GDF cutoff includes only the home cell.

• "ewald-jk-fallback": Ewald-3D \(J\) + real-space \(K\). Used for 3D systems with image cells inside the GDF cutoff.

run_periodic_job and the k-point GDF drivers call this driver only for dim < 3. A bulk Γ job runs the GDF Hamiltonian of every other mesh (fitted \(J\), exxdiv='ewald' exchange, Sundararaman, R. and Arias, T. A. (2013). *Regularization of the Coulomb singularity in exact exchange by Wigner-Seitz truncated interactions: Towards chemical accuracy in nontrivial systems.* Physical Review B, **87**, 165122. doi:[10.1103/PhysRevB.87.165122](https://doi.org/10.1103/PhysRevB.87.165122)) through run_pbc_gdf_rhf or the k-mesh engine at (1, 1, 1); see vibeqc.periodic_gdf_route.

The SR/LR source builder shares reciprocal AO-pair preparation across groups of at most eight ket k points, then applies each ket’s Bloch phase. This uses the factorization in Eq. 11 of Ye, Hong-Zhou and Berkelbach, Timothy C. (2021). *Fast periodic Gaussian density fitting by range separation.* Journal of Chemical Physics, **154**(13), 131104. doi:[10.1063/5.0046617](https://doi.org/10.1063/5.0046617) without changing the image cutoff, reciprocal mesh or integral tolerance. Fixed-size Fourier workspaces and the caller’s transient cap bound the group size; a tight budget uses smaller groups or the scalar path. Real-space integrals are still built separately for each transfer vector. This optimization makes no change to the one-electron image domain or its convergence bounds.

3.2.1 RSGDF all-FT-Bloch path, chemical-accuracy milestone (v0.9)

For tight ionic 3D crystals (LiH, MgO and similar), the compensated-cell path in vibeqc.pbc_gdf now runs an all-FT-Bloch route in its gdf_method="rsgdf" mode (Ye, Hong-Zhou and Berkelbach, Timothy C. (2021). *Fast periodic Gaussian density fitting by range separation.* Journal of Chemical Physics, **154**(13), 131104. doi:[10.1063/5.0046617](https://doi.org/10.1063/5.0046617)): the 3c tensor is assembled on a dense FFT mesh from Bloch-summed AO pair-FTs, rather than the older sparse-mesh SR + LR split that conflated single-orbital lattice sums (libint short range) with molecular pair-FT long range. The unified all-FT mesh handles both the vacuum-box limit and the tight ionic limit with the same algorithm.

What landed across commit chain 018784c0..d7f4b3bd brings the PBC GDF route to within chemical accuracy on LiH:

System	Pre-session	Post-session
H₂ / 12-bohr / def2-svp-jk	0.6 mHa	0.04 µHa @ ke=400
LiH primitive FCC / sto-3g / def2-svp-jk	+11 Ha (broken)	−2.3 mHa @ ke=200

The three contributing fixes:

1. All-FT-Bloch path (018784c0), dense FFT mesh + Bloch-summed pair-FT for the 3c side, replacing the sparse SR+LR split inside the rsgdf route.

2. modrho-direct fix (c949d748), build_lpq_native_fft now takes the modrho aux basis directly, so the per-L 1/√(4π/(2L+1)) factor baked into the modrho coefficients is not double-applied. Pre-fix Lpq was 330× too large on LiH.

3. Proper Ewald-Madelung (d7f4b3bd, Ewald, Paul Peter (1921). *Die Berechnung optischer und elektrostatischer Gitterpotentiale.* Annalen der Physik, **369**(3), 253-287. doi:[10.1002/andp.19213690304](https://doi.org/10.1002/andp.19213690304)), replaced the cubic Wigner shortcut α_M / L (exact only for cubic Bravais cells, 1.77 % off for primitive FCC) with a proper Ewald self-energy sum over real-space images + reciprocal G-vectors. Bit-exact to PySCF’s pbc.tools.madelung on cubic, FCC, and intermediate cells. See ewald.md for the user- facing surface and python/vibeqc/madelung.py for the implementation.

The remaining −2.3 mHa LiH residue is not a GDF bug, it is located in the V_ne (nuclear-electron) one-electron Hcore: compute_nuclear_lattice_ewald builds V_long via grid quadrature whose Becke partitioning breaks the (1,1,1) 3-fold rotational symmetry of the FCC primitive cell (Li 2p_x ≠ 2p_y ≠ 2p_z by ~1.9 mHa per AO). A prototype analytic FT replacement preserves the symmetry to machine precision; productionising it is the periodic-1e chat’s job, not GDF.

3.2.2 Cauchy-Schwarz screening of the three-centre fit

The multi-k drivers accept fit_screen_threshold (default 0.0 = off, exact) on the gdf_method="rsgdf" path:

r = vq.run_krhf_periodic_gdf(
    system, basis, (2, 2, 2),
    fit_screen_threshold=1e-10,
)

AO pairs \(\mu\nu\) whose Schwarz bound \(\max_P \sqrt{(P|P)}\,\sqrt{(\mu\nu|\mu\nu)}\) on the fit mesh stays below the threshold are dropped from every per-\((k_i,k_j)\) three-centre fit tensor (Neese, Wennmohs, Hansen and Becker, Chem. Phys. 356, 98 (2009), § 3.1, Neese, Frank, Wennmohs, Frank, et al. (2009). *Efficient, approximate and parallel Hartree-Fock and hybrid DFT calculations. A ‘chain-of-spheres’ algorithm for the Hartree-Fock exchange.* Chemical Physics, **356**(1–3), 98–109. doi:[10.1016/j.chemphys.2008.10.036](https://doi.org/10.1016/j.chemphys.2008.10.036): the RI analogue of Häser-Ahlrichs direct-SCF screening). The pair-norm factor is an analytic upper bound, so the screen never underestimates a pair; 1e-10-class thresholds reproduce the unscreened energy to well below SCF accuracy. The kept/dropped pair counts are reported in the SCF log (no silent truncation). The CCM RI routes (run_ccm_rhf_gdf / run_ccm_rks_gdf) forward driver kwargs, so they inherit the option unchanged. Screened fits compose with compute_gradient=True (Γ and multi-k, since 2026-07-30): the analytic gradient differentiates the screened objective at the SCF’s fixed pair mask. Non-rsgdf fit methods raise NotImplementedError rather than silently ignoring the threshold.

The Γ-point drivers (run_pbc_gdf_rhf / run_pbc_gdf_uhf / run_pbc_gdf_uks, and the Γ (1,1,1) fast path of the multi-k drivers) accept the same option. There it additionally selects the memory-lean build: the fit sweeps its reciprocal mesh in G-chunks and the Ewald-3D \(V_{ne}\) Fourier transform streams its own chunked AO-pair FT, so the dense \((n_{\mathrm{AO}}, n_{\mathrm{AO}}, n_G)\) pair-FT tensor (the dominant memory of a Γ-point GDF build on production-basis cells) is never materialised (at the cost of running the pair FT once per consumer instead of once shared).

An explicit fock_mixing / fmixing_percent / level_shift or smearing request on a default-Γ closed-shell job runs the k-mesh engine at (1, 1, 1), which implements it on the same Hamiltonian. The convergence="auto" profile knobs (e.g. the ionic-insulator FMIXING 30% that MgO-class cells resolve) are capability-filtered to zero instead, so the run keeps the faster exact-Γ run_pbc_gdf_rhf evaluation and its DIIS/accelerator stack. The strategy block in the .out states the filtered knob and the reason. Called directly, the exact-Γ drivers raise NotImplementedError for a control they do not implement (run_pbc_gdf_rhf: smearing, Fock mixing, level shift; run_pbc_gdf_uhf / run_pbc_gdf_uks: level shift) instead of dropping it.

3.2.3 Space-group reduction of the exchange build

Exact exchange is the one \(n_k^2\) term in a multi-k GDF SCF: \(K(k_i)\) needs a sum over the whole Brillouin zone for every bra \(k_i\). The space-group saving is therefore on the bra index alone: build \(K(k_i)\) only at the irreducible-wedge representatives and symmetry-transport it across each star (Dovesi, Roberto (1986). *On the role of symmetry in the ab initio Hartree-Fock linear-combination-of-atomic-orbitals treatment of periodic systems.* International Journal of Quantum Chemistry, **29**, 1755–1774. doi:[10.1002/qua.560290608](https://doi.org/10.1002/qua.560290608)). Reducing the ket sum as well is not the same operation and is wrong (measured \(+1.389\) Ha on MgO).

Request it from the production entry point with symmetry_reduce_k:

r = vq.run_periodic_job(
    system, basis,
    method="RHF",
    jk_method="gdf",
    kpoints=(4, 4, 4),
    symmetry="attach",          # the wedge needs a symmetry model
    symmetry_reduce_k=True,
    aux_basis="def2-svp-jk",
    output="nacl_reduced",
)

The Lpq cderi cache then holds \(n_{\mathrm{IBZ}} \times n_k\) blocks plus the \(n_k - n_{\mathrm{IBZ}}\) Hartree diagonals outside the wedge, instead of \(n_k^2\), and the memory preflight is charged at that reduced count. Measured 2026-09-24 on NaCl primitive / def2-SVP with aux_basis="def2-svp-jk" (33 AOs, 240 aux) at \((4,4,4)\), \(n_{\mathrm{IBZ}} = 8\) of 64: 568 stored pairs against 4096, so the Lpq term drops from 15.95 GB to 2.21 GB (7.2x) and the total preflight from 26.32 GB to 5.71 GB (the fit source, the nuclear Fourier workspace, the one-electron lattice blocks and the per-\(k\) buffers do not reduce). The .out states the reduction explicitly (symmetry_reduce_k = exchange bras 64 -> 8 (irreducible wedge)) directly above the [memory] block, so the estimate can be accounted for.

Everything else is untouched (diagonalisation at every \(k\), occupations, the energy expression, and the result shape), so bands / DOS / COOP / QVF consumers see no change. The transported \(K\) equals the one the same fit would build directly only for operations that map the executed exchange source onto itself, so the driver admits operations against that source before any integral is computed:

Exchange source	Finite support	Admitted operations
SR/LR fit (gdf_method="rsgdf" on a 3-D cell, the default)	Pair-centred: AO-pair images \(\lvert A - B - R\rvert \le\) cutoff, auxiliary images by distance from the pair, and the physical \(\lvert q + G\rvert\) sphere	Every operation of the space group, glides and screws included
Classic fit (fit_source="classic", default classic_support_policy="physical_pairs")	Pair-centred: AO-pair images \(\lvert A - B - R\rvert \le\) cutoff_bohr and the same bounded \(\lvert q + G\rvert\) sphere; the classic fit has no auxiliary image domain	Every operation of the space group, glides and screws included
Historical classic replay (classic_support_policy="historical_index_ball")	One origin-centred cell ball \(\lvert R\rvert \le\) cutoff_bohr for every shell pair	Only operations that shift every atom by the same lattice vector (the shift of a glide or screw is computed from its full Seitz operation)
compcell, MDF, a screened classic fit, private cache adapters	Not qualified	None: refused

The star map then uses admitted operations only, and among those only operations that map the whole Monkhorst-Pack mesh onto itself: a finite-mesh exchange is covariant under an operation only when that operation permutes the mesh points, so on an anisotropic mesh of a cubic cell an operation that merely sends one representative onto a mesh point is not used. When the usable operations do not reach every point of the wedge’s stars, the run is refused with the cause named: the source when its admission is the limit (rock-salt LiH on the historical classic replay admits 8 of 48 operations), the mesh when every operation is admitted but the mesh breaks the symmetry (diamond on the SR/LR fit at \((2,3,4)\)). Coverage depends on the mesh as well as the cell: a historical-replay cell whose admitted operations cover the \((2,2,2)\) stars can still be refused at \((3,3,3)\), whose stars need other operations, and diamond on the historical classic replay (12 of 48 operations admitted, 6 of them glide-type) is refused at \((2,2,2)\) and \((3,3,3)\) but resolves at \((2,2,3)\).

Admission also checks the quadrature actually passed to the driver. Every point must match the declared full mesh in native order, modulo reciprocal lattice vectors, and every weight must equal the normalized uniform weight within floating-point roundoff. Reordered, duplicated, displaced or nonfinite points and custom weights are refused before integral work. Finding only the wedge representatives is insufficient: the exchange ket sum uses every point. Malformed array shapes and complex-valued points or weights are refused before input normalization can reshape them or discard their imaginary components. When the input is already an API-produced irreducible mesh, the driver checks its original representatives, weights and ir_mapping against the current system and declared mesh before expanding it. Reciprocal wraps and a consistent reordering of representatives, weights and mapping indices are accepted. Nonfinite or displaced representatives, altered orbit weights and inconsistent mapping entries are refused instead of being discarded by expansion. Rebuild the reduction for the current system or supply its full mesh explicitly. If every orbit is a singleton, the validated reduction is still normalized to native full-grid order before either driver uses it. All points are then computed, so this case saves no exchange rows. This normalization requires a valid reduction mapping; an explicit full mesh without that mapping must already satisfy the full-grid ordering check above. Native BlochKMesh.is_shift and Python KPoints.shift half-step flags are preserved when constructing both full and reduced meshes. Tuple meshes remain Gamma-centred. A shifted mesh must still pass the same source and whole-mesh operation admission; a shift does not override a failed symmetry gate.

Both drivers reserve the aggregate mesh-binding and star-planning numerical storage and work before expanding an input wedge or scanning coordinates. This reservation covers the full point count, both native meshes and matching tables. It excludes Python/container overhead and spglib’s internal workspace; integral and SCF memory admission remain separate. Inconsistent inputs are not reordered, averaged or renormalized. Supply the declared mesh or run with symmetry_reduce_k=False to use the unreduced path.

A glide or screw \(\{R \mid \tau\}\) maps a home-cell atom to \(x_b = R x_a + \tau + L_a\) with \(L_a\) a lattice vector, and in the AO-centred Bloch sum \(\phi_{\mu k} = \sum_T e^{i k \cdot T} \chi_\mu(r - x_a - T)\) it acts as \(e^{i k' \cdot L_a} D(R)\) at \(k' = R k\): the fractional translation is carried entirely by the atom-image shift \(L_a\), and no separate \(e^{-i k \cdot \tau}\) factor exists in this convention (that factor belongs to the atom-centred convention, where it is a global phase and cancels in the transport of an operator). The driver’s transport computes \(L_a\) from the full Seitz operation (Zicovich-Wilson, C. M. and Dovesi, Roberto (1998). *On the use of symmetry-adapted crystalline orbitals in SCF-LCAO periodic calculations. I. The construction of the symmetrized orbitals.* International Journal of Quantum Chemistry, **67**(5), 299–309. doi:[10.1002/(SICI)1097-461X(1998)67:5<299::AID-QUA3>3.0.CO;2-Q](https://doi.org/10.1002/(SICI)1097-461X(1998)67:5<299::AID-QUA3>3.0.CO;2-Q), Eqs. 21-26), so a non-symmorphic group is admitted like any other on the SR/LR source. Measured on diamond \(Fd\bar{3}m\) (48 operations, 24 with a fractional translation) and hcp \(P6_3/mmc\) (24, 12) with explicit 8-bohr SR/LR sources, pure s/p/d shells and nonsymmetric complex trial densities, \(K\) transported through every operation, proper and improper, glide and screw, matches the directly built \(K\) to \(6.7\times10^{-15}\) and \(2.0\times10^{-15}\); dropping the per-atom phases misses by \(0.37\) and \(3.9\times10^{-2}\) (tests/test_periodic_k_gdf_ibz_nonsymmorphic.py). Before this admission the same case iterated and was then refused as a “variational symmetry break” of a closed-shell state, because \(J\) and \(K\) from the historical classic fit are not covariant under the other 40 operations. Measured on a two-atom \(Pm\bar{3}m\) cell with an explicit 8-bohr finite source and a nonsymmetric complex trial density, \(K\) transported through each of the 48 operations matches the directly built \(K\) to \(3.6\times10^{-15}\) relative on the SR/LR fit and to \(2.3\times10^{-15}\) on the physical classic fit. On the historical classic replay the 6 shift-constant operations match to \(1.8\times10^{-15}\) and the 42 others miss by \(1.3\times10^{-2}\) to \(1.8\times10^{-2}\) (tests/test_periodic_k_gdf_ibz_admission.py, re-measured 2026-10-06). On diamond and hcp the physical classic fit transports through every operation, glides and screws included, to \(8.5\times10^{-15}\) and \(2.1\times10^{-15}\) (tests/test_periodic_k_gdf_ibz_nonsymmorphic.py). The admitted transport is therefore exact for its finite approximation, independent of cutoff_bohr; the former cell-list warning is gone.

The electronic state is a separate precondition. Every SCF iteration compares the density with its transported wedge (each star member against its representative, and each representative against its own image under every admitted operation that fixes it, time reversal included), against a floor calibrated on the overlap, and refuses a density that breaks the wedge’s symmetry either way. A state that broke its little group at a representative while being propagated consistently to the star members used to pass the member comparison at roundoff; it is now refused with the little group named. That gate is a numerical threshold, not a proof of state invariance.

The occupations of a wedge-native run are filled once per star: the zero-temperature or smearing fill is evaluated on the wedge representatives with each star’s total weight and copied to the star members, so star-related \(k\) points carry identical occupations by construction. Symmetry-degenerate levels are not exactly degenerate in a periodic calculation (the one-electron lattice sums split them by about \(10^{-12}\) Ha on LiH, a decade above the fill’s roundoff tolerance), so at each representative any levels split by no more than the measured star spectral spread, or by less than one nanohartree, are treated as one degenerate group and share their electrons equally (the nanohartree floor covers a little-group-degenerate pair at a \(k\) point that is its own star, which has no partner to measure its noise on). A star-degenerate frontier therefore receives the symmetric ensemble occupation (the triplet HCORE guess on rock-salt LiH at \((2,2,2)\) puts \(1/6\) on each of its six L-point frontier states, where the plain fill put \(1/2\) on two of them and was refused by the gate at a residual of \(0.51\)). The full-BZ run keeps the plain fill, so at such a frontier the two runs are different states by design: the wedge-native one is the symmetry-adapted ensemble, the full-BZ one whichever basin the noise selected. bz_integration="gilat" cannot be combined with the wedge-native run.

Every precondition fails closed rather than silently running full-BZ, because a performance flag that looks effective and is not is worse than a refusal:

• jk_method="gdf", dim=3, a true multi-k Monkhorst-Pack mesh, and method in RHF / RKS / UHF / UKS (spin-restricted open-shell GDF has no wedge surface);

• an attached symmetry model: pass symmetry="attach", or attach validated operations to the PeriodicSystem yourself;

• exact exchange, because on a pure functional there is no \(n_k^2\) term to reduce;

• an exchange source whose admitted operations cover every star, as tabulated above;

• no geometry optimization: the wedge-native SCF converges the k-transported exchange, whose gradient is not the captured single-point objective.

The underlying driver keyword is ibz_native=True on run_krhf_periodic_gdf / run_krks_periodic_gdf / run_kuhf_periodic_gdf / run_kuks_periodic_gdf, for callers driving the drivers directly.

3.3 Advantages

• Best closed-shell scaling. With aux size \(M \approx 3N\), per-iteration cost is \(\mathcal{O}(N^2 M) = \mathcal{O}(N^3)\), compared to \(\mathcal{O}(N^4)\) for the direct-K piece in BIPOLE.

• Production multi-k for closed-shell RHF / RKS in 3D, plus a verified closed-shell energy + analytic-gradient route in 2D. Slab GDF uses a separate signed truncated metric (Spencer, James and Alavi, Ali (2008). *Efficient calculation of the exact exchange energy in periodic systems using a truncated Coulomb potential.* Physical Review B, **77**, 193110. doi:[10.1103/PhysRevB.77.193110](https://doi.org/10.1103/PhysRevB.77.193110)) and rigorous V_ne(k) / nuclear gauge; it is not the old bulk builder with a normal-axis vacuum cell. jk_method="auto" still picks the direct SLAB_EWALD_2D route for dim=2; request "gdf" explicitly.

• Hybrid DFT included. Exact-exchange coefficients are contracted out of the same Lpq factor that builds \(J\).

• Same kernel as molecular RIJK. Auxiliary-basis selection, Cholesky robustness, and linear-dependence handling are all shared with the molecular DFJKBuilder.

3.4 Disadvantages

• No 1-D wire. The bare 1-D \(V_{ne}\) / \(E_{nn}\) lattice sums are only conditionally convergent, and the neutral three-centre fit collapses to a transverse-uniform sheet term on a vacuum-padded cell, so the three Coulomb-bearing pieces cannot share one gauge below three dimensions. The GDF drivers therefore raise NotImplementedError for dim=1, at Gamma and with a k mesh, instead of returning a cutoff-dependent number. Model the wire as a dim=3 cell with vacuum on the two transverse axes and converge the energy with respect to that vacuum; the cost grows with the box volume.

• AUTO routing is dimensional. In 3D, AUTO keeps BIPOLE as the UHF/UKS default and explicit jk_method="gdf" opts into native open-shell GDF. In 1D, AUTO raises NotImplementedError before any setup, for every method, and names the dim=3 vacuum-box model: no public J/K builder has a consistent 1-D Coulomb gauge, and vibe-qc does not rewrite a dim=1 system as dim=3 for you, because the vacuum is part of the model. An open-shell wire has no route: the open-shell GDF drivers accept only dim=3 at Gamma, and raise on a dim=1 k mesh. Native open-shell GDF in 3D ships at Gamma (run_pbc_gdf_uhf / run_pbc_gdf_uks, including per-spin Fermi-Dirac smearing) and multi-k (run_kuhf_periodic_gdf / run_kuks_periodic_gdf).

• The analytic periodic GDF gradient has route-specific envelopes. The original compute_gdf_gradient implementation covers 3D Gamma gdf_method="compcell" for closed-shell RHF and RKS (run_pbc_gdf_rhf(functional=..., compute_gradient=True), milestone 2: general-\(a_x\) J/K machinery plus the lattice-summed XC Pulay primitive on the SCF’s own quadrature) and, since milestone 3b (2026-07-29), UHF (run_pbc_gdf_uhf(compute_gradient=True)), on both AFT settings: the 2c/3c AFT reciprocal-space corrections are differentiated analytically, so the production apply_aft_correction=True fit no longer needs to be disabled for forces. Its Ewald \(V_{ne}\) derivative is fully analytic, it reuses the SCF’s resolved 2c/3c cutoffs, and it differentiates the threshold-truncated fit eigenspace. The same milestone identified and fixed the long-standing flat \(3.4\times10^{-3}\) Ha/bohr FD residual (a truncated-direct-sum nuclear gradient paired with the Ewald-gauge \(E_{nn}\)): the end-to-end H2/ STO-3G FD regressions now hold every route (RHF, LDA/PBE/PBE0, and the UHF triplet) to \(10^{-6}\) Ha/bohr (observed \(\sim 10^{-8}\)). UKS followed as milestone 4 (run_pbc_gdf_uks(functional=..., compute_gradient=True)): the per-spin assembly takes the functional’s exact-exchange fraction and adds the spin-polarised XC Pulay; UKS(M=1) collapses onto the RKS gradient at machine precision and the triplet passes the same AFT-on FD gate. The production rsgdf fit gained the same full Gamma ladder as milestone 6 (gdf_method="rsgdf", compute_gradient=True, RHF/RKS/UHF/UKS): its reciprocal-space metric and 3c tensor carry atom positions only in exp(-iG.R) phases, so their weighted centre derivatives are exact, and the full-SCF FD sits at ~4e-9 even at the ke=200 production default. Since the open-shell drivers default to rsgdf, run_pbc_gdf_uhf/uks(compute_gradient=True) needs no extra knobs. Tailed dense-core rsgdf fits are differentiated too: the long-standing dense-core analytic-vs-FD mismatch (~2.5e-5 Ha/bohr on MgO) was an e_nuc real-space truncation artefact in the energy, fixed 2026-07-29 by switching the 3D GDF drivers to the converged Ewald nuclear term. The molecular-grid Gamma compcell RKS/UKS path (use_periodic_becke=False) includes grid-point and partition-weight motion through the existing fixed-density XC central difference. This completes the XC contribution, but a separate compcell Coulomb derivative residual remains for some near-null kept metric modes: the recorded LiH/STO-3G/PBE 6-bohr cube retains about \(1.2\times10^{-5}\) Ha/bohr after the grid response, above the \(10^{-6}\) full-gradient target. The periodic-Becke route and its independent acceptance are unchanged. Multi-k GDF gradients landed 2026-07-30 (G-PBC-002 Item 4): run_krhf_periodic_gdf / run_krks_periodic_gdf / run_kuhf_periodic_gdf / run_kuks_periodic_gdf accept compute_gradient=True and return the analytic nuclear gradient on result.gradient, differentiating the shared-q rsgdf fit k-resolvedly (per-q metric/3c centre derivatives, DF-K over every momentum-transfer group, the exxdiv W-shift, and the XC response on the SCF’s own inverse-Bloch folded density). The KS response combines the analytic fixed-grid AO/GGA derivative with a fixed-density central difference of the moving atom-centred quadrature points and Becke weights; this is an XC-only correction, not a 6N-SCF gradient. Full-SCF FD gates on H2 (2,1,1) at h=2e-4: KRHF \(1.1\times10^{-8}\), KUHF triplet \(1.6\times10^{-8}\), KRKS LDA \(4.1\times10^{-10}\) / PBE0 \(1.5\times10^{-8}\) Ha/bohr. A genuine-spin LiH triplet/PBE KUKS gate matches FD to about \(7\times10^{-9}\) Ha/bohr and conserves net force, with the KUKS(M=1) collapse onto KRKS retained at machine precision. Supported multi-k envelope: 3D cells, gdf_method="rsgdf" on the cached-Lpq SCF (use_compcell=True; every 3D rsgdf run promotes to it, pure DFT included), Γ-centered uniform full-BZ meshes, gapped integer T = 0 Aufbau occupations, pure and global-hybrid functionals. A band-overlap T = 0 global-Aufbau ensemble is supported for energies but fails closed for analytic gradients until its sharp-Fermi response is validated. The ASE backend="gdf" calculator routes kpts != (1,1,1) through these drivers (rsgdf only; the tuple kpts is Γ-centered there, unlike the ewald backend’s shifted Monkhorst-Pack on even meshes). Schwarz-screened fits are differentiated too since 2026-07-30 (Γ and multi-k): the gradient cache re-derives the SCF’s per-q pair mask from the same mask builder and differentiates the screened objective at that FIXED mask (masked pairs are hard zeros of the fit, so their derivative is exactly zero; the mask itself is a discrete quantity with no derivative). Screened full-SCF FD gates: Γ H2 \(2.7\times10^{-9}\), KRHF (2,1,1) \(1.2\times10^{-8}\) Ha/bohr. Finite-temperature (Fermi-Dirac) smearing is differentiated too since 2026-07-30 (Γ open-shell + multi-k): the analytic gradient is \(dA/dR\) of the driver’s reported Mermin free energy \(A = E - T S\): by Mermin’s variational principle (Phys. Rev. 137, A1441 (1965)) the occupation- and \(\mu\)-response terms vanish at self-consistency, so the \(T = 0\) assemblers run on the fractional-occupation density and energy-weighted density (the standard smeared-force statement, Marzari & Vanderbilt, PRL 82, 3296 (1999), Marzari, Nicola, Vanderbilt, David, et al. (1999). *Thermal Contraction and Disordering of the Al(110) Surface.* Physical Review Letters, **82**(16), 3296–3299. doi:[10.1103/PhysRevLett.82.3296](https://doi.org/10.1103/PhysRevLett.82.3296)). Smeared full-SCF FD gates against the free energy: Γ UHF \(3.8\times10^{-10}\) (compcell) / \(1.9\times10^{-10}\) (rsgdf), KRHF (2,1,1) \(3.5\times10^{-10}\) Ha/bohr at genuinely fractional occupations (\(f \approx 1.86/0.15\)), plus a \(T \to 0\) consistency pin at \(\le 10^{-15}\). FD-differencing result.energy instead of result.free_energy disagrees by the \(T\,dS/dR\) term (~0.1 Ha/bohr on the gate fixture); the free energy is the force surface. bz_integration="smearing" is accepted (it is literally the same occupation path). Restriction: Fermi-Dirac only for now: the MP/MV generalized entropies are conjugate to their occupations, but their GDF force paths do not yet have dedicated full-SCF gates, and MP occupations can be negative while the present fractional-block assembly takes square roots. Fixed-geometry smeared GDF gradients remain available and differentiate the free energy. High-level smeared GDF optimization now fails closed: although the force is the matching \(dA/dR\), the optimization, trajectory, and QVF schemas do not yet distinguish the Mermin objective \(A\) from internal energy \(E\). Slab (dim=2) GDF gradients landed 2026-07-30 (G-PBC-002 § 6): on a dim=2 system, run_krhf_periodic_gdf / run_krks_periodic_gdf with compute_gradient=True assemble the analytic gradient of the slab route’s two-gauge composition: the bare Parry / de Leeuw-Perram 2D-Ewald nuclear term (Parry, D. E. (1975). *The electrostatic potential in the surface region of an ionic crystal.* Surface Science, **49**, 433–440. doi:[10.1016/0039-6028(75)90362-3](https://doi.org/10.1016/0039-6028(75)90362-3), de Leeuw, S. W. and Perram, J. W. (1979). *Electrostatic lattice sums for semi-infinite lattices.* Molecular Physics, **37**(4), 1313–1322. doi:[10.1080/00268977900100951](https://doi.org/10.1080/00268977900100951)), the slab \(V_{ne}\) derivative (z-resolved reciprocal blocks included), the S/T lattice folds, the signed Sundararaman-Arias truncated-metric DF-J/K fit derivative (the indefinite metric’s Fréchet response, finite negative \(K(0)\) zero mode and all), and the BvK probe-charge exxdiv overlap response, plus the periodic-Becke XC Pulay for RKS. Full-SCF FD on the compact H2 (2,2,1) slab fixture at h=2e-4: KRHF \(\le 2.2\times10^{-10}\) Ha/bohr per component (gate 1e-6, slab normal included; Γ (1,1,1) reduction \(2.1\times10^{-10}\)). KRKS is gated at \(2\times10^{-4}\): the XC term is exact at fixed grid (\(9\times10^{-12}\), asserted), and the full-SCF residual (PBE/PBE0 \(\sim 10^{-5}\) to \(8\times10^{-5}\) on that deliberately compact cell) is the moving periodic-Becke partition’s becke_image_radius_bohr truncation derivative; it falls to \(1.4\times10^{-6}\) at a 20-bohr image radius and is the same quadrature-limited class as the documented Γ KS gates. run_periodic_job(optimize=True, jk_method="gdf") on a slab relaxes on this objective through the same optimizer capture as the bulk routes; optimize_cell stays fail-closed (no slab GDF analytic stress). Slab open shell, smearing, custom meshes, and restarts remain fail-closed on the slab route itself; the ASE calculator still accepts only fully periodic (pbc=(True, True, True)) cells, so ASE slab forces remain closed at the boundary. Not yet differentiated (pinned pre-dry-run errors naming the reason): MDF/compcell multi-k fits, non-Fermi-Dirac smearing flavors, bz_integration="gilat", custom-weight k-point sets, pure-functional bulk multi-k RKS without the cached-Lpq fit, range-separated functionals, k_exchange="cosx", DFT+U, wire (dim=1) meshes. Optimizer wiring landed 2026-07-30: run_periodic_job(optimize=True) with a GDF-routed SCF (Γ and multi-k) relaxes on this GDF analytic-gradient objective – the exact driver call the SCF dispatched, re-run with compute_gradient=True at each candidate geometry – rather than switching to the BIPOLE force surface, so the relaxed geometry is a stationary point of the energy the job reports. The .out optimization block states the active force objective; unsupported gradient envelopes fail during runner preflight before dry-run/SCF, and periodic dispersion and positive-temperature requests fail before SCF because the optimizer does not add D3 per geometry and cannot yet label a Mermin free-energy objective distinctly; optimize_cell=True fails closed on the GDF and BIPOLE routes because neither has a certified production stress and coupled cell optimizer. Periodic GDF KS also fails closed on VV10/nonlocal-correlation and double-hybrid functionals because those drivers do not add the nonlocal VV10 potential or the perturbative MP2 correlation term. Molecular RIJK has separate analytic-gradient coverage; see density_fitting.md.

• gdf_method="compcell" on tight ionic cells (fixed 2026-09-29, GitLab #172). Until then the Gamma-only compcell Hartree converged to non-physical fixed points on compact cells with Z > 1 atoms (rocksalt LiH/STO-3G: +579.8 Ha at HF, +1172.6 Ha at UKS-PBE, against PySCF’s -8.24 Ha), and the drivers warned and tagged the result. Four defects combined: the Gaussian three-centre integrals kept the auxiliary function in the home cell, \(\sum_R (P_0|\mu_0\nu_R)\), while the periodic tensor needs its images, \(\sum_{T,R} (P_T|\mu_0\nu_R)\) (Sun et al. 2017, Eq. 15); the finite \(G = 0\) term of the compensated real-space integrals was kept; the reciprocal-space correction scaled p-AO pairs wrongly; and the Gamma driver used a flat one-electron domain for compcell and MDF. With all four fixed, compcell lands 6e-8 Ha from PySCF GDF on H₂/STO-3G in a 6-bohr cube (before: 162 mHa), 4.8 µHa on LiH/STO-3G and 27 µHa on MgO/STO-3G, where the remaining spread is the choice of near-null metric modes (the SR/LR fit moves by 3.5 µHa between cuts there). MDF lands 2e-8 Ha from PySCF MDF on the H₂ cube and 0.9 µHa on LiH. The multi-k builders now resolve the ket k point of each AO pair as well as the momentum transfer. The open-shell Gamma drivers have defaulted to gdf_method="rsgdf" since 2026-07-09, and run_pbc_gdf_rhf since 2026-09, like every other GDF entry point. A KS run through run_pbc_gdf_uks with options=None defaults to PeriodicKSOptions, so the periodic-Becke-grid/torus-density XC pairing engages by default (the old PeriodicRHFOptions fallback evaluated XC in the v0.8.x molecular-grid convention: -7.9638 Ha on the LiH fixture).

• gdf_method="mdf" on compact dense-core cells (admitted since 2026-09-29, GitLab #172). From 2026-07-29 all five GDF drivers refused the class (Z >= 8 in < 250 bohr³/atom): on MgO primitive FCC/STO-3G, Gamma MDF was non-convergent and trial-dependent at mdf_ke_cutoff 40/60 (+1703/+2143/+418 Ha), falsely converged at 80 to -1651 Ha, and multi-k (1,1,2) landed at -11741 Ha; later, with the metric built accurately, it was bounded but 92 mHa off. Both failures were the auxiliary images its Gaussian three-centre block omitted. With them summed, and the plane-wave-dressed metric’s floor recalibrated from 1e-4 to 1e-5 of its largest eigenvalue, MgO/STO-3G lands 3.4 µHa from PySCF MDF at Gamma and 70 µHa at (2,2,2); LiH/STO-3G 0.9 µHa and 2.8 µHa. MDF on all-electron cores depends on the plane-wave mesh by about 1e-4 Ha, in PySCF as in vibe-qc (Ne/STO-3G in a 10-bohr box moves by -65 µHa from mdf_ke_cutoff=60 to 80, and PySCF’s own MDF by -130 µHa from its default mesh to 50³), so compare MDF results at a fixed mdf_ke_cutoff. The default gdf_method="rsgdf" stays the production route on these cells (+2.9 µHa from PySCF GDF on MgO).

• An explicit gdf_method on a pure-DFT k-mesh job (fixed 2026-09-29, found on GitLab #172). A closed-shell pure-functional dim=3 k-mesh run with gdf_method="mdf" or "compcell" and the default use_compcell=False kept the Ewald-3D Hartree, which builds no density fit, so the requested method was accepted and silently ignored. This covered run_krks_periodic_gdf and run_periodic_job(method="rks", gdf_method=...) with a k-mesh. It now builds the cached fit with the requested method, as the default "rsgdf" already did, and logs the promotion. On H₂/STO-3G in a 6-bohr cube at (2,1,1) with PBE, both methods used to return the Ewald-J energy -1.1995011057 Ha. They now return -1.1995012305 Ha (mdf) and -1.1995485915 Ha (compcell, 0.13 µHa from rsgdf). The open-shell drivers always built the requested fit and are unchanged.

• rsgdf high-|G| tail on dense-core cells. The default gdf_method="rsgdf" on a dim=3 cell uses the SR/LR fit source at Gamma and on every k-mesh, closed- and open-shell. It resolves the tight-core AO products without an exponent-scaled reciprocal tail (MgO/STO-3G: 2.9 µHa from PySCF GDF at Gamma, 3.2 µHa at (1,1,2)), is never tagged +PARITY_HELD, and ignores rsgdf_tail_ke_cutoff with a DeprecationWarning (rsgdf_g_precision controls its raw integral error). The tail survives on the classic multi-k source, fit_source="classic". Untailed on the dense-core class it sits -0.50 Ha from PySCF on MgO/STO-3G, so the multi-k drivers warn and tag +PARITY_HELD until rsgdf_tail_ke_cutoff reaches 10x the steepest AO primitive exponent, at n_k^2 extra pair builds. A tail on a multi-k gdf_method="mdf" or "compcell" run raises, because only the rsgdf builders carry one. Before 2026-09 a pure-DFT multi-k job took an unfitted Ewald-3D Hartree, which accepted the knob and silently ignored it (IID 146) and later auto-sized it (IID 518); that Hartree was removed from the k-mesh driver (GitLab #415).

• Auxiliary basis coverage. Standard def2-*-jkfit aux bases (Weigend, Florian (2008). *Hartree-Fock exchange fitting basis sets for H to Rn.* Journal of Computational Chemistry, **29**(2), 167–175. doi:[10.1002/jcc.20702](https://doi.org/10.1002/jcc.20702)) work cleanly. pob-tzvp does not yet have a matching aux: designing a pob-tzvp-jk is a paper-worthy item on its own; the pragmatic workaround is def2-tzvp-jk (aux size scales with orbital size, not orbital identity) or falling back to BIPOLE on pob-TZVP work.

• Aux conditioning is real. Diffuse aux primitives produce a near-singular metric; vibe-qc uses pivoted Cholesky (Lehtola, Susi (2019). *Curing basis set overcompleteness with pivoted Cholesky decompositions.* Journal of Chemical Physics, **151**(24), 241102. doi:[10.1063/1.5139948](https://doi.org/10.1063/1.5139948)) with a linear_dep_thr knob (default 1e-9) and the modrho renormalisation, but you can still see warnings on very diffuse composite bases. Cut diffuse exponents with drop_eta if needed.

Mixed density fitting (MDF) and the compact-cell fit

gdf_method selects how the periodic three-index tensor is fitted, on run_periodic_job(..., jk_method="gdf", gdf_method=...) and on every GDF driver:

• "rsgdf" (the default on every entry point since 2026-09): range-separated GDF. The Coulomb kernel is split into an erfc part fitted analytically in real space and an erf part evaluated on a reciprocal mesh, so no lattice sum is truncated in the fit. The SR/LR fit source used at Gamma and on every k mesh resolves the tight-core AO products without an exponent-scaled reciprocal tail; the tail knob survives only on the classic multi-k source.

• "compcell": the compensated-cell Gaussian fit. Each auxiliary function carries a compensating charge so that the real-space two- and three-centre Coulomb sums converge, and a reciprocal-space correction restores the long-range part; the compensating width is compcell_eta.

• "mdf": mixed density fitting (Sun, Berkelbach, McClain and Chan 2017, Sun, Qiming, Berkelbach, Timothy C., et al. (2017). *Gaussian and plane-wave mixed density fitting for periodic systems.* Journal of Chemical Physics, **147**(16), 164119. doi:[10.1063/1.4998644](https://doi.org/10.1063/1.4998644)). Plane waves are added to the Gaussian auxiliary set: the plane-wave block (up to mdf_ke_cutoff, default 40 Ha) resolves the smooth long-range density and the Gaussians the cusp region, and the plane-wave-dressed metric is orthogonalised under a linear-dependence threshold (gdf_linear_dep_threshold, an absolute cut, PySCF’s convention).

What #172 fixed (2026-09-29). Both Gaussian fits kept the auxiliary function in the home cell, \(\sum_R (P_0|\mu_0\nu_R)\), where the periodic tensor needs its images, \(\sum_{T,R}(P_T|\mu_0\nu_R)\) (Sun et al. 2017, Eq. 15). With the images summed, compcell lands 6e-8 Ha from PySCF GDF and MDF 2e-8 Ha from PySCF MDF on H₂/STO-3G in a 6-bohr cube (before: 162 mHa and 33 mHa), and 4.8 µHa and 0.9 µHa on LiH/STO-3G. compcell also drops the \(G = 0\) term of its compensated real-space integrals, as MDF and PySCF do, so its energy no longer depends on compcell_eta (on H₂ in a 16-bohr cube it moved by +0.6 mHa onto the SR/LR fit); its reciprocal-space three-centre correction now applies the solid-harmonic factor to the AO-pair transform, so cells with p or higher AOs are right (MgO/STO-3G: +3.2e5 Ha before, 0.03 mHa from PySCF GDF after); the exact-Gamma drivers resolve the automatic one-electron AO domain for every fitting method (compcell on LiH was 4.7 mHa off its own (1,1,1) run); and the analytic Gamma compcell gradient differentiates the auxiliary images and the \(G = 0\) term, matching finite differences to 2e-9 Ha/bohr on the 6-bohr cube. The multi-k builders resolve the ket k point of each AO pair as well as the momentum transfer. Measured on the vibe-qc 0.17.x development build (Apple M2, 41 s for all three runs):

import numpy as np
import vibeqc as vq

# H2/STO-3G in a 6-bohr cube: the #172 fixture.
system = vq.PeriodicSystem(
    dim=3, lattice=np.eye(3) * 6.0,
    unit_cell=[vq.Atom(1, [2.3, 3.0, 3.0]), vq.Atom(1, [3.7, 3.0, 3.0])],
)
basis = vq.BasisSet(system.unit_cell_molecule(), "sto-3g")

for fit in ("rsgdf", "compcell", "mdf"):
    r = vq.run_periodic_job(
        system, basis, method="RHF", jk_method="gdf", gdf_method=fit,
        output=f"h2-cube-{fit}",
    )
    print(f"{fit:9s} E = {r.energy:.8f} Ha  converged = {r.converged}")
# rsgdf     E = -1.23855629 Ha  converged = True
# compcell  E = -1.23855635 Ha  converged = True
# mdf       E = -1.23853107 Ha  converged = True

MDF on compact dense-core cells. Since 2026-09-29 (GitLab #172) the compact dense-core class (Z >= 8 in less than 250 bohr³ per atom: MgO, NaCl and the like) is admitted for gdf_method="mdf". From 2026-07-29 it had been refused, because the Gaussian three-centre block omitted the auxiliary images; with the images summed and the floor of the plane-wave-dressed metric recalibrated from 1e-4 to 1e-5 of its largest eigenvalue, MgO/STO-3G lands 3.4 µHa from PySCF MDF at Gamma and 70 µHa at (2,2,2), and LiH/STO-3G 0.9 µHa and 2.8 µHa. MDF on all-electron cores depends on the plane-wave mesh by about 1e-4 Ha, in PySCF as in vibe-qc, so compare MDF results at a fixed mdf_ke_cutoff. The default "rsgdf" stays the production route on these cells (+2.9 µHa from PySCF GDF on MgO). The behaviour of the classic fit source on the dense-core class, which still warns and tags the backend +PARITY_HELD on some multi-k routes, is described in the list of known limitations above.

3.5 When to use GDF

• You are doing closed-shell RHF / RKS on a dim=3 bulk crystal, or a closed-shell dim=2 slab on a full Gamma-centered tuple mesh (energies and analytic forces). A wire is run as a dim=3 cell with transverse vacuum, since the GDF drivers raise on dim=1.

• You want the cheapest production route on bulk for hybrid DFT and large basis (def2-SVP through def2-TZVP class).

• You want an open-shell UHF / UKS single-point on the density-fitting metric (Gamma or multi-k) in 3D. Pass jk_method="gdf" explicitly because AUTO routes open-shell to BIPOLE. There is no open-shell route for a dim=1 wire.

• You want analytic forces on the periodic path: the GDF gradient ladder covers Gamma, multi-k, and closed-shell dim=2 slabs.

3.6 Quick start

import vibeqc as vq

# 2D graphene sheet, RKS / PBE, Gamma-only.
# AUTO resolves a dim=2 slab to the direct SLAB_EWALD_2D route. Select GDF
# explicitly only for the bounded closed-shell fitted energy route.
sysp  = build_graphene()                       # dim=2 PeriodicSystem
basis = vq.BasisSet(sysp.unit_cell_molecule(), "def2-svp")

result = vq.run_periodic_job(
    sysp, basis,
    method="RKS", functional="pbe",
    jk_method="gdf",                           # signed slab-truncated fit
    kpoints=(2, 2, 1),
    output="graphene-pbe",
)
print(f"E_total = {result.energy:+.6f} Ha / cell")

Multi-k closed-shell example:

import vibeqc as vq
from vibeqc import monkhorst_pack
from vibeqc.periodic_k_gdf import run_krhf_periodic_gdf

sysp  = build_mgo()
basis = vq.BasisSet(sysp.unit_cell_molecule(), "def2-svp")
kmesh = monkhorst_pack(sysp, [4, 4, 4])

opts = vq.PeriodicRHFOptions()
opts.lattice_opts.cutoff_bohr = 12.0
opts.use_diis = True

result = run_krhf_periodic_gdf(sysp, basis, kmesh, opts,
                               aux_basis="def2-svp-jk")
print(f"E_total = {result.energy:+.6f} Ha   k-mesh = 4x4x4")

The full runner-level workflow is at examples/periodic/input-mgo-rocksalt-rhf.py and the unified-runner variants under examples/periodic/input-mgo-*.

3.7 Citations

• J. L. Whitten, “Coulombic potential energy integrals and approximations”, J. Chem. Phys. 58, 4496 (1973).

• B. I. Dunlap, J. W. D. Connolly, J. R. Sabin, “On some approximations in applications of X-alpha theory”, J. Chem. Phys. 71, 3396 (1979).

• O. Vahtras, J. Almlof, M. W. Feyereisen, “Integral approximations for LCAO-SCF calculations”, Chem. Phys. Lett. 213, 514 (1993).

• K. Eichkorn, O. Treutler, H. Ohm, M. Haser, R. Ahlrichs, “Auxiliary basis sets to approximate Coulomb potentials”, Chem. Phys. Lett. 240, 283 (1995).

• F. Weigend, “Accurate Coulomb-fitting basis sets for H to Rn” (the def2/J Coulomb-fitting sets), Phys. Chem. Chem. Phys. 8, 1057 (2006). DOI 10.1039/b515623h.

• Q. Sun, T. C. Berkelbach, J. D. McClain, G. K.-L. Chan, “Gaussian and plane-wave mixed density fitting for periodic systems”, J. Chem. Phys. 147, 164119 (2017). The reference algorithm for the modrho-compensated periodic GDF that vibe-qc reproduces.

3.8 Periodic ROHF and ROKS: restricted open-shell SCF

For open-shell systems that need a single set of doubly-occupied spatial orbitals plus a defined set of singly-occupied orbitals (rather than the independent α / β orbitals of UHF), vibe-qc ships periodic restricted open-shell Hartree-Fock drivers on the Ewald-3D gauge and Gamma-point GPW. The Gamma-point GPW route also provides the corresponding spin-pure ROKS method with periodic XC:

import vibeqc as vq
from vibeqc import run_periodic_rohf_gpw, run_periodic_roks_gpw
from vibeqc.periodic_rohf_ewald import run_rohf_periodic_gamma_ewald3d
from vibeqc.periodic_rohf_multi_k_ewald import run_rohf_periodic_multi_k_ewald3d

run_rohf_periodic_gamma_ewald3d runs the Gamma-point ROHF SCF; run_rohf_periodic_multi_k_ewald3d runs the multi-k variant over a Monkhorst-Pack mesh. Both share the Ewald-3D one-electron gauge and J / K build with the UHF Ewald path, so the total energy is gauge-consistent with the other 3D routes.

run_periodic_rohf_gpw is the 3D Gamma-point HF alternative. It combines the smooth-grid GPW Hartree J with per-spin exact K and the periodic Ewald one-electron/nuclear gauge, then uses the shared Roothaan loop (Roothaan, Clemens C. J. (1960). *Self-Consistent Field Theory for Open Shells of Electronic Systems.* Reviews of Modern Physics, **32**(2), 179–185. doi:[10.1103/RevModPhys.32.179](https://doi.org/10.1103/RevModPhys.32.179)) to produce one restricted spatial-orbital set with integer closed/open/virtual occupations of 2/1/0. The same maintained-preview route is available through vq.run_periodic_job(system, basis, method="ROHF", jk_method="gpw").

The corrected-Ewald-exchange multi-k engine is also wired into the public runner on the BIPOLE route: vq.run_periodic_job(system, basis, method="ROHF", jk_method="bipole") (the AUTO default when only method="ROHF" is given) runs run_rohf_periodic_multi_k_ewald3d at Gamma (a (1,1,1) mesh) or on a full Monkhorst-Pack mesh via kpoints=[n1, n2, n3], with the ordinary .out / .system / QVF / citation output. The LiH+/STO-3G (3,1,1) anchor of that engine is validated out of process against PySCF KROHF/GDF to 0.062 mHa. BIPOLE-specific knobs the engine does not honour (MOM, ODA, multipole far-field, ewald_omega overrides, sr_image_precision, Fock mixing) fail closed rather than being silently ignored.

run_periodic_roks_gpw(..., functional="pbe") uses the same restricted orbital and Roothaan contracts while adding the established spin-polarised GPW XC potential and energy. It is also available through vq.run_periodic_job(system, basis, method="ROKS", functional="pbe", jk_method="gpw"). LDA, GGA, meta-GGA, and global-hybrid functionals share the Gamma GPW UKS envelope. Range-separated and direct double-hybrid ROKS are not implemented on this backend.

ROKS also runs on the BIPOLE route, on the same corrected-Ewald-exchange multi-k engine as ROHF above: vq.run_periodic_job(system, basis, method="ROKS", functional="pbe", jk_method="bipole") reaches run_roks_periodic_multi_k_ewald3d at Gamma (a (1,1,1) mesh) or on a full Monkhorst-Pack mesh via kpoints=[n1, n2, n3]. Each per-spin Fock is \(F_\sigma(k) = H_\mathrm{core}(k) + J(k) - c_\mathrm{full} K_\sigma(k) + V_{xc}^\sigma(k)\), and the two are combined by the same Roothaan coupling ROHF uses. Supported: LDA, GGA, meta-GGA and global hybrids. A pure functional builds no exchange at all: neither the exchange caches nor the \(O(n_k^2)\) long-range \(K\) loop, which is what makes dense meshes affordable there. Screened hybrids (hse06 and friends) need an erfc-attenuated exchange arm this engine does not carry and fail closed, as do long-range-corrected functionals and double hybrids. The Li/PBE/STO-3G (1,1,2) anchor agrees with an out-of-process PySCF KROKS().density_fit() reference to 0.11 mHa, and the closed-shell PBE0 limit reproduces the BIPOLE RKS route to 0.02 mHa.

Note

Like the ROHF route, this driver keeps the cell’s alpha and beta counts at every k point. Its ROKS-specific occupation rule also detects when the one open electron selects an exactly degenerate orbital set. It locks that set as one restricted-open shell and shares the electron equally over its members; for example, a twofold shell has occupations 0.5/0.5. This is the fixed, zero-temperature ensemble for a symmetry-degenerate shell, not electronic smearing. The Li/PBE/STO-3G 6-bohr-box (2,2,1) case now converges in 6 iterations instead of 87-200 cycles with an arbitrary integer SOMO.

General metallic filling, finite-temperature open-shell smearing, and multiple partially occupied open shells are not inferred by this path. Use a supported unrestricted route when the system needs those occupation models.

ROHF on the native GDF route. run_periodic_job(method="ROHF", jk_method="gdf") is the normal output-producing entry point; run_krohf_periodic_gdf (periodic_rohf_gdf.py) is the density-fitted alternative to the Ewald-3D engine above. Gamma is the (1,1,1) mesh of the same code path, so there is only one driver to reason about:

from vibeqc import run_periodic_job

result = run_periodic_job(
    system, basis, method="ROHF", jk_method="gdf",
    kpoints=(3, 1, 1),                    # omit for Gamma
    aux_basis="def2-svp-jk", gdf_method="rsgdf",
)
result.energy, result.s_squared, result.mo_occupations[0]

It reuses the multi-k GDF machinery verbatim – the per-\((k_i, k_j)\) cderi cache, Hartree \(J\) from the BZ-summed total density, the exxdiv='ewald' BvK-supercell Madelung shift applied once per spin, and the converged Ewald nuclear repulsion – and differs from run_kuhf_periodic_gdf only in the orbital step, where the two spin Focks are combined by Roothaan’s effective Fock. At multiplicity 1 it reproduces run_krhf_periodic_gdf to machine precision; Li/STO-3G in a 10-bohr box at Gamma sits \(+2.86\times10^{-5}\) Ha from out-of-process PySCF 2.13.1 KROHF/GDF, and that residual collapses to \(+2.2\times10^{-7}\) Ha as the rsgdf mesh is refined – the shared fit-truncation ladder, not a convention difference.

Exchange is the only \(O(n_k^2)\) term, so it is built by contracting the cderi cache against the occupied MO blocks rather than the \(n_\mathrm{bf}\times n_\mathrm{bf}\) densities: the same operator at \(n_\mathrm{occ}/n_\mathrm{bf}\) of the cost, as two GEMM calls per k-pair (so the work stays inside threaded BLAS). Since 2026-08-02 every multi-k GDF driver builds exchange that way, not just this one. The rsgdf cderi setup rides the existing MPI k-point farming described in § 6.

Note

That contraction is an exact reassociation, not an approximation, and it does not change what any run reports: converged energies move by at most \(2.5\times10^{-11}\) Ha, with identical iteration counts. It also does not speed up a small job noticeably, because on test-scale cells the exchange build is a negligible fraction of the wall time (the cderi build dominates); the saving is aimed at production meshes and realistic bases, where exchange is what dominates.

ROKS on GDF is a separate increment and is not implemented: a functional raises rather than silently running something else. Smearing and 2D slabs raise as well. Analytic gradients, COSX exchange and symmetry-reduced (ibz_native) exchange are not part of this driver’s argument surface at all, so there is no way to ask for them. The public route adds .out, .system, Molden/QVF sidecars when applicable, and the ROHF plus GDF references; it does not change the standalone driver’s energy or gauge.

For both methods, GAPW, 1D/2D cells, multi-k GPW, gradients, restart densities, and electronic smearing fail closed. The standalone drivers’ smearing_alpha argument controls nuclear-charge smoothing in the one-electron potential, not electronic Fermi smearing.

ROHF avoids the spin contamination of UHF at the cost of the more involved coupling operator; reach for it when a clean \(\langle S^2\rangle\) matters more than variational freedom.

4. GPW and GAPW: the plane-wave / augmented route

4.1 Theory

The Gaussian + plane-wave (GPW, Lippert, Gerald, Hutter, Jürg, and Parrinello, Michele (1997). *A hybrid Gaussian and plane wave density functional scheme.* Molecular Physics, **92**(3), 477-487. doi:[10.1080/00268979709482119](https://doi.org/10.1080/00268979709482119)) and Gaussian-augmented plane-wave (GAPW) methods, due to Lippert and Hutter and implemented as the core CP2K periodic-DFT engine (VandeVondele, Joost, Krack, Matthias, et al. (2005). *Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach.* Computer Physics Communications, **167**(2), 103-128. doi:[10.1016/j.cpc.2004.12.014](https://doi.org/10.1016/j.cpc.2004.12.014)), replace the four-centre Coulomb build with a Poisson solve on a uniform real-space grid:

1. Collocate the Gaussian-built electron density onto a uniform 3D grid sampled in fractional cell coordinates, \(\rho(\mathbf{r}) = \sum_{\mu\nu} D_{\mu\nu} \chi_\mu(\mathbf{r})\chi_\nu(\mathbf{r})\).

2. Solve Poisson in reciprocal space: \(V_H(\mathbf{G}) = 4\pi\,\rho(\mathbf{G})/|\mathbf{G}|^2\), with \(V_H(\mathbf{G}=0)\) pinned to zero. A neutral cell handles the long-range tail through this gauge; an unbalanced charge gets the standard uniform-background self-image shift, which is the same Madelung convention as the rest of vibe-qc.

3. Project back to the AO basis: \(J_{\mu\nu} = \int \chi_\mu(\mathbf{r}) \chi_\nu(\mathbf{r}) V_H(\mathbf{r})\,\mathrm{d}\mathbf{r} \approx \Delta V \sum_g \chi_\mu(\mathbf{r}_g) \chi_\nu(\mathbf{r}_g) V_H(\mathbf{r}_g)\).

The cost of the FFT is \(\mathcal{O}(N_g \log N_g)\) where \(N_g = N_x N_y N_z\) is set by the plane-wave cutoff \(|\mathbf{G}_{\max}|^2 = 2 E_{\text{cut}}\). For a smooth pseudo- density that is much cheaper than building the full four-centre ERI tensor.

The catch is that core densities are sharply peaked, and a grid fine enough to resolve them is enormous. The two responses to this are:

• GPW: use a pseudopotential to remove the core electrons. The remaining valence pseudo-density is smooth and a moderate \(E_{\text{cut}}\) (tens to a few hundred Ry) suffices.

• GAPW (Lippert and Hutter, 1999, Lippert, Gerald, Hutter, Jürg, and Parrinello, Michele (1999). *The Gaussian and augmented-plane-wave density functional method for ab initio molecular dynamics simulations.* Theoretical Chemistry Accounts, **103**(2), 124-140. doi:[10.1007/s002140050523](https://doi.org/10.1007/s002140050523); all-electron energy expression: Krack, Matthias and Parrinello, Michele (2000). *All-electron ab-initio molecular dynamics.* Physical Chemistry Chemical Physics, **2**(10), 2105-2112. doi:[10.1039/b001167n](https://doi.org/10.1039/b001167n)): keep the all-electron problem, but split the density into a smooth part \(\tilde\rho\) that lives on the coarse plane-wave grid plus per-atom hard / soft corrections \(\rho_a - \tilde\rho_a\) that live on a fine radial grid around each nucleus. The Hartree potential decomposes as

  \[ V_H[\rho] = V_H[\tilde\rho] + \sum_a \Big( V_H[\rho_a] - V_H[\tilde\rho_a] \Big), \]

  where the smooth piece reuses the GPW Poisson solve and the atomic pieces are local 1D radial Poisson solves with no FFT. The augmentation restores all-electron accuracy without paying the cost of resolving the core on the smooth grid. The PAW generalisation (Blochl 1994, Blöchl, P. E. (1994). *Projector augmented-wave method.* Physical Review B, **50**(24), 17953–17979. doi:[10.1103/PhysRevB.50.17953](https://doi.org/10.1103/PhysRevB.50.17953)) takes the same dual-grid idea further by transforming all-electron orbitals through a per-atom projector basis.

4.2 Status in vibe-qc

The GPW / GAPW track grew out of the v0.10.x acceleration program and now drives a full periodic SCF. The design doc with the full milestone history is at design_periodic_gapw.md.

Both routes are reached with jk_method="gpw" or jk_method="gapw" on run_periodic_job. GPW is production within its shipped envelope: Gamma RHF / UHF / RKS / UKS, multi-k pure-DFT RKS / ROKS / UKS, Gamma analytic gradients, DFT+U (Dudarev, S. L., Botton, G. A., et al. (1998). *Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study.* Physical Review B, **57**(3), 1505–1509. doi:[10.1103/PhysRevB.57.1505](https://doi.org/10.1103/PhysRevB.57.1505)) and D3-BJ (Grimme, Stefan, Ehrlich, Stephan, and Goerigk, Lars (2011). *Effect of the damping function in dispersion corrected density functional theory.* Journal of Computational Chemistry, **32**(7), 1456–1465. doi:[10.1002/jcc.21759](https://doi.org/10.1002/jcc.21759)). The .out route description ends with an experimental tag (J/K method: GPW -- ...; experimental) while the GPW / GAPW family is completed. run_periodic_job calls the GPW drivers with quiet=True, so the drivers’ own opt-in warning is silent, but the grid and J-builder infrastructure still emit GAPWExperimentalWarning once per build (measured 2026-09-29 on a Gamma RHF run: PlaneWaveGrid built at cutoff=... and GPW J build: GpwJBuilder ...); read that warning as a scope reminder, not as a defect flag. GAPW is experimental, and GPW ROHF / ROKS are a maintained preview. GAPW RHF/UHF additionally requires gapw_molecular_limit=True; without that explicit isolated-cell declaration it fails closed because its promoted analytic one-centre correction is not a dense-crystal image-summed functional. What runs today:

• Gamma-point SCF for RHF, RKS, UHF, and UKS via run_periodic_rhf_gpw / run_periodic_uhf_gpw / run_periodic_uks_gpw (GPW) and run_periodic_rhf_gapw / run_periodic_uhf_gapw / run_periodic_uks_gapw (GAPW). RKS / UKS require a functional=.

• Maintained-preview 3D Gamma-point ROHF and ROKS via run_periodic_rohf_gpw / run_periodic_roks_gpw, or the corresponding run_periodic_job(method="ROHF", jk_method="gpw") or run_periodic_job(method="ROKS", functional=..., jk_method="gpw") route. Both use GPW J and the shared Roothaan loop with integer 2/1/0 occupations; ROHF adds per-spin exact K and ROKS adds spin-polarised XC plus functional-scaled K. They are not GAPW, multi-k, or electronic-smearing routes.

• Multi-k RKS (pure DFT) via run_periodic_rks_gpw_multi_k and run_periodic_rks_gapw_multi_k, including LDA/GGA and self-regularising meta-GGAs such as r2SCAN (Furness, James W., Kaplan, Aaron D., et al. (2020). *Accurate and Numerically Efficient r²SCAN Meta-Generalized Gradient Approximation.* Journal of Physical Chemistry Letters, **11**(19), 8208–8215. doi:[10.1021/acs.jpclett.0c02405](https://doi.org/10.1021/acs.jpclett.0c02405)). Hybrid / HF multi-k is not yet wired on this route and raises.

• Forces, +U, Fermi-Dirac smearing (Mermin, N. D. (1965). *Thermal Properties of the Inhomogeneous Electron Gas.* Physical Review, **137**(5A), A1441–A1443. doi:[10.1103/PhysRev.137.A1441](https://doi.org/10.1103/PhysRev.137.A1441)), and D3(BJ) dispersion. GPW/GAPW cover RHF/RKS/UHF/UKS +U at Gamma and RKS +U at multi-k; multi-k open-shell +U is not wired. GPW retains the high-level finite-difference optimization and Hessian plumbing. GAPW RHF/UHF instead rejects high-level optimize, optimize_cell, and hessian requests so a follow-on workflow cannot silently change the analytic one-centre energy surface. Fixed-cell GAPW forces remain available through VibeqcGAPW; block-mode DFT retains the existing derivative plumbing.

The milestone trail that landed this:

Milestone	Status	Deliverable
M1	landed	PlaneWaveGrid infrastructure + Gaussian-on-grid collocation
M2	landed	GPW Hartree-J via FFT-Poisson + full Gamma-point SCF
M3	landed	GAPW dual-grid augmentation; multi-k RKS; open-shell UHF / UKS
M4	in progress	Metallic GAPW refinements with smearing
M5 (stretch)	scheduled	ACE-decorated K for hybrids + on-demand PAW dataset fetcher

The grid kernels live in periodic_gapw_j.py (GPW) and periodic_gapw_augment.py (GAPW augmentation), with the runner dispatch in periodic_runner.py. The GpwJBuilder orchestrates density collocation -> FFT-Poisson -> AO-basis projection into the Hartree-J matrix; the GAPW driver adds the per-atom augmentation on top.

Note

GPW is production within the envelope above and GAPW is experimental. The experimental tag on the GPW route line, and the infrastructure warning, are there because the route is younger than GDF / BIPOLE and its parity sweep against CP2K is still broadening, not because it is unimplemented. For tight-cell hybrid work the mature choice is still GDF or BIPOLE. The label vocabulary is defined in Methods: theory and status.

4.3 Advantages

• Linear-log scaling Hartree-J. FFT-Poisson is \(\mathcal{O}(N_g \log N_g)\), independent of basis size. For large supercells with modest valence densities this is the decisive scaling win.

• Natural fit for large cells. Slabs with hundreds of atoms and metallic systems where smearing matters benefit most.

• GAPW gives all-electron accuracy without bundled pseudopotentials, by carrying the core augmentation on per-atom radial grids.

• Pluggable exchange. GAPW replaces the J build only; exchange goes through whichever K-builder the existing periodic SCF would use. ACE decoration (Lin Lin 2016, Lin, Lin (2016). *Adaptively Compressed Exchange Operator.* Journal of Chemical Theory and Computation, **12**(5), 2242–2249. doi:[10.1021/acs.jctc.6b00092](https://doi.org/10.1021/acs.jctc.6b00092)) is the planned hybrid-DFT accelerator and is being designed so it can be shared with the locality-exploitation track (ADMM, Guidon, Manuel, Hutter, Jürg, and VandeVondele, Joost (2010). *Auxiliary Density Matrix Methods for Hartree-Fock Exchange Calculations.* Journal of Chemical Theory and Computation, **6**(8), 2348–2364. doi:[10.1021/ct1002225](https://doi.org/10.1021/ct1002225); OT, VandeVondele, Joost and Hutter, Jürg (2003). *An efficient orbital transformation method for electronic structure calculations.* Journal of Chemical Physics, **118**(10), 4365–4369. doi:[10.1063/1.1543154](https://doi.org/10.1063/1.1543154)).

4.4 Disadvantages (intrinsic and current)

• Younger than GDF / BIPOLE. GPW runs a full SCF and is production within its envelope, but its CP2K parity sweep is still broadening, which is why the .out route line keeps an experimental tag; GAPW is experimental. For a production tight-cell hybrid, prefer GDF or BIPOLE today.

• Multi-k is pure-DFT RKS only. Hybrid and HF multi-k are not yet wired on the GPW / GAPW route; multi-k hybrids stay on BIPOLE.

• Plane-wave cutoff is an extra knob. Sensible defaults are derived from the tightest primitive exponent and a CP2K-style REL_CUTOFF safety factor (calibrated against REL_CUTOFF = 60 Ry).

• Pseudopotential licensing is non-trivial. The v0.10.x design decision (see design_periodic_gapw.md section 9) is to bundle no PAW dataset and instead build a vqfetch-style on-demand fetcher that filters out GPL / proprietary records. GBRV (CC-BY-SA), ABINIT-JTH (mixed), GPAW datasets (GPL), PSlibrary (mixed), and VASP POTCARs (proprietary) all fail the bundle test under MPL-2.0 (see docs/license.md for the full inventory).

• 3D-focused. The route is exercised mainly on 3D periodic cells; 2D slab and 1D wire coverage is thinner than on GDF, which remains the go-to for low-dimensional systems.

• FFTW3 dependency is GPL-v2. The combined binary linking FFTW3 (Frigo, Matteo and Johnson, Steven G. (2005). *The Design and Implementation of FFTW3.* Proceedings of the IEEE, **93**(2), 216–231. doi:[10.1109/JPROC.2004.840301](https://doi.org/10.1109/JPROC.2004.840301)) is effectively GPL-v2. A future FFT-backend abstraction (VIBEQC_FFT_BACKEND={fftw3, pocketfft, kissfft}) is on the roadmap to allow commercial-friendly binaries.

4.5 When to use GPW / GAPW

Reach for GPW / GAPW when you want plane-wave Hartree-J on a large smooth-density cell: large-cell periodic DFT (slabs and bulk supercells with hundreds of atoms) where the BIPOLE direct-K cost and the GDF aux-basis cost both become uncomfortable, and where the \(\mathcal{O}(N_g \log N_g)\) FFT-Poisson build pays off as \(N\) grows. GAPW adds all-electron accuracy without pseudopotential files. Both run a full \(\Gamma\)-point SCF (RHF / RKS / UHF / UKS) plus multi-k RKS today; GPW is production within that envelope, GAPW is experimental. For tight-cell hybrids and for production CRYSTAL parity, stay on GDF / BIPOLE until the GPW / GAPW parity sweep and hybrid multi-k land.

4.6 Citations

• G. Lippert, J. Hutter, M. Parrinello, “A hybrid Gaussian and plane-wave density functional scheme”, Mol. Phys. 92, 477 (1997).

• G. Lippert, J. Hutter, M. Parrinello, “The Gaussian and augmented-plane-wave density functional method for ab initio molecular dynamics simulations”, Theor. Chem. Acc. 103, 124 (1999). The GAPW formulation.

• M. Krack, M. Parrinello, “All-electron ab initio molecular dynamics”, Phys. Chem. Chem. Phys. 2, 2105 (2000).

• J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing, J. Hutter, “Quickstep: fast and accurate density functional calculations using a mixed Gaussian and plane waves approach”, Comput. Phys. Commun. 167, 103 (2005). The CP2K Quickstep reference.

• P. E. Blochl, “Projector augmented-wave method”, Phys. Rev. B 50, 17953 (1994). The PAW formulation that the M5 stretch goal targets.

• L. Lin, “Adaptively compressed exchange operator”, J. Chem. Theory Comput. 12, 2242 (2016). The ACE acceleration of exact exchange that M5 wraps around GAPW for hybrid-DFT scaling.

4.7 RSGAPW: range-separated GAPW for screened hybrids

Range-separated hybrids (HSE06, ωB97X, ωB97M (Mardirossian, Narbe and Head-Gordon, Martin (2016). *ωB97M-V: A combinatorially optimized, range-separated hybrid, meta-GGA density functional with VV10 nonlocal correlation.* Journal of Chemical Physics, **144**(21), 214110. doi:[10.1063/1.4952647](https://doi.org/10.1063/1.4952647)), CAM-B3LYP) split the Coulomb operator with the error function, \(1/r = \mathrm{erf}(\omega r)/r + \mathrm{erfc}(\omega r)/r\), and apply a different fraction of exact exchange to each piece. vibe-qc carries these on the plane-wave grid through a dedicated range-separated GAPW path:

from vibeqc.periodic_gapw_range_sep import (
    run_periodic_rhf_rsgapw,   # range-separated HF exchange
    run_periodic_rks_rsgapw,   # range-separated hybrid KS
)

The long-range (smooth) Coulomb piece is solved by FFT-Poisson on the plane-wave grid; the short-range piece uses libint erfc-screened integrals. The GAPW per-atom augmentation sees only the total (unscreened) Coulomb potential, so the \(\omega\)-split applies to the smooth-grid J and the K builders alone. omega_for_functional(name) returns the \((\omega, \text{HF}_{SR}, \text{HF}_{LR})\) triple for the supported functionals.

This is a dedicated entry point, not a jk_method route: you call run_periodic_rhf_rsgapw / run_periodic_rks_rsgapw directly rather than selecting it through run_periodic_job(jk_method=...). Citations: Heyd, Scuseria, Ernzerhof, J. Chem. Phys. 118, 8207 (2003) (HSE, Heyd, Jochen, Scuseria, Gustavo E., and Ernzerhof, Matthias (2003). *Hybrid functionals based on a screened Coulomb potential.* Journal of Chemical Physics, **118**(18), 8207–8215. doi:[10.1063/1.1564060](https://doi.org/10.1063/1.1564060); HSE06 screening parameter: Krukau, Aliaksandr V., Vydrov, Oleg A., et al. (2006). *Influence of the exchange screening parameter on the performance of screened hybrid functionals.* Journal of Chemical Physics, **125**(22), 224106. doi:[10.1063/1.2404663](https://doi.org/10.1063/1.2404663)); Chai and Head-Gordon, J. Chem. Phys. 128, 084106 (2008) (ωB97X, Chai, Jeng-Da and Head-Gordon, Martin (2008). *Systematic optimization of long-range corrected hybrid density functionals.* Journal of Chemical Physics, **128**(8), 084106. doi:[10.1063/1.2834918](https://doi.org/10.1063/1.2834918)); Yanai, Tew, Handy, Chem. Phys. Lett. 393, 51 (2004) (CAM-B3LYP, Yanai, Takeshi, Tew, David P., and Handy, Nicholas C. (2004). *A new hybrid exchange-correlation functional using the Coulomb-attenuating method (CAM-B3LYP).* Chemical Physics Letters, **393**(1–3), 51–57. doi:[10.1016/j.cplett.2004.06.011](https://doi.org/10.1016/j.cplett.2004.06.011)).

4.8 RIJCOSX: density-fitted J with chain-of-spheres K

RIJCOSX pairs a GDF-fitted Coulomb build with a seminumerical chain-of-spheres (COSX) exchange build on a periodic Becke grid (Becke, A. D. (1988). *A multicenter numerical integration scheme for polyatomic molecules.* Journal of Chemical Physics, **88**(4), 2547–2553. doi:[10.1063/1.454033](https://doi.org/10.1063/1.454033)). It is the periodic analogue of the molecular RIJCOSX accelerator (Neese 2009, Neese, Frank, Wennmohs, Frank, et al. (2009). *Efficient, approximate and parallel Hartree-Fock and hybrid DFT calculations. A ‘chain-of-spheres’ algorithm for the Hartree-Fock exchange.* Chemical Physics, **356**(1–3), 98–109. doi:[10.1016/j.chemphys.2008.10.036](https://doi.org/10.1016/j.chemphys.2008.10.036)) and is selected with jk_method="rijcosx":

result = vq.run_periodic_job(sysp, basis, method="RHF",
                             jk_method="rijcosx")

What runs today:

• Gamma RHF. The dedicated run_periodic_rijcosx_rhf driver keeps the original periodic RIJCOSX path for one-point closed-shell HF. On a tight 3-D cell its gdf_method="rsgdf" fit is the classic fit on physical AO pairs; classic_support_policy="historical_index_ball" replays the former origin-centred ball (GitLab #466). Its two Gamma JK builders iterate on the smooth seminumerical exchange surface, then apply the analytic one-center replacement to the final Fock and orbitals after convergence. The reported energy remains on the iterated seminumerical surface, matching the molecular RIJCOSX contract.

• True multi-k RHF / RKS / UHF / UKS. run_periodic_job(..., jk_method="rijcosx", kpoints=(...)) uses the native multi-k GDF loop for RI-J and switches exact exchange to the composed periodic COSX backend (k_exchange="cosx"). Closed-shell density mixers, closed-shell READ restart, closed-shell +U, and closed-shell rsgdf_tail_ke_cutoff follow the same support envelope as the GDF multi-k route. After convergence the multi-k drivers apply the analytic one-center replacement to the returned Fock matrices and orbitals, evaluated with the same short-range kernel as the bridge’s real-space exchange; the reported energy remains on the iterated seminumerical surface, matching the molecular RIJCOSX and Gamma-builder contract.

• One-cell RKS on vacuum-padded cells. run_periodic_job(..., method="RKS", jk_method="rijcosx", kpoints=(1, 1, 1)) runs the Gamma RIJCOSX RKS driver when every active lattice length is at least 12 bohr and the cell is not tight: molecular DF J plus periodic COSX K through the dedicated Gamma builder, XC integrated on the molecular Becke grid (exact in that envelope), and the post-convergence one-center replacement applied to the returned Kohn-Sham matrix. Hybrids route their exact-exchange fraction through COSX; pure functionals never build K.

• One-cell UHF/UKS on vacuum-padded cells. The open-shell Gamma RIJCOSX drivers follow the same envelope rule as RKS: molecular DF J from the total density, per-spin periodic COSX K, per-spin post-convergence one-center replacement, and (UKS) XC on the molecular Becke grid. Spin counts follow the system multiplicity.

• Home-cell electrostatics on the vacuum-padded envelope. Because the molecular DF J of these one-cell drivers contains no image-cell electrons, the nuclear attraction, the nuclear repulsion and the COSX exchange are evaluated over the home cell as well, regardless of cutoff_bohr and nuclear_cutoff_bohr. In periodic boundary conditions only the total electrostatic energy of a neutral cell is well defined (Makov and Payne 1995, Makov, G. and Payne, M. C. (1995). *Periodic boundary conditions in ab initio calculations.* Physical Review B, **51**(7), 4014–4022. doi:[10.1103/PhysRevB.51.4014](https://doi.org/10.1103/PhysRevB.51.4014)); keeping image nuclei or the image-ket exchange class without the matching image electrons leaves the bare monopole sums in the energy. Before this was enforced (issue #707), a one-electron H atom in a 12-bohr box converged to -1.638 Ha instead of -0.467 Ha, and H2 RHF in a 24-bohr box was 0.5 Ha too low. Every box inside the envelope now reproduces the molecular RIJCOSX result to the COSX quadrature floor; the neglected image terms are the neutral-cell multipole tail that the molecular DF J already neglects.

• One-cell tight-cell RKS/UHF/UKS fail closed. Tight or short-axis cells need the periodic-density XC grid and the Gamma-folded exchange classes, so vibe-qc refuses to label those one-cell calculations as RIJCOSX.

When to use it: explore the periodic COSX approximation for HF and hybrid DFT on compact periodic cells, with GDF as the independent reference route. For pure functionals such as PBE, RIJCOSX and RI share the same zero-exact- exchange KS Fock build; use a hybrid such as PBE0 when you specifically want to exercise COSX exchange. Citation: F. Neese, F. Wennmohs, A. Hansen, U. Becker, “Efficient, approximate and parallel Hartree-Fock and hybrid DFT calculations. A ‘chain-of-spheres’ algorithm for the Hartree-Fock exchange”, Chem. Phys. 356, 98 (2009).

5. Side-by-side comparison

5.1 Dimensionality and capability matrix

	1D wires	2D slabs	3D bulk	RHF	UHF	RKS pure	RKS hybrid	UKS	Multi-k	Analytic gradient	Cell opt
SLAB_EWALD_2D	no	yes	no	yes	no	yes	yes	yes	yes	no	no
BIPOLE	no (direct-API diagnostic only)	raises	yes	yes	yes	yes	yes	yes	yes	preview	no
GDF (native)	no (raises; use dim=3 with vacuum)	closed-shell	yes	yes	yes (Γ + multi-k)	yes	yes	yes (Γ + multi-k)	yes	yes (Γ / multi-k / closed-shell slab)	no
GPW	limited	raises	yes	yes (Γ)	yes (Γ)	yes (Γ)	yes (Γ)	yes (Γ)	yes (RKS)	FD	no
GAPW	limited	raises	yes	yes (Γ)	yes (Γ)	yes (Γ)	yes (Γ)	yes (Γ)	yes (RKS)	block only	no
RIJCOSX	yes (multi-k)	raises	yes	yes	yes (multi-k)	yes (multi-k)	yes (multi-k)	yes (multi-k)	yes	no	no
RSGAPW	limited	raises	yes	yes (Γ)	no	yes (Γ)	yes (Γ, RS)	no	no	FD	FD
ROHF / ROKS	no	no	yes	ROHF (Ewald; GDF Γ + multi-k; GPW Γ)	n/a	ROKS (Ewald; GPW Γ)	ROKS (Ewald global hybrids; GPW Γ)	no	yes (ROHF + ROKS Ewald)	no	no

Slabs have two routes. SLAB_EWALD_2D is what jk_method="auto" selects for dim=2, and the explicit slab-GDF route is the closed-shell alternative that also carries analytic gradients (optimize=True, 2026-07-30). Every bulk builder raises on dim=2 rather than returning an a3-dependent energy. Fermi-Dirac smearing and open-shell HF still need a dim=3 cell (vq.slab(..., periodic_z=True)); so do freeze-mask relaxation and NEB, which run on the BIPOLE force stack.

The direct run_pbc_bipole_* APIs retain a 1D DIRECT_TRUNCATED molecular-limit diagnostic. run_periodic_job does not expose it: explicit BIPOLE rejects dim=1, and so does AUTO. Run a wire as a dim=3 cell with transverse vacuum.

Notes on the matrix: GDF open-shell (UHF / UKS) ships at Gamma and multi-k in 3D only; a dim=1 wire has no open-shell route. In 3D, AUTO routes open shell to BIPOLE by design, so you opt into GDF with an explicit jk_method="gdf"; in 1D, AUTO refuses. GPW / GAPW multi-k is pure-DFT RKS only (hybrid / HF multi-k raise on that route). GPW forces and fixed-cell gradients use the finite-difference path the runner shares; cell optimization fails closed. GAPW direct analytic gradients cover block mode, while high-level HF optimization and Hessian requests fail closed; block-mode DFT retains its existing derivative plumbing, but cell optimization also fails closed. VibeqcGAPW uses analytic block-mode forces or full-SCF finite differences for the molecular-limit analytic mode. RIJCOSX is Gamma RHF plus true multi-k RHF/RKS/UHF/UKS through the GDF/COSX backend; one-cell RKS/UHF/UKS requests fail closed so they are not mislabeled. RSGAPW (range-separated GAPW) remains a dedicated entry point. ROHF has dedicated Ewald entry points, public and standalone native-GDF entry points (run_periodic_job(method="ROHF", jk_method="gdf") and run_krohf_periodic_gdf, Gamma + multi-k), and a public jk_method="gpw" route at 3D Gamma (the ROHF “RHF” column reads as restricted open-shell HF on the Ewald-3D gauge at Gamma and multi-k, on the GDF route at Gamma and multi-k, plus GPW at Gamma). RSGAPW handles screened hybrids such as HSE06 / ωB97X / CAM-B3LYP.

5.2 Scaling and cost intuition

	Per-iter J	Per-iter K	Memory	Notes
BIPOLE	\(\mathcal{O}(N^2)\) per cell + reciprocal \(J^{LR}\)	\(\mathcal{O}(N^4)\) direct with Schwarz	\(\mathcal{O}(N^2)\)	\(K\) dominates at large cutoff. The quartet multipole prototype is fail-closed and is not a supported shortcut.
GDF	\(\mathcal{O}(N^2 M)\) via \(L_{P,\mu\nu}\) on fitted branches; \(\mathcal{O}(N_k N^2 N_g N_c)\) batched AFT on the 3D pure-DFT Ewald-J fallback	\(\mathcal{O}(N^2 M)\) via \(L_{P,\mu\nu}\); absent for pure DFT	\(\mathcal{O}(N^2 M)\) for fitted \(L\); bounded FT batch plus \(\mathcal{O}(N_c N^2)\) J blocks on the pure-DFT fallback	Aux size \(M \approx 3N\); the pure-DFT fallback never retains the all-cell/all-\(G\) FT tensor.
GPW / GAPW	\(\mathcal{O}(N_g \log N_g)\) FFT + \(\mathcal{O}(N N_g)\) collocation + projection	inherits whichever K builder is composed with it (target: ACE)	\(\mathcal{O}(N_g)\) for the grid	\(N_g\) depends on the plane-wave cutoff, not on \(N\); the win shows up as \(N\) grows.

N is the AO basis size per unit cell, M is the auxiliary basis size, N_k is the number of k-points, N_g is the reciprocal or plane-wave grid count, and N_c is the number of real-space lattice cells.

5.3 Method-selection chart

Periodic SCF
  |
  +-- dim=2 slab
  |     +-- AUTO: SLAB_EWALD_2D (energies)
  |     +-- explicit closed-shell GDF (energies + forces)
  |
  +-- dim=1 wire
  |     +-- AUTO raises for every method; use a dim=3 cell with transverse vacuum
  |
  +-- dim=3
        +-- UHF/UKS: AUTO selects BIPOLE
        +-- closed shell: GDF by default
              +-- BIPOLE for exact-route FD forces and fixed-cell optimization
              +-- GPW/GAPW for large smooth-density systems

The slab branch comes first because it is exclusive: no bulk builder accepts dim=2. Explicit closed-shell slab GDF supplies analytic forces. For smearing, open shell, freeze masks, or periodic Gaussian NEB, build the cell as dim=3 with vq.slab(..., periodic_z=True) instead. These are separate capability envelopes: finite-temperature BIPOLE single points are supported, but BIPOLE atom optimization and periodic Gaussian NEB require zero-temperature occupations.

The open-shell branch points at BIPOLE because that is the AUTO default; open-shell GDF (UHF / UKS, Gamma + multi-k) also ships and is reachable with an explicit jk_method="gdf".

5.4 Coexistence

The three routes share one infrastructure stack and may be mixed in the same session without re-installation:

• S, T, V_{ne} lattice integrals live in C++ behind compute_overlap_lattice, compute_kinetic_lattice, and the shared Ewald-gauged \(V_{ne}\) helper. All three drivers consume them.

• nuclear_repulsion_per_cell is the same Ewald-3D nuclear sum for everyone on a 3D system.

• BasisSet and Molecule constructors are shared with the molecular SCF; the same auxiliary-basis manager (vibeqc.aux_basis.make_aux_basis_set) backs both molecular RIJK and periodic GDF.

• The .system manifest records which J/K route was active, plus its parameters (Ewald omega for BIPOLE, aux basis name + linear_dep_thr for GDF, plane-wave cutoff and grid shape for GPW / GAPW). A future re-run reproduces the exact route.

5.5 Stress and variable-cell relaxation

The stress tensor is the strain derivative of the energy per volume (Nielsen, O. H. and Martin, Richard M. (1985). *Quantum-mechanical theory of stress and force.* Physical Review B, **32**(6), 3780–3791. doi:[10.1103/PhysRevB.32.3780](https://doi.org/10.1103/PhysRevB.32.3780)). With the symmetric strain \(\boldsymbol\epsilon\) applied to lattice vectors and atomic positions alike, \(\mathbf a_i \to (\mathbf 1 + \boldsymbol\epsilon)\mathbf a_i\) and \(\mathbf R_A \to (\mathbf 1 + \boldsymbol\epsilon)\mathbf R_A\),

\[ \sigma_{ij} = \frac{1}{V}\,\frac{\partial E}{\partial \epsilon_{ij}} \approx \frac{1}{V}\,\frac{E(+h\,\mathbf e_{ij}) - E(-h\,\mathbf e_{ij})}{2h}, \]

and vibe-qc has no analytic periodic stress on any ab initio route: every stress it reports is this central difference of the total energy, with the basis rebuilt on each strained cell so that Hellmann-Feynman (Feynman, R. P. (1939). *Forces in Molecules.* Physical Review, **56**(4), 340–343. doi:[10.1103/PhysRev.56.340](https://doi.org/10.1103/PhysRev.56.340)) and Pulay (Pulay, P. (1969). *Ab initio calculation of force constants and equilibrium geometries in polyatomic molecules. I. Theory.* Molecular Physics, **17**(2), 197–204. doi:[10.1080/00268976900100941](https://doi.org/10.1080/00268976900100941)) contributions arrive together.

• GPW / GAPW: compute_stress_gpw(system, basis, result, *, basis_name, ..., eps=1e-4) and compute_stress_gapw(...) in vibeqc.periodic_gapw_stress return the \(3\times 3\) tensor in Ha/bohr³ from finite differences of the converged GPW or GAPW energy.

• Semiempirical: compute_stress_fd(system, energy_fn, h=0.001) in vibeqc.semiempirical.periodic differentiates any energy callable; for dim < 3 only the first dim rows and columns are meaningful.

• optimize_cell=True on run_periodic_job (it requires optimize=True) is honoured only on the full-k DFTB0 (Porezag, D., Frauenheim, Th., et al. (1995). *Construction of tight-binding-like potentials on the basis of density-functional theory: Application to carbon.* Physical Review B, **51**, 12947. doi:[10.1103/PhysRevB.51.12947](https://doi.org/10.1103/PhysRevB.51.12947)) and SCC-DFTB (Elstner, M., Porezag, D., et al. (1998). *Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties.* Physical Review B, **58**, 7260. doi:[10.1103/PhysRevB.58.7260](https://doi.org/10.1103/PhysRevB.58.7260)) semiempirical routes, where vibeqc.semiempirical.periodic.optimize_cell relaxes positions and lattice together through ASE (Larsen, Ask Hjorth, Mortensen, Jens Jørgen, et al. (2017). *The atomic simulation environment - a Python library for working with atoms.* Journal of Physics: Condensed Matter, **29**(27), 273002. doi:[10.1088/1361-648X/aa680e](https://doi.org/10.1088/1361-648X/aa680e); isotropic=True restricts the cell to a volume scaling). On BIPOLE it raises (BIPOLE variable-cell optimization is unavailable. The historical strain convention and coupled atom/cell convergence were not certified), on GDF (optimize_cell=True is not implemented on the GDF route -- there is no GDF analytic stress), on the slab-GDF route (no slab GDF analytic stress), on AICCM (variable-cell relaxation is not implemented) and on external-XC routes. Use a converged equation-of-state scan (vibeqc.eos) for equilibrium volumes and bulk moduli on those routes.

• Not a stress: compute_stress_tensor(system, gradient) in vibeqc.bipole_gradient is the atomic-force virial \(s_{ij} = -(1/V)\sum_A R_{A,i} F_{A,j}\). It omits the explicit lattice and Ewald strain dependence and the Gaussian-basis Pulay stress, and on H₂/STO-3G it comes out opposite in sign to \(\partial E/\partial\epsilon\). It is kept only as a force-virial diagnostic; do not relax a cell with it or report it as a stress. The BIPOLE cell relaxers relax_cell, relax_cell_gradient and relax_full fail closed for the same reason.

See Methods: phonons, stress and variable cells and the periodic geometry optimisation tutorial.

6. Worked examples by bonding situation and dimensionality

The examples below intentionally cover different bonding regimes because the right route is rarely a function of dimensionality alone. Diffuseness of the valence density, charge transfer, band-gap, and the strength of the long-range Coulomb tail all push you toward different J/K kernels.

All systems use the canonical lattice parameters from examples/periodic/_systems.py (Springer Materials).

6.1 Ionic, closed-shell 3D: MgO and NaCl (rocksalt)

Ionic rocksalt crystals are the textbook stress test for the periodic Coulomb gauge: charge transfer of order one electron between cation and anion, Madelung sums that diverge under naive truncation, and a wide gap so the SCF itself is well-behaved. Both BIPOLE (Saunders, V. R., Freyria-Fava, C., et al. (1992). *On the electrostatic potential in crystalline systems where the charge density is expanded in Gaussian functions.* Molecular Physics, **77**, 629–665. doi:[10.1080/00268979200102671](https://doi.org/10.1080/00268979200102671)) and GDF (Sun, Qiming, Berkelbach, Timothy C., et al. (2017). *Gaussian and plane-wave mixed density fitting for periodic systems.* Journal of Chemical Physics, **147**(16), 164119. doi:[10.1063/1.4998644](https://doi.org/10.1063/1.4998644)) are production routes here.

import numpy as np
import vibeqc as vq

ANG2BOHR = 1.0 / 0.529177210903

# MgO FCC primitive (1 Mg + 1 O), a = 4.211 A.
a = 4.211 * ANG2BOHR
lattice = (a / 2.0) * np.array(
    [[0.0, 1.0, 1.0],
     [1.0, 0.0, 1.0],
     [1.0, 1.0, 0.0]]
)
mgo = vq.PeriodicSystem(
    3, lattice,
    [vq.Atom(12, [0.0, 0.0, 0.0]),
     vq.Atom(8,  [a / 2.0, a / 2.0, a / 2.0])],
)
basis = vq.BasisSet(mgo.unit_cell_molecule(), "sto-3g")

# Route A: BIPOLE, multi-k 2x2x2 (matches CRYSTAL14 SHRINK 2 baseline).
res_bipole = vq.run_periodic_job(
    mgo, basis, method="RHF",
    jk_method="bipole",
    kpoints=(2, 2, 2),
    output="mgo-bipole",
)

# Route B: GDF Gamma-only with the universal JK aux.
res_gdf = vq.run_periodic_job(
    mgo, basis, method="RHF",
    jk_method="gdf",
    aux_basis="def2-universal-jkfit",
    output="mgo-gdf",
)

print(f"BIPOLE 2x2x2 : {res_bipole.energy:+.6f} Ha")
print(f"GDF    Gamma : {res_gdf.energy:+.6f} Ha")

NaCl (a = 5.64 A) has a much more diffuse valence density on Cl and so probes aux conditioning more aggressively. Swap the geometry helpers and rerun; the GDF route may need drop_eta on the auxiliary basis if the metric is near-singular.

Pick BIPOLE for fixed-cell finite-difference forces, or GDF where its analytic gradient envelope applies. Use neither for variable-cell work.

6.2 Covalent network 3D: diamond and silicon (Fd-3m)

Diamond C and silicon are the canonical covalent network crystals: every atom sp3, narrow density between bonded atoms, no charge transfer, indirect band gap of order eV. The lattice vectors are FCC primitive (60-deg angles between the column vectors, not orthogonal axes).

import numpy as np
import vibeqc as vq

ANG2BOHR = 1.0 / 0.529177210903

# Si diamond primitive (2 Si in the rhombohedral cell), a = 5.431 A.
a = 5.431 * ANG2BOHR
lattice = (a / 2.0) * np.array(
    [[0.0, 1.0, 1.0],
     [1.0, 0.0, 1.0],
     [1.0, 1.0, 0.0]]
)
si = vq.PeriodicSystem(
    3, lattice,
    [vq.Atom(14, [0.0, 0.0, 0.0]),
     vq.Atom(14, [a / 4.0, a / 4.0, a / 4.0])],
)
basis = vq.BasisSet(si.unit_cell_molecule(), "sto-3g")

# GDF, multi-k 4x4x4 (closed-shell RHF or RKS / PBE).
res_rks = vq.run_periodic_job(
    si, basis,
    method="RKS", functional="pbe",
    jk_method="gdf",
    aux_basis="def2-universal-jkfit",
    kpoints=(4, 4, 4),
    output="si-pbe",
)
print(f"Si bulk PBE : {res_rks.energy:+.6f} Ha / cell")

For diamond C, swap to Atom(6, ...) and a = 3.567 A. Diamond’s parity reference at STO-3G with BIPOLE on SHRINK 8 is \(-74.8771393842\,\text{Ha}\) (delta vs CRYSTAL14: \(-0.145\,\text{mHa}\)).

Covalent crystals are also a clean target for the GPW / GAPW route (Lippert, Gerald, Hutter, Jürg, and Parrinello, Michele (1997). *A hybrid Gaussian and plane wave density functional scheme.* Molecular Physics, **92**(3), 477-487. doi:[10.1080/00268979709482119](https://doi.org/10.1080/00268979709482119)), which runs a full SCF today: the smooth pseudo-density on the FFT grid represents the bonding region efficiently, and this is where the plane-wave scaling win shows up first as the cell grows.

6.3 Metallic 3D: fcc aluminium with Fermi-Dirac smearing

Metals are the second classical stress test for periodic SCF: partial occupations at the Fermi level, no gap to anchor DIIS, and integer-occupation aufbau diverges. vibe-qc handles this with finite-temperature Fermi-Dirac smearing (Mermin, N. D. (1965). *Thermal Properties of the Inhomogeneous Electron Gas.* Physical Review, **137**(5A), A1441–A1443. doi:[10.1103/PhysRev.137.A1441](https://doi.org/10.1103/PhysRev.137.A1441)) on the multi-k BIPOLE RKS / RHF path. The driver auto-routes through BIPOLE when smearing_metallic=True.

import numpy as np
import vibeqc as vq

ANG2BOHR = 1.0 / 0.529177210903

# Al fcc primitive (1 Al), a = 4.05 A.
a = 4.05 * ANG2BOHR
lattice = (a / 2.0) * np.array(
    [[0.0, 1.0, 1.0],
     [1.0, 0.0, 1.0],
     [1.0, 1.0, 0.0]]
)
al = vq.PeriodicSystem(
    3, lattice, [vq.Atom(13, [0.0, 0.0, 0.0])],
)
basis = vq.BasisSet(al.unit_cell_molecule(), "sto-3g")

result = vq.run_periodic_job(
    al, basis,
    method="RKS", functional="pbe",
    jk_method="bipole",
    kpoints=(8, 8, 8),
    smearing_temperature="auto",   # ~0.01 Ha for metallic systems
    smearing_metallic=True,
    output="al-fcc-pbe",
)
print(f"Al fcc PBE   : E   = {result.energy:+.6f} Ha")
print(f"               F   = {result.free_energy:+.6f} Ha (= E - T*S)")
print(f"               E_f = {result.fermi_level:+.6f} Ha")

Two things to note. First, the meaningful thermodynamic quantity on a smeared metal is the free energy \(F = E - T S\), not the total energy; both are reported. Second, smearing_metallic=True relaxes the gap-closure safety guards that the standard SCF inserts; on a wide-gap insulator it should stay False (the default) so an accidental run on the wrong template fails loudly. For UKS open-shell metals, use the BIPOLE UKS route when you need per-spin Fermi-Dirac smearing today. The non-BIPOLE multi-k Ewald-UKS path still raises on positive smearing_temperature so spin-polarized metallic inputs do not run with silently integer occupations.

6.4 Open-shell metallic 3D: bcc lithium (BIPOLE UHF)

Light alkali metals are the simplest place to exercise the open-shell periodic SCF. BIPOLE is the AUTO open-shell route and the one that carries open-shell Fermi-Dirac smearing on the multi-k UKS / RKS path; open-shell GDF and GPW / GAPW also run UHF / UKS (GDF at Gamma and multi-k, GPW / GAPW at Gamma) when selected explicitly.

PeriodicSystem defaults to multiplicity 1. When an odd-electron primitive cell reaches open-shell SCF with that default, run_periodic_job selects the lowest parity-compatible multiplicity (normally a doublet); the same rule is applied after conventional-to-primitive reduction. An explicit multiplicity is preserved. Restricted RHF and RKS instead reject an odd-electron final cell and direct the calculation to UHF or UKS.

The run_periodic_job(method="ROHF", jk_method="gpw") route and its standalone run_periodic_rohf_gpw driver are available for a 3D Gamma-point spin-pure HF reference, but they keep integer 2/1/0 occupations and are not electronic-smearing routes. The standalone driver’s smearing_alpha option is nuclear smoothing, not a metal occupation temperature.

import numpy as np
import vibeqc as vq
from vibeqc.pbc_bipole_uhf import run_pbc_bipole_uhf

a = 6.6   # bohr, bcc Li (conventional cubic edge)
sysp = vq.PeriodicSystem(
    3, (a / 2) * np.array([[-1, 1, 1], [1, -1, 1], [1, 1, -1]]),   # bcc primitive
    [vq.Atom(3, [0.0, 0.0, 0.0])],
)
basis = vq.BasisSet(sysp.unit_cell_molecule(), "sto-3g")
kmesh = vq.monkhorst_pack(sysp, [2, 2, 2], use_symmetry=False)

opts = vq.PeriodicRHFOptions()
opts.lattice_opts.cutoff_bohr = 20.0   # 14 bohr fails the BIPOLE fold check
opts.use_diis = True
opts.max_iter = 40
opts.initial_guess = vq.InitialGuess.SAD

result = run_pbc_bipole_uhf(sysp, basis, kmesh, opts,
                            ewald_precision=1e-6)
print(f"E_total = {result.energy:+.6f} Ha   <S^2> = {result.s_squared:.4f}")

A molecular-limit open-shell BIPOLE input (an isolated Li atom in a 10-bohr box) lives at examples/periodic/input-bipole-li-uhf.py.

6.5 Mixed ionic-covalent 3D: ZnO wurtzite (hexagonal cell)

ZnO wurtzite mixes polar ionic character (Zn-O charge transfer) with covalent sp3 bonds in a non-orthogonal hexagonal cell (gamma = 120 deg). It is the canonical test that the periodic machinery does not secretly assume cubic axes.

import numpy as np
import vibeqc as vq

ANG2BOHR = 1.0 / 0.529177210903

# ZnO wurtzite (P6_3 m c, #186). a = 3.2495 A, c = 5.2069 A, u = 0.382.
a, c, u = 3.2495 * ANG2BOHR, 5.2069 * ANG2BOHR, 0.382
lattice = np.array([
    [a,            0.0, 0.0],
    [-a / 2.0,     a * np.sqrt(3.0) / 2.0, 0.0],
    [0.0,          0.0, c],
]).T

# Two formula units in the conventional hex cell.
unit_cell = [
    vq.Atom(30, lattice @ np.array([1/3, 2/3, 0.0])),     # Zn
    vq.Atom(30, lattice @ np.array([2/3, 1/3, 0.5])),
    vq.Atom( 8, lattice @ np.array([1/3, 2/3, u])),       # O
    vq.Atom( 8, lattice @ np.array([2/3, 1/3, 0.5 + u])),
]
zno = vq.PeriodicSystem(3, lattice, unit_cell)
basis = vq.BasisSet(zno.unit_cell_molecule(), "pob-tzvp")

# GDF would need a matching pob-tzvp-jk aux which is not yet shipped;
# BIPOLE handles pob-TZVP natively without an aux basis.
result = vq.run_periodic_job(
    zno, basis,
    method="RKS", functional="pbe",
    jk_method="bipole",
    kpoints=(4, 4, 4),
    output="zno-wurtzite-pbe",
)
print(f"ZnO wurtzite PBE : {result.energy:+.6f} Ha")

For wurtzite-class systems on pob-TZVP, BIPOLE is the operational default until the pob-TZVP JK aux basis is published (see density_fitting.md).

6.6 Van der Waals 3D: neon fcc (RKS / PBE + D3)

Closed-shell rare-gas solids are dominated by dispersion. The HF / pure-GGA part is mostly a sanity check on the periodic SCF; the cohesion comes from the D3 / D4 correction added on top (Grimme, Stefan, Antony, Jens, et al. (2010). *A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu.* Journal of Chemical Physics, **132**(15), 154104. doi:[10.1063/1.3382344](https://doi.org/10.1063/1.3382344), Grimme, Stefan, Ehrlich, Stephan, and Goerigk, Lars (2011). *Effect of the damping function in dispersion corrected density functional theory.* Journal of Computational Chemistry, **32**(7), 1456–1465. doi:[10.1002/jcc.21759](https://doi.org/10.1002/jcc.21759); Caldeweyher, Eike, Ehlert, Sebastian, et al. (2019). *A generally applicable atomic-charge dependent London dispersion correction.* Journal of Chemical Physics, **150**(15), 154122. doi:[10.1063/1.5090222](https://doi.org/10.1063/1.5090222)). Both BIPOLE and GDF work here; GDF is cheaper since the system is closed-shell.

import numpy as np
import vibeqc as vq

ANG2BOHR = 1.0 / 0.529177210903
a = 4.43 * ANG2BOHR   # Ne fcc conventional cubic
lattice = (a / 2.0) * np.array(
    [[0.0, 1.0, 1.0],
     [1.0, 0.0, 1.0],
     [1.0, 1.0, 0.0]]
)
ne = vq.PeriodicSystem(3, lattice, [vq.Atom(10, [0.0, 0.0, 0.0])])
basis = vq.BasisSet(ne.unit_cell_molecule(), "def2-svp")

result = vq.run_periodic_job(
    ne, basis,
    method="RKS",
    functional="pbe-d3bj",          # PBE + D3(BJ) dispersion
    jk_method="gdf",
    aux_basis="def2-svp-jk",
    kpoints=(2, 2, 2),
    output="ne-fcc-pbe-d3bj",
)
print(f"Ne fcc PBE+D3 : {result.energy:+.6f} Ha")

vdW solids are also where the opt-in nature of GAPW’s all- electron augmentation matters most: relative cohesive energies on the meV scale are easy to get wrong with too-coarse a smooth grid, so a CP2K parity sweep on Ne / Ar is part of the M3 acceptance set in the GAPW design doc.

6.7 1D: trans-polyene chain (Peierls distortion)

A 1D chain shows why dimensionality matters: Peierls distortion (R. E. Peierls, Quantum Theory of Solids, Oxford University Press (1955)) opens a gap in the otherwise metallic uniform chain. vibe-qc has no route with a consistent 1-D Coulomb gauge yet, so a PeriodicSystem(dim=1, ...) raises under GDF and under AUTO. Model the chain as a dim=3 cell with vacuum on the two transverse axes and sample k along the chain only. The uniform and the distorted geometries run on the same code path; check how the energy moves with the vacuum before comparing it with another code.

import numpy as np
import vibeqc as vq

# Trans-polyacetylene-like H-H chain. Column 0 is the chain axis;
# columns 1 and 2 are transverse vacuum, part of the model.
a       = 4.0    # bohr, lattice period
delta   = 0.20   # bohr, Peierls displacement (=0 for the uniform chain)
pad     = 15.0   # bohr, transverse vacuum

sysp = vq.PeriodicSystem(
    dim=3,
    lattice=np.diag([a, pad, pad]),
    unit_cell=[vq.Atom(1, [0.0,           0.0, 0.0]),
               vq.Atom(1, [1.4 + delta,   0.0, 0.0])],
)
basis = vq.BasisSet(sysp.unit_cell_molecule(), "sto-3g")

result = vq.run_periodic_job(
    sysp, basis,
    method="RHF",
    jk_method="gdf",
    kpoints=(8, 1, 1),
    output="h-chain-peierls",
)
print(f"H-chain (delta = {delta:+.2f} bohr) : {result.energy:+.6f} Ha")

The matching bandstructure / DOS workflow lives at examples/periodic/input-h-chain-bands.py; the uniform vs Peierls comparison at examples/periodic/input-h-chain-peierls.py.

6.8 2D: graphene sheet (semi-metal) and a graphene nanoribbon

Graphene at half-filling is the textbook 2D Dirac semimetal: the Fermi level sits exactly at the K-point band crossing, so any finite-temperature smearing > 0 introduces partial occupation in the pi* band. A Gamma-only single-point will still run cleanly because the cone is not at Gamma. On this route the (6, 6, 1) tuple below builds the classical shifted Monkhorst-Pack mesh, which does not contain K either; a Γ-centred mesh whose in-plane divisions are multiples of 3 does.

import numpy as np
import vibeqc as vq

# Graphene sheet, primitive hexagonal cell, a = 2.46 A in-plane.
# slab_2d takes only the two in-plane vectors: no vacuum, real-z atoms.
ANG2BOHR = 1.0 / 0.529177210903
a = 2.46 * ANG2BOHR

graphene = vq.slab_2d(
    [a,       0.0,                    0.0],
    [a / 2.0, a * np.sqrt(3.0) / 2.0, 0.0],
    [vq.Atom(6, [0.0,     0.0,                    0.0]),
     vq.Atom(6, [a / 2.0, a * np.sqrt(3.0) / 6.0, 0.0])],
)
basis = vq.BasisSet(graphene.unit_cell_molecule(), "def2-svp")

result = vq.run_periodic_job(
    graphene, basis,
    method="RKS", functional="pbe",
    jk_method="auto",                 # dim=2 -> SLAB_EWALD_2D
    kpoints=(6, 6, 1),                # in-plane mesh; 1 along the normal
    output="graphene-pbe",
)
print(f"Graphene PBE : {result.energy:+.6f} Ha / cell")

For an armchair graphene nanoribbon (effectively a 1D system), use a dim=3 cell with vacuum on the two transverse axes and sample k along the ribbon axis only, as for the chain in section 6.7; a dim=1 cell has no route.

Low-dimensional cells get no AUTO convergence profile

The AUTO convergence profiler’s covalent / ionic / metallic signals are 3D-only: they are gated on a tight-cell test that requires dim == 3. Every dim=1 and dim=2 cell therefore reports

profile: unknown
  - no low-dimensional branch (dim=2, ...) -- the covalent/ionic/metallic
    signals are 3D-only, so this cell is unclassified by construction,
    not judged inconclusive; a semimetallic or small-gap sheet needs an
    explicit smearing choice

and resolves smearing_temperature = 0. That default is wrong for a semimetal: graphene at half-filling cannot converge on integer occupations, and a small-gap sheet may not either. Pass smearing_method and smearing_temperature explicitly for such systems on routes that support them. Note that the explicit closed-shell slab GDF route currently rejects smearing outright (NotImplementedError), so a 2D semimetal has no smeared path there yet.

Band occupations on the multi-k slab GDF route are filled under one global Fermi level across the whole mesh, not per k point. When the bands overlap, that fill is fractional at the Fermi-degenerate group and the .out reports occupation class = fractionally occupied / smeared rather than a zero-temperature band gap.

A negative band gap is printed, not clamped

gap (indirect) is reported as CBM - VBM with no floor applied. If the value comes out negative, the band ordering across the mesh is inverted and the line carries a negative: band ordering inverted across the mesh annotation. Earlier releases clamped this at 0.000000 Ha, which made an inverted ordering indistinguishable from a genuine zero gap. Treat a negative gap as a diagnostic to investigate, not as a physical result.

6.9 Slab: MgO(001) 5-layer surface

A slab is a genuinely 2D periodic system: the in-plane lattice vectors are the surface mesh, and there is nothing periodic along the normal. Give vq.slab_2d those two vectors and the atoms at their real z. No vacuum is needed, and none is used: the third lattice column is synthesized for AO and spglib bookkeeping, and the SCF energy does not depend on it.

import numpy as np
import vibeqc as vq

ANG2BOHR = 1.0 / 0.529177210903
a = 4.211 * ANG2BOHR

# MgO(001) 5-layer slab, primitive 1x1 surface cell. No vacuum: the two
# in-plane vectors are the whole periodicity, and the layers keep their
# real z.
in_plane = a / np.sqrt(2.0)

layers = []
for k, layer_z in enumerate(np.linspace(-a, a, 5)):   # layer spacing a/2
    if k % 2 == 0:
        layers.append(vq.Atom(12, [0.0,      0.0,      layer_z]))
        layers.append(vq.Atom( 8, [in_plane / 2.0, in_plane / 2.0, layer_z]))
    else:
        layers.append(vq.Atom( 8, [0.0,      0.0,      layer_z]))
        layers.append(vq.Atom(12, [in_plane / 2.0, in_plane / 2.0, layer_z]))

slab = vq.slab_2d([in_plane, 0.0, 0.0], [0.0, in_plane, 0.0], layers)
basis = vq.BasisSet(slab.unit_cell_molecule(), "pob-tzvp")

result = vq.run_periodic_job(
    slab, basis,
    method="RKS", functional="pbe",
    jk_method="auto",                 # dim=2 -> SLAB_EWALD_2D
    kpoints=(4, 4, 1),                # in-plane mesh; 1 along the normal
    output="mgo-001-slab",
)
print(f"MgO(001) 5-layer : {result.energy:+.6f} Ha / cell")

Warning

The pre-v0.16 bulk GDF path was wrong on dim=2: it summed the Coulomb interaction as if the layers were stacked a3 apart, so the nuclear repulsion depended on the k-mesh and the total energy could be off by thousands of Hartree. The current explicit slab-GDF route is a different implementation with a finite 2D Coulomb metric and rigorous Parry-gauge (Parry, D. E. (1975). *The electrostatic potential in the surface region of an ionic crystal.* Surface Science, **49**, 433–440. doi:[10.1016/0039-6028(75)90362-3](https://doi.org/10.1016/0039-6028(75)90362-3)) one-electron terms. It is closed-shell only; analytic gradients ship since 2026-07-30.

Surface slabs are the use case where GPW / GAPW’s scaling tells: the smooth-grid Hartree-J cost is set by the cell volume, which scales linearly with slab thickness, while the BIPOLE direct-K cost grows quadratically with the number of in-plane atoms. GPW / GAPW run a full SCF today (Gamma RHF / ROHF / ROKS / RKS / UHF / UKS, multi-k pure-DFT RKS / ROKS / UKS), but GPW/GAPW do not accept a dim=2 slab. The general slab route is SLAB_EWALD_2D; closed-shell energy and relaxation work may explicitly use the slab-truncated GDF route (analytic gradients since 2026-07-30).

6.10 Method-choice summary by bonding regime

System type	Recommended first choice	Notes
Ionic 3D bulk (MgO, NaCl, LiH)	GDF (closed-shell) or BIPOLE (forces / opt)	Both cleanly handle the Madelung gauge.
Covalent 3D bulk (diamond, Si)	GDF for cost, BIPOLE for CRYSTAL parity	GPW / GAPW is a natural smooth-density fit and runs full SCF today.
Metallic 3D bulk (Al, Cu)	BIPOLE + smearing (smearing_metallic=True)	RKS/UHF/UKS BIPOLE support Fermi-Dirac; BIPOLE RHF remains gated.
Open-shell 3D (Li bcc, magnetic oxides)	BIPOLE UHF / UKS (AUTO)	GDF UHF / UKS (Gamma + multi-k) and GPW / GAPW (Gamma) also run open-shell; ROHF and ROKS are available on Ewald-3D (Gamma + multi-k), while maintained-preview GPW provides Gamma ROHF and Gamma or multi-k pure-DFT ROKS with integer occupations.
Mixed ionic-covalent (ZnO, TiO2, Al2O3)	BIPOLE on pob-TZVP, GDF on def2-*-jk	pob-TZVP-jk aux is on the roadmap.
Van der Waals 3D (Ne, Ar, molecular crystals)	GDF + dispersion correction	GPW / GAPW carry D3(BJ) too; their meV-cohesion CP2K sign-off is still broadening.
1D wires (H chain, polyacetylene)	GDF on a dim=3 cell with transverse vacuum	AUTO and the GDF drivers raise on dim=1 at Gamma and with a k mesh; public BIPOLE rejects dim=1.
2D sheets (graphene, h-BN)	SLAB_EWALD_2D (AUTO); explicit GDF for closed-shell energy scaling	Both are vacuum-free and a3-invariant. Slab GDF is full-tuple and zero-temperature; closed-shell analytic forces ship (optimize=True).
2D surface slabs	SLAB_EWALD_2D (AUTO); explicit GDF for closed-shell energies + relaxation	No vacuum needed. Plain closed-shell relaxation: explicit GDF with optimize=True. Freeze-mask/NEB/smearing workflows: vq.slab(..., periodic_z=True) and run dim=3.

7. Surface reactions: adsorption, TS search, vibrational analysis

Surface chemistry (adsorption energies, reaction barriers, transition states, intermediates, reaction-path scans) is the flagship workflow vibe-qc is built around. Bulk single-points are infrastructure for it; molecular SCF is infrastructure for it. This section walks the end-to-end stack: slab setup, adsorbate placement, geometry relaxation, NEB reaction path, vibrational analysis at the transition state.

7.1 Slab setup: dim = 2 is the physical model

There are two ways to make a slab look “non-periodic in the surface normal direction”. They are not equally correct.

• dim = 2 (the physical model, and the default). Give the two in-plane lattice vectors and put the atoms at their real Cartesian z. There is no vacuum gap. Build it with vq.slab_2d(a1, a2, atoms) or vq.slab(...), which returns dim=2. jk_method="auto" resolves to CoulombMethod.SLAB_EWALD_2D, the rigorous vacuum-free 2D Ewald gauge (Parry 1975, Parry, D. E. (1975). *The electrostatic potential in the surface region of an ionic crystal.* Surface Science, **49**, 433–440. doi:[10.1016/0039-6028(75)90362-3](https://doi.org/10.1016/0039-6028(75)90362-3); de Leeuw and Perram 1979, de Leeuw, S. W. and Perram, J. W. (1979). *Electrostatic lattice sums for semi-infinite lattices.* Molecular Physics, **37**(4), 1313–1322. doi:[10.1080/00268977900100951](https://doi.org/10.1080/00268977900100951)). The total energy is provably invariant to the synthesized third lattice column and the nuclear repulsion is k-mesh independent. Available for RHF / RKS / UKS at Gamma and multi-k. Analytic gradients in dim = 2 ship on the closed-shell slab-GDF route (2026-07-30): jk_method="gdf" with optimize=True (or the drivers’ compute_gradient=True) relaxes atoms on the slab analytic-gradient objective. Still not available in dim = 2: gradients on the direct slab_ewald_2d route itself, NEB / FD Hessians / freeze masks, Fermi-Dirac smearing, UHF, and the bulk J/K routes (bipole, gpw, gapw, rijcosx), all of which raise rather than return a wrong number.

• dim = 3 with thick vacuum (the plane-wave-code convention). Make the c-axis large enough that the slab images do not see each other (20 to 30 bohr of vacuum on each side), and run the whole calculation in 3D. Build it with vq.slab(..., periodic_z=True). Production finite-difference forces work through BIPOLE; analytic gradients remain a gated preview. Smearing is available, and the periodic ASE (Larsen, Ask Hjorth, Mortensen, Jens Jørgen, et al. (2017). *The atomic simulation environment - a Python library for working with atoms.* Journal of Physics: Condensed Matter, **29**(27), 273002. doi:[10.1088/1361-648X/aa680e](https://doi.org/10.1088/1361-648X/aa680e)) adapter (vibeqc.ase_periodic) accepts the geometry without special-casing. The cost is real: you pay for one extra Bloch sum in the normal direction, and a residual slab-slab dipole-image interaction shifts the absolute energy by an amount that grows with the surface dipole and shrinks only as the vacuum grows. It cancels in differences such as binding energies only to the extent that the dipole is the same on both sides.

Warning

Do not build a dim=2 cell and then hand it to a bulk J/K builder. A slab run through a 3D Coulomb gauge is treated as a crystal of sheets stacked a3 apart, so its energy depends on a3 and on the k-mesh. vibe-qc now raises rather than letting that happen. Before v0.15.32 it did not, and jk_method="auto" silently chose GDF.

Which to use. For single points and reference energies, use dim = 2: it is the correct physics and costs less. Explicit closed-shell slab GDF also supports fixed-cell relaxation. For freeze masks, periodic Gaussian NEB, open shell, or metallic-smearing single points, use dim = 3 + thick vacuum today and accept the image error. Do not combine positive smearing with BIPOLE optimization or periodic Gaussian NEB; both routes fail closed until their result schemas represent the Mermin free-energy objective. Periodic BIPOLE Hessians fail closed in every dimension.

import numpy as np
import vibeqc as vq

ANG2BOHR = 1.0 / 0.529177210903

def build_pt111_slab(n_layers=4, n_xy=2, vacuum_bohr=20.0):
    """Pt(111) n_xy x n_xy x n_layers slab. dim = 3 with thick vacuum."""
    a_pt = 3.92 * ANG2BOHR / np.sqrt(2.0)        # nearest-neighbour, bohr
    in_plane = a_pt * n_xy
    layer_spacing = a_pt * np.sqrt(2.0 / 3.0)    # FCC (111) interlayer

    # Hex 2D mesh in-plane, c-axis = slab + vacuum.
    lattice = np.array([
        [in_plane,            0.0,                      0.0],
        [in_plane / 2.0,      in_plane * np.sqrt(3.0) / 2.0, 0.0],
        [0.0,                 0.0,                      n_layers * layer_spacing + 2 * vacuum_bohr],
    ]).T

    atoms = []
    for layer in range(n_layers):
        z = layer * layer_spacing + vacuum_bohr
        # ABCABC stacking offset per layer.
        shift = (layer % 3) * np.array([in_plane / 3.0,
                                        in_plane * np.sqrt(3.0) / 6.0, 0.0])
        for i in range(n_xy):
            for j in range(n_xy):
                pos = (i * np.array([a_pt, 0, 0])
                       + j * np.array([a_pt / 2.0,
                                       a_pt * np.sqrt(3.0) / 2.0, 0])
                       + shift
                       + np.array([0, 0, z]))
                atoms.append(vq.Atom(78, list(pos)))

    return vq.PeriodicSystem(dim=3, lattice=lattice, unit_cell=atoms)

slab = build_pt111_slab(n_layers=4, n_xy=2, vacuum_bohr=20.0)
print(f"{len(slab.unit_cell)} Pt atoms in 2x2x4 slab")

For convenience you can also build a slab through ASE and convert:

from ase.build import fcc111
from vibeqc.ase_periodic import atoms_to_periodic_system

ase_slab = fcc111("Pt", size=(2, 2, 4), vacuum=10.0, a=3.92)
ase_slab.pbc = (True, True, True)   # the dim = 3 + vacuum convention
slab = atoms_to_periodic_system(ase_slab)

atoms_to_periodic_system currently requires pbc.sum() == 3, which is exactly the convention above. Slabs declared with pbc = (True, True, False) still fail closed at the boundary even though the dim = 2 slab-GDF gradient now exists (2026-07-30): mapping ASE’s 2D pbc onto the slab route’s synthesized-normal convention is a boundary-convention decision that has not been made, so the constraint stays until it is. Native dim = 2 relaxation goes through run_periodic_job(jk_method="gdf", optimize=True) instead.

7.2 Adsorbate placement

vibe-qc ships a native, pure-Python slab + adsorbate builder (vibeqc.slab + vibeqc.place_adsorbate) that does not depend on ASE at runtime. See slabs_and_adsorbates.md for the full reference. Quick form:

import vibeqc as vq

slab, info = vq.slab("Pt", "fcc", (1, 1, 1), n_layers=4, vacuum=10.0,
                     supercell=(2, 2), periodic_z=True)
# periodic_z=True -> dim=3 with a real vacuum gap. Freeze masks, NEB,
# open shell, and smeared single points use the dim=3 BIPOLE stack, but
# positive smearing cannot be combined with BIPOLE optimization or NEB.
# dim=2 analytic forces exist only on the closed-shell slab-GDF route
# (run_periodic_job(jk_method="gdf", optimize=True)).
slab_with_h2 = vq.place_adsorbate(slab, "H2", site="top")

Built-in adsorbate library: H, H2, N2, O2, CO, OH, NH, NH3, H2O, CH4. Site names: "top" / "bridge" / "hollow" / "fcc-hollow" / "hcp-hollow" / "long-bridge" / "short-bridge", plus explicit (x, y) coordinates. info.bottom_layer_indices(n) is the convenience accessor that pairs with the freeze mask in § 7.3.

The ASE bridge is still available if you prefer ASE builders:

from ase.build import add_adsorbate, fcc111, molecule

ase_slab = fcc111("Pt", size=(2, 2, 4), vacuum=10.0, a=3.92)
add_adsorbate(ase_slab, molecule("H2"), height=2.0, position="ontop")
ase_slab.pbc = (True, True, True)
slab_with_h2 = vq.ase_periodic.atoms_to_periodic_system(ase_slab)

For charged adsorbates, set charge= and multiplicity= on atoms_to_periodic_system (or, with the native builder, on the PeriodicSystem constructor directly).

7.3 Geometry relaxation

Fixed-cell BIPOLE atom relaxation goes through vibeqc.bipole_optimize.relax_atoms, using finite-difference forces by default. A freeze_indices= mask is available for fully 3D periodic systems, and VibeQCPeriodic supplies the ASE Calculator surface for the same fully periodic boundary. Explicit closed-shell slab GDF instead supports unmasked dim=2 analytic-gradient optimization through run_periodic_job(jk_method="gdf", optimize=True).

Do not use the former Pt/POB example as a runnable recipe. The bundled POB coverage for Pt is ECP-paired and the public BIPOLE all-electron guard rejects that combination. Also, relax_atoms has no dispersion objective, so a compound functional spelling such as "pbe-d3bj" is not a substitute for a separately implemented dispersion force. Choose a basis with a complete all-electron record for every element and a BIPOLE-supported functional. Long-range-corrected range-separated hybrids, VV10/nonlocal correlation, and double hybrids fail closed. Converge the k-mesh and lattice cutoffs for the actual system.

7.4 Reaction-path search via climbing-image NEB

vibe-qc ships a native Nudged Elastic Band driver (Henkelman, Graeme and Jónsson, Hannes (2000). *Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points.* Journal of Chemical Physics, **113**(22), 9978-9985. doi:[10.1063/1.1323224](https://doi.org/10.1063/1.1323224)), vibeqc.run_neb, that handles molecular and bounded periodic systems through one entry point. See neb.md for the executable light-element periodic example.

Periodic Gaussian NEB fails closed before image evaluation unless the system is dim=3, every element has a supported all-electron basis, and no periodic dispersion correction is requested. POB basis names, ECP metadata or paired bases, dim=1/dim=2, and dispersion_params are rejected. These guards mean the Pt/POB/D3 sketch formerly shown here was not a runnable production route.

Two flavours of NEB matter for surface chemistry:

• Plain NEB (climbing_image=False). Relaxes the chain onto the MEP but does not pin the highest-energy image to the saddle; the TS is interpolated from the band, which is approximate.

• Climbing-image NEB (climbing_image=True, Henkelman, Uberuaga, Jonsson 2000, Henkelman, Graeme, Uberuaga, Blas P., and Jónsson, Hannes (2000). *A climbing image nudged elastic band method for finding saddle points and minimum energy paths.* Journal of Chemical Physics, **113**(22), 9901-9904. doi:[10.1063/1.1329672](https://doi.org/10.1063/1.1329672)). The highest-energy image is pushed along the tangent to climb up to the saddle point exactly. Use this for any production barrier number.

The ASE NEB on top of the molecular VibeQC calculator is still the right pattern for molecule-in-vacuum reaction paths (precedent: examples/workflows/input-nh3-umbrella-neb.py for NH3 umbrella inversion). For periodic systems, use native run_neb only inside its declared all-electron, no-dispersion, 3D Gaussian envelope.

7.5 Transition-state characterisation: imaginary modes

A true first-order saddle is normally certified by one imaginary normal mode along the reaction coordinate. Periodic BIPOLE Hessian requests now fail closed before SCF because the former construction differentiated an isolated unit-cell Hamiltonian rather than the executed periodic energy. Do not use a dim=3 vacuum cell or a frozen subset as a workaround. At present, obtain a periodic Hessian from a separately validated route before attaching ZPE or thermal corrections to a BIPOLE NEB barrier.

7.6 Pt(111) workflow status

The Pt(111) geometry builders remain useful for constructing and inspecting candidate structures. There is currently no supported end-to-end native Pt/POB/D3 recipe: the basis/ECP guard rejects that Gaussian combination, periodic Gaussian NEB rejects dispersion, and periodic BIPOLE Hessians fail closed. Treat these as explicit method boundaries, not knobs to bypass.

7.7 Recommended defaults for surface-reactions runs

Knob	Default	Reason
jk_method	"bipole" (dim=3 force workflows) / "auto" (dim=2 single points)	BIPOLE is the production FD force route for closed- and open-shell periodic systems, and it needs dim=3. On a dim=2 slab, AUTO selects SLAB_EWALD_2D.
dim	2 for physical slab single points and supported closed-shell GDF relaxation; 3 with converged vacuum for smearing, freeze masks, and periodic Gaussian NEB	Slab GDF analytic gradients ship. Periodic Gaussian NEB remains 3D-only.
basis	A basis with a complete all-electron record for every element	POB and ECP-paired bases are rejected by periodic Gaussian NEB; converge the chosen basis for the target chemistry.
functional and dispersion	A BIPOLE-supported semilocal, global-hybrid, or validated screened-hybrid functional for relaxation/NEB; dispersion only on routes that document matching forces	Long-range-corrected range-separated hybrids, VV10/nonlocal correlation, and double hybrids fail closed. functional="pbe-d3bj" is not a libxc functional, and periodic Gaussian NEB rejects dispersion_params before image evaluation.
kpoints	(N, N, 1) with N chosen so \(\lvert\mathbf{k}\rvert\cdot \lvert\mathbf{a}_\parallel\rvert \approx 2\pi / (10\,\text{bohr})\)	Standard slab Brillouin-zone sampling; one k-point in the normal direction.
cutoff_bohr	Converge for the chosen basis and system	The required real-space extent is basis and geometry dependent; do not inherit the former Pt/POB sketch’s cutoff.
smearing_temperature	0 for BIPOLE relaxation and periodic Gaussian NEB; "auto" only for fixed-geometry metallic single points	The SCF drivers support smeared metals, but the BIPOLE optimizer and periodic Gaussian NEB reject positive temperature until their outputs distinguish the Mermin objective from internal energy.
relax_atoms.conv_tol_grad	\(5 \cdot 10^{-4}\) Ha/bohr (production: \(1 \cdot 10^{-4}\))	Standard tightness; ASE’s BFGS fmax = 0.05 eV/A equivalent.
NEB images	7 (5 intermediate)	Enough resolution to resolve a single saddle; expand to 9-11 for multi-step paths.
NEB optimiser	quick-min (default in vibeqc.run_neb) + climbing_image=True	Quick-min is the velocity-Verlet-style integrator the native NEB driver uses; robust on noisy SCF gradients and parallelises trivially per image.
Reaction-coord scan (cheap first pass before NEB)	vibeqc.relaxed_scan(system, basis, ("bond", i, j), values, freeze_indices=...)	Single-coordinate sweep with per-step warm-start; resolves the rough barrier shape in O(N_points) SCFs without an NEB band. Use it to choose n_images and identify which atoms move.

7.8 Current caveats and gaps

Items marked LANDED below reflect the historical snapshot. New work belongs in the current core repository.

1. VibeQCPeriodic Calculator wrapper. LANDED. Use atoms.calc = VibeQCPeriodic(...) for fully periodic ASE Atoms. Direct pbc.sum() != 3 ASE conversion remains outside this wrapper.

2. dim = 2 production forces. LANDED for explicit closed-shell GDF. run_periodic_job(jk_method="gdf", optimize=True) uses the analytic slab gradient. The direct SLAB_EWALD_2D route itself, slab open shell, freeze masks, and NEB remain unsupported.

3. dim = 2 Fermi-Dirac smearing. Open. smearing_temperature is not supported by the slab-GDF envelope, while GPW / GAPW / RIJCOSX do not accept a slab. A metallic dim = 2 slab therefore cannot be smeared today. Use dim = 3 + thick vacuum for metal surfaces that need it.

4. Atom-freeze mask in relax_atoms. LANDED 2026-05-25 in commit eb83c6c7. Use relax_atoms(..., freeze_indices=[...]) (see § 7.3); SlabInfo.bottom_layer_indices(n) is the convenience accessor.

5. Periodic NEB native to vibe-qc. LANDED 2026-05-25 in commits d47a3b55 … 0395159f (NEB increments 1-5). Public API: vibeqc.run_neb(reactant, product, basis, ..., kpoints=..., freeze_indices=..., climbing_image=True) returning NEBResult with .write_qvf(path) for vibe-view rendering. Reference: neb.md.

6. Periodic BIPOLE Hessian. Open. Requests fail closed because the former construction did not differentiate the executed periodic Hamiltonian; a frozen-index option cannot repair that method boundary.

7. Surface-dipole correction. Open (lowest priority). Polar slabs (terminations with net dipole, e.g. polar ZnO surfaces) need a dipole correction to cancel the spurious slab-slab image interaction in the dim = 3 + vacuum convention. The workaround is symmetric (non-polar) slab termination.

The remaining open gaps include slab smearing/open shell, broader NEB basis and dispersion support, a correct periodic Hessian, and polar-slab dipole corrections. The design_periodic_gapw.md plane-wave route in v0.10.x is the other half of the long-term answer because slab supercells with hundreds of atoms are where GPW / GAPW’s \(\mathcal{O}(N_g \log N_g)\) Hartree-J scaling overtakes both BIPOLE and GDF.

8. Reproducibility and provenance

vibe-qc records every routing choice in the run output so a paper or a later re-run can reconstruct the calculation without relying on AUTO heuristic versions:

• The .out banner names the resolved method (never the literal "auto").

• The .system manifest records Ewald omega (BIPOLE), aux basis name + linear-dependence threshold (GDF), and the plane-wave cutoff plus grid shape (GPW / GAPW).

• The .references and .bibtex siblings carry the per-method citations through the citation database. Pisani-Dovesi-Roetti 1988 / Saunders 1992 / Dovesi 2014 / Erba 2023 fire on the registered BIPOLE route; Whitten 1973 / Dunlap 1979 / Eichkorn 1995 / Weigend 2006 / Sun 2017 fire on the GDF route; Lippert-Hutter 1999 / VandeVondele 2005 / Krack-Parrinello 2000 / Blochl 1994 / Lin 2016 fire on the GPW / GAPW route.

If you publish work that used vibe-qc, the references block in the .out file (plus the BibTeX sibling) is the copy-paste surface; cite the method paper for the route you actually ran, not just the catch-all CRYSTAL or CP2K citation.

9. Roadmap pointers

• BIPOLE Phase 6a / 6b / 7: redesign the quartet far-pair branch around the exact three-translation Fock domain before any activation, complete the higher-l compute_ext_el_spheropole, and pin the numerical parity sign-off on the 15-system demo suite. Status in bipole.md.

• GDF open-shell: native UHF / UKS drivers have landed at Gamma (run_pbc_gdf_uhf / run_pbc_gdf_uks) and multi-k (run_kuhf_periodic_gdf / run_kuks_periodic_gdf). Analytic gradients cover the production rsgdf Γ ladder (RHF/RKS/UHF/UKS), the multi-k KRHF/KRKS/KUHF/KUKS drivers, and run_periodic_job(optimize=True) relaxes GDF-routed jobs on that same GDF objective; remaining gradient work is stress/variable-cell plus the guarded envelopes (screened fits, smearing, IBZ, range separation).

• GPW / GAPW: M2 (Gamma SCF) and M3 (GAPW dual grid, multi-k RKS, open-shell) have landed; M4 metallic refinements with smearing and M5 ACE-decorated K + on-demand PAW dataset fetcher remain. Track in design_periodic_gapw.md.

• RIJCOSX: Gamma RHF, vacuum-padded Gamma RKS/UHF/UKS, and true multi-k RHF / RKS / UHF / UKS ship, including the post-convergence one-center replacement on every route, weighted full-BZ quadratures, and HSE-type screened hybrids (hse06-class, closed-shell) on the multi-k backend. The screened exchange uses the physical erfc kernel with its finite G=0 mode included (the VASP/CRYSTAL convention shared with the BIPOLE screened route); no exxdiv shift applies. Remaining work includes tight-cell parity for the dedicated Gamma K, open-shell multi-k screened hybrids, and the periodic analytic gradient.

• Top-level RSGDF and CFMM: enum reserved; the standalone jk_method="rsgdf" / "cfmm" routes are not implemented. (The internal gdf_method="rsgdf" backend that the GDF route uses for tight ionic cells is shipped and is the production default, see § 3.2.1.)

See also

• bipole.md: BIPOLE driver reference (phases, options, gradients, optimisation).

• density_fitting.md: molecular RIJ / RIJK / RIJCOSX and the JKBuilder dispatch that GDF reuses.

• multi_k_scf.md: multi-k SCF surface (BIPOLE and GDF), smearing for metals, and the multi-k GDF roadmap.

• ewald.md: the Ewald primitives that the one-electron and the long-range J builds share.

• k_points.md: Monkhorst-Pack mesh generation and BZ sampling.

• crystal_lattices.md: the 14 Bravais lattices with worked examples.

• citations.md: how method routing produces the per-run references block.

• design_native_gdf.md: design walkthrough for the native periodic GDF Lpq construction (modrho compensation, aux conditioning, Cholesky robustness).

• design_periodic_gapw.md: GAPW design doc, milestone breakdown, PAW dataset license analysis.

---

[1]

Γ-CCM and χ-CCM are distinct CCM construction approaches: union-and-weight/Wigner-Seitz integral weighting and finite-translation-group characters, respectively. A declared common exchange-q=0 convention is a comparison constraint, not an identification of the constructions. Their distinction is not a choice of Coulomb kernels, and equality for a specified operator and route must be established rather than inferred from the labels.