"""Periodic GDF SCF driver via compensated charges (compcell) + exxdiv='ewald'.
This is the production GDF driver for vibe-qc's periodic SCF. It pairs
:func:`vibeqc.aux_basis.build_lpq_compcell` (Sun, *J. Chem. Phys.* **147**,
164119 (2017); PySCF ``_CCGDFBuilder``) with PySCF's ``exxdiv='ewald'``
Madelung K-matrix shift (McClain, Sun, Chan, Berkelbach, *J. Chem.
Theory Comput.* **13**, 1209 (2017)) to match PySCF's
``KRHF(cell).density_fit()`` SCF total energy to µHa.
The driver is intentionally narrow in this initial landing:
* Closed-shell RHF only (UKS / open-shell not yet wired).
* Γ-only path (single k-point at the BZ origin); multi-k is the
next milestone (see ``run_krhf_periodic_gdf`` for the multi-k
scaffolding that this will plug into).
* Coordinated with the parallel BIPOLE chat through the
:class:`PBCMethod` enum (defined here for this branch; once the
BIPOLE driver lands the enum + dispatcher will move to a shared
``pbc_options.py``).
The Coulomb (J) and exchange (K) matrices are both built from the
compcell Lpq factor -- the "true GDF" path. This differs from
:func:`run_rhf_periodic_gamma_gdf` which on dim=3 short-circuits J to
the Ewald-3D composed builder and K to the molecular-limit real-space
kernel (see Pitfall 3 in the GDF-chat handover). For diffuse aux
(every standard JKfit aux) this distinction matters: only compcell
keeps the Lpq factor cutoff-stable, which keeps J/K stable too.
.. note:: **AFT correction is applied by default**
(``apply_aft_correction=True``); η is largely a residual knob.
The AFT long-range correction (PySCF ``_CCGDFBuilder.get_2c2e``
j2c_p subtraction) conditions the compcell metric and makes the
answer largely η-independent. With both the 2c and 3c AFT pieces
in the matched ``libint`` convention (the default, η ≈ 1.0) the
compcell path reaches sub-mHa PySCF parity on H2/vacuum-box.
For µHa parity use ``gdf_method='rsgdf'`` (range-separated GDF;
sub-µHa on H2, -0.5 µHa on LiH primitive FCC) -- the recommended
µHa path. The ``aft_ft_convention='libcint'`` combination is a
known +178 mHa foot-gun and emits a runtime warning; keep the
default ``'libint'``.
Historical (pre-AFT): without the correction, parity was
η-tunable at η ≈ 0.25-0.3 on loose-aux vacuum boxes and the
compensated metric was ill-conditioned on tight ionic cells.
References
----------
* Sun, *J. Comput. Chem.* **38**, 2399 (2017), DOI 10.1002/jcc.24890
-- periodic GDF formulation.
* Sun et al., *J. Chem. Phys.* **147**, 164119 (2017),
DOI 10.1063/1.4998644 -- eigendecomposition + threshold protocol.
* Mintmire, Sabin, Trickey, *Phys. Rev. A* **25**, 88 (1982);
Dunlap, Connolly, Sabin, *J. Chem. Phys.* **71**, 3396 (1979);
Whitten, *J. Chem. Phys.* **58**, 4496 (1973) -- modrho theory.
* McClain, Sun, Chan, Berkelbach, *J. Chem. Theory Comput.* **13**,
1209 (2017), DOI 10.1021/acs.jctc.6b01184 -- exxdiv='ewald'.
* Ye, Berkelbach, *J. Chem. Phys.* **154**, 131104 (2021),
DOI 10.1063/5.0046617 -- RSGDF alternative to compcell.
"""
from __future__ import annotations
import enum
import warnings
from dataclasses import dataclass, field
from typing import List, NamedTuple, Optional, Sequence, Tuple, Union
import numpy as np
from ._vibeqc_core import (
BasisSet,
CoulombMethod,
InitialGuess,
LatticeSumOptions,
PeriodicRHFOptions,
PeriodicSystem,
SCFIteration,
bloch_sum,
compute_kinetic_lattice,
compute_overlap_lattice,
direct_lattice_cells,
nuclear_repulsion_per_cell,
)
from .aux_basis import (
build_lpq_compcell,
build_lpq_mdf,
build_lpq_native,
build_lpq_native_fft,
default_aux_for,
make_aux_basis_set,
make_modrho_aux_basis,
)
from .guess import initial_densities_open_shell, initial_density_closed_shell
from .linear_dependence import scf_preflight_overlap_check
from .madelung import (
apply_exxdiv_ewald_to_K,
exxdiv_ewald_energy_shift,
madelung_constant_for_cell,
)
from .occupations import aufbau_occupations_per_k as _aufbau_occupations_per_k
from .periodic_rhf_ewald import _canonical_orthogonalizer
from .periodic_scf_accelerators import DynamicDamping, PeriodicSCFAccelerator
from .periodic_v_ne import compute_nuclear_lattice_dispatch
from .periodic_screened_exchange import reject_unscreened_range_separated
from .progress import ProgressLogger, resolve_progress
__all__ = [
"PBCMethod",
"PBCExxDiv",
"PBCGDFResult",
"run_pbc_gdf_rhf",
]
class PBCMethod(enum.Enum):
"""Periodic SCF algorithm selector.
Coordinated with the BIPOLE chat for the v0.8.0 release --
``GDF`` is owned by this branch, ``BIPOLE`` by the sibling
chat. The enum will move to ``vibeqc.pbc_options`` once both
drivers land.
"""
GDF = "gdf"
BIPOLE = "bipole"
class PBCExxDiv(enum.Enum):
"""Treatment of the G=0 self-image divergence in HF exchange.
* ``EWALD`` -- PySCF's default. Adds the Ewald-Madelung K-matrix
shift ``K(k) += ξ . S(k).D(k).S(k)`` per k-point, where
``ξ = a_M / L`` is the cell Madelung constant. The shift is
the standard fix for the O(1/N_k) bias in HF exchange on a
finite Monkhorst-Pack mesh (McClain et al. 2017). REQUIRED
for µHa parity with PySCF ``exxdiv='ewald'``.
* ``NONE`` -- Skip the shift. Reproduces PySCF
``exxdiv=None`` / vibe-qc's pre-2026-05-17 behaviour. Useful
for cross-comparison and for systems where the user wants
the unshifted K (e.g. for explicit charge-correction sweeps).
"""
EWALD = "ewald"
NONE = "none"
[docs]
@dataclass
class PBCGDFResult:
"""SCF result from :func:`run_pbc_gdf_rhf` / :func:`run_pbc_gdf_rks`."""
energy: float
e_electronic: float
e_nuclear: float
e_coulomb: float
e_hf_exchange: float
e_exxdiv: float
n_iter: int
converged: bool
mo_energies: np.ndarray
mo_coeffs: np.ndarray
density: np.ndarray
fock: np.ndarray
overlap: np.ndarray
hcore: np.ndarray = field(default_factory=lambda: np.empty((0, 0)))
scf_trace: List[SCFIteration] = field(default_factory=list)
aux_basis_name: str = ""
n_aux: int = 0
n_fit: int = 0
madelung_constant: float = 0.0
exxdiv: str = "ewald"
compcell_eta: float = 0.2
backend: str = "pbc-gdf-compcell"
gradient: Optional[np.ndarray] = None
e_xc: float = 0.0
functional: str = ""
def _density_from_orbitals_and_occupations(
C: np.ndarray,
occupations: np.ndarray,
) -> np.ndarray:
D = (C * np.asarray(occupations, dtype=float)[None, :]) @ C.T
return 0.5 * (D + D.T)
def _open_shell_gamma_occupy(
C_alpha,
eps_alpha,
C_beta,
eps_beta,
n_alpha,
n_beta,
smear_opts,
):
"""Γ-point per-spin occupations + densities for the open-shell GDF drivers.
``smear_opts`` disabled (T = 0) -> integer hard cutoff, **bit-identical** to
the pre-smearing path (``D_s = C_s[:, :n_s] C_s[:, :n_s]ᵀ``; occupations
returned empty). ``T > 0`` -> per-spin Fermi-Dirac with independent global
mu_a, mu_b at the single Γ k-point (:func:`apply_smearing_open_shell`,
weights = [1]); the density is built from the fractional occupations.
Returns ``(D_a, D_b, occ_a, occ_b, mu_a, mu_b, entropy)``."""
if not smear_opts.enabled:
D_a = C_alpha[:, :n_alpha] @ C_alpha[:, :n_alpha].T
D_b = (
C_beta[:, :n_beta] @ C_beta[:, :n_beta].T
if n_beta > 0
else np.zeros_like(D_a)
)
return D_a, D_b, np.empty(0), np.empty(0), 0.0, 0.0, 0.0
from .smearing import apply_smearing_open_shell
a_res, b_res = apply_smearing_open_shell(
[np.asarray(eps_alpha)],
[np.asarray(eps_beta)],
weights=[1.0],
n_alpha=n_alpha,
n_beta=n_beta,
smearing=smear_opts,
)
occ_a = np.asarray(a_res.occupations_per_k[0], dtype=float)
occ_b = np.asarray(b_res.occupations_per_k[0], dtype=float)
D_a = _density_from_orbitals_and_occupations(C_alpha, occ_a)
D_b = (
_density_from_orbitals_and_occupations(C_beta, occ_b)
if n_beta > 0
else np.zeros_like(D_a)
)
return (
D_a,
D_b,
occ_a,
occ_b,
float(a_res.mu),
float(b_res.mu),
float(a_res.entropy + b_res.entropy),
)
def _gamma_open_shell_s_squared(
n_alpha,
n_beta,
C_alpha,
C_beta,
S,
occ_a,
occ_b,
smear_opts,
) -> float:
"""<S^2> for a Γ open-shell determinant. T = 0 -> the exact integer-occupation
value (:func:`_spin_squared`, **bit-identical** to the pre-smearing path);
T > 0 -> the fractional-occupation (ensemble-UHF) value via
:func:`_multi_k_s_squared` at the single Γ k-point (weights = [1])."""
from .periodic_uhf_ewald import _spin_squared
if not smear_opts.enabled:
return float(_spin_squared(n_alpha, n_beta, C_alpha, C_beta, S))
from .periodic_k_gdf import _multi_k_s_squared
return float(
_multi_k_s_squared(
n_alpha,
n_beta,
[C_alpha],
[C_beta],
[S],
[1.0],
occ_alpha_k=[occ_a],
occ_beta_k=[occ_b],
)
)
def _gauge_lat_opts_ewald_3d(
src: LatticeSumOptions,
system: PeriodicSystem,
) -> LatticeSumOptions:
"""Return a clone of ``src`` with coulomb_method forced to EWALD_3D
for the V_ne and nuclear-repulsion lattice sums on 3D-periodic
systems. Matches PySCF's exxdiv='ewald' / CRYSTAL14 / vibe-qc's
own EWALD_3D direct path. For dim<3 this is a passthrough.
See ``periodic_rhf_gdf._gauge_lat_opts_for_v_ne_and_e_nuc`` for
the original rationale (2026-05-13 gauge-regression fix).
"""
if int(system.dim) != 3:
return src
dst = LatticeSumOptions()
for attr in dir(src):
if attr.startswith("_"):
continue
try:
value = getattr(src, attr)
except Exception:
continue
if callable(value):
continue
try:
setattr(dst, attr, value)
except Exception:
pass
dst.coulomb_method = CoulombMethod.EWALD_3D
return dst
# Dimensionless tail-resolution threshold for lifting the dense-core Γ GDF
# parity hold: rsgdf_tail_ke_cutoff / zeta_max >= this ratio (zeta_max = the
# steepest AO primitive exponent). Calibrated on the exact P01 MgO/STO-3G Γ
# RHF cell (zeta_max = 299.24) against the PySCF GDF target
# -271.049457612534 Ha:
# tail/zeta_max ~ 5.3 (1600 Ha) -> -1.7 mHa
# tail/zeta_max ~ 10.7 (3200 Ha) -> +1.3 uHa (PySCF GDF<->RSDF ~ 1 uHa)
# so >= 10 marks the mesh regime where the tight-core AO-pair FT is resolved
# and absolute parity holds.
_RSGDF_PARITY_TAIL_RATIO = 10.0
# High-|G| tail completion is dominated by same-cell tight-core AO products.
# A conservative primitive-pair FT shell/cell screen keeps dense-core
# production tails from spending hours on exponentially dead diffuse/inter-cell
# blocks while leaving the exact unscreened builder available for direct
# validation.
_RSGDF_TAIL_PAIR_FT_SCREEN = 1.0e-10
def _gamma_dense_core_gdf_parity_held(
system: PeriodicSystem,
gdf_method: str,
ao_basis: Optional[BasisSet] = None,
tail_ke_cutoff: Optional[float] = None,
rsgdf_ke_cutoff: float = 200.0,
) -> bool:
"""Return True for the Γ-only tight-core GDF parity-hold class.
H2/vacuum-box and LiH/STO-3G Γ RSGDF have PySCF parity coverage. The
release-paper P01 MgO/STO-3G Γ cell was different: at the default
``rsgdf_ke_cutoff`` the tight Mg/O core AO products are unresolved on
the dense reciprocal mesh and the electronic terms sat ~0.5 Ha from PySCF
while the Ewald nuclear term agreed. The same tight-core failure can appear
in sparse molecular-limit boxes (for example OH/STO-3G in a 20 bohr cube):
it is a basis reciprocal-resolution issue, not only a cell-density issue.
The tight-core Γ class is tagged held -- UNLESS the high-|G| tail
completion is requested with enough reciprocal support to resolve the
steepest AO pair (``tail_ke_cutoff >= _RSGDF_PARITY_TAIL_RATIO x
zeta_max``): with the tail, P01 reaches +1.3 uHa vs PySCF GDF at
``rsgdf_tail_ke_cutoff=3200`` (the validated fix; see the calibration
note on ``_RSGDF_PARITY_TAIL_RATIO``).
"""
if str(gdf_method) not in {"rsgdf", "mdf"}:
return False
if int(system.dim) != 3 or not system.unit_cell:
return False
zeta_max = _max_ao_primitive_exponent(ao_basis) if ao_basis is not None else None
tight_basis_requires_tail = (
zeta_max is not None
and zeta_max > 0.0
and _RSGDF_PARITY_TAIL_RATIO * zeta_max
> float(rsgdf_ke_cutoff) + 1e-12
)
compact_dense_core = False
if not tight_basis_requires_tail:
z_max = max(int(atom.Z) for atom in system.unit_cell)
if z_max < 8:
return False
try:
cell_vol = float(abs(np.linalg.det(np.asarray(system.lattice, dtype=float))))
except Exception:
return False
cell_vol_per_atom = cell_vol / max(len(system.unit_cell), 1)
compact_dense_core = cell_vol_per_atom < 250.0
if not compact_dense_core:
return False
# Tail completion sized to the steepest AO primitive lifts the hold
# (rsgdf only -- the mdf path has no tail plumbing).
if (
tail_ke_cutoff is not None
and zeta_max is not None
and str(gdf_method) == "rsgdf"
):
if (
zeta_max is not None
and zeta_max > 0.0
and float(tail_ke_cutoff) >= _RSGDF_PARITY_TAIL_RATIO * zeta_max
):
return False
return bool(tight_basis_requires_tail or compact_dense_core)
def _max_ao_primitive_exponent(basis: BasisSet) -> Optional[float]:
try:
return max(float(max(shell.exponents)) for shell in basis.shells())
except Exception:
return None
def _auto_rsgdf_tail_ke_cutoff(
system: PeriodicSystem,
gdf_method: str,
basis: BasisSet,
tail_ke_cutoff: Optional[float],
rsgdf_ke_cutoff: float = 200.0,
) -> Optional[float]:
"""Resolve the production default for tight-core Gamma RSGDF tails.
``None`` means production behavior: if this is the tight-core Gamma hold
class, choose a tail sized past the calibrated parity threshold. The
classifier is basis-driven so sparse molecular-limit boxes with O/Ne/Mg
STO-3G cores get the same correction as compact ionic cells. Explicit
numeric values are preserved, including ``0`` for fast diagnostic runs that
intentionally keep the parity hold.
"""
if tail_ke_cutoff is not None or str(gdf_method) != "rsgdf":
return tail_ke_cutoff
if not _gamma_dense_core_gdf_parity_held(
system,
gdf_method,
basis,
None,
rsgdf_ke_cutoff=float(rsgdf_ke_cutoff),
):
return None
zeta_max = _max_ao_primitive_exponent(basis)
if zeta_max is None or zeta_max <= 0.0:
return None
return float(
max(float(rsgdf_ke_cutoff), 1.1 * _RSGDF_PARITY_TAIL_RATIO * zeta_max)
)
def _warn_gamma_dense_core_gdf_parity_hold(
system: PeriodicSystem,
gdf_method: str,
plog: ProgressLogger,
ao_basis: Optional[BasisSet] = None,
tail_ke_cutoff: Optional[float] = None,
rsgdf_ke_cutoff: float = 200.0,
) -> bool:
if not _gamma_dense_core_gdf_parity_held(
system,
gdf_method,
ao_basis,
tail_ke_cutoff,
rsgdf_ke_cutoff=float(rsgdf_ke_cutoff),
):
return False
msg = (
"run_pbc_gdf: Gamma-only "
f"{gdf_method} absolute-energy parity is HELD for tight-core "
"basis/cell combinations at this reciprocal-mesh resolution: the "
"tight core AO-pair FT is unresolved (the P01 MgO/STO-3G audit "
"measured a ~0.5 Ha electronic offset vs PySCF at the 200 Ha default; "
"OH/STO-3G molecular-limit boxes show the same basis-resolution "
"failure while the Ewald nuclear term matches). Remediation: pass "
"rsgdf_tail_ke_cutoff >= 10 x (steepest AO primitive exponent) "
"to enable the high-|G| tail completion -- P01 reaches ~1 uHa "
"parity at rsgdf_tail_ke_cutoff=3200 -- or use a separately "
"validated multi-k route."
)
warnings.warn(msg, RuntimeWarning, stacklevel=3)
plog.info(" WARNING: " + msg)
return True
def _warn_gamma_compcell_tight_ionic_cell(
system: PeriodicSystem,
gdf_method: str,
driver: str,
) -> bool:
"""Warn when the Γ-only compcell Hartree cannot be trusted.
Γ-only compcell GDF is not sufficient for tight cells with heavy
atoms (Z > 1) where AO-pair images overlap: the Γ-only Hartree
(G=0 dropped) cannot resolve the overlap between periodic images
and the SCF converges to a non-physical fixed point (LiH FCC
primitive: +579.8 Ha at HF vs PySCF -8 Ha; the energy-sanity
guard rejects it downstream). The multi-k path is the production
route for these cells.
Heuristic: Z > 1 in a cell < 500 bohr³ total is tight enough that
AO-pair images overlap. H₂ in a 12-bohr box (1728 bohr³) is safe;
LiH FCC primitive (115 bohr³) is not.
"""
if int(system.dim) != 3 or str(gdf_method) != "compcell":
return False
if not system.unit_cell:
return False
z_max = max(atom.Z for atom in system.unit_cell)
cell_vol = float(abs(np.linalg.det(np.asarray(system.lattice))))
if z_max <= 1 or cell_vol >= 500.0:
return False
warnings.warn(
f"{driver}: Γ-only compcell GDF may be unreliable "
f"for this tight ionic cell (max Z={z_max}, cell volume "
f"{cell_vol:.0f} bohr³). The Γ-only Hartree (G=0 "
"dropped) cannot fully resolve AO-pair overlap between "
"periodic images. Prefer the multi-k path: "
"vibeqc.run_krhf_periodic_gdf(..., use_compcell=True). "
"Validated at µHa parity on LiH FCC at kmesh=(2,2,2).",
stacklevel=3,
)
return True
def _gdf_backend_with_parity_hold(backend: str, parity_held: bool) -> str:
if not parity_held or "+PARITY_HELD" in backend:
return backend
return f"{backend}+PARITY_HELD"
def _build_j_from_lpq(Lpq: np.ndarray, D: np.ndarray) -> np.ndarray:
# W_muν,κl = S_L Lpq[L,muν].conj(Lpq[L,κl]); J_muν = S_κl W_muν,κl D_κl.
# The complex combined MDF cderi [L_gauss; cderi_pw] needs the conj()
# (real D at Γ => J is real symmetric); the real compcell/rsgdf Lpq is
# left untouched so the path stays BIT-identical (conj on a real view
# perturbs einsum's contraction order at ~1e-15).
cplx = np.iscomplexobj(Lpq)
Lc = Lpq.conj() if cplx else Lpq
rho = np.einsum("Lij,ij->L", Lc, D, optimize=True)
J = np.einsum("L,Lij->ij", rho, Lpq, optimize=True)
if cplx:
J = np.real(J)
return 0.5 * (J + J.T)
def _build_k_from_lpq(Lpq: np.ndarray, D: np.ndarray) -> np.ndarray:
# K_muν = S_κl W_muκ,νl D_κl = S_L (Lpq[L].D.Lpq[L]+)_muν. As in
# _build_j_from_lpq, the real path is left bit-identical; the complex
# MDF cderi takes the conjugate transpose.
cplx = np.iscomplexobj(Lpq)
Lc = Lpq.conj() if cplx else Lpq
K = np.einsum("Lmk,kl,Lnl->mn", Lpq, D, Lc, optimize=True)
if cplx:
K = np.real(K)
return 0.5 * (K + K.T)
class _PbcGdfGammaSetup(NamedTuple):
"""Shared Γ-GDF setup outputs for ``run_pbc_gdf_{rhf,uhf,uks}``."""
S: np.ndarray
Hcore: np.ndarray
X: np.ndarray
n_kept: int
Lpq: np.ndarray
aux: BasisSet
madelung: float
e_nuc: float
gauge_lat_opts: LatticeSumOptions
def _pbc_gdf_gamma_setup(
system: PeriodicSystem,
basis: BasisSet,
lat_opts: LatticeSumOptions,
*,
aux_name: str,
aux_drop_eta: float,
exxdiv: PBCExxDiv,
gdf_method: str,
compcell_eta: float,
apply_aft_correction: bool,
aft_precision: float,
aft_ft_convention: str,
rsgdf_ke_cutoff: float,
rsgdf_tail_ke_cutoff: Optional[float],
mdf_ke_cutoff: float,
rcut_strategy: Optional[object],
rcut_precision: float,
gdf_linear_dep_threshold: float,
linear_dep_threshold: float,
fit_screen_threshold: float = 0.0,
plog: ProgressLogger,
) -> "_PbcGdfGammaSetup":
"""Shared Γ-point GDF setup for the closed- and open-shell pbc_gdf
drivers (``run_pbc_gdf_{rhf,uhf,uks}``).
Builds the one-electron integrals (``S`` / ``Hcore`` at Γ; ``V_ne``
and ``e_nuc`` forced through the Ewald-3D gauge to match PySCF's
``exxdiv='ewald'``), the canonical orthogonaliser ``X``, the ``Lpq``
cderi (compcell or rsgdf -- sharing the dense AO-pair FT with the V_ne
build on the rsgdf path), and the exxdiv Madelung constant. The
occupation-specific ``n_kept`` adequacy check stays in each driver
(RHF needs ``n_occ`` directions; UHF/UKS need ``max(na, nb)``).
``fit_screen_threshold > 0`` selects the memory-lean rsgdf fit: the
Schwarz-screened, G-chunked Γ build (build_lpq_native_fft delegating
to the Bloch builder -- handovers/HANDOVER_GDF_FIT_SCREENING.md) and
the streamed V_ne FT. The shared dense pair-FT bundle is NOT built
in that mode (it is exactly the ``(n_ao, n_ao, n_G)`` materialisation
the screen exists to avoid), so the pair FT runs once per consumer
(chunked) instead of once shared (dense).
"""
if float(fit_screen_threshold) < 0.0:
raise ValueError(
"pbc_gdf: fit_screen_threshold must be >= 0; "
f"got {fit_screen_threshold}"
)
if float(fit_screen_threshold) > 0.0 and gdf_method != "rsgdf":
# Loud, not silent: the Schwarz fit screen lives in the rsgdf
# builder; a threshold on compcell/mdf would be ignored.
raise NotImplementedError(
"pbc_gdf: fit_screen_threshold is implemented for "
f"gdf_method='rsgdf' only (got {gdf_method!r})."
)
# V_ne and the nuclear-nuclear repulsion are forced through the
# Ewald-3D gauge for 3D-periodic systems so they match PySCF's
# exxdiv='ewald' convention. ``lat_opts.coulomb_method`` continues
# to control J/K routing (compcell GDF here).
gauge_lat_opts = _gauge_lat_opts_ewald_3d(lat_opts, system)
if int(system.dim) == 3 and gauge_lat_opts is not lat_opts:
plog.info(
"V_ne / e_nuc gauge: Ewald-3D "
"(forced for 3D-periodic systems regardless of "
f"lat_opts.coulomb_method={lat_opts.coulomb_method!r})"
)
with plog.stage(
"integrals_lattice",
detail=f"S/T/V at cutoff {lat_opts.cutoff_bohr:.2f} bohr",
):
S_lat = compute_overlap_lattice(basis, system, lat_opts)
T_lat = compute_kinetic_lattice(basis, system, lat_opts)
# Dense-mesh AO-pair FT r̂_muν(G), shared between the V_ne FT
# (here) and the rsgdf cderi build (below) so the dominant
# ``ao_pair_fourier_transform_bloch`` cost runs ONCE per SCF
# instead of twice (rsgdf+dim-3 only). ``_resolve_pair_ft_shared``
# guards the (basis, ke_cutoff, cutoff_bohr) provenance.
pair_ft_bundle = None
if gdf_method == "rsgdf" and int(system.dim) == 3:
from .periodic_v_ne import compute_v_ne_ewald_3d_ft_gamma
if float(fit_screen_threshold) > 0.0:
# Memory-lean mode: no shared dense bundle (see the
# docstring note above); V_ne streams its own chunked
# pair FT and the screened fit build below chunks its
# own sweep.
V = compute_v_ne_ewald_3d_ft_gamma(
basis,
system,
gauge_lat_opts,
ke_cutoff=float(rsgdf_ke_cutoff),
stream_pair_ft=True,
)
else:
from .aux_basis import _rsgdf_dense_pair_ft
pair_ft_bundle = _rsgdf_dense_pair_ft(
basis, system, float(rsgdf_ke_cutoff), lat_opts
)
V = compute_v_ne_ewald_3d_ft_gamma(
basis,
system,
gauge_lat_opts,
ke_cutoff=float(rsgdf_ke_cutoff),
pair_ft_shared=pair_ft_bundle,
)
V_lat = None # Γ-only path: V already bloch-summed
else:
V_lat = compute_nuclear_lattice_dispatch(basis, system, gauge_lat_opts)
V = None
k_gamma = np.zeros(3)
S = np.real(bloch_sum(S_lat, k_gamma))
T = np.real(bloch_sum(T_lat, k_gamma))
if V is None:
V = np.real(bloch_sum(V_lat, k_gamma))
Hcore = 0.5 * ((T + V) + (T + V).T)
S = 0.5 * (S + S.T)
scf_preflight_overlap_check(S, plog=plog, label="S(Γ)", basis=basis)
X, n_kept = _canonical_orthogonalizer(S, linear_dep_threshold)
mol = system.unit_cell_molecule()
with plog.stage("aux_basis", detail=aux_name):
aux = make_aux_basis_set(
mol,
aux_name=aux_name,
drop_eta=float(aux_drop_eta),
)
plog.info(f"aux basis: {aux_name} ({aux.nbasis} BFs / {aux.nshells} shells)")
if gdf_method == "rsgdf":
with plog.stage(
"rsgdf_lpq",
detail=(
f"rsgdf-fft aux={aux.nbasis}, ke={rsgdf_ke_cutoff:g} Ha, "
f"fit_thr={gdf_linear_dep_threshold:.1e}"
),
):
# All-FT path: dense FFT mesh + Bloch-summed pair-FT (reusing
# the shared pair_ft_bundle). H₂ sub-µHa; LiH ~18 mHa.
aux_modrho = make_modrho_aux_basis(aux, mol)
Lpq = build_lpq_native_fft(
system,
basis,
aux_modrho,
ke_cutoff=float(rsgdf_ke_cutoff),
tail_ke_cutoff=(
float(rsgdf_tail_ke_cutoff)
if rsgdf_tail_ke_cutoff is not None
else None
),
tail_pair_ft_screen=(
# The Schwarz fit screen drops negligible pairs in
# base AND tail sweeps; the native tail pair screen
# is the unscreened path's knob only (the screened
# Bloch-delegated build rejects it loudly).
_RSGDF_TAIL_PAIR_FT_SCREEN
if (
rsgdf_tail_ke_cutoff is not None
and float(fit_screen_threshold) == 0.0
)
else 0.0
),
lat_opts=lat_opts,
linear_dep_thr=float(gdf_linear_dep_threshold),
pair_ft_shared=pair_ft_bundle,
fit_screen_threshold=float(fit_screen_threshold),
progress=plog,
)
elif gdf_method == "compcell":
with plog.stage(
"compcell_lpq",
detail=(
f"compcell aux={aux.nbasis}, eta={compcell_eta:g}, "
f"fit_thr={gdf_linear_dep_threshold:.1e}"
),
):
Lpq = build_lpq_compcell(
system,
basis,
aux,
molecule=mol,
lat_opts=lat_opts,
linear_dep_thr=float(gdf_linear_dep_threshold),
compcell_eta=float(compcell_eta),
apply_aft_correction=bool(apply_aft_correction),
aft_precision=float(aft_precision),
aft_ft_convention=str(aft_ft_convention),
rcut_strategy=rcut_strategy,
rcut_precision=float(rcut_precision),
)
elif gdf_method == "mdf":
with plog.stage(
"mdf_lpq",
detail=(
f"mdf aux={aux.nbasis}, eta={compcell_eta:g}, "
f"mdf_ke={mdf_ke_cutoff:g} Ha, fit_thr={gdf_linear_dep_threshold:.1e}"
),
):
# Mixed Density Fitting (Sun-Berkelbach 2017): compensated
# Gaussian fit (steep cores exact) + plane-wave residual. The
# builder returns the two parts separately; combine into one
# complex cderi [L_gauss (real->complex); cderi_pw] so the
# (complex-aware) _build_{j,k}_from_lpq handle both uniformly
# via W = S_L cderi[L].conj(cderi[L]).
_mdf = build_lpq_mdf(
system,
basis,
aux,
molecule=mol,
lat_opts=lat_opts,
linear_dep_thr=float(gdf_linear_dep_threshold),
compcell_eta=float(compcell_eta),
mdf_ke_cutoff=float(mdf_ke_cutoff),
rcut_strategy=rcut_strategy,
rcut_precision=float(rcut_precision),
)
Lpq = np.concatenate(
[_mdf.L_gauss.astype(np.complex128), _mdf.cderi_pw], axis=0
)
plog.info(
f"MDF cderi: {_mdf.n_kept_gauss} Gaussian + {_mdf.n_pw} PW fit vectors"
)
else:
raise ValueError(
f"_pbc_gdf_gamma_setup: gdf_method must be 'compcell', 'rsgdf', "
f"or 'mdf'; got {gdf_method!r}"
)
plog.info(
f"Lpq: {Lpq.shape[0]} fit vectors, "
f"shape=({Lpq.shape[0]}, {Lpq.shape[1]}, {Lpq.shape[2]})"
)
madelung = madelung_constant_for_cell(system) if exxdiv is PBCExxDiv.EWALD else 0.0
if exxdiv is PBCExxDiv.EWALD:
plog.info(f"exxdiv='ewald': K-shift xi = {madelung:.6f} / bohr (alpha_M/L)")
e_nuc = nuclear_repulsion_per_cell(system, gauge_lat_opts)
plog.info(f"E_nuc = {e_nuc:.10f} Ha")
return _PbcGdfGammaSetup(
S=S,
Hcore=Hcore,
X=X,
n_kept=int(n_kept),
Lpq=Lpq,
aux=aux,
madelung=float(madelung),
e_nuc=float(e_nuc),
gauge_lat_opts=gauge_lat_opts,
)
def run_pbc_gdf_rhf(
system: PeriodicSystem,
basis: BasisSet,
options: Optional[PeriodicRHFOptions] = None,
*,
functional: Optional[str] = None,
kmesh: Sequence[int] = (1, 1, 1),
aux_basis: Optional[str] = None,
aux_drop_eta: float = 0.0,
exxdiv: Union[PBCExxDiv, str] = PBCExxDiv.EWALD,
gdf_method: str = "compcell",
compcell_eta: float = 1.0,
apply_aft_correction: bool = True,
aft_precision: float = 1e-10,
aft_ft_convention: str = "libint",
rsgdf_omega: float = 0.4,
rsgdf_g_precision: float = 1e-10,
rsgdf_ke_cutoff: float = 200.0,
rsgdf_tail_ke_cutoff: Optional[float] = None,
mdf_ke_cutoff: float = 40.0,
rcut_strategy: Optional[object] = "pyscf_auto",
rcut_precision: float = 1e-8,
linear_dep_threshold: float = 1e-7,
gdf_linear_dep_threshold: float = 1e-9,
fit_screen_threshold: float = 0.0,
compute_gradient: bool = False,
progress: Union[bool, ProgressLogger, None] = None,
verbose: Optional[int] = None,
) -> PBCGDFResult:
"""Closed-shell periodic RHF/RKS via compcell GDF.
Γ-only in this initial landing. For ``kmesh != (1, 1, 1)``,
raises ``NotImplementedError`` -- the multi-k path is the next
milestone (see ``run_krhf_periodic_gdf`` which will pull
compcell Lpq + exxdiv shift in once this Γ path is validated).
Parameters
----------
system, basis
Periodic system and orbital basis (same conventions as the
rest of the periodic stack).
options
:class:`PeriodicRHFOptions`. If ``None``, defaults are used.
Relevant fields: ``max_iter``, ``conv_tol_energy``,
``conv_tol_grad``, ``use_diis``, ``diis_start_iter``,
``diis_subspace_size``, ``damping``, ``initial_guess``,
``lattice_opts``.
functional
Optional libxc functional name (e.g. ``"pbe"``, ``"pbe0"``) --
the closed-shell KS sibling of :func:`run_pbc_gdf_uks`'s
``functional``. ``None`` falls back to ``options.functional``
and, if that is empty too, the driver runs as plain RHF. For
hybrids the exact-exchange channel is the a_x-scaled GDF ``K``
with the same ``exxdiv='ewald'`` Madelung shift as RHF (the
identical convention to :func:`apply_exxdiv_ewald_to_K` on the
multi-k path); pure DFT (a_x = 0) skips ``K`` entirely. KS runs
with ``options=None`` default to :class:`PeriodicKSOptions`,
whose ``use_periodic_becke=True`` selects the periodic-Becke
grid + Γ-torus density pairing (b3f74aa9). Range-separated
hybrids raise ``NotImplementedError`` (no erf-attenuated Lpq).
kmesh
``(n1, n2, n3)`` Monkhorst-Pack mesh. Γ-only ``(1,1,1)`` in
this landing.
aux_basis
Aux basis name. Defaults to ``default_aux_for(basis.name)``.
aux_drop_eta
Per-primitive aux cull threshold (currently a no-op in
``make_aux_basis_set``).
exxdiv
:class:`PBCExxDiv` or its string value (``'ewald'`` /
``'none'``). Default ``EWALD`` (PySCF parity).
gdf_method
Lpq builder selector -- ``'compcell'`` (default) for the Sun-2017
compensated-charge path
(:func:`vibeqc.aux_basis.build_lpq_compcell`), or ``'rsgdf'`` for
the Ye-Berkelbach 2021 range-separated path
(:func:`vibeqc.aux_basis.build_lpq_native` with
``algorithm='rsgdf'``). RSGDF reaches sub-µHa PySCF parity on H₂
vacuum-box (~0.3 µHa vs the compcell+AFT 0.6 mHa); it is the
recommended path when ``compcell_eta`` tuning is awkward. The
compcell-only kwargs (``compcell_eta``, ``apply_aft_correction``,
``aft_*``, ``rcut_*``) are ignored under ``gdf_method='rsgdf'``;
the rsgdf-only kwargs (``rsgdf_omega``, ``rsgdf_g_precision``)
are ignored under ``gdf_method='compcell'``.
compcell_eta
Smooth-Gaussian exponent for the compensating basis (default
1.0; with AFT on, the result is largely η-independent). See
:func:`vibeqc.aux_basis.build_lpq_compcell`.
Only used when ``gdf_method='compcell'``.
rsgdf_omega
Range-separation parameter w (bohr⁻¹) for the
``gdf_method='rsgdf'`` path. Default ``0.4``.
SR + LR are independent of w at convergence -- the value
balances real-space cost against G-mesh cost.
rsgdf_g_precision
G-mesh truncation precision for the RSGDF LR
reciprocal-space sum (default ``1e-10``). Only used when
``gdf_method='rsgdf'``.
linear_dep_threshold
Overlap eigenvalue floor for canonical orthogonalisation.
gdf_linear_dep_threshold
Aux metric eigenvalue floor for the GDF fit
(eigendecomposition-with-threshold).
fit_screen_threshold
Cauchy-Schwarz screen on the three-centre fit (rsgdf only;
``handovers/HANDOVER_GDF_FIT_SCREENING.md``). ``0.0`` (default)
keeps the historical dense build with the shared V_ne/cderi
pair-FT bundle. Positive values select the memory-lean mode:
the fit build drops AO pairs whose Schwarz bound is below the
threshold and sweeps the G-mesh in chunks, and the V_ne FT
streams its own chunked pair FT -- the dense
``(n_ao, n_ao, n_G)`` pair-FT tensor is never materialised.
``1e-10`` reproduces the unscreened energy to < 1e-8 Ha/cell on
the pinned controls. Raises for compcell/mdf (the screen lives
in the rsgdf builder). Same option as on the multi-k
``run_krhf/krks/kuhf_periodic_gdf`` drivers.
compute_gradient
If ``True``, compute the convenience analytic/numerical GDF
gradient after SCF convergence. The default is ``False`` so
ordinary energy calculations do not pay the expensive post-SCF
gradient path.
progress, verbose
Live progress logging passthrough.
Returns
-------
PBCGDFResult
"""
kmesh = tuple(int(x) for x in kmesh)
if kmesh != (1, 1, 1):
raise NotImplementedError(
"run_pbc_gdf_rhf: only kmesh=(1,1,1) (Γ-only) is implemented "
f"in this driver; got kmesh={kmesh}. For multi-k periodic "
"SCF, see `vibeqc.run_krhf_periodic_gdf` (uses the legacy "
"bare-aux Lpq path; works but has the bare-aux divergence on "
"diffuse JKfit auxiliaries). Multi-k compcell GDF is a "
"v0.9.0 milestone (will patch run_krhf_periodic_gdf to use "
"build_lpq_compcell + apply_exxdiv_ewald_to_K)."
)
if isinstance(exxdiv, str):
exxdiv = PBCExxDiv(exxdiv)
if options is None:
# KS run (functional given here or via options) -> PeriodicKSOptions,
# so the periodic-Becke grid + torus density pairing engages by
# default (b3f74aa9); plain-RHF fallback keeps PeriodicRHFOptions.
from ._vibeqc_core import PeriodicKSOptions
options = PeriodicKSOptions() if functional else PeriodicRHFOptions()
opts = options
lat_opts: LatticeSumOptions = opts.lattice_opts
plog = resolve_progress(progress, verbose=verbose)
func_name = functional or str(getattr(opts, "functional", "") or "")
is_ks = bool(func_name)
func = None
alpha = 1.0
if is_ks:
from ._vibeqc_core import Functional
func = Functional(func_name, 1)
alpha = float(func.hf_exchange_fraction)
if bool(getattr(func, "is_range_separated", False)):
raise NotImplementedError(
"run_pbc_gdf_rhf: range-separated hybrids need an "
f"erf-attenuated Lpq cderi, which the Γ GDF path does not "
f"build yet (functional={func_name!r}). Use a global hybrid "
"(e.g. pbe0) or a pure functional."
)
if is_ks and compute_gradient:
raise NotImplementedError(
"run_pbc_gdf_rhf: compute_gradient is wired for the RHF "
"energy expression only; the KS gradient (V_xc derivative "
"terms) is not implemented on this driver."
)
# Hardening: detect obvious misuse before launching expensive
# integrals. Each check fails fast with an actionable message.
if int(system.dim) != 3:
raise NotImplementedError(
f"run_pbc_gdf_rhf: only dim=3 (full 3D periodic) is "
f"implemented; got system.dim={system.dim}. For molecular-"
"in-vacuum-box (dim=3 with large vacuum spacing) the path "
"works fine -- make sure system.dim is set to 3."
)
n_elec = system.n_electrons()
if n_elec % 2 != 0:
raise ValueError(
"run_pbc_gdf_rhf: closed-shell RHF requires even electron "
f"count; got {n_elec}. For open-shell systems, UHF/UKS "
"compcell GDF is a v0.9.0 milestone."
)
if system.multiplicity != 1:
raise ValueError(
"run_pbc_gdf_rhf: closed-shell RHF requires multiplicity=1; "
f"got {system.multiplicity}. For open-shell systems, UHF/UKS "
"compcell GDF is a v0.9.0 milestone."
)
# Charge-neutrality check: compcell construction assumes the unit
# cell is charge-neutral. For charged cells (defect calculations,
# isolated ions in vacuum boxes), the Madelung cancellation between
# nuclear and electronic G=0 contributions breaks and the SCF
# energy includes a divergent G=0 self-image term.
Q_nuc = float(sum(atom.Z for atom in system.unit_cell))
if abs(Q_nuc - n_elec) > 0.5:
raise ValueError(
"run_pbc_gdf_rhf: unit cell is not charge-neutral "
f"(Q_nuclei={Q_nuc:.0f}, n_electrons={n_elec}; net charge "
f"= {Q_nuc - n_elec:+.0f}). The compcell construction + "
"exxdiv='ewald' shift both assume neutrality; charged "
"cells get a divergent G=0 contribution that this driver "
"does NOT correct. For neutral cells with an explicit "
"charge state, use the molecular limit driver "
"(`vibeqc.run_rhf`) in a large vacuum box, or wait for "
"v0.9.0's charged-cell support."
)
# ---- Ionic-system guard: \u0393-only compcell is not sufficient
# for tight cells with heavy atoms (Z > 1) where AO-pair images
# overlap. The multi-k path is the production route for these.
_warn_gamma_compcell_tight_ionic_cell(system, gdf_method, "run_pbc_gdf_rhf")
if not apply_aft_correction and gdf_method == "compcell":
# The AFT long-range correction is now ON by default (v0.12.0).
# Disabling it reverts to the bare compcell path which requires
# manual \u03b7 tuning and diverges for tight ionic cells.
warnings.warn(
"run_pbc_gdf_rhf: apply_aft_correction=False disables the "
"AFT long-range correction (ON by default since v0.12.0). "
"Without AFT, \u03b7 is a tuning knob and tight ionic cells "
"diverge. Prefer the default (AFT on, \u03b7=1.0, "
"rcut_strategy='pyscf_auto') for production use.",
stacklevel=2,
)
if gdf_method == "compcell" and aft_ft_convention == "libcint":
warnings.warn(
"run_pbc_gdf_rhf: aft_ft_convention='libcint' uses PySCF's "
"FT convention, which is inconsistent with vibe-qc's own "
"bare lattice sum convention (libint). This produces ~+200 mHa "
"errors. Use the default aft_ft_convention='libint' instead.",
stacklevel=2,
)
n_occ = n_elec // 2
aux_name = aux_basis or default_aux_for(basis.name)
rsgdf_tail_ke_cutoff = _auto_rsgdf_tail_ke_cutoff(
system,
gdf_method,
basis,
rsgdf_tail_ke_cutoff,
rsgdf_ke_cutoff=float(rsgdf_ke_cutoff),
)
plog.info(
f"PBC-GDF {'RKS ' + func_name if is_ks else 'RHF'} / "
f"aux={aux_name}, exxdiv={exxdiv.value}, "
+ (f"alpha={alpha:g}, " if is_ks else "")
+ f"eta={compcell_eta:g}, kmesh={kmesh}, "
f"cutoff={lat_opts.cutoff_bohr:.2f} bohr"
)
if rsgdf_tail_ke_cutoff is not None and gdf_method == "rsgdf":
plog.info(
"RSGDF high-|G| tail cutoff: "
f"{float(rsgdf_tail_ke_cutoff):g} Ha"
)
parity_held = _warn_gamma_dense_core_gdf_parity_hold(
system,
gdf_method,
plog,
ao_basis=basis,
tail_ke_cutoff=rsgdf_tail_ke_cutoff,
rsgdf_ke_cutoff=float(rsgdf_ke_cutoff),
)
plog.info(f"basis: {basis.name} ({basis.nbasis} BFs / {basis.nshells} shells)")
plog.info(
"lattice cells: "
f"one-electron/GDF cutoff -> "
f"{len(direct_lattice_cells(system, lat_opts.cutoff_bohr))}, "
"nuclear cutoff -> "
f"{len(direct_lattice_cells(system, lat_opts.nuclear_cutoff_bohr))}"
)
# ---- One-electron integrals + Lpq cderi + exxdiv (shared setup) ---
setup = _pbc_gdf_gamma_setup(
system,
basis,
lat_opts,
aux_name=aux_name,
aux_drop_eta=aux_drop_eta,
exxdiv=exxdiv,
gdf_method=gdf_method,
compcell_eta=compcell_eta,
apply_aft_correction=apply_aft_correction,
aft_precision=aft_precision,
aft_ft_convention=aft_ft_convention,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
rsgdf_tail_ke_cutoff=rsgdf_tail_ke_cutoff,
mdf_ke_cutoff=mdf_ke_cutoff,
rcut_strategy=rcut_strategy,
rcut_precision=rcut_precision,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
linear_dep_threshold=linear_dep_threshold,
fit_screen_threshold=fit_screen_threshold,
plog=plog,
)
S, Hcore, X = setup.S, setup.Hcore, setup.X
Lpq, aux = setup.Lpq, setup.aux
madelung, e_nuc = setup.madelung, setup.e_nuc
gauge_lat_opts = setup.gauge_lat_opts
mol = system.unit_cell_molecule()
if n_occ > setup.n_kept:
raise RuntimeError(
"run_pbc_gdf_rhf: canonical orthogonalisation dropped too "
f"many directions (n_occ={n_occ}, n_kept={setup.n_kept})"
)
use_davidson = getattr(opts, "use_davidson", False)
dav_opts = getattr(opts, "davidson", None)
dav_dim = getattr(opts, "davidson_min_dim", 100)
use_dav = use_davidson and S.shape[0] >= dav_dim
if use_dav and dav_opts is None:
from vibeqc._vibeqc_core import DavidsonOptions
dav_opts = DavidsonOptions()
# ---- XC grid + density-set template (KS only) ----------------------
grid = None
D_set = None
_set_xc_density = None
if is_ks:
from ._vibeqc_core import GridOptions, build_grid, build_xc_periodic
from .periodic_grid import build_periodic_becke_grid
from .periodic_rhf_gdf import _density_set_gamma
grid_options = getattr(opts, "grid", None) or GridOptions()
_set_xc_density = _density_set_gamma
if bool(getattr(opts, "use_periodic_becke", False)):
grid = build_periodic_becke_grid(
system,
grid_options=grid_options,
image_radius_bohr=float(getattr(opts, "becke_image_radius_bohr", 0.0)),
)
# Periodic-Becke grid pairs with the Γ-torus density (every
# lattice block populated -> build_xc_periodic cross-cell
# mode); the home-cell-only set is the molecular-limit
# density that pairs with the molecular grid below (the
# 2026-07-09 KRKS finding class, b3f74aa9 --
# HANDOVER_AICCM_DIRECT_TORUS.md §4).
from .periodic_rhf_gdf import _density_set_torus_gamma
_set_xc_density = _density_set_torus_gamma
else:
grid = build_grid(system.unit_cell_molecule(), grid_options)
D_set = compute_overlap_lattice(basis, system, lat_opts)
k_gamma = np.zeros(3) # Γ-point Bloch phase for the V_xc fold below
def _xc_at(D_in: np.ndarray) -> Tuple[float, np.ndarray]:
"""E_xc + folded, symmetrised V_xc(Γ) for the current density."""
_set_xc_density(D_set, D_in)
xc = build_xc_periodic(basis, system, grid, func, D_set, lat_opts)
V = np.real(bloch_sum(xc.V_xc, k_gamma))
return float(xc.e_xc), 0.5 * (V + V.T)
# ---- SCF loop -----------------------------------------------------
def diagonalise(F: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
Fp = X.T @ F @ X
Fp = 0.5 * (Fp + Fp.T)
if use_dav and dav_opts is not None:
from vibeqc._vibeqc_core import davidson_solve
if dav_opts.n_eig == 0:
dav_opts.n_eig = Fp.shape[0]
if dav_opts.guess_vectors is not None:
pass # already set from previous iteration
dres = davidson_solve(Fp, dav_opts)
if not dres.converged:
raise RuntimeError(
f"Davidson did not converge after {dres.n_iter} iters"
)
eps, Cp = dres.eigenvalues, dres.eigenvectors
dav_opts.guess_vectors = Cp
else:
eps, Cp = np.linalg.eigh(Fp)
return X @ Cp, eps
def occupations_from_eps(eps: np.ndarray) -> np.ndarray:
return np.asarray(_aufbau_occupations_per_k([eps], n_occ)[0], dtype=float)
guess = getattr(opts, "initial_guess", InitialGuess.HCORE)
D_engine = initial_density_closed_shell(
mol,
basis,
n_occ,
guess,
is_periodic=True,
periodic_system=system,
lattice_opts=lat_opts,
# READ restart (Γ-only): prior g=0 cell density (pre-resolved from
# read_from, or read + projected from read_path). Ignored unless READ.
read_density=getattr(opts, "read_density", None),
read_path=getattr(opts, "read_path", ""),
)
if D_engine is not None:
plog.info(f"initial guess: {guess.name} (density via GuessEngine)")
D = D_engine
C0, eps0 = diagonalise(Hcore)
else:
plog.info(f"initial guess: {guess.name} (Hcore-diagonalise)")
C0, eps0 = diagonalise(Hcore)
occ0 = occupations_from_eps(eps0)
D = _density_from_orbitals_and_occupations(C0, occ0)
D_prev = D.copy()
damping = float(opts.damping)
if not (0.0 <= damping < 1.0):
raise ValueError(f"run_pbc_gdf_rhf: damping must be in [0, 1); got {damping}")
damper: Optional[DynamicDamping] = None
if bool(getattr(opts, "dynamic_damping", False)):
damper = DynamicDamping(
initial_alpha=damping,
alpha_min=float(getattr(opts, "dynamic_damping_min", 0.0)),
alpha_max=float(getattr(opts, "dynamic_damping_max", 0.95)),
)
use_diis = bool(opts.use_diis)
diis_start_iter = int(opts.diis_start_iter)
accel: Optional[PeriodicSCFAccelerator] = (
PeriodicSCFAccelerator(opts) if use_diis else None
)
max_iter = int(opts.max_iter)
scf_trace: List[SCFIteration] = []
result = PBCGDFResult(
energy=0.0,
e_electronic=0.0,
e_nuclear=float(e_nuc),
e_coulomb=0.0,
e_hf_exchange=0.0,
e_exxdiv=0.0,
n_iter=0,
converged=False,
mo_energies=np.empty(0),
mo_coeffs=np.empty((0, 0)),
density=D.copy(),
fock=np.empty((0, 0)),
overlap=S,
hcore=Hcore,
scf_trace=scf_trace,
aux_basis_name=aux_name,
n_aux=int(aux.nbasis),
n_fit=int(Lpq.shape[0]),
madelung_constant=float(madelung),
exxdiv=exxdiv.value,
compcell_eta=float(compcell_eta),
backend=_gdf_backend_with_parity_hold(
f"pbc-gdf-{gdf_method}" + ("-rks" if is_ks else ""), parity_held
),
functional=func_name,
)
plog.banner(
f"SCF (PBC-GDF {gdf_method}{' RKS ' + func_name if is_ks else ''})"
)
plog.info(" iter energy (Ha) dE ||[F,DS]|| DIIS")
E_prev = 0.0
C_final = C0
eps_final = eps0
for iter_idx in range(1, max_iter + 1):
diis_active = use_diis and iter_idx >= diis_start_iter
D_used = (
D
if (iter_idx == 1 or damping == 0.0 or diis_active)
else damping * D_prev + (1.0 - damping) * D
)
J = _build_j_from_lpq(Lpq, D_used)
K = _build_k_from_lpq(Lpq, D_used) if alpha != 0.0 else None
# exxdiv='ewald' K-shift (per Γ-point as a 1-element list). For
# hybrids the shift rides the exact-exchange channel, scaled by
# a_x through the -a_x/2 K contraction below -- the identical
# convention to the multi-k apply_exxdiv_ewald_to_K route and
# the real-Γ direct route's ξ_N seam. Pure DFT (a_x = 0) has no
# exact-exchange channel, so no shift enters at all.
e_exx = 0.0
if K is not None and exxdiv is PBCExxDiv.EWALD:
K = apply_exxdiv_ewald_to_K([K], [S], [D_used], madelung)[0]
e_exx = exxdiv_ewald_energy_shift(
[D_used],
[S],
madelung,
hf_exchange_fraction=alpha,
weights=[1.0],
)
E_xc = 0.0
V_xc = 0.0
if is_ks:
E_xc, V_xc = _xc_at(D_used)
F = Hcore + J - (0.5 * alpha * K if K is not None else 0.0) + V_xc
F = 0.5 * (F + F.T)
E_core = float(np.einsum("ij,ij->", D_used, Hcore))
E_J = 0.5 * float(np.einsum("ij,ij->", D_used, J))
# The exxdiv shift is already inside K above, so E_hf includes it.
E_K = (
-0.25 * alpha * float(np.einsum("ij,ij->", D_used, K))
if K is not None
else 0.0
)
E_elec = E_core + E_J + E_K + E_xc
E_total = E_elec + float(e_nuc)
FDS = F @ D_used @ S
grad = FDS - FDS.T
grad_norm = float(np.linalg.norm(grad))
dE = E_total - E_prev
converged = (
iter_idx > 1
and abs(dE) < float(opts.conv_tol_energy)
and grad_norm < float(opts.conv_tol_grad)
)
scf_trace.append(
SCFIteration(
iter=iter_idx,
energy=float(E_total),
delta_e=float(dE if iter_idx > 1 else 0.0),
grad_norm=float(grad_norm),
diis_subspace=(accel.subspace_size if accel is not None else 0),
)
)
plog.iteration(
iter_idx,
energy=float(E_total),
dE=float(dE if iter_idx > 1 else 0.0),
grad=float(grad_norm),
diis=(accel.subspace_size if accel is not None else 0),
)
if accel is not None:
F_ex = accel.extrapolate_rhf(
F,
error=grad,
density=D_used,
energy=E_total,
mo_coeffs=C_final,
mo_energies=eps_final,
n_occ=n_occ,
)
if diis_active:
F = F_ex
C_new, eps_new = diagonalise(F)
occ = occupations_from_eps(eps_new)
D_prev = D_used
D = _density_from_orbitals_and_occupations(C_new, occ)
C_final = C_new
eps_final = eps_new
if damper is not None:
damper.update(E_total)
E_prev = E_total
result.energy = E_total
result.e_electronic = E_elec
result.e_coulomb = E_J
result.e_hf_exchange = E_K
result.e_xc = E_xc
result.e_exxdiv = e_exx
result.n_iter = iter_idx
result.mo_energies = eps_new
result.mo_coeffs = C_new
result.density = D_used
result.fock = F
if converged:
# One final, clean energy evaluation on the converged D.
J_f = _build_j_from_lpq(Lpq, D)
K_f = _build_k_from_lpq(Lpq, D) if alpha != 0.0 else None
e_exx_f = 0.0
if K_f is not None and exxdiv is PBCExxDiv.EWALD:
K_f = apply_exxdiv_ewald_to_K([K_f], [S], [D], madelung)[0]
e_exx_f = exxdiv_ewald_energy_shift(
[D],
[S],
madelung,
hf_exchange_fraction=alpha,
weights=[1.0],
)
E_xc_f = 0.0
V_xc_f = 0.0
if is_ks:
E_xc_f, V_xc_f = _xc_at(D)
F_f = (
Hcore
+ J_f
- (0.5 * alpha * K_f if K_f is not None else 0.0)
+ V_xc_f
)
F_f = 0.5 * (F_f + F_f.T)
C_f, eps_f = diagonalise(F_f)
E_core_f = float(np.einsum("ij,ij->", D, Hcore))
E_J_f = 0.5 * float(np.einsum("ij,ij->", D, J_f))
E_K_f = (
-0.25 * alpha * float(np.einsum("ij,ij->", D, K_f))
if K_f is not None
else 0.0
)
E_elec_f = E_core_f + E_J_f + E_K_f + E_xc_f
result.energy = E_elec_f + float(e_nuc)
result.e_electronic = E_elec_f
result.e_coulomb = E_J_f
result.e_hf_exchange = E_K_f
result.e_xc = E_xc_f
result.e_exxdiv = e_exx_f
result.mo_energies = eps_f
result.mo_coeffs = C_f
result.density = D
result.fock = F_f
result.converged = True
if compute_gradient:
# Compute analytic gradient at the converged geometry only
# when requested; the V_ne numerical piece is too expensive
# to run unconditionally after plain energy calculations.
try:
from .periodic_gdf_gradient import compute_gdf_gradient
result.gradient = compute_gdf_gradient(
system,
basis,
result,
aux_basis_name=aux_name,
compcell_eta=float(compcell_eta),
)
except Exception:
result.gradient = None
else:
result.gradient = None
_check_energy_sanity(result, system, plog)
plog.converged(
n_iter=result.n_iter,
energy=result.energy,
converged=True,
)
return result
result.mo_coeffs = C_final
result.mo_energies = eps_final
result.converged = False
_check_energy_sanity(result, system, plog)
plog.converged(
n_iter=result.n_iter,
energy=result.energy,
converged=False,
)
return result
def run_pbc_gdf_rks(
system: PeriodicSystem,
basis: BasisSet,
options: Optional[PeriodicRHFOptions] = None,
*,
functional: Optional[str] = None,
**kwargs,
) -> PBCGDFResult:
"""Γ-only closed-shell periodic RKS via rsgdf/compcell GDF.
The closed-shell KS sibling of :func:`run_pbc_gdf_uks`: a thin
wrapper over :func:`run_pbc_gdf_rhf` that requires a ``functional``
(from the keyword or ``options.functional``). Hybrids use the
a_x-scaled GDF exchange with the ``exxdiv='ewald'`` Madelung shift;
pure functionals skip ``K`` entirely. See :func:`run_pbc_gdf_rhf`
for the full parameter list.
"""
func = functional or str(getattr(options, "functional", "") or "")
if not func:
raise ValueError("run_pbc_gdf_rks requires functional=...")
return run_pbc_gdf_rhf(
system, basis, options, functional=str(func), **kwargs
)
@dataclass
class PBCGDFUHFResult:
"""Result of a Γ-only open-shell UHF compcell/rsgdf GDF SCF.
Mirrors :class:`PBCGDFResult` (closed-shell) with per-spin a/b
orbital/density/Fock blocks and the ``<S^2>`` spin-contamination
diagnostic, matching :class:`vibeqc.periodic_uhf_ewald.PeriodicUHFEwaldResult`.
"""
energy: float
e_electronic: float
e_nuclear: float
e_coulomb: float
e_hf_exchange: float
e_exxdiv: float
n_iter: int
converged: bool
s_squared: float
s_squared_ideal: float
# a spin
mo_energies_alpha: np.ndarray
mo_coeffs_alpha: np.ndarray
density_alpha: np.ndarray
fock_alpha: np.ndarray
# b spin
mo_energies_beta: np.ndarray
mo_coeffs_beta: np.ndarray
density_beta: np.ndarray
fock_beta: np.ndarray
overlap: np.ndarray
hcore: np.ndarray = field(default_factory=lambda: np.empty((0, 0)))
scf_trace: List[SCFIteration] = field(default_factory=list)
aux_basis_name: str = ""
n_aux: int = 0
n_fit: int = 0
madelung_constant: float = 0.0
exxdiv: str = "ewald"
compcell_eta: float = 1.0
backend: str = "pbc-gdf-compcell-uhf"
# Open-shell Fermi-Dirac smearing (Γ, per-spin global mu_a/mu_b). All zero /
# empty at T = 0 (the integer-Aufbau default), so existing callers are
# behavior-neutral.
smearing_temperature: float = 0.0
fermi_level_alpha: float = 0.0
fermi_level_beta: float = 0.0
entropy: float = 0.0
free_energy: float = 0.0
occupations_alpha: np.ndarray = field(default_factory=lambda: np.empty(0))
occupations_beta: np.ndarray = field(default_factory=lambda: np.empty(0))
fock_mixing: float = 0.0
def run_pbc_gdf_uhf(
system: PeriodicSystem,
basis: BasisSet,
options: Optional[PeriodicRHFOptions] = None,
*,
aux_basis: Optional[str] = None,
aux_drop_eta: float = 0.0,
exxdiv: Union[PBCExxDiv, str] = PBCExxDiv.EWALD,
gdf_method: str = "rsgdf",
compcell_eta: float = 1.0,
apply_aft_correction: bool = True,
aft_precision: float = 1e-10,
aft_ft_convention: str = "libint",
rsgdf_omega: float = 0.4,
rsgdf_g_precision: float = 1e-10,
rsgdf_ke_cutoff: float = 200.0,
rsgdf_tail_ke_cutoff: Optional[float] = None,
mdf_ke_cutoff: float = 40.0,
rcut_strategy: Optional[object] = "pyscf_auto",
rcut_precision: float = 1e-8,
linear_dep_threshold: float = 1e-7,
gdf_linear_dep_threshold: float = 1e-9,
fit_screen_threshold: float = 0.0,
check_energy_sanity: bool = True,
progress: Union[bool, ProgressLogger, None] = None,
verbose: Optional[int] = None,
) -> PBCGDFUHFResult:
"""Γ-only open-shell periodic UHF via rsgdf/compcell GDF.
The open-shell sibling of :func:`run_pbc_gdf_rhf`: same Lpq cderi
(spin-independent), Hartree ``J`` from the total density
``D_a + D_b``, per-spin exchange ``K_s`` from the shared Lpq, and the
``exxdiv='ewald'`` Madelung K-shift applied per spin. a/b occupations
follow the molecule's ``multiplicity``
(``n_a = (n_e + mult - 1) // 2``, ``n_b = (n_e - mult + 1) // 2``).
At ``multiplicity = 1`` it reproduces :func:`run_pbc_gdf_rhf` to SCF
tolerance (the closed-shell-limit gate).
``gdf_method`` defaults to ``'rsgdf'`` (2026-07-09; the open-shell
drivers used to default to ``'compcell'`` like the RHF sibling). The
Γ-only compcell Hartree cannot resolve overlapping AO-pair images on
tight ionic cells -- its q-only cderi is exact only in the vacuum-box
limit -- and the SCF converges to a non-physical fixed point
(rocksalt LiH/STO-3G: +579.8 Ha at UHF, +1172.6 Ha at UKS-PBE, both
``converged=True``; the rsgdf lane lands at the sane -8.33 / -8.23 Ha
on the same cell). ``gdf_method='compcell'`` stays selectable for its
own development; the tight-ionic warning + the energy-sanity guard
still fence it.
Parameters are identical to :func:`run_pbc_gdf_rhf`; see its docstring.
Additionally, ``check_energy_sanity`` (default ``True``) rejects a
*converged* non-physical total energy with a ``RuntimeError`` instead
of returning it -- the CLAUDE.md §7 guard; see
:func:`_check_energy_sanity`. Pass ``False`` to bypass (diagnostics
only). Returns a :class:`PBCGDFUHFResult` with per-spin a/b blocks and
the ``<S^2>`` diagnostic.
"""
from .periodic_uhf_ewald import _spin_squared
if options is None:
options = PeriodicRHFOptions()
opts = options
# SPINLOCK. SPIN_SCHEDULE (two-phase: lock spin then release) is delegated
# before any SCF setup; the sub-runs re-enter with spinlock OFF. PATTERN_HOLD
# (MOM-hold the seeded occupied set, below) runs inside the SCF loop.
from ._vibeqc_core import SpinlockMode
from .spinlock_periodic import check_spinlock_support, run_spin_schedule
if (
getattr(opts, "spinlock_mode", SpinlockMode.OFF) == SpinlockMode.SPIN_SCHEDULE
and int(getattr(opts, "spinlock_iterations", 0)) > 0
):
return run_spin_schedule(
lambda sysx, o: run_pbc_gdf_uhf(
sysx,
basis,
o,
aux_basis=aux_basis,
aux_drop_eta=aux_drop_eta,
exxdiv=exxdiv,
gdf_method=gdf_method,
compcell_eta=compcell_eta,
apply_aft_correction=apply_aft_correction,
aft_precision=aft_precision,
aft_ft_convention=aft_ft_convention,
rsgdf_omega=rsgdf_omega,
rsgdf_g_precision=rsgdf_g_precision,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
rsgdf_tail_ke_cutoff=rsgdf_tail_ke_cutoff,
mdf_ke_cutoff=mdf_ke_cutoff,
rcut_strategy=rcut_strategy,
rcut_precision=rcut_precision,
linear_dep_threshold=linear_dep_threshold,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
fit_screen_threshold=fit_screen_threshold,
check_energy_sanity=check_energy_sanity,
progress=progress,
verbose=verbose,
),
system,
opts,
)
check_spinlock_support(
opts,
{SpinlockMode.PATTERN_HOLD, SpinlockMode.SPIN_SCHEDULE},
"the GDF UHF driver",
)
lat_opts = opts.lattice_opts
plog = resolve_progress(progress, verbose=verbose)
if not isinstance(exxdiv, PBCExxDiv):
exxdiv = PBCExxDiv(exxdiv)
if int(system.dim) != 3:
raise NotImplementedError(
"run_pbc_gdf_uhf: only dim=3 periodic systems are supported."
)
if str(gdf_method) not in ("compcell", "rsgdf", "mdf"):
raise ValueError(
f"run_pbc_gdf_uhf: gdf_method must be 'compcell', 'rsgdf', or "
f"'mdf'; got {gdf_method!r}"
)
n_elec = system.n_electrons()
mult = int(system.multiplicity)
if mult < 1:
raise ValueError(f"run_pbc_gdf_uhf: multiplicity must be >= 1, got {mult}")
if (n_elec + mult - 1) % 2 != 0:
raise ValueError(
f"run_pbc_gdf_uhf: n_electrons={n_elec} and multiplicity={mult} "
"cannot be split into integer a/b occupations."
)
n_alpha = (n_elec + mult - 1) // 2
n_beta = (n_elec - mult + 1) // 2
if n_alpha < 0 or n_beta < 0:
raise ValueError(
f"run_pbc_gdf_uhf: invalid occupations n_a={n_alpha}, "
f"n_b={n_beta} for n_e={n_elec}, mult={mult}."
)
Q_nuc = float(sum(a.Z for a in system.unit_cell))
if abs(Q_nuc - n_elec) > 0.5:
raise ValueError(
"run_pbc_gdf_uhf: cell is not charge-neutral "
f"(Q_nuclei={Q_nuc:.0f}, n_electrons={n_elec})."
)
aux_name = aux_basis or default_aux_for(basis.name)
rsgdf_tail_ke_cutoff = _auto_rsgdf_tail_ke_cutoff(
system,
gdf_method,
basis,
rsgdf_tail_ke_cutoff,
rsgdf_ke_cutoff=float(rsgdf_ke_cutoff),
)
plog.info(
f"PBC-GDF UHF / aux={aux_name}, exxdiv={exxdiv.value}, "
f"n_alpha={n_alpha}, n_beta={n_beta} (mult={mult})"
)
_warn_gamma_compcell_tight_ionic_cell(system, gdf_method, "run_pbc_gdf_uhf")
if rsgdf_tail_ke_cutoff is not None and gdf_method == "rsgdf":
plog.info(
"RSGDF high-|G| tail cutoff: "
f"{float(rsgdf_tail_ke_cutoff):g} Ha"
)
parity_held = _warn_gamma_dense_core_gdf_parity_hold(
system,
gdf_method,
plog,
ao_basis=basis,
tail_ke_cutoff=rsgdf_tail_ke_cutoff,
rsgdf_ke_cutoff=float(rsgdf_ke_cutoff),
)
# ---- One-electron integrals + Lpq cderi + exxdiv (shared setup) ---
setup = _pbc_gdf_gamma_setup(
system,
basis,
lat_opts,
aux_name=aux_name,
aux_drop_eta=aux_drop_eta,
exxdiv=exxdiv,
gdf_method=gdf_method,
compcell_eta=compcell_eta,
apply_aft_correction=apply_aft_correction,
aft_precision=aft_precision,
aft_ft_convention=aft_ft_convention,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
rsgdf_tail_ke_cutoff=rsgdf_tail_ke_cutoff,
mdf_ke_cutoff=mdf_ke_cutoff,
rcut_strategy=rcut_strategy,
rcut_precision=rcut_precision,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
linear_dep_threshold=linear_dep_threshold,
fit_screen_threshold=fit_screen_threshold,
plog=plog,
)
S, Hcore, X = setup.S, setup.Hcore, setup.X
Lpq, aux = setup.Lpq, setup.aux
madelung, e_nuc = setup.madelung, setup.e_nuc
if max(n_alpha, n_beta) > setup.n_kept:
raise RuntimeError(
f"run_pbc_gdf_uhf: orthogonalisation kept {setup.n_kept} directions; "
f"need >= {max(n_alpha, n_beta)} (n_a={n_alpha}, n_b={n_beta})."
)
use_davidson = getattr(opts, "use_davidson", False)
dav_opts = getattr(opts, "davidson", None)
dav_dim = getattr(opts, "davidson_min_dim", 100)
use_dav = use_davidson and S.shape[0] >= dav_dim
if use_dav and dav_opts is None:
from vibeqc._vibeqc_core import DavidsonOptions
dav_opts = DavidsonOptions()
# ---- SCF loop -----------------------------------------------------
def diagonalise(F: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
Fp = X.T @ F @ X
Fp = 0.5 * (Fp + Fp.T)
if use_dav and dav_opts is not None:
from vibeqc._vibeqc_core import davidson_solve
if dav_opts.n_eig == 0:
dav_opts.n_eig = Fp.shape[0]
if dav_opts.guess_vectors is not None:
pass # already set from previous iteration
dres = davidson_solve(Fp, dav_opts)
if not dres.converged:
raise RuntimeError(
f"Davidson did not converge after {dres.n_iter} iters"
)
eps, Cp = dres.eigenvalues, dres.eigenvectors
dav_opts.guess_vectors = Cp
else:
eps, Cp = np.linalg.eigh(Fp)
return X @ Cp, eps
# Open-shell Fermi-Dirac smearing (Γ, per-spin global mu_a/mu_b). T = 0 keeps
# the exact integer-Aufbau path (bit-identical to the pre-smearing driver).
smear_T = float(getattr(opts, "smearing_temperature", 0.0) or 0.0)
if smear_T < 0.0:
raise ValueError("run_pbc_gdf_uhf: smearing_temperature must be >= 0")
from .smearing import SmearingOptions as _SmearingOptions
smear_opts = _SmearingOptions.from_legacy_kwarg(smear_T)
C_alpha, eps_alpha = diagonalise(Hcore)
C_beta, eps_beta = C_alpha.copy(), eps_alpha.copy()
(D_alpha, D_beta, occ_a, occ_b, fermi_a, fermi_b, entropy_cur) = (
_open_shell_gamma_occupy(
C_alpha, eps_alpha, C_beta, eps_beta, n_alpha, n_beta, smear_opts
)
)
D_alpha_prev, D_beta_prev = D_alpha.copy(), D_beta.copy()
# ATOMSPIN broken-symmetry seed / READ restart. The GDF open-shell guess
# above is a spin-degenerate Hcore start (C_beta = C_alpha), which cannot
# break to an AFM / ferrimagnetic solution and carries no prior density.
# When a per-atom spin pattern (ATOMSPIN) or a READ restart is requested,
# overwrite the g=0 guess with the engine's per-spin density (ATOMSPIN
# requires the SAD guess; READ injects the prior per-spin g=0 density).
# Only engaged for those two cases, so the default GDF guess is unchanged.
_guess = getattr(opts, "initial_guess", InitialGuess.SAD)
_atomic_spins = getattr(opts, "atomic_spins", None) or None
if (
_atomic_spins is not None
or _guess == InitialGuess.READ
or _guess in {
InitialGuess.SAP,
InitialGuess.PATOM,
InitialGuess.HUECKEL,
InitialGuess.MINAO,
}
):
_seed_guess = InitialGuess.SAD if _guess == InitialGuess.PATOM else _guess
_split = initial_densities_open_shell(
system.unit_cell_molecule(),
basis,
n_alpha,
n_beta,
_seed_guess,
is_periodic=True,
periodic_system=system,
lattice_opts=lat_opts,
atomic_spins=_atomic_spins,
read_density_alpha=getattr(opts, "read_density_alpha", None),
read_density_beta=getattr(opts, "read_density_beta", None),
read_path=getattr(opts, "read_path", ""),
)
if _split is not None:
D_alpha, D_beta = _split
D_alpha_prev, D_beta_prev = D_alpha.copy(), D_beta.copy()
if _guess == InitialGuess.PATOM:
plog.info("initial guess: PATOM (SAD + one GDF in-field step)")
J = _build_j_from_lpq(Lpq, D_alpha + D_beta)
K_alpha = _build_k_from_lpq(Lpq, D_alpha)
K_beta = _build_k_from_lpq(Lpq, D_beta)
if exxdiv is PBCExxDiv.EWALD:
K_alpha = apply_exxdiv_ewald_to_K(
[K_alpha], [S], [D_alpha], madelung
)[0]
K_beta = apply_exxdiv_ewald_to_K(
[K_beta], [S], [D_beta], madelung
)[0]
F_alpha_seed = 0.5 * ((Hcore + J - K_alpha) + (Hcore + J - K_alpha).T)
F_beta_seed = 0.5 * ((Hcore + J - K_beta) + (Hcore + J - K_beta).T)
C_alpha, eps_alpha = diagonalise(F_alpha_seed)
C_beta, eps_beta = diagonalise(F_beta_seed)
(
D_alpha,
D_beta,
occ_a,
occ_b,
fermi_a,
fermi_b,
entropy_cur,
) = _open_shell_gamma_occupy(
C_alpha,
eps_alpha,
C_beta,
eps_beta,
n_alpha,
n_beta,
smear_opts,
)
D_alpha_prev, D_beta_prev = D_alpha.copy(), D_beta.copy()
damping = float(opts.damping)
fock_mixing_value = float(getattr(opts, "fock_mixing", 0.0) or 0.0)
if not (0.0 <= fock_mixing_value < 1.0):
raise ValueError(
"run_pbc_gdf_uhf: fock_mixing must be in [0, 1); "
f"got {fock_mixing_value}"
)
if fock_mixing_value != 0.0:
plog.info(f"fock mixing: {100.0 * fock_mixing_value:.1f}% previous Fock")
damper: Optional[DynamicDamping] = None
if bool(getattr(opts, "dynamic_damping", False)):
damper = DynamicDamping(
initial_alpha=damping,
alpha_min=float(getattr(opts, "dynamic_damping_min", 0.0)),
alpha_max=float(getattr(opts, "dynamic_damping_max", 0.95)),
)
use_diis = bool(opts.use_diis)
diis_start_iter = int(opts.diis_start_iter)
accel: Optional[PeriodicSCFAccelerator] = (
PeriodicSCFAccelerator(opts) if use_diis else None
)
scf_trace: List[SCFIteration] = []
result = PBCGDFUHFResult(
energy=0.0,
e_electronic=0.0,
e_nuclear=float(e_nuc),
e_coulomb=0.0,
e_hf_exchange=0.0,
e_exxdiv=0.0,
n_iter=0,
converged=False,
s_squared=0.0,
s_squared_ideal=0.25 * (mult - 1) * (mult + 1),
mo_energies_alpha=np.empty(0),
mo_coeffs_alpha=np.empty((0, 0)),
density_alpha=D_alpha.copy(),
fock_alpha=np.empty((0, 0)),
mo_energies_beta=np.empty(0),
mo_coeffs_beta=np.empty((0, 0)),
density_beta=D_beta.copy(),
fock_beta=np.empty((0, 0)),
overlap=S,
hcore=Hcore,
scf_trace=scf_trace,
aux_basis_name=aux_name,
n_aux=int(aux.nbasis),
n_fit=int(Lpq.shape[0]),
madelung_constant=float(madelung),
exxdiv=exxdiv.value,
compcell_eta=float(compcell_eta),
backend=_gdf_backend_with_parity_hold(
f"pbc-gdf-{gdf_method}-uhf", parity_held
),
fock_mixing=fock_mixing_value,
)
plog.banner(f"SCF (PBC-GDF {gdf_method} UHF)")
E_prev = 0.0
C_alpha_f, eps_alpha_f = C_alpha, eps_alpha
C_beta_f, eps_beta_f = C_beta, eps_beta
# SPINLOCK PATTERN_HOLD: hold the seeded broken-symmetry occupied set by
# maximum overlap (MOM) with the previous cycle for cycles
# 2..spinlock_iterations, then release -- protects an ATOMSPIN seed from
# collapsing to the symmetric solution. The no-smearing occupy is
# column-order (D_s = C_s[:, :n_s] C_s[:, :n_s]ᵀ), so reordering the held
# occupied to the front makes the fill pick them up. Mirrors the Γ UHF
# Ewald driver and the C++ molecular path.
from .mom import reorder_occupied_by_max_overlap as _mom_reorder
_pattern_hold = (
getattr(opts, "spinlock_mode", SpinlockMode.OFF) == SpinlockMode.PATTERN_HOLD
and int(getattr(opts, "spinlock_iterations", 0)) > 0
)
_spinlock_iters = int(getattr(opts, "spinlock_iterations", 0))
_ca_occ_prev = None
_cb_occ_prev = None
F_alpha_prev_mixed: Optional[np.ndarray] = None
F_beta_prev_mixed: Optional[np.ndarray] = None
for iter_idx in range(1, int(opts.max_iter) + 1):
if damper is not None:
damping = damper.alpha
# SPINLOCK PATTERN_HOLD: the accelerator (DIIS / EDIIS / ADIIS /
# KDIIS -- whatever PeriodicSCFAccelerator resolved) is suspended
# (no history recorded, no extrapolation, damping stays live)
# while the hold is active. Fock extrapolation across held-window
# iterates steers the SCF toward the symmetric attractor by
# continuous orbital rotation -- a collapse the occupation-selecting
# MOM hold cannot see -- and poisons the post-release history with
# out-of-basin iterates. The history starts fresh at release.
hold_active = _pattern_hold and iter_idx <= _spinlock_iters
diis_active = (
use_diis and iter_idx >= diis_start_iter and not hold_active
)
if iter_idx == 1 or damping == 0.0 or diis_active:
D_alpha_used, D_beta_used = D_alpha, D_beta
else:
D_alpha_used = damping * D_alpha_prev + (1.0 - damping) * D_alpha
D_beta_used = damping * D_beta_prev + (1.0 - damping) * D_beta
J = _build_j_from_lpq(Lpq, D_alpha_used + D_beta_used)
K_alpha = _build_k_from_lpq(Lpq, D_alpha_used)
K_beta = _build_k_from_lpq(Lpq, D_beta_used)
e_exx = 0.0
if exxdiv is PBCExxDiv.EWALD:
K_alpha = apply_exxdiv_ewald_to_K([K_alpha], [S], [D_alpha_used], madelung)[
0
]
K_beta = apply_exxdiv_ewald_to_K([K_beta], [S], [D_beta_used], madelung)[0]
e_exx = exxdiv_ewald_energy_shift(
[D_alpha_used], [S], madelung, hf_exchange_fraction=1.0, weights=[1.0]
) + exxdiv_ewald_energy_shift(
[D_beta_used], [S], madelung, hf_exchange_fraction=1.0, weights=[1.0]
)
F_alpha = 0.5 * ((Hcore + J - K_alpha) + (Hcore + J - K_alpha).T)
F_beta = 0.5 * ((Hcore + J - K_beta) + (Hcore + J - K_beta).T)
# UHF energy: 1/2Tr[(Da+Db)Hcore] + 1/2Tr[Da Fa] + 1/2Tr[Db Fb]
# (the exxdiv shift is already inside the K_s folded into F_s).
E_elec = (
0.5 * float(np.einsum("ij,ij->", D_alpha_used + D_beta_used, Hcore))
+ 0.5 * float(np.einsum("ij,ij->", D_alpha_used, F_alpha))
+ 0.5 * float(np.einsum("ij,ij->", D_beta_used, F_beta))
)
E_total = E_elec + float(e_nuc)
FDS_a = F_alpha @ D_alpha_used @ S
FDS_b = F_beta @ D_beta_used @ S
grad_a = FDS_a - FDS_a.T
grad_b = FDS_b - FDS_b.T
grad_norm = float(
np.sqrt(np.linalg.norm(grad_a) ** 2 + np.linalg.norm(grad_b) ** 2)
)
dE = E_total - E_prev
converged = (
iter_idx > 1
and abs(dE) < float(opts.conv_tol_energy)
and grad_norm < float(opts.conv_tol_grad)
)
scf_trace.append(
SCFIteration(
iter=iter_idx,
energy=float(E_total),
delta_e=float(dE if iter_idx > 1 else 0.0),
grad_norm=float(grad_norm),
diis_subspace=(accel.subspace_size if accel is not None else 0),
)
)
plog.iteration(
iter_idx,
energy=float(E_total),
dE=float(dE if iter_idx > 1 else 0.0),
grad=float(grad_norm),
diis=(accel.subspace_size if accel is not None else 0),
)
result.energy = E_total
result.e_electronic = E_elec
result.e_coulomb = 0.5 * float(
np.einsum("ij,ij->", D_alpha_used + D_beta_used, J)
)
result.e_hf_exchange = -0.5 * (
float(np.einsum("ij,ij->", D_alpha_used, K_alpha))
+ float(np.einsum("ij,ij->", D_beta_used, K_beta))
)
result.e_exxdiv = e_exx
result.n_iter = iter_idx
result.mo_energies_alpha = eps_alpha_f
result.mo_coeffs_alpha = C_alpha_f
result.density_alpha = D_alpha_used
result.fock_alpha = F_alpha
result.mo_energies_beta = eps_beta_f
result.mo_coeffs_beta = C_beta_f
result.density_beta = D_beta_used
result.fock_beta = F_beta
result.smearing_temperature = smear_T
result.entropy = float(entropy_cur)
result.free_energy = E_total - smear_T * float(entropy_cur)
result.fermi_level_alpha = float(fermi_a)
result.fermi_level_beta = float(fermi_b)
result.occupations_alpha = np.asarray(occ_a, dtype=float)
result.occupations_beta = np.asarray(occ_b, dtype=float)
if converged:
result.converged = True
result.s_squared = _gamma_open_shell_s_squared(
n_alpha, n_beta, C_alpha_f, C_beta_f, S, occ_a, occ_b, smear_opts
)
if check_energy_sanity:
_check_energy_sanity(
result,
system,
plog,
driver="run_pbc_gdf_uhf",
raise_if_converged=True,
)
plog.converged(n_iter=iter_idx, energy=E_total, converged=True)
return result
# Skipped entirely (not even recorded) while the PATTERN_HOLD window
# is active; see the hold_active note at the loop head.
if accel is not None and not hold_active:
F_alpha_ex, F_beta_ex = accel.extrapolate_uhf(
F_alpha,
F_beta,
error_alpha=grad_a,
error_beta=grad_b,
density_alpha=D_alpha_used,
density_beta=D_beta_used,
energy=E_total,
mo_coeffs_alpha=C_alpha_f,
mo_coeffs_beta=C_beta_f,
mo_energies_alpha=eps_alpha_f,
mo_energies_beta=eps_beta_f,
n_alpha=n_alpha,
n_beta=n_beta,
)
if diis_active:
F_alpha, F_beta = F_alpha_ex, F_beta_ex
if fock_mixing_value != 0.0:
if F_alpha_prev_mixed is not None and F_beta_prev_mixed is not None:
F_alpha_mix = (
(1.0 - fock_mixing_value) * F_alpha
+ fock_mixing_value * F_alpha_prev_mixed
)
F_beta_mix = (
(1.0 - fock_mixing_value) * F_beta
+ fock_mixing_value * F_beta_prev_mixed
)
F_alpha = 0.5 * (F_alpha_mix + F_alpha_mix.T)
F_beta = 0.5 * (F_beta_mix + F_beta_mix.T)
F_alpha_prev_mixed = F_alpha.copy()
F_beta_prev_mixed = F_beta.copy()
C_alpha_f, eps_alpha_f = diagonalise(F_alpha)
C_beta_f, eps_beta_f = diagonalise(F_beta)
if _pattern_hold and 1 < iter_idx <= _spinlock_iters:
if n_alpha > 0 and _ca_occ_prev is not None:
C_alpha_f, eps_alpha_f = _mom_reorder(
C_alpha_f, eps_alpha_f, S, _ca_occ_prev, n_alpha
)
if n_beta > 0 and _cb_occ_prev is not None:
C_beta_f, eps_beta_f = _mom_reorder(
C_beta_f, eps_beta_f, S, _cb_occ_prev, n_beta
)
if _pattern_hold and iter_idx <= _spinlock_iters:
_ca_occ_prev = (
np.asarray(C_alpha_f[:, :n_alpha]).copy() if n_alpha > 0 else None
)
_cb_occ_prev = (
np.asarray(C_beta_f[:, :n_beta]).copy() if n_beta > 0 else None
)
D_alpha_prev, D_beta_prev = D_alpha_used, D_beta_used
(D_alpha, D_beta, occ_a, occ_b, fermi_a, fermi_b, entropy_cur) = (
_open_shell_gamma_occupy(
C_alpha_f,
eps_alpha_f,
C_beta_f,
eps_beta_f,
n_alpha,
n_beta,
smear_opts,
)
)
if damper is not None:
damper.update(E_total)
E_prev = E_total
result.s_squared = float(
_gamma_open_shell_s_squared(
n_alpha, n_beta, C_alpha_f, C_beta_f, S, occ_a, occ_b, smear_opts
)
)
result.converged = False
if check_energy_sanity:
_check_energy_sanity(
result,
system,
plog,
driver="run_pbc_gdf_uhf",
raise_if_converged=True,
)
plog.converged(n_iter=result.n_iter, energy=result.energy, converged=False)
return result
@dataclass
class PBCGDFUKSResult:
"""Result of a Γ-only open-shell UKS compcell/rsgdf GDF SCF.
As :class:`PBCGDFUHFResult` plus the XC energy ``e_xc`` and the
``functional`` name. ``e_hf_exchange`` is the (a-scaled) HF-exchange
energy for hybrids; ``0`` for pure DFT.
"""
energy: float
e_electronic: float
e_nuclear: float
e_coulomb: float
e_hf_exchange: float
e_xc: float
e_exxdiv: float
n_iter: int
converged: bool
s_squared: float
s_squared_ideal: float
functional: str
mo_energies_alpha: np.ndarray
mo_coeffs_alpha: np.ndarray
density_alpha: np.ndarray
fock_alpha: np.ndarray
mo_energies_beta: np.ndarray
mo_coeffs_beta: np.ndarray
density_beta: np.ndarray
fock_beta: np.ndarray
overlap: np.ndarray
hcore: np.ndarray = field(default_factory=lambda: np.empty((0, 0)))
scf_trace: List[SCFIteration] = field(default_factory=list)
aux_basis_name: str = ""
n_aux: int = 0
n_fit: int = 0
madelung_constant: float = 0.0
exxdiv: str = "ewald"
compcell_eta: float = 1.0
backend: str = "pbc-gdf-compcell-uks"
# Open-shell Fermi-Dirac smearing (Γ, per-spin global mu_a/mu_b); zero/empty
# at T = 0, so existing callers are behavior-neutral.
smearing_temperature: float = 0.0
fermi_level_alpha: float = 0.0
fermi_level_beta: float = 0.0
entropy: float = 0.0
free_energy: float = 0.0
occupations_alpha: np.ndarray = field(default_factory=lambda: np.empty(0))
occupations_beta: np.ndarray = field(default_factory=lambda: np.empty(0))
fock_mixing: float = 0.0
def run_pbc_gdf_uks(
system: PeriodicSystem,
basis: BasisSet,
options: Optional[PeriodicRHFOptions] = None,
*,
functional: Optional[str] = None,
aux_basis: Optional[str] = None,
aux_drop_eta: float = 0.0,
exxdiv: Union[PBCExxDiv, str] = PBCExxDiv.EWALD,
gdf_method: str = "rsgdf",
compcell_eta: float = 1.0,
apply_aft_correction: bool = True,
aft_precision: float = 1e-10,
aft_ft_convention: str = "libint",
rsgdf_omega: float = 0.4,
rsgdf_g_precision: float = 1e-10,
rsgdf_ke_cutoff: float = 200.0,
rsgdf_tail_ke_cutoff: Optional[float] = None,
mdf_ke_cutoff: float = 40.0,
rcut_strategy: Optional[object] = "pyscf_auto",
rcut_precision: float = 1e-8,
linear_dep_threshold: float = 1e-7,
gdf_linear_dep_threshold: float = 1e-9,
fit_screen_threshold: float = 0.0,
check_energy_sanity: bool = True,
progress: Union[bool, ProgressLogger, None] = None,
verbose: Optional[int] = None,
) -> PBCGDFUKSResult:
"""Γ-only open-shell periodic UKS via rsgdf/compcell GDF.
The DFT sibling of :func:`run_pbc_gdf_uhf`: same spin-independent Lpq,
Hartree ``J`` from ``D_a + D_b``, native per-spin libxc ``V_xc`` from
:func:`build_xc_periodic_uks` on the periodic Becke grid, and -- for
hybrids -- the a-scaled per-spin exchange ``K_s`` from the shared Lpq
with the ``exxdiv='ewald'`` Madelung shift. The HF-exchange fraction
``a`` is the functional's; pure DFT (``a = 0``) skips ``K`` entirely.
``functional`` selects the XC functional (e.g. ``"pbe"``, ``"b3lyp"``,
``"pbe0"``); ``None`` falls back to ``options.functional`` and, if that
is empty too, the driver runs as plain UHF (``a = 1``, no ``V_xc``) --
so ``run_pbc_gdf_uks(functional=None)`` reproduces
:func:`run_pbc_gdf_uhf`/:func:`run_pbc_gdf_rhf`. a/b occupations follow
``multiplicity``. ``check_energy_sanity`` (default ``True``) rejects a
*converged* non-physical total energy with a ``RuntimeError`` instead
of returning it -- the CLAUDE.md §7 guard; see
:func:`_check_energy_sanity`. Pass ``False`` to bypass (diagnostics
only). Returns a :class:`PBCGDFUKSResult`.
``gdf_method`` defaults to ``'rsgdf'`` (2026-07-09; see
:func:`run_pbc_gdf_uhf` -- the Γ-only compcell Hartree returned a
converged +1172.6 Ha on rocksalt LiH/STO-3G at UKS-PBE). KS runs with
``options=None`` default to :class:`PeriodicKSOptions`, whose
``use_periodic_becke=True`` selects the periodic-Becke grid + Γ-torus
density pairing (b3f74aa9); the old ``PeriodicRHFOptions`` fallback
silently evaluated XC in the v0.8.x molecular-grid convention
(-7.9638 vs -8.2339 Ha on the same LiH fixture).
"""
from ._vibeqc_core import (
Functional,
GridOptions,
PeriodicKSOptions,
build_grid,
build_xc_periodic_uks,
)
from .periodic_grid import build_periodic_becke_grid
from .periodic_rhf_gdf import _density_set_gamma
from .periodic_uhf_ewald import _spin_squared
if options is None:
# KS run (functional given here or via options) -> PeriodicKSOptions,
# so the periodic-Becke grid + torus density pairing engages by
# default; plain-UHF fallback keeps PeriodicRHFOptions.
options = (
PeriodicKSOptions() if functional else PeriodicRHFOptions()
)
opts = options
# SPINLOCK. SPIN_SCHEDULE (two-phase) is delegated before SCF setup;
# PATTERN_HOLD (MOM-hold the seeded occupied set) runs inside the SCF loop.
from ._vibeqc_core import SpinlockMode
from .spinlock_periodic import check_spinlock_support, run_spin_schedule
if (
getattr(opts, "spinlock_mode", SpinlockMode.OFF) == SpinlockMode.SPIN_SCHEDULE
and int(getattr(opts, "spinlock_iterations", 0)) > 0
):
return run_spin_schedule(
lambda sysx, o: run_pbc_gdf_uks(
sysx,
basis,
o,
functional=functional,
aux_basis=aux_basis,
aux_drop_eta=aux_drop_eta,
exxdiv=exxdiv,
gdf_method=gdf_method,
compcell_eta=compcell_eta,
apply_aft_correction=apply_aft_correction,
aft_precision=aft_precision,
aft_ft_convention=aft_ft_convention,
rsgdf_omega=rsgdf_omega,
rsgdf_g_precision=rsgdf_g_precision,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
rsgdf_tail_ke_cutoff=rsgdf_tail_ke_cutoff,
mdf_ke_cutoff=mdf_ke_cutoff,
rcut_strategy=rcut_strategy,
rcut_precision=rcut_precision,
linear_dep_threshold=linear_dep_threshold,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
fit_screen_threshold=fit_screen_threshold,
check_energy_sanity=check_energy_sanity,
progress=progress,
verbose=verbose,
),
system,
opts,
)
check_spinlock_support(
opts,
{SpinlockMode.PATTERN_HOLD, SpinlockMode.SPIN_SCHEDULE},
"the GDF UKS driver",
)
lat_opts = opts.lattice_opts
plog = resolve_progress(progress, verbose=verbose)
if not isinstance(exxdiv, PBCExxDiv):
exxdiv = PBCExxDiv(exxdiv)
if int(system.dim) != 3:
raise NotImplementedError("run_pbc_gdf_uks: only dim=3 is supported.")
if str(gdf_method) not in ("compcell", "rsgdf", "mdf"):
raise ValueError(
f"run_pbc_gdf_uks: gdf_method must be 'compcell', 'rsgdf', or "
f"'mdf'; got {gdf_method!r}"
)
func_name = functional or str(getattr(opts, "functional", "") or "")
is_ks = bool(func_name)
func = Functional(func_name, 2) if is_ks else None # spin-polarized
# Same fail-closed rule as run_pbc_gdf_rhf's erf-attenuated-Lpq
# guard: the fitted K is full-range only.
reject_unscreened_range_separated(func, where="run_pbc_gdf_uks")
alpha = float(func.hf_exchange_fraction) if func is not None else 1.0
n_elec = system.n_electrons()
mult = int(system.multiplicity)
if mult < 1:
raise ValueError(f"run_pbc_gdf_uks: multiplicity must be >= 1, got {mult}")
if (n_elec + mult - 1) % 2 != 0:
raise ValueError(
f"run_pbc_gdf_uks: n_electrons={n_elec} and multiplicity={mult} "
"cannot be split into integer a/b occupations."
)
n_alpha = (n_elec + mult - 1) // 2
n_beta = (n_elec - mult + 1) // 2
Q_nuc = float(sum(a.Z for a in system.unit_cell))
if abs(Q_nuc - n_elec) > 0.5:
raise ValueError(
f"run_pbc_gdf_uks: cell is not charge-neutral "
f"(Q_nuclei={Q_nuc:.0f}, n_electrons={n_elec})."
)
aux_name = aux_basis or default_aux_for(basis.name)
rsgdf_tail_ke_cutoff = _auto_rsgdf_tail_ke_cutoff(
system,
gdf_method,
basis,
rsgdf_tail_ke_cutoff,
rsgdf_ke_cutoff=float(rsgdf_ke_cutoff),
)
plog.info(
f"PBC-GDF UKS / func={func_name or '(none->UHF)'}, alpha={alpha:g}, "
f"aux={aux_name}, n_alpha={n_alpha}, n_beta={n_beta} (mult={mult})"
)
_warn_gamma_compcell_tight_ionic_cell(system, gdf_method, "run_pbc_gdf_uks")
if rsgdf_tail_ke_cutoff is not None and gdf_method == "rsgdf":
plog.info(
"RSGDF high-|G| tail cutoff: "
f"{float(rsgdf_tail_ke_cutoff):g} Ha"
)
parity_held = _warn_gamma_dense_core_gdf_parity_hold(
system,
gdf_method,
plog,
ao_basis=basis,
tail_ke_cutoff=rsgdf_tail_ke_cutoff,
rsgdf_ke_cutoff=float(rsgdf_ke_cutoff),
)
# ---- One-electron integrals + Lpq cderi + exxdiv (shared setup) ---
setup = _pbc_gdf_gamma_setup(
system,
basis,
lat_opts,
aux_name=aux_name,
aux_drop_eta=aux_drop_eta,
exxdiv=exxdiv,
gdf_method=gdf_method,
compcell_eta=compcell_eta,
apply_aft_correction=apply_aft_correction,
aft_precision=aft_precision,
aft_ft_convention=aft_ft_convention,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
rsgdf_tail_ke_cutoff=rsgdf_tail_ke_cutoff,
mdf_ke_cutoff=mdf_ke_cutoff,
rcut_strategy=rcut_strategy,
rcut_precision=rcut_precision,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
linear_dep_threshold=linear_dep_threshold,
fit_screen_threshold=fit_screen_threshold,
plog=plog,
)
S, Hcore, X = setup.S, setup.Hcore, setup.X
Lpq, aux = setup.Lpq, setup.aux
madelung, e_nuc = setup.madelung, setup.e_nuc
if max(n_alpha, n_beta) > setup.n_kept:
raise RuntimeError(
f"run_pbc_gdf_uks: orthogonalisation kept {setup.n_kept} directions; "
f"need >= {max(n_alpha, n_beta)} (n_a={n_alpha}, n_b={n_beta})."
)
use_davidson = getattr(opts, "use_davidson", False)
dav_opts = getattr(opts, "davidson", None)
dav_dim = getattr(opts, "davidson_min_dim", 100)
use_dav = use_davidson and S.shape[0] >= dav_dim
if use_dav and dav_opts is None:
from vibeqc._vibeqc_core import DavidsonOptions
dav_opts = DavidsonOptions()
# ---- XC grid + per-spin density-set templates (KS only) -----------
grid = None
D_alpha_set = D_beta_set = None
_set_xc_density = _density_set_gamma
if is_ks:
grid_options = getattr(opts, "grid", None) or GridOptions()
if bool(getattr(opts, "use_periodic_becke", False)):
grid = build_periodic_becke_grid(
system,
grid_options=grid_options,
image_radius_bohr=float(getattr(opts, "becke_image_radius_bohr", 0.0)),
)
# Periodic-Becke grid pairs with the Γ-torus density (every
# lattice block populated -> build_xc_periodic_uks cross-cell
# mode); the home-cell-only set is the molecular-limit
# density that pairs with the molecular grid below (the
# 2026-07-09 KRKS finding class, spin-polarised flavour --
# HANDOVER_AICCM_DIRECT_TORUS.md §4).
from .periodic_rhf_gdf import _density_set_torus_gamma
_set_xc_density = _density_set_torus_gamma
else:
grid = build_grid(system.unit_cell_molecule(), grid_options)
D_alpha_set = compute_overlap_lattice(basis, system, lat_opts)
D_beta_set = compute_overlap_lattice(basis, system, lat_opts)
k_gamma = np.zeros(3) # Γ-point Bloch phase for the V_xc fold below
# ---- SCF loop -----------------------------------------------------
def diagonalise(F: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
Fp = X.T @ F @ X
Fp = 0.5 * (Fp + Fp.T)
if use_dav and dav_opts is not None:
from vibeqc._vibeqc_core import davidson_solve
if dav_opts.n_eig == 0:
dav_opts.n_eig = Fp.shape[0]
if dav_opts.guess_vectors is not None:
pass # already set from previous iteration
dres = davidson_solve(Fp, dav_opts)
if not dres.converged:
raise RuntimeError(
f"Davidson did not converge after {dres.n_iter} iters"
)
eps, Cp = dres.eigenvalues, dres.eigenvectors
dav_opts.guess_vectors = Cp
else:
eps, Cp = np.linalg.eigh(Fp)
return X @ Cp, eps
# Open-shell Fermi-Dirac smearing (Γ, per-spin global mu_a/mu_b). T = 0 keeps
# the exact integer-Aufbau path (bit-identical to the pre-smearing driver).
smear_T = float(getattr(opts, "smearing_temperature", 0.0) or 0.0)
if smear_T < 0.0:
raise ValueError("run_pbc_gdf_uks: smearing_temperature must be >= 0")
from .smearing import SmearingOptions as _SmearingOptions
smear_opts = _SmearingOptions.from_legacy_kwarg(smear_T)
C_alpha, eps_alpha = diagonalise(Hcore)
C_beta, eps_beta = C_alpha.copy(), eps_alpha.copy()
(D_alpha, D_beta, occ_a, occ_b, fermi_a, fermi_b, entropy_cur) = (
_open_shell_gamma_occupy(
C_alpha, eps_alpha, C_beta, eps_beta, n_alpha, n_beta, smear_opts
)
)
D_alpha_prev, D_beta_prev = D_alpha.copy(), D_beta.copy()
# ATOMSPIN broken-symmetry seed / READ restart. See run_pbc_gdf_uhf -- the
# GDF open-shell guess is a spin-degenerate Hcore start; a per-atom spin
# pattern (ATOMSPIN) or a prior per-spin density (READ) must come from the
# engine. Engaged only for those cases, so the default GDF guess is unchanged.
_guess = getattr(opts, "initial_guess", InitialGuess.SAD)
_atomic_spins = getattr(opts, "atomic_spins", None) or None
if (
_atomic_spins is not None
or _guess == InitialGuess.READ
or _guess in {
InitialGuess.SAP,
InitialGuess.PATOM,
InitialGuess.HUECKEL,
InitialGuess.MINAO,
}
):
_seed_guess = InitialGuess.SAD if _guess == InitialGuess.PATOM else _guess
_split = initial_densities_open_shell(
system.unit_cell_molecule(),
basis,
n_alpha,
n_beta,
_seed_guess,
is_periodic=True,
periodic_system=system,
lattice_opts=lat_opts,
atomic_spins=_atomic_spins,
read_density_alpha=getattr(opts, "read_density_alpha", None),
read_density_beta=getattr(opts, "read_density_beta", None),
read_path=getattr(opts, "read_path", ""),
)
if _split is not None:
D_alpha, D_beta = _split
D_alpha_prev, D_beta_prev = D_alpha.copy(), D_beta.copy()
if _guess == InitialGuess.PATOM:
plog.info("initial guess: PATOM (SAD + one GDF in-field step)")
J = _build_j_from_lpq(Lpq, D_alpha + D_beta)
K_alpha = _build_k_from_lpq(Lpq, D_alpha)
K_beta = _build_k_from_lpq(Lpq, D_beta)
if exxdiv is PBCExxDiv.EWALD:
K_alpha = apply_exxdiv_ewald_to_K(
[K_alpha], [S], [D_alpha], madelung
)[0]
K_beta = apply_exxdiv_ewald_to_K(
[K_beta], [S], [D_beta], madelung
)[0]
F_alpha_seed = 0.5 * ((Hcore + J - K_alpha) + (Hcore + J - K_alpha).T)
F_beta_seed = 0.5 * ((Hcore + J - K_beta) + (Hcore + J - K_beta).T)
C_alpha, eps_alpha = diagonalise(F_alpha_seed)
C_beta, eps_beta = diagonalise(F_beta_seed)
(
D_alpha,
D_beta,
occ_a,
occ_b,
fermi_a,
fermi_b,
entropy_cur,
) = _open_shell_gamma_occupy(
C_alpha,
eps_alpha,
C_beta,
eps_beta,
n_alpha,
n_beta,
smear_opts,
)
D_alpha_prev, D_beta_prev = D_alpha.copy(), D_beta.copy()
damping = float(opts.damping)
fock_mixing_value = float(getattr(opts, "fock_mixing", 0.0) or 0.0)
if not (0.0 <= fock_mixing_value < 1.0):
raise ValueError(
"run_pbc_gdf_uks: fock_mixing must be in [0, 1); "
f"got {fock_mixing_value}"
)
if fock_mixing_value != 0.0:
plog.info(f"fock mixing: {100.0 * fock_mixing_value:.1f}% previous Fock")
damper: Optional[DynamicDamping] = None
if bool(getattr(opts, "dynamic_damping", False)):
damper = DynamicDamping(
initial_alpha=damping,
alpha_min=float(getattr(opts, "dynamic_damping_min", 0.0)),
alpha_max=float(getattr(opts, "dynamic_damping_max", 0.95)),
)
use_diis = bool(opts.use_diis)
diis_start_iter = int(opts.diis_start_iter)
accel: Optional[PeriodicSCFAccelerator] = (
PeriodicSCFAccelerator(opts) if use_diis else None
)
scf_trace: List[SCFIteration] = []
result = PBCGDFUKSResult(
energy=0.0,
e_electronic=0.0,
e_nuclear=float(e_nuc),
e_coulomb=0.0,
e_hf_exchange=0.0,
e_xc=0.0,
e_exxdiv=0.0,
n_iter=0,
converged=False,
s_squared=0.0,
s_squared_ideal=0.25 * (mult - 1) * (mult + 1),
functional=func_name,
mo_energies_alpha=np.empty(0),
mo_coeffs_alpha=np.empty((0, 0)),
density_alpha=D_alpha.copy(),
fock_alpha=np.empty((0, 0)),
mo_energies_beta=np.empty(0),
mo_coeffs_beta=np.empty((0, 0)),
density_beta=D_beta.copy(),
fock_beta=np.empty((0, 0)),
overlap=S,
hcore=Hcore,
scf_trace=scf_trace,
aux_basis_name=aux_name,
n_aux=int(aux.nbasis),
n_fit=int(Lpq.shape[0]),
madelung_constant=float(madelung),
exxdiv=exxdiv.value,
compcell_eta=float(compcell_eta),
backend=_gdf_backend_with_parity_hold(
f"pbc-gdf-{gdf_method}-{'uks' if is_ks else 'uhf'}",
parity_held,
),
fock_mixing=fock_mixing_value,
)
plog.banner(f"SCF (PBC-GDF {gdf_method} {'UKS ' + func_name if is_ks else 'UHF'})")
E_prev = 0.0
C_alpha_f, eps_alpha_f = C_alpha, eps_alpha
C_beta_f, eps_beta_f = C_beta, eps_beta
# SPINLOCK PATTERN_HOLD: see run_pbc_gdf_uhf. Hold the seeded broken-symmetry
# occupied set by MOM for cycles 2..spinlock_iterations, then release. The
# no-smearing occupy is column-order, so reordering the held occupied to the
# front makes the fill pick them up.
from .mom import reorder_occupied_by_max_overlap as _mom_reorder
_pattern_hold = (
getattr(opts, "spinlock_mode", SpinlockMode.OFF) == SpinlockMode.PATTERN_HOLD
and int(getattr(opts, "spinlock_iterations", 0)) > 0
)
_spinlock_iters = int(getattr(opts, "spinlock_iterations", 0))
_ca_occ_prev = None
_cb_occ_prev = None
F_alpha_prev_mixed: Optional[np.ndarray] = None
F_beta_prev_mixed: Optional[np.ndarray] = None
for iter_idx in range(1, int(opts.max_iter) + 1):
if damper is not None:
damping = damper.alpha
# SPINLOCK PATTERN_HOLD: the accelerator (DIIS / EDIIS / ADIIS /
# KDIIS -- whatever PeriodicSCFAccelerator resolved) is suspended
# (no history recorded, no extrapolation, damping stays live)
# while the hold is active. Fock extrapolation across held-window
# iterates steers the SCF toward the symmetric attractor by
# continuous orbital rotation -- a collapse the occupation-selecting
# MOM hold cannot see -- and poisons the post-release history with
# out-of-basin iterates. The history starts fresh at release.
hold_active = _pattern_hold and iter_idx <= _spinlock_iters
diis_active = (
use_diis and iter_idx >= diis_start_iter and not hold_active
)
if iter_idx == 1 or damping == 0.0 or diis_active:
D_alpha_used, D_beta_used = D_alpha, D_beta
else:
D_alpha_used = damping * D_alpha_prev + (1.0 - damping) * D_alpha
D_beta_used = damping * D_beta_prev + (1.0 - damping) * D_beta
J = _build_j_from_lpq(Lpq, D_alpha_used + D_beta_used)
K_alpha = _build_k_from_lpq(Lpq, D_alpha_used) if alpha != 0.0 else None
K_beta = _build_k_from_lpq(Lpq, D_beta_used) if alpha != 0.0 else None
e_exx = 0.0
if alpha != 0.0 and exxdiv is PBCExxDiv.EWALD:
K_alpha = apply_exxdiv_ewald_to_K([K_alpha], [S], [D_alpha_used], madelung)[
0
]
K_beta = apply_exxdiv_ewald_to_K([K_beta], [S], [D_beta_used], madelung)[0]
e_exx = exxdiv_ewald_energy_shift(
[D_alpha_used], [S], madelung, hf_exchange_fraction=alpha, weights=[1.0]
) + exxdiv_ewald_energy_shift(
[D_beta_used], [S], madelung, hf_exchange_fraction=alpha, weights=[1.0]
)
F_HF_alpha = J - alpha * K_alpha if K_alpha is not None else J
F_HF_beta = J - alpha * K_beta if K_beta is not None else J
E_xc = 0.0
V_xc_alpha = V_xc_beta = 0.0
if is_ks:
_set_xc_density(D_alpha_set, D_alpha_used)
_set_xc_density(D_beta_set, D_beta_used)
xc = build_xc_periodic_uks(
basis, system, grid, func, D_alpha_set, D_beta_set, lat_opts
)
V_xc_alpha = np.real(bloch_sum(xc.V_alpha, k_gamma))
V_xc_beta = np.real(bloch_sum(xc.V_beta, k_gamma))
V_xc_alpha = 0.5 * (V_xc_alpha + V_xc_alpha.T)
V_xc_beta = 0.5 * (V_xc_beta + V_xc_beta.T)
E_xc = float(xc.e_xc)
F_alpha = Hcore + F_HF_alpha + V_xc_alpha
F_beta = Hcore + F_HF_beta + V_xc_beta
F_alpha = 0.5 * (F_alpha + F_alpha.T)
F_beta = 0.5 * (F_beta + F_beta.T)
E_core = float(np.einsum("ij,ij->", D_alpha_used + D_beta_used, Hcore))
E_J = 0.5 * float(np.einsum("ij,ij->", D_alpha_used + D_beta_used, J))
E_K = 0.0
if K_alpha is not None:
E_K = (
-0.5
* alpha
* (
float(np.einsum("ij,ij->", D_alpha_used, K_alpha))
+ float(np.einsum("ij,ij->", D_beta_used, K_beta))
)
)
E_elec = E_core + E_J + E_K + E_xc
E_total = E_elec + float(e_nuc)
FDS_a = F_alpha @ D_alpha_used @ S
FDS_b = F_beta @ D_beta_used @ S
grad_a = FDS_a - FDS_a.T
grad_b = FDS_b - FDS_b.T
grad_norm = float(
np.sqrt(np.linalg.norm(grad_a) ** 2 + np.linalg.norm(grad_b) ** 2)
)
dE = E_total - E_prev
converged = (
iter_idx > 1
and abs(dE) < float(opts.conv_tol_energy)
and grad_norm < float(opts.conv_tol_grad)
)
scf_trace.append(
SCFIteration(
iter=iter_idx,
energy=float(E_total),
delta_e=float(dE if iter_idx > 1 else 0.0),
grad_norm=float(grad_norm),
diis_subspace=(accel.subspace_size if accel is not None else 0),
)
)
plog.iteration(
iter_idx,
energy=float(E_total),
dE=float(dE if iter_idx > 1 else 0.0),
grad=float(grad_norm),
diis=(accel.subspace_size if accel is not None else 0),
)
result.energy = E_total
result.e_electronic = E_elec
result.e_coulomb = E_J
result.e_hf_exchange = E_K
result.e_xc = E_xc
result.e_exxdiv = e_exx
result.n_iter = iter_idx
result.mo_energies_alpha = eps_alpha_f
result.mo_coeffs_alpha = C_alpha_f
result.density_alpha = D_alpha_used
result.fock_alpha = F_alpha
result.mo_energies_beta = eps_beta_f
result.mo_coeffs_beta = C_beta_f
result.density_beta = D_beta_used
result.fock_beta = F_beta
result.smearing_temperature = smear_T
result.entropy = float(entropy_cur)
result.free_energy = E_total - smear_T * float(entropy_cur)
result.fermi_level_alpha = float(fermi_a)
result.fermi_level_beta = float(fermi_b)
result.occupations_alpha = np.asarray(occ_a, dtype=float)
result.occupations_beta = np.asarray(occ_b, dtype=float)
if converged:
result.converged = True
result.s_squared = _gamma_open_shell_s_squared(
n_alpha, n_beta, C_alpha_f, C_beta_f, S, occ_a, occ_b, smear_opts
)
if check_energy_sanity:
_check_energy_sanity(
result,
system,
plog,
driver="run_pbc_gdf_uks",
raise_if_converged=True,
)
plog.converged(n_iter=iter_idx, energy=E_total, converged=True)
return result
# Skipped entirely (not even recorded) while the PATTERN_HOLD window
# is active; see the hold_active note at the loop head.
if accel is not None and not hold_active:
F_alpha_ex, F_beta_ex = accel.extrapolate_uhf(
F_alpha,
F_beta,
error_alpha=grad_a,
error_beta=grad_b,
density_alpha=D_alpha_used,
density_beta=D_beta_used,
energy=E_total,
mo_coeffs_alpha=C_alpha_f,
mo_coeffs_beta=C_beta_f,
mo_energies_alpha=eps_alpha_f,
mo_energies_beta=eps_beta_f,
n_alpha=n_alpha,
n_beta=n_beta,
)
if diis_active:
F_alpha, F_beta = F_alpha_ex, F_beta_ex
if fock_mixing_value != 0.0:
if F_alpha_prev_mixed is not None and F_beta_prev_mixed is not None:
F_alpha_mix = (
(1.0 - fock_mixing_value) * F_alpha
+ fock_mixing_value * F_alpha_prev_mixed
)
F_beta_mix = (
(1.0 - fock_mixing_value) * F_beta
+ fock_mixing_value * F_beta_prev_mixed
)
F_alpha = 0.5 * (F_alpha_mix + F_alpha_mix.T)
F_beta = 0.5 * (F_beta_mix + F_beta_mix.T)
F_alpha_prev_mixed = F_alpha.copy()
F_beta_prev_mixed = F_beta.copy()
C_alpha_f, eps_alpha_f = diagonalise(F_alpha)
C_beta_f, eps_beta_f = diagonalise(F_beta)
if _pattern_hold and 1 < iter_idx <= _spinlock_iters:
if n_alpha > 0 and _ca_occ_prev is not None:
C_alpha_f, eps_alpha_f = _mom_reorder(
C_alpha_f, eps_alpha_f, S, _ca_occ_prev, n_alpha
)
if n_beta > 0 and _cb_occ_prev is not None:
C_beta_f, eps_beta_f = _mom_reorder(
C_beta_f, eps_beta_f, S, _cb_occ_prev, n_beta
)
if _pattern_hold and iter_idx <= _spinlock_iters:
_ca_occ_prev = (
np.asarray(C_alpha_f[:, :n_alpha]).copy() if n_alpha > 0 else None
)
_cb_occ_prev = (
np.asarray(C_beta_f[:, :n_beta]).copy() if n_beta > 0 else None
)
D_alpha_prev, D_beta_prev = D_alpha_used, D_beta_used
(D_alpha, D_beta, occ_a, occ_b, fermi_a, fermi_b, entropy_cur) = (
_open_shell_gamma_occupy(
C_alpha_f,
eps_alpha_f,
C_beta_f,
eps_beta_f,
n_alpha,
n_beta,
smear_opts,
)
)
if damper is not None:
damper.update(E_total)
E_prev = E_total
result.s_squared = float(
_gamma_open_shell_s_squared(
n_alpha, n_beta, C_alpha_f, C_beta_f, S, occ_a, occ_b, smear_opts
)
)
result.converged = False
if check_energy_sanity:
_check_energy_sanity(
result,
system,
plog,
driver="run_pbc_gdf_uks",
raise_if_converged=True,
)
plog.converged(n_iter=result.n_iter, energy=result.energy, converged=False)
return result
# Positive-energy slack for the Γ-GDF energy-sanity guard (Ha). A bound
# neutral cell has E_total < 0, but cramped/artificial Bravais smoke-test
# cells can converge just above zero (the 8-bohr hexagonal H₂ coverage cell
# lands at +0.067 Ha); a positive energy beyond this slack is unphysical.
# Same convention as vibeqc.periodic_k_gdf.POSITIVE_E_SLACK_HA.
POSITIVE_E_SLACK_HA = 1.0
def _check_energy_sanity(
result: PBCGDFResult,
system: PeriodicSystem,
plog: ProgressLogger,
*,
driver: str = "run_pbc_gdf_rhf",
raise_if_converged: bool = False,
) -> None:
"""Post-condition: SCF total energy is in a physically sensible
range relative to the unit cell's atomic content.
Two failure signatures are rejected, mirroring the multi-k guard
(:func:`vibeqc.periodic_k_gdf._check_energy_sanity`):
* **Runaway** -- ``|E_total|`` beyond ``max(10.S Z^2, 100)`` Ha (the
loose hydrogenic bound on absolute binding, factor-10 slack).
Seen on tight ionic crystals when the Γ-only compcell Hartree is
inconsistent (the LiH FCC primitive cell converges to +579.8 Ha
at RHF/UHF and +1172.6 Ha at UKS-PBE vs PySCF's -8.24 Ha) or the
AFT convention is mismatched (1e4-1e6 Ha range) -- the SCF
"converges" to a numerical fixed point of the broken Fock that
has no physical meaning.
* **Unbound** -- ``E_total`` positive beyond
``POSITIVE_E_SLACK_HA``. A bound neutral cell has
``E_total < 0``; small positive energies are tolerated because
cramped/artificial Bravais smoke-test cells can legitimately
converge just above zero.
Callers with pinned warn-and-tag behaviour (``run_pbc_gdf_rhf``;
see ``test_pbc_gdf_rhf_production_defaults_lih_ionic``) keep the
default ``raise_if_converged=False``: warn loudly + tag the
backend ``+SANITY_FAILED``, but hand the value back for parity
diagnostics. The open-shell/KS drivers pass
``raise_if_converged=True``: a non-physical energy reported as
*converged* is the silent-corruption pattern CLAUDE.md §7 warns
about, so they RAISE instead of returning it (bypass with
``check_energy_sanity=False``). A non-converged run warns + tags
either way (``converged=False`` already signals failure).
"""
E = float(result.energy)
z_sum_sq = sum(atom.Z**2 for atom in system.unit_cell)
sane_bound = max(10.0 * z_sum_sq, 100.0) # loose hydrogenic + floor
runaway = abs(E) > sane_bound
unbound = E > POSITIVE_E_SLACK_HA
if not (runaway or unbound):
return
reason = "runaway divergence" if runaway else "positive (unbound) total energy"
msg = (
f"{driver}: SCF total energy {E:.6e} Ha is non-physical "
f"({reason}). A bound, neutral cell has E_total < 0 (a positive "
f"energy beyond {POSITIVE_E_SLACK_HA:g} Ha slack is unbound) and "
f"|E_total| < {sane_bound:.2e} Ha (loose hydrogenic bound on "
"S_atoms Z^2/2). The SCF has converged to a numerical fixed point "
"of a broken Fock (CLAUDE.md Sec.7 -- not a convergence-aid "
"problem). Known trigger: Gamma-only gdf_method='compcell' on a "
"tight ionic cell, where the Gamma-only Hartree (G=0 dropped) "
"cannot resolve overlapping AO-pair images -- use the multi-k "
"route (run_krhf/run_kuhf/run_kuks_periodic_gdf, validated at uHa "
"parity on LiH FCC at kmesh=(2,2,2)) or gdf_method='rsgdf' at "
"Gamma. Pass check_energy_sanity=False to bypass this guard "
"(diagnostics only)."
)
try:
result.backend = result.backend + "+SANITY_FAILED"
except Exception:
pass
if raise_if_converged and result.converged:
raise RuntimeError(msg)
plog.info(" WARNING: " + msg)