"""Gamma-only periodic RHF/RKS with native Gaussian density fitting.
This is the first self-hosted GDF SCF driver in vibe-qc's periodic
stack. It uses the native C++ periodic 2c/3c auxiliary-integral kernels
through :func:`vibeqc.aux_basis.build_lpq_native`, then contracts the
resulting ``Lpq`` tensor in Python to build J and K.
Scope
-----
* Γ-only RHF and closed-shell RKS.
* Native vibe-qc integrals only; PySCF and CRYSTAL remain out-of-process
validation references.
* Multi-k GDF is still pending.
The implementation is deliberately explicit and chatty in the progress
log because periodic GDF setup can spend noticeable wall time in the
2c/3c lattice sums before the first SCF row appears.
"""
from __future__ import annotations
from dataclasses import dataclass, field
from typing import List, Optional, Tuple, Union
import numpy as np
from ._vibeqc_core import (
BasisSet,
CoulombMethod,
Functional,
GridOptions,
InitialGuess,
LatticeMatrixSet,
LatticeSumOptions,
PeriodicKSOptions,
PeriodicRHFOptions,
PeriodicSystem,
SCFIteration,
bloch_sum,
build_grid,
build_jk_gamma_molecular_limit,
build_xc_periodic,
compute_kinetic_lattice,
compute_overlap_lattice,
direct_lattice_cells,
nuclear_repulsion_per_cell,
)
from .aux_basis import build_lpq_native, default_aux_for, make_aux_basis_set
from .ewald_composed import build_j_ewald_3d
from .ewald_j import auto_grid
from .guess import initial_density_closed_shell
from .linear_dependence import scf_preflight_overlap_check
from .occupations import (
aufbau_occupations_per_k as _aufbau_occupations_per_k,
)
from .occupations import (
fermi_dirac_occupations_per_k as _fermi_dirac_occupations_per_k,
)
from .occupations import (
hartree_to_kelvin_temperature as _hartree_to_kelvin_temperature,
)
from .options_dump import dump_active_settings
from .periodic_grid import build_periodic_becke_grid
from .periodic_rhf_ewald import (
_canonical_orthogonalizer,
_PulayDIIS,
_reject_unsupported_python_accelerator,
)
from .progress import ProgressLogger, resolve_progress
from .symmetry_integrals_reduced import (
compute_kinetic_lattice_reduced,
compute_overlap_lattice_reduced,
)
__all__ = [
"PeriodicRHFGDFResult",
"run_rhf_periodic_gamma_gdf",
"run_rks_periodic_gamma_gdf",
]
[docs]
@dataclass
class PeriodicRHFGDFResult:
"""Result of :func:`run_rhf_periodic_gamma_gdf`.
The ``backend`` field indicates the actual J/K mechanism:
- ``"native-gamma-gdf"`` — true density fitting via Lpq
(1D/2D systems; 3D molecular-limit cells whose GDF cutoff
includes only the home cell).
- ``"ewald-jk-fallback"`` — Ewald-3D J + real-space K
(3D systems with image cells in the GDF cutoff; see
:mod:`vibeqc.pbc_gdf` for the compensated-cell GDF path)."""
energy: float
e_electronic: float
e_nuclear: float
n_iter: int
converged: bool
mo_energies: np.ndarray
mo_coeffs: np.ndarray
density: np.ndarray
fock: np.ndarray
overlap: np.ndarray
scf_trace: List[SCFIteration] = field(default_factory=list)
hcore: np.ndarray = field(default_factory=lambda: np.empty((0, 0)))
aux_basis_name: str = ""
n_aux: int = 0
n_fit: int = 0
linear_dep_threshold: float = 1e-9
gdf_algorithm: str = "bare"
e_xc: float = 0.0
e_coulomb: float = 0.0
e_hf_exchange: float = 0.0
functional: str = ""
fock_mixing: float = 0.0
level_shift: float = 0.0
level_shift_warmup_cycles: int = 0
smearing_temperature: float = 0.0
fermi_level: float = 0.0
entropy: float = 0.0
free_energy: float = 0.0
occupations: np.ndarray = field(default_factory=lambda: np.empty(0))
backend: str = "native-gamma-gdf"
# --- Ewald-J build parameters (Step 2 of the 2026-05-13 gauge fix) ---
# Defaults match `periodic_rks_ewald.py` (`omega = 0.5 / bohr`, FFT
# Poisson spacing `0.3 bohr`); the composed Ewald-3D J is ω-invariant
# up to numerical precision (validated in `validate_ewald_identity`), so
# tuning matters only for cost. Exposed here as constants so the choice
# is discoverable; can be promoted to function-level kwargs in a follow-up
# if user-tuning becomes needed.
_J_EWALD_OMEGA: float = 0.0
_J_EWALD_SPACING_BOHR: float = 0.3
def _gauge_lat_opts_for_v_ne_and_e_nuc(
src: LatticeSumOptions,
system: PeriodicSystem,
) -> LatticeSumOptions:
"""Return a clone of ``src`` with the Coulomb gauge forced to Ewald-3D
for the V_ne lattice sum and the nuclear-nuclear lattice sum, when
the system is fully 3D-periodic.
Rationale (gauge-regression fix, 2026-05-13):
The native gamma GDF driver's three Coulomb-bearing pieces
(V_ne, e_nuc, J / K via Lpq) must share a single gauge for the
total energy to be physically meaningful. PySCF
(``exxdiv='ewald'``), CRYSTAL14, and vibe-qc's own ``EWALD_3D``
direct path all use the Madelung/Ewald gauge — that's the
convention the now-sealed v0.8.0 STO-3G parity baseline
(``examples/regression/crystal_parity/baseline_sto3g/PARITY_TABLE.md``)
is calibrated against.
The user-supplied ``lat_opts.coulomb_method`` controls **J/K**
routing (DIRECT_TRUNCATED vs EWALD_3D vs DF backends), not
one-electron / nuclear gauge — those need to be Ewald-3D
unconditionally on 3D-periodic systems.
For ``dim < 3`` (vacuum-padded slab / wire) the bare lattice
sum is correct, so the passthrough avoids a no-op detour
through Ewald-3D which is 3D-only on the C++ side anyway.
See ``docs/handover_periodic_scf_gauge_regression_2026-05-13.md``
§ "Fix path — three steps" / Step 1.
"""
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:
"""Closed-shell Coulomb matrix from a fitted periodic ERI factor."""
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:
"""Closed-shell exchange matrix from a fitted periodic ERI factor.
``D`` is the total RHF density, including the factor of two for
doubly occupied orbitals. The RHF Fock then uses ``J - 0.5 K``.
"""
K = np.einsum("Lmk,kl,Lnl->mn", Lpq, D, Lpq, optimize=True)
return 0.5 * (K + K.T)
def _density_set_gamma(
template: LatticeMatrixSet,
D: np.ndarray,
) -> LatticeMatrixSet:
"""Return a Γ-only density set with ``D`` in the home-cell block."""
if D.shape[0] != template.nbf:
raise ValueError(
f"_density_set_gamma: D has nbf={D.shape[0]} but template "
f"has nbf={template.nbf}"
)
zero = np.zeros_like(D)
for i in range(len(template)):
template.set_block(i, D if i == 0 else zero)
return template
def _resolve_fock_mixing(options, override: Optional[float]) -> float:
"""Return previous-Fock weight for CRYSTAL-style FMIXING."""
value = override
if value is None:
value = getattr(options, "fock_mixing", None)
if value is None:
percent = getattr(options, "fmixing_percent", None)
if percent is not None:
value = float(percent) / 100.0
if value is None:
value = 0.0
value = float(value)
if not (0.0 <= value < 1.0):
raise ValueError(
f"run_rhf_periodic_gamma_gdf: fock_mixing must be in [0, 1); got {value}"
)
return value
def _resolve_level_shift_warmup_cycles(
options,
*,
level_shift: float,
max_iter: int,
override: Optional[int] = None,
) -> int:
"""Return the number of shifted startup cycles.
``None`` or ``-1`` means "auto": if a nonzero level shift is active,
use up to five shifted startup cycles while leaving at least one
final unshifted tail cycle. ``0`` keeps the historic persistent
level-shift behaviour. Positive integers request an explicit warm-up
length and are capped to preserve the unshifted tail cycle.
"""
max_iter = int(max_iter)
if abs(float(level_shift)) == 0.0 or max_iter <= 1:
return 0
raw = override
if raw is None:
raw = getattr(options, "level_shift_warmup_cycles", None)
if raw is None:
raw = -1
raw = int(raw)
if raw < -1:
raise ValueError(
"run_rhf_periodic_gamma_gdf: level_shift_warmup_cycles must "
f"be >= 0 or -1 for auto; got {raw}"
)
if raw < 0:
raw = 5
if raw == 0:
return 0
return min(raw, max_iter - 1)
def _density_from_orbitals_and_occupations(
C: np.ndarray,
occupations: np.ndarray,
) -> np.ndarray:
"""One-particle density for closed-shell occupations in ``[0, 2]``."""
D = (C * np.asarray(occupations, dtype=float)[None, :]) @ C.T
return 0.5 * (D + D.T)
[docs]
def run_rhf_periodic_gamma_gdf(
system: PeriodicSystem,
basis: BasisSet,
options: Optional[Union[PeriodicRHFOptions, PeriodicKSOptions]] = None,
*,
functional: Optional[str] = None,
aux_basis: Optional[str] = None,
aux_drop_eta: float = 0.0,
linear_dep_threshold: float = 1e-7,
gdf_linear_dep_threshold: float = 1e-9,
apply_modrho: bool = True,
gdf_algorithm: str = "bare",
rsgdf_omega: float = 0.4,
fock_mixing: Optional[float] = None,
level_shift_warmup_cycles: Optional[int] = None,
symmetry_stabilize: bool = False,
symmetry_reduce_fock: bool = False,
progress: Union[bool, ProgressLogger, None] = None,
verbose: Optional[int] = None,
) -> PeriodicRHFGDFResult:
"""Run Γ-only closed-shell periodic RHF/RKS via native GDF.
Parameters mirror the older target API, but the implementation is
now entirely vibe-qc owned. If ``functional`` is provided, the
driver adds native libxc ``V_xc`` on the periodic/molecular Becke
grid and uses the functional's exact-exchange fraction for hybrids.
"""
opts = options if options is not None else PeriodicRHFOptions()
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 = Functional(func_name, 1) if is_ks else None
alpha = float(func.hf_exchange_fraction) if func is not None else 1.0
fock_mixing_value = _resolve_fock_mixing(opts, fock_mixing)
level_shift = float(getattr(opts, "level_shift", 0.0))
max_iter = int(opts.max_iter)
warmup_cycles = _resolve_level_shift_warmup_cycles(
opts,
level_shift=level_shift,
max_iter=max_iter,
override=level_shift_warmup_cycles,
)
smearing_T = float(getattr(opts, "smearing_temperature", 0.0))
if smearing_T < 0.0:
raise ValueError(
"run_rhf_periodic_gamma_gdf: smearing_temperature must be >= 0"
)
n_elec = system.n_electrons()
if n_elec % 2 != 0:
raise ValueError(
"run_rhf_periodic_gamma_gdf: closed-shell RHF/RKS requires even "
f"electron count; got {n_elec}"
)
if system.multiplicity != 1:
raise ValueError(
"run_rhf_periodic_gamma_gdf: closed-shell RHF/RKS requires "
f"multiplicity=1; got {system.multiplicity}"
)
n_occ = n_elec // 2
aux_name = aux_basis or default_aux_for(basis.name)
label = f"RKS {func_name}" if is_ks else "RHF"
plog.info(
f"{label} Gamma native GDF / "
f"aux={aux_name}, algorithm={gdf_algorithm}, "
f"cutoff={lat_opts.cutoff_bohr:.2f} bohr"
)
plog.info(f"basis: {basis.name} ({basis.nbasis} BFs / {basis.nshells} shells)")
dim = int(system.dim)
active_lengths = [
float(np.linalg.norm(np.asarray(system.lattice, dtype=float)[:, i]))
for i in range(dim)
]
plog.info(
"periodicity: "
f"dim={dim}D, active lengths="
+ ", ".join(f"{x:.3f}" for x in active_lengths)
+ " bohr"
)
n_int_cells = len(direct_lattice_cells(system, lat_opts.cutoff_bohr))
n_nuc_cells = len(direct_lattice_cells(system, lat_opts.nuclear_cutoff_bohr))
plog.info(
"lattice cells: "
f"one-electron/GDF cutoff -> {n_int_cells}, "
f"nuclear cutoff -> {n_nuc_cells}"
)
dump_active_settings(
plog,
[
("PeriodicKSOptions" if is_ks else "PeriodicRHFOptions", opts),
("LatticeSumOptions", lat_opts),
(
"GDF kwargs",
{
"functional": func_name or None,
"hf_exchange_fraction": alpha,
"fock_mixing": fock_mixing_value,
"fmixing_percent": 100.0 * fock_mixing_value,
"level_shift": level_shift,
"level_shift_warmup_cycles": warmup_cycles,
"smearing_temperature": smearing_T,
"aux_basis": aux_name,
"aux_drop_eta": float(aux_drop_eta),
"linear_dep_threshold": float(linear_dep_threshold),
"gdf_linear_dep_threshold": float(gdf_linear_dep_threshold),
"apply_modrho": bool(apply_modrho),
"gdf_algorithm": gdf_algorithm,
"rsgdf_omega": float(rsgdf_omega),
},
),
],
)
# ---- Functional + grid --------------------------------------------
grid = None
if is_ks:
grid_options = getattr(opts, "grid", None)
if grid_options is None:
grid_options = 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)),
)
else:
grid = build_grid(system.unit_cell_molecule(), grid_options)
molecular_limit_gdf = int(system.dim) == 3 and n_int_cells == 1
use_ewald_jk = int(system.dim) == 3 and not molecular_limit_gdf
# ---- One-electron integrals at Γ ----------------------------------
#
# Force V_ne (and e_nuc below) through the Ewald-3D gauge regardless
# of the user-supplied coulomb_method. See
# ``_gauge_lat_opts_for_v_ne_and_e_nuc`` for rationale + the
# 2026-05-13 gauge-regression handover doc. ``coulomb_method`` on
# the user-facing ``lat_opts`` still controls J/K routing.
gauge_lat_opts = (
_gauge_lat_opts_for_v_ne_and_e_nuc(lat_opts, system)
if use_ewald_jk
else lat_opts
)
if use_ewald_jk 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",
):
_use_sym = system.symmetry is not None and system.symmetry.operations
if _use_sym:
ops = system.symmetry.operations
plog.info(
f"S/T integrals: symmetry-reduced path "
f"(SG {system.symmetry.international_symbol}, "
f"{system.symmetry.order} ops)"
)
_, S_blocks = compute_overlap_lattice_reduced(
basis,
system,
lat_opts,
ops,
)
S_lat = compute_overlap_lattice(basis, system, lat_opts)
for i in range(len(S_lat)):
S_lat.set_block(i, S_blocks[i])
_, T_blocks = compute_kinetic_lattice_reduced(
basis,
system,
lat_opts,
ops,
)
T_lat = compute_kinetic_lattice(basis, system, lat_opts)
for i in range(len(T_lat)):
T_lat.set_block(i, T_blocks[i])
else:
S_lat = compute_overlap_lattice(basis, system, lat_opts)
T_lat = compute_kinetic_lattice(basis, system, lat_opts)
from .periodic_v_ne import compute_nuclear_lattice_dispatch
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_rhf_periodic_gamma_gdf: canonical orthogonalisation "
f"dropped too many directions (n_occ={n_occ}, n_kept={n_kept})"
)
# ---- Native GDF factor --------------------------------------------
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)")
with plog.stage(
"native_lpq",
detail=(
f"2c/3c lattice ERIs, aux={aux.nbasis}, "
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=bool(apply_modrho),
algorithm=gdf_algorithm,
rsgdf_omega=float(rsgdf_omega),
)
plog.info(
f"Lpq: {Lpq.shape[0]} fit vectors, "
f"shape=({Lpq.shape[0]}, {Lpq.shape[1]}, {Lpq.shape[2]})"
)
e_nuc = nuclear_repulsion_per_cell(system, gauge_lat_opts)
# --- Steps 2 + 3 of the 2026-05-13 gauge fix: Ewald-3D J + K -------
#
# Step 2: for dim == 3 the Hartree J is built via the composed
# Ewald-3D builder (short-range erfc/r + long-range FFT-Poisson of
# erf/r), bringing J into the same Madelung/Ewald gauge as V_ne /
# e_nuc (Step 1) and CRYSTAL14 / PySCF `exxdiv='ewald'`.
#
# Step 3 (revised): the K build *also* moves off the native modrho-
# Lpq factor for dim == 3. The 2026-05-14 forensic K-vs-K comparison
# on LiH/STO-3G showed `_build_k_from_lpq` gives ||K|| ~5x too large
# on tight ionic cells — `build_lpq_native`'s native Lpq tensor is
# wrong for tight cells where AO images overlap. So for v0.8.0 K is
# built via the same full-range real-space path the working EWALD_3D
# driver uses (`build_jk_gamma_molecular_limit` at omega=0).
#
# The 2026-05-14 handover originally framed this molecular-limit
# kernel as a builder bug ("misses g_λ ≠ g_σ contributions"); the
# 2026-05-15 cutoff-convergence sweep (see
# `docs/handover_periodic_jk_revert_2026-05-15.md`) showed the
# alternative "homogeneous-D triple-cell-sum" build diverges with
# lattice cutoff and the molecular-limit kernel is in fact the
# correct Γ-only periodic Fock build (the inverse-Bloch convention
# P(g) = P_Γ · δ_{g,0}). The −0.9 Ha gap to CRYSTAL14 SHRINK 8 8
# on LiH primitive is BZ-sampling residual, not a builder bug —
# multi-k sampling is the fix for tight ionic crystals.
#
# Consequence: in v0.8.0 this "GDF" driver does NOT use density
# fitting for J or K on tight dim == 3 cells — it is Ewald-J +
# real-space-K, identical in operation to
# `run_rhf_periodic_gamma_ewald3d`. In the molecular-limit regime
# where the GDF cutoff contains only the home cell, the native Lpq
# tensor remains the intended backend and exercises arbitrary 3D
# lattice metrics without the tight-cell image-overlap pathology.
#
# Pre-compute the FFT grid once; it depends only on lattice geometry.
j_ewald_grid_shape: Optional[Tuple[int, int, int]] = None
if use_ewald_jk:
lat_arr = np.asarray(system.lattice, dtype=float)
j_ewald_grid_shape = tuple(
int(x) for x in auto_grid(lat_arr, _J_EWALD_SPACING_BOHR)
)
plog.info(
"J gauge: Ewald-3D "
f"(omega={'auto' if _J_EWALD_OMEGA <= 0.0 else f'{_J_EWALD_OMEGA:.3f}'} / bohr, "
f"FFT grid {j_ewald_grid_shape[0]}x"
f"{j_ewald_grid_shape[1]}x{j_ewald_grid_shape[2]})"
)
plog.info(
"K gauge: full-range real-space "
"(molecular-limit kernel = correct Γ-only build; "
"tight-ionic accuracy needs multi-k)"
)
elif molecular_limit_gdf:
plog.info(
"J/K backend: native Gamma GDF "
"(3D molecular-limit cell; GDF cutoff contains home cell only)"
)
def _build_j(D_in: np.ndarray) -> np.ndarray:
"""Coulomb J for the current density. For dim==3 use the
Ewald-3D composed builder (gauge-correct); fall back to the
Lpq-based J for dim<3 (vacuum-padded; bare gauge is fine)."""
if use_ewald_jk:
_omega = _J_EWALD_OMEGA
if _omega <= 0.0:
from .bipole_ext_el_pole import crystal_default_ewald_alpha
V_cell = float(
abs(np.linalg.det(np.asarray(system.lattice, dtype=float)))
)
_omega = crystal_default_ewald_alpha(V_cell)
J_out = build_j_ewald_3d(
basis,
system,
D_in,
omega=_omega,
lattice_opts=lat_opts,
grid_shape=j_ewald_grid_shape,
origin=None,
spacing_bohr=_J_EWALD_SPACING_BOHR,
)
else:
J_out = _build_j_from_lpq(Lpq, D_in)
return 0.5 * (J_out + J_out.T)
def _build_k(D_in: np.ndarray) -> np.ndarray:
"""Exchange K for the current density. For dim==3 use the
full-range real-space molecular-limit builder (the correct
Γ-only build per the 2026-05-15 cutoff-convergence finding);
fall back to the Lpq-based K for dim<3 (vacuum-padded — the
native Lpq tensor is fine when AO images don't overlap)."""
if use_ewald_jk:
jk = build_jk_gamma_molecular_limit(
basis,
system,
lat_opts,
D_in,
0.0,
)
K_out = np.asarray(jk.K)
else:
K_out = _build_k_from_lpq(Lpq, D_in)
return 0.5 * (K_out + K_out.T)
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,
) -> Tuple[np.ndarray, float, float]:
if smearing_T <= 0.0:
occ = _aufbau_occupations_per_k([eps], n_occ)[0]
return np.asarray(occ, dtype=float), 0.0, 0.0
occ_list, mu, entropy = _fermi_dirac_occupations_per_k(
[eps],
[1.0],
float(n_elec),
smearing_T,
)
return np.asarray(occ_list[0], dtype=float), float(mu), float(entropy)
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()
occ, fermi_level, entropy = occupations_from_eps(eps0)
if smearing_T > 0.0:
plog.info(
"smearing: Fermi-Dirac kBT = "
f"{smearing_T:.6g} Ha "
f"({_hartree_to_kelvin_temperature(smearing_T):.1f} K)"
)
damping = float(opts.damping)
if not (0.0 <= damping < 1.0):
raise ValueError(
f"run_rhf_periodic_gamma_gdf: damping must be in [0, 1); got {damping}"
)
if fock_mixing_value != 0.0:
plog.info(
"fock mixing: CRYSTAL FMIXING "
f"{100.0 * fock_mixing_value:.1f}% "
"(previous Fock/KS matrix weight)"
)
_reject_unsupported_python_accelerator(opts, "run_rhf_periodic_gamma_gdf")
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
)
if level_shift != 0.0:
if warmup_cycles > 0:
cycle_word = "cycle" if warmup_cycles == 1 else "cycles"
plog.info(
"level-shift warm-up: "
f"{warmup_cycles} {cycle_word} at {level_shift:.3f} Ha; "
"then restart unshifted"
)
else:
plog.info(
f"level shift: {level_shift:.3f} Ha applied at each diagonalization"
)
D_set = compute_overlap_lattice(basis, system, lat_opts) if is_ks else None
plog.banner(f"SCF ({label} Gamma, native GDF)")
plog.info(" iter energy (Ha) dE ||[F,DS]|| DIIS")
scf_trace: List[SCFIteration] = []
result = PeriodicRHFGDFResult(
energy=0.0,
e_electronic=0.0,
e_nuclear=float(e_nuc),
e_xc=0.0,
e_coulomb=0.0,
e_hf_exchange=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,
scf_trace=scf_trace,
hcore=Hcore,
aux_basis_name=aux_name,
n_aux=int(aux.nbasis),
n_fit=int(Lpq.shape[0]),
linear_dep_threshold=float(gdf_linear_dep_threshold),
gdf_algorithm=gdf_algorithm,
functional=func_name,
fock_mixing=fock_mixing_value,
level_shift=level_shift,
level_shift_warmup_cycles=warmup_cycles,
smearing_temperature=smearing_T,
fermi_level=float(fermi_level),
entropy=float(entropy),
free_energy=0.0,
occupations=np.asarray(occ, dtype=float),
backend=("ewald-jk-fallback" if use_ewald_jk else "native-gamma-gdf"),
)
E_prev = 0.0
C_final = C0
eps_final = eps0
F_prev_mixed: Optional[np.ndarray] = None
P_cache: list = []
if symmetry_stabilize and system.symmetry is not None:
from .symmetry_integrals import symmorphic_operations as _so
from .symmetry_scf import build_ao_permutation_cache as _bpc
ops = _so(system.symmetry.operations)
if ops:
P_cache = _bpc(system, basis, ops)
plog.info(f"symmetry_stabilize: {len(P_cache)} symmorphic ops")
for iter_idx in range(1, max_iter + 1):
if warmup_cycles > 0 and iter_idx == warmup_cycles + 1:
if diis is not None:
diis = _PulayDIIS(
max_subspace=int(opts.diis_subspace_size),
)
F_prev_mixed = None
plog.info("restart: unshifted Fock with fresh DIIS history")
active_level_shift = (
level_shift
if (
level_shift != 0.0 and (warmup_cycles == 0 or iter_idx <= warmup_cycles)
)
else 0.0
)
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
)
if symmetry_reduce_fock and system.symmetry is not None:
from .symmetry_fock_reduced import compute_jk_gamma_reduced as _jk_red
J, K = _jk_red(
basis,
system,
lat_opts_2e,
D_used,
system.symmetry.operations,
)
else:
J = _build_j(D_used)
K = _build_k(D_used) if alpha != 0.0 else None
if symmetry_stabilize and P_cache:
from .symmetry_scf import symmetrize_matrix as _sym
J = _sym(J, P_cache)
if K is not None:
K = _sym(K, P_cache)
F_HF_part = J - 0.5 * alpha * K if K is not None else J
E_xc = 0.0
V_xc = 0.0
if is_ks:
_density_set_gamma(D_set, D_used)
xc_contrib = build_xc_periodic(
basis,
system,
grid,
func,
D_set,
lat_opts,
)
V_xc = np.real(bloch_sum(xc_contrib.V_xc, k_gamma))
V_xc = 0.5 * (V_xc + V_xc.T)
E_xc = float(xc_contrib.e_xc)
F = Hcore + F_HF_part + V_xc
F = 0.5 * (F + F.T)
E_core = float(np.einsum("ij,ij->", D_used, Hcore))
E_hf_trace = 0.5 * float(np.einsum("ij,ij->", D_used, F_HF_part))
E_elec = E_xc + E_core + E_hf_trace
E_total = E_elec + float(e_nuc)
E_coulomb = 0.5 * float(np.einsum("ij,ij->", D_used, J))
E_hf_K = (
-0.25 * alpha * float(np.einsum("ij,ij->", D_used, K))
if K is not None
else 0.0
)
free_energy = E_total - smearing_T * entropy
FDS = F @ D_used @ S
grad = FDS - FDS.T
grad_norm = float(np.linalg.norm(grad))
dE = free_energy - E_prev
converged = (
iter_idx > 1
and (warmup_cycles == 0 or iter_idx > warmup_cycles)
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(free_energy),
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(free_energy),
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),
)
plog.energy_decomposition(
iter_idx,
E_kin=float(np.einsum("ij,ij->", D_used, T)),
E_ne=float(np.einsum("ij,ij->", D_used, V)),
E_J=E_coulomb,
E_xc=E_xc if is_ks else None,
E_K=E_hf_K,
E_nuc=float(e_nuc),
)
if diis is not None:
F_ex = diis.extrapolate(F, grad)
if diis_active:
F = F_ex
if fock_mixing_value != 0.0:
if F_prev_mixed is not None:
F_mixed = (
1.0 - fock_mixing_value
) * F + fock_mixing_value * F_prev_mixed
F = 0.5 * (F_mixed + F_mixed.T)
F_prev_mixed = F.copy()
if active_level_shift != 0.0:
F_diag = (
F
+ active_level_shift * S
- (active_level_shift / 2.0) * (S @ D_used @ S)
)
F_diag = 0.5 * (F_diag + F_diag.T)
else:
F_diag = F
C_new, eps_new = diagonalise(F_diag)
occ, fermi_level, entropy = 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 = free_energy
result.energy = E_total
result.e_electronic = E_elec
result.e_xc = E_xc
result.e_coulomb = E_coulomb
result.e_hf_exchange = E_hf_K
result.n_iter = iter_idx
result.mo_energies = eps_new
result.mo_coeffs = C_new
result.density = D_used
result.fock = F
result.fermi_level = float(fermi_level)
result.entropy = float(entropy)
result.free_energy = float(free_energy)
result.occupations = np.asarray(occ, dtype=float)
if converged:
J_f = _build_j(D)
K_f = _build_k(D) if alpha != 0.0 else None
F_HF_f = J_f - 0.5 * alpha * K_f if K_f is not None else J_f
E_xc_f = 0.0
V_xc_f = 0.0
if is_ks:
_density_set_gamma(D_set, D)
xc_f = build_xc_periodic(
basis,
system,
grid,
func,
D_set,
lat_opts,
)
V_xc_f = np.real(bloch_sum(xc_f.V_xc, k_gamma))
V_xc_f = 0.5 * (V_xc_f + V_xc_f.T)
E_xc_f = float(xc_f.e_xc)
F_f = 0.5 * ((Hcore + F_HF_f + V_xc_f) + (Hcore + F_HF_f + V_xc_f).T)
C_f, eps_f = diagonalise(F_f)
occ_f, fermi_level_f, entropy_f = occupations_from_eps(eps_f)
E_core_f = float(np.einsum("ij,ij->", D, Hcore))
E_hf_trace_f = 0.5 * float(np.einsum("ij,ij->", D, F_HF_f))
E_elec_f = E_xc_f + E_core_f + E_hf_trace_f
result.energy = E_elec_f + float(e_nuc)
result.e_electronic = E_elec_f
result.e_xc = E_xc_f
result.e_coulomb = 0.5 * float(np.einsum("ij,ij->", D, J_f))
result.e_hf_exchange = (
-0.25 * alpha * float(np.einsum("ij,ij->", D, K_f))
if K_f is not None
else 0.0
)
result.mo_energies = eps_f
result.mo_coeffs = C_f
result.density = D
result.fock = F_f
result.fermi_level = float(fermi_level_f)
result.entropy = float(entropy_f)
result.free_energy = float(result.energy - smearing_T * entropy_f)
result.occupations = np.asarray(occ_f, dtype=float)
result.converged = True
# Diagnostic: HOMO-LUMO gap check at convergence. Mirrors
# CRYSTAL14's `POSSIBLY CONDUCTING STATE - EFERMI(AU)`
# warning that fires when the gap between the highest
# occupied and lowest unoccupied MO is below threshold.
# For a closed-shell Γ-only calc on a true insulator the
# gap should be tens of mHa or larger; below ~1 mHa the
# system is effectively metallic at this k-point and
# SCF convergence may be fragile / sensitive to start
# guess. This is purely informational — we don't fail
# the SCF. The 1e-3 Ha (~27 meV) threshold matches
# CRYSTAL14's default reporting cutoff.
_CONDUCTING_GAP_THRESHOLD_HA = 1e-3
if n_occ > 0 and n_occ < eps_f.shape[0]:
homo = float(eps_f[n_occ - 1])
lumo = float(eps_f[n_occ])
gap = lumo - homo
if gap < _CONDUCTING_GAP_THRESHOLD_HA:
plog.info(
"POSSIBLY CONDUCTING STATE at convergence: "
f"HOMO-LUMO gap = {gap:.6f} Ha "
f"({gap * 27.2114:.3f} eV) "
f"below the {_CONDUCTING_GAP_THRESHOLD_HA:.0e} Ha threshold "
f"(HOMO={homo:.6f}, LUMO={lumo:.6f}). "
"For metals / small-gap systems consider Fermi smearing "
"(opts.smearing_temperature > 0) or denser k-mesh; "
"Γ-only RHF on a true metal can be ill-defined."
)
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
plog.converged(
n_iter=result.n_iter,
energy=result.energy,
converged=False,
)
return result
# Closed-shell RKS via GDF — the RHF/GDF driver already dispatches to
# the KS branch when ``options`` is :class:`PeriodicKSOptions` (or
# ``functional=`` is set); this alias gives the call site the right
# name for what it does.
run_rks_periodic_gamma_gdf = run_rhf_periodic_gamma_gdf