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.

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. Column 2 is vacuum normal to the sheet.
sys_2d = PeriodicSystem(
    dim=2,
    lattice=np.diag([2.5, 2.5, 30.0]),
    unit_cell=[Atom(6, [0, 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:

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

jk_method="auto" (the default) picks 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:

`jk_method=`	What it dispatches to	Status
`"gdf"`	Native Gaussian density fitting, Gamma-only RHF / RKS / hybrids in 1D / 2D / 3D; multi-k GDF in 1D / 2D / 3D via `run_krhf_periodic_gdf` / `run_krks_periodic_gdf`.	Production
`"bipole"`	CRYSTAL-gauge Ewald J-split; multi-k RHF / UHF / RKS / UKS; analytic gradients and atom + cell optimisation; 3D only on the Ewald-J split, 2D / 1D fall back to direct-truncated.	Production (3D)
`"direct"`	4-centre real-space lattice sum with Schwarz screening. Valid only for vacuum-padded molecular-limit cells.	Diagnostic
`"fft_poisson"`	Legacy Ewald-split FFT-Poisson; known broken on ionic crystals (MgO RHF off by +241 Ha at v0.7.0).	Regression-only
`"gpw"`	Gaussian-on-grid collocation + FFT-Poisson Hartree-J on a smooth real-space grid. Pseudopotential-only.	Experimental (M2 of v0.10.x; J builder lands, SCF wiring in progress)
`"gapw"`	GPW plus per-atom Gaussian augmentation on a fine radial grid for all-electron accuracy.	Roadmap (M3 of v0.10.x)
`"rsgdf"`	Range-separated GDF (Ye and Berkelbach, 2021).	Roadmap
`"cfmm"`	Continuous fast multipole.	Roadmap

The rest of this page is the side-by-side treatment of the three methods that matter day-to-day: BIPOLE, GDF, and the emerging GPW / GAPW route.

2. BIPOLE: CRYSTAL-gauge Ewald J-split¶

2.1 Theory¶

BIPOLE is vibe-qc’s CRYSTAL-style 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; CRYSTAL mirrors this through COMMON/VRSMAD/. 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 CRYSTAL-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 -> CRYSTAL-style 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.
Analytic atomic gradients for all four flavours, plus atomic and cell relaxation (force-based stress tensor, L-BFGS-B over the 6 Voigt strain components). 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. The only vibe-qc route with multi-k SCF on UHF / UKS as well as RHF / RKS.
Hybrid DFT works out of the box. Pass any libxc hybrid in the RKS / UKS driver; HF exchange shares the direct-K machinery.
Analytic gradients + cell optimisation. Atom optimisation in fractional coordinates, lattice optimisation via the force-based stress tensor.
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, $\text{dim} < 3$ runs stay on the legacy direct-only branch (valid for vacuum-padded cells); 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 analytic gradients 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:

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.
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.
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 / 2D / 3D via run_rhf_periodic_gamma_gdf (and the runner-level run_periodic_job(..., jk_method="gdf")).
Multi-k RHF / RKS in 1D / 2D / 3D 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 / 2D systems 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.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, 2D, and 3D. The same code path drives wires, slabs, and bulk.
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 / 2D 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¶

No open-shell. UHF / UKS native-GDF drivers are not in v0.8 / v0.9; use BIPOLE.
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.
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 periodic system in any of 1D, 2D, 3D.
You want the cheapest production route on bulk for hybrid DFT and large basis (def2-SVP through def2-TZVP class).
You do not need analytic forces on the periodic path (yet).

3.6 Quick start¶

import vibeqc as vq

# 2D graphene-like sheet, RKS / PBE, Gamma-only GDF.
sysp  = build_graphene()                       # 2D PeriodicSystem
basis = vq.BasisSet(sysp.unit_cell_molecule(), "def2-svp")

result = vq.run_periodic_job(
    sysp, basis,
    method="RKS", functional="pbe",
    jk_method="gdf",
    aux_basis="def2-svp-jk",
    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.

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:

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})$.
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.
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 is the v0.10.x acceleration program. The design doc with the full milestone breakdown is at design_periodic_gapw.md. Today (2026-05-24):

Milestone	Status	Deliverable
M1	landed	`PlaneWaveGrid` infrastructure + Gaussian-on-grid collocation, behind an experimental flag
M2	in progress (J builder landed, SCF wiring in progress)	GPW Hartree-J via FFT-Poisson, Gamma-only, closed-shell LDA
M3	scheduled	GAPW dual-grid augmentation
M4	scheduled	Metallic GAPW with smearing (composes with the scf-mix chat)
M5 (stretch)	scheduled	ACE-decorated K for hybrids + PAW pseudopotentials with vqfetch-style on-demand fetcher

What is callable today lives in periodic_gapw_j.py: the GpwJBuilder orchestrates density collocation -> FFT-Poisson -> AO-basis projection into a Hartree-J matrix on the smooth real-space grid. It is unit-tested against the molecular Hartree energy and the Madelung-shift expectation but does not yet drive a full SCF: jk_method="gpw" and jk_method="gapw" raise NotImplementedError until M2c lands the SCF wiring and a CP2K subprocess oracle gives sub-mHa parity on Si bulk. Track progress at .release-status/v0.10.x/gapw.md.

4.3 Advantages (target, M2 onwards)¶

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)¶

Not user-facing yet. GPW SCF wiring is in flight; GAPW is one milestone behind.
Plane-wave cutoff is an extra knob. Sensible defaults will be derived from the tightest primitive exponent and a CP2K-style REL_CUTOFF safety factor (calibrating against REL_CUTOFF = 60 Ry is the M2 target).
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 first. The first wave is 3D periodic only; 2D slab and 1D wire GPW lands once dual-grid is solid.
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¶

Today: only for development of the route itself, or for J-only parity work against CP2K. Once M2 lands, the target use case is 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.

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.

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
BIPOLE	direct only	direct only	yes	yes	yes	yes	yes	yes	yes	yes	yes
GDF (native)	yes	yes	yes	yes	no	yes	yes	no	yes (closed-shell)	no	no
GPW	scheduled	scheduled	wiring	LDA target	no	LDA target	(M5 via ACE)	no	scheduled	scheduled	scheduled
GAPW	scheduled	scheduled	scheduled	scheduled	no	scheduled	(M5 via ACE)	no	scheduled	scheduled	scheduled

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

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? |
                       +-------------+--------------+
                                     |
                  +------------------+-----------------+
                  |                                    |
              Open-shell?                         Closed-shell
              (UHF / UKS)
                  |                                    |
              **BIPOLE**                  +------------+------------+
                                          |                         |
                                  1D / 2D system               3D system
                                          |                         |
                                     **GDF**                  +-----+------+
                                                              |            |
                                                       Multi-k or hybrid   Large supercell
                                                              |            (hundreds of atoms)
                                                       **GDF**             |
                                                       (or **BIPOLE**      **GPW / GAPW**
                                                       if you need         (once M2/M3 land;
                                                       analytic forces     today: BIPOLE or
                                                       or cell opt)        GDF, whichever fits)

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 the cleanest target for the upcoming GPW / GAPW route: the smooth pseudo-density on the FFT grid represents the bonding region efficiently. Once M2 lands, this is where the plane-wave win should first show up at scale.

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, the multi-k Ewald-UKS smearing path is not implemented yet; use the closed-shell route, or accept the Gamma-only UHF/UKS BIPOLE driver without smearing for debugging.

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 only native route that dispatches UHF / UKS today.

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.
ANG2BOHR = 1.0 / 0.529177210903
a = 2.46 * ANG2BOHR
lattice = np.array([
    [a,            0.0, 0.0],
    [a / 2.0,      a * np.sqrt(3.0) / 2.0, 0.0],
    [0.0,          0.0, 30.0],
]).T

graphene = vq.PeriodicSystem(
    dim=2, lattice=lattice,
    unit_cell=[
        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="gdf",
    aux_basis="def2-svp-jk",
    kpoints=(6, 6, 1),
    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¶

Surfaces are where the 1D / 2D vacuum-padding convention earns its keep. A slab is a 2D periodic system: in-plane lattice vectors are the surface mesh, the out-of-plane column is vacuum large enough that the slab images do not see each other (10 to 20 bohr of vacuum on each side is the usual rule).

import numpy as np
import vibeqc as vq

ANG2BOHR = 1.0 / 0.529177210903
a = 4.211 * ANG2BOHR
vacuum = 20.0   # bohr on each side of the slab

# MgO(001) 5-layer slab, primitive 1x1 surface cell.
in_plane = a / np.sqrt(2.0)
lattice = np.array([
    [in_plane, 0.0, 0.0],
    [0.0,      in_plane, 0.0],
    [0.0,      0.0,      a * 2.0 + 2.0 * vacuum],
]).T

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.PeriodicSystem(dim=2, lattice=lattice, unit_cell=layers)
basis = vq.BasisSet(slab.unit_cell_molecule(), "pob-tzvp")

result = vq.run_periodic_job(
    slab, basis,
    method="RKS", functional="pbe",
    jk_method="gdf",                  # 2D: no 3D Madelung problem
    aux_basis="def2-tzvp-jk",         # fallback aux for pob-TZVP
    kpoints=(4, 4, 1),
    output="mgo-001-slab",
)
print(f"MgO(001) 5-layer : {result.energy:+.6f} Ha / cell")

Surface slabs are the use case where GPW / GAPW will eventually dominate once it lands: 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. Until then, GDF in 2D is the production choice.

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 will be the natural fit once it ships.
Metallic 3D bulk (Al, Cu)	BIPOLE + smearing (`smearing_metallic=True`)	UKS metallic smearing not yet implemented.
Open-shell 3D (Li bcc, magnetic oxides)	BIPOLE UHF / UKS	Only native open-shell route today.
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 will dominate once meV-cohesion sign-off lands.
1D wires (H chain, polyacetylene)	GDF (any zeta)	No 3D Madelung problem; BIPOLE direct only.
2D sheets (graphene, h-BN)	GDF (any zeta)	Same as wires.
2D surface slabs	GDF (small to medium), GPW / GAPW (large; future)	Vacuum padding decouples slab images.

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 and the vacuum-padding convention¶

There are two operationally valid ways to make a slab look “non-periodic in the surface normal direction”:

dim = 2. Make the third lattice column purely vacuum and set dim=2. The GDF builder treats only the two in-plane axes as periodic, so the 3D Madelung gauge does not apply. Cheap, clean, and the right physical model. No analytic gradients in v0.8 / v0.9: BIPOLE’s Ewald-J split is 3D-only and falls back to DIRECT_TRUNCATED in 2D; GDF has no periodic gradient driver yet. Use this for single-points and reference energies.
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. Analytic gradients work through the BIPOLE driver, and the existing periodic ASE adapter (vibeqc.ase_periodic) accepts the geometry without any special-casing. The cost is that you pay for one extra Bloch sum in the normal direction even though nothing is happening there, and a residual slab-slab dipole-image interaction shifts the absolute energy by a small constant (cancels in differences such as binding energies).

For force-driven workflows (relaxation, NEB, FD Hessian) today, the dim = 3 + thick-vacuum convention is the operational recommendation. Once GDF periodic gradients land, dim = 2 will become 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.

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", "fcc111", n_layers=4, vacuum=10.0,
                     supercell=(2, 2, 1))
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 analytic gradients¶

Atom relaxation goes through vibeqc.bipole_optimize.relax_atoms, which uses BIPOLE analytic forces 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", "fcc111", n_layers=4, vacuum=10.0,
                     supercell=(2, 2, 1))
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", "fcc111", n_layers=4, vacuum=10.0,
                     supercell=(2, 2, 1))
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,                          # joblib per-image parallel SCFs
)

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", "fcc111", n_layers=4, vacuum=10.0,
                     supercell=(2, 2, 1))
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 (analytic BIPOLE gradients + 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.fun - res_a.fun) * 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 v0.5 placeholder examples/ase_workflows/surface-h2-pt111-singlepoint.py is now obsolete; this skeleton is the canonical native pattern.

7.7 Recommended defaults for surface-reactions runs¶

Knob	Default	Reason
`jk_method`	`"bipole"`	Only route with analytic forces today; closed- and open-shell.
`dim`	3 (with 20-30 bohr vacuum)	ASE adapter + analytic gradients require it.
`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 $	\mathbf{k}
`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 HANDOVER_SURFACE_REACTIONS.md; items marked LANDED below have been crossed off there too.

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.
dim = 2 analytic gradients. Open. GDF has no periodic gradient driver yet; BIPOLE’s Ewald J-split is 3D-only and falls back to DIRECT_TRUNCATED in 2D, which gives correct forces only for well-padded slabs. The dim = 3 + thick-vacuum convention is the operational workaround.
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.
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.
Hessian frozen_indices option for slab TS characterisation. Open. Right now the 6N displacement count is the bottleneck on large slabs.
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, once landed, 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 once it ships.

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 and handover_bipole_v0_8_2026_05_18.md.
GDF multi-k UHF / UKS: open-shell native-GDF drivers are scoped for v0.9.x.
GPW / GAPW: M2 SCF wiring + CP2K parity sign-off, M3 GAPW dual grid, M4 metallic with smearing, M5 ACE-decorated K and on-demand PAW dataset fetcher. Track in design_periodic_gapw.md.
RSGDF and CFMM: enum reserved; implementation not started.

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}\)

	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.