"""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.
.. warning:: **AFT correction not yet applied — η is a tuning knob.**
The AFT long-range correction (PySCF
``_CCGDFBuilder.get_2c2e`` j2c_p subtraction) makes the
compcell answer η-independent and conditions the metric for
tight ionic systems. Without it:
* **H2 / vacuum-box / loose-aux scenarios**: mHa-scale PySCF
parity at η ≈ 0.25-0.3. Tunable. Functional.
* **Tight ionic crystals (LiH, MgO, NaCl, …)**: the compensated
metric is numerically singular; SCF diverges. Use
:func:`run_rhf_periodic_gamma_ewald3d` (the current
production path for ionic crystals) until the AFT correction
ships.
AFT correction is the next milestone for full µHa PySCF parity.
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
from dataclasses import dataclass, field
from typing import List, 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_native,
default_aux_for,
make_aux_basis_set,
)
from .guess import 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 _PulayDIIS, _canonical_orthogonalizer
from .periodic_v_ne import compute_nuclear_lattice_dispatch
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
``ξ = α_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`."""
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"
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 _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
def _build_j_from_lpq(Lpq: np.ndarray, D: np.ndarray) -> np.ndarray:
rho = np.einsum("Lij,ij->L", Lpq, D, optimize=True)
J = np.einsum("L,Lij->ij", rho, Lpq, optimize=True)
return 0.5 * (J + J.T)
def _build_k_from_lpq(Lpq: np.ndarray, D: np.ndarray) -> np.ndarray:
K = np.einsum("Lmk,kl,Lnl->mn", Lpq, D, Lpq, optimize=True)
return 0.5 * (K + K.T)
def run_pbc_gdf_rhf(
system: PeriodicSystem,
basis: BasisSet,
options: Optional[PeriodicRHFOptions] = 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 = 0.2,
apply_aft_correction: bool = False,
aft_precision: float = 1e-10,
aft_ft_convention: str = "libcint",
rsgdf_omega: float = 0.4,
rsgdf_g_precision: float = 1e-10,
rcut_strategy: Optional[object] = None,
rcut_precision: float = 1e-8,
linear_dep_threshold: float = 1e-7,
gdf_linear_dep_threshold: float = 1e-9,
progress: Union[bool, ProgressLogger, None] = None,
verbose: Optional[int] = None,
) -> PBCGDFResult:
"""Closed-shell periodic RHF 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``.
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
0.2). See :func:`vibeqc.aux_basis.build_lpq_compcell`.
Only used when ``gdf_method='compcell'``.
rsgdf_omega
Range-separation parameter ω (bohr⁻¹) for the
``gdf_method='rsgdf'`` path. Default ``0.4``.
SR + LR are independent of ω 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).
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)
opts = options if options is not None else PeriodicRHFOptions()
lat_opts: LatticeSumOptions = opts.lattice_opts
plog = resolve_progress(progress, verbose=verbose)
# 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."
)
n_occ = n_elec // 2
aux_name = aux_basis or default_aux_for(basis.name)
plog.info(
f"PBC-GDF RHF / aux={aux_name}, exxdiv={exxdiv.value}, "
f"eta={compcell_eta:g}, kmesh={kmesh}, "
f"cutoff={lat_opts.cutoff_bohr:.2f} bohr"
)
plog.info(
f"basis: {basis.name} "
f"({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 at Γ ----------------------------------
#
# 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)
V_lat = compute_nuclear_lattice_dispatch(basis, system, gauge_lat_opts)
k_gamma = np.zeros(3)
S = np.real(bloch_sum(S_lat, k_gamma))
T = np.real(bloch_sum(T_lat, k_gamma))
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)
if n_occ > n_kept:
raise RuntimeError(
"run_pbc_gdf_rhf: canonical orthogonalisation dropped too "
f"many directions (n_occ={n_occ}, n_kept={n_kept})"
)
# ---- Compcell Lpq -------------------------------------------------
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} "
f"({aux.nbasis} BFs / {aux.nshells} shells)"
)
if gdf_method == "rsgdf":
with plog.stage(
"rsgdf_lpq",
detail=(
f"rsgdf aux={aux.nbasis}, ω={rsgdf_omega:g}, "
f"fit_thr={gdf_linear_dep_threshold:.1e}"
),
):
Lpq = build_lpq_native(
system, basis, aux,
lat_opts=lat_opts,
linear_dep_thr=float(gdf_linear_dep_threshold),
apply_modrho=True,
algorithm="rsgdf",
rsgdf_omega=float(rsgdf_omega),
rsgdf_g_precision=float(rsgdf_g_precision),
)
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),
)
else:
raise ValueError(
f"run_pbc_gdf_rhf: gdf_method must be 'compcell' or 'rsgdf'; "
f"got {gdf_method!r}"
)
plog.info(
f"Lpq: {Lpq.shape[0]} fit vectors, "
f"shape=({Lpq.shape[0]}, {Lpq.shape[1]}, {Lpq.shape[2]})"
)
# ---- exxdiv setup -------------------------------------------------
madelung = (
madelung_constant_for_cell(system) if exxdiv is PBCExxDiv.EWALD else 0.0
)
if exxdiv is PBCExxDiv.EWALD:
plog.info(
"exxdiv='ewald': K-shift xi = "
f"{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")
# ---- SCF loop -----------------------------------------------------
def diagonalise(F: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
Fp = X.T @ F @ X
Fp = 0.5 * (Fp + Fp.T)
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,
)
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}"
)
use_diis = bool(opts.use_diis)
diis_start_iter = int(opts.diis_start_iter)
diis = _PulayDIIS(
max_subspace=int(opts.diis_subspace_size),
) 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=f"pbc-gdf-{gdf_method}",
)
plog.banner(f"SCF (PBC-GDF {gdf_method})")
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)
# exxdiv='ewald' K-shift (per Γ-point as a 1-element list).
e_exx = 0.0
if exxdiv is PBCExxDiv.EWALD:
K_list = apply_exxdiv_ewald_to_K([K], [S], [D_used], madelung)
K = K_list[0]
e_exx = exxdiv_ewald_energy_shift(
[D_used], [S], madelung,
hf_exchange_fraction=1.0, weights=[1.0],
)
F = Hcore + J - 0.5 * K
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 * float(np.einsum("ij,ij->", D_used, K))
E_elec = E_core + E_J + E_K
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=(diis.subspace_size if diis 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=(diis.subspace_size if diis is not None else 0),
)
if diis is not None:
F_ex = diis.extrapolate(F, grad)
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
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_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)
e_exx_f = 0.0
if exxdiv is PBCExxDiv.EWALD:
K_list = apply_exxdiv_ewald_to_K([K_f], [S], [D], madelung)
K_f = K_list[0]
e_exx_f = exxdiv_ewald_energy_shift(
[D], [S], madelung,
hf_exchange_fraction=1.0, weights=[1.0],
)
F_f = 0.5 * ((Hcore + J_f - 0.5 * K_f)
+ (Hcore + J_f - 0.5 * K_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 * float(np.einsum("ij,ij->", D, K_f))
E_elec_f = E_core_f + E_J_f + E_K_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_exxdiv = e_exx_f
result.mo_energies = eps_f
result.mo_coeffs = C_f
result.density = D
result.fock = F_f
result.converged = True
_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 _check_energy_sanity(
result: PBCGDFResult,
system: PeriodicSystem,
plog: ProgressLogger,
) -> None:
"""Post-condition: SCF total energy is in a physically sensible
range relative to the unit cell's atomic content.
For a bound closed-shell system, ``|E_total|`` should be of order
``Σ_atoms Z²/2`` (the loose hydrogenic upper bound on absolute
binding energy) — say a factor of 10 generous slack. Energies
that exceed this by orders of magnitude (1e4–1e6 Ha range) are
the runaway pattern we see on tight ionic crystals when the AFT
convention is mismatched — the SCF "converges" to a numerical
fixed point of the broken Fock that has no physical meaning.
Surface this loudly. Users who don't see the warning might use
the nonsensical result downstream.
"""
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
if abs(result.energy) > sane_bound:
plog.info(
f" WARNING: SCF total energy = {result.energy:.4e} Ha is "
f"orders of magnitude beyond the physically sensible "
f"bound (|E| < {sane_bound:.2e} Ha for this cell's atomic "
"content). This is the runaway-divergence pattern seen on "
"tight ionic crystals with mismatched AFT convention. "
"Verify: (a) apply_aft_correction=False at default for "
"ionic cells, or (b) aft_ft_convention paired correctly "
"with rcut_strategy. Do NOT trust this E_total as physical."
)
# Don't raise — some callers (parity diagnostics, debug scripts)
# legitimately want the nonsense value back for analysis. Just
# tag the result.
try:
result.backend = result.backend + "+SANITY_FAILED"
except Exception:
pass