Plane-wave Hartree-J via GPW / GAPW¶

Feature status — production + experimental gates

The Gaussian + plane-wave (GPW/GAPW) route serves two tiers:

Production (ungated) — ship as of v0.12:

RHF / RKS / UHF / UKS at Γ via run_periodic_job(jk_method="gpw") or the standalone run_periodic_*_gpw entries.
Multi-k pure-DFT via run_periodic_rks_gpw_multi_k.
DIIS, density damping, analytic forces, finite-difference Hessians, DOS / band-path helpers, .npz restart I/O, and the VibeqcGPW ASE calculator.
Both V_ne conventions ("ewald" default / "smeared_erfc").
Reproduces the molecular RHF limit to sub-µHa (0.9 µHa vs pyscf RHF on H₂/STO-3G, run out of process).

Experimental (opt-in, gated by GAPWExperimentalWarning):

GAPW augmentation — all-electron per-atom radial × Lebedev Hartree/XC correction (\`jk_method="gapw"\```, periodic_gapw_augment``). Not yet at CP2K parity for tight ionic crystals; the atomic-grid infrastructure ships but the augmentation SCF path + PAW atomic-data fetch are open roadmap items.
Range-separated hybrids — HSE06 / ωB97X / CAM-B3LYP on the GPW/GAPW path (periodic_gapw_range_sep).
Orbital Transformation (OT) solver — direct energy minimisation avoiding diagonalisation (periodic_gapw_ot).
Multi-k GAPW — all-electron multi-k with augmentation correction.
MPI z-slab grid overlay — distributed-memory Hartree-J build via GpwJBuilder(mpi_aware=True) (z-slab domain decomposition over mpi4py). The multi-rank build was broken for every world size > 1 through v0.12.0 (the v0.12 post-release audit caught it; fixed and pinned by a simulated-rank parity suite, tests/test_mpi_gpw_parity.py). Remains experimental until validated on a real mpirun multi-rank run. The serial GPW route is unaffected — world size 1 always degenerated to the correct serial build.

For tight all-electron crystals today, use the GDF or BIPOLE routes.

Known limitation — boundary-atom density quadrature (∫ρ ≠ N):

The periodic-AO image sum in collocate_density_on_grid sums nearest-neighbour images (image_radius=1, 27 cells) to fold escaping Gaussian tails back into the cell. This is exact for compact (valence) bases in cells with side ≥ ~6 bohr, but diffuse bases in very tight cells may still leak — the quadrature error shows as a small under-integration of the density (∫ρ < N). GPW cubes, ELF maps, and post-SCF grid-based properties derived from the collocated density are affected. The SCF energy and wavefunction are not affected (they use the rigorous lattice-summed AO integrals). Users with diffuse bases in tight cells can increase image_radius via the _AO_IMAGE_RADIUS module constant (not yet a public option).

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, the GAPW route adds a per-atom radial-grid augmentation (the “GA” in GAPW) that restores the hard cores — this is wired and functional (jk_method="gapw") but remains experimental (guarded by GAPWExperimentalWarning).

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.

V_ne convention switch (M3a vs. M3b)¶

Both periodic V_ne conventions ship and are physically equivalent on neutral cells. The default "ewald" is the M3a path that uses compute_nuclear_lattice_ewald (G = 0 dropped, v_bg jellium shift applied) — production-tested, no extra parameter. The opt-in "smeared_erfc" is the CP2K-style decomposition

V_ne = V_ne_long  +  V_ne_short  +  c_α · S
     = (FFT-Poisson on Σ_a Z_a · (α/√π)³ exp(-α²|r-R_a|²))
     + ⟨χ_μ | -Σ_a Z_a · erfc(α · r_a) / r_a | χ_ν⟩
     + π · (Σ_a Z_a) / (α² · V_cell) · S    # PW-DFT α-correction

routing every Hartree-class quantity through the same Poisson solver — the convention the GAPW augmentation work plugs into.

# M3a default (Ewald):
r_ewald = vq.run_periodic_rhf_gpw(system, basis, cutoff_ha=300.0)

# M3b CP2K-style erfc/erf split:
r_smear = vq.run_periodic_rhf_gpw(
    system, basis, cutoff_ha=300.0,
    v_ne_convention="smeared_erfc",
    smearing_alpha=2.0,             # bohr⁻¹; None ⇒ auto from grid
)

# Both converge to the same SCF total within < 1 mHa on neutral
# cells (parity tests pin this at < 50 µHa on He STO-3G across
# α ∈ {1.0, 1.5, 2.0} and at < 1 µHa on H2 STO-3G at α = 2.0).

Atomic radial × Lebedev grids¶

The per-atom quadrature grids for the GAPW hard/soft Hartree decomposition are available in vibeqc.periodic_gapw_atomic_grid. These grids are now wired into the GAPW SCF — the augmentation uses radial Poisson solves + Lebedev integration inside each atomic sphere, with a multipole-compensating ρ₀ on the smooth FFT grid.

What the module provides:

AtomicRadialGrid — a frozen dataclass that owns a 1D Mura-Knowles radial mesh (r, w_r, with the r² Jacobian pre-folded into w_r) and a Lebedev-Laikov angular mesh (angular_xyz, w_a, with Σ_l w_l = 4π). Helpers: cartesian_points() returns the product grid as (n_radial, n_angular, 3); combined_weights() returns the (n_radial, n_angular) w_k · w_l; integrate(values) evaluates Σ_kl values[k,l] · w_k · w_l directly.
AtomicRadialGrid.build(centre_bohr, n_radial=..., alpha=..., lebedev_order=...) and AtomicRadialGrid.from_element( centre_bohr, Z, ...) are the two constructors. from_element pulls α from default_alpha_for_element(Z) and emits a GAPWExperimentalWarning (suppressible via quiet=True).
default_alpha_for_element(Z) — per-element radial scale picked to hug the nucleus and the diffuse tail equally; the CP2K grid_atom family’s choice.
lebedev_supported_orders() — the Lebedev orders SciPy exposes (scipy.integrate.lebedev_rule); the grid constructor validates against this list.

import numpy as np
import vibeqc as vq

grid_O = vq.periodic_gapw_atomic_grid.AtomicRadialGrid.from_element(
    centre_bohr=np.array([0.0, 0.0, 0.0]),
    Z=8,
    n_radial=80,
    lebedev_order=29,
    quiet=True,
)
# grid_O.cartesian_points() has shape (n_radial, n_angular, 3)
# grid_O.combined_weights() has shape (n_radial, n_angular)
# Σ_k w_r[k] = ∫_0^∞ r² · 1 dr  truncated to the mesh range
# Σ_l w_a[l] = 4π

The radial rule is Mura-Knowles (1996) — r(x) = -α · ln(1 - x³) on uniform x ∈ (0, 1); numerically equivalent to the log-transformed Gauss-Chebyshev rule the CP2K source uses for the smooth integrands GAPW sees. The implementation is paraphrased from the published paper (no GPL source reproduced).

Note

The grids are unpartitioned — no Becke weighting. GAPW only needs them inside an augmentation sphere around each nucleus, so the molecular Becke grid (vibeqc._vibeqc_core.build_grid) is the wrong tool here. The augmentation work in M3b-aug-B applies the soft-cut at the augmentation radius on top of these raw atomic shells.

See docs/design_periodic_gapw.md § Q2 / D1 / D2 for the audit reference against CP2K’s grid_atom.F.

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

As of the M3b dispatch wiring, the runner serialises GPW SCF output through the standard .out / .system / .molden artefacts. Usage:

result = vq.run_periodic_job(
    system, basis,
    method="RHF",
    jk_method="gpw",          # M3b: now reachable
    output="my_gpw_calc",
)
print(f"E_total = {result.energy:.10f} Ha")

The runner wraps the standalone GpwScfResult with the :class:vibeqc.periodic_gapw_runner_adapter.GpwRunnerResult adapter so the downstream output-rendering / manifest plumbing sees the same shape the GDF / BIPOLE routes produce (energy / e_electronic / e_nuclear / e_coulomb / e_hf_exchange / n_iter / converged / mo_coeffs / mo_energies / density / fock / overlap / occupations / scf_trace).

Constraints (current): RHF / RKS / UHF / UKS. Γ-only for RHF, UHF, UKS; multi-k RKS (pure DFT) available via kpoints. Neutral cells (RKS / UKS require a functional=). jk_method="auto" continues to pick GDF (closed-shell) / BIPOLE (open-shell); GPW/GAPW are opt-in experimental routes.

Tutorials¶

The walk-throughs below mirror GPAW’s tutorial style: a short motivation, a single self-contained Python block, and a pointer at the API reference earlier on this page so you can jump back to the underlying knob when you want to tweak it. Every example uses H₂ in a vacuum-padded cube — small enough to converge in seconds on a laptop, big enough that the periodic answer matches the molecular reference. Drop each block into a .py file next to your checkout and run with the project’s .venv/bin/python.

Tutorial 1: Your first GPW calculation¶

The fastest way to convince yourself the plane-wave-on-Gaussian machinery is wired up is to run a single GPW SCF and compare its total energy against a textbook number. H₂ at 1.4 bohr in a 15-bohr cube is a good first system: it is closed-shell, the vacuum-padded cell carries no Madelung shift (see The Madelung-shift convention), and the molecular RHF/STO-3G energy is the well-known −1.1167 Ha.

What you’re seeing in the block below: we build a 3D PeriodicSystem with a cubic lattice, attach an STO-3G basis to its two H atoms, and call vibeqc.run_periodic_rhf_gpw. The function returns a GpwScfResult carrying the converged energy, iteration count, and the per-term breakdown documented in the API section above. The e_hartree and e_nuclear_attraction fields both live in the Ewald gauge at this milestone, so the total is directly comparable against the molecular RHF number.

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

L = 15.0
d = 1.4  # H–H bond, bohr
system = core.PeriodicSystem()
system.dim = 3
system.lattice = np.eye(3) * L
system.unit_cell = [
    core.Atom(1, [L / 2 - d / 2, L / 2, L / 2]),
    core.Atom(1, [L / 2 + d / 2, L / 2, L / 2]),
]
basis = vq.BasisSet(vq.Molecule(list(system.unit_cell), 0, 1), "sto-3g")

warnings.simplefilter("ignore", vq.GAPWExperimentalWarning)

result = vq.run_periodic_rhf_gpw(
    system, basis,
    cutoff_ha=300.0,
    max_iter=50,
    conv_tol_energy=1e-9,
)
print(f"E_HF (periodic GPW) = {result.energy:.6f} Ha")
print(f"molecular reference  ≈ -1.1167 Ha (H2 / STO-3G)")

Next: Tutorial 2: Periodic DFT on a hydrogen molecule.

Tutorial 2: Periodic DFT on a hydrogen molecule¶

GPW’s J kernel is method-agnostic: once the Hartree-J is in the Ewald gauge, swapping the molecular RHF on top for a pure-DFT RKS only changes the K piece. The M3d milestone wired libxc into the GPW SCF, so you reach RKS-LDA on the periodic side with a single keyword swap.

Why does the periodic LDA number land on the molecular LDA number? Because the H₂ cell is vacuum-padded and neutral. With 12+ bohr of empty space between periodic images, the nearest H₂ replica sits beyond the practical range of the Hartree kernel; on a neutral cell the Ewald V_ne / E_Madelung cross-term cancels the electronic self-image shift (see The Madelung-shift convention); and LDA’s exchange-correlation is fully local so it can’t pick up any inter-image contribution at all. The takeaway: if you want a molecular DFT answer out of the periodic stack, pad the cell.

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

L, d = 15.0, 1.4
system = core.PeriodicSystem()
system.dim = 3
system.lattice = np.eye(3) * L
system.unit_cell = [
    core.Atom(1, [L / 2 - d / 2, L / 2, L / 2]),
    core.Atom(1, [L / 2 + d / 2, L / 2, L / 2]),
]
basis = vq.BasisSet(vq.Molecule(list(system.unit_cell), 0, 1), "sto-3g")

warnings.simplefilter("ignore", vq.GAPWExperimentalWarning)

result = vq.run_periodic_rks_gpw(
    system, basis,
    functional="lda",
    cutoff_ha=300.0,
    conv_tol_energy=1e-8,
)
print(f"E_LDA (periodic GPW) = {result.energy:.6f} Ha")
print("(matches molecular RKS-LDA on a vacuum-padded cell)")

Next: Tutorial 3: Multi-k Brillouin-zone sampling.

Tutorial 3: Multi-k Brillouin-zone sampling¶

For metallic and dispersive systems the SCF total energy depends on a Brillouin-zone integral over the occupied bands, and Γ-only sampling is not enough. The M3e milestone added vibeqc.run_periodic_rks_gpw_multi_k — a closed-shell, pure-DFT route that Bloch-sums the core operators (T, S, V_ne) on a Monkhorst-Pack mesh while sharing a single FFT-grid density (and therefore a single Hartree-J and V_xc) across all k-points.

On the H₂ / STO-3G vacuum-padded cube the answer is independent of kmesh — the bands are flat (no dispersion) because each H₂ image is electronically isolated, so every k-point sees the same Fock spectrum. That’s the diagnostic: if your kmesh sweep shows < µHa drift, you’re correctly in the molecular limit and Γ is sufficient. The moment you shrink the cell or pack atoms densely enough that bands disperse, the kmesh result will diverge from Γ and you’ll need to converge it (start at (2,2,2), double until ΔE < your threshold).

The block below runs the same H₂ cell at a (2,2,2) Monkhorst- Pack mesh and prints the energy alongside Γ-only for comparison.

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

L, d = 15.0, 1.4
system = core.PeriodicSystem()
system.dim = 3
system.lattice = np.eye(3) * L
system.unit_cell = [
    core.Atom(1, [L / 2 - d / 2, L / 2, L / 2]),
    core.Atom(1, [L / 2 + d / 2, L / 2, L / 2]),
]
basis = vq.BasisSet(vq.Molecule(list(system.unit_cell), 0, 1), "sto-3g")

warnings.simplefilter("ignore", vq.GAPWExperimentalWarning)

r_gamma = vq.run_periodic_rks_gpw(
    system, basis, functional="lda", cutoff_ha=300.0,
)
kmesh = vq.monkhorst_pack(system, [2, 2, 2])
r_mk = vq.run_periodic_rks_gpw_multi_k(
    system, basis, kmesh,
    functional="lda",
    cutoff_ha=300.0,
)
print(f"E_LDA (Γ-only)    = {r_gamma.energy:.8f} Ha")
print(f"E_LDA (2x2x2 MP)  = {r_mk.energy:.8f} Ha")
print(f"|ΔE|              = {abs(r_gamma.energy - r_mk.energy):.2e} Ha")

Next: Tutorial 4: Choosing the V_ne convention.

Tutorial 4: Choosing the V_ne convention (Ewald vs CP2K-style)¶

GPW’s nuclear-attraction operator V_ne is gauge-invariant on a neutral cell — two conventions ship and both give the same SCF total to better than 1 mHa. The default v_ne_convention="ewald" builds V_ne in the same Ewald gauge as the Hartree-J directly. The opt-in v_ne_convention="smeared_erfc" is the CP2K-style decomposition that splits V_ne into a smooth long-range piece (routed through the FFT-Poisson solver as a Gaussian-smeared nuclear charge density) plus a short-range erfc/r correction in AO basis (see V_ne convention switch (M3a vs. M3b) on this page).

For routine vacuum-padded RHF / RKS jobs the default Ewald convention is the right pick: one less parameter, no smearing_alpha, and no spurious Madelung shift to track. The smeared-erfc convention earns its keep when you head toward GAPW augmentation: it routes every Hartree-class quantity through the same FFT-Poisson solver, which is the convention the per-atom radial-grid hard/soft cancellation will exploit. If you want to start instrumenting against the augmentation work, develop with smeared_erfc; if you just want an energy out, stick with Ewald.

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

L, d = 15.0, 1.4
system = core.PeriodicSystem()
system.dim = 3
system.lattice = np.eye(3) * L
system.unit_cell = [
    core.Atom(1, [L / 2 - d / 2, L / 2, L / 2]),
    core.Atom(1, [L / 2 + d / 2, L / 2, L / 2]),
]
basis = vq.BasisSet(vq.Molecule(list(system.unit_cell), 0, 1), "sto-3g")

warnings.simplefilter("ignore", vq.GAPWExperimentalWarning)

r_ewald = vq.run_periodic_rhf_gpw(
    system, basis, cutoff_ha=300.0, v_ne_convention="ewald",
)
r_smear = vq.run_periodic_rhf_gpw(
    system, basis, cutoff_ha=300.0,
    v_ne_convention="smeared_erfc", smearing_alpha=2.0,
)
print(f"E (Ewald V_ne)        = {r_ewald.energy:.8f} Ha")
print(f"E (smeared-erfc V_ne) = {r_smear.energy:.8f} Ha")
print(f"|ΔE|                  = {abs(r_ewald.energy - r_smear.energy):.2e} Ha")

Next: Tutorial 5: Visualising the density.

Tutorial 5: Visualising the density¶

A periodic GPW job collocates the density on the FFT grid as its very first SCF step, so dropping it to disk is essentially free. Driving the calculation through the runner — via run_periodic_job with write_density=True — emits the standard .out / .system / .molden triple alongside an .xsf volumetric file. Open the .xsf in VESTA or any volumetric- file viewer (vibe-view, VMD, Avogadro) to inspect the bonding density.

The block below runs the same H₂ vacuum-padded cube under RKS-LDA and asks for the density dump. The runner picks the GPW J/K backend (jk_method="gpw") and routes the rest of the periodic-SCF plumbing through the normal run_periodic_job flow, so the output layout matches every other periodic job’s. Future milestones may add a dedicated CUBE writer (e.g. a write_cube_density_periodic helper) and a QVF volume.density section for the GAPW runner; for now the XSF file is the one to open.

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

L, d = 15.0, 1.4
system = core.PeriodicSystem()
system.dim = 3
system.lattice = np.eye(3) * L
system.unit_cell = [
    core.Atom(1, [L / 2 - d / 2, L / 2, L / 2]),
    core.Atom(1, [L / 2 + d / 2, L / 2, L / 2]),
]
basis = vq.BasisSet(vq.Molecule(list(system.unit_cell), 0, 1), "sto-3g")

warnings.simplefilter("ignore", vq.GAPWExperimentalWarning)

result = vq.run_periodic_job(
    system, basis,
    method="RKS",
    functional="lda",
    jk_method="gpw",
    write_density=True,
    output="h2_gpw_density",
)
print(f"E_LDA = {result.energy:.8f} Ha")
print("open h2_gpw_density.xsf in VESTA or vibe-view")

Next: Tutorial 6: Saving and resuming a GPW calculation.

Tutorial 6: Saving and resuming a GPW calculation¶

A converged GPW SCF carries everything a downstream post-SCF analysis needs — density, MO coefficients + energies, the per-term breakdown, and the lattice / basis / atoms that built it. The v0.12-prep restart format ships three helpers that round-trip the lot through a single compressed .npz archive: vibeqc.save_gpw_result(path, result, basis, system) writes the archive, vibeqc.load_gpw_result(path) reads it back into a plain dict, and vibeqc.describe_gpw_result(path) returns a short human-readable summary string for quick inspection.

This is useful any time you’d like to avoid re-running the SCF — running a band path against a converged density, running a DOS sweep over different smearing widths, sharing a result with a collaborator, or just resuming an analysis after a long SCF without having to wait the SCF out again. Both GpwScfResult (Γ-only) and GpwMultiKScfResult (multi-k) serialise through the same pair of calls.

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

L, d = 15.0, 1.4
system = core.PeriodicSystem()
system.dim = 3
system.lattice = np.eye(3) * L
system.unit_cell = [
    core.Atom(1, [L / 2 - d / 2, L / 2, L / 2]),
    core.Atom(1, [L / 2 + d / 2, L / 2, L / 2]),
]
basis = vq.BasisSet(vq.Molecule(list(system.unit_cell), 0, 1), "sto-3g")

warnings.simplefilter("ignore", vq.GAPWExperimentalWarning)

result = vq.run_periodic_rks_gpw(
    system, basis, functional="lda", cutoff_ha=300.0,
)
written = vq.save_gpw_result("h2_gpw.npz", result, basis, system)
print(f"wrote {written} ({written.stat().st_size} bytes)")
print(vq.describe_gpw_result(written))

# Resume later in a fresh process:
data = vq.load_gpw_result(written)
print(f"resumed E = {data['energy']:.8f} Ha")
print(f"resumed density shape = {data['density'].shape}")

Next: Tutorial 7: Computing a density of states.

Tutorial 7: Computing a density of states¶

The Γ-only spectrum is one diagonalisation; a multi-k spectrum is a stack of diagonalisations weighted by the Monkhorst-Pack weights. Either way, what you want next is usually a Gaussian-broadened density of states — and a quick HOMO / LUMO read-out to anchor the energy axis. The v0.12-prep post-SCF helpers in vibeqc.periodic_gapw_postscf (re-exported on the top-level vibeqc.* namespace) cover both: gaussian_dos(mo_energies, ...) is the raw broadening kernel, compute_dos_from_result(result, sigma=..., ...) pulls eigenvalues straight off a GpwScfResult or GpwMultiKScfResult and weights k-points for you, and compute_homo_lumo(result) returns the (e_homo, e_lumo, gap) triple from the converged spectrum.

If you also need a band path between high-symmetry k-points, band_path_eigenvalues(result, basis, system, k_path, ...) diagonalises F(k) C = ε S(k) C at each k on a user-supplied path — useful for plotting band structures once you have a converged density.

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

L, d = 15.0, 1.4
system = core.PeriodicSystem()
system.dim = 3
system.lattice = np.eye(3) * L
system.unit_cell = [
    core.Atom(1, [L / 2 - d / 2, L / 2, L / 2]),
    core.Atom(1, [L / 2 + d / 2, L / 2, L / 2]),
]
basis = vq.BasisSet(vq.Molecule(list(system.unit_cell), 0, 1), "sto-3g")

warnings.simplefilter("ignore", vq.GAPWExperimentalWarning)

result = vq.run_periodic_rks_gpw(
    system, basis, functional="lda", cutoff_ha=300.0,
)
e_homo, e_lumo, gap = vq.compute_homo_lumo(result, system=system)
print(f"HOMO = {e_homo:+.4f} Ha   LUMO = {e_lumo:+.4f} Ha   gap = {gap:.4f} Ha")

e_grid, dos = vq.compute_dos_from_result(
    result, system=system, sigma=0.01, e_range=(e_homo - 0.5, e_lumo + 0.5),
)
print(f"DOS array shape = {dos.shape};  max DOS = {dos.max():.2f} states/Ha")

Next: Tutorial 8: Geometry optimisation with ASE.

Tutorial 8: Geometry optimisation with ASE¶

The v0.12 R1 cycle added ASE-compatible numerical forces to the VibeqcGPW Calculator — central-difference around each Cartesian coordinate, exposed through the standard atoms.get_forces() interface. That is the only piece ASE’s geometry-optimisation drivers need: drop the calculator on an ase.Atoms object and any optimiser in ase.optimize will relax it.

H₂ at HF/STO-3G has a well-known equilibrium bond length near 1.346 bohr (~0.712 Å). The block below starts from a stretched 1.6-bohr bond and runs BFGS for a handful of steps; the bond contracts toward the optimum and the final length lands inside the noise floor of the numerical-force step size. The same recipe drives any closed-shell GPW geometry: pick a basis + cutoff that converge your system’s total energy, attach the calculator, and hand the Atoms off to BFGS / LBFGS / FIRE.

import warnings
import numpy as np
from ase import Atoms
from ase.units import Bohr
from ase.optimize import BFGS
import vibeqc as vq

warnings.simplefilter("ignore", vq.GAPWExperimentalWarning)

L_bohr = 12.0
L_ang = L_bohr * Bohr
r0_bohr = 1.6  # stretched start
p1 = np.array([L_bohr / 2 - r0_bohr / 2, L_bohr / 2, L_bohr / 2]) * Bohr
p2 = np.array([L_bohr / 2 + r0_bohr / 2, L_bohr / 2, L_bohr / 2]) * Bohr
atoms = Atoms("H2", positions=[p1, p2], pbc=True)
atoms.set_cell(np.eye(3) * L_ang)
atoms.calc = vq.VibeqcGPW(basis="sto-3g", cutoff_ha=200.0)

r_initial = atoms.get_distance(0, 1)
BFGS(atoms, logfile=None).run(fmax=0.05, steps=10)
r_final = atoms.get_distance(0, 1)
print(f"initial bond = {r_initial / Bohr:.3f} bohr")
print(f"final   bond = {r_final / Bohr:.3f} bohr (HF/STO-3G opt ≈ 1.346 bohr)")

Next: Tutorial 9: Metallic systems with smearing.

Tutorial 9: Metallic systems with smearing¶

Metallic and small-gap systems oscillate under sharp Aufbau occupations because the Fermi level passes through a band of states the SCF can’t decide between. The fix is Fermi-Dirac smearing: spread the occupations smoothly across an energy window of width kT, converge the resulting variationally correct free energy F = E − T·S, then extrapolate to T=0 if needed. The v0.12 R1 cycle wires this onto the multi-k GPW entry run_periodic_rks_gpw_multi_k via the smearing_temperature kwarg (in Hartree).

The block below runs the H₂ vacuum-padded cube at a tiny finite kT just to show what the surface looks like — on a gapped molecular reference the entropy stays small and F sits just below E by T·S. On a real metal you’d typically pick kT in the 0.001–0.01 Ha range (≈ 300–3000 K), monitor the smearing_entropy term to make sure it isn’t ballooning, and report free_energy (not energy) as the converged total.

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

L, d = 15.0, 1.4
system = core.PeriodicSystem()
system.dim = 3
system.lattice = np.eye(3) * L
system.unit_cell = [
    core.Atom(1, [L / 2 - d / 2, L / 2, L / 2]),
    core.Atom(1, [L / 2 + d / 2, L / 2, L / 2]),
]
basis = vq.BasisSet(vq.Molecule(list(system.unit_cell), 0, 1), "sto-3g")

warnings.simplefilter("ignore", vq.GAPWExperimentalWarning)

kmesh = vq.monkhorst_pack(system, [2, 2, 2])
result = vq.run_periodic_rks_gpw_multi_k(
    system, basis, kmesh,
    functional="lda", cutoff_ha=300.0,
    smearing_temperature=0.005,
)
print(f"E (internal)   = {result.energy:.8f} Ha")
print(f"F = E - T·S    = {result.free_energy:.8f} Ha")
print(f"smearing T     = {result.smearing_temperature:.4f} Ha")
print(f"entropy S      = {result.smearing_entropy:.4e}")

Next: Tutorial 10: Periodic UHF on a hydrogen atom.

Tutorial 10: Periodic UHF on a hydrogen atom¶

The v0.12 R2 cycle wires the open-shell GPW route: run_periodic_uhf_gpw (pure HF) and run_periodic_uks_gpw (DFT with a libxc functional) both take per-spin occupations n_alpha / n_beta, run a Γ-only joint α + β DIIS-accelerated SCF, and return a GpwUhfScfResult / GpwUksScfResult with per-spin densities, MO coefficients, and MO energies. They reduce bit-for-bit to the closed-shell run_periodic_rhf_gpw / run_periodic_rks_gpw results when n_alpha == n_beta — that parity is the regression test pin.

The block below runs a single H atom in a vacuum-padded box at (n_alpha, n_beta) = (1, 0) — the spin-polarised reference for the 1s electron. The reported energy is the standard UHF total; spin_squared (<S²>) tracks the spin contamination and lands near 0.75 (the doublet exact value) when the open-shell SCF converges cleanly.

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

L = 12.0
system = core.PeriodicSystem()
system.dim = 3
system.lattice = np.eye(3) * L
system.unit_cell = [core.Atom(1, [L / 2, L / 2, L / 2])]
# Single H atom = 1 electron → multiplicity 2 (doublet); mult 1 would raise
# "n_electrons and multiplicity are inconsistent".
basis = vq.BasisSet(vq.Molecule(list(system.unit_cell), 0, 2), "sto-3g")

warnings.simplefilter("ignore", vq.GAPWExperimentalWarning)

result = vq.run_periodic_uhf_gpw(
    system, basis,
    n_alpha=1, n_beta=0,
    cutoff_ha=200.0,
)
print(f"E (UHF/STO-3G) = {result.energy:.8f} Ha")
print(f"n_alpha        = {result.n_alpha}")
print(f"n_beta         = {result.n_beta}")
print(f"converged      = {result.converged} ({result.n_iter} iter)")

Next: Tutorial 11: Geometry optimisation with analytic forces.

Tutorial 11: Geometry optimisation with analytic forces¶

The v0.12 R2 cycle replaces Tutorial 8’s central-difference numerical-force loop with an analytic-where-possible gradient driver (vibeqc.periodic_gapw_gradient.compute_gradient_gpw): the Pulay overlap-Lagrangian uses the C++ periodic gradient primitive, and the combined Hellmann-Feynman piece (T + V_ne + Hartree-on-grid + XC-on-grid + E_nn at fixed density) is a central-difference on the energy evaluator only — no extra SCF iteration per displacement. Cost per atom per axis: two fixed-density energy evaluations. On STO-3G H₂ this is ~6× faster than the 6N-SCF path from Tutorial 8 and produces forces that agree to <1% of the FD step size.

VibeqcGPW.calculate(properties=("energy", "forces")) now defaults to the analytic path on the Γ-only route, so ASE’s optimisers pick it up automatically. The block below repeats the Tutorial 8 H₂ BFGS relaxation with the analytic gradient selected by default; pass use_numerical_forces=True to fall back to the legacy 6N-SCF path for cross-checks.

import warnings
import numpy as np
from ase import Atoms
from ase.units import Bohr
from ase.optimize import BFGS
import vibeqc as vq

warnings.simplefilter("ignore", vq.GAPWExperimentalWarning)

L_bohr = 12.0
L_ang = L_bohr * Bohr
r0_bohr = 1.6  # stretched start
p1 = np.array([L_bohr / 2 - r0_bohr / 2, L_bohr / 2, L_bohr / 2]) * Bohr
p2 = np.array([L_bohr / 2 + r0_bohr / 2, L_bohr / 2, L_bohr / 2]) * Bohr
atoms = Atoms("H2", positions=[p1, p2], pbc=True)
atoms.set_cell(np.eye(3) * L_ang)
atoms.calc = vq.VibeqcGPW(basis="sto-3g", cutoff_ha=200.0)

r_initial = atoms.get_distance(0, 1)
BFGS(atoms, logfile=None).run(fmax=0.05, steps=10)
r_final = atoms.get_distance(0, 1)
print(f"initial bond = {r_initial / Bohr:.3f} bohr")
print(f"final   bond = {r_final / Bohr:.3f} bohr (HF/STO-3G opt ≈ 1.346 bohr)")

What’s coming¶

The implemented route — RHF / RKS / UHF / UKS at Γ, multi-k pure-DFT, DIIS, analytic forces, Hessians, DOS / band helpers, restart I/O, the VibeqcGPW ASE calculator, the C++ SCF host (run_rhf_scf_gpw_cpp), and the PeriodicJKMethod.GPW runner wiring — is in place. The open roadmap items are the all-electron and hybrid extensions:

GAPW augmentation — the per-atom radial-grid all-electron correction on top of the smooth-grid GPW J that closes the parity gap against CP2K on tight ionic crystals. The atomic radial × Lebedev grids already ship (periodic_gapw_atomic_grid); the multipole-compensating ρ_0 and its wiring through the SCF are the remaining work. M3a’s Ewald V_ne (the default) and M3b’s CP2K-style v_ne_convention="smeared_erfc" split are the conventions the augmentation plugs into.
CP2K parity sweep via the examples/regression/core/runner_cp2k.py subprocess runner — gated on the augmentation above.
Metallic GAPW — Fermi-Dirac smearing already ships on the multi-k GPW route (smearing_temperature=); pairing it with the augmentation for genuine all-electron metals is the open piece.
ACE-decorated K for hybrids, and PAW pseudopotentials via an on-demand paw_fetcher.py.

The full milestone map is in docs/design_periodic_gapw.md.