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.

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. Two independent experimental ab-initio CCM lines ship in the current release: Γ-CCM / aiccm2026dev-a (direct-torus, HF→CCSD(T)/DLPNO) and χ-CCM / aiccm2026dev-b (finite-character, SCF + finite-torus correlation). Both are under active head-to-head comparison on a 28-system benchmark. 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.

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
from vibeqc import Atom, PeriodicSystem

# 1D atomic wire (H-H along x). Columns 1 and 2 are vacuum.
sys_1d = PeriodicSystem(
    dim=1,
    lattice=np.diag([4.0, 30.0, 30.0]),
    unit_cell=[Atom(1, [0, 0, 0]),
               Atom(1, [0, 0, 1.4])],
)

# 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, conventional cubic cell, bohr).
sys_3d = PeriodicSystem(
    dim=3,
    lattice=np.diag([5.4, 5.4, 5.4]),
    unit_cell=[Atom(14, [0, 0, 0]),
               Atom(14, [0.25, 0.25, 0.25])],
)

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. BIPOLE 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).

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.

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)

# Independent finite cyclic-cluster RHF/RKS (experimental):
result = vq.run_periodic_job(sysp, basis, method="RKS", functional="pbe0",
                             jk_method="aiccm2026dev-b",
                             kpoints=(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, for RHF / RKS / UKS. Every bulk J/K builder raises on a slab rather than returning an a3-dependent energy.

• dim=1 / 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 representation labels for the direct-torus and finite-character (Γ-centred character-mesh) forms of one declared finite-BvK-torus Hamiltonian. They are not interchangeable unless the exchange-q=0 convention is matched.

jk_method=	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, no smearing, no UHF yet. Aliases: "slab", "ewald_2d".	Production (energies)
"gdf"	Native Gaussian density fitting. Closed-shell RHF / RKS / hybrids at Gamma in 1D / 3D and multi-k via run_krhf_periodic_gdf / run_krks_periodic_gdf; open-shell UHF / UKS at Gamma (run_pbc_gdf_uhf / run_pbc_gdf_uks, with per-spin Fermi-Dirac smearing) and multi-k (run_kuhf_periodic_gdf / run_kuks_periodic_gdf). Raises on dim=2: it is a bulk builder and would treat a slab as a crystal of sheets. AUTO still picks BIPOLE for open-shell by design.	Production
"aiccm2026dev-b"	χ-CCM, the finite-character (Γ-centred character-mesh) CCM implementation. aiccm_lattice_extension=(N1,N2,N3) defines the real cyclic dimensions and derives the equivalent character net; aiccm_wigner_seitz_shells=s requests odd dimensions 2s+1. aiccm_backend= selects "four_center", "ri", or "rijcosx". See the derivation.	Experimental, opt-in
"bipole"	an Ewald J-split; multi-k RHF / UHF / RKS / UKS; production FD forces plus atom + cell optimisation; analytic gradients remain a gated preview. 3D only on the Ewald-J split; dim=1 falls back to direct-truncated. Raises on dim=2: use "auto" / "slab_ewald_2d", or build a dim=3-with-vacuum cell (vq.slab(..., periodic_z=True)) when you need forces.	Production (3D)
"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. Runs a full SCF: Gamma RHF / RKS / UHF / UKS and multi-k RKS, with forces, FD Hessian, closed-shell DFT+U (RHF/RKS, including multi-k RKS), smearing, and D3(BJ). Open-shell +U uses BIPOLE today. Emits GAPWExperimentalWarning.	Production (opt-in warning)
"gapw"	GPW plus per-atom Gaussian augmentation on a fine radial grid for all-electron accuracy. Same SCF coverage as GPW (Gamma RHF / RKS / UHF / UKS, multi-k RKS, forces, FD Hessian, closed-shell DFT+U, smearing, D3(BJ)). Open-shell +U uses BIPOLE today. Emits GAPWExperimentalWarning.	Production (opt-in warning)
"rijcosx"	RI-J plus seminumerical chain-of-spheres exchange K (Neese 2009) 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.	Implemented, experimental periodic COSX
"rsgdf"	Top-level route not implemented (enum-reserved): range-separated GDF as a standalone jk_method (Ye and Berkelbach, 2021). 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. Enum-reserved; not implemented.	Roadmap

Two further periodic routes ship today but are reached through dedicated entry points rather than jk_method:

• 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 (restricted open-shell HF): 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. See § 3.8.

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) and is cross-validated out-of-process against a program from its own method family: GDF ↔ PySCF, BIPOLE ↔ CRYSTAL, GPW / GAPW ↔ CP2K / GPAW. 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. In vibe-qc the equivalent is the Ewald-3D J-split

  \[ 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 via shell-pair Schwarz screening (the same machinery that the molecular four-index path uses).

The CRYSTAL paper trail traces back to Saunders, Dovesi, Roetti et al. (1992) for the original BIELET / BIPOLE multipole-far-pair construction, and to Dovesi et al. (2014) and Erba et al. (2023) for the modern reference of the EXT EL-POLE / EXT EL-SPHEROPOLE conventions vibe-qc is targeting. 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)\).

• 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 -> optional level-shift -> diagonalisation -> optional MOM occupied-subspace reorder -> density rebuild.

• Multi-k acceleration in the dense-k regime. Symmetry-reduced Monkhorst-Pack inputs are accepted and internally expanded back to the full mesh until true IBZ orbit expansion lands.

• Production finite-difference atomic forces for all four flavours, plus atomic and cell relaxation. The gradient-based cell optimizer finite-differences the SCF energy over the 6 Voigt strain components; 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}\)

The 15-system CRYSTAL parity demo suite lives at tests/demos/ and the matching CRYSTAL14 reference inputs ship beside the vibe-qc inputs.

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; BIPOLE remains the AUTO open-shell default and the multi-k hybrid route.)

• Hybrid DFT works out of the box. Pass any libxc hybrid in the RKS / UKS driver; HF exchange shares the direct-K machinery.

• Production FD forces + cell optimisation. Atom optimisation in fractional coordinates, lattice optimisation via an FD strain-gradient of the SCF energy.

• Direct CRYSTAL parity. The shared-Ewald-\(\alpha\) gauge means every CRYSTAL14 reference number can be reproduced inside one per-mHa target without conversion fudge factors.

• 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 on the Ewald-J path. With the auto default, a dim=1 wire stays on the legacy direct-only branch (valid for vacuum-padded cells), and a dim=2 slab raises: use SLAB_EWALD_2D. Passing use_ewald_j_split=True on a 1D wire or 2D slab raises.

• Cost scales with the real-space cutoff. Direct \(K\) at large cutoff is the dominant per-iteration cost. The Saunders-Dovesi multipole-far-pair branch that would amortise this is wired through Phase 5b but not yet auto-engaged inside the Fock build (Phase 6a in the BIPOLE phase status).

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

2.7 Citations

• V. R. Saunders, R. Dovesi, C. Roetti, M. Causa, N. M. Harrison, R. Orlando, C. M. Zicovich-Wilson, CRYSTAL92 User’s Manual, University of Torino, 1992. Original BIELET / BIPOLE multipole-far-pair 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; Dunlap 1979; Eichkorn et al. 1995; 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: 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, 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 1D / 3D (a dim=2 slab raises) via run_rhf_periodic_gamma_gdf (and the runner-level run_periodic_job(..., jk_method="gdf")).

• Multi-k RHF / RKS in 1D / 3D (a dim=2 slab raises) via run_krhf_periodic_gdf(..., kmesh, ...) and run_krks_periodic_gdf(..., kmesh, ...). A (1, 1, 1) mesh delegates to the Gamma-GDF path; any other Monkhorst-Pack mesh runs the full multi-k loop with k-dependent 3c / 2c blocks, Bloch-phase assembly, and multi-k DIIS. The multi-k Ewald-gauge fix brings a multi-k mesh to within 0.000 mHa of the Gamma energy in the molecular limit.

• 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 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 1D wires and 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. See vibeqc.pbc_gdf for the compensated-cell GDF path that will eventually subsume this branch.

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: 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), 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. See HANDOVER_GDF_V0_11_2026_05_29.md § 2026-05-29 for the diagnostic chain.

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: 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. 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). Γ k-meshes that fall back to the legacy molecular-limit GDF driver (open options such as level_shift / fock_mixing / smearing) raise NotImplementedError rather than silently dropping the threshold.

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 1D and 3D. The same code path drives wires and bulk. Slabs (dim=2) are not on this route: GDF is a bulk builder, and it now raises on a slab rather than treating it as a crystal of sheets. Use SLAB_EWALD_2D (what jk_method="auto" picks for dim=2).

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

• Dimensionality clean in 1D. For a wire there is no 3D Madelung problem to invent gauge for; the compensated lattice sum converges in real space and the route is genuinely single-shot.

3.4 Disadvantages

• AUTO routes open-shell to BIPOLE. Native open-shell GDF does ship (run_pbc_gdf_uhf / run_pbc_gdf_uks at Gamma, with per-spin Fermi-Dirac smearing; run_kuhf_periodic_gdf / run_kuks_periodic_gdf multi-k) and is reachable with an explicit jk_method="gdf" on a UHF / UKS job. AUTO keeps BIPOLE as the open-shell default by design, so you opt into GDF open-shell deliberately.

• No analytic gradients on the periodic GDF path yet. Molecular RIJK has had its analytic gradient bug closed since v0.8.0 (see density_fitting.md), but the periodic gradient is still BIPOLE-only.

• Gamma-only gdf_method="compcell" cannot do tight ionic cells. On compact cells with Z > 1 atoms (e.g. rocksalt LiH primitive FCC, 115 bohr³) the Gamma-only compcell Hartree cannot resolve AO-pair overlap between periodic images and the SCF converges to a non-physical fixed point (LiH/STO-3G: +579.8 Ha at HF, +1172.6 Ha at UKS-PBE, vs PySCF’s -8.24 Ha). The drivers warn on entry and the energy-sanity guard rejects the converged garbage: run_pbc_gdf_uhf / run_pbc_gdf_uks raise a RuntimeError (bypass with check_energy_sanity=False for diagnostics); run_pbc_gdf_rhf warns and tags the result backend +SANITY_FAILED. Because of this, the open-shell Gamma drivers default to gdf_method="rsgdf" since 2026-07-09 (on the broken LiH fixture the rsgdf defaults land at -8.33 Ha UHF / -8.23 Ha UKS-PBE, at the documented few-mHa diffuse-Li fitting floor from PySCF); compcell stays selectable for development and keeps the warning + guard. A KS run through run_pbc_gdf_uks with options=None also defaults to PeriodicKSOptions now, so the periodic-Becke-grid/torus-density XC pairing engages by default (the old PeriodicRHFOptions fallback silently evaluated XC in the v0.8.x molecular-grid convention: -7.9638 Ha on the same fixture). run_pbc_gdf_rhf keeps its compcell default (its warn-and-tag contract is pinned); use the multi-k route (validated at µHa parity on LiH FCC at kmesh=(2,2,2)) or gdf_method="rsgdf" at Gamma for such cells. Note the production runner already dispatches Gamma GDF jobs with gdf_method="rsgdf"; the compcell default only concerns direct API calls to the closed-shell driver.

• Auxiliary basis coverage. Standard def2-*-jkfit aux bases 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 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.

3.5 When to use GDF

• You are doing closed-shell RHF / RKS on a dim=1 wire or a dim=3 bulk crystal. (A dim=2 slab goes to SLAB_EWALD_2D instead.)

• 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); pass jk_method="gdf" explicitly, since AUTO routes open-shell to BIPOLE.

• You do not need analytic forces on the periodic path (yet).

3.6 Quick start

import vibeqc as vq

# 2D graphene sheet, RKS / PBE, Gamma-only.
# A dim=2 slab does NOT go through GDF: leave jk_method at "auto" and it
# resolves to SLAB_EWALD_2D, the vacuum-free 2D Ewald gauge. Passing
# jk_method="gdf" here raises NotImplementedError.
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="auto",                          # -> SLAB_EWALD_2D
    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 (def2-jkfit family)”, 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: restricted open-shell HF

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 a periodic restricted open-shell Hartree-Fock driver on the Ewald-3D gauge:

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. 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) and Gaussian-augmented plane-wave (GAPW) methods, due to Lippert and Hutter and implemented as the core CP2K periodic-DFT engine, 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): 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) 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, and both emit a GAPWExperimentalWarning when you opt in (the maturity flag, not a NotImplementedError). 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=.

• 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. Hybrid / HF multi-k is not yet wired on this route and raises.

• Forces, the finite-difference Hessian, closed-shell DFT+U, Fermi-Dirac smearing, and D3(BJ) dispersion, threaded through the same runner plumbing the other periodic routes use. For +U, GPW/GAPW cover RHF/RKS at Gamma and multi-k RKS; open-shell UHF/UKS +U uses the BIPOLE route today.

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 and GAPW carry an opt-in GAPWExperimentalWarning 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.

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) is the planned hybrid-DFT accelerator and is being designed so it can be shared with the locality-exploitation track (ADMM, OT).

4.4 Disadvantages (intrinsic and current)

• Younger than GDF / BIPOLE. The route runs a full SCF but its CP2K parity sweep is still broadening, so it carries the opt-in GAPWExperimentalWarning. 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 is effectively GPL-v2 (CLAUDE.md section 1). 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. 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, 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); Chai and Head-Gordon, J. Chem. Phys. 128, 084106 (2008) (ωB97X); Yanai, Tew, Handy, Chem. Phys. Lett. 393, 51 (2004) (CAM-B3LYP).

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. It is the periodic analogue of the molecular RIJCOSX accelerator (Neese 2009) 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.

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

• One-cell RKS/UHF/UKS fail closed. The Gamma RKS/GDF path is not pair-resolved and the COSX bridge is a multi-k backend, so vibe-qc refuses to label a one-cell RKS/UHF/UKS calculation 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	direct only	raises	yes	yes	yes	yes	yes	yes	yes	yes	yes
GDF (native)	yes	raises	yes	yes	yes (Γ + multi-k)	yes	yes	yes (Γ + multi-k)	yes	no	no
GPW	limited	raises	yes	yes (Γ)	yes (Γ)	yes (Γ)	yes (Γ)	yes (Γ)	yes (RKS)	FD	FD
GAPW	limited	raises	yes	yes (Γ)	yes (Γ)	yes (Γ)	yes (Γ)	yes (Γ)	yes (RKS)	FD	FD
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	no	no	yes	ROHF	n/a	no	no	no	yes	no	no

Slabs have exactly one route. SLAB_EWALD_2D is the only method that accepts dim=2, and jk_method="auto" selects it for you. Every bulk builder raises there rather than returning an a3-dependent energy. The price today is that a slab gets energies but not analytic gradients, and no Fermi-Dirac smearing; for force-driven or metallic work, build a dim=3 cell with vq.slab(..., periodic_z=True).

“Direct only” for BIPOLE on 1D means the driver auto-falls back to DIRECT_TRUNCATED; this is adequate for vacuum-padded molecular-limit wires (~10 bohr of vacuum or more) but not for tight ones. For tight wires today, use GDF.

Notes on the matrix: GDF open-shell (UHF / UKS) ships at Gamma and multi-k but AUTO routes open-shell to BIPOLE by design, so you opt into it with an explicit jk_method="gdf". GPW / GAPW multi-k is pure-DFT RKS only (hybrid / HF multi-k raise on that route); their gradients and cell optimisation go through the finite-difference path the runner shares. 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) and ROHF are dedicated entry points, not jk_method routes (RSGAPW handles screened hybrids such as HSE06 / ωB97X / CAM-B3LYP; the ROHF “RHF” column reads as restricted open-shell HF on the Ewald-3D gauge, Gamma and multi-k).

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; multipole far-pair branch (Phase 6a) is the planned amortiser.
GDF	\(\mathcal{O}(N^2 M)\) via \(L_{P,\mu\nu}\)	\(\mathcal{O}(N^2 M)\) via \(L_{P,\mu\nu}\)	\(\mathcal{O}(N^2 M)\) for the \(L\) tensor	Aux size \(M \approx 3N\); cost is Cholesky-bound for the first SCF call, contraction-bound thereafter.
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_g is the plane-wave grid count.

5.3 Method-selection chart

                       +----------------------------+
                       | Periodic SCF, which route? |
                       +-------------+--------------+
                                     |
                            Is it a slab (dim=2)?
                                     |
                        +------------+------------+
                        | yes                     | no
                        |                         |
                **SLAB_EWALD_2D**       +---------+---------+
                (the only route;        |                   |
                 AUTO picks it;     Open-shell?        Closed-shell
                 energies only,     (UHF / UKS)             |
                 no gradients,          |          +--------+--------+
                 no smearing)       **BIPOLE**     |                 |
                                              1D wire          3D system
                                                  |                 |
                                             **GDF**          +-----+------+
                                                              |            |
                                                       Multi-k or hybrid   Large supercell
                                                              |            (hundreds of atoms)
                                                       **GDF**             |
                                                       (or **BIPOLE**      **GPW / GAPW**
                                                       if you need         (plane-wave J;
                                                       FD forces           BIPOLE/GDF still
                                                       or cell opt)        fine for hybrids)

The slab branch comes first because it is exclusive: no bulk builder accepts dim=2. If you need forces or smearing on a surface, do not look for another jk_method, build the cell as dim=3 with vq.slab(..., periodic_z=True) instead.

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.

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 and GDF 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 forces / cell relaxation, GDF for the cheapest single-point at large orbital basis.

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

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

a = 6.6   # bohr, bcc Li
sysp = vq.PeriodicSystem(
    3, np.eye(3) * a,
    [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 = 14.0
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}")

The full open-shell input 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. 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 is the cleanest demonstration of why dimensionality matters: the 3D Madelung problem vanishes, GDF runs natively without any compensating-charge machinery, and Peierls distortion opens a gap in the otherwise metallic uniform chain. Both the uniform and the distorted geometries should run on the same code path.

import numpy as np
import vibeqc as vq

# Trans-polyacetylene-like H-H chain. Column 0 is the periodic
# axis; columns 1 and 2 are vacuum (>= 25 bohr decouples images).
a       = 4.0    # bohr, lattice period
delta   = 0.20   # bohr, Peierls displacement (=0 for the uniform chain)

sysp = vq.PeriodicSystem(
    dim=1,
    lattice=np.diag([a, 30.0, 30.0]),
    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.

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 embedded in a 2D sheet of vacuum), the same code path holds with dim=1 and only the ribbon-axis lattice column physical.

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(-2.0 * a, 2.0 * a, 5)):
    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

This example previously read jk_method="gdf", # 2D: no 3D Madelung problem. That was wrong, and it is the bug this convention change fixes. GDF is a bulk builder: handed a dim=2 slab 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 was off by thousands of Hartree. GDF now raises on dim=2.

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 / RKS / UHF / UKS, multi-k RKS), but none of them accept a dim=2 slab: the production slab route is SLAB_EWALD_2D. A 2D-truncated Coulomb kernel for GDF and GPW, which would let large slabs use density fitting or the grid, is on the roadmap.

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 available on the Ewald-3D gauge.
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 (any zeta)	No 3D Madelung problem; BIPOLE direct only.
2D sheets (graphene, h-BN)	SLAB_EWALD_2D (what AUTO picks)	Vacuum-free 2D Ewald. Energies only: no gradients, no smearing.
2D surface slabs	SLAB_EWALD_2D (what AUTO picks)	No vacuum needed; the total is a3-invariant. For forces or smearing, use 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; de Leeuw and Perram 1979). 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. Not yet available in dim = 2: analytic gradients (so no relaxation, NEB, or FD Hessian), Fermi-Dirac smearing, UHF, and the bulk J/K routes (gdf, 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). Analytic gradients work through the BIPOLE driver, smearing is available, and the periodic ASE 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.16.0 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. For force-driven workflows (relaxation, NEB, FD Hessian) and for metallic slabs that need smearing, use dim = 3 + thick vacuum today, and accept the image error. Slab analytic gradients and slab smearing are the two gaps that will make dim = 2 viable end-to-end.

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) will be supported when the dim = 2 gradient driver lands. The ASE bridge exists to feed force-driven workflows, and forces on a dim = 2 slab do not ship yet, so there is nothing to gain from relaxing the constraint before then.

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. Force-driven
# workflows need it: slab analytic gradients in dim=2 are not shipped.
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 with FD BIPOLE forces

Atom relaxation goes through vibeqc.bipole_optimize.relax_atoms, which uses the finite-difference BIPOLE force path by default and L-BFGS-B in fractional coordinates (stable on general lattices). It supports all four SCF flavours (RHF / UHF / RKS / UKS) and any libxc functional.

import vibeqc as vq
from vibeqc.bipole_optimize import relax_atoms

slab_with_h2 = build_pt111_h2_initial()    # see § 7.2
kmesh = vq.monkhorst_pack(slab_with_h2, [3, 3, 1])

result = relax_atoms(
    slab_with_h2,
    basis_name="pob-tzvp",
    kmesh=kmesh,
    method="RKS",
    functional="pbe-d3bj",            # include dispersion for adsorbates
    cutoff_bohr=12.0,
    max_iter=60,
    conv_tol_grad=5e-4,               # Ha/bohr; tighten for production
)
print(f"Adsorbed E = {result.fun:+.6f} Ha   ({result.nit} steps)")
relaxed_slab = result.system          # PeriodicSystem at the optimum

For a slab + adsorbate, freezing the bottom layers during relaxation is the standard convention (bottom of slab represents the bulk and should keep its bulk geometry). relax_atoms(..., freeze_indices=[...]) natively does this: the (fixed, fixed) L-BFGS-B box-bounds implementation zeroes gradient components on frozen atoms so the reported |grad| reflects free DOFs only.

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. Force-driven
# workflows need it: slab analytic gradients in dim=2 are not shipped.
slab_with_h2 = vq.place_adsorbate(slab, "H2", site="top")

# Freeze the bottom 2 substrate layers; relax the top 2 + adsorbate.
frozen = info.bottom_layer_indices(2)

result = relax_atoms(
    slab_with_h2,
    basis_name="pob-tzvp",
    kmesh=vq.monkhorst_pack(slab_with_h2, [3, 3, 1]),
    method="RKS",
    functional="pbe-d3bj",
    freeze_indices=frozen,
    cutoff_bohr=12.0,
    max_iter=60,
    conv_tol_grad=5e-4,
)

The ASE bridge + FixAtoms constraint route still works (and is the right path if you are doing a mixed ASE workflow with custom optimisers), but the native freeze_indices= is the canonical form for the surface-catalysis pattern.

The “full ASE Calculator” wrapper for periodic systems is being built; until then periodic_forces is the helper you call by hand inside a custom calculator subclass (see § 7.8).

7.4 Reaction-path search via climbing-image NEB

vibe-qc ships a native Nudged Elastic Band driver, vibeqc.run_neb, that handles both molecular and periodic systems through one unified entry point (the periodic dispatch turns on when you pass kpoints=). The full reference is at neb.md; the surface-reactions short form:

import vibeqc as vq
from vibeqc.bipole_optimize import relax_atoms

# Endpoint A: H2 above Pt(111), relaxed (see § 7.3).
# Endpoint B: 2 H atoms in adjacent fcc / hcp hollow sites, relaxed.
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. Force-driven
# workflows need it: slab analytic gradients in dim=2 are not shipped.
endpoint_a = vq.place_adsorbate(slab, "H2", site="top")
endpoint_b = vq.place_adsorbate(
    vq.place_adsorbate(slab, "H", site="fcc-hollow"),
    "H", site="hcp-hollow",
)
frozen = info.bottom_layer_indices(2)
# (Endpoint relaxations elided; see § 7.3.)

result = vq.run_neb(
    endpoint_a, endpoint_b,
    basis="pob-tzvp",
    n_images=5,                         # intermediate; 7 total with endpoints
    method="RKS", functional="pbe",
    dispersion_params="d3bj",
    kpoints=(3, 3, 1),
    freeze_indices=frozen,
    interpolation="idpp",               # image-dependent pair potential
    climbing_image=True,
    conv_tol_force=1e-3,                # ~ 0.05 eV/A
    n_jobs=1,                           # periodic images are memory-heavy
)

result.write_qvf("h2-pt111-neb.qvf")    # vibe-view renders the band

print(f"TS image = {result.climbing_image_index}")
print(f"Barrier  = "
      f"{(result.energies[result.climbing_image_index]"
      f" - result.energies[0]) * 27.2114:.2f} eV")

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). 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, the native run_neb is the production path.

7.5 Transition-state characterisation: imaginary modes

A true first-order saddle has exactly one imaginary normal-mode frequency along the reaction coordinate. vibe-qc’s periodic runner ships a finite-difference Hessian behind the hessian=True flag:

result = vq.run_periodic_job(
    ts_slab, basis,                   # ts_slab from the climbing-image NEB
    method="RKS", functional="pbe-d3bj",
    jk_method="bipole",
    kpoints=(3, 3, 1),
    hessian=True,                     # 6N FD displacements
    output="ts-hessian",
)
print(f"Imaginary modes: {sum(1 for f in result.hessian.freqs_cm if f.imag)}")
print(f"Largest imaginary freq: {min(f.imag for f in result.hessian.freqs_cm):+.1f} cm^-1")

For a TS this should report exactly one imaginary mode; for a minimum, zero. The eigenvector of the imaginary mode is the reaction coordinate at the saddle. The zero-point energy and the thermodynamic free-energy correction (via the rigid-rotor / harmonic-oscillator partition function) come out of the same Hessian result, and are the right correction to apply when turning a bare \(\Delta E^\ddag\) into a \(\Delta G^\ddag\) that compares with kinetics experiment.

For large slabs the full 6N FD Hessian is expensive (one SCF + gradient per displacement, factor 6N total). The conventional shortcut is to freeze the bottom layers when building the Hessian, restricting the FD displacements to the adsorbate plus the top-layer atoms that participate in the reaction. The analytic-gradient FD Hessian already amortises this when the frozen atoms are excluded from the displacement list (an explicit hessian_options.frozen_indices argument is on the roadmap; the operational workaround today is to build the Hessian on a sub-system with only the moving atoms via a custom driver call).

7.6 Worked example: H2 dissociative adsorption on Pt(111)

End-to-end native skeleton. Build the slab, place H2, relax the initial state, relax the dissociated final state, run climbing- image NEB, characterise the TS. Everything uses the native vibe-qc surface-chemistry stack (no ASE-runtime dependency for the calculation itself).

import vibeqc as vq
from vibeqc.bipole_optimize import relax_atoms

# 1. Slab + endpoint geometries.
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. Force-driven
# workflows need it: slab analytic gradients in dim=2 are not shipped.
endpoint_a = vq.place_adsorbate(slab, "H2", site="top")
endpoint_b = vq.place_adsorbate(
    vq.place_adsorbate(slab, "H", site="fcc-hollow"),
    "H", site="hcp-hollow",
)
frozen = info.bottom_layer_indices(2)
kmesh  = vq.monkhorst_pack(endpoint_a, [3, 3, 1])

# 2. Relax both endpoints (FD BIPOLE forces + freeze mask).
res_a = relax_atoms(endpoint_a, "pob-tzvp", kmesh,
                    method="RKS", functional="pbe-d3bj",
                    freeze_indices=frozen,
                    cutoff_bohr=12.0, max_iter=80, conv_tol_grad=5e-4)
res_b = relax_atoms(endpoint_b, "pob-tzvp", kmesh,
                    method="RKS", functional="pbe-d3bj",
                    freeze_indices=frozen,
                    cutoff_bohr=12.0, max_iter=80, conv_tol_grad=5e-4)

print(f"Reaction energy (adsorbed -> dissociated): "
      f"{(res_b.energy - res_a.energy) * 27.2114:+.2f} eV")

# 3. Climbing-image NEB (native, single call).
neb = vq.run_neb(
    res_a.system, res_b.system,
    basis="pob-tzvp",
    n_images=5,
    method="RKS", functional="pbe",
    dispersion_params="d3bj",
    kpoints=(3, 3, 1),
    freeze_indices=frozen,
    interpolation="idpp",
    climbing_image=True,
    conv_tol_force=1e-3,
    n_jobs=1,
)
neb.write_qvf("h2-pt111-neb.qvf")
ts_slab = neb.path[neb.climbing_image_index]
barrier_eV = (neb.energies[neb.climbing_image_index]
              - neb.energies[0]) * 27.2114
print(f"Forward barrier (no ZPE) : {barrier_eV:+.2f} eV")

# 4. TS Hessian (FD; for large slabs use freeze_indices once Gap 5 lands).
ts_res = vq.run_periodic_job(
    ts_slab, vq.BasisSet(ts_slab.unit_cell_molecule(), "pob-tzvp"),
    method="RKS", functional="pbe",
    dispersion="d3bj",
    jk_method="bipole",
    kpoints=(3, 3, 1),
    hessian=True,
    output="h2-pt111-ts",
)
imag = [f for f in ts_res.hessian.freqs_cm if f.imag]
print(f"Imaginary modes at TS : {len(imag)} (expect 1)")
print(f"Reaction-coord freq   : {imag[0].imag:+.1f} cm^-1")

The older geometry-only fixture examples/ase_workflows/surface-h2-pt111-singlepoint.py remains useful as a committed Pt(111) + H₂ setup; this skeleton is the canonical native periodic-driver pattern for production runs.

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 single points; 3 (with 20-30 bohr vacuum, via periodic_z=True) for force-driven or smeared workflows	dim=2 is the correct physics and is a3-invariant, but slab analytic gradients and slab smearing are not shipped yet, and the BIPOLE FD optimizer route is 3D.
basis	pob-tzvp (oxides / metals) or def2-svp	pob-TZVP is the Bredow-group periodic standard; def2-SVP is the cheap default when GDF + an aux is needed.
functional + dispersion="d3bj"	pbe + d3bj via the dftd3 backend (commit 5dbfbfc3)	PBE handles the bulk of slab DFT literature; D3(BJ) is required for adsorption energies on transition metals. run_periodic_job(..., dispersion="d3bj", dispersion_backend="dftd3") is the production wiring; the inline functional="pbe-d3bj" form also works.
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	12-14	Tight enough for sub-mHa SCF convergence on pob-TZVP; matches the CRYSTAL14 parity baseline.
smearing_temperature	0 (insulator slabs) or "auto" (metallic slabs)	Pt(111), Cu(111), Ni(111), Au(111) all need smearing.
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

The honest list of what still needs work. Tracked in handovers/HANDOVER_SURFACE_REACTIONS.md; items marked LANDED below have been crossed off there too.

1. VibeQCPeriodic Calculator wrapper. Open. ASE’s Calculator subclass for periodic systems (atoms.calc = VibeQCPeriodic(...)) is the missing single-line ergonomics piece. Today the periodic ASE bridge is the function-style periodic_forces(atoms, basis, kpts=, functional=); wrapping it inside a Calculator subclass is straightforward custom code, and is the natural follow-on once the atoms_to_periodic_system pbc.sum() == 3 constraint is relaxed for dim = 2.

2. dim = 2 production forces. Open. dim = 2 energies ship through the vacuum-free SLAB_EWALD_2D gauge, but its analytic gradient is not written yet, so forces, relaxation, NEB, and the FD Hessian all raise on a slab. Bulk builders are no longer a fallback here: they raise on dim = 2 rather than return an a3-dependent number. The dim = 3 + thick-vacuum convention (vq.slab(..., periodic_z=True)) is the operational workaround.

3. dim = 2 Fermi-Dirac smearing. Open. smearing_temperature is wired through GDF / GPW / GAPW / multi-k RIJCOSX, none of which accept a slab, so a metallic dim = 2 slab 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. Hessian frozen_indices option for slab TS characterisation. Open. Right now the 6N displacement count is the bottleneck on large slabs.

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 (1, 2, 5, 6) are what the surface- reactions focus elevates on the roadmap; 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. Saunders 1992 / Dovesi 2014 / Erba 2023 will fire when the BIPOLE route is registered into routes.methods["bipole"]; 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: wire the multipole far-pair branch into the Fock build for asymptotic acceleration, complete the higher-l compute_ext_el_spheropole, and pin the CRYSTAL14 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); remaining work is the periodic analytic gradient on the GDF path.

• 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 ships; bringing the tight-cell COSX exchange K to Γ-folded GDF parity (and RKS / multi-k) is the remaining work.

• 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 representation labels at a declared exchange-q=0 convention, currently exchange_q0="bvk-ewald" for the production χ-CCM path. They are not two theories and not two Coulomb kernels.