Plane-wave Hartree-J via GPW (experimental)¶

Feature status — experimental (M2-full)

The Gaussian + plane-wave (GPW) Hartree-J route is the v0.10.x acceleration program’s first deliverable (design doc). At M2-full the infrastructure is in place and unit-tested end-to-end, but not yet wired into run_periodic_job — the production runner-level dispatch (where users naturally reach via jk_method="gpw") hard-errors with a pointer to the standalone entry described below.

For routine periodic calculations, prefer:

run_rhf_periodic_gamma_gdf for Γ-only closed-shell hybrid / pure DFT,
run_pbc_bipole_* for multi-k RHF / UHF / RKS / UKS via the CRYSTAL-gauge Ewald J-split.

What GPW exposes today is diagnostic — a closed-shell RHF SCF on simple test cells that demonstrates the M2 chain works end-to-end and pins the cubic-cell Madelung-shift convention. Production GPW (DIIS, level-shift, dynamic damping, runner-level dispatch, all-electron GAPW augmentation, CP2K parity sweep) is M3 work.

What GPW is¶

The Gaussian + plane-wave (GPW) method (Lippert & Hutter 1997, Mol. Phys. 92, 477) replaces the Gaussian-Gaussian four-index Coulomb integral that dominates the periodic Fock build’s J piece with three smaller steps on a smooth real-space grid:

Density collocation. Evaluate ρ(r) = Σ_μν D_μν · χ_μ(r) · χ_ν(r) on a uniform 3D grid spanning the unit cell — pairs the Gaussian basis with a plane-wave-style sampling.
FFT-Poisson. Solve V(r) = ∫ 1/|r − r'| ρ(r') dr' on the same grid via a forward / multiply-by-4π/G² / backward FFT cycle. The G=0 component is pinned to zero (the cell carries an implicit uniform neutralising background).
AO-basis projection. Compute J_μν = ∫ χ_μ(r) χ_ν(r) V(r) dr by quadrature on the same grid.

The cost is O(n_basis² · N_grid) for the collocation + projection plus O(N_grid · log N_grid) for the FFT. For a large enough basis the FFT becomes irrelevant and only the AO work matters; that’s the regime where GPW wins against direct- space periodic Coulomb. CP2K is the canonical implementation the v0.10.x track is reproducing.

GPW assumes pseudopotentials (the density is “smooth” — no sharp 1s cores). For all-electron periodic systems M3 adds a per-atom radial-grid augmentation (the “GA” in GAPW) that restores the hard cores; that work isn’t yet in vibe-qc.

What today’s implementation does¶

The M2-full stack ships these surfaces, all on the public vibeqc.* namespace:

vibeqc.PlaneWaveGrid — frozen dataclass holding a uniform 3D grid (nx, ny, nz, lattice, cutoff). Lazy fractional + Cartesian coordinate generators; reciprocal- lattice vectors with FFT wrap-around convention.
vibeqc.make_grid — construct a grid sized from a plane-wave cutoff (Hartree). Rounds each extent up to a {2ᵃ · 3ᵇ · 5ᶜ} value so FFTW3 hits its fast paths.
vibeqc.collocate_density_on_grid — Gaussian density on the grid via the existing C++ evaluate_ao accessor. Integrates to tr(D · S) (the total number of electrons for a closed-shell density) to finite-grid quadrature precision.
vibeqc.collocate_point_charges_on_grid — point charges as normalised Gaussians of width σ on the grid, summing periodic images out to image_shells neighbours. Used by M1e’s Madelung-constant Hartree test.
vibeqc.GpwJBuilder — caches the AO-on-grid values across SCF iterations; build_J(D) returns the Hartree-J matrix.
vibeqc.evaluate_gpw_energy — single-point periodic-energy breakdown at a given density. Returns vibeqc.GpwEnergyBreakdown with the per-term decomposition (e_kinetic, e_nuclear_attraction, e_hartree, e_hf_exchange, e_nuclear_repulsion, e_total).
vibeqc.run_periodic_rhf_gpw — minimal iterative closed-shell RHF SCF using GPW for J and the existing molecular kernel for K. Returns vibeqc.GpwScfResult with energy, density, MO coefficients + energies, and the per-term breakdown.

Minimal worked example — He STO-3G¶

import warnings
import numpy as np
import vibeqc as vq
from vibeqc import _vibeqc_core as core

# He atom at the centre of a 16-bohr cubic cell.
L = 16.0
system = core.PeriodicSystem()
system.dim = 3
system.lattice = np.eye(3) * L
system.unit_cell = [core.Atom(2, [L / 2, L / 2, L / 2])]
basis = vq.BasisSet(
    vq.Molecule(list(system.unit_cell), 0, 1),
    "sto-3g",
)

# Suppress the experimental warning during opt-in use.
warnings.simplefilter("ignore", vq.GAPWExperimentalWarning)

result = vq.run_periodic_rhf_gpw(
    system, basis,
    cutoff_ha=300.0,        # ~600 Ry; CP2K's production default
    max_iter=50,
    conv_tol_energy=1e-9,
    conv_tol_density=1e-7,
)

print(f"converged: {result.converged}")
print(f"n_iter: {result.n_iter}")
print(f"E_total: {result.energy:.10f} Ha")
print()
print("breakdown:")
for field in ("e_kinetic", "e_nuclear_attraction",
              "e_hartree", "e_hf_exchange",
              "e_nuclear_repulsion"):
    print(f"  {field}: {getattr(result.breakdown, field):+.6f}")

The standalone diagnostic evaluate_gpw_energy(system, basis, D, grid=grid) does the same per-term decomposition at a density D you supply — useful when you want to evaluate the GPW energy at the converged molecular density for comparison.

The Madelung-shift convention¶

GPW’s FFT-Poisson kernel runs on a periodically-repeated cell. Even a single neutral atom in a finite box (He, H₂, …) is treated as one cell of an infinite periodic lattice — the electrostatic energy you get out is the lattice’s energy, not the isolated-molecule energy. The difference is the cubic- cell Madelung shift: each point charge interacts with its own periodic images plus the uniform neutralising background. The shift per unit charge in a cubic cell of side L is

ΔE_Madelung = q² · ξ_NaCl / (2 · L)

with ξ_NaCl = 2.837297… the Madelung constant for the simple cubic lattice.

At the M2 J-only level, this shift falls between the periodic GPW Hartree and the molecular Hartree as a clean single-term correction. The vibeqc.evaluate_gpw_energy result’s e_hartree field carries the electronic-side Madelung shift on a charged cell.

At the M3a SCF level, V_ne is built in the same Ewald gauge as the Hartree-J (compute_nuclear_lattice_dispatch with CoulombMethod.EWALD_3D — G = 0 reciprocal mode dropped, v_bg jellium shift applied), so the three pieces all sit in the same gauge. On a neutral cell the cross-term cancels both the electronic and nuclear self-image shifts exactly:

E_periodic_M3a ≈ E_molecular  (vacuum-padded neutral cell;
                              < 1 mHa residual from finite
                              cell + finite grid)

For He STO-3G in a 16-bohr cube, periodic E_HF and molecular E_HF now agree to better than 2·10⁻⁵ Ha (verified by tests/test_periodic_gapw_j.py:: test_rhf_gpw_he_recovers_molecular_in_vacuum_padded_cell).

Before the M3a lift, V_ne was at the molecular limit while the Hartree-J and nuclear-repulsion both lived in the Ewald gauge, leaving a known “2× per-charge Madelung shift” in the SCF total energy (0.71 Ha for He STO-3G at L = 16). That is gone: any residual cell-size dependence is now the physical finite-cell / finite-grid contribution, not a gauge artefact.

At the M2 J-only level (evaluate_gpw_energy(...).e_hartree in isolation), the periodic-vs-molecular Hartree difference still carries the electronic self-image shift q² · ξ_NaCl / (2L). That is the FFT-Poisson convention of the J build itself and is the M2 internal-validation anchor that pairs with M1e’s neutral-cell sub-µHa NaCl test. The e_hartree field thus continues to read out a Madelung-shifted Hartree energy on charged cells; the total-energy cancellation kicks in only when V_ne and E_nn are summed in.

What the runner reports if you ask for it via `run_periodic_job`¶

Until M3 lands the GpwScfResult → PeriodicSCFResult adapter that lets the runner-level dispatcher serialise GPW SCF output into the standard .out / .system / .bibtex artefacts, the runner refuses:

>>> vq.run_periodic_job(system, basis, jk_method="gpw", ...)
NotImplementedError: PeriodicJKMethod.GPW is not yet plumbed
through `run_periodic_job`. Use
`vibeqc.run_periodic_rhf_gpw(system, basis, *, grid=None,
cutoff_ha=300)` instead — the standalone GPW SCF entry.
Closed-shell RHF only at M2-full; DFT / open-shell are M3
work. See `docs/design_periodic_gapw.md` § GPW.

That’s intentional — the pointer at run_periodic_rhf_gpw guides users at the working entry without silently dropping through to a half-wired code path (CLAUDE.md § 7).

What’s coming¶

The next gapw-chat session picks up the M3 chain:

C++ JKBuilder trampoline for GpwJBuilder so the production C++ SCF (run_rhf_scf_with_jk — DIIS, dynamic damping, level-shift) can host the GPW J kernel.
PeriodicJKMethod.GPW runner wiring — adapt GpwScfResult into the periodic-SCF result shape the run_periodic_job plumbing expects.
GAPW augmentation (M3b) — per-atom radial-grid all- electron correction on top of the smooth-grid GPW J. Closes the all-electron parity gap against CP2K. (The convention half of M3 — lifting V_ne to the periodic Ewald gauge — landed at M3a; the augmentation piece is what’s still pending.)
CP2K parity sweep (M3c) via the M1f examples/regression/core/runner_cp2k.py subprocess runner.
Metallic GAPW + smearing (composes with the v0.10.x density-mixing flagship).
ACE-decorated K for hybrids; PAW pseudopotentials via an on-demand paw_fetcher.py (the design-doc decision to pull the fetcher forward into M5).

The full milestone map is in docs/design_periodic_gapw.md; the per-session drop-box at .release-status/v0.10.x/gapw.md (gitignored) is where the in-flight work is tracked.