"""High-level "run a periodic job" -- companion to :func:`run_job`.

Mirrors the molecular ``run_job`` API for periodic SCFs:

  - ``output.out``       text log (banner, system, basis, SCF trace,
                         energies, properties)
  - ``output.system``    runtime manifest (CPU, OS, libs, wall-time)
  - ``output.molden``    Γ-point MOs in MOLDEN format (when MOs exist)
  - ``output.xsf``       SCF density on a primitive-cell grid
                         (when ``write_density=True``)

Usage::

    from vibeqc import PeriodicSystem, BasisSet
    from vibeqc.periodic_runner import run_periodic_job

    sys_p = PeriodicSystem(...)
    basis = BasisSet(sys_p.unit_cell_molecule(), "sto-3g")
    run_periodic_job(
        sys_p, basis,
        method="RHF",
        output="output",
    )

Status: native GDF covers Gamma and multi-k RHF/RKS/UHF/UKS. The
public RIJCOSX route keeps Gamma RHF on the dedicated periodic COSX
driver and routes true multi-k RHF/RKS/UHF/UKS through the native GDF
loop with ``k_exchange="cosx"``. The earlier PySCF-backed GDF spike is
retired because PySCF and CRYSTAL are external reference programs, not
in-process vibe-qc backends. Explicit ``jk_method="fft_poisson"`` is
retired as a public route; ``AUTO`` resolves to native GDF for
closed-shell jobs and BIPOLE for open-shell jobs.
"""

from __future__ import annotations

import os
import time
import warnings
from contextvars import ContextVar
from functools import wraps
from pathlib import Path
from types import SimpleNamespace
from typing import TYPE_CHECKING, Any, List, Optional, Sequence, Tuple, Union

import numpy as np

if TYPE_CHECKING:
    from ._vibeqc_core import BlochKMesh
    from .bands import BandStructure
    from .kpoints import KPoints

from ._vibeqc_core import (
    Atom,
    BasisSet,
    Crystal,
    D3BJParams,
    Functional,
    InitialGuess,
    PeriodicKSOptions,
    PeriodicRHFOptions,
    PeriodicSystem,
    SpaceGroup,
    SpinlockMode,
    attach_symmetry,
    direct_lattice_cells,
    get_num_threads,
    to_primitive,
)
from .banner import banner, enforce_runtime_pin_from_env, library_versions
from .occupations import (
    hartree_to_kelvin_temperature,
    resolve_smearing_temperature,
)
from .output import (
    HeaderlessBlock,
    HeaderlessColumn,
    OutputChannel,
    OutputPlan,
    OutputWriter,
    active_policy,
    dry_run_manifest,
    flush,
    is_dry_run_estimate_requested,
    is_dry_run_requested,
    render_duration,
    render_energy_labeled,
    section_header,
    warn,
    write,
)
from .output.formats.scf_log import write_scf_trace
from .output._errors import (
    OutputFailureKind,
    warn_output_failure,
    warn_writer_failure,
)
from .output.citations import (
    citation_manifest_rows,
    format_references_block,
    load_default_database,
    write_references_block,
)
from .pbc_gdf import (
    _gamma_dense_core_gdf_parity_held,
    run_pbc_gdf_rhf,
    run_pbc_gdf_uhf,
    run_pbc_gdf_uks,
)
from .periodic_aiccm2026dev_b import (
    _density_blocks_per_k,
    _spin_density_blocks_per_k,
    cyclic_lattice_extension,
    inverse_bloch_transform,
    run_aiccm2026dev_b_rhf,
    run_aiccm2026dev_b_rks,
    run_aiccm2026dev_b_uhf,
    run_aiccm2026dev_b_uks,
)
from .periodic_convergence_auto import (
    ConvergenceStrategy,
    KnobResolution,
    classify_periodic_system,
    insulator_smearing_warning,
    resolve_convergence_strategy,
)
from .periodic_jk_method import (
    PeriodicJKMethod,
    describe_jk_method,
    pick_jk_method,
    validate_jk_method,
)
from .periodic_k_gdf import (
    _gamma_kmesh_info,
    run_krhf_periodic_gdf,
    run_krks_periodic_gdf,
    run_kuhf_periodic_gdf,
    run_kuks_periodic_gdf,
)


_PERIODIC_OUTPUT_WRITER: ContextVar[OutputWriter | None] = ContextVar(
    "vibeqc_periodic_output_writer",
    default=None,
)


def _periodic_output_lifecycle(func):
    """Finish or crash the job-scoped writer at the public API boundary."""

    @wraps(func)
    def wrapped(*args, **kwargs):
        token = _PERIODIC_OUTPUT_WRITER.set(None)
        try:
            result = func(*args, **kwargs)
        except Exception:
            writer = _PERIODIC_OUTPUT_WRITER.get()
            if writer is not None and writer.status == "running":
                writer.crash()
            raise
        else:
            writer = _PERIODIC_OUTPUT_WRITER.get()
            if writer is not None and writer.status == "running":
                writer.finish()
            return result
        finally:
            _PERIODIC_OUTPUT_WRITER.reset(token)

    return wrapped


from .periodic_rhf_gdf import (
    PeriodicRHFGDFResult,
    run_rhf_periodic_gamma_gdf,
)
from .progress import ProgressLogger, resolve_progress
from .smearing import SmearingOptions

# Sentinel distinguishing "smearing_temperature not given" (auto strategy
# may fill it) from the documented explicit values 0.0 / None (both mean
# "smearing off, by user choice").
_SMEARING_UNSET = object()

# Sentinel distinguishing "bz_integration not given" from explicit None
# (which means "use smearing").
_BZ_INTEGRATION_UNSET = object()

_BOHR_TO_ANGSTROM = 0.529177210903

__all__ = ["run_periodic_job"]


def _periodic_xc_gradient_dry_run_estimate_bytes(
    system: PeriodicSystem,
    basis: BasisSet,
    *,
    method_upper: str,
    functional: str | None,
    lattice_cutoff_bohr: float | None = None,
) -> int | None:
    """Best-effort dry-run estimate for periodic RKS/UKS gradients."""
    if method_upper not in ("RKS", "UKS"):
        return None

    from .memory import estimate_periodic_xc_gradient

    opts = PeriodicKSOptions()
    grid = opts.grid
    n_atoms = len(system.unit_cell)
    n_grid_points = (
        n_atoms
        * int(getattr(grid, "n_radial", 75))
        * int(getattr(grid, "n_theta", 17))
        * int(getattr(grid, "n_phi", 36))
    )
    # RKS/UKS Ewald paths may enlarge the density cutoff to 18 bohr before
    # calling the analytic gradient. Charging at least that many cells keeps
    # vq placement from seeing the cheaper default 15-bohr stencil.
    cutoff_base = (
        float(lattice_cutoff_bohr)
        if lattice_cutoff_bohr is not None
        else float(getattr(opts.lattice_opts, "cutoff_bohr", 15.0))
    )
    cutoff = max(cutoff_base, 18.0)
    n_cells = len(direct_lattice_cells(system, cutoff))
    spin = 2 if method_upper == "UKS" else 1
    func = Functional(str(functional), spin)
    kind = getattr(getattr(func, "kind", None), "name", None)
    if kind is None:
        kind = str(getattr(func, "kind", ""))

    return estimate_periodic_xc_gradient(
        n_basis=int(basis.nbasis),
        n_atoms=n_atoms,
        n_grid_points=n_grid_points,
        n_cells=n_cells,
        functional_kind=str(kind),
        open_shell=method_upper == "UKS",
    ).total_bytes


def _basis_primitive_count(basis: BasisSet) -> int | None:
    try:
        return int(sum(len(sh.exponents) for sh in basis.shells()))
    except Exception:
        return None


def _periodic_gpw_gapw_estimate(
    system: PeriodicSystem,
    basis: BasisSet,
    *,
    resolved_jk: PeriodicJKMethod,
    method_upper: str,
    functional: str | None,
    cutoff_ha: float,
    kpoints: object,
):
    """Best-effort estimate for GPW/GAPW FFT-grid SCFs."""
    if resolved_jk not in (PeriodicJKMethod.GPW, PeriodicJKMethod.GAPW):
        return None

    from .memory import estimate_periodic_gpw_gapw
    from .periodic_gapw_grid import nx_for_axis

    lattice = np.asarray(system.lattice, dtype=float)
    n_grid_points = 1
    for axis_length in np.linalg.norm(lattice, axis=0):
        n_grid_points *= nx_for_axis(float(axis_length), float(cutoff_ha))

    functional_kind: str | None = None
    if method_upper in ("RKS", "UKS") and functional is not None:
        func = Functional(str(functional), 2 if method_upper == "UKS" else 1)
        kind = getattr(getattr(func, "kind", None), "name", None)
        functional_kind = str(kind if kind is not None else getattr(func, "kind", ""))

    try:
        n_kpoints = _bloch_kmesh_size(_runner_bloch_kmesh(system, kpoints))
    except Exception:
        n_kpoints = 1

    compact_multik = False
    if n_kpoints > 1:
        try:
            from .periodic_gapw_j import _multik_gpw_is_molecular_limit

            compact_multik = not _multik_gpw_is_molecular_limit(system)
        except Exception:
            compact_multik = False

    n_soft_basis: int | None = None
    augmentation_active: bool | None = False
    if resolved_jk == PeriodicJKMethod.GAPW:
        n_soft_basis = int(basis.nbasis)
        augmentation_active = True
        try:
            from .periodic_gapw_augment import softened_basis

            with warnings.catch_warnings():
                warnings.simplefilter("ignore")
                soft_basis = softened_basis(basis, system)
            n_soft_basis = int(soft_basis.nbasis)
            full_prim = _basis_primitive_count(basis)
            soft_prim = _basis_primitive_count(soft_basis)
            if full_prim is not None and soft_prim is not None:
                augmentation_active = soft_prim != full_prim
        except Exception:
            augmentation_active = True

    return estimate_periodic_gpw_gapw(
        n_basis=int(basis.nbasis),
        n_grid_points=n_grid_points,
        route=resolved_jk.value,
        functional_kind=functional_kind,
        open_shell=method_upper in ("UHF", "UKS"),
        n_kpoints=n_kpoints,
        compact_multik=compact_multik,
        n_soft_basis=n_soft_basis,
        n_atoms=len(system.unit_cell),
        augmentation_active=augmentation_active,
    )


def _periodic_gpw_gapw_dry_run_estimate_bytes(
    system: PeriodicSystem,
    basis: BasisSet,
    *,
    resolved_jk: PeriodicJKMethod,
    method_upper: str,
    functional: str | None,
    cutoff_ha: float,
    kpoints: object,
) -> int | None:
    """Best-effort dry-run estimate for GPW/GAPW FFT-grid SCFs."""
    est = _periodic_gpw_gapw_estimate(
        system,
        basis,
        resolved_jk=resolved_jk,
        method_upper=method_upper,
        functional=functional,
        cutoff_ha=cutoff_ha,
        kpoints=kpoints,
    )
    return None if est is None else est.total_bytes


def _periodic_functional_needs_exchange(
    functional: str | None,
    *,
    open_shell: bool,
) -> bool:
    if functional is None:
        return False
    try:
        func = Functional(str(functional), 2 if open_shell else 1)
    except Exception:
        return False
    return bool(getattr(func, "is_hybrid", False)) or bool(
        getattr(func, "is_range_separated", False)
    )


def _periodic_gdf_aux_basis_size(
    system: PeriodicSystem,
    basis: BasisSet,
    aux_basis: str | None,
) -> int:
    if aux_basis:
        try:
            from .aux_basis import make_aux_basis_set

            return int(
                make_aux_basis_set(
                    system.unit_cell_molecule(),
                    aux_name=aux_basis,
                ).nbasis
            )
        except Exception:
            pass
    return max(1, 3 * int(basis.nbasis))


def _periodic_gdf_estimate(
    system: PeriodicSystem,
    basis: BasisSet,
    *,
    resolved_jk: PeriodicJKMethod,
    method_upper: str,
    functional: str | None,
    kpoints: object,
    aux_basis: str | None,
):
    if resolved_jk not in (PeriodicJKMethod.GDF, PeriodicJKMethod.RIJCOSX):
        return None

    from .memory import estimate_periodic_multik_gdf

    try:
        n_kpoints = _bloch_kmesh_size(_runner_bloch_kmesh(system, kpoints))
    except Exception:
        n_kpoints = 1
    n_aux = _periodic_gdf_aux_basis_size(system, basis, aux_basis)
    open_shell = method_upper in ("UHF", "UKS")
    need_k_pairs = method_upper in ("RHF", "UHF")
    if method_upper in ("RKS", "UKS"):
        need_k_pairs = _periodic_functional_needs_exchange(
            functional,
            open_shell=open_shell,
        )
    if resolved_jk == PeriodicJKMethod.RIJCOSX:
        need_k_pairs = True
    label_parts = [method_upper]
    if functional:
        label_parts.append(str(functional))
    label_parts.append(resolved_jk.value)
    route_label = " ".join(label_parts)
    estimate = estimate_periodic_multik_gdf(
        n_basis=int(basis.nbasis),
        n_aux=n_aux,
        n_kpoints=n_kpoints,
        need_k_pairs=need_k_pairs,
        open_shell=open_shell,
    )
    return SimpleNamespace(
        estimate=estimate,
        n_kpoints=n_kpoints,
        route_label=route_label,
    )


# ============================================================
# Helpers -- symmetry reduction
# ============================================================

# ============================================================
# Helpers -- ECP auto-attach
# ============================================================

_POB_SOURCE_DIR_BY_BASIS = {
    "pob-tzvp": "pob-TZVP",
    "pob-tzvp-rev2": "pob-TZVP-rev2",
    "pob-dzvp-rev2": "pob-DZVP-rev2",
}


def _system_atomic_numbers(system) -> set[int]:
    return {int(atom.Z) for atom in system.unit_cell}


def _format_atomic_numbers(zs: Sequence[int]) -> str:
    try:
        from .basis_crystal import _ELEMENT_SYMBOLS as _symbols
    except Exception:
        _symbols = []

    labels: list[str] = []
    for z in sorted(int(z) for z in zs):
        sym = _symbols[z] if 0 <= z < len(_symbols) else f"Z={z}"
        labels.append(f"{sym}(Z={z})")
    return ", ".join(labels)


def _bundled_pob_source_atoms(name: str) -> list:
    """Return bundled per-element POB source records, if available.

    The checked-in source files cover all-electron POB elements (for example
    Ni/O in P16) and let runtime ECP resolution avoid any network access. The
    heavier POB ECP records are not bundled as source files in this tree; those
    still require the legacy Bredow fetcher and therefore fail closed below if
    they cannot be resolved.
    """
    source_dir_name = _POB_SOURCE_DIR_BY_BASIS.get(name.lower())
    if source_dir_name is None:
        return []
    source_dir = (
        Path(__file__).resolve().parent
        / "basis_library"
        / "sources"
        / source_dir_name
    )
    if not source_dir.is_dir():
        return []

    from .basis_crystal import parse_crystal_atom_basis_file

    atoms = []
    for path in sorted(source_dir.iterdir()):
        if not path.is_file() or path.name.startswith(".") or "_" not in path.name:
            continue
        atoms.append(parse_crystal_atom_basis_file(path))
    return atoms


def _resolve_ecp_data(system, basis) -> tuple:
    """Detect pob-TZVP-REV2 inline ECPs and build C++-ready data.

    When ``basis`` is a CRYSTAL-format basis (pob-TZVP-REV2 family)
    with heavy-element ECP blocks embedded in the basis files themselves,
    convert the parsed :class:`CrystalECP` data to the inline-primitive
    format the C++ periodic SCF drivers expect.

    Returns (ecp_primitive_blocks, ecp_home_centers, effective_charges,
    total_ncore) -- all empty/zero only when the requested POB atoms are
    genuinely all-electron. ECP-bearing POB atoms fail closed if the ECP
    records cannot be resolved; paper routes must not silently continue as
    all-electron calculations after missing runtime data.
    """
    name = getattr(basis, "name", "") or ""
    name_key = name.lower()
    if not name_key.startswith("pob"):
        return [], [], [], 0

    requested_zs = _system_atomic_numbers(system)
    local_atoms = _bundled_pob_source_atoms(name_key)
    local_by_z = {int(atom.Z): atom for atom in local_atoms}
    if requested_zs and requested_zs.issubset(local_by_z):
        from .basis_crystal import build_periodic_ecp_data

        return build_periodic_ecp_data(system, local_atoms)

    missing_local_zs = sorted(requested_zs.difference(local_by_z))
    try:
        from .basis_crystal import (
            build_periodic_ecp_data,
            fetch_bredow_basis_sets,
        )

        _, atoms_dict = fetch_bredow_basis_sets(
            names=[name_key], verbose=False, return_atoms=True
        )
        if not atoms_dict:
            raise RuntimeError("Bredow fetcher returned no parsed atom records")

        atom_list = atoms_dict.get(name_key, [])
        if not atom_list:
            raise RuntimeError(
                f"Bredow fetcher returned no records for basis {name!r}"
            )
        by_z = {int(atom.Z): atom for atom in atom_list}
        missing_fetched_zs = sorted(requested_zs.difference(by_z))
        if missing_fetched_zs:
            raise RuntimeError(
                "Bredow records do not include requested element(s): "
                f"{_format_atomic_numbers(missing_fetched_zs)}"
            )
        missing_ecp_zs = [
            z for z in missing_local_zs
            if not (by_z[z].has_ecp and by_z[z].ecp is not None)
        ]
        if missing_ecp_zs:
            raise RuntimeError(
                "requested POB heavy element record(s) did not include inline "
                f"ECP data: {_format_atomic_numbers(missing_ecp_zs)}"
            )

        return build_periodic_ecp_data(system, atom_list)
    except Exception as exc:
        raise RuntimeError(
            f"_resolve_ecp_data: basis {name!r} requires ECP source data for "
            f"{_format_atomic_numbers(missing_local_zs)}, but the data could "
            "not be resolved. The periodic SCF will not continue as an "
            f"all-electron calculation. Underlying error: {exc}"
        ) from exc


def _validate_smearing_dispatch(
    *,
    method: str,
    jk_method: PeriodicJKMethod,
    smearing_temperature: float,
    kpoints: object = None,
) -> None:
    """Fail fast when a selected backend cannot honour smearing.

    ``kpoints`` is the runner's k-mesh argument. The runner exposes
    open-shell smearing through the GDF drivers; the low-level multi-k Ewald
    UHF/UKS drivers also support per-spin smearing, but there is no current
    ``run_periodic_job`` Ewald dispatch route for them.
    """
    if float(smearing_temperature) <= 0.0:
        return
    if jk_method == PeriodicJKMethod.BIPOLE:
        if method in ("RKS", "UHF", "UKS"):
            return
        raise NotImplementedError(
            "run_periodic_job: smearing_temperature > 0 is implemented "
            "for BIPOLE RKS/UHF/UKS, but BIPOLE RHF still requires integer "
            f"occupations. Got method={method!r}."
        )
    if method in ("UHF", "UKS"):
        # Open-shell smearing is wired on the runner's GDF drivers. The
        # exported multi-k Ewald UHF/UKS drivers also have per-spin smearing,
        # but run_periodic_job no longer exposes an Ewald dispatch branch.
        # BIPOLE UHF is accepted in the BIPOLE branch above; BIPOLE RHF
        # remains integer-occupation only.
        if jk_method in (PeriodicJKMethod.GDF, PeriodicJKMethod.RIJCOSX):
            return
        raise NotImplementedError(
            "run_periodic_job: open-shell (UHF/UKS) smearing_temperature > 0 "
            "is wired through run_periodic_job on the GDF drivers "
            "(jk_method='gdf', Gamma or multi-k) and on the multi-k "
            "RIJCOSX/COSX route. The exported multi-k "
            "Ewald UHF/UKS drivers support spin-resolved smearing directly; "
            "BIPOLE UHF supports spin-resolved smearing directly through "
            "this runner; BIPOLE RHF and the other runner routes do not. "
            f"Got method={method!r}, jk_method={jk_method.value!r}, "
            f"kpoints={'set' if kpoints is not None else 'None'}."
        )
    if jk_method == PeriodicJKMethod.RIJCOSX:
        if kpoints is not None:
            return
        raise NotImplementedError(
            "run_periodic_job: RIJCOSX smearing_temperature > 0 is wired "
            "only on the true multi-k COSX route. Pass kpoints= with at "
            "least two k-points, or omit smearing for the Gamma RIJCOSX "
            "RHF driver."
        )
    if jk_method not in (
        PeriodicJKMethod.GDF,
        PeriodicJKMethod.GPW,
        PeriodicJKMethod.GAPW,
    ):
        raise NotImplementedError(
            "run_periodic_job: smearing_temperature > 0 is currently "
            "wired through the GDF, GPW, GAPW, and multi-k RIJCOSX "
            "closed-shell RHF/RKS "
            f"paths; selected J/K method is {jk_method.value!r}."
        )


def _validate_density_mixer_dispatch(
    *,
    density_mixer: Optional[str],
    density_mixer_depth: int,
    density_mixer_beta: float,
    density_mixer_kerker: bool,
    kerker_k0: float,
    kerker_strength: float,
    kerker_cutoff_ha: float,
    jk_method: Optional["PeriodicJKMethod"] = None,
    method: Optional[str] = None,
    kpoints: Optional[object] = None,
) -> None:
    """Validate high-level density-mixer requests and gate unsupported routes.

    Anderson/Broyden/Kerker are wired in the lower-level EWALD_3D RKS driver
    (`run_rks_periodic_scf` / `run_rks_periodic_multi_k_ewald3d`) and -- since
    the prompt-22 exposure -- in the closed-shell multi-k GDF driver
    (`run_krhf_periodic_gdf`), which `run_periodic_job` reaches via
    `jk_method="gdf"` + `method="RHF"/"RKS"` or via the public multi-k
    RIJCOSX route (`jk_method="rijcosx"`, same GDF loop with
    `k_exchange="cosx"`). All other runner routes (GPW/GAPW, BIPOLE,
    AICCM, open-shell GDF/RIJCOSX) still fail closed: they cannot honour
    these knobs without silently changing the SCF update scheme.
    """
    non_default_parameter = (
        int(density_mixer_depth) != 8
        or float(density_mixer_beta) != 0.5
        or float(kerker_k0) != 1.5
        or float(kerker_strength) != 1.0
        or float(kerker_cutoff_ha) != 120.0
    )
    if (
        density_mixer is None
        and not density_mixer_kerker
        and not non_default_parameter
    ):
        return

    mixer_key = None if density_mixer is None else str(density_mixer).strip().lower()
    if mixer_key in (None, "", "none"):
        if density_mixer_kerker or non_default_parameter:
            raise ValueError(
                "run_periodic_job: density_mixer_* and kerker_* options require "
                "density_mixer='anderson' or 'broyden'."
            )
        return
    if mixer_key == "diis":
        if density_mixer_kerker or non_default_parameter:
            raise ValueError(
                "run_periodic_job: density_mixer='diis' is the default Fock-DIIS "
                "route and does not accept density_mixer_* or kerker_* options."
            )
        return
    if mixer_key not in {"anderson", "broyden"}:
        raise ValueError(
            "run_periodic_job: density_mixer must be None, 'diis', 'anderson', "
            "or 'broyden'. Use density_mixer_kerker=True to add Kerker "
            "preconditioning to Anderson/Broyden."
        )
    method_u = (method or "").upper()
    if (
        jk_method is not None
        and jk_method in (PeriodicJKMethod.GDF, PeriodicJKMethod.RIJCOSX)
        and method_u in ("RHF", "RKS")
        and kpoints is not None
    ):
        # Supported: the closed-shell multi-k GDF driver honours
        # Anderson/Broyden density mixing with optional Kerker
        # preconditioning (ported from the EWALD_3D driver). Gamma-only
        # requests must pass kpoints=(1,1,1) explicitly -- the runner's
        # default-Gamma GDF path calls run_pbc_gdf_rhf, which has no
        # mixer plumbing, and silently dropping the knob is worse than
        # asking for the Nk=1 mesh.
        return
    raise NotImplementedError(
        "run_periodic_job: density_mixer="
        f"{density_mixer!r} is supported on the closed-shell GDF route "
        "(jk_method='gdf' or multi-k 'rijcosx', method='RHF'/'RKS', "
        "explicit kpoints=) and the "
        "lower-level multi-k EWALD_3D RKS drivers "
        "(vibeqc.run_rks_periodic_scf / run_rks_periodic_multi_k_ewald3d). "
        f"Got method={method!r}, jk_method="
        f"{jk_method.value if jk_method is not None else None!r}, "
        f"kpoints={'set' if kpoints is not None else 'None'}. For a "
        "Gamma-only mixed run pass kpoints=(1,1,1) explicitly. Use "
        "DIIS/FMIXING controls for the GPW/GAPW/BIPOLE/AICCM and "
        "open-shell routes."
    )


_COMPACT_MGGA_DENSITY_MIXER_PROFILES = {
    "covalent-insulator",
    "ionic-insulator",
    "metallic-candidate",
}


def _functional_is_scan_family(functional: Optional[str]) -> bool:
    key = (
        str(functional or "")
        .strip()
        .lower()
        .replace("-", "")
        .replace("_", "")
    )
    return key == "scan" or key.startswith("r2scan")


def _runner_bloch_kmesh(
    system: PeriodicSystem,
    kpoints: Optional[
        Union[Tuple[int, int, int], List[int], int, "KPoints", "BlochKMesh"]
    ],
):
    """Materialize the high-level runner's k-point input as a BlochKMesh."""
    from ._vibeqc_core import monkhorst_pack as _mp

    if kpoints is None:
        return _mp(system, [1, 1, 1])
    if hasattr(kpoints, "to_bloch_kmesh") or (
        hasattr(kpoints, "kpoints") and hasattr(kpoints, "weights")
    ):
        from .kpoints import as_bloch_kmesh

        return as_bloch_kmesh(kpoints)
    if isinstance(kpoints, (list, tuple)):
        mesh = list(kpoints)
    else:
        mesh = [kpoints, kpoints, kpoints]
    return _mp(system, [int(n) for n in mesh])


def _bloch_kmesh_size(kmesh) -> int:
    for attr in ("n_kpoints", "nkpts", "num_kpoints"):
        n = getattr(kmesh, attr, None)
        if isinstance(n, int):
            return int(n)
    for attr in ("kpoints_cart", "kpoints", "kpoints_frac"):
        pts = getattr(kmesh, attr, None)
        if pts is not None:
            return int(len(pts))
    return int(len(kmesh))


def _is_multik_kpoints(kpoints) -> bool:
    """Best-effort: does this ``kpoints`` request more than the Γ-point?

    Used to gate the Γ-only READ restart in run_periodic_job. The multi-k
    drivers carry their own defensive READ gate, so a miss here is caught
    downstream; this exists to give the clearer up-front error message.
    """
    if kpoints is None:
        return False
    if isinstance(kpoints, bool):
        return False
    if isinstance(kpoints, int):
        return kpoints > 1
    if isinstance(kpoints, (list, tuple)):
        vals = list(kpoints)
        if len(vals) == 3 and all(isinstance(v, int) for v in vals):
            return int(vals[0]) * int(vals[1]) * int(vals[2]) > 1
        return len(vals) > 1
    for attr in ("n_kpoints", "nkpts", "num_kpoints"):
        n = getattr(kpoints, attr, None)
        if isinstance(n, int):
            return n > 1
    kpts = getattr(kpoints, "kpoints", None)
    if kpts is not None:
        try:
            return len(kpts) > 1
        except TypeError:
            pass
    return True  # unknown explicit object -- assume multi-k (fail closed)


def _reduce_system_to_primitive(
    system: PeriodicSystem,
    *,
    symprec: float = 1e-4,
) -> tuple[PeriodicSystem, SpaceGroup]:
    """Reduce a ``PeriodicSystem`` to its primitive cell via spglib.

    Returns ``(primitive_system, space_group)`` where ``space_group``
    is the symmetry analysis of the **original** cell. The primitive
    system carries its own symmetry analysis as well.

    Raises ``ValueError`` if the primitive cell is identical to the
    input (i.e. the input is already primitive), since no reduction
    is possible.
    """
    # Attach symmetry to the original cell first.
    attach_symmetry(system, symprec=symprec)
    sg_original = system.symmetry
    if sg_original is None:
        raise RuntimeError("attach_symmetry did not populate system.symmetry")

    # Build a Crystal from the PeriodicSystem so we can call to_primitive.
    L = np.asarray(system.lattice, dtype=float, order="F")
    n_atoms_in = len(system.unit_cell)

    # Fractional coordinates: r_frac = L^{-1} . r_cart
    inv_L = np.linalg.inv(L)
    frac_coords = np.empty((3, n_atoms_in), dtype=float, order="F")
    species = []
    for i, atom in enumerate(system.unit_cell):
        r_cart = np.array(atom.xyz, dtype=float)
        frac = inv_L @ r_cart
        frac_coords[:, i] = frac
        species.append(int(atom.Z))

    crystal_in = Crystal(L, frac_coords, species)
    crystal_prim = to_primitive(crystal_in, symprec=symprec)

    n_atoms_out = crystal_prim.n_atoms
    if n_atoms_out == n_atoms_in:
        raise ValueError(
            "reduce_to_primitive: input cell is already primitive "
            f"({n_atoms_in} atoms, space group "
            f"{sg_original.international_symbol} "
            f"(No. {sg_original.number})). "
            "Set reduce_to_primitive=False to skip reduction."
        )

    # Build a new PeriodicSystem from the primitive Crystal.
    prim_lattice = np.asarray(crystal_prim.lattice, dtype=float, order="F")
    prim_atoms = []
    for col in range(n_atoms_out):
        r_cart = prim_lattice @ np.asarray(
            crystal_prim.fractional_coords[:, col], dtype=float
        )
        prim_atoms.append(Atom(int(crystal_prim.species[col]), r_cart.tolist()))

    system_prim = PeriodicSystem(
        dim=system.dim,
        lattice=prim_lattice,
        unit_cell=prim_atoms,
        charge=system.charge,
        multiplicity=system.multiplicity,
    )
    # Attach symmetry to the primitive cell too.
    attach_symmetry(system_prim, symprec=symprec)

    return system_prim, sg_original


def _build_primitive_summary(
    system_in: PeriodicSystem,
    system_prim: PeriodicSystem,
    sg: SpaceGroup,
) -> str:
    """Human-readable summary of the cell reduction."""
    n_in = len(system_in.unit_cell)
    n_out = len(system_prim.unit_cell)
    ratio = n_in / n_out if n_out > 0 else 1.0
    L_in = np.asarray(system_in.lattice, dtype=float)
    L_out = np.asarray(system_prim.lattice, dtype=float)
    vol_in = float(abs(np.linalg.det(L_in)))
    vol_out = float(abs(np.linalg.det(L_out)))

    lines = [
        "  Cell reduction (reduce_to_primitive=True)",
        "  " + "-" * 56,
        f"    space group      = {sg.international_symbol} "
        f"(No. {sg.number}), point group {sg.point_group}",
        f"    symmetry order   = {sg.order}",
        f"    input cell       = {n_in} atoms, volume = {vol_in:.4f} bohr^3",
        f"    primitive cell   = {n_out} atoms, volume = {vol_out:.4f} bohr^3",
        f"    reduction factor = {ratio:.1f}x fewer atoms",
    ]
    if sg.equivalent_atoms:
        lines.append(
            f"    inequivalent     = {len(set(sg.equivalent_atoms))} atom type(s)"
        )
    lines.append("")
    return "\n".join(lines)


# ============================================================
# Helpers -- text output sections
# ============================================================


def _system_summary(system: PeriodicSystem) -> str:
    """Print lattice + atoms in bohr."""
    L = np.asarray(system.lattice, dtype=float)
    dim = int(system.dim)
    if dim == 1:
        measure_label = "periodic length"
        measure_unit = "bohr"
        measure = float(np.linalg.norm(L[:, 0]))
    elif dim == 2:
        measure_label = "periodic area"
        measure_unit = "bohr^2"
        measure = float(np.linalg.norm(np.cross(L[:, 0], L[:, 1])))
    else:
        measure_label = "cell volume"
        measure_unit = "bohr^3"
        measure = float(abs(np.linalg.det(L)))
    periodicity = HeaderlessBlock(
        "Periodicity",
        [HeaderlessColumn("<")],
        body_indent=2,
    )
    # These are label:value rows, but their historical contract is a
    # *minimum* 14-character label field rather than a table-wide aligned
    # column.  A longer runtime label ("periodic length") must not push the
    # equals signs on the other rows to the right.
    periodicity.add_row(f"{'dimensionality':<14s} = {dim}D")
    active_axes = ", ".join(f"a{i + 1}" for i in range(dim))
    periodicity.add_row(f"{'active axes':<14s} = {active_axes}")
    periodicity.add_row(
        f"{measure_label:<14s} = {measure:.4f} {measure_unit}"
    )

    lattice = HeaderlessBlock(
        "Lattice (bohr)",
        [
            HeaderlessColumn("<", min_width=2),
            HeaderlessColumn("<", min_width=1),
            HeaderlessColumn(">", min_width=14),
            HeaderlessColumn(">", min_width=14),
            HeaderlessColumn(">", min_width=14),
        ],
        body_indent=2,
        gutter=1,
    )
    # ``system.lattice`` columns are the Cartesian lattice vectors (C++
    # ``PeriodicSystem.lattice``: "Columns = Cartesian lattice vectors"), so
    # a{i+1} is column i -- ``L[:, i]``. Printing the rows (``L[i, :]``) would
    # display the transpose, which agrees with the true vectors only for a
    # symmetric lattice matrix and is wrong for a skewed / triclinic cell. The
    # measure lines above (periodic length/area) already read columns.
    for i in range(3):
        lattice.add_row(
            f"a{i + 1}",
            "=",
            f"{L[0, i]:14.8f}",
            f"{L[1, i]:14.8f}",
            f"{L[2, i]:14.8f}",
        )
    if dim < 3:
        supercell_volume = float(abs(np.linalg.det(L)))
        lattice.footer(f"embedding volume = {supercell_volume:.4f} bohr^3")

    atoms = HeaderlessBlock(
        f"Atoms (bohr) -- {len(system.unit_cell)} in unit cell",
        [
            HeaderlessColumn(">", min_width=4),
            HeaderlessColumn(">", min_width=5, gutter_after=3),
            HeaderlessColumn(">", min_width=14),
            HeaderlessColumn(">", min_width=14),
            HeaderlessColumn(">", min_width=14),
        ],
        body_indent=2,
        gutter=2,
    )
    for i, atom in enumerate(system.unit_cell, start=1):
        x, y, z = atom.xyz
        atoms.add_row(
            i,
            f"Z={atom.Z:3d}",
            f"{x:14.8f}",
            f"{y:14.8f}",
            f"{z:14.8f}",
        )
    atoms.footer(
        f"n_electrons = {system.n_electrons()}  "
        f"multiplicity = {system.multiplicity}"
    )
    return "\n\n".join(
        (
            periodicity.render(active_policy()),
            lattice.render(active_policy()),
            atoms.render(active_policy()),
        )
    ) + "\n"


def _basis_summary(basis: BasisSet) -> str:
    return (
        "  Basis\n"
        "  " + "-" * 56 + "\n"
        f"    name    = {basis.name}\n"
        f"    nbasis  = {basis.nbasis}\n"
        f"    nshells = {basis.nshells}\n"
        "\n"
    )


def _smearing_summary(result) -> str:
    """Finite-temperature (smearing) block, appended after the shared SCF
    trace + energy-component breakdown.

    Only the smearing-specific quantities live here now: kBT, electronic
    temperature, entropy, Helmholtz free energy, and the Fermi level. The
    iteration table, the ``converged in N iterations`` line, and the
    energy-component breakdown are produced by the shared
    :func:`vibeqc.output.formats.scf_log.format_scf_trace`, identical to the
    molecular path. Returns ``""`` for a zero-temperature (non-smeared)
    result, which is the common case.
    """
    smearing_T = float(getattr(result, "smearing_temperature", 0.0))
    if smearing_T <= 0.0:
        return ""
    lines = [
        "  Finite-temperature (smearing)",
        "  " + "-" * 56,
        f"    {'kBT_smearing (Ha)':>20s} = {smearing_T:20.10f}",
        f"    {'T_elec (K)':>20s} = "
        f"{hartree_to_kelvin_temperature(smearing_T):20.3f}",
        f"    {'entropy S/kB':>20s} = "
        f"{float(getattr(result, 'entropy', 0.0)):20.10f}",
        f"    {'free_energy (Ha)':>20s} = "
        f"{float(getattr(result, 'free_energy', result.energy)):20.10f}",
        f"    {'fermi_level (Ha)':>20s} = "
        f"{float(getattr(result, 'fermi_level', 0.0)):20.10f}",
        "",
    ]
    return "\n".join(lines) + "\n"


def _band_summary(result) -> str:
    """Print band extrema and gap for multi-k periodic jobs.

    Scans all k-point eigenstates and reports the valence-band
    maximum (VBM), conduction-band minimum (CBM), and the
    associated direct/indirect band gaps across the mesh.
    """
    mo_energies = getattr(
        result,
        "mo_energies",
        getattr(result, "mo_energies_alpha", None),
    )
    if mo_energies is None:
        return ""
    if not isinstance(mo_energies, (list, tuple)) or len(mo_energies) <= 1:
        return ""
    occ_list = getattr(result, "occupations", None)
    if not isinstance(occ_list, (list, tuple)) or len(occ_list) == 0:
        # Fallback for result objects without occupation metadata
        # (e.g. PBCBipoleRHFResult). Use energy ordering to infer n_occ.
        eps0 = np.asarray(mo_energies[0])
        n_occ_est = int(np.sum(eps0 < 0))
        if n_occ_est == 0:
            n_e = getattr(result, "n_electrons", 0)
            if hasattr(result, "mo_coeffs") and isinstance(
                result.mo_coeffs, (list, tuple)
            ):
                n_occ_est = result.mo_coeffs[0].shape[1] // 2
            elif n_e:
                n_occ_est = n_e // 2
            else:
                return ""
        # Build synthetic occupation lists: 2.0 for occupied, 0.0 for virtual
        occ_list = []
        for eps_arr in mo_energies:
            n_bf = len(np.asarray(eps_arr))
            occ_k = np.zeros(n_bf, dtype=float)
            occ_k[:n_occ_est] = 2.0
            occ_list.append(occ_k)

    n_k = len(mo_energies)
    # Gather occupied/virtual indices per k-point from occupation thresholds.
    vbm = -1e300  # highest occupied energy across all k
    cbm = +1e300  # lowest unoccupied energy across all k
    vbm_k = -1
    cbm_k = -1
    direct_gap_min = +1e300
    direct_gap_min_k = -1

    for k_idx in range(n_k):
        eps = np.asarray(mo_energies[k_idx])
        occ = np.asarray(occ_list[k_idx], dtype=float)
        occ_mask = occ > 1e-8
        virt_mask = occ < 2.0 - 1e-8
        if not np.any(occ_mask) or not np.any(virt_mask):
            continue
        ho = float(np.max(eps[occ_mask]))
        lu = float(np.min(eps[virt_mask]))
        if ho > vbm:
            vbm = ho
            vbm_k = k_idx
        if lu < cbm:
            cbm = lu
            cbm_k = k_idx
        direct_gap = lu - ho
        if direct_gap < direct_gap_min:
            direct_gap_min = direct_gap
            direct_gap_min_k = k_idx

    if vbm_k < 0 or cbm_k < 0:
        return ""

    indirect_gap = cbm - vbm
    ev = 27.211386245988  # Ha -> eV
    lines = [
        "  Band extrema (multi-k)",
        "  " + "-" * 56,
    ]
    k_ref = getattr(result, "kpoints", None)
    if k_ref is not None and isinstance(k_ref, (list, tuple)):
        _fmt_k = lambda k: ", ".join(f"{x:+.4f}" for x in k)
        if 0 <= vbm_k < len(k_ref):
            lines.append(
                f"    VBM     = {vbm:14.10f} Ha  ({vbm * ev:8.3f} eV)"
                f"  at k = [{_fmt_k(k_ref[vbm_k])}]"
            )
        else:
            lines.append(
                f"    VBM     = {vbm:14.10f} Ha  ({vbm * ev:8.3f} eV)"
                f"  at k_idx = {vbm_k}"
            )
        if 0 <= cbm_k < len(k_ref):
            lines.append(
                f"    CBM     = {cbm:14.10f} Ha  ({cbm * ev:8.3f} eV)"
                f"  at k = [{_fmt_k(k_ref[cbm_k])}]"
            )
        else:
            lines.append(
                f"    CBM     = {cbm:14.10f} Ha  ({cbm * ev:8.3f} eV)"
                f"  at k_idx = {cbm_k}"
            )
    else:
        lines.append(
            f"    VBM     = {vbm:14.10f} Ha  ({vbm * ev:8.3f} eV)  (k_idx = {vbm_k})"
        )
        lines.append(
            f"    CBM     = {cbm:14.10f} Ha  ({cbm * ev:8.3f} eV)  (k_idx = {cbm_k})"
        )
    gap_type = "direct" if vbm_k == cbm_k else "indirect"
    lines.append(
        f"    gap ({gap_type}) = {indirect_gap:10.6f} Ha  ({indirect_gap * ev:8.3f} eV)"
    )
    lines.append(
        f"    direct gap (min)"
        f" = {direct_gap_min:10.6f} Ha  ({direct_gap_min * ev:8.3f} eV)"
        f"  (k_idx = {direct_gap_min_k})"
    )
    lines.append("")
    return "\n".join(lines)


def _mo_summary(result, n_show: int = 20) -> str:
    """Print HOCO/LUCO (crystal orbitals) + nearest energies.

    This function is only called for periodic jobs, so the frontier
    orbitals are crystalline orbitals, not molecular orbitals.
    """
    # UHF / UKS results store per-spin MO energies; fall back to the
    # alpha spin for the summary (beta is symmetric in closed-shell-
    # equivalent cases and informational anyway).
    mo_energies = getattr(
        result,
        "mo_energies",
        getattr(result, "mo_energies_alpha", None),
    )
    if mo_energies is None:
        return ""
    # Multi-k: list of per-k arrays; use the first k-point's energies.
    if isinstance(mo_energies, (list, tuple)):
        if len(mo_energies) == 0:
            return ""
        eps = np.asarray(mo_energies[0])
        occ_list = getattr(result, "occupations", None)
        occ = (
            np.asarray(occ_list[0], dtype=float)
            if isinstance(occ_list, (list, tuple)) and len(occ_list) > 0
            else np.empty(0)
        )
        is_multi_k = True
    else:
        eps = np.asarray(mo_energies)
        occ = np.asarray(getattr(result, "occupations", np.empty(0)), dtype=float)
        is_multi_k = False
    if occ.shape == eps.shape and occ.size:
        occ_mask = occ > 1e-8
        virt_mask = occ < 2.0 - 1e-8
        homo_idx = int(np.where(occ_mask)[0][-1]) if np.any(occ_mask) else -1
        lumo_idx = int(np.where(virt_mask)[0][0]) if np.any(virt_mask) else eps.size
        centre = max(homo_idx + 1, 0) if homo_idx >= 0 else min(lumo_idx, eps.size)
    else:
        # Fallback for older result objects without occupation metadata.
        n_occ = int(np.sum(eps < 0))
        homo_idx = n_occ - 1
        lumo_idx = n_occ
        centre = n_occ
    if is_multi_k:
        n_k = len(mo_energies)
        lines = [
            f"  Crystal orbital energies (Ha) at k-point 0"
            f" (of {n_k}) -- sorted low -> high",
            "  " + "-" * 56,
        ]
    else:
        lines = [
            "  Crystal orbital energies (Ha) -- sorted low -> high",
            "  " + "-" * 56,
        ]
    half = n_show // 2
    lo = max(0, centre - half)
    hi = min(eps.size, centre + half)
    for i in range(lo, hi):
        marker = " HOCO" if i == homo_idx else " LUCO" if i == lumo_idx else ""
        if occ.shape == eps.shape and occ.size:
            lines.append(
                f"    {i + 1:4d}   {eps[i]:18.10f}   occ={occ[i]:8.5f}{marker}"
            )
        else:
            lines.append(f"    {i + 1:4d}   {eps[i]:18.10f}{marker}")
    lines.append("")
    return "\n".join(lines)


def _is_lattice_matrix_set_like(obj) -> bool:
    return (
        hasattr(obj, "cells") and hasattr(obj, "blocks") and hasattr(obj, "set_block")
    )


def _real_block_for_periodic_output(
    block,
    *,
    label: str,
    abs_tol: float = 1.0e-10,
    rel_tol: float = 1.0e-7,
) -> np.ndarray:
    """Return a real matrix block for real scalar output artifacts."""
    arr = np.asarray(block)
    if np.iscomplexobj(arr):
        if arr.size:
            max_imag = float(np.max(np.abs(arr.imag)))
            max_real = float(np.max(np.abs(arr.real)))
        else:
            max_imag = 0.0
            max_real = 0.0
        limit = max(float(abs_tol), float(rel_tol) * max(max_real, 1.0))
        if max_imag > limit:
            raise ValueError(
                f"{label} has a non-negligible imaginary component "
                f"(max|Im|={max_imag:.2e}, max|Re|={max_real:.2e}, "
                f"tolerance={limit:.2e}); periodic density output is "
                "real-only and the k-point density fold is not "
                "time-reversal consistent"
            )
        arr = arr.real
    return np.ascontiguousarray(arr, dtype=float)


def _system_with_valid_unit_cell_multiplicity(
    system: PeriodicSystem,
) -> PeriodicSystem:
    """Return ``system`` or an output-only copy with a valid unit-cell spin."""
    unit = system.unit_cell_molecule()
    charge = int(getattr(unit, "charge", system.charge))
    multiplicity = int(getattr(unit, "multiplicity", system.multiplicity))
    if charge == int(system.charge) and multiplicity == int(system.multiplicity):
        return system
    atoms = [
        Atom(int(atom.Z), [float(x) for x in atom.xyz])
        for atom in system.unit_cell
    ]
    return PeriodicSystem(
        int(system.dim),
        np.asarray(system.lattice, dtype=float),
        atoms,
        charge=charge,
        multiplicity=multiplicity,
    )


def _fold_per_k_density_for_periodic_output(
    basis: BasisSet,
    system: PeriodicSystem,
    D_per_k: Sequence[np.ndarray],
    kpoints_cart,
    weights,
    lat_opts,
):
    """Inverse-Bloch fold per-k AO density matrices for real output grids."""
    from ._vibeqc_core import compute_overlap_lattice

    D_set = compute_overlap_lattice(basis, system, lat_opts)
    kpts = np.asarray(kpoints_cart, dtype=float).reshape(-1, 3)
    w_arr = np.asarray(weights, dtype=float).reshape(-1)
    if len(D_per_k) == 0:
        raise ValueError("SCF result has an empty density list")
    if kpts.shape[0] != len(D_per_k) or w_arr.shape[0] != len(D_per_k):
        raise ValueError(
            "per-k density output requires aligned density, k-point, and "
            f"weight lists; got {len(D_per_k)} densities, {kpts.shape[0]} "
            f"k-points, {w_arr.shape[0]} weights"
        )
    total_w = float(w_arr.sum())
    if total_w <= 0.0:
        raise ValueError("per-k density output requires positive k-point weights")
    if not np.isclose(total_w, 1.0, rtol=1.0e-9, atol=1.0e-12):
        w_arr = w_arr / total_w

    terms = _time_reversal_completed_density_terms(
        D_per_k,
        kpts,
        w_arr,
    )
    for cell_idx, cell in enumerate(D_set.cells):
        r_cart = np.asarray(cell.r_cart, dtype=float).reshape(3)
        P_g = np.zeros((basis.nbasis, basis.nbasis), dtype=np.complex128)
        for weight, k_cart, D in terms:
            phase = np.exp(-1j * float(np.dot(k_cart, r_cart)))
            P_g += float(weight) * phase * D
        D_set.set_block(
            cell_idx,
            _real_block_for_periodic_output(
                P_g,
                label=f"periodic density block g={tuple(cell.index)}",
            ),
        )
    return D_set


def _wrap_reciprocal_fractional(frac: np.ndarray) -> np.ndarray:
    """Wrap reciprocal fractional coordinates into [-0.5, 0.5)."""
    arr = np.asarray(frac, dtype=float)
    return arr - np.floor(arr + 0.5)


def _fractional_kpoints_for_output(
    system: PeriodicSystem,
    kpoints_cart: np.ndarray,
) -> np.ndarray:
    """Return reciprocal fractional k coordinates for output folding."""
    B = np.asarray(system.reciprocal_lattice(), dtype=float)
    return (np.linalg.pinv(B) @ np.asarray(kpoints_cart, dtype=float).T).T


def _time_reversal_completed_density_terms(
    D_per_k: Sequence[np.ndarray],
    kpoints_cart: np.ndarray,
    weights: np.ndarray,
) -> list[tuple[float, np.ndarray, np.ndarray]]:
    """Return density-fold terms symmetrized over time reversal.

    Every k point is folded as the explicit conjugate pair

        w/2 exp(-i k.R) D(k) + w/2 exp(+i k.R) D(k)^*

    so each lattice block of the fold is 2 Re(...) -- real to machine
    precision by construction. This both materializes the implicit -k
    partner of a symmetry-reduced mesh and time-reversal-symmetrizes an
    explicit +/-k mesh: the converged per-k densities satisfy
    D(-k) = D(k)^* only up to the per-k lattice-truncation asymmetry of
    the Fock build (~1e-6 relative; see
    handovers/HANDOVER_GDF_FIT_SCREENING.md "Open findings"), so trusting
    the explicit partner leaves a spurious imaginary residue of that size
    in the folded blocks. The symmetrized fold equals the real part of
    the pristine fold, which is exact for the time-reversal-symmetric
    Hamiltonians vibe-qc ships (no magnetic field / spin-orbit coupling).
    """
    kpts = np.asarray(kpoints_cart, dtype=float).reshape(-1, 3)
    w_arr = np.asarray(weights, dtype=float).reshape(-1)
    terms: list[tuple[float, np.ndarray, np.ndarray]] = []

    for weight, k_cart, D_k in zip(w_arr, kpts, D_per_k):
        D = np.asarray(D_k, dtype=np.complex128)
        D = 0.5 * (D + D.conj().T)
        half = 0.5 * float(weight)
        k = np.asarray(k_cart, dtype=float)
        terms.append((half, k, D))
        terms.append((half, -k, D.conj()))
    return terms


def _result_kpoints_cart(result) -> Optional[np.ndarray]:
    """Return result k-points in Cartesian bohr^-1 coordinates, if present."""
    kpts = getattr(result, "kpoints_cart", None)
    if kpts is None:
        kpts = getattr(result, "kpoints", None)
    if kpts is None:
        kmesh = getattr(result, "kmesh", None)
        if kmesh is not None:
            kpts = getattr(kmesh, "kpoints_cart", None)
            if kpts is None:
                kpts = getattr(kmesh, "kpoints", None)
    if kpts is None:
        return None
    arr = np.asarray(kpts, dtype=float)
    if arr.size == 0:
        return None
    return arr.reshape(-1, 3)


def _result_kpoint_weights(result) -> Optional[np.ndarray]:
    """Return result k-point weights, if present."""
    weights = getattr(result, "kpoint_weights", None)
    if weights is None:
        weights = getattr(result, "weights", None)
    if weights is None:
        kmesh = getattr(result, "kmesh", None)
        if kmesh is not None:
            weights = getattr(kmesh, "weights", None)
    if weights is None:
        return None
    arr = np.asarray(weights, dtype=float)
    if arr.size == 0:
        return None
    return arr.reshape(-1)


def _density_proxy_with_k_metadata(result, density) -> SimpleNamespace:
    """Proxy one spin density while preserving k metadata for output folds."""
    payload: dict[str, object] = {"density": density}
    for name in ("kpoints_cart", "kpoint_weights", "kpoints", "weights", "kmesh"):
        if hasattr(result, name):
            payload[name] = getattr(result, name)
    return SimpleNamespace(**payload)


def _gamma_index_for_multi_k(result, n_items: int) -> int:
    """Find the actual Gamma point for Gamma-only writer fallbacks."""
    if n_items <= 0:
        raise ValueError("multi-k result has no k-point data")
    if n_items == 1:
        return 0
    kpts = _result_kpoints_cart(result)
    if kpts is None:
        raise ValueError(
            "multi-k Gamma-only output requires k-point metadata; "
            "cannot assume the first k-point is Gamma"
        )
    if kpts.shape[0] != n_items:
        raise ValueError(
            "multi-k Gamma-only output requires aligned k-point metadata; "
            f"got {n_items} data blocks and {kpts.shape[0]} k-points"
        )
    norms = np.linalg.norm(kpts, axis=1)
    gamma_idx = int(np.argmin(norms))
    if float(norms[gamma_idx]) > 1.0e-10:
        raise ValueError(
            "multi-k Gamma-only output requires an explicit Gamma k-point; "
            f"closest |k|={float(norms[gamma_idx]):.2e} bohr^-1"
        )
    return gamma_idx


def _density_lattice_set_for_output(
    basis: BasisSet,
    system: PeriodicSystem,
    result,
    lat_opts,
):
    """Return a LatticeMatrixSet density for periodic grid/DOS writers."""
    D_attr = getattr(result, "density", None)
    if D_attr is None:
        raise ValueError("SCF result has no .density attribute")
    if _is_lattice_matrix_set_like(D_attr):
        from ._vibeqc_core import compute_overlap_lattice

        # SCF drivers may converge a density on a shorter lattice cutoff
        # than the fixed output grid/DOS template. Build that target first,
        # then zero-extend the converged density by lattice-cell key. Never
        # use positional alignment: equal-radius cells have no ordering
        # contract across independently constructed lattice sets.
        D_set = compute_overlap_lattice(basis, system, lat_opts)
        target_cells = list(D_set.cells)
        target_blocks = list(D_set.blocks)
        if len(target_cells) != len(target_blocks):
            raise ValueError(
                "density output template received misaligned cells and blocks"
            )
        target_keys = [
            tuple(int(value) for value in cell.index) for cell in target_cells
        ]
        if len(set(target_keys)) != len(target_keys):
            raise ValueError("density output template received duplicate cells")

        source_cells = list(D_attr.cells)
        source_blocks = list(D_attr.blocks)
        if len(source_cells) != len(source_blocks):
            raise ValueError(
                "density lattice-set output received misaligned cells and blocks"
            )

        block_by_cell: dict[tuple[int, int, int], np.ndarray] = {}
        for cell, block in zip(source_cells, source_blocks):
            key = tuple(int(value) for value in cell.index)
            if key in block_by_cell:
                raise ValueError(
                    "density lattice-set output received duplicate cell "
                    f"{key}"
                )
            real_block = _real_block_for_periodic_output(
                block,
                label=f"periodic density block g={key}",
            )
            expected_shape = (int(D_set.nbf), int(D_set.nbf))
            if real_block.shape != expected_shape:
                raise ValueError(
                    "density lattice-set output block "
                    f"g={key} has shape {real_block.shape}, expected "
                    f"{expected_shape}"
                )
            block_by_cell[key] = real_block

        source_only = set(block_by_cell).difference(target_keys)
        if source_only:
            missing = sorted(source_only)
            raise ValueError(
                "density output template omits converged lattice cells: "
                f"{missing}"
            )

        zero = np.zeros((int(D_set.nbf), int(D_set.nbf)), dtype=float)
        for i, key in enumerate(target_keys):
            D_set.set_block(i, block_by_cell.get(key, zero))
        return D_set

    from ._vibeqc_core import compute_overlap_lattice

    if isinstance(D_attr, (list, tuple)):
        kpoints_cart = _result_kpoints_cart(result)
        weights = _result_kpoint_weights(result)
        has_complex_blocks = any(np.iscomplexobj(np.asarray(d)) for d in D_attr)
        if has_complex_blocks and (kpoints_cart is None or weights is None):
            raise ValueError(
                "complex per-k density output requires k-point and weight "
                "metadata so the density can be inverse-Bloch folded; "
                "refusing to real-project individual k blocks"
            )
        if kpoints_cart is not None and weights is not None:
            return _fold_per_k_density_for_periodic_output(
                basis,
                system,
                D_attr,
                kpoints_cart,
                weights,
                lat_opts,
            )
        if len(D_attr) == 0:
            raise ValueError("SCF result has an empty density list")
        D_set = compute_overlap_lattice(basis, system, lat_opts)
        D_home = sum(
            _real_block_for_periodic_output(
                d,
                label=f"periodic density k-point {i}",
            )
            for i, d in enumerate(D_attr)
        ) / len(D_attr)
    else:
        D_set = compute_overlap_lattice(basis, system, lat_opts)
        D_home = _real_block_for_periodic_output(
            D_attr,
            label="periodic density matrix",
        )
    zero = np.zeros_like(D_home)
    for i in range(len(D_set)):
        D_set.set_block(i, D_home if i == 0 else zero)
    return D_set


def _sum_lattice_density_sets_for_output(
    basis: BasisSet,
    system: PeriodicSystem,
    lat_opts,
    *density_sets,
    label: str,
):
    """Add spin-resolved lattice density sets for real output artifacts."""
    if not density_sets:
        raise ValueError("no density lattice sets supplied")
    ref = density_sets[0]
    ref_cells = list(ref.cells)
    ref_blocks = list(ref.blocks)
    if len(ref_cells) != len(ref_blocks):
        raise ValueError("spin density lattice set has misaligned cells and blocks")
    ref_keys = [tuple(int(value) for value in cell.index) for cell in ref_cells]
    for other in density_sets[1:]:
        other_cells = list(other.cells)
        other_blocks = list(other.blocks)
        if len(other_cells) != len(other_blocks):
            raise ValueError(
                "spin density lattice set has misaligned cells and blocks"
            )
        other_keys = [
            tuple(int(value) for value in cell.index) for cell in other_cells
        ]
        if other_keys != ref_keys:
            raise ValueError(
                "spin density lattice sets have different keyed cell order"
            )

    from vibeqc.pbc_bipole_common import _copy_lattice_with_blocks

    blocks = []
    for i, parts in enumerate(zip(*(list(ds.blocks) for ds in density_sets))):
        cell = ref_cells[i]
        combined = sum(np.asarray(part) for part in parts)
        blocks.append(
            _real_block_for_periodic_output(
                combined,
                label=f"{label} block g={tuple(cell.index)}",
            )
        )
    return _copy_lattice_with_blocks(
        basis,
        system,
        lat_opts,
        ref.cells,
        blocks,
    )


def _evaluate_density_matrix_on_lattice_grid(
    density: np.ndarray,
    basis: BasisSet,
    lattice_bohr: np.ndarray,
    shape: tuple[int, int, int],
    *,
    system: PeriodicSystem | None = None,
    ao_image_radius: int = 1,
) -> tuple[np.ndarray, np.ndarray]:
    """Evaluate a one-cell AO density matrix on a lattice-spanning grid.

    ``PeriodicSystem.lattice`` stores lattice vectors as columns. QVF grids
    store per-voxel vectors as rows, so the returned ``voxel_vectors`` are
    ``lattice.T / shape`` and ``shape[i] * voxel_vectors[i]`` exactly spans
    lattice vector ``i``.
    """
    L_bohr = np.asarray(lattice_bohr, dtype=float)
    if L_bohr.shape != (3, 3):
        raise ValueError("periodic QVF lattice grid requires a 3x3 lattice")
    nx, ny, nz = (int(shape[0]), int(shape[1]), int(shape[2]))
    if min(nx, ny, nz) < 1:
        raise ValueError("periodic QVF lattice grid shape must be positive")
    voxel_vectors = L_bohr.T / np.array([nx, ny, nz], dtype=float)[:, None]
    D = _real_block_for_periodic_output(density, label="periodic QVF density")
    if D.ndim != 2:
        raise ValueError(
            "periodic QVF density matrix must be rank 2; "
            f"got shape {D.shape!r}"
        )

    n_points = nx * ny * nz
    rho_flat = np.empty(n_points, dtype=float)
    chunk_size = 200_000
    for start in range(0, n_points, chunk_size):
        stop = min(start + chunk_size, n_points)
        linear = np.arange(start, stop, dtype=np.int64)
        ix = linear // (ny * nz)
        iy = (linear // nz) % ny
        iz = linear % nz
        frac = np.column_stack(
            [
                ix.astype(float) / float(nx),
                iy.astype(float) / float(ny),
                iz.astype(float) / float(nz),
            ]
        )
        points = frac @ L_bohr.T
        chi = _evaluate_ao_for_periodic_qvf(
            basis,
            system,
            points,
            image_radius=ao_image_radius,
        )
        rho_flat[start:stop] = np.einsum("mi,ij,mj->m", chi, D, chi)
    return rho_flat.reshape((nx, ny, nz)), voxel_vectors


def _evaluate_ao_for_periodic_qvf(
    basis: BasisSet,
    system: PeriodicSystem | None,
    points: np.ndarray,
    *,
    image_radius: int,
) -> np.ndarray:
    """Evaluate AO values with image sums on periodic axes for QVF grids."""
    from ._vibeqc_core import evaluate_ao

    if system is None or int(image_radius) <= 0:
        return evaluate_ao(basis, points)

    import itertools

    from .ewald_j import get_shifted_basis

    dim = max(0, min(3, int(getattr(system, "dim", 3))))
    ranges = [
        range(-int(image_radius), int(image_radius) + 1) if axis < dim else (0,)
        for axis in range(3)
    ]
    chi = evaluate_ao(basis, points)
    for shift in itertools.product(*ranges):
        if shift == (0, 0, 0):
            continue
        shifted = get_shifted_basis(basis, system, shift)
        chi += evaluate_ao(shifted, points)
    return chi


def _evaluate_orbital_on_lattice_grid(
    coefficients: np.ndarray,
    basis: BasisSet,
    lattice_bohr: np.ndarray,
    shape: tuple[int, int, int],
    *,
    system: PeriodicSystem | None = None,
    ao_image_radius: int = 1,
) -> tuple[np.ndarray, np.ndarray]:
    """Evaluate one AO coefficient vector on a lattice-spanning grid."""
    L_bohr = np.asarray(lattice_bohr, dtype=float)
    if L_bohr.shape != (3, 3):
        raise ValueError("periodic QVF lattice grid requires a 3x3 lattice")
    nx, ny, nz = (int(shape[0]), int(shape[1]), int(shape[2]))
    if min(nx, ny, nz) < 1:
        raise ValueError("periodic QVF lattice grid shape must be positive")
    voxel_vectors = L_bohr.T / np.array([nx, ny, nz], dtype=float)[:, None]
    coeff = np.asarray(coefficients, dtype=np.complex128).reshape(-1)
    if coeff.shape[0] != int(basis.nbasis):
        raise ValueError(
            "periodic QVF orbital coefficient length "
            f"{coeff.shape[0]} does not match basis size {int(basis.nbasis)}"
        )

    n_points = nx * ny * nz
    psi_flat = np.empty(n_points, dtype=np.complex128)
    chunk_size = 200_000
    for start in range(0, n_points, chunk_size):
        stop = min(start + chunk_size, n_points)
        linear = np.arange(start, stop, dtype=np.int64)
        ix = linear // (ny * nz)
        iy = (linear // nz) % ny
        iz = linear % nz
        frac = np.column_stack(
            [
                ix.astype(float) / float(nx),
                iy.astype(float) / float(ny),
                iz.astype(float) / float(nz),
            ]
        )
        points = frac @ L_bohr.T
        chi = _evaluate_ao_for_periodic_qvf(
            basis,
            system,
            points,
            image_radius=ao_image_radius,
        )
        psi_flat[start:stop] = chi @ coeff
    psi = _real_block_for_periodic_output(
        psi_flat.reshape((nx, ny, nz)),
        label="periodic QVF orbital grid",
    )
    return psi, voxel_vectors


def _aiccm_b_qvf_translations(mesh: tuple[int, int, int]) -> list[tuple[int, int, int]]:
    return [
        (int(i), int(j), int(k))
        for i in range(int(mesh[0]))
        for j in range(int(mesh[1]))
        for k in range(int(mesh[2]))
    ]


def _aiccm_b_qvf_supercell_system(
    system: PeriodicSystem,
    mesh: tuple[int, int, int],
) -> PeriodicSystem:
    """Build the visual BvK supercell used by χ-CCM-B QVF output."""
    lattice = np.asarray(system.lattice, dtype=float)
    mesh_arr = np.asarray(mesh, dtype=float)
    super_lattice = lattice * mesh_arr[None, :]
    atoms: list[Atom] = []
    for translation in _aiccm_b_qvf_translations(mesh):
        shift = lattice @ np.asarray(translation, dtype=float)
        for atom in system.unit_cell:
            position = np.asarray(atom.xyz, dtype=float) + shift
            atoms.append(Atom(int(atom.Z), position.tolist()))
    return PeriodicSystem(
        int(system.dim),
        super_lattice,
        atoms,
        charge=int(system.charge) * int(np.prod(mesh_arr)),
        multiplicity=1,
    )


def _aiccm_b_residue_density_blocks_for_qvf(
    density_like,
    mesh: tuple[int, int, int],
    *,
    label: str,
) -> dict[tuple[int, int, int], np.ndarray]:
    """Fold lattice-set density aliases to one block per cyclic residue."""
    residues = _aiccm_b_qvf_translations(mesh)
    if not _is_lattice_matrix_set_like(density_like):
        if tuple(mesh) != (1, 1, 1):
            raise ValueError(
                f"{label} is not a lattice matrix set for χ-CCM-B mesh {mesh!r}"
            )
        return {
            (0, 0, 0): _real_block_for_periodic_output(
                density_like,
                label=label,
            )
        }

    buckets: dict[tuple[int, int, int], list[np.ndarray]] = {
        residue: [] for residue in residues
    }
    mesh_arr = np.asarray(mesh, dtype=int)
    for cell, block in zip(density_like.cells, density_like.blocks):
        index = np.asarray(cell.index, dtype=int)
        residue = tuple(int(x) for x in np.mod(index, mesh_arr))
        if residue in buckets:
            buckets[residue].append(
                _real_block_for_periodic_output(
                    block,
                    label=f"{label} residue {residue}",
                )
            )
    missing = [residue for residue, blocks in buckets.items() if not blocks]
    if missing:
        raise ValueError(
            f"{label} is missing cyclic density residues {missing!r}"
        )
    return {
        residue: np.mean(np.stack(blocks, axis=0), axis=0)
        for residue, blocks in buckets.items()
    }


def _aiccm_b_density_alias_can_fold_directly(
    density_like,
    mesh: tuple[int, int, int],
) -> bool:
    if _is_lattice_matrix_set_like(density_like):
        return True
    return tuple(mesh) == (1, 1, 1) and np.asarray(density_like).ndim == 2


def _aiccm_b_effective_electrons_for_qvf(result) -> int:
    diagnostics = getattr(result, "aiccm2026dev_b", None)
    value = getattr(result, "effective_n_electrons", None)
    if value is None and diagnostics is not None:
        value = getattr(diagnostics, "effective_electron_count", None)
    if value is None:
        raise TypeError(
            "χ-CCM-B QVF density folding requires an effective electron count"
        )
    return int(value)


def _aiccm_b_spin_occupations_for_qvf(
    system: PeriodicSystem,
    effective_electrons: int,
) -> tuple[int, int]:
    two_s = int(system.multiplicity) - 1
    alpha_twice = int(effective_electrons) + two_s
    beta_twice = int(effective_electrons) - two_s
    if alpha_twice < 0 or beta_twice < 0 or alpha_twice % 2 or beta_twice % 2:
        raise ValueError(
            "χ-CCM-B QVF density folding requires the effective electron "
            "count and multiplicity to define integer alpha/beta occupations"
        )
    return alpha_twice // 2, beta_twice // 2


def _aiccm_b_residue_density_blocks_from_k_for_qvf(
    result,
    system: PeriodicSystem,
    mesh: tuple[int, int, int],
) -> dict[tuple[int, int, int], np.ndarray]:
    """Fold stored B k-density blocks to one block per cyclic residue."""
    effective_electrons = _aiccm_b_effective_electrons_for_qvf(result)
    if hasattr(result, "density_alpha") and hasattr(result, "density_beta"):
        n_alpha, n_beta = _aiccm_b_spin_occupations_for_qvf(
            system,
            effective_electrons,
        )
        alpha = _spin_density_blocks_per_k(result, "alpha", n_alpha)
        beta = _spin_density_blocks_per_k(result, "beta", n_beta)
        if len(alpha) != len(beta):
            raise ValueError("χ-CCM-B QVF alpha/beta density block counts differ")
        density_k = [a + b for a, b in zip(alpha, beta)]
    else:
        if int(system.multiplicity) != 1 or int(effective_electrons) % 2:
            raise ValueError(
                "χ-CCM-B QVF restricted density folding requires a singlet "
                "record with an even effective electron count"
            )
        density_k = _density_blocks_per_k(result, effective_electrons)

    kpoints_frac = np.asarray(getattr(result, "kpoints_frac"), dtype=float).reshape(
        -1,
        3,
    )
    weights = np.asarray(getattr(result, "kpoint_weights"), dtype=float).reshape(-1)
    if len(density_k) != kpoints_frac.shape[0] or weights.shape != (
        kpoints_frac.shape[0],
    ):
        raise ValueError(
            "χ-CCM-B QVF density folding requires matching density, k-point, "
            "and weight counts"
        )
    residues = _aiccm_b_qvf_translations(mesh)
    blocks = inverse_bloch_transform(density_k, kpoints_frac, residues, weights)
    imaginary_residual = float(np.max(np.abs(blocks.imag))) if blocks.size else 0.0
    if imaginary_residual > 1.0e-7:
        raise NotImplementedError(
            "χ-CCM-B QVF density folding would discard a non-real density "
            f"residue ({imaginary_residual:.3e})"
        )
    return {
        residue: np.ascontiguousarray(block, dtype=float)
        for residue, block in zip(residues, blocks.real)
    }


def _aiccm_b_full_density_matrix_for_qvf(
    result,
    system: PeriodicSystem,
    mesh: tuple[int, int, int],
) -> np.ndarray:
    """Expand χ-CCM-B residue density blocks to a full BvK AO matrix."""

    def full_from_blocks(blocks_by_residue: dict[tuple[int, int, int], np.ndarray]):
        residues = _aiccm_b_qvf_translations(mesh)
        nbf = int(next(iter(blocks_by_residue.values())).shape[0])
        full = np.zeros((len(residues) * nbf, len(residues) * nbf), dtype=float)
        for origin_index, origin in enumerate(residues):
            row = slice(origin_index * nbf, (origin_index + 1) * nbf)
            for target_index, target in enumerate(residues):
                delta = tuple(
                    (int(target[axis]) - int(origin[axis])) % int(mesh[axis])
                    for axis in range(3)
                )
                col = slice(target_index * nbf, (target_index + 1) * nbf)
                full[row, col] = blocks_by_residue[delta]
        return full

    if hasattr(result, "density_alpha") and hasattr(result, "density_beta"):
        if _aiccm_b_density_alias_can_fold_directly(
            result.density_alpha,
            mesh,
        ):
            alpha = _aiccm_b_residue_density_blocks_for_qvf(
                result.density_alpha,
                mesh,
                label="periodic QVF alpha density",
            )
            beta = _aiccm_b_residue_density_blocks_for_qvf(
                result.density_beta,
                mesh,
                label="periodic QVF beta density",
            )
            return full_from_blocks(alpha) + full_from_blocks(beta)
        return full_from_blocks(
            _aiccm_b_residue_density_blocks_from_k_for_qvf(
                result,
                system,
                mesh,
            )
        )
    density_like = getattr(result, "density", None)
    if _aiccm_b_density_alias_can_fold_directly(density_like, mesh):
        density = _aiccm_b_residue_density_blocks_for_qvf(
            density_like,
            mesh,
            label="periodic QVF density",
        )
    else:
        density = _aiccm_b_residue_density_blocks_from_k_for_qvf(
            result,
            system,
            mesh,
        )
    return full_from_blocks(density)


def _aiccm_b_qvf_gamma_orbital_grids(
    result,
    super_basis: BasisSet,
    super_system: PeriodicSystem,
    lattice_bohr: np.ndarray,
    shape: tuple[int, int, int],
    mesh: tuple[int, int, int],
) -> list[dict[str, object]]:
    """Return Γ-character HOMO/LUMO grids in the χ-CCM-B BvK supercell."""
    mo_coeffs = getattr(result, "mo_coeffs", None)
    if not isinstance(mo_coeffs, (list, tuple)) or not mo_coeffs:
        return []
    if hasattr(result, "mo_coeffs_alpha") or hasattr(result, "mo_coeffs_beta"):
        return []

    kfrac = getattr(result, "kpoints_frac", None)
    if kfrac is None:
        return []
    kfrac_arr = np.asarray(kfrac, dtype=float).reshape(-1, 3)
    if kfrac_arr.shape[0] != len(mo_coeffs):
        raise ValueError(
            "aiccm2026dev-b QVF orbital grids require aligned kpoints_frac "
            f"({kfrac_arr.shape[0]}) and mo_coeffs ({len(mo_coeffs)})"
        )
    wrapped = _wrap_reciprocal_fractional(kfrac_arr)
    gamma_idx = int(np.argmin(np.linalg.norm(wrapped, axis=1)))
    if float(np.linalg.norm(wrapped[gamma_idx])) > 1.0e-12:
        raise ValueError(
            "aiccm2026dev-b QVF orbital grids require an explicit Gamma "
            "character in the finite mesh"
        )

    C_gamma = np.asarray(mo_coeffs[gamma_idx], dtype=np.complex128)
    if C_gamma.ndim != 2:
        raise ValueError(
            "aiccm2026dev-b QVF Gamma coefficient block must be rank 2; "
            f"got shape {C_gamma.shape!r}"
        )
    n_ao_cell, n_bands = C_gamma.shape
    n_cells = int(np.prod(np.asarray(mesh, dtype=int)))
    if int(super_basis.nbasis) != n_cells * n_ao_cell:
        raise ValueError(
            "aiccm2026dev-b QVF supercell basis size does not match "
            f"mesh*AO count: {int(super_basis.nbasis)} != {n_cells}*{n_ao_cell}"
        )

    diag = getattr(result, "aiccm2026dev_b", None)
    n_electrons = getattr(result, "effective_n_electrons", None)
    if n_electrons is None and diag is not None:
        n_electrons = getattr(diag, "effective_electron_count", None)
    try:
        n_occ = int(n_electrons) // 2
    except Exception:
        n_occ = 0
    indices: list[int] = []
    if 0 < n_occ <= n_bands:
        indices.append(n_occ - 1)
    if n_occ < n_bands:
        indices.append(n_occ)
    if not indices and n_bands:
        indices.append(0)

    energies = getattr(result, "mo_energies", None)
    if isinstance(energies, (list, tuple)) and len(energies) > gamma_idx:
        eps_gamma = np.asarray(energies[gamma_idx], dtype=float).reshape(-1)
    elif energies is not None:
        eps_gamma = np.asarray(energies, dtype=float).reshape(-1)
    else:
        eps_gamma = np.zeros(n_bands, dtype=float)

    translations = _aiccm_b_qvf_translations(mesh)
    norm = 1.0 / np.sqrt(float(n_cells))
    out: list[dict[str, object]] = []
    for band_index in indices:
        column = np.array(C_gamma[:, band_index], dtype=np.complex128, copy=True)
        if column.size:
            pivot = int(np.argmax(np.abs(column)))
            if abs(column[pivot]) > 0.0:
                column *= np.exp(-1j * float(np.angle(column[pivot])))
        coeff = np.zeros(int(super_basis.nbasis), dtype=np.complex128)
        for cell_index, _translation in enumerate(translations):
            start = cell_index * n_ao_cell
            stop = start + n_ao_cell
            coeff[start:stop] = norm * column
        psi, voxel_vectors = _evaluate_orbital_on_lattice_grid(
            coeff,
            super_basis,
            lattice_bohr,
            shape,
            system=super_system,
        )
        if band_index == n_occ - 1:
            role = "HOMO"
        elif band_index == n_occ:
            role = "LUMO"
        else:
            role = f"band {band_index}"
        out.append(
            {
                "label": f"chi-CCM Gamma {role}",
                "data": psi,
                "origin": np.zeros(3, dtype=float),
                "span": voxel_vectors,
                "band_index": int(band_index),
                "energy_eh": (
                    float(eps_gamma[band_index])
                    if band_index < eps_gamma.shape[0]
                    else 0.0
                ),
                "occupation": 2.0 if band_index < n_occ else 0.0,
                "spin": "both",
                "component": "real",
            }
        )
    return out


def _json_safe_qvf_value(value):
    """Convert small diagnostic payloads to JSON-safe Python values."""
    from dataclasses import asdict, is_dataclass

    if is_dataclass(value):
        return _json_safe_qvf_value(asdict(value))
    if isinstance(value, np.ndarray):
        return _json_safe_qvf_value(value.tolist())
    if isinstance(value, np.generic):
        return value.item()
    if isinstance(value, dict):
        return {
            str(key): _json_safe_qvf_value(item)
            for key, item in value.items()
        }
    if isinstance(value, (list, tuple)):
        return [_json_safe_qvf_value(item) for item in value]
    return value


def _aiccm_b_qvf_vendor_sections(result) -> list[dict[str, object]]:
    """Return first-party QVF vendor metadata sections for χ-CCM-B output."""
    diagnostics = getattr(result, "aiccm2026dev_b", None)
    convention = getattr(diagnostics, "finite_torus_convention", None)
    if convention is None:
        return []
    payload = {
        "schema_version": 1,
        "method_selector": "aiccm2026dev-b",
        "prose_name": "chi-CCM",
        "ccm_approach": convention.ccm_approach,
        "ccm_construction": convention.ccm_construction,
        "evaluation_representation": convention.evaluation_representation,
        # Retained verbatim for vendor-schema-v1 readers.  The exact machine
        # identity is the evaluation_representation field above.
        "representation": "finite-character Gamma-centred character mesh",
        "family": "variational finite-BvK-torus CCM",
        "result_backend": getattr(result, "backend", None),
        "backend": getattr(diagnostics, "backend", None),
        "electronic_method": getattr(diagnostics, "electronic_method", None),
        "mesh": _json_safe_qvf_value(getattr(diagnostics, "mesh", None)),
        "n_cyclic_cells": getattr(diagnostics, "n_cyclic_cells", None),
        "finite_torus_convention": _json_safe_qvf_value(convention),
    }
    return [
        {
            "id": "x_vibeqc_aiccm2026dev_b_convention",
            "kind": "x_vibeqc.aiccm2026dev_b_convention",
            "member": "convention",
            "label": "chi-CCM-B finite-torus convention",
            "payload": _json_safe_qvf_value(payload),
        }
    ]


def _aiccm_b_qvf_wannier_center_sections(
    localization_result,
) -> list[dict[str, object]]:
    """Return vibe-view Wannier-centre overlay sections for χ-CCM-B output."""

    entries: list[dict[str, object]] = []
    bohr2_to_ang2 = _BOHR_TO_ANGSTROM * _BOHR_TO_ANGSTROM

    def append_block(block, *, prefix: str = "") -> None:
        centers = np.asarray(getattr(block, "centers_bohr"), dtype=float)
        spreads = np.asarray(getattr(block, "spreads_bohr2"), dtype=float)
        centers = centers.reshape((-1, 3))
        spreads = spreads.reshape((-1,))
        if centers.shape[0] != spreads.shape[0]:
            raise ValueError(
                "χ-CCM-B Wannier centre count does not match spread count"
            )
        for index, (center_bohr, spread_bohr2) in enumerate(zip(centers, spreads)):
            label = (
                f"{prefix}Wannier {index + 1}"
                if prefix
                else f"Wannier {index + 1}"
            )
            entries.append(
                {
                    "center": [
                        float(value) * _BOHR_TO_ANGSTROM
                        for value in center_bohr
                    ],
                    "spread": float(spread_bohr2) * bohr2_to_ang2,
                    "label": f"χ-CCM-B {label}",
                }
            )

    if hasattr(localization_result, "alpha"):
        append_block(localization_result.alpha, prefix="alpha ")
        beta = getattr(localization_result, "beta", None)
        if beta is not None:
            append_block(beta, prefix="beta ")
    else:
        append_block(localization_result)

    if not entries:
        return []
    return [
        {
            "id": "x_ccm_wannier_centers",
            "kind": "x_ccm.wannier_centers",
            "member": "centers",
            "label": "χ-CCM-B Wannier centres",
            "payload": {"centers": _json_safe_qvf_value(entries)},
        }
    ]


def _qvf_extensions_with(
    extensions: dict[str, object] | None,
    namespace: str,
) -> dict[str, object]:
    """Return a QVF extension map containing ``namespace``."""
    merged = dict(extensions or {})
    merged.setdefault(namespace, {"version": "1.0", "critical": False})
    return merged


# ============================================================
# Main entry point
# ============================================================


def _has_valid_mo_coeffs(result) -> bool:
    """Check whether result.mo_coeffs is non-empty (array or list)."""
    mc = getattr(result, "mo_coeffs", None)
    if mc is None:
        return False
    if isinstance(mc, (list, tuple)):
        return len(mc) > 0 and hasattr(mc[0], "size") and mc[0].size > 0
    return hasattr(mc, "size") and mc.size > 0


def _mark_legacy_gamma_gdf_parity_hold(
    result,
    system: PeriodicSystem,
    plog: ProgressLogger,
) -> bool:
    """Tag dense-core Gamma GDF fallback results as held for PySCF parity."""
    if not _gamma_dense_core_gdf_parity_held(system, "rsgdf"):
        return False
    backend = str(getattr(result, "backend", "") or "")
    if "+PARITY_HELD" not in backend:
        result.backend = (
            f"{backend}+PARITY_HELD" if backend else "PARITY_HELD"
        )
    msg = (
        "run_periodic_job: Gamma-only GDF legacy fallback absolute-energy "
        "parity is HELD for dense-core ionic cells. The P01 MgO/STO-3G "
        "audit found matched Ewald nuclear terms but a large electronic "
        "GDF gauge offset versus PySCF; keep this row held until the "
        "electronic gauge is resolved."
    )
    warnings.warn(msg, RuntimeWarning, stacklevel=3)
    plog.info("  WARNING: " + msg)
    return True


def _gamma_proxy_for_multi_k(result) -> _GammaProxy:
    """Wrap a multi-k result to expose Γ-point (k=0) MOs for molden/etc."""
    mo_coeffs = result.mo_coeffs
    if isinstance(mo_coeffs, (list, tuple)):
        gamma_idx = _gamma_index_for_multi_k(result, len(mo_coeffs))
    else:
        gamma_idx = 0
    return _kpoint_proxy_for_multi_k(result, gamma_idx)


def _qvf_wavefunction_proxy_for_multi_k(result, system) -> tuple[object, list[float]]:
    """Return a result proxy and fractional k point for QVF wavefunction.gto.

    QVF can carry one selected complex Bloch wavefunction. Prefer Γ when the
    SCF k list contains it, preserving older archives; otherwise use the first
    k point and record its reciprocal fractional coordinate.
    """
    mo_coeffs = getattr(result, "mo_coeffs", None)
    if not isinstance(mo_coeffs, (list, tuple)):
        return result, [0.0, 0.0, 0.0]
    n_items = len(mo_coeffs)
    if n_items <= 0:
        raise ValueError("multi-k result has no k-point data")
    kpts = _result_kpoints_cart(result)
    if kpts is None:
        if n_items == 1:
            return _kpoint_proxy_for_multi_k(result, 0), [0.0, 0.0, 0.0]
        raise ValueError(
            "multi-k QVF wavefunction output requires k-point metadata"
        )
    if kpts.shape[0] != n_items:
        raise ValueError(
            "multi-k QVF wavefunction output requires aligned k-point metadata; "
            f"got {n_items} data blocks and {kpts.shape[0]} k-points"
        )
    norms = np.linalg.norm(kpts, axis=1)
    gamma_idx = int(np.argmin(norms))
    idx = gamma_idx if float(norms[gamma_idx]) <= 1.0e-10 else 0
    frac = _wrap_reciprocal_fractional(
        _fractional_kpoints_for_output(system, kpts)
    )
    return _kpoint_proxy_for_multi_k(result, idx), [float(x) for x in frac[idx]]


def _kpoint_proxy_for_multi_k(result, idx: int) -> _GammaProxy:
    """Wrap one k-point block of a multi-k result in single-k result shape."""
    mo_coeffs = result.mo_coeffs
    mo_energies = result.mo_energies
    occupations = getattr(result, "occupations", None)
    if isinstance(occupations, (list, tuple)):
        if len(occupations) == 0:
            occupations = None
        elif idx < len(occupations):
            occupations = occupations[idx]
        else:
            raise ValueError(
                "multi-k result has occupation metadata for fewer k-points "
                f"than orbital blocks; requested index {idx}, got "
                f"{len(occupations)} occupation blocks"
            )
    density = result.density
    return _GammaProxy(
        mo_coeffs=mo_coeffs[idx]
        if isinstance(mo_coeffs, (list, tuple))
        else mo_coeffs,
        mo_energies=mo_energies[idx]
        if isinstance(mo_energies, (list, tuple))
        else mo_energies,
        occupations=occupations,
        density=density[idx]
        if isinstance(density, (list, tuple))
        else density,
        overlap=getattr(result, "overlap", None),
    )


class _GammaProxy:
    """Duck-typed result wrapping Γ-point (k=0) of a multi-k result."""

    __slots__ = ("mo_coeffs", "mo_energies", "occupations", "density", "overlap")

    def __init__(self, *, mo_coeffs, mo_energies, occupations, density, overlap):
        self.mo_coeffs = mo_coeffs
        self.mo_energies = mo_energies
        self.occupations = occupations
        self.density = density
        self.overlap = overlap


class _GapwMultiKRunnerProxy:
    """Duck-typed proxy wrapping a ``GpwMultiKScfResult`` for the runner.

    ``GpwMultiKScfResult`` stores per-k data under ``mo_coeffs_k`` /
    ``mo_energies_k`` / ``occupations_k``, but the runner's output code
    (molden, MO summary, density) expects ``mo_coeffs`` / ``mo_energies`` /
    ``occupations`` as per-k list-shaped attributes. This proxy delegates
    those names to the ``*_k`` equivalents and passes all other attribute
    access through to the original result.
    """

    __slots__ = ("_inner",)

    def __init__(self, inner):
        self._inner = inner

    @property
    def mo_coeffs(self):
        return list(self._inner.mo_coeffs_k)

    @property
    def mo_energies(self):
        return list(self._inner.mo_energies_k)

    @property
    def occupations(self):
        return list(self._inner.occupations_k)

    @property
    def e_electronic(self) -> float:
        bd = getattr(self._inner, "breakdown", None)
        if bd is not None:
            return float(
                getattr(bd, "e_kinetic", 0.0)
                + getattr(bd, "e_nuclear_attraction", 0.0)
                + getattr(bd, "e_hartree", 0.0)
                + getattr(bd, "e_xc", 0.0)
                + getattr(bd, "e_hf_exchange", 0.0)
            )
        return 0.0

    @property
    def e_nuclear(self) -> float:
        bd = getattr(self._inner, "breakdown", None)
        if bd is not None:
            return float(getattr(bd, "e_nuclear_repulsion", 0.0))
        return 0.0

    @property
    def e_xc(self) -> float:
        bd = getattr(self._inner, "breakdown", None)
        if bd is not None:
            return float(getattr(bd, "e_xc", 0.0))
        return 0.0

    def __getattr__(self, name):
        # Delegate all other attribute lookups to the inner result.
        # Properties (mo_coeffs etc.) are handled by the class descriptors
        # and never reach __getattr__.
        return getattr(self._inner, name)


def _plan_periodic_semiempirical_method(
    method: str,
    system: PeriodicSystem,
    *,
    kpoints: object | None,
):
    """Return a validated semiempirical route plan, or ``None``."""
    from vibeqc.semiempirical.routes import (
        BOUNDARY_PERIODIC_GAMMA,
        BOUNDARY_PERIODIC_K,
        is_semiempirical_method,
        plan_periodic_semiempirical_route,
    )

    if not is_semiempirical_method(method):
        return None
    boundary = (
        BOUNDARY_PERIODIC_K
        if kpoints is not None
        else BOUNDARY_PERIODIC_GAMMA
    )
    return plan_periodic_semiempirical_route(
        method,
        system,
        boundary=boundary,
    )


def _run_periodic_semiempirical_engine(
    system: PeriodicSystem,
    route_plan,
    *,
    max_iter: int,
    conv_tol: float,
):
    """Run one basis-free Gamma-point periodic semiempirical backend."""
    method_key = route_plan.method_key
    if method_key in ("dftb0", "scc_dftb"):
        from vibeqc._vibeqc_core import semiempirical as _se
        from vibeqc.semiempirical.parameters import default_parameters

        params = default_parameters()
        open_shell = route_plan.spin == "unrestricted"
        if method_key == "dftb0":
            if open_shell:
                return _se.run_udftb0_gamma(system, params)
            return _se.run_dftb0_gamma(system, params)

        opts = _se.PeriodicSCCOptions()
        opts.max_iter = int(max_iter)
        opts.conv_tol_charge = float(conv_tol)
        if open_shell:
            return _se.run_uscc_dftb_gamma(system, params, opts)
        return _se.run_scc_dftb_gamma(system, params, opts)

    if method_key == "gfn2_xtb":
        from vibeqc._vibeqc_core.semiempirical.xtb import run_gfn2_xtb_gamma
        from vibeqc.semiempirical.methods.gfn2_params import load_gfn2_params

        return run_gfn2_xtb_gamma(system, load_gfn2_params())

    if method_key == "pm6":
        from vibeqc.semiempirical.methods.periodic_pm6 import run_pm6_gamma

        return run_pm6_gamma(
            system,
            max_iter=int(max_iter),
            conv_tol=float(conv_tol),
        )

    if method_key in ("om1", "om2", "om3"):
        from vibeqc.semiempirical.methods.periodic_omx import run_omx_gamma

        return run_omx_gamma(
            system,
            variant=method_key,
            max_iter=int(max_iter),
            conv_tol=float(conv_tol),
        )

    raise ValueError(f"unknown periodic semiempirical method {method_key!r}")


def _periodic_semiempirical_attr(result: object, *names: str) -> object | None:
    for name in names:
        value = getattr(result, name, None)
        if value is not None:
            return value
    return None


def _run_periodic_semiempirical_job(
    system: PeriodicSystem,
    *,
    route_plan,
    output: Union[str, os.PathLike],
    dry_run: bool,
    record_hostname: bool,
    citations: bool,
    write_xyz_file: bool,
    write_poscar_file: bool,
    write_xsf_structure_file: bool,
    write_cif_file: bool,
    write_molden_file: bool,
    write_density: bool,
    write_population_file: bool,
    output_qvf: bool,
    jk_method: Union[str, "PeriodicJKMethod"],
    aux_basis: Optional[str],
    smearing: Optional[SmearingOptions],
    smearing_temperature: Union[float, str, None],
    smearing_metallic: Optional[bool],
    smearing_band_gap_hartree: Optional[float],
    dispersion: Optional[Union[str, bool, "D3BJParams"]],
    optimize: bool,
    optimize_cell_requested: bool,
    hessian: bool,
    tddft: bool,
    coop_cohp: bool,
    band_structure: object | None,
    qvf_wannier_centers: bool,
    kpoints: object | None,
    checkpoint_qvf: Union[str, os.PathLike, None],
    dft_plus_u: object | None,
    atomic_spins: object | None,
    spinlock: str | None,
    read_from: object | None,
    restart_from: Union[str, os.PathLike, None],
    functional: str | None,
    max_iter: int,
    conv_tol_energy: float,
):
    """Public ``run_periodic_job`` branch for basis-free SE methods."""
    method_key = route_plan.method_key
    if functional is not None:
        raise ValueError(
            f"run_periodic_job: functional={functional!r} is not used with "
            f"basis-free semiempirical method={method_key!r}."
        )
    unsupported: list[str] = []
    if optimize:
        unsupported.append("optimize")
    if optimize_cell_requested:
        unsupported.append("optimize_cell")
    if hessian:
        unsupported.append("hessian")
    if tddft:
        unsupported.append("tddft")
    if coop_cohp:
        unsupported.append("coop_cohp")
    if band_structure is not None:
        unsupported.append("band_structure")
    if qvf_wannier_centers:
        unsupported.append("qvf_wannier_centers")
    if write_density:
        unsupported.append("write_density")
    jk_label = (
        jk_method.value if isinstance(jk_method, PeriodicJKMethod) else str(jk_method)
    ).strip().lower()
    if jk_label not in ("auto", ""):
        unsupported.append("jk_method")
    if aux_basis is not None:
        unsupported.append("aux_basis")
    smearing_requested = (
        smearing is not None
        or smearing_temperature is not _SMEARING_UNSET
        or smearing_metallic is not None
        or smearing_band_gap_hartree is not None
    )
    if smearing_requested:
        unsupported.append("smearing")
    if dispersion not in (None, False):
        unsupported.append("dispersion")
    if kpoints is not None:
        unsupported.append("kpoints")
    if checkpoint_qvf is not None:
        unsupported.append("checkpoint_qvf")
    if dft_plus_u:
        unsupported.append("dft_plus_u")
    if atomic_spins is not None:
        unsupported.append("atomic_spins")
    if spinlock is not None:
        unsupported.append("spinlock")
    if read_from is not None:
        unsupported.append("read_from")
    if restart_from is not None:
        unsupported.append("restart_from")
    if unsupported:
        opts = ", ".join(sorted(unsupported))
        raise NotImplementedError(
            "run_periodic_job: basis-free periodic semiempirical methods are "
            f"Gamma-point single-point routes; unsupported option(s): {opts}."
        )

    output_stem = Path(os.fspath(output))
    output_stem.parent.mkdir(parents=True, exist_ok=True)
    out_path = output_stem.with_suffix(".out")
    basis_label = "<basis-free>"
    plan = OutputPlan.from_run_job_kwargs(
        output=output_stem,
        method=method_key,
        basis=basis_label,
        functional=None,
        write_molden_file=False,
        write_xyz=write_xyz_file,
        write_poscar=write_poscar_file,
        write_xsf_structure=write_xsf_structure_file,
        write_cif=write_cif_file,
        write_density_xsf=False,
        write_population=False,
        citations=citations,
        crash_dump=False,
        output_qvf=False,
        job_kind="periodic_scf",
    )
    if dry_run or is_dry_run_requested():
        dry_run_manifest(plan, record_hostname=record_hostname)
        return None

    t_start = time.perf_counter()
    _output_writer = OutputWriter(plan, record_hostname=record_hostname)
    _PERIODIC_OUTPUT_WRITER.set(_output_writer)
    try:
        result = _run_periodic_semiempirical_engine(
            system,
            route_plan,
            max_iter=max_iter,
            conv_tol=conv_tol_energy,
        )
        t_total = time.perf_counter() - t_start

        cite_block_text = ""
        cite_manifest_rows: list[dict[str, Any]] = []
        if citations:
            try:
                refs = load_default_database().assemble(
                    method=method_key,
                    basis="",
                    functional=None,
                    periodic=True,
                    uses_integrals=False,
                )
                cite_manifest_rows = citation_manifest_rows(refs)
                _output_writer.dispatch_role(
                    "citations",
                    citations=refs,
                    raise_on_error=True,
                )
                cite_block_text = format_references_block(refs)
            except Exception as exc:
                warn(
                    f"citation emission failed: {type(exc).__name__}: {exc}",
                )

        with OutputChannel.to_file(out_path):
            write(banner() + "\n\n")
            write(
                f"  Job: PERIODIC {method_key.upper()}  "
                "basis=<basis-free>\n"
            )
            write("  Route: Gamma-point periodic semiempirical\n\n")
            write(_system_summary(system))
            write(section_header("Semiempirical options", width=56))
            write(f"    method              = {method_key}\n")
            write(f"    max_iter            = {int(max_iter)}\n")
            write(f"    conv_tol            = {float(conv_tol_energy)}\n\n")
            skipped = []
            if write_molden_file:
                skipped.append("Molden orbitals")
            if write_population_file:
                skipped.append("population analysis")
            if output_qvf:
                skipped.append("QVF")
            if skipped:
                write(
                    "  Basis-free route: skipped Gaussian-basis artefacts "
                    f"({', '.join(skipped)}).\n\n"
                )

            write(section_header("Results", width=56))
            energy = float(getattr(result, "energy"))
            write(
                "    total energy        = "
                f"{render_energy_labeled(energy, width=0, precision=10, sign=True)}\n"
            )
            converged = _periodic_semiempirical_attr(result, "converged")
            if converged is not None:
                write(f"    converged           = {bool(converged)}\n")
            n_iter = _periodic_semiempirical_attr(result, "n_iter", "iterations")
            if n_iter is not None:
                write(f"    iterations          = {int(n_iter)}\n")
            n_basis = _periodic_semiempirical_attr(result, "n_basis")
            if n_basis is not None:
                write(f"    semiempirical basis = {int(n_basis)} functions\n")
            write("\n")

            if citations and cite_block_text:
                write_references_block(block=cite_block_text)
                flush()

        _output_writer.record(out_path, wall_time_s=t_total)

        if write_xyz_file:
            _output_writer.dispatch_role(
                "geometry",
                only_format="extended-xyz",
                system=system,
                energy_ha=float(getattr(result, "energy")),
                comment=f"vibe-qc periodic {method_key}",
                raise_on_error=True,
            )
        if write_poscar_file:
            _output_writer.dispatch_role(
                "geometry",
                only_format="poscar",
                system=system,
                comment=f"vibe-qc periodic {method_key}",
                raise_on_error=True,
            )
        if write_xsf_structure_file:
            _output_writer.dispatch_role(
                "geometry",
                only_format="xsf",
                system=system,
                raise_on_error=True,
            )
        if write_cif_file:
            _output_writer.dispatch_role(
                "geometry",
                only_format="cif",
                system=system,
                comment=f"vibe-qc periodic {method_key}",
                raise_on_error=True,
            )

        if citations and cite_block_text:
            _output_writer.set_citations(cite_manifest_rows)
        _output_writer.finish(wall_seconds=t_total)
        return result
    except Exception:
        _output_writer.crash(wall_seconds=time.perf_counter() - t_start)
        raise


@_periodic_output_lifecycle
def run_periodic_job(
    system: PeriodicSystem,
    basis: Optional[BasisSet],
    *,
    method: str = "RHF",
    functional: Optional[str] = None,
    jk_method: Union[str, "PeriodicJKMethod"] = "auto",
    # Plane-wave grid cutoff (Hartree) for GPW / GAPW routes.
    # Default 300 Ha ≈ 600 Ry gives ~10⁻⁴ Ha grid convergence on
    # compact contracted Gaussians. Ignored for GDF / BIPOLE routes.
    cutoff_ha: float = 300.0,
    aux_basis: Optional[str] = None,
    # GDF backend selection for jk_method="gdf". ``None`` selects rsgdf
    # (the per-path default for the multi-k + Γ open-shell drivers AND, since
    # 2026-06-15, the Γ closed-shell RHF path: a plain dim=3 neutral Γ RHF
    # now routes through run_pbc_gdf_rhf with exxdiv='ewald' -- PySCF µHa
    # parity -- instead of the legacy molecular-limit gamma driver, which is
    # kept only for RKS / dim<3 / charged / smeared / symmetry runs). Set to
    # "compcell" / "rsgdf" / "mdf" to choose explicitly; "mdf" (Mixed
    # Density Fitting) closes the all-electron Gaussian-DF floor. All three
    # explicit values route Γ closed-shell RHF through run_pbc_gdf_rhf.
    # rsgdf_ke_cutoff controls the reciprocal-space auxiliary mesh used by the
    # range-separated GDF backend.  mdf_ke_cutoff is the modest residual-mesh
    # cutoff used only by gdf_method="mdf".
    gdf_method: Optional[str] = None,
    rsgdf_ke_cutoff: float = 200.0,
    rsgdf_tail_ke_cutoff: Optional[float] = None,
    mdf_ke_cutoff: float = 40.0,
    # Electron-repulsion representation for jk_method="aiccm2026dev-b":
    # corrected-gauge direct four-centre, pair-resolved RI, or RIJCOSX.
    aiccm_backend: str = "four_center",
    # Primary B-stream finite-size control: number of primitive lattice
    # vectors in each cyclic supercell direction.  The equivalent k net is
    # derived, never chosen independently.  ``aiccm_wigner_seitz_shells=s``
    # is the odd-cluster shorthand N=2s+1 (s complete layers per side).
    aiccm_lattice_extension: Optional[Union[int, Sequence[int]]] = None,
    aiccm_wigner_seitz_shells: Optional[Union[int, Sequence[int]]] = None,
    # B-stream-only space-group mode. ``diagnostic`` constructs and checks
    # the finite-torus subgroup and k orbits without changing the SCF;
    # ``integrals`` fails closed until petite-list parity is established.
    aiccm_symmetry: str = "off",
    aiccm_symmetry_require_full_group: bool = False,
    output: Union[str, os.PathLike] = "output",
    use_diis: bool = True,
    # Diagonalisation solver for the SCF procedure.
    # "dense"  -- NumPy/ScaLAPACK dense eigh (default).
    # "davidson" -- block-Davidson iterative solver.
    # "lobpcg" -- LOBPCG iterative solver (handled by Python SCF loop).
    solver: str = "dense",
    # Convergence strategy: "auto" fills any convergence knob the user
    # did not set, from a cheap pre-SCF classification of the system
    # (ionic / covalent / metallic / molecular-limit); "off" keeps the
    # plain defaults. Omitted (None) -> auto unless any explicit knob
    # below is given. The chosen strategy, its per-knob values, and the
    # classification reasons are stated in the .out file. v1 scope:
    # jk_method="bipole" and jk_method="gdf" (other routes fall back
    # to plain defaults).
    convergence: Optional[str] = None,
    # None = not given (auto strategy may choose); a float (incl. 0.0)
    # is an explicit user choice that auto never overrides.
    damping: Optional[float] = None,
    fmixing_percent: Optional[float] = None,
    fock_mixing: Optional[float] = None,
    density_mixer: Optional[str] = None,
    density_mixer_depth: int = 8,
    density_mixer_beta: float = 0.5,
    density_mixer_kerker: bool = False,
    kerker_k0: float = 1.5,
    kerker_strength: float = 1.0,
    kerker_cutoff_ha: float = 120.0,
    smearing: Optional[SmearingOptions] = None,
    smearing_temperature: Union[float, str, None] = _SMEARING_UNSET,  # type: ignore[assignment]
    smearing_unit: str = "hartree",
    smearing_method: str = "fermi-dirac",
    smearing_metallic: Optional[bool] = None,
    smearing_band_gap_hartree: Optional[float] = None,
    # BZ integration backend. None / "smearing" (temperature-broadening,
    # the default), or "gilat" (parameter-free Gilat-Raubenheimer net,
    # T=0 tetrahedron-family integrator). When not given (sentinel),
    # auto-reads from a KPoints.recommend() result's .bz_integration
    # attribute. Explicit user args win over KPoints metadata.
    bz_integration: Optional[str] = _BZ_INTEGRATION_UNSET,
    diis_start_iter: int = 2,
    diis_subspace_size: int = 8,
    max_iter: int = 80,
    conv_tol_energy: float = 1e-7,
    initial_guess: str = "SAD",
    write_molden_file: bool = True,
    write_density: bool = False,
    density_spacing_bohr: float = 0.2,
    write_xyz_file: bool = True,
    write_poscar_file: bool = False,
    write_xsf_structure_file: bool = True,
    write_cif_file: bool = True,
    write_population_file: bool = True,
    citations: bool = True,
    dry_run: bool = False,
    memory_override: bool = False,
    record_hostname: bool = True,
    progress: Union[bool, ProgressLogger, None] = None,
    verbose: Optional[int] = None,
    # --- Dispersion ---------------------------------------------------
    dispersion: Optional[Union[str, bool, "D3BJParams"]] = None,
    dispersion_backend: str = "auto",
    dispersion_cutoff_bohr: float = 50.0,
    # --- BIPOLE-specific ---------------------------------------------
    # Direct-lattice cutoff (bohr) for BIPOLE Fock and nuclear sums.
    # ``cutoff_ha`` is a plane-wave grid cutoff and remains GPW/GAPW-only.
    bipole_cutoff_bohr: Optional[float] = None,
    bipole_nuclear_cutoff_bohr: Optional[float] = None,
    # M5 production default for the erfc SR internal ket-image ball.
    # ``None`` restores the historical unpadded traversal.
    sr_image_precision: Optional[float] = 1e-6,
    # M4b (pair-resolved truncation): QQR-style separation-aware
    # screening for the BIPOLE SR erfc J/K build
    # (LatticeSumOptions.sr_range_screening). The M5 precision path
    # enables it automatically; this explicit flag also permits screened
    # historical/unpadded diagnostics. Ignored by non-BIPOLE routes.
    sr_range_screening: bool = False,
    ewald_omega: Optional[float] = None,
    ewald_precision: float = 1e-8,
    use_oda: bool = False,
    oda_trust_lambda_max: float = 1.0,
    use_mom: bool = False,
    # None = not given (auto strategy may choose); a float (incl. 0.0)
    # is an explicit user choice that auto never overrides.
    level_shift: Optional[float] = None,
    # --- Multi-k (GDF, RIJCOSX, and BIPOLE) -------------------------
    kpoints: Optional[
        Union[Tuple[int, int, int], List[int], int, "KPoints", "BlochKMesh"]
    ] = None,
    # --- Cell reduction -----------------------------------------------
    reduce_to_primitive: bool = False,
    symmetry_precision: float = 1e-4,
    symmetry: Union[bool, str] = False,
    symmetry_stabilize: bool = False,
    symmetry_reduce_fock: bool = False,
    # --- BIPOLE multipole far-field ----------------------------------
    use_multipole_far_field: bool = False,
    multipole_l_max: int = 2,
    # --- BIPOLE exchange convention (option (b), 2026-06-10) ----------
    # None = auto: the Ewald exchange split (corrected gauge -- full
    # Bloch density, split K with the exxdiv G=0 correction, no
    # spheropole term) is ON for 3D Γ-only RHF BIPOLE runs and OFF
    # otherwise (multi-k pends the q!=0 LR-exchange channels). Pass
    # False to force the legacy Γ-locality gauge (needed e.g. to
    # combine the multipole far-field branch with Γ sampling, or for
    # the analytic-gradient preview). 'exchange_exxdiv' picks the K=0
    # convention: 'ewald' (probe-charge Madelung; PySCF-equivalent
    # default) or 'none'.
    use_exchange_ewald_split: Optional[bool] = None,
    exchange_exxdiv: str = "ewald",
    # --- Geometry optimization ----------------------------------------
    optimize: bool = False,
    optimize_max_iter: int = 30,
    optimize_conv_tol_grad: float = 1e-4,
    optimize_cell: bool = False,
    # QVF visualisation archive (v1).
    output_qvf: bool = True,
    # Optional vibe-view overlay: localize χ-CCM-B occupied orbitals and
    # embed x_ccm.wannier_centers in the QVF archive.
    qvf_wannier_centers: bool = False,
    # Opt-in live QVF checkpointing for vibe-view hot-reload. When
    # ``checkpoint_qvf`` is set, a running snapshot is atomically
    # (re)written there as the SCF climbs -- every ``checkpoint_every``
    # cycles -- carrying ``provenance.run_status="running"`` and a
    # monotonic ``provenance.checkpoint.seq``; the terminal snapshot is
    # labeled ``"converged"`` / ``"failed"``. ``checkpoint_every=0``
    # keeps the start + end frames but no per-iteration cadence.
    checkpoint_qvf: Optional[Union[str, os.PathLike]] = None,
    checkpoint_every: int = 0,
    # Harmonic vibrational frequencies via finite-difference Hessian.
    hessian: bool = False,
    # Partial Hessian: atom indices to hold fixed (e.g. bulk-like bottom
    # slab layers) so only the unfrozen atoms are displaced -- 6M instead
    # of 6N SCFs. Frozen atoms anchor the cell, so the frequencies are
    # vibrational-only (no gas-phase trans/rot). See HessianFDOptions.
    hessian_frozen_indices: Optional[List[int]] = None,
    # Pre-computed band structure to include in QVF (vibe-view bands plot).
    band_structure: Optional["BandStructure"] = None,
    # Compute COOP/COHP bonding analysis (requires output_qvf=True).
    coop_cohp: bool = False,
    # TD-DFT excited states (Gamma-point only, TDA).
    tddft: bool = False,
    tddft_n_states: int = 5,
    # DOS/PDOS/COOP k-mesh dimensions for QVF output (default [8,8,8]).
    # For 1D/2D systems, consider e.g. [32,1,1] or [12,12,1].
    dos_kmesh: Optional[Sequence[int]] = None,
    # --- DFT+U (Dudarev rotationally-invariant) ------------------------
    # Shipped on main (see docs/user_guide/dft_plus_u.md):
    #   - RHF +U: Γ-only, via vibeqc.run_rhf_periodic_gamma
    #     (multi-k RHF +U raises NotImplementedError).
    #   - RKS +U: Γ and multi-k both supported, via vibeqc.run_rks_periodic.
    #     Multi-k (Increment 4c) uses the k-averaged AO occupation matrix
    #     and adds S(k) V_AO S(k) per k in cpp/src/periodic_scf.cpp.
    #   - UHF / UKS +U: Γ GPW and Γ/multi-k BIPOLE are supported via
    #     their open-shell drivers.
    #   - GDF +U: queued (GDF Γ-only driver has no +U hook yet).
    dft_plus_u: Optional[List["HubbardSite"]] = None,
    # --- ATOMSPIN: broken-symmetry magnetic seed (UHF/UKS only) -------
    # Per-atom spin tag in unit-cell atom order: +1 (majority alpha),
    # -1 (majority beta), 0 (unpolarised). Seeds an AFM / ferrimagnetic
    # g=0 spin pattern (Bloch-sums to a broken-symmetry D(k)). Supported
    # on the Γ UHF/UKS Ewald / GDF / BIPOLE drivers and the multi-k UHF/UKS
    # Ewald drivers.
    # Requires the SAD guess (the default). See docs/roadmap.md Sec.G2.
    atomic_spins: Optional[List[int]] = None,
    # --- READ: restart from a prior SCF -------------------------------
    # Gamma periodic restarts use the prior g=0 cell density, projected onto
    # this cell's basis. Closed-shell multi-k GDF/GPW/GAPW restarts accept an
    # in-memory prior result and rebuild the per-k Bloch density blocks from
    # its coefficients/occupations. File-backed multi-k READ stays fail-closed
    # until QVF carries all-k wavefunction data. See docs/roadmap.md Sec.G2.
    read_from: Optional[object] = None,
    # --- SPINLOCK: broken-symmetry magnetic convergence (UHF/UKS only) -
    # "pattern_hold" holds the seeded (ATOMSPIN) occupied set by maximum
    # overlap (MOM) for ``spinlock_iterations`` cycles, then releases --
    # protecting an AFM seed from collapsing to the symmetric solution.
    # "spin_schedule" runs a two-phase SCF: converge at locked
    # n_alpha-n_beta = ``spinlock_value`` for ``spinlock_iterations`` cycles,
    # then release to the multiplicity target. PATTERN_HOLD is supported on
    # the Γ-Ewald, GDF, BIPOLE and multi-k UKS Ewald drivers; SPIN_SCHEDULE
    # on the Γ-Ewald, GDF and BIPOLE drivers (unsupported (driver, mode)
    # pairs fail closed). See docs/roadmap.md Sec.G2.
    spinlock: Optional[str] = None,
    spinlock_value: int = 0,
    spinlock_iterations: int = 0,
    # Restart from a previous GPW/GAPW calculation (.npz file written
    # by ``save_gpw_result``). Warm-starts the SCF from the saved density.
    restart_from: Optional[Union[str, os.PathLike]] = None,
):
    """Run a periodic SCF job and write the standard output files.

    Mirrors :func:`vibeqc.run_job` but for periodic systems.
    Gamma RHF/RKS GDF is the default when no k-mesh is specified.
    Multi-k RHF/RKS GDF, multi-k RHF/RKS/UHF/UKS RIJCOSX, and all four
    BIPOLE methods are available via ``kpoints``. These routes accept a
    mesh tuple/list, a scalar mesh size, :class:`KPoints`, or a native
    :class:`BlochKMesh`.

    Parameters
    ----------
    system, basis
        Periodic system + AO basis.  For basis-free periodic semiempirical
        methods (`dftb0`, `scc_dftb`, `gfn2_xtb`, `pm6`, `om1`, `om2`,
        `om3`), pass `basis=None`; the runner dispatches before any
        Gaussian-basis setup.
    method
        ``"RHF"``, ``"UHF"``, ``"RKS"`` or ``"UKS"``. Closed-shell
        RHF / RKS default to the Gamma or multi-k GDF path depending
        on ``kpoints``. With ``jk_method="bipole"``, all four methods
        dispatch through the BIPOLE Gamma or multi-k route. With
        ``jk_method="rijcosx"``, Gamma RHF uses the dedicated RIJCOSX
        driver and true multi-k meshes use the GDF/COSX backend for
        RHF/RKS/UHF/UKS.
    functional
        XC functional for ``method="RKS"`` or ``method="UKS"``.
    output
        Path stem; produces ``{output}.out``, ``{output}.system``,
        ``{output}.molden``, ``{output}.xsf`` (when ``write_density``).
    band_structure
        Optional :class:`BandStructure` pre-computed by
        :func:`vibeqc.band_structure` (or ``_hcore``).
        When given together with ``output_qvf=True``, the band
        structure is embedded in the QVF archive so vibe-view can
        render an interactive Plotly band-structure plot.
        Compute it before calling this function -- the same workflow
        used for matplotlib plotting with
        :func:`vibeqc.plot.band_structure_figure`.
    coop_cohp
        When ``True`` together with ``output_qvf=True``, compute
        COOP and COHP bonding analysis and embed ``dos.coop`` +
        ``dos.cohp`` sections in the QVF archive. Uses the same
        DOS k-mesh as the total/projected DOS (``[8,8,8]``).
        The Hcore matrix (T + V) needed for COHP is computed
        independently within the DOS/COHP k-mesh block; no
        additional user input is required.
    tddft
        When True, compute TD-DFT vertical excitation energies via
        the Tamm-Dancoff approximation (TDA) at the Gamma point.
        Writes excitation energies, oscillator strengths, and
        dominant transitions to the .out file. Requires
        ``_has_valid_mo_coeffs(result)`` (true for Γ-only and
        multi-k results). Not embedded in QVF yet.
    tddft_n_states
        Number of excited states to compute when ``tddft=True``.
        Default 5.
    dos_kmesh
        Override the DOS/PDOS/COOP k-mesh dimensions. Default
        ``[8, 8, 8]``. For 1D systems use e.g. ``[32, 1, 1]``;
        for 2D use ``[12, 12, 1]``. Only used when
        ``output_qvf=True``.
    qvf_wannier_centers
        When True for ``jk_method="aiccm2026dev-b"``, localize the occupied
        finite-torus space with the B-owned Wannier gauge and embed an
        ``x_ccm.wannier_centers`` vendor section in the QVF archive for
        vibe-view's centre overlay. Requires ``output_qvf=True``.
    aux_basis
        Optional auxiliary basis for ``jk_method="gdf"`` and the RI-J part
        of ``jk_method="rijcosx"``. If omitted, vibe-qc chooses the current
        native-GDF default for ``basis.name``.
    bipole_cutoff_bohr, bipole_nuclear_cutoff_bohr
        Direct-lattice cutoff radii in bohr for ``jk_method="bipole"``.
        ``bipole_cutoff_bohr`` controls the electronic BIPOLE Fock sums; when
        ``bipole_nuclear_cutoff_bohr`` is omitted, the BIPOLE route keeps the
        nuclear/Ewald real-space cutoff no longer than the electronic J/K
        cutoff so neutral-cell cancellation stays in the corrected Ewald
        gauge. These are separate from
        ``cutoff_ha``, which is a GPW/GAPW plane-wave grid cutoff.
    convergence
        Convergence-strategy selector. ``"auto"`` classifies the system
        from cheap pre-SCF signals (composition electronegativity
        spread, cell volume, vacuum axes, electron parity) into a
        profile -- ionic-insulator, covalent-insulator,
        metallic-candidate or molecular-limit -- and fills every
        convergence knob the user did not set (Fermi-Dirac smearing,
        FMIXING, level shift, damping) with profile defaults grounded
        in measured behaviour. MgO-class ionic cells get FMIXING 30 %
        with integer occupations; smearing is never selected automatically
        for an insulating profile because it can converge a wrong-energy
        metallic basin.
        ``"off"``/``"none"`` keeps the plain defaults. Omitted
        (``None``): auto applies **only when no explicit convergence
        knob is given** -- any explicit ``damping=`` /
        ``fmixing_percent=`` / ``fock_mixing=`` / ``level_shift=`` / smearing input
        switches to fully-manual mode and nothing is auto-filled.
        Either way the ``.out`` file carries a "Convergence strategy"
        block stating the mode (AUTO default / AUTO requested / manual
        / off), the per-knob values with their provenance, and the
        classification reasons. Explicit knobs are never overridden.
        Scope: applied on ``jk_method="bipole"`` and ``"gdf"``;
        other routes run with plain defaults and label the block
        accordingly.
    use_diis, damping, fmixing_percent, fock_mixing, diis_start_iter,
    diis_subspace_size, max_iter, conv_tol_energy
        SCF controls forwarded to the periodic driver.
        ``fmixing_percent`` mirrors CRYSTAL's ``FMIXING`` keyword:
        the percentage of the previous Fock/KS matrix mixed into the
        matrix diagonalised on the next cycle. It is separate from
        density damping.
        ``fock_mixing`` is the same knob on the fractional 0.0-1.0
        scale; pass only one spelling.
        ``density_mixer`` / ``density_mixer_kerker`` expose the periodic
        Anderson/Broyden/Kerker API surface on the closed-shell multi-k GDF
        and RIJCOSX routes. The lower-level Ewald RKS drivers also support
        these mixers directly. Other routes fail closed on active requests
        instead of silently ignoring them or changing electrostatic gauge.
        Compact 3D closed-shell RKS/GDF SCAN/r2SCAN-family jobs with an
        explicit ``kpoints=`` mesh select Anderson density mixing by default
        unless ``density_mixer=`` or ``convergence="off"`` is set.
        ``smearing`` accepts the new :class:`vibeqc.SmearingOptions`
        surface. The legacy ``smearing_temperature`` may be a numeric
        electronic ``k_B T``
        (interpreted via ``smearing_unit``), ``"auto"``, ``"metal"``,
        ``"small-gap"``, ``"debug"``, ``"none"`` / ``"off"``, or
        ``None``. ``smearing_method`` selects ``"fermi-dirac"`` (default),
        ``"methfessel-paxton"`` or ``"marzari-vanderbilt"`` (all
        implemented). ``smearing_metallic`` and ``smearing_band_gap_hartree``
        guide the conservative ``"auto"`` guess.
    initial_guess
        ``"SAD"`` (default) or ``"HCORE"``.
    write_molden_file
        Emit ``{output}.molden`` of the Γ-point MOs (using the unit-cell
        molecule + basis as the molecular target).
    write_density
        Emit ``{output}.xsf`` with the SCF density on a primitive-cell
        grid (XSF works for any lattice; cube is orthorhombic-only).
    density_spacing_bohr
        Grid spacing for the XSF density. Default 0.2 bohr.
    solver
        Diagonalisation solver.
        ``"dense"`` (default) uses NumPy/ScaLAPACK dense eigh.
        ``"davidson"`` uses the block-Davidson iterative solver
        (``opts.use_davidson = True``). ``"lobpcg"`` also sets
        ``use_davidson = True``; the Python SCF loop then detects
        the LOBPCG preference and dispatches through the Python
        solver stack instead of dense or Davidson diagonalisation.
    hessian
        When True, compute harmonic vibrational frequencies for
        the unit-cell molecule via finite-difference Hessian.
        Frequencies and IR intensities printed to .out and
        embedded in QVF for vibe-view. Default False.
        Cost: ~6N SCF evaluations for the unit cell.
    """
    enforce_runtime_pin_from_env()

    # Fail fast on a molecular system. Passing a Molecule here used to die
    # deep in the setup with "'Molecule' object has no attribute 'lattice'";
    # say what the user actually needs instead.
    if not hasattr(system, "lattice") or not hasattr(system, "unit_cell"):
        raise TypeError(
            f"run_periodic_job: system must be a PeriodicSystem (a cell "
            f"with a lattice); got {type(system).__name__}. Molecules have "
            f"no lattice or Brillouin zone -- k-points, band structures, "
            f"and DOS meshes are periodic-boundary-condition features. Use "
            f"vibeqc.run_job(...) for molecular calculations, or build a "
            f"PeriodicSystem (set .lattice / .unit_cell / .dim) for a "
            f"crystal."
        )

    semiempirical_route = _plan_periodic_semiempirical_method(
        method,
        system,
        kpoints=kpoints,
    )
    if semiempirical_route is not None:
        return _run_periodic_semiempirical_job(
            system,
            route_plan=semiempirical_route,
            output=output,
            dry_run=dry_run,
            record_hostname=record_hostname,
            citations=citations,
            write_xyz_file=write_xyz_file,
            write_poscar_file=write_poscar_file,
            write_xsf_structure_file=write_xsf_structure_file,
            write_cif_file=write_cif_file,
            write_molden_file=write_molden_file,
            write_density=write_density,
            write_population_file=write_population_file,
            output_qvf=output_qvf,
            jk_method=jk_method,
            aux_basis=aux_basis,
            smearing=smearing,
            smearing_temperature=smearing_temperature,
            smearing_metallic=smearing_metallic,
            smearing_band_gap_hartree=smearing_band_gap_hartree,
            dispersion=dispersion,
            optimize=optimize,
            optimize_cell_requested=optimize_cell,
            hessian=hessian,
            tddft=tddft,
            coop_cohp=coop_cohp,
            band_structure=band_structure,
            qvf_wannier_centers=qvf_wannier_centers,
            kpoints=kpoints,
            checkpoint_qvf=checkpoint_qvf,
            dft_plus_u=dft_plus_u,
            atomic_spins=atomic_spins,
            spinlock=spinlock,
            read_from=read_from,
            restart_from=restart_from,
            functional=functional,
            max_iter=max_iter,
            conv_tol_energy=conv_tol_energy,
        )

    method_upper = method.upper()
    if method_upper not in ("RHF", "RKS", "UHF", "UKS"):
        raise NotImplementedError(
            f"run_periodic_job: method={method!r} not supported. "
            f"Supported: RHF, RKS, UHF, UKS."
        )
    if basis is None:
        raise TypeError(
            f"run_periodic_job: method={method_upper!r} requires a BasisSet. "
            "Pass basis=None only for basis-free periodic semiempirical "
            "methods such as 'scc_dftb' or 'pm6'."
        )
    if method_upper in ("RHF", "UHF") and functional is not None:
        raise ValueError(
            f"functional={functional!r} given but method={method_upper}; "
            f"use method='RKS' or 'UKS' for KS calculations"
        )
    if method_upper in ("RKS", "UKS") and functional is None:
        raise ValueError(f"method={method_upper!r} requires functional=...")
    # Odd electron count with a closed-shell method (RHF/RKS) can never be
    # a restricted single determinant. The individual periodic drivers each
    # reject this deep in dispatch with a terse "requires even electron
    # count"; catch it here with an actionable message (odd-Z metals like
    # Cu need spin-polarised UHF/UKS). Charged cells shift the parity, so
    # this uses the actual electron count, not sum(Z).
    if method_upper in ("RHF", "RKS"):
        try:
            _n_elec = int(system.n_electrons())
        except Exception:
            _n_elec = None
        if _n_elec is not None and _n_elec % 2 != 0:
            _ks = method_upper == "RKS"
            raise ValueError(
                f"method={method_upper!r} is a closed-shell (restricted) "
                f"method but the cell has an odd electron count ({_n_elec}); "
                "a restricted single determinant needs paired electrons. Use "
                f"method='{'UKS' if _ks else 'UHF'}' (spin-polarised) — for an "
                "open-shell metal such as an odd-Z transition metal, add "
                "atomic_spins=... to seed the moments. If the cell should be "
                "neutral and even, check the charge / composition."
            )
    if atomic_spins is not None and method_upper not in ("UHF", "UKS"):
        # ATOMSPIN is an open-shell broken-symmetry seed; reject the
        # closed-shell misuse early (before the dry-run manifest / setup).
        raise ValueError(
            "atomic_spins (ATOMSPIN) is an open-shell broken-symmetry seed; "
            f"it requires method='UHF' or 'UKS', got {method_upper!r}."
        )
    # SPINLOCK: resolve the string mode to a SpinlockMode (open-shell only).
    # The selected driver fails closed if it does not implement the mode.
    _spinlock_mode = SpinlockMode.OFF
    if spinlock is not None:
        if method_upper not in ("UHF", "UKS"):
            raise ValueError(
                "spinlock is an open-shell broken-symmetry convergence aid; "
                f"it requires method='UHF' or 'UKS', got {method_upper!r}."
            )
        _spinlock_map = {
            "spin_schedule": SpinlockMode.SPIN_SCHEDULE,
            "pattern_hold": SpinlockMode.PATTERN_HOLD,
            "off": SpinlockMode.OFF,
        }
        _spinlock_key = str(spinlock).strip().lower().replace("-", "_")
        if _spinlock_key not in _spinlock_map:
            raise ValueError(
                f"unknown spinlock={spinlock!r} (valid: 'spin_schedule', "
                "'pattern_hold')."
            )
        _spinlock_mode = _spinlock_map[_spinlock_key]
        if _spinlock_mode != SpinlockMode.OFF and int(spinlock_iterations) <= 0:
            raise ValueError(
                "spinlock requires spinlock_iterations > 0 (the number of "
                "locked/held SCF cycles before release)."
            )
    # READ restart validations, early (before the dry-run manifest / setup).
    _is_read_request = initial_guess.strip().upper() in ("READ", "MOREAD", "COREAD")
    _read_multik_request = _is_read_request and _is_multik_kpoints(kpoints)
    if read_from is not None and not _is_read_request:
        raise ValueError(
            "read_from is only used with initial_guess='read' (or 'moread' / "
            "'coread'); set initial_guess to request a READ restart."
        )
    if _read_multik_request:
        if method_upper not in ("RHF", "RKS"):
            raise NotImplementedError(
                "multi-k periodic READ restart is currently closed-shell "
                "only (RHF/RKS). Open-shell multi-k READ needs per-spin "
                "per-k Bloch densities."
            )
        if read_from is None:
            raise ValueError(
                "multi-k periodic initial_guess='read' needs an in-memory "
                "multi-k source result or a .qvf archive with all-k Bloch "
                "restart data."
            )
        if isinstance(read_from, (str, os.PathLike)):
            _read_suffix = Path(os.fspath(read_from)).suffix.lower()
            if _read_suffix != ".qvf":
                raise NotImplementedError(
                    "multi-k periodic READ restart from non-QVF file sources "
                    "is not implemented: a multi-k restart needs all per-k "
                    "complex Bloch coefficients and occupations."
                )

    # Resolve the J/K method (AUTO -> concrete pick) and validate the
    # combination of method x lattice x basis. This is intentionally
    # done before any expensive setup so user errors fire early.
    resolved_jk = pick_jk_method(
        jk_method,
        lattice=np.asarray(system.lattice, dtype=float),
        basis_name=basis.name,
        n_atoms=len(system.unit_cell),
        scf_method=method_upper,
        dim=int(system.dim),
    )
    validate_jk_method(
        resolved_jk,
        lattice=np.asarray(system.lattice, dtype=float),
        basis_name=basis.name,
    )
    _requested_bloch_kmesh = _runner_bloch_kmesh(system, kpoints)
    _requested_kmesh_size = _bloch_kmesh_size(_requested_bloch_kmesh)
    _requested_true_multik = kpoints is not None and _requested_kmesh_size > 1
    _requested_gamma_only = (
        kpoints is None
        or _gamma_kmesh_info(system, _requested_bloch_kmesh) is not None
    )
    if _read_multik_request and resolved_jk not in (
        PeriodicJKMethod.GDF,
        PeriodicJKMethod.RIJCOSX,
        PeriodicJKMethod.GPW,
        PeriodicJKMethod.GAPW,
    ):
        raise NotImplementedError(
            "multi-k periodic READ from in-memory results is wired for "
            "closed-shell GDF, RIJCOSX, GPW, and GAPW routes only; the "
            f"selected jk_method={resolved_jk.value!r} has no per-k density "
            "restart hook yet."
        )
    if resolved_jk == PeriodicJKMethod.AICCM2026DEV_B:
        from .periodic.exchange_convention import BVK_EWALD, exchange_q0_label

        _requested_b_exchange_q0 = exchange_q0_label(exchange_exxdiv)
        if _requested_b_exchange_q0 != BVK_EWALD:
            raise ValueError(
                "jk_method='aiccm2026dev-b' fixes the finite-torus "
                "exchange_q0 convention to 'bvk-ewald'; "
                f"exchange_exxdiv={exchange_exxdiv!r} requests "
                f"{_requested_b_exchange_q0!r}, which the χ-CCM-B drivers "
                "do not implement. Omit exchange_exxdiv or pass 'ewald'."
            )
        if method_upper not in ("RHF", "RKS", "UHF", "UKS"):
            raise NotImplementedError(
                "jk_method='aiccm2026dev-b' implements RHF, RKS, UHF, and "
                f"UKS; got method={method_upper!r}"
            )
        if dft_plus_u:
            raise NotImplementedError(
                "jk_method='aiccm2026dev-b' does not yet implement DFT+U"
            )
        if optimize:
            raise NotImplementedError(
                "jk_method='aiccm2026dev-b' does not yet implement analytic "
                "χ-CCM-B nuclear gradients, so periodic geometry optimization "
                "is disabled. Use aiccm2026dev_b_gradient_status(result) to "
                "inspect the declared finite-torus convention and open "
                "gradient terms."
            )
        if hessian:
            raise NotImplementedError(
                "jk_method='aiccm2026dev-b' does not yet implement χ-CCM-B "
                "force constants. The generic periodic Hessian path would "
                "differentiate a non-B unit-cell model, so it is disabled."
            )
    if resolved_jk == PeriodicJKMethod.BIPOLE and int(system.dim) < 3:
        # Fail fast at the API level: the bipolar-expansion Fock build relies on a
        # 3-D Ewald/Madelung lattice sum (probe_charge_madelung_supercell, the
        # shifted reciprocal lattice) that is undefined for a vacuum-padded 1-D/2-D
        # cell -- otherwise this surfaces only deep inside the SCF setup as a terse
        # "requires dim=3" ValueError.
        raise NotImplementedError(
            f"jk_method='bipole' requires a 3-D periodic system (dim=3); got "
            f"dim={int(system.dim)}. The bipolar-expansion Fock build uses a 3-D "
            f"Ewald/Madelung lattice sum that is not defined for a vacuum-padded "
            f"1-D/2-D cell. For 1-D/2-D use jk_method='gdf' (Γ-only real-space GDF) "
            f"or the EWALD_3D / four-center references; low-D bipole support needs a "
            f"low-dimensional Ewald sum (periodic-SCF roadmap)."
        )
    _bipole_cutoff_bohr = (
        None if bipole_cutoff_bohr is None else float(bipole_cutoff_bohr)
    )
    _bipole_nuclear_cutoff_bohr = (
        None
        if bipole_nuclear_cutoff_bohr is None
        else float(bipole_nuclear_cutoff_bohr)
    )
    if _bipole_cutoff_bohr is not None or _bipole_nuclear_cutoff_bohr is not None:
        if resolved_jk != PeriodicJKMethod.BIPOLE:
            raise NotImplementedError(
                "run_periodic_job: bipole_cutoff_bohr and "
                "bipole_nuclear_cutoff_bohr apply only to "
                "jk_method='bipole'. Use cutoff_ha for GPW/GAPW grid "
                "routes, rsgdf_ke_cutoff/rsgdf_tail_ke_cutoff for GDF, "
                "or select jk_method='bipole'."
            )
        if _bipole_cutoff_bohr is not None and _bipole_cutoff_bohr <= 0.0:
            raise ValueError(
                "run_periodic_job: bipole_cutoff_bohr must be positive "
                f"(bohr); got {bipole_cutoff_bohr!r}."
            )
        if (
            _bipole_nuclear_cutoff_bohr is not None
            and _bipole_nuclear_cutoff_bohr <= 0.0
        ):
            raise ValueError(
                "run_periodic_job: bipole_nuclear_cutoff_bohr must be "
                f"positive (bohr); got {bipole_nuclear_cutoff_bohr!r}."
            )
        if _bipole_nuclear_cutoff_bohr is None:
            _bipole_nuclear_cutoff_bohr = _bipole_cutoff_bohr
    _rijcosx_closed_shell_multik = (
        resolved_jk == PeriodicJKMethod.RIJCOSX
        and method_upper in ("RHF", "RKS")
        and _requested_true_multik
    )
    if (
        rsgdf_tail_ke_cutoff is not None
        and resolved_jk != PeriodicJKMethod.GDF
        and not _rijcosx_closed_shell_multik
    ):
        raise NotImplementedError(
            "run_periodic_job: rsgdf_tail_ke_cutoff is currently supported only "
            "with jk_method='gdf' or closed-shell true multi-k "
            "jk_method='rijcosx'. Other J/K routes would ignore the high-|G| "
            "RSGDF tail correction, so this combination fails closed."
        )
    _jk_requested_label = (
        jk_method.value
        if isinstance(jk_method, PeriodicJKMethod)
        else str(jk_method)
    ).strip().lower()
    _jk_method_explicit = _jk_requested_label not in ("auto", "")
    _executed_jk_method = resolved_jk.value
    _dft_plus_u_route = "none"
    _dftu_native_routes = (
        PeriodicJKMethod.BIPOLE,
        PeriodicJKMethod.GPW,
        PeriodicJKMethod.GAPW,
    )
    _dftu_gdf_native_route = (
        dft_plus_u
        and resolved_jk == PeriodicJKMethod.GDF
        and method_upper in ("RHF", "RKS")
        and _requested_true_multik
    )
    _dftu_rijcosx_native_route = (
        dft_plus_u
        and resolved_jk == PeriodicJKMethod.RIJCOSX
        and method_upper in ("RHF", "RKS")
        and _requested_true_multik
    )
    if (
        dft_plus_u
        and _jk_method_explicit
        and resolved_jk not in _dftu_native_routes
        and not _dftu_gdf_native_route
        and not _dftu_rijcosx_native_route
    ):
        raise NotImplementedError(
            "run_periodic_job: dft_plus_u with explicit "
            f"jk_method={_jk_requested_label!r} is not wired. Supported "
            "explicit +U routes are jk_method='bipole', Gamma "
            "jk_method='gpw' for RHF/RKS/UHF/UKS, multi-k GPW RKS, "
            "closed-shell jk_method='gapw', and closed-shell multi-k "
            "jk_method='gdf'/'rijcosx'. Omit jk_method to use the legacy "
            "Gamma/direct +U route, or choose a supported explicit backend."
        )
    # --- Cell reduction (symmetry) ------------------------------------
    system_original_info: Optional[str] = None

    # Resolve the 'symmetry' kwarg into a reduction flag.
    _sym_val = str(symmetry).lower() if isinstance(symmetry, str) else ""
    _do_reduce = (
        reduce_to_primitive or symmetry is True or _sym_val in ("auto", "reduce")
    )
    _do_attach = _do_reduce or _sym_val == "attach"

    if _do_attach and system.symmetry is None:
        attach_symmetry(system, symprec=symmetry_precision)

    if _do_reduce:
        try:
            system_prim, sg_input = _reduce_system_to_primitive(
                system,
                symprec=symmetry_precision,
            )
        except ValueError:
            # Already primitive -- just attach symmetry and continue
            if system.symmetry is None:
                attach_symmetry(system, symprec=symmetry_precision)
        else:
            if system.symmetry is None:
                attach_symmetry(system, symprec=symmetry_precision)
            system_original_info = _build_primitive_summary(
                system,
                system_prim,
                sg_input,
            )
            basis = BasisSet(system_prim.unit_cell_molecule(), basis.name)
            system = system_prim

    if fmixing_percent is not None and fock_mixing is not None:
        raise ValueError(
            "run_periodic_job: pass either fmixing_percent= or fock_mixing=, "
            "not both"
        )
    if fock_mixing is not None:
        fock_mixing_value = float(fock_mixing)
        if not (0.0 <= fock_mixing_value < 1.0):
            raise ValueError(
                "run_periodic_job: fock_mixing must be in [0, 1); "
                f"got {fock_mixing}"
            )
        fmixing_percent = 100.0 * fock_mixing_value

    fock_mixing = 0.0
    if fmixing_percent is not None:
        fock_mixing = float(fmixing_percent) / 100.0
        if not (0.0 <= fock_mixing < 1.0):
            raise ValueError(
                "run_periodic_job: fmixing_percent must be in [0, 100); "
                f"got {fmixing_percent}"
            )

    output_stem = Path(os.fspath(output))
    out_path = output_stem.with_suffix(".out")
    molden_path = output_stem.with_suffix(".molden")
    xsf_path = output_stem.with_suffix(".xsf")
    out_path.parent.mkdir(parents=True, exist_ok=True)

    if qvf_wannier_centers:
        if not output_qvf:
            raise ValueError(
                "qvf_wannier_centers=True requires output_qvf=True so the "
                "x_ccm.wannier_centers overlay has a container to write into"
            )
        if resolved_jk != PeriodicJKMethod.AICCM2026DEV_B:
            raise NotImplementedError(
                "qvf_wannier_centers=True is currently implemented only for "
                "jk_method='aiccm2026dev-b', where the finite-torus occupied "
                "localization convention is explicit"
            )

    _population_variant = (
        "bipole"
        if resolved_jk == PeriodicJKMethod.BIPOLE
        else (
            "aiccm2026dev-b"
            if resolved_jk == PeriodicJKMethod.AICCM2026DEV_B
            else "standard"
        )
    )
    real_plan = OutputPlan.from_run_job_kwargs(
        output=output_stem,
        method=method_upper,
        basis=basis.name,
        functional=functional,
        write_molden_file=write_molden_file,
        write_xyz=write_xyz_file,
        write_poscar=write_poscar_file,
        write_xsf_structure=write_xsf_structure_file,
        write_density_xsf=write_density,
        write_cif=write_cif_file,
        write_population=write_population_file,
        population_variant=_population_variant,
        citations=citations,
        crash_dump=False,
        output_qvf=output_qvf,
        job_kind="periodic_scf",
    )

    # Dry-run short-circuit (Phase O5). Mirrors the molecular run_job
    # path: build the OutputPlan from current kwargs, write a one-shot
    # ``{output}.system`` with ``[outputs].status = "dry_run"``, print
    # the declared-artefacts summary, and return None without running
    # the SCF. Honours both the ``dry_run=True`` kwarg and the
    # ``VIBEQC_DRY_RUN=1`` env var that vq's submit pre-flight sets.
    if dry_run or is_dry_run_requested():
        _estimate_bytes: Optional[int] = None
        if is_dry_run_estimate_requested():
            try:
                _estimate_bytes = _periodic_gpw_gapw_dry_run_estimate_bytes(
                    system,
                    basis,
                    resolved_jk=resolved_jk,
                    method_upper=method_upper,
                    functional=functional,
                    cutoff_ha=cutoff_ha,
                    kpoints=kpoints,
                )
            except Exception:
                _estimate_bytes = None
            if _estimate_bytes is None:
                try:
                    _gdf_memory = _periodic_gdf_estimate(
                        system,
                        basis,
                        resolved_jk=resolved_jk,
                        method_upper=method_upper,
                        functional=functional,
                        kpoints=kpoints,
                        aux_basis=aux_basis,
                    )
                    if _gdf_memory is not None:
                        _estimate_bytes = _gdf_memory.estimate.total_bytes
                except Exception:
                    _estimate_bytes = None
            if optimize:
                try:
                    _gradient_estimate_bytes = (
                        _periodic_xc_gradient_dry_run_estimate_bytes(
                            system,
                            basis,
                            method_upper=method_upper,
                            functional=functional,
                            lattice_cutoff_bohr=_bipole_cutoff_bohr,
                        )
                    )
                    if _gradient_estimate_bytes is not None:
                        _estimate_bytes = max(
                            int(_estimate_bytes or 0),
                            int(_gradient_estimate_bytes),
                        )
                except Exception:
                    pass
        dry_run_manifest(
            real_plan,
            record_hostname=record_hostname,
            estimate_bytes=_estimate_bytes,
        )
        return None

    from .periodic.exchange_convention import exchange_q0_label

    _output_writer = OutputWriter(
        real_plan,
        record_hostname=record_hostname,
        extra_run_fields={
            "exchange_q0": exchange_q0_label(exchange_exxdiv),
            "jk_method_requested": _jk_requested_label,
            "jk_method_resolved": resolved_jk.value,
            "jk_method_executed": _executed_jk_method,
            "dft_plus_u": bool(dft_plus_u),
            "dft_plus_u_route": _dft_plus_u_route,
        },
    )
    _PERIODIC_OUTPUT_WRITER.set(_output_writer)

    plog = resolve_progress(progress, verbose=verbose)

    # --- Opt-in live QVF checkpointing (vibe-view hot-reload) -----------
    # Build the checkpointer + its QVF-only plan up front. The job's real
    # OutputPlan already owns the final archive; this narrower plan is only a
    # type-gate because checkpoint sections come from per-snapshot context.
    # When enabled, wrap ``plog`` so each SCF cycle
    # that lands a ``plog.iteration(...)`` also refreshes the checkpoint
    # QVF on the configured cadence. This covers the routes that stream
    # per-iteration through the shared logger -- the Ewald, GDF, BIPOLE, and
    # GPW drivers all take ``progress=plog`` and call ``plog.iteration(...)``
    # per SCF cycle. The GPW route (``periodic_gapw_j.run_periodic_rhf_gpw`` /
    # ``run_periodic_rks_gpw_multi_k`` + the open-shell UHF/UKS/multi-k
    # siblings in ``periodic_gapw_open_shell``) now threads the same
    # ``progress=`` handle into its Python SCF loop, so GPW jobs get
    # per-iteration cadence too. Only the molecular compiled-C++ SCF still
    # runs without a per-iteration Python hook, so molecular jobs remain
    # start + terminal frames only.
    from .output.checkpoint import (
        QvfCheckpointer as _QvfCheckpointer,
        wrap_progress_for_checkpoints as _wrap_progress_for_checkpoints,
    )

    _checkpointer = _QvfCheckpointer(
        checkpoint_qvf if output_qvf else None,
        checkpoint_every,
        plan=OutputPlan.from_run_job_kwargs(
            output=output_stem,
            method=method_upper,
            basis=basis.name,
            functional=functional,
            output_qvf=True,
            job_kind="periodic_scf",
        ),
    )
    if _checkpointer.enabled and checkpoint_every > 0:

        def _periodic_scf_checkpoint(_n: int, _fields: dict) -> None:
            _checkpointer.maybe_snapshot(
                _n,
                energy_eh=_fields.get("energy"),
                system=system,
                method=method_upper,
                basis=basis.name,
                functional=functional,
            )

        plog = _wrap_progress_for_checkpoints(plog, _periodic_scf_checkpoint)

    # --- Auto-read smearing / bz_integration from KPoints metadata -----
    # When kpoints is a KPoints.recommend() result with .smearing /
    # .bz_integration set, and the user did not explicitly pass those
    # args, auto-apply them.  Explicit user args always win (mirrors the
    # periodic_convergence_auto "explicit wins" contract).
    from .kpoints import KPoints as _KPoints

    _kpts_smearing = None
    _kpts_bz = None
    _kpts_uses_ml = False
    # [routes.numerics] keys naming the published k-point constructions this
    # job used, so the references block cites them (CLAUDE.md § 8). The SCF
    # mesh contributes its own key; an attached band structure contributes its
    # path convention. Plain Monkhorst-Pack and explicit lists contribute none.
    _numerics: List[str] = []
    if isinstance(kpoints, _KPoints):
        _kpts_smearing = getattr(kpoints, "smearing", None)
        _kpts_bz = getattr(kpoints, "bz_integration", None)
        _kpts_uses_ml = getattr(kpoints, "uses_ml_predictor", False)
        _numerics.extend(getattr(kpoints, "citation_numerics", ()) or ())
    _numerics.extend(
        getattr(getattr(band_structure, "kpath", None), "citation_numerics", ()) or ()
    )
    # Resolve bz_integration: explicit arg wins, else KPoints metadata.
    if bz_integration is _BZ_INTEGRATION_UNSET:
        bz_integration = _kpts_bz  # None / "smearing" / "gilat"
    # Validate bz_integration early.
    if bz_integration is not None:
        bz_integration = str(bz_integration).strip().lower()
        if bz_integration not in ("smearing", "gilat"):
            raise ValueError(
                "run_periodic_job: bz_integration must be None, 'smearing', "
                f"or 'gilat'; got {bz_integration!r}."
            )

    # Build the option object expected by the selected native driver.
    opts = (
        PeriodicKSOptions() if method_upper in ("RKS", "UKS") else PeriodicRHFOptions()
    )
    if method_upper in ("RKS", "UKS"):
        opts.functional = str(functional)
    if resolved_jk == PeriodicJKMethod.BIPOLE:
        if sr_range_screening:
            # M4b (pair-resolved truncation): QQR-style separation-aware
            # screening for the SR erfc J/K build; cited via the
            # bipole_sr_range route below.
            opts.lattice_opts.sr_range_screening = True
        if _bipole_cutoff_bohr is not None:
            opts.lattice_opts.cutoff_bohr = _bipole_cutoff_bohr
        if _bipole_nuclear_cutoff_bohr is not None:
            opts.lattice_opts.nuclear_cutoff_bohr = _bipole_nuclear_cutoff_bohr
        elif opts.lattice_opts.nuclear_cutoff_bohr > opts.lattice_opts.cutoff_bohr:
            # BIPOLE's corrected Ewald gauge relies on neutral-cell
            # cancellation between V_ne, E_nn, and the matching electronic
            # J/K cell set.  The generic LatticeSumOptions default has a
            # longer nuclear cutoff (25 bohr) than electronic cutoff
            # (15 bohr), which overbinds molecular-limit hybrid RKS rows by
            # mHa.  Keep the public BIPOLE default coherent unless the user
            # explicitly requested a separate nuclear cutoff.
            opts.lattice_opts.nuclear_cutoff_bohr = opts.lattice_opts.cutoff_bohr
    # --- Auto-attach ECPs from pob-TZVP-REV2 CRYSTAL ECP blocks -------
    # When the basis carries inline ECP data (detected via Z+200 header
    # in CRYSTAL-format basis files), convert to libecpint inline-primitive
    # blocks and attach to the options object (Phase 14g).
    _ecp_blocks, _ecp_centers, _eff_z, _total_ncore = _resolve_ecp_data(system, basis)
    if _ecp_blocks:
        opts.ecp_primitive_blocks = _ecp_blocks
        opts.ecp_home_centers = _ecp_centers
        opts.ecp_effective_charges = _eff_z
        opts.ecp_total_ncore = _total_ncore
    opts.use_diis = bool(use_diis)
    # --- Diagonalisation solver -----------------------------------------
    if solver not in ("dense", "davidson", "lobpcg"):
        raise ValueError(
            f"run_periodic_job: solver='{solver}' is not recognised. "
            f"Supported: 'dense', 'davidson', 'lobpcg'."
        )
    if solver != "dense":
        opts.use_davidson = True
    # --- Smearing: resolve any user-engaged smearing input first ------
    smearing_engaged = (
        smearing is not None
        or smearing_temperature is not _SMEARING_UNSET
        or smearing_metallic is not None
        or smearing_band_gap_hartree is not None
    )
    # Auto-read smearing from KPoints metadata when not explicitly given.
    if not smearing_engaged and _kpts_smearing is not None:
        smearing = _kpts_smearing
        smearing_engaged = True
    _smearing_temperature_arg = (
        0.0 if smearing_temperature is _SMEARING_UNSET else smearing_temperature
    )
    if smearing is not None:
        if _smearing_temperature_arg not in (0.0, None):
            raise ValueError(
                "run_periodic_job: pass either smearing= or "
                "smearing_temperature=, not both"
            )
        if not isinstance(smearing, SmearingOptions):
            raise TypeError(
                "run_periodic_job: smearing must be a SmearingOptions instance"
            )
        if smearing.enabled and smearing.flavor not in (
            "fermi-dirac",
            "methfessel-paxton",
            "marzari-vanderbilt",
        ):
            raise NotImplementedError(
                "run_periodic_job: smearing flavor "
                f"{smearing.flavor!r} is not implemented"
            )
        smearing_temperature_hartree = float(smearing.temperature)
        smearing_method_label = smearing.flavor
        smearing_source = smearing.source
        smearing_reason = smearing.reason
    else:
        smearing_resolution = resolve_smearing_temperature(
            _smearing_temperature_arg,
            unit=smearing_unit,
            method=smearing_method,
            metallic=smearing_metallic,
            band_gap_hartree=smearing_band_gap_hartree,
            n_electrons=system.n_electrons(),
        )
        smearing_temperature_hartree = float(smearing_resolution.temperature)
        smearing_method_label = smearing_resolution.method
        smearing_source = smearing_resolution.source
        smearing_reason = smearing_resolution.reason

    # --- Automatic convergence strategy (transparency contract) -------
    # Explicit user knobs are never overridden; auto fills only unset
    # knobs, and the .out states the mode, the classification, and the
    # per-knob reasons. v1 applies auto on the BIPOLE and GDF routes --
    # other routes run with mode "off" unless the user set knobs
    # (mode "manual" labels them explicitly there too).
    _auto_supported = resolved_jk in (
        PeriodicJKMethod.BIPOLE,
        PeriodicJKMethod.GDF,
    )
    _requested_convergence = convergence
    if (
        convergence is not None
        and str(convergence).strip().lower() == "auto"
        and not _auto_supported
    ):
        _requested_convergence = "off"
    elif convergence is None and not _auto_supported:
        _requested_convergence = "off"
    convergence_strategy = resolve_convergence_strategy(
        system,
        method=method_upper,
        convergence=_requested_convergence,
        # Only the BIPOLE KS drivers carry an in-driver FMIXING-30%
        # default for DFT functionals, and since the 2026-07-13 Gap-B
        # validation only when DIIS is off (under DIIS the mixing is
        # redundant damping that measurably slows convergence -- see
        # pbc_bipole_rks.py). The floor mirrors the driver so the
        # printed strategy matches what actually runs; GDF honours the
        # resolved value verbatim, so no floor there.
        ks_driver_fock_mixing_floor=(
            resolved_jk == PeriodicJKMethod.BIPOLE and not bool(use_diis)
        ),
        explicit={
            "fock_mixing": fock_mixing if fmixing_percent is not None else None,
            "level_shift": (float(level_shift) if level_shift is not None else None),
            "damping": float(damping) if damping is not None else None,
            "smearing_temperature": (
                smearing_temperature_hartree if smearing_engaged else None
            ),
        },
    )
    _gdf_capability_knobs = dict(convergence_strategy.knobs)
    _gdf_capability_changed = False
    if resolved_jk == PeriodicJKMethod.GDF:
        # Mirror of the _gamma_default_pure_gdf_ok domain gate in the Γ GDF
        # dispatch below, minus the knob values this filter itself resolves.
        # When the default-Γ closed-shell RHF run qualifies for the
        # PySCF-µHa-validated run_pbc_gdf_rhf in every non-knob respect, an
        # AUTO-resolved fock-mixing / level-shift (e.g. the ionic-insulator
        # profile's FMIXING 30% on MgO-class cells) must not force the run
        # onto the legacy molecular-limit fallback, whose dense-core absolute
        # energies are PARITY_HELD (G-GDF-001 in HANDOVER_GATED_ITEMS.md).
        # Zero the auto knobs instead -- the parity route carries its own
        # convergence handling (DIIS + accelerators + the auto-sized
        # high-|G| tail). Explicit user knobs keep the legacy fallback.
        _gamma_pure_gdf_capable = (
            method_upper == "RHF"
            and kpoints is None
            and gdf_method is None
            and int(system.dim) == 3
            and int(system.n_electrons()) % 2 == 0
            and int(system.multiplicity) == 1
            and abs(
                float(sum(atom.Z for atom in system.unit_cell))
                - int(system.n_electrons())
            )
            <= 0.5
            and float(convergence_strategy.knobs["smearing_temperature"].value)
            <= 0.0
            and not symmetry_stabilize
            and not symmetry_reduce_fock
        )
        _level_shift_resolution = convergence_strategy.knobs["level_shift"]
        _gdf_level_shift_filter_reason = None
        if (
            _level_shift_resolution.source == "auto"
            and _level_shift_resolution.value != 0.0
        ):
            if method_upper in ("UHF", "UKS"):
                _gdf_level_shift_filter_reason = (
                    "capability-filtered to zero: open-shell GDF drivers do "
                    "not implement the level-shift operator"
                )
            elif _requested_gamma_only and gdf_method is not None:
                _gdf_level_shift_filter_reason = (
                    "capability-filtered to zero: the exact-Gamma GDF driver "
                    "selected by an explicit gdf_method does not implement "
                    "the level-shift operator"
                )
            elif _gamma_pure_gdf_capable:
                _gdf_level_shift_filter_reason = (
                    "capability-filtered to zero: the default-Gamma "
                    "closed-shell RHF run stays on the PySCF-parity "
                    "run_pbc_gdf_rhf driver, which does not implement the "
                    "level-shift operator (the legacy fallback's dense-core "
                    "absolute energies are parity-held)"
                )
        if _gdf_level_shift_filter_reason is not None:
            _gdf_capability_knobs["level_shift"] = KnobResolution(
                0.0,
                "auto",
                _gdf_level_shift_filter_reason,
            )
            _gdf_capability_changed = True

        _fock_mixing_resolution = convergence_strategy.knobs["fock_mixing"]
        _gdf_fock_mixing_filter_reason = None
        if (
            _fock_mixing_resolution.source == "auto"
            and _fock_mixing_resolution.value != 0.0
        ):
            if method_upper in ("UHF", "UKS") and kpoints is not None:
                _gdf_fock_mixing_filter_reason = (
                    "capability-filtered to zero: open-shell k-mesh GDF "
                    "drivers do not implement Fock mixing"
                )
            elif (
                method_upper == "RHF"
                and _requested_gamma_only
                and gdf_method is not None
            ):
                _gdf_fock_mixing_filter_reason = (
                    "capability-filtered to zero: the closed-shell "
                    "exact-Gamma GDF driver selected by an explicit "
                    "gdf_method does not implement Fock mixing"
                )
            elif _gamma_pure_gdf_capable:
                _gdf_fock_mixing_filter_reason = (
                    "capability-filtered to zero: the default-Gamma "
                    "closed-shell RHF run stays on the PySCF-parity "
                    "run_pbc_gdf_rhf driver, which does not implement Fock "
                    "mixing (the legacy fallback's dense-core absolute "
                    "energies are parity-held)"
                )
        if _gdf_fock_mixing_filter_reason is not None:
            _gdf_capability_knobs["fock_mixing"] = KnobResolution(
                0.0,
                "auto",
                _gdf_fock_mixing_filter_reason,
            )
            _gdf_capability_changed = True
    if _gdf_capability_changed:
        # Auto selection must respect the capabilities of the chosen driver.
        # Leave unset/default requests at zero and state why in the strategy
        # block. Explicit unsupported requests remain fail-closed errors below.
        convergence_strategy = ConvergenceStrategy(
            mode=convergence_strategy.mode,
            classification=convergence_strategy.classification,
            knobs=_gdf_capability_knobs,
        )
    _conv_auto_note = None
    if (
        convergence is not None
        and str(convergence).strip().lower() == "auto"
        and not _auto_supported
    ):
        _conv_auto_note = (
            'convergence="auto" is wired for jk_method="bipole" and '
            'jk_method="gdf" only in this version; using plain defaults '
            "for this route"
        )

    _density_mixer_auto_note = None
    _density_mixer_params_default = (
        int(density_mixer_depth) == 8
        and float(density_mixer_beta) == 0.5
        and not bool(density_mixer_kerker)
        and float(kerker_k0) == 1.5
        and float(kerker_strength) == 1.0
        and float(kerker_cutoff_ha) == 120.0
    )
    _convergence_off_requested = (
        convergence is not None
        and str(convergence).strip().lower() in ("off", "none")
    )
    if (
        density_mixer is None
        and _density_mixer_params_default
        and not _convergence_off_requested
        and resolved_jk == PeriodicJKMethod.GDF
        and method_upper == "RKS"
        and kpoints is not None
        and _functional_is_scan_family(functional)
        and int(system.dim) == 3
        and int(system.multiplicity) == 1
        and int(system.n_electrons()) % 2 == 0
    ):
        _mgga_cls = classify_periodic_system(system)
        if _mgga_cls.profile in _COMPACT_MGGA_DENSITY_MIXER_PROFILES:
            density_mixer = "anderson"
            density_mixer_beta = 0.35
            _density_mixer_auto_note = (
                "compact periodic SCAN/r2SCAN RKS/GDF profile "
                f"'{_mgga_cls.profile}' uses density_mixer='anderson' "
                "with beta=0.35 by default; the GDF driver disables "
                "Fock-DIIS and FMIXING while density mixing is active. "
                "Pass density_mixer='diis' or convergence='off' to keep "
                "the Fock-DIIS route"
            )

    opts.damping = convergence_strategy.value("damping")
    opts.fock_mixing = convergence_strategy.value("fock_mixing")
    if resolved_jk in (PeriodicJKMethod.BIPOLE, PeriodicJKMethod.GDF):
        opts.level_shift = convergence_strategy.value("level_shift")
    _smear_res = convergence_strategy.knobs["smearing_temperature"]
    if not smearing_engaged and _smear_res.source == "auto" and _smear_res.value > 0.0:
        # Auto strategy turned smearing on: surface it through the same
        # .out lines the explicit smearing path uses.
        smearing_method_label = "fermi-dirac"
        smearing_source = "auto-strategy"
        smearing_reason = _smear_res.reason
    opts.smearing_temperature = convergence_strategy.value("smearing_temperature")
    _insulator_smearing_note = insulator_smearing_warning(
        system,
        opts.smearing_temperature,
        band_gap_hartree=smearing_band_gap_hartree,
        metallic=smearing_metallic,
    )
    if _insulator_smearing_note is not None:
        warnings.warn(_insulator_smearing_note, UserWarning, stacklevel=2)
    _validate_smearing_dispatch(
        method=method_upper,
        jk_method=resolved_jk,
        smearing_temperature=opts.smearing_temperature,
        kpoints=kpoints,
    )
    _validate_density_mixer_dispatch(
        density_mixer=density_mixer,
        density_mixer_depth=density_mixer_depth,
        density_mixer_beta=density_mixer_beta,
        density_mixer_kerker=density_mixer_kerker,
        kerker_k0=kerker_k0,
        kerker_strength=kerker_strength,
        kerker_cutoff_ha=kerker_cutoff_ha,
        jk_method=resolved_jk,
        method=method_upper,
        kpoints=kpoints,
    )
    opts.diis_start_iter = int(diis_start_iter)
    opts.diis_subspace_size = int(diis_subspace_size)
    opts.max_iter = int(max_iter)
    opts.conv_tol_energy = float(conv_tol_energy)
    # pybind11 enums expose ``__members__`` rather than supporting
    # ``Enum[name]`` subscripting directly. MOREAD / COREAD are accepted
    # ORCA-style spellings of the READ restart guess.
    _guess_aliases = {"MOREAD": "READ", "COREAD": "READ"}
    _guess_name = _guess_aliases.get(initial_guess.upper(), initial_guess.upper())
    try:
        opts.initial_guess = InitialGuess.__members__[_guess_name]
    except KeyError as exc:
        valid = ", ".join(InitialGuess.__members__.keys())
        raise ValueError(
            f"unknown initial_guess={initial_guess!r} (valid: {valid})"
        ) from exc

    # ATOMSPIN broken-symmetry seed (UHF/UKS only; the closed-shell misuse is
    # rejected above). Carried on the options struct so the open-shell
    # driver's guess step assembles the per-atom broken-symmetry density. The
    # GuessEngine validates the per-atom tag count against the SCF cell (so it
    # stays correct under cell reduction).
    if atomic_spins is not None:
        opts.atomic_spins = [int(s) for s in atomic_spins]

    # SPINLOCK (UHF/UKS only; resolved + validated above). Carried on the
    # options struct so the open-shell driver runs the two-phase schedule
    # (SPIN_SCHEDULE) or MOM-holds the seeded pattern (PATTERN_HOLD).
    if _spinlock_mode != SpinlockMode.OFF:
        opts.spinlock_mode = _spinlock_mode
        opts.spinlock_value = int(spinlock_value)
        opts.spinlock_iterations = int(spinlock_iterations)

    # READ restart. Gamma restarts resolve the prior g=0 cell density and put
    # it on the options struct. Closed-shell multi-k restarts keep the native
    # per-k density list separate so the selected multi-k driver can inject
    # D(k) directly instead of degrading to a Gamma density.
    _read_density_k_closed = None
    if opts.initial_guess == InitialGuess.READ:
        _read_path = opts.read_path
        _read_obj = read_from
        if isinstance(read_from, (str, os.PathLike)):
            _read_path = os.fspath(read_from)
            _read_obj = None
        opts.read_path = _read_path

        if _read_multik_request:
            from .guess_read import resolve_periodic_read_density_k_closed

            _read_density_k_closed = resolve_periodic_read_density_k_closed(
                read_from=_read_obj,
                read_path=_read_path,
                expected_n_k=_bloch_kmesh_size(_runner_bloch_kmesh(system, kpoints)),
                n_basis=basis.nbasis,
            )
        else:
            from .guess_read import (
                resolve_periodic_read_densities_open,
                resolve_periodic_read_density_closed,
            )

            if method_upper in ("UHF", "UKS"):
                _da, _db = resolve_periodic_read_densities_open(
                    basis, read_path=_read_path, read_from=_read_obj
                )
                opts.read_density_alpha = _da
                opts.read_density_beta = _db
            else:
                opts.read_density = resolve_periodic_read_density_closed(
                    basis, read_path=_read_path, read_from=_read_obj
                )

    label = f"{method_upper}"
    if functional:
        label = f"{label} / {functional}"

    _periodic_memory = None
    try:
        if resolved_jk in (PeriodicJKMethod.GPW, PeriodicJKMethod.GAPW):
            _gpw_estimate = _periodic_gpw_gapw_estimate(
                system,
                basis,
                resolved_jk=resolved_jk,
                method_upper=method_upper,
                functional=functional,
                cutoff_ha=cutoff_ha,
                kpoints=kpoints,
            )
            if _gpw_estimate is not None:
                _periodic_memory = SimpleNamespace(
                    estimate=_gpw_estimate,
                    kind="generic",
                )
        elif resolved_jk in (PeriodicJKMethod.GDF, PeriodicJKMethod.RIJCOSX):
            _gdf_estimate = _periodic_gdf_estimate(
                system,
                basis,
                resolved_jk=resolved_jk,
                method_upper=method_upper,
                functional=functional,
                kpoints=kpoints,
                aux_basis=aux_basis,
            )
            if _gdf_estimate is not None:
                _gdf_estimate.kind = "gdf"
                _periodic_memory = _gdf_estimate
    except Exception:
        _periodic_memory = None

    t_job_start = time.perf_counter()

    hessian_result = None  # populated by hessian=True; consumed by QVF writer

    with OutputChannel.to_file(out_path):
        # --- Banner ---------------------------------------------------
        write(banner() + "\n\n")
        libs = library_versions()
        write(f"  Job: PERIODIC {label}  basis={basis.name}\n")
        write(f"  J/K method: {describe_jk_method(resolved_jk)}\n")
        if jk_method != "auto" and jk_method != PeriodicJKMethod.AUTO:
            write(f"    (user-requested: {jk_method!r})\n")
        else:
            write(f"    (resolved from AUTO)\n")
        write("\n")
        write(_system_summary(system))
        if system_original_info is not None:
            write(system_original_info)
        elif system.symmetry is not None:
            sg = system.symmetry
            write(
                f"  Symmetry: {sg.international_symbol} (No. {sg.number}), "
                f"point group {sg.point_group}, order {sg.order}\n\n"
            )
        write(_basis_summary(basis))
        write(section_header("SCF options", width=56))
        write(f"    use_diis            = {opts.use_diis}\n")
        write(f"    damping             = {opts.damping}\n")
        if fmixing_percent is not None:
            write(f"    fmixing_percent     = {float(fmixing_percent)}\n")
        if opts.smearing_temperature > 0.0 or smearing_source != "explicit":
            write(f"    smearing_method      = {smearing_method_label}\n")
            write(f"    smearing_source      = {smearing_source}\n")
            if smearing_reason:
                write(f"    smearing_reason      = {smearing_reason}\n")
            write(f"    smearing_temperature = {opts.smearing_temperature}\n")
            if opts.smearing_temperature > 0.0:
                write(
                    "    smearing_temperature_K = "
                    f"{hartree_to_kelvin_temperature(opts.smearing_temperature)}\n"
                )
        if bz_integration is not None:
            write(f"    bz_integration      = {bz_integration}\n")
        write(f"    diis_start_iter     = {opts.diis_start_iter}\n")
        write(f"    diis_subspace_size  = {opts.diis_subspace_size}\n")
        write(f"    max_iter            = {opts.max_iter}\n")
        write(f"    conv_tol_energy     = {opts.conv_tol_energy}\n")
        write(f"    initial_guess       = {initial_guess.upper()}\n")
        if opts.fock_mixing != 0.0 and fmixing_percent is None:
            write(f"    fock_mixing         = {opts.fock_mixing}\n")
        if resolved_jk == PeriodicJKMethod.BIPOLE:
            write(
                "    bipole_cutoff_bohr  = "
                f"{float(opts.lattice_opts.cutoff_bohr)}\n"
            )
            write(
                "    bipole_nuclear_cutoff_bohr = "
                f"{float(opts.lattice_opts.nuclear_cutoff_bohr)}\n"
            )
        if resolved_jk == PeriodicJKMethod.BIPOLE and (
            convergence_strategy.value("level_shift") != 0.0
        ):
            write(
                "    level_shift         = "
                f"{convergence_strategy.value('level_shift')}\n"
            )
        if resolved_jk == PeriodicJKMethod.GDF or (
            resolved_jk == PeriodicJKMethod.AICCM2026DEV_B
            and aiccm_backend.strip().lower().replace("-", "_") != "four_center"
        ):
            write(f"    aux_basis           = {aux_basis or '<auto>'}\n")
            write(f"    gdf_method          = {gdf_method or 'rsgdf'}\n")
            write(f"    rsgdf_ke_cutoff     = {float(rsgdf_ke_cutoff)}\n")
            if rsgdf_tail_ke_cutoff is not None:
                write(
                    "    rsgdf_tail_ke_cutoff = "
                    f"{float(rsgdf_tail_ke_cutoff)}\n"
                )
            write(f"    mdf_ke_cutoff       = {float(mdf_ke_cutoff)}\n")
            write(
                "    OpenMP threads      = "
                f"{get_num_threads()}  (max for native parallel regions)\n"
            )
            write("\n")
        # --- Convergence-strategy transparency block -------------------
        # States whether the convergence aids were chosen automatically
        # (by default or by convergence="auto"), set manually, or left
        # plain -- with the classification evidence for auto choices.
        write(section_header("Convergence strategy", width=56))
        for _line in convergence_strategy.log_lines():
            write(f"    {_line}\n")
        if _conv_auto_note:
            write(f"    note: {_conv_auto_note}\n")
        if _density_mixer_auto_note:
            write(f"    note: {_density_mixer_auto_note}\n")
        if _insulator_smearing_note:
            # A real warning: normalised to the canonical "  WARNING: ..."
            # surface (survives --quiet, fires a structured event) rather
            # than the 4-space in-block echo it used to be.
            warn(_insulator_smearing_note, role="insulator_smearing")
        write("\n")
        plog.info(
            "convergence strategy: "
            f"{convergence_strategy.mode}"
            + (
                f" ({convergence_strategy.classification.profile})"
                if convergence_strategy.classification is not None
                else ""
            )
        )
        if _density_mixer_auto_note:
            plog.info(_density_mixer_auto_note)
        if _insulator_smearing_note:
            plog.warn(_insulator_smearing_note)

        if _periodic_memory is not None:
            from .memory import (
                check_memory,
                check_periodic_gdf_memory,
                format_memory_report,
            )

            _estimate = _periodic_memory.estimate
            _gdf_override = bool(os.environ.get("VIBEQC_GDF_MEMORY_OVERRIDE"))
            _override_requested = memory_override or (
                getattr(_periodic_memory, "kind", "") == "gdf" and _gdf_override
            )
            write(
                "  "
                + format_memory_report(
                    _estimate,
                    override_requested=_override_requested,
                ).replace("\n", "\n  ")
                + "\n\n"
            )
            flush()
            if getattr(_periodic_memory, "kind", "") == "gdf":
                check_periodic_gdf_memory(
                    _estimate,
                    n_kpoints=int(_periodic_memory.n_kpoints),
                    route_label=str(_periodic_memory.route_label),
                    allow_exceed=_override_requested,
                )
            else:
                check_memory(_estimate, allow_exceed=memory_override)

        # --- Restart from previous GPW/GAPW calculation ---------------
        restart_density = None
        if restart_from is not None:
            from .periodic_gapw_restart import load_gpw_result

            restart_path = Path(os.fspath(restart_from))
            if not restart_path.exists():
                raise FileNotFoundError(
                    f"restart_from={restart_from!r}: file not found"
                )
            data = load_gpw_result(str(restart_path))
            kind = data.get("kind", "")
            if kind not in ("gpw_scf", "gpw_multi_k_scf"):
                raise ValueError(
                    f"restart_from={restart_from!r}: unsupported kind "
                    f"{kind!r} (expected 'gpw_scf' or 'gpw_multi_k_scf')"
                )
            restart_density = np.asarray(data["density"], dtype=float)
            plog.info(
                f"Restart from {restart_path} ({kind}, "
                f"E = {render_energy_labeled(float(data.get('energy', 0)), width=0, precision=6)})"
            )
            write(f"    restart_from        = {restart_from}\n")
            write(
                f"    restart_energy      = {render_energy_labeled(float(data.get('energy', 0)), width=20, precision=10)}\n"
            )

        # --- SCF (dispatch on resolved jk_method) --------------------
        t0 = time.perf_counter()

        plog.banner(f"run_periodic_job  PERIODIC {label}  basis={basis.name}")
        plog.info(f"Output file: {out_path}")

        # Initial checkpoint frame: the input geometry, so a live viewer
        # has a structure to show the instant the SCF starts. Per-cycle
        # cadence frames follow via the wrapped ``plog`` (if enabled).
        if _checkpointer.enabled:
            _checkpointer.snapshot(
                system=system,
                method=method_upper,
                basis=basis.name,
                functional=functional,
            )

        def _record_dft_plus_u_route(route: str, executed_jk: str) -> None:
            nonlocal _dft_plus_u_route, _executed_jk_method
            _dft_plus_u_route = route
            _executed_jk_method = executed_jk
            write(f"    dft_plus_u_route    = {route}\n")

        # --- DFT+U interception (Increment 4d) ----------------------
        # When dft_plus_u is set, route to the appropriate +U-capable
        # driver (Increment 4d-bipole for UHF/UKS; Increment 4c for
        # multi-k RKS; Γ-only for RHF).  Closed-shell true multi-k GDF is
        # native below; FFT_POISSON still has no +U Fock-build hook.
        if (
            dft_plus_u
            and resolved_jk not in _dftu_native_routes
            and not _dftu_gdf_native_route
            and not _dftu_rijcosx_native_route
        ):
            from . import HubbardSite as _HubbardSite  # noqa: F401
            from . import (
                run_rhf_periodic_gamma as _run_rhf_periodic_gamma,
            )
            from . import (
                run_rks_periodic as _run_rks_periodic,
            )

            if method_upper == "UKS":
                # Open-shell UKS +U via the BIPOLE driver (Increment
                # 4d-bipole UKS). Same per-spin pattern as UHF, plus
                # the UKS XC contribution.
                from .pbc_bipole_uks import run_pbc_bipole_uks

                kmesh = _runner_bloch_kmesh(system, kpoints)
                _record_dft_plus_u_route(
                    "legacy_auto_bipole_uks_"
                    + ("multi_k" if kpoints is not None else "gamma"),
                    "bipole",
                )
                result = run_pbc_bipole_uks(
                    system,
                    basis,
                    kmesh,
                    opts,
                    functional=functional,
                    linear_dep_threshold=1e-7,
                    use_ewald_j_split=True,
                    ewald_omega=ewald_omega,
                    ewald_precision=ewald_precision,
                    use_oda=use_oda,
                    oda_trust_lambda_max=oda_trust_lambda_max,
                    use_mom=use_mom,
                    use_multipole_far_field=use_multipole_far_field,
                    multipole_l_max=multipole_l_max,
                    use_exchange_ewald_split=use_exchange_ewald_split,
                    exchange_exxdiv=exchange_exxdiv,
                    sr_image_precision=sr_image_precision,
                    progress=plog,
                    dft_plus_u=dft_plus_u,
                )
            elif method_upper == "UHF":
                # Open-shell UHF +U via the BIPOLE driver (Increment
                # 4d-bipole). Route directly to run_pbc_bipole_uhf
                # with the user-supplied kmesh.
                from .pbc_bipole_uhf import run_pbc_bipole_uhf

                kmesh = _runner_bloch_kmesh(system, kpoints)
                _record_dft_plus_u_route(
                    "legacy_auto_bipole_uhf_"
                    + ("multi_k" if kpoints is not None else "gamma"),
                    "bipole",
                )
                result = run_pbc_bipole_uhf(
                    system,
                    basis,
                    kmesh,
                    opts,
                    linear_dep_threshold=1e-7,
                    use_ewald_j_split=True,
                    ewald_omega=ewald_omega,
                    ewald_precision=ewald_precision,
                    use_oda=use_oda,
                    oda_trust_lambda_max=oda_trust_lambda_max,
                    use_mom=use_mom,
                    use_multipole_far_field=use_multipole_far_field,
                    multipole_l_max=multipole_l_max,
                    use_exchange_ewald_split=use_exchange_ewald_split,
                    exchange_exxdiv=exchange_exxdiv,
                    sr_image_precision=sr_image_precision,
                    progress=plog,
                    dft_plus_u=dft_plus_u,
                )
            elif method_upper == "RHF":
                if kpoints is not None and tuple(kpoints) != (1, 1, 1):
                    raise NotImplementedError(
                        f"DFT+U on multi-k periodic RHF (kpoints="
                        f"{kpoints!r}) is not yet wired. Increment 4c "
                        "covers the closed-shell-DFT (RKS) multi-k "
                        "path through cpp/src/periodic_scf.cpp; for "
                        "RHF you can run Γ-only today via "
                        "kpoints=None or kpoints=(1,1,1)."
                    )
                _record_dft_plus_u_route("legacy_auto_direct_rhf_gamma", "direct")
                # PeriodicRHFOptions field-by-field copy of the salient
                # convergence knobs; the +U-via-DIRECT path uses its
                # own Coulomb backend so jk_method-specific options
                # don't apply.
                rhf_opts = PeriodicRHFOptions()
                rhf_opts.max_iter = opts.max_iter
                rhf_opts.conv_tol_energy = opts.conv_tol_energy
                rhf_opts.conv_tol_grad = float(getattr(opts, "conv_tol_grad", 1e-6))
                rhf_opts.damping = float(getattr(opts, "damping", 0.5))
                rhf_opts.use_diis = opts.use_diis
                rhf_opts.diis_start_iter = opts.diis_start_iter
                rhf_opts.diis_subspace_size = opts.diis_subspace_size
                rhf_opts.lattice_opts = opts.lattice_opts
                rhf_opts.use_davidson = getattr(opts, "use_davidson", False)
                result = _run_rhf_periodic_gamma(
                    system,
                    basis,
                    rhf_opts,
                    dft_plus_u=dft_plus_u,
                )
            else:
                # RKS via the multi-k DIRECT_TRUNCATED driver. Now
                # supports arbitrary kmesh (Increment 4c) -- uses
                # the k-averaged AO occupation matrix and adds
                # S(k) V_AO S(k) per k.
                from .kpoints import KPoints

                _record_dft_plus_u_route(
                    "legacy_auto_direct_rks_"
                    + ("gamma" if kpoints is None else "multi_k"),
                    "direct",
                )
                ks_opts = PeriodicKSOptions()
                ks_opts.functional = functional or "lda"
                ks_opts.max_iter = opts.max_iter
                ks_opts.conv_tol_energy = opts.conv_tol_energy
                ks_opts.conv_tol_grad = float(getattr(opts, "conv_tol_grad", 1e-6))
                ks_opts.damping = float(getattr(opts, "damping", 0.5))
                ks_opts.use_diis = opts.use_diis
                ks_opts.diis_start_iter = opts.diis_start_iter
                ks_opts.diis_subspace_size = opts.diis_subspace_size
                ks_opts.lattice_opts = opts.lattice_opts
                ks_opts.use_davidson = getattr(opts, "use_davidson", False)
                if kpoints is None:
                    kmesh = KPoints.monkhorst_pack(system, (1, 1, 1))
                else:
                    kp = (
                        list(kpoints)
                        if isinstance(kpoints, (list, tuple))
                        else [kpoints, kpoints, kpoints]
                    )
                    kmesh = KPoints.monkhorst_pack(system, tuple(kp))
                result = _run_rks_periodic(
                    system,
                    basis,
                    kmesh,
                    ks_opts,
                    dft_plus_u=dft_plus_u,
                )
        elif resolved_jk == PeriodicJKMethod.AICCM2026DEV_B:
            if bz_integration not in (None, "smearing"):
                raise NotImplementedError(
                    "aiccm2026dev-b uses the exact finite cyclic-group sum; "
                    "alternative Brillouin-zone integration is not applicable"
                )
            if (
                aiccm_lattice_extension is not None
                or aiccm_wigner_seitz_shells is not None
            ):
                if kpoints is not None:
                    raise ValueError(
                        "aiccm2026dev-b accepts either the real-space "
                        "aiccm_lattice_extension/aiccm_wigner_seitz_shells "
                        "control or the legacy kpoints mesh alias, not both"
                    )
                _extension = cyclic_lattice_extension(
                    system,
                    aiccm_lattice_extension,
                    wigner_seitz_shells=aiccm_wigner_seitz_shells,
                )
                aiccm_mesh = _extension.repetitions
            elif kpoints is None:
                _extension = cyclic_lattice_extension(system)
                aiccm_mesh = _extension.repetitions
            elif isinstance(kpoints, (int, list, tuple)):
                _extension = cyclic_lattice_extension(system, mesh=kpoints)
                aiccm_mesh = _extension.repetitions
            elif getattr(kpoints, "mesh", None) is not None:
                if tuple(getattr(kpoints, "shift", (0, 0, 0))) != (0, 0, 0):
                    raise ValueError(
                        "aiccm2026dev-b requires a Gamma-centred k mesh "
                        "(shift=(0,0,0)) because it represents a cyclic cluster"
                    )
                _extension = cyclic_lattice_extension(system, mesh=tuple(kpoints.mesh))
                aiccm_mesh = _extension.repetitions
            else:
                raise TypeError(
                    "aiccm2026dev-b kpoints must be a cyclic mesh size, a mesh "
                    "tuple/list, or a Gamma-centred KPoints object"
                )
            write(f"    lattice_extension   = {aiccm_mesh}\n")
            write(
                "    WS half-extent      = "
                f"{_extension.wigner_seitz_half_extent} lattice vectors\n"
            )
            write(f"    equivalent k net    = {aiccm_mesh} (Gamma-centred)\n")
            write(f"    aiccm_backend       = {aiccm_backend}\n")
            write(f"    aiccm_symmetry      = {aiccm_symmetry}\n")
            if method_upper == "RHF":
                result = run_aiccm2026dev_b_rhf(
                    system,
                    basis,
                    aiccm_mesh,
                    opts,
                    backend=aiccm_backend,
                    aux_basis=aux_basis,
                    gdf_method=(gdf_method or "rsgdf"),
                    rsgdf_ke_cutoff=rsgdf_ke_cutoff,
                    mdf_ke_cutoff=mdf_ke_cutoff,
                    fock_mixing=opts.fock_mixing,
                    symmetry_mode=aiccm_symmetry,
                    symmetry_precision=symmetry_precision,
                    symmetry_require_full_group=(aiccm_symmetry_require_full_group),
                    progress=plog,
                )
            elif method_upper == "RKS":
                result = run_aiccm2026dev_b_rks(
                    system,
                    basis,
                    functional or "pbe",
                    aiccm_mesh,
                    opts,
                    backend=aiccm_backend,
                    aux_basis=aux_basis,
                    gdf_method=(gdf_method or "rsgdf"),
                    rsgdf_ke_cutoff=rsgdf_ke_cutoff,
                    mdf_ke_cutoff=mdf_ke_cutoff,
                    fock_mixing=opts.fock_mixing,
                    symmetry_mode=aiccm_symmetry,
                    symmetry_precision=symmetry_precision,
                    symmetry_require_full_group=(aiccm_symmetry_require_full_group),
                    progress=plog,
                )
            elif method_upper == "UHF":
                result = run_aiccm2026dev_b_uhf(
                    system,
                    basis,
                    aiccm_mesh,
                    opts,
                    backend=aiccm_backend,
                    aux_basis=aux_basis,
                    gdf_method=(gdf_method or "rsgdf"),
                    rsgdf_ke_cutoff=rsgdf_ke_cutoff,
                    mdf_ke_cutoff=mdf_ke_cutoff,
                    fock_mixing=opts.fock_mixing,
                    symmetry_mode=aiccm_symmetry,
                    symmetry_precision=symmetry_precision,
                    symmetry_require_full_group=(aiccm_symmetry_require_full_group),
                    progress=plog,
                )
            else:  # UKS
                result = run_aiccm2026dev_b_uks(
                    system,
                    basis,
                    functional or "pbe",
                    aiccm_mesh,
                    opts,
                    backend=aiccm_backend,
                    aux_basis=aux_basis,
                    gdf_method=(gdf_method or "rsgdf"),
                    rsgdf_ke_cutoff=rsgdf_ke_cutoff,
                    mdf_ke_cutoff=mdf_ke_cutoff,
                    fock_mixing=opts.fock_mixing,
                    symmetry_mode=aiccm_symmetry,
                    symmetry_precision=symmetry_precision,
                    symmetry_require_full_group=(aiccm_symmetry_require_full_group),
                    progress=plog,
                )
        elif resolved_jk == PeriodicJKMethod.GDF:
            # Closed-shell multi-k GDF and the legacy closed-shell Gamma
            # driver apply the shared Saunders-Hillier operator.  The pure
            # PBC-GDF Gamma driver and all open-shell GDF drivers do not yet
            # implement it, so reject those combinations instead of silently
            # accepting and dropping an explicit or auto-resolved shift.
            if (
                float(opts.level_shift) != 0.0
                and method_upper in ("UHF", "UKS")
            ):
                raise NotImplementedError(
                    "run_periodic_job: level_shift is not implemented for "
                    f"open-shell {method_upper}/GDF drivers. Use "
                    "jk_method='bipole' "
                    "for a shifted open-shell periodic SCF, or pass "
                    "level_shift=0.0."
                )
            if (
                float(opts.level_shift) != 0.0
                and _requested_gamma_only
                and gdf_method is not None
            ):
                raise NotImplementedError(
                    "run_periodic_job: Gamma-only GDF with an explicit gdf_method "
                    f"({gdf_method!r}) does not implement level_shift. Omit "
                    "gdf_method to use the level-shift-capable closed-shell "
                    "Gamma fallback, provide a non-Gamma k-point mesh, or pass "
                    "level_shift=0.0."
                )
            if (
                float(opts.fock_mixing) != 0.0
                and method_upper in ("UHF", "UKS")
                and kpoints is not None
            ):
                raise NotImplementedError(
                    "run_periodic_job: fock_mixing is not implemented for "
                    f"open-shell {method_upper}/GDF k-mesh drivers. Omit "
                    "kpoints only if Gamma is intended, or pass "
                    "fock_mixing=0.0."
                )
            if (
                float(opts.fock_mixing) != 0.0
                and method_upper == "RHF"
                and _requested_gamma_only
                and gdf_method is not None
            ):
                raise NotImplementedError(
                    "run_periodic_job: closed-shell Gamma-only RHF/GDF with "
                    f"an explicit gdf_method ({gdf_method!r}) does not "
                    "implement fock_mixing. Omit gdf_method to use the "
                    "Fock-mixing-capable Gamma fallback, provide a non-Gamma "
                    "k-point mesh, or pass fock_mixing=0.0."
                )
            # Default-Γ closed-shell RHF (no explicit gdf_method) now routes
            # through the PySCF-µHa-validated run_pbc_gdf_rhf (exxdiv='ewald'),
            # the same driver as the explicit-gdf_method and open-shell Γ
            # UHF/UKS paths -- but ONLY when the cell is in its supported
            # domain. The legacy run_rhf_periodic_gamma_gdf (molecular limit /
            # exxdiv=None) stays the fallback for everything it cannot do:
            # RKS, dim<3, charged cells, finite-T smearing, and the symmetry /
            # Fock-mixing convergence aids it threads. The gate mirrors
            # run_pbc_gdf_rhf's own preconditions (pbc_gdf.py ~l.563-601) and
            # the knobs that driver honours (DIIS + damping, NOT fock_mixing /
            # level_shift / smearing) so we never route a cell it would reject
            # or silently drop a convergence aid. It also matches the
            # kpoints=(1,1,1) routing in run_krhf so default-Γ and an explicit
            # Γ k-mesh never disagree. Pre-2026-06-15 the default fell to the
            # legacy driver, which computes the molecular limit and disagreed
            # with PySCF / the explicit path by the finite-size Madelung shift
            # (~5.8 mHa on H2/sto-3g/12-bohr).
            _gamma_q_nuc = float(sum(atom.Z for atom in system.unit_cell))
            _gamma_n_elec = int(system.n_electrons())
            _gamma_default_pure_gdf_ok = (
                method_upper == "RHF"
                and int(system.dim) == 3
                and _gamma_n_elec % 2 == 0
                and int(system.multiplicity) == 1
                and abs(_gamma_q_nuc - _gamma_n_elec) <= 0.5
                and float(opts.smearing_temperature) <= 0.0
                and float(opts.fock_mixing) == 0.0
                and float(getattr(opts, "level_shift", 0.0)) == 0.0
                and not symmetry_stabilize
                and not symmetry_reduce_fock
            )
            # Closed-shell Γ RKS can use the native pure-GDF KS engine via the
            # spin-unrestricted implementation in its singlet limit. This is
            # the only Γ KS path here whose Hartree J is the Lpq/GDF build; the
            # legacy run_rhf_periodic_gamma_gdf fallback uses the molecular-
            # limit Ewald-J bridge and is not a PySCF-GDF parity route for
            # condensed cells. Keep the gate narrow so we never drop knobs that
            # run_pbc_gdf_uks does not honour. The Γ UKS/GDF driver applies
            # Fock mixing itself, so ionic auto profiles can remain on the
            # pure-GDF path instead of falling back to the legacy bridge.
            _gamma_default_rks_gdf_ok = (
                method_upper == "RKS"
                and int(system.dim) == 3
                and _gamma_n_elec % 2 == 0
                and int(system.multiplicity) == 1
                and abs(_gamma_q_nuc - _gamma_n_elec) <= 0.5
                and float(opts.smearing_temperature) <= 0.0
                and float(getattr(opts, "level_shift", 0.0)) == 0.0
                and not symmetry_stabilize
                and not symmetry_reduce_fock
            )
            # Multi-k GDF dispatch when a k-mesh is requested.
            if kpoints is not None:
                plog.info(f"kmesh = {kpoints}")
                write(f"    kpoints            = {kpoints}\n")
                if bz_integration is not None:
                    write(f"    bz_integration     = {bz_integration}\n")
                if method_upper == "UKS":
                    if bz_integration == "gilat":
                        raise NotImplementedError(
                            "run_periodic_job: bz_integration='gilat' is "
                            "wired for closed-shell multi-k GDF (RHF/RKS) "
                            "only; open-shell GDF needs per-spin "
                            "Gilat-Raubenheimer occupations."
                        )
                    if rsgdf_tail_ke_cutoff is not None:
                        raise NotImplementedError(
                            "run_periodic_job: multi-k UKS/GDF does not yet "
                            "accept rsgdf_tail_ke_cutoff. Use the Gamma GDF "
                            "path or omit kpoints until the open-shell "
                            "multi-k tail route is wired."
                        )
                    result = run_kuks_periodic_gdf(
                        system,
                        basis,
                        kpoints,
                        opts,
                        functional=functional,
                        aux_basis=aux_basis,
                        gdf_method=(gdf_method or "rsgdf"),
                        rsgdf_ke_cutoff=rsgdf_ke_cutoff,
                        mdf_ke_cutoff=mdf_ke_cutoff,
                        progress=plog,
                    )
                elif method_upper == "UHF":
                    if bz_integration == "gilat":
                        raise NotImplementedError(
                            "run_periodic_job: bz_integration='gilat' is "
                            "wired for closed-shell multi-k GDF (RHF/RKS) "
                            "only; open-shell GDF needs per-spin "
                            "Gilat-Raubenheimer occupations."
                        )
                    if rsgdf_tail_ke_cutoff is not None:
                        raise NotImplementedError(
                            "run_periodic_job: multi-k UHF/GDF does not yet "
                            "accept rsgdf_tail_ke_cutoff. Use the Gamma GDF "
                            "path or omit kpoints until the open-shell "
                            "multi-k tail route is wired."
                        )
                    result = run_kuhf_periodic_gdf(
                        system,
                        basis,
                        kpoints,
                        opts,
                        functional=None,
                        aux_basis=aux_basis,
                        gdf_method=(gdf_method or "rsgdf"),
                        rsgdf_ke_cutoff=rsgdf_ke_cutoff,
                        mdf_ke_cutoff=mdf_ke_cutoff,
                        progress=plog,
                    )
                elif method_upper == "RKS":
                    if dft_plus_u:
                        _record_dft_plus_u_route("gdf_rks_multi_k", "gdf")
                    result = run_krks_periodic_gdf(
                        system,
                        basis,
                        kpoints,
                        opts,
                        functional=functional,
                        aux_basis=aux_basis,
                        gdf_method=(gdf_method or "rsgdf"),
                        rsgdf_ke_cutoff=rsgdf_ke_cutoff,
                        rsgdf_tail_ke_cutoff=rsgdf_tail_ke_cutoff,
                        mdf_ke_cutoff=mdf_ke_cutoff,
                        bz_integration=bz_integration,
                        fock_mixing=opts.fock_mixing,
                        density_mixer=density_mixer,
                        density_mixer_depth=density_mixer_depth,
                        density_mixer_beta=density_mixer_beta,
                        density_mixer_kerker=density_mixer_kerker,
                        kerker_k0=kerker_k0,
                        kerker_strength=kerker_strength,
                        kerker_cutoff_ha=kerker_cutoff_ha,
                        dft_plus_u_sites=dft_plus_u,
                        initial_density_k=_read_density_k_closed,
                        progress=plog,
                    )
                else:
                    if dft_plus_u:
                        _record_dft_plus_u_route("gdf_rhf_multi_k", "gdf")
                    result = run_krhf_periodic_gdf(
                        system,
                        basis,
                        kpoints,
                        opts,
                        functional=None,
                        aux_basis=aux_basis,
                        gdf_method=(gdf_method or "rsgdf"),
                        rsgdf_ke_cutoff=rsgdf_ke_cutoff,
                        rsgdf_tail_ke_cutoff=rsgdf_tail_ke_cutoff,
                        mdf_ke_cutoff=mdf_ke_cutoff,
                        bz_integration=bz_integration,
                        fock_mixing=opts.fock_mixing,
                        density_mixer=density_mixer,
                        density_mixer_depth=density_mixer_depth,
                        density_mixer_beta=density_mixer_beta,
                        density_mixer_kerker=density_mixer_kerker,
                        kerker_k0=kerker_k0,
                        kerker_strength=kerker_strength,
                        kerker_cutoff_ha=kerker_cutoff_ha,
                        dft_plus_u_sites=dft_plus_u,
                        initial_density_k=_read_density_k_closed,
                        progress=plog,
                    )
            elif method_upper == "UKS":
                # Γ open-shell UKS -- the pure-GDF driver on the rsgdf
                # path (µHa-validated; consistent with the multi-k route).
                result = run_pbc_gdf_uks(
                    system,
                    basis,
                    opts,
                    functional=functional,
                    aux_basis=aux_basis,
                    gdf_method=(gdf_method or "rsgdf"),
                    rsgdf_ke_cutoff=rsgdf_ke_cutoff,
                    rsgdf_tail_ke_cutoff=rsgdf_tail_ke_cutoff,
                    mdf_ke_cutoff=mdf_ke_cutoff,
                    progress=plog,
                )
            elif method_upper == "UHF":
                result = run_pbc_gdf_uhf(
                    system,
                    basis,
                    opts,
                    aux_basis=aux_basis,
                    gdf_method=(gdf_method or "rsgdf"),
                    rsgdf_ke_cutoff=rsgdf_ke_cutoff,
                    rsgdf_tail_ke_cutoff=rsgdf_tail_ke_cutoff,
                    mdf_ke_cutoff=mdf_ke_cutoff,
                    progress=plog,
                )
            elif _gamma_default_rks_gdf_ok:
                result = run_pbc_gdf_uks(
                    system,
                    basis,
                    opts,
                    functional=functional,
                    aux_basis=aux_basis,
                    gdf_method=(gdf_method or "rsgdf"),
                    rsgdf_ke_cutoff=rsgdf_ke_cutoff,
                    rsgdf_tail_ke_cutoff=rsgdf_tail_ke_cutoff,
                    mdf_ke_cutoff=mdf_ke_cutoff,
                    progress=plog,
                )
            elif gdf_method is not None or _gamma_default_pure_gdf_ok:
                # Closed-shell Γ RHF -> the PySCF-µHa-validated run_pbc_gdf_rhf
                # (exxdiv='ewald'). Fires for an explicit gdf_method (e.g.
                # 'mdf') AND for the plain default (gdf_method=None) once the
                # domain gate above passes. Γ RKS with an explicit gdf_method
                # is wired through the singlet UKS/GDF path above when that
                # path can honour all requested knobs; otherwise keep the
                # historical fail-closed behaviour for explicit gdf_method.
                if method_upper == "RKS":
                    raise NotImplementedError(
                        "run_periodic_job: Γ RKS with an explicit gdf_method "
                        f"({gdf_method!r}, e.g. 'mdf') cannot be combined with "
                        "the requested convergence/symmetry options. Drop "
                        "unsupported knobs such as fmixing_percent/fock_mixing "
                        "or use method='UKS'."
                    )
                result = run_pbc_gdf_rhf(
                    system,
                    basis,
                    opts,
                    aux_basis=aux_basis,
                    gdf_method=(gdf_method or "rsgdf"),
                    rsgdf_ke_cutoff=rsgdf_ke_cutoff,
                    rsgdf_tail_ke_cutoff=rsgdf_tail_ke_cutoff,
                    mdf_ke_cutoff=mdf_ke_cutoff,
                    exxdiv="ewald",
                    progress=plog,
                )
            else:
                # Legacy Γ GDF driver (molecular limit / exxdiv=None): the
                # fallback for the cases run_pbc_gdf_rhf does not support -- RKS,
                # dim<3, charged cells, finite-T smearing, and symmetry /
                # Fock-mixing convergence aids. (Open-shell Γ already routed to
                # run_pbc_gdf_u{hf,ks} above.)
                if rsgdf_tail_ke_cutoff is not None:
                    raise NotImplementedError(
                        "run_periodic_job: rsgdf_tail_ke_cutoff requires the "
                        "pure PBC-GDF route. The requested Γ job falls back to "
                        "the legacy molecular-limit GDF driver, which has no "
                        "high-|G| tail correction."
                    )
                result = run_rhf_periodic_gamma_gdf(
                    system,
                    basis,
                    opts,
                    functional=(functional if method_upper == "RKS" else None),
                    aux_basis=aux_basis,
                    fock_mixing=opts.fock_mixing,
                    symmetry_stabilize=symmetry_stabilize,
                    symmetry_reduce_fock=symmetry_reduce_fock,
                    progress=plog,
                )
                _mark_legacy_gamma_gdf_parity_hold(result, system, plog)
        elif resolved_jk == PeriodicJKMethod.SLAB_EWALD_2D:
            # dim=2 vacuum-free slab (rigorous Parry / de Leeuw-Perram gauge).
            # Route through lattice_opts.coulomb_method=SLAB_EWALD_2D: the
            # RHF/RKS dispatchers split Gamma vs multi-k on the kmesh
            # internally; UKS uses its slab drivers directly. The slab-normal
            # axis is non-periodic, so the k-mesh along it is always a single
            # Gamma point. See handovers/HANDOVER_SLAB_EWALD_2D.md.
            from ._vibeqc_core import CoulombMethod as _CoulombMethod
            from .periodic_rhf_dispatch import run_rhf_periodic_scf as _run_rhf_slab
            from .periodic_ks_dispatch import run_rks_periodic_scf as _run_rks_slab
            from .kpoints import KPoints, as_bloch_kmesh as _as_bloch

            if int(system.dim) != 2:
                raise NotImplementedError(
                    "jk_method='slab_ewald_2d' requires a dim=2 slab; got "
                    f"dim={int(system.dim)}."
                )
            opts.lattice_opts.coulomb_method = _CoulombMethod.SLAB_EWALD_2D
            if kpoints is None:
                _kp = [1, 1, 1]
            elif isinstance(kpoints, int):
                _kp = [int(kpoints), int(kpoints), 1]
            elif isinstance(kpoints, (list, tuple)):
                _seq = list(kpoints)
                _kp = (_seq + [1, 1, 1])[:3]
                _kp[2] = 1  # slab normal is non-periodic
            else:
                _kp = None  # already a KPoints / BlochKMesh
            slab_kmesh = (
                kpoints if _kp is None
                else KPoints.monkhorst_pack(system, tuple(_kp))
            )
            if method_upper == "RHF":
                result = _run_rhf_slab(system, basis, slab_kmesh, opts, progress=plog)
            elif method_upper == "RKS":
                result = _run_rks_slab(system, basis, slab_kmesh, opts, progress=plog)
            elif method_upper == "UKS":
                from .periodic_uks_ewald import (
                    run_uks_periodic_gamma_ewald2d as _uks_gamma_slab,
                )
                from .periodic_uks_multi_k_ewald import (
                    run_uks_periodic_multi_k_ewald3d as _uks_multik_slab,
                )
                _bm = _as_bloch(slab_kmesh)
                if len(_bm.kpoints) == 1 and np.allclose(_bm.kpoints[0], 0.0):
                    result = _uks_gamma_slab(system, basis, opts, progress=plog)
                else:
                    result = _uks_multik_slab(system, basis, _bm, opts, progress=plog)
            else:
                # UHF-on-slab is a follow-on; pick_jk_method already blocks it,
                # so this is defense-in-depth (CLAUDE.md Sec. 7 fail-closed).
                raise NotImplementedError(
                    "SLAB_EWALD_2D open-shell HF (UHF) on slabs is a follow-on; "
                    "use RHF, RKS, or UKS."
                )
        elif resolved_jk == PeriodicJKMethod.FFT_POISSON:
            # Unreachable: validate_jk_method (above) raises on FFT_POISSON,
            # retired as a user route (v0.13.0) because Γ-only EWALD_3D is
            # wrong on dense ionic crystals. Kept as a defense-in-depth
            # fail-closed should a future caller bypass validation. The
            # internal Γ-only drivers (run_r{h,k}f_periodic_gamma_ewald3d)
            # remain for dilute periodic + mechanics, fail-closed on dense
            # cells (CLAUDE.md Sec.7).
            raise ValueError(
                "jk_method='fft_poisson' (Γ-only EWALD_3D) is retired "
                "(v0.13.0). Use jk_method='gdf' (default), 'bipole', or "
                "'gpw'."
            )
        elif resolved_jk == PeriodicJKMethod.RIJCOSX:
            if not _requested_true_multik:
                if method_upper != "RHF":
                    raise NotImplementedError(
                        "run_periodic_job: Gamma RIJCOSX is implemented for "
                        "RHF only. RKS/UHF/UKS RIJCOSX use the true multi-k "
                        "GDF/COSX backend; pass a mesh with at least two "
                        f"k-points (got kpoints={kpoints!r})."
                    )
                if float(getattr(opts, "smearing_temperature", 0.0) or 0.0) > 0.0:
                    raise NotImplementedError(
                        "run_periodic_job: Gamma RIJCOSX RHF uses integer "
                        "occupations. Pass a true multi-k mesh to use the "
                        "GDF/COSX smearing path."
                    )
                if density_mixer not in (None, "", "none", "diis"):
                    raise NotImplementedError(
                        "run_periodic_job: density_mixer is wired on the "
                        "true multi-k RIJCOSX route only. Pass a mesh with "
                        "at least two k-points, or use the Gamma RIJCOSX RHF "
                        "driver with DIIS."
                    )
                if opts.initial_guess == InitialGuess.READ:
                    raise NotImplementedError(
                        "run_periodic_job: READ restart is wired on the "
                        "true multi-k RIJCOSX route only. The Gamma RIJCOSX "
                        "RHF driver currently starts from its default guess."
                    )
                from .periodic_rijcosx import run_periodic_rijcosx_rhf

                result = run_periodic_rijcosx_rhf(
                    system,
                    basis,
                    opts,
                    aux_basis=aux_basis,
                    gdf_method=(gdf_method or "compcell"),
                    rsgdf_ke_cutoff=rsgdf_ke_cutoff,
                    progress=plog,
                )
            else:
                plog.info(f"kmesh = {kpoints}")
                write(f"    kpoints            = {kpoints}\n")
                write("    k_exchange         = cosx\n")
                if bz_integration is not None:
                    write(f"    bz_integration     = {bz_integration}\n")
                rijcosx_gdf_method = gdf_method or "rsgdf"
                write(f"    aux_basis           = {aux_basis or '<auto>'}\n")
                write(f"    gdf_method          = {rijcosx_gdf_method}\n")
                write(f"    rsgdf_ke_cutoff     = {float(rsgdf_ke_cutoff)}\n")
                if rsgdf_tail_ke_cutoff is not None:
                    write(
                        "    rsgdf_tail_ke_cutoff = "
                        f"{float(rsgdf_tail_ke_cutoff)}\n"
                    )
                write(f"    mdf_ke_cutoff       = {float(mdf_ke_cutoff)}\n")
                if method_upper == "UKS":
                    if bz_integration == "gilat":
                        raise NotImplementedError(
                            "run_periodic_job: bz_integration='gilat' is "
                            "wired for closed-shell multi-k GDF/RIJCOSX "
                            "(RHF/RKS) only; open-shell RIJCOSX needs "
                            "per-spin Gilat-Raubenheimer occupations."
                        )
                    if rsgdf_tail_ke_cutoff is not None:
                        raise NotImplementedError(
                            "run_periodic_job: multi-k UKS/RIJCOSX does not "
                            "yet accept rsgdf_tail_ke_cutoff. The open-shell "
                            "multi-k GDF tail route must be wired first."
                        )
                    result = run_kuks_periodic_gdf(
                        system,
                        basis,
                        kpoints,
                        opts,
                        functional=functional,
                        aux_basis=aux_basis,
                        gdf_method=rijcosx_gdf_method,
                        rsgdf_ke_cutoff=rsgdf_ke_cutoff,
                        mdf_ke_cutoff=mdf_ke_cutoff,
                        k_exchange="cosx",
                        progress=plog,
                    )
                elif method_upper == "UHF":
                    if bz_integration == "gilat":
                        raise NotImplementedError(
                            "run_periodic_job: bz_integration='gilat' is "
                            "wired for closed-shell multi-k GDF/RIJCOSX "
                            "(RHF/RKS) only; open-shell RIJCOSX needs "
                            "per-spin Gilat-Raubenheimer occupations."
                        )
                    if rsgdf_tail_ke_cutoff is not None:
                        raise NotImplementedError(
                            "run_periodic_job: multi-k UHF/RIJCOSX does not "
                            "yet accept rsgdf_tail_ke_cutoff. The open-shell "
                            "multi-k GDF tail route must be wired first."
                        )
                    result = run_kuhf_periodic_gdf(
                        system,
                        basis,
                        kpoints,
                        opts,
                        functional=None,
                        aux_basis=aux_basis,
                        gdf_method=rijcosx_gdf_method,
                        rsgdf_ke_cutoff=rsgdf_ke_cutoff,
                        mdf_ke_cutoff=mdf_ke_cutoff,
                        k_exchange="cosx",
                        progress=plog,
                    )
                elif method_upper == "RKS":
                    if dft_plus_u:
                        _record_dft_plus_u_route("rijcosx_rks_multi_k", "rijcosx")
                    result = run_krks_periodic_gdf(
                        system,
                        basis,
                        kpoints,
                        opts,
                        functional=functional,
                        aux_basis=aux_basis,
                        use_compcell=True,
                        k_exchange="cosx",
                        gdf_method=rijcosx_gdf_method,
                        rsgdf_ke_cutoff=rsgdf_ke_cutoff,
                        rsgdf_tail_ke_cutoff=rsgdf_tail_ke_cutoff,
                        mdf_ke_cutoff=mdf_ke_cutoff,
                        bz_integration=bz_integration,
                        fock_mixing=opts.fock_mixing,
                        density_mixer=density_mixer,
                        density_mixer_depth=density_mixer_depth,
                        density_mixer_beta=density_mixer_beta,
                        density_mixer_kerker=density_mixer_kerker,
                        kerker_k0=kerker_k0,
                        kerker_strength=kerker_strength,
                        kerker_cutoff_ha=kerker_cutoff_ha,
                        dft_plus_u_sites=dft_plus_u,
                        initial_density_k=_read_density_k_closed,
                        progress=plog,
                    )
                else:
                    if dft_plus_u:
                        _record_dft_plus_u_route("rijcosx_rhf_multi_k", "rijcosx")
                    result = run_krhf_periodic_gdf(
                        system,
                        basis,
                        kpoints,
                        opts,
                        functional=None,
                        aux_basis=aux_basis,
                        use_compcell=True,
                        k_exchange="cosx",
                        gdf_method=rijcosx_gdf_method,
                        rsgdf_ke_cutoff=rsgdf_ke_cutoff,
                        rsgdf_tail_ke_cutoff=rsgdf_tail_ke_cutoff,
                        mdf_ke_cutoff=mdf_ke_cutoff,
                        bz_integration=bz_integration,
                        fock_mixing=opts.fock_mixing,
                        density_mixer=density_mixer,
                        density_mixer_depth=density_mixer_depth,
                        density_mixer_beta=density_mixer_beta,
                        density_mixer_kerker=density_mixer_kerker,
                        kerker_k0=kerker_k0,
                        kerker_strength=kerker_strength,
                        kerker_cutoff_ha=kerker_cutoff_ha,
                        dft_plus_u_sites=dft_plus_u,
                        initial_density_k=_read_density_k_closed,
                        progress=plog,
                    )
        elif resolved_jk == PeriodicJKMethod.DIRECT:
            # Dispatch DIRECT through the periodic_jk_direct wrapper.
            # AUTO never picks DIRECT (it diverges on tight ionic
            # crystals); the user has explicitly opted in here.
            # NOTE: this is a placeholder dispatch -- there's no DIRECT
            # SCF driver yet that builds J + K via jk_via_direct each
            # iter. It's a wrapper one would call directly outside the
            # periodic_runner SCF loop. Track the SCF-driver port in
            # docs/design_native_gdf.md (DIRECT loop variant).
            raise NotImplementedError(
                "DIRECT SCF driver is not wired into run_periodic_job "
                "yet. Use vibeqc.periodic_jk_direct.jk_via_direct(...) "
                "for one-shot J/K builds. AUTO picks GDF for SCF; opt "
                "in to DIRECT only for vacuum-padded debug studies."
            )
        elif resolved_jk == PeriodicJKMethod.GPW:
            # M2-full / M3a / M3b / M3d / M3e GPW SCF entry.
            # Γ-only: RHF, RKS, UHF, UKS via run_periodic_rhf_gpw / _uhf_gpw / _uks_gpw.
            # Multi-k: pure-DFT RKS via run_periodic_rks_gpw_multi_k (M3e).
            if method_upper not in ("RHF", "RKS", "UHF", "UKS"):
                raise NotImplementedError(
                    f"PeriodicJKMethod.GPW currently supports RHF, "
                    f"RKS, UHF, and UKS; got method={method!r}."
                )
            from .periodic_gapw_j import run_periodic_rhf_gpw
            from .periodic_gapw_runner_adapter import (
                gpw_result_to_runner_shape,
                gpw_uhf_result_to_runner_shape,
                gpw_uks_result_to_runner_shape,
            )

            is_dft = method_upper in ("RKS", "UKS")
            xc_for_gpw = functional if is_dft else None
            if is_dft and not xc_for_gpw:
                raise ValueError(
                    f"PeriodicJKMethod.GPW + method={method_upper!r} "
                    "requires a functional= argument (e.g. 'lda', "
                    "'pbe', 'b3lyp'). Got functional=None."
                )
            def _gpw_explicit_gamma_mesh(value) -> bool:
                if value is None:
                    return False
                if isinstance(value, (int, np.integer)):
                    return int(value) == 1
                if isinstance(value, (list, tuple)):
                    return tuple(int(x) for x in value) == (1, 1, 1)
                return False

            gpw_use_gamma_branch = (
                method_upper != "RKS"
                and _gpw_explicit_gamma_mesh(kpoints)
            )
            if kpoints is not None and not gpw_use_gamma_branch:
                # Multi-k GPW: pure-DFT RKS and UKS (LDA/GGA/meta-GGA;
                # hybrids raise from the drivers). Build a BlochKMesh and
                # dispatch.
                if method_upper not in ("RKS", "UKS"):
                    raise NotImplementedError(
                        f"PeriodicJKMethod.GPW multi-k through "
                        f"run_periodic_job supports RKS and UKS (pure "
                        f"DFT); got method={method!r}. Multi-k UHF / "
                        f"hybrids need per-k exact exchange and are not "
                        f"wired."
                    )
                if method_upper == "UKS" and dft_plus_u:
                    raise NotImplementedError(
                        "run_periodic_job: dft_plus_u is not wired on the "
                        "multi-k UKS GPW driver yet. Use the Gamma UKS GPW "
                        "route or multi-k BIPOLE for open-shell +U."
                    )
                from ._vibeqc_core import monkhorst_pack as _mp
                from .periodic_gapw_j import run_periodic_rks_gpw_multi_k

                kp = (
                    list(kpoints)
                    if isinstance(kpoints, (list, tuple))
                    else [kpoints, kpoints, kpoints]
                )
                # Use symmetry-reduced k-mesh when the system has
                # symmetry attached (attach_symmetry was called
                # upstream by the reduce_to_primitive / symmetry
                # kwarg resolution).
                if system.symmetry is not None:
                    from .kpoints import KPoints

                    kmesh_obj = KPoints.monkhorst_pack(system, tuple(kp), symmetry=True)
                    kmesh = kmesh_obj.to_bloch_kmesh()
                    plog.info(
                        f"  GPW multi-k: mesh = {kp} "
                        f"(symmetry-reduced: {len(kmesh.kpoints)} k)"
                    )
                    write(
                        f"    kpoints            = {kp} "
                        f"(symmetry-reduced: {len(kmesh.kpoints)} k)\n"
                    )
                else:
                    kmesh = _mp(system, list(kp))
                    plog.info(f"  GPW multi-k: mesh = {kp}")
                    write(f"    kpoints            = {kp}\n")

                if method_upper == "UKS":
                    # Pure-DFT multi-k UKS (LDA/GGA/meta-GGA; hybrids
                    # raise in the driver). Smearing for open-shell
                    # non-GDF routes is gated upstream; +U gated above.
                    from .periodic_gapw_open_shell import (
                        run_periodic_uks_gpw_multi_k,
                    )
                    from .periodic_gapw_runner_adapter import (
                        gpw_uks_multik_result_to_runner_shape,
                    )

                    gpw_result = run_periodic_uks_gpw_multi_k(
                        system,
                        basis,
                        kmesh,
                        functional=functional,
                        cutoff_ha=cutoff_ha,
                        max_iter=opts.max_iter,
                        conv_tol_energy=opts.conv_tol_energy,
                        conv_tol_density=1e-7,
                        quiet=True,
                        progress=plog,
                    )
                    if (not opts.use_diis) or opts.damping != 0.0:
                        write(
                            "    note: use_diis/damping are not honored "
                            "by the GPW multi-k UKS driver (internal "
                            "per-k Pulay DIIS)\n"
                        )
                    result = gpw_uks_multik_result_to_runner_shape(
                        gpw_result, system, basis
                    )
                else:
                    if dft_plus_u:
                        _record_dft_plus_u_route("gpw_rks_multi_k", "gpw")
                    gpw_result = run_periodic_rks_gpw_multi_k(
                        system,
                        basis,
                        kmesh,
                        functional=functional,
                        cutoff_ha=cutoff_ha,
                        max_iter=opts.max_iter,
                        conv_tol_energy=opts.conv_tol_energy,
                        conv_tol_density=1e-7,
                        smearing_temperature=opts.smearing_temperature,
                        smearing_method=smearing_method_label,
                        fock_mixing=(
                            opts.fock_mixing if fmixing_percent is not None else None
                        ),
                        dft_plus_u_sites=dft_plus_u,
                        initial_density_k=_read_density_k_closed,
                        quiet=True,
                        progress=plog,
                    )
                    if getattr(gpw_result, "fock_mixing", 0.0) != 0.0:
                        write(
                            "    fock_mixing         = "
                            f"{float(gpw_result.fock_mixing):.3g} "
                            "(GPW multi-k)\n"
                        )
                    if (not opts.use_diis) or opts.damping != 0.0:
                        write(
                            "    note: use_diis/damping are not honored by "
                            "the GPW multi-k driver (internal scheme)\n"
                        )
                    result = _GapwMultiKRunnerProxy(gpw_result)
            elif method_upper in ("RHF", "RKS"):
                if dft_plus_u:
                    _record_dft_plus_u_route(
                        f"gpw_{method_upper.lower()}_gamma",
                        "gpw",
                    )
                gpw_result = run_periodic_rhf_gpw(
                    system,
                    basis,
                    cutoff_ha=cutoff_ha,
                    max_iter=opts.max_iter,
                    conv_tol_energy=opts.conv_tol_energy,
                    conv_tol_density=1e-7,
                    functional=xc_for_gpw,
                    use_diis=opts.use_diis,
                    damping=opts.damping,
                    diis_subspace_size=opts.diis_subspace_size,
                    diis_start_iter=opts.diis_start_iter,
                    smearing_temperature=opts.smearing_temperature,
                    smearing_method=smearing_method_label,
                    dft_plus_u_sites=dft_plus_u,
                    quiet=True,
                    progress=plog,
                    **(
                        {}
                        if restart_density is None
                        else {"initial_density": restart_density}
                    ),
                )
                result = gpw_result_to_runner_shape(
                    gpw_result,
                    system,
                    basis,
                )
            elif method_upper == "UHF":
                from .periodic_gapw_open_shell import (
                    run_periodic_uhf_gpw,
                )

                if dft_plus_u:
                    _record_dft_plus_u_route("gpw_uhf_gamma", "gpw")
                gpw_result = run_periodic_uhf_gpw(
                    system,
                    basis,
                    cutoff_ha=cutoff_ha,
                    max_iter=opts.max_iter,
                    conv_tol_energy=opts.conv_tol_energy,
                    conv_tol_density=1e-7,
                    use_diis=opts.use_diis,
                    damping=opts.damping,
                    diis_subspace_size=opts.diis_subspace_size,
                    diis_start_iter=opts.diis_start_iter,
                    dft_plus_u_sites=dft_plus_u,
                    quiet=True,
                    progress=plog,
                    **(
                        {}
                        if restart_density is None
                        else {"initial_density": restart_density}
                    ),
                )
                result = gpw_uhf_result_to_runner_shape(
                    gpw_result,
                    system,
                    basis,
                )
            else:  # method_upper == "UKS"
                from .periodic_gapw_open_shell import (
                    run_periodic_uks_gpw,
                )

                if dft_plus_u:
                    _record_dft_plus_u_route("gpw_uks_gamma", "gpw")
                gpw_result = run_periodic_uks_gpw(
                    system,
                    basis,
                    functional=xc_for_gpw,
                    cutoff_ha=cutoff_ha,
                    max_iter=opts.max_iter,
                    conv_tol_energy=opts.conv_tol_energy,
                    conv_tol_density=1e-7,
                    use_diis=opts.use_diis,
                    damping=opts.damping,
                    diis_subspace_size=opts.diis_subspace_size,
                    diis_start_iter=opts.diis_start_iter,
                    dft_plus_u_sites=dft_plus_u,
                    quiet=True,
                    progress=plog,
                    **(
                        {}
                        if restart_density is None
                        else {"initial_density": restart_density}
                    ),
                )
                result = gpw_uks_result_to_runner_shape(
                    gpw_result,
                    system,
                    basis,
                )
        elif resolved_jk == PeriodicJKMethod.GAPW:
            # M3c GAPW all-electron SCF entry. Same dispatch structure
            # as the GPW block, but uses the per-atom augmentation
            # correction for all-electron accuracy. Multi-k RKS is
            # available via :func:`run_periodic_rks_gapw_multi_k`.
            #
            # The GAPW all-electron augmentation route's two open correctness
            # bugs are FIXED (2026-06-26): the Hartree Fock is now the exact
            # derivative of the energy (analytic J = dE_H/dD), and the bonded-
            # molecule overlapping-augmentation double-count is root-caused
            # (own-atom compensator + partition-of-unity). The run-level
            # GAPWExperimentalWarning that flagged those bugs is therefore
            # retired for jk_method='gapw'. The remaining caveats are accuracy,
            # not correctness -- the per-atom augmentation has a grid-convergent
            # absolute residual (hundreds of mHa for 2nd-row atoms at the
            # default N=24 grid), and pure GAPW-HF on light elements carries the
            # intrinsic soft-basis bond residual -- and are documented in
            # docs/user_guide/gapw.md (the .out "(experimental)" label stays).
            if method_upper not in ("RHF", "RKS", "UHF", "UKS"):
                raise NotImplementedError(
                    f"PeriodicJKMethod.GAPW currently supports RHF, "
                    f"RKS, UHF, and UKS; got method={method!r}."
                )
            def _gapw_explicit_gamma_mesh(value) -> bool:
                if value is None:
                    return False
                if isinstance(value, (int, np.integer)):
                    return int(value) == 1
                if isinstance(value, (list, tuple)):
                    return tuple(int(x) for x in value) == (1, 1, 1)
                return False

            gapw_use_gamma_branch = (
                method_upper != "RKS"
                and _gapw_explicit_gamma_mesh(kpoints)
            )
            if kpoints is not None and not gapw_use_gamma_branch:
                # Multi-k GAPW: pure DFT only (RKS).
                if method_upper != "RKS":
                    raise NotImplementedError(
                        f"PeriodicJKMethod.GAPW multi-k only supports "
                        f"RKS (pure DFT) at v0.10.x; hybrid / HF "
                        f"multi-k is not wired. Got method={method!r}."
                    )
                from ._vibeqc_core import monkhorst_pack as _mp
                from .periodic_gapw_augment import (
                    run_periodic_rks_gapw_multi_k,
                )

                kp = (
                    list(kpoints)
                    if isinstance(kpoints, (list, tuple))
                    else [kpoints, kpoints, kpoints]
                )
                # Use symmetry-reduced k-mesh when the system has
                # symmetry attached.
                if system.symmetry is not None:
                    from .kpoints import KPoints

                    kmesh_obj = KPoints.monkhorst_pack(system, tuple(kp), symmetry=True)
                    kmesh = kmesh_obj.to_bloch_kmesh()
                    plog.info(
                        f"  GAPW multi-k: mesh = {kp} "
                        f"(symmetry-reduced: {len(kmesh.kpoints)} k)"
                    )
                    write(
                        f"    kpoints            = {kp} "
                        f"(symmetry-reduced: {len(kmesh.kpoints)} k)\n"
                    )
                else:
                    kmesh = _mp(system, list(kp))
                    plog.info(f"  GAPW multi-k: mesh = {kp}")
                    write(f"    kpoints            = {kp}\n")

                if dft_plus_u:
                    _record_dft_plus_u_route("gapw_rks_multi_k", "gapw")
                gapw_result = run_periodic_rks_gapw_multi_k(
                    system,
                    basis,
                    kmesh,
                    functional=functional,
                    cutoff_ha=cutoff_ha,
                    max_iter=opts.max_iter,
                    conv_tol_energy=opts.conv_tol_energy,
                    conv_tol_density=1e-7,
                    dft_plus_u_sites=dft_plus_u,
                    initial_density_k=_read_density_k_closed,
                    quiet=True,
                )
                # Adapt the multi-k result to the runner duck-type.
                # GpwMultiKScfResult carries per-k data under
                # ``mo_coeffs_k`` / ``mo_energies_k`` / ``occupations_k``;
                # the runner output code looks for ``mo_coeffs``,
                # ``mo_energies``, and ``occupations`` as per-k lists.
                result = _GapwMultiKRunnerProxy(gapw_result)
            else:
                from .periodic_gapw_augment import (
                    run_periodic_rhf_gapw,
                    run_periodic_rks_gapw,
                    run_periodic_uhf_gapw,
                    run_periodic_uks_gapw,
                )
                from .periodic_gapw_runner_adapter import (
                    gpw_result_to_runner_shape,
                    gpw_uhf_result_to_runner_shape,
                    gpw_uks_result_to_runner_shape,
                )

                is_dft = method_upper in ("RKS", "UKS")
                xc_for_gapw = functional if is_dft else None
                if is_dft and not xc_for_gapw:
                    raise ValueError(
                        f"PeriodicJKMethod.GAPW + method={method_upper!r} "
                        "requires a functional= argument (e.g. 'lda', "
                        "'pbe', 'b3lyp'). Got functional=None."
                    )
                if method_upper in ("RHF", "RKS"):
                    if dft_plus_u:
                        _record_dft_plus_u_route(
                            f"gapw_{method_upper.lower()}_gamma",
                            "gapw",
                        )
                    gapw_result = run_periodic_rhf_gapw(
                        system,
                        basis,
                        cutoff_ha=cutoff_ha,
                        max_iter=opts.max_iter,
                        conv_tol_energy=opts.conv_tol_energy,
                        conv_tol_density=1e-7,
                        functional=xc_for_gapw,
                        use_diis=opts.use_diis,
                        damping=opts.damping,
                        diis_subspace_size=opts.diis_subspace_size,
                        diis_start_iter=opts.diis_start_iter,
                        smearing_temperature=opts.smearing_temperature,
                        smearing_method=smearing_method_label,
                        dft_plus_u_sites=dft_plus_u,
                        quiet=True,
                        **(
                            {}
                            if restart_density is None
                            else {"initial_density": restart_density}
                        ),
                    )
                    result = gpw_result_to_runner_shape(
                        gapw_result,
                        system,
                        basis,
                    )
                elif method_upper == "UHF":
                    if dft_plus_u:
                        _record_dft_plus_u_route("gapw_uhf_gamma", "gapw")
                    gapw_result = run_periodic_uhf_gapw(
                        system,
                        basis,
                        cutoff_ha=cutoff_ha,
                        max_iter=opts.max_iter,
                        conv_tol_energy=opts.conv_tol_energy,
                        conv_tol_density=1e-7,
                        use_diis=opts.use_diis,
                        damping=opts.damping,
                        diis_subspace_size=opts.diis_subspace_size,
                        diis_start_iter=opts.diis_start_iter,
                        dft_plus_u_sites=dft_plus_u,
                        quiet=True,
                        **(
                            {}
                            if restart_density is None
                            else {"initial_density": restart_density}
                        ),
                    )
                    result = gpw_uhf_result_to_runner_shape(
                        gapw_result,
                        system,
                        basis,
                    )
                else:  # method_upper == "UKS"
                    if dft_plus_u:
                        _record_dft_plus_u_route("gapw_uks_gamma", "gapw")
                    gapw_result = run_periodic_uks_gapw(
                        system,
                        basis,
                        functional=xc_for_gapw,
                        cutoff_ha=cutoff_ha,
                        max_iter=opts.max_iter,
                        conv_tol_energy=opts.conv_tol_energy,
                        conv_tol_density=1e-7,
                        use_diis=opts.use_diis,
                        damping=opts.damping,
                        diis_subspace_size=opts.diis_subspace_size,
                        diis_start_iter=opts.diis_start_iter,
                        dft_plus_u_sites=dft_plus_u,
                        quiet=True,
                        **(
                            {}
                            if restart_density is None
                            else {"initial_density": restart_density}
                        ),
                    )
                    result = gpw_uks_result_to_runner_shape(
                        gapw_result,
                        system,
                        basis,
                    )
        elif resolved_jk == PeriodicJKMethod.BIPOLE:
            # CRYSTAL-gauge Ewald J-split -- all four method flavours.
            # Use user-provided k-mesh if given, else default to Γ-only
            kmesh = _runner_bloch_kmesh(system, kpoints)
            if kpoints is not None:
                plog.info(f"  BIPOLE kmesh: {kpoints}")
            if bz_integration == "gilat" and method_upper != "RKS":
                raise NotImplementedError(
                    "run_periodic_job: bz_integration='gilat' is wired for "
                    "BIPOLE RKS only; BIPOLE RHF/UHF/UKS need route-specific "
                    "fractional/per-spin Gilat-Raubenheimer density support."
                )
            if bz_integration is not None and kpoints is not None:
                write(f"    bz_integration     = {bz_integration}\n")
            if dft_plus_u:
                _record_dft_plus_u_route(
                    f"bipole_{method_upper.lower()}_"
                    + ("multi_k" if kpoints is not None else "gamma"),
                    "bipole",
                )

            if method_upper == "RHF":
                from .pbc_bipole import run_pbc_bipole_rhf

                result = run_pbc_bipole_rhf(
                    system,
                    basis,
                    kmesh,
                    opts,
                    linear_dep_threshold=1e-7,
                    use_ewald_j_split=True,
                    ewald_omega=ewald_omega,
                    ewald_precision=ewald_precision,
                    use_oda=use_oda,
                    oda_trust_lambda_max=oda_trust_lambda_max,
                    use_mom=use_mom,
                    use_multipole_far_field=use_multipole_far_field,
                    multipole_l_max=multipole_l_max,
                    use_exchange_ewald_split=use_exchange_ewald_split,
                    exchange_exxdiv=exchange_exxdiv,
                    sr_image_precision=sr_image_precision,
                    progress=plog,
                    dft_plus_u=dft_plus_u,
                    bz_integration=bz_integration,
                )
            elif method_upper == "UHF":
                from .pbc_bipole_uhf import run_pbc_bipole_uhf

                result = run_pbc_bipole_uhf(
                    system,
                    basis,
                    kmesh,
                    opts,
                    linear_dep_threshold=1e-7,
                    use_ewald_j_split=True,
                    ewald_omega=ewald_omega,
                    ewald_precision=ewald_precision,
                    use_oda=use_oda,
                    oda_trust_lambda_max=oda_trust_lambda_max,
                    use_mom=use_mom,
                    use_multipole_far_field=use_multipole_far_field,
                    multipole_l_max=multipole_l_max,
                    use_exchange_ewald_split=use_exchange_ewald_split,
                    exchange_exxdiv=exchange_exxdiv,
                    sr_image_precision=sr_image_precision,
                    progress=plog,
                    dft_plus_u=dft_plus_u,
                )
            elif method_upper == "RKS":
                from .pbc_bipole_rks import run_pbc_bipole_rks

                result = run_pbc_bipole_rks(
                    system,
                    basis,
                    kmesh,
                    opts,
                    functional=functional,
                    linear_dep_threshold=1e-7,
                    use_ewald_j_split=True,
                    ewald_omega=ewald_omega,
                    ewald_precision=ewald_precision,
                    use_oda=use_oda,
                    oda_trust_lambda_max=oda_trust_lambda_max,
                    use_mom=use_mom,
                    use_multipole_far_field=use_multipole_far_field,
                    multipole_l_max=multipole_l_max,
                    use_exchange_ewald_split=use_exchange_ewald_split,
                    exchange_exxdiv=exchange_exxdiv,
                    sr_image_precision=sr_image_precision,
                    progress=plog,
                    dft_plus_u=dft_plus_u,
                    # RKS is the one BIPOLE driver implementing Gilat-net BZ
                    # integration (the guard above rejects gilat elsewhere);
                    # without this forwarding a requested "gilat" was recorded
                    # in the .out but silently ignored by the SCF.
                    bz_integration=bz_integration,
                )
            elif method_upper == "UKS":
                from .pbc_bipole_uks import run_pbc_bipole_uks

                result = run_pbc_bipole_uks(
                    system,
                    basis,
                    kmesh,
                    opts,
                    functional=functional,
                    linear_dep_threshold=1e-7,
                    use_ewald_j_split=True,
                    ewald_omega=ewald_omega,
                    ewald_precision=ewald_precision,
                    use_oda=use_oda,
                    oda_trust_lambda_max=oda_trust_lambda_max,
                    use_mom=use_mom,
                    use_multipole_far_field=use_multipole_far_field,
                    multipole_l_max=multipole_l_max,
                    use_exchange_ewald_split=use_exchange_ewald_split,
                    exchange_exxdiv=exchange_exxdiv,
                    sr_image_precision=sr_image_precision,
                    progress=plog,
                    dft_plus_u=dft_plus_u,
                )
            else:
                # Defensive: method_upper is validated to RHF/RKS/UHF/UKS at
                # the top of run_periodic_job and never reassigned, so this is
                # unreachable today. The guard keeps the BIPOLE branch self-
                # consistent with the GPW/GAPW branches (which guard locally)
                # so widening the allowed-method set upstream can never silently
                # fall through to an UnboundLocalError on `result` below.
                raise NotImplementedError(
                    "run_periodic_job: PeriodicJKMethod.BIPOLE dispatch: "
                    f"method={method_upper!r} is not supported "
                    "(expected RHF, RKS, UHF, or UKS)."
                )
        else:
            raise NotImplementedError(
                f"Periodic JK method {resolved_jk.value!r} dispatch "
                f"is not yet implemented in v0.7.1-spike."
            )
        t_scf = time.perf_counter() - t0

        # --- Write SCF trace + energies ------------------------------
        # Shared with the molecular path: iteration table, the
        # "converged in N iterations" line, and the energy-component
        # breakdown all come from format_scf_trace so periodic and
        # molecular .out files read identically. Periodic-only smearing
        # quantities follow. ``label`` (RHF/RKS/...) is already in the
        # "Job: PERIODIC <label>" header above, so it is not repeated here.
        write_scf_trace(
            result,
            include_banner=False,
            include_properties=False,
            energy_label="Total energy",
            trailing="\n\n",
        )
        write(_smearing_summary(result))
        write(_band_summary(result))
        write(_mo_summary(result))
        if resolved_jk == PeriodicJKMethod.AICCM2026DEV_B:
            _aiccm_diag = result.aiccm2026dev_b
            write(section_header("AICCM2026DEV-B invariants", width=56))
            write(f"    cyclic mesh         = {_aiccm_diag.mesh}\n")
            write(f"    cyclic cells        = {_aiccm_diag.n_cyclic_cells}\n")
            write(f"    electronic method   = {_aiccm_diag.electronic_method}\n")
            write(f"    integral backend    = {_aiccm_diag.backend}\n")
            write(f"    CCM approach        = {_aiccm_diag.ccm_approach}\n")
            write(f"    CCM construction    = {_aiccm_diag.ccm_construction}\n")
            write(
                "    evaluation repr.    = "
                f"{_aiccm_diag.evaluation_representation}\n"
            )
            write(
                "    coulomb kernel      = "
                f"{getattr(_aiccm_diag, 'coulomb_kernel', 'not-recorded')}\n"
            )
            write(
                "    exchange q=0        = "
                f"{getattr(_aiccm_diag, 'exchange_q0', 'not-recorded')}\n"
            )
            write(
                "    q=0 applicability   = "
                f"{getattr(_aiccm_diag, 'exchange_q0_applicability', 'not-recorded')}\n"
            )
            write(
                "    boundary model      = "
                f"{getattr(_aiccm_diag, 'boundary_model', 'not-recorded')}\n"
            )
            _aiccm_convention = getattr(_aiccm_diag, "finite_torus_convention", None)
            if _aiccm_convention is not None:
                write(
                    "    BvK Madelung cell   = "
                    f"{_aiccm_convention.bvk_madelung_supercell_repetitions}\n"
                )
            write(
                "    WS partition error  = "
                f"{_aiccm_diag.wigner_seitz_partition_error:.3e}\n"
            )
            write(
                "    idempotency error   = "
                f"{_aiccm_diag.density_idempotency_error:.3e}\n"
            )
            write(
                f"    electron-count err  = {_aiccm_diag.electron_count_error:.3e}\n"
            )
            write(
                "    inverse-Bloch Im    = "
                f"{_aiccm_diag.inverse_bloch_imaginary_residual:.3e}\n\n"
            )
            _symmetry_diag = getattr(result, "aiccm2026dev_b_symmetry", None)
            if _symmetry_diag is not None:
                _symmetry_plan = _symmetry_diag.plan
                write(section_header("AICCM2026DEV-B space-group diagnostic", width=56))
                write(
                    "    space group         = "
                    f"{_symmetry_plan.international_symbol} "
                    f"(No. {_symmetry_plan.space_group_number})\n"
                )
                write(
                    "    compatible ops      = "
                    f"{_symmetry_plan.n_operations_compatible}/"
                    f"{_symmetry_plan.n_operations_full}\n"
                )
                write(
                    "    full/irreducible k  = "
                    f"{_symmetry_plan.n_kpoints_full}/"
                    f"{_symmetry_plan.n_kpoints_irreducible}\n"
                )
                write("    acceleration        = diagnostic only (not applied)\n")
                write(
                    "    shell pairs         = "
                    f"{_symmetry_diag.n_unique_shell_pairs}/"
                    f"{_symmetry_diag.n_shell_pairs} unique\n"
                )
                if _symmetry_diag.n_shell_quartets is not None:
                    write(
                        "    shell quartets      = "
                        f"{_symmetry_diag.n_unique_shell_quartets}/"
                        f"{_symmetry_diag.n_shell_quartets} unique\n"
                    )
                else:
                    write("    shell quartets      = count skipped (>24 shells)\n")
                if _symmetry_diag.gamma_fock_residual is not None:
                    write(
                        "    Gamma Fock residual = "
                        f"{_symmetry_diag.gamma_fock_residual:.3e}\n"
                    )
                if _symmetry_diag.gamma_density_residual is not None:
                    write(
                        "    Gamma dens residual = "
                        f"{_symmetry_diag.gamma_density_residual:.3e}\n"
                    )
                write("\n")

        # --- Convergence-strategy post-SCF check ----------------------
        # The pre-SCF classification is a guess; the converged spectrum
        # is the truth. On auto-mode runs, compare and say so when they
        # disagree (generalizes the GDF conducting-state warning).
        if convergence_strategy.mode.startswith("auto") and bool(
            getattr(result, "converged", False)
        ):
            from .periodic_convergence_auto import (
                converged_gap_hartree,
                post_scf_profile_check,
            )

            _n_elec_chk = int(round(system.n_electrons()))
            _mult_chk = int(getattr(system, "multiplicity", 1) or 1)
            _n_a_chk = (_n_elec_chk + (_mult_chk - 1)) // 2
            _n_b_chk = _n_elec_chk - _n_a_chk
            _gap_chk = converged_gap_hartree(
                result,
                n_alpha=(
                    _n_a_chk if method_upper in ("UHF", "UKS") else _n_elec_chk // 2
                ),
                n_beta=_n_b_chk if method_upper in ("UHF", "UKS") else None,
            )
            _profile_check = post_scf_profile_check(convergence_strategy, _gap_chk)
            if _profile_check is not None:
                _lvl, _msg = _profile_check
                write(section_header("Convergence strategy check", width=56))
                write(f"    {_lvl}: {_msg}\n\n")
                if _lvl == "warning":
                    warnings.warn(
                        f"run_periodic_job: {_msg}",
                        UserWarning,
                        stacklevel=2,
                    )
                else:
                    plog.info(f"convergence strategy check: {_msg}")

        # --- Post-SCF dispersion (D3-BJ) -----------------------------
        # Mirrors the molecular runner's _DispersionAugmented wrapping
        # (see python/vibeqc/runner.py): the SCF result is left
        # untouched (energy = pure SCF) and a wrapper exposes the
        # dispersion piece via .e_dispersion / .energy_total. The
        # periodic lattice sum is handled by
        # vibeqc.dispersion_periodic.compute_d3bj_periodic.
        if dispersion is not None and dispersion is not False:
            from .dispersion_periodic import compute_d3bj_periodic
            from .runner import _DispersionAugmented, _resolve_dispersion

            d3_params = _resolve_dispersion(
                dispersion,
                functional if method_upper in ("RKS", "UKS") else None,
            )
            if d3_params is not None:
                disp = compute_d3bj_periodic(
                    system,
                    d3_params,
                    cutoff_bohr=dispersion_cutoff_bohr,
                    backend=dispersion_backend,
                )
                e_scf = float(getattr(result, "energy", 0.0))
                e_total = e_scf + float(disp.energy)
                write(
                    "\n  Dispersion correction (D3-BJ, periodic)\n"
                    "  " + "-" * 52 + "\n"
                    f"  {'backend':>10s} {disp.backend!s:>14s}\n"
                    f"  {'supercell':>10s} {disp.supercell!s:>14s}\n"
                    f"  {'s6':>10s} {d3_params.s6:14.6f}\n"
                    f"  {'s8':>10s} {d3_params.s8:14.6f}\n"
                    f"  {'a1':>10s} {d3_params.a1:14.6f}\n"
                    f"  {'a2':>10s} {d3_params.a2:14.6f}\n"
                    f"  {'E_disp':>10s} {render_energy_labeled(disp.energy, width=14, precision=8)}"
                    f"  ({disp.energy * 627.5094740631:+.4f} kcal/mol)\n"
                    f"  {'E_SCF':>10s} {render_energy_labeled(e_scf, width=14, precision=8)}\n"
                    f"  {'E_total':>10s} {render_energy_labeled(e_total, width=14, precision=8)}\n"
                )
                plog.info(
                    f"  E_disp = {disp.energy:+.6e} Ha "
                    f"({disp.backend}, supercell={disp.supercell})"
                )
                result = _DispersionAugmented(result, float(disp.energy), d3_params)

        # --- Molden ---------------------------------------------------
        _qvf_wf = None
        _qvf_bloch_wf = None
        if write_molden_file and _has_valid_mo_coeffs(result):
            try:
                from .runner import _DispersionAugmented  # noqa: F401

                mol = system.unit_cell_molecule()
                # Multi-k: use Γ-point (k=0) MOs for molden output.
                molden_result = result
                if isinstance(result.mo_coeffs, (list, tuple)):
                    molden_result = _gamma_proxy_for_multi_k(result)
                _output_writer.dispatch_role(
                    "orbitals",
                    only_format="molden",
                    result=molden_result,
                    basis=basis,
                    molecule=mol,
                    title=output_stem.name,
                    raise_on_error=True,
                )
                write(f"  Molecular orbitals written to {molden_path.name}\n")
            except Exception as exc:
                warn(
                    f"molden write failed: {type(exc).__name__}: {exc}",
                    role="molden",
                )
                warn_writer_failure(
                    exc,
                    molden_path,
                    role="molden",
                    category=OutputFailureKind.optional_artifact,
                    writer=_output_writer,
                )
        if output_qvf and _has_valid_mo_coeffs(result):
            # Periodic QVF consumers render density/orbitals from precomputed
            # torus-periodic grids. A molecular ``wavefunction.gto`` section
            # would be evaluated without image sums and clips face-straddling
            # crystalline orbitals, so the periodic runner deliberately omits
            # it. Molden output above remains the text/orbital export path.
            _qvf_wf = None
            if (
                method_upper in ("RHF", "RKS")
                and resolved_jk != PeriodicJKMethod.AICCM2026DEV_B
            ):
                try:
                    from vibeqc.output.formats.qvf import qvf_bloch_wf_data

                    if isinstance(getattr(result, "mo_coeffs", None), (list, tuple)):
                        _kpts_cart = _result_kpoints_cart(result)
                        if _kpts_cart is None:
                            raise ValueError(
                                "multi-k QVF READ payload requires k-point metadata"
                            )
                        _kpts_frac = _wrap_reciprocal_fractional(
                            _fractional_kpoints_for_output(system, _kpts_cart)
                        )
                    else:
                        _kpts_frac = np.zeros((1, 3), dtype=float)
                    _qvf_bloch_wf = qvf_bloch_wf_data(
                        result,
                        basis,
                        system.unit_cell_molecule(),
                        k_points=_kpts_frac,
                    )
                except Exception as _qvf_bloch_exc:
                    warn_output_failure(
                        _qvf_bloch_exc,
                        output_stem.with_suffix(".qvf"),
                        role="qvf_bloch_wavefunction_prep",
                        category=OutputFailureKind.compatibility_fallback,
                    )

        # --- Population dump (Phase O6, periodic) ---------------------
        _qvf_pop = None
        _qvf_bond_orders = None
        _qvf_dipole = None
        if write_population_file and _has_valid_mo_coeffs(result):
            try:
                from ._vibeqc_core import Molecule as _Mol  # noqa: F401

                mol_p = system.unit_cell_molecule()
                pop_result = result
                if resolved_jk == PeriodicJKMethod.BIPOLE:
                    from vibeqc.output.formats.population import (
                        compute_bipole_population_summary,
                        unsupported_population_summary,
                    )

                    try:
                        _bipole_pop_kmesh = (
                            kmesh
                            if "kmesh" in locals()
                            else _runner_bloch_kmesh(system, kpoints)
                        )
                        _qvf_pop = compute_bipole_population_summary(
                            result,
                            basis,
                            mol_p,
                            system,
                            lattice_options=opts.lattice_opts,
                            kmesh=_bipole_pop_kmesh,
                        )
                    except Exception as _bipole_pop_exc:
                        warn_output_failure(
                            _bipole_pop_exc,
                            output_stem.with_suffix(".population.txt"),
                            role="bipole_population_summary",
                            category=OutputFailureKind.compatibility_fallback,
                        )
                        _qvf_pop = unsupported_population_summary(
                            "periodic BIPOLE population properties are "
                            "unavailable because lattice-summed Mulliken "
                            "evaluation failed"
                        )
                elif resolved_jk == PeriodicJKMethod.AICCM2026DEV_B:
                    from vibeqc.output.formats.population import (
                        compute_aiccm2026dev_b_population_summary,
                        unsupported_population_summary,
                    )

                    try:
                        _qvf_pop = compute_aiccm2026dev_b_population_summary(
                            result,
                            basis,
                            mol_p,
                            system,
                        )
                    except Exception as _aiccm_pop_exc:
                        warn_output_failure(
                            _aiccm_pop_exc,
                            output_stem.with_suffix(".population.txt"),
                            role="aiccm2026dev_b_population_summary",
                            category=OutputFailureKind.compatibility_fallback,
                        )
                        _qvf_pop = unsupported_population_summary(
                            "periodic AICCM2026DEV-B population properties "
                            "are unavailable because finite-torus population "
                            "evaluation failed"
                        )
                else:
                    if isinstance(result.mo_coeffs, (list, tuple)):
                        pop_result = _gamma_proxy_for_multi_k(result)
                _output_writer.dispatch_role(
                    "population",
                    result=pop_result,
                    basis=basis,
                    molecule=mol_p,
                    population_summary=_qvf_pop,
                    raise_on_error=True,
                )
                write(
                    f"  Population dump written to "
                    f"{output_stem.name}.population.txt + "
                    f"{output_stem.name}.population.json\n"
                )
                # Also compute summary for QVF atom_properties section.
                if output_qvf:
                    try:
                        if _qvf_pop is None:
                            from vibeqc.output.formats.population import (
                                compute_population_summary,
                            )

                            _qvf_pop = compute_population_summary(
                                pop_result,
                                basis,
                                mol_p,
                            )
                        # Extract bond-order table for bond_orders QVF section.
                        if _qvf_pop.mayer_bonds:
                            import math

                            _pos = [(a.xyz[0], a.xyz[1], a.xyz[2]) for a in mol_p.atoms]
                            _qvf_bond_orders = {
                                "method": "mayer",
                                "pairs": [
                                    {
                                        "i": int(i),
                                        "j": int(j),
                                        "order": float(order),
                                        "symbol_i": si,
                                        "symbol_j": sj,
                                        "distance_ang": float(
                                            math.dist(_pos[i], _pos[j]) * 0.529177210903
                                        ),
                                    }
                                    for (i, j, si, sj, order) in _qvf_pop.mayer_bonds
                                ],
                            }
                        # Extract dipole moment for manifest root metadata.
                        if _qvf_pop.dipole is not None:
                            _qvf_dipole = {
                                "total_debye": float(_qvf_pop.dipole["total_debye"]),
                                "vector_debye": [
                                    float(_qvf_pop.dipole["x_ebohr"]) * 2.541746473,
                                    float(_qvf_pop.dipole["y_ebohr"]) * 2.541746473,
                                    float(_qvf_pop.dipole["z_ebohr"]) * 2.541746473,
                                ],
                                "origin": str(
                                    _qvf_pop.dipole.get("origin_bohr", "origin")
                                ),
                            }
                    except Exception as _qvf_pop_exc:
                        warn_output_failure(
                            _qvf_pop_exc,
                            output_stem.with_suffix(".qvf"),
                            role="qvf_population_prep",
                            category=OutputFailureKind.compatibility_fallback,
                        )
            except Exception as exc:
                warn(
                    f"population dump failed: {type(exc).__name__}: {exc}",
                    role="population",
                )
                warn_writer_failure(
                    exc,
                    output_stem.with_suffix(".population.txt"),
                    role="population_summary",
                    category=OutputFailureKind.optional_artifact,
                    writer=_output_writer,
                )

        # --- Geometry siblings (Phase O5) -----------------------------
        # Extended-XYZ (lattice in comment line), POSCAR, and the
        # structure-only XSF -- every viewer / chem-toolkit / fellow
        # QC code can consume at least one of these. Default-on with
        # opt-out via the matching kwarg; failures are best-effort
        # warnings so a finished SCF never gets dragged down by a
        # geometry writer.
        if write_xyz_file:
            xyz_path = output_stem.with_suffix(".xyz")
            try:
                energy_ha = float(getattr(result, "energy", float("nan")))
                if energy_ha != energy_ha:  # NaN sentinel
                    energy_ha = None
                _output_writer.dispatch_role(
                    "geometry",
                    only_format="extended-xyz",
                    system=system,
                    energy_ha=energy_ha,
                    raise_on_error=True,
                )
                write(
                    f"  Final geometry written to {xyz_path.name} "
                    f"(extended XYZ with lattice)\n"
                )
            except Exception as exc:
                warn(
                    f"extended-xyz write failed: {type(exc).__name__}: {exc}",
                    role="extended_xyz",
                )
                warn_writer_failure(
                    exc,
                    xyz_path,
                    role="extended_xyz",
                    category=OutputFailureKind.optional_artifact,
                    writer=_output_writer,
                )

        if write_poscar_file:
            poscar_path = output_stem.with_suffix(".POSCAR")
            try:
                _output_writer.dispatch_role(
                    "geometry",
                    only_format="poscar",
                    system=system,
                    comment=f"vibe-qc periodic {label} basis={basis.name}",
                    raise_on_error=True,
                )
                write(f"  Structure written to {poscar_path.name} (VASP-5 POSCAR)\n")
            except Exception as exc:
                warn(
                    f"POSCAR write failed: {type(exc).__name__}: {exc}",
                    role="poscar",
                )
                warn_writer_failure(
                    exc,
                    poscar_path,
                    role="poscar",
                    category=OutputFailureKind.optional_artifact,
                    writer=_output_writer,
                )

        if write_cif_file:
            cif_path = output_stem.with_suffix(".cif")
            try:
                _output_writer.dispatch_role(
                    "geometry",
                    only_format="cif",
                    system=system,
                    comment=f"vibe-qc periodic {label} basis={basis.name}",
                    raise_on_error=True,
                )
                write(f"  Structure written to {cif_path.name} (CIF, P 1)\n")
            except Exception as exc:
                warn(f"CIF write failed: {type(exc).__name__}: {exc}", role="cif")
                warn_writer_failure(
                    exc,
                    cif_path,
                    role="cif",
                    category=OutputFailureKind.optional_artifact,
                    writer=_output_writer,
                )

        # --- Citations (Phase O5b) ------------------------------------
        # Assemble the citation list for this periodic job (software +
        # libint + libxc-if-DFT + basis-set + functional + DIIS +
        # spglib + ECP-if-used) and emit {stem}.bibtex / .references
        # siblings plus a "## References" block in the .out so the
        # text log is self-contained. Failures are non-fatal -- the
        # SCF result is the load-bearing artefact, citations are
        # observability.
        cite_block_text: str | None = None
        # Full assembled provenance rows destined for the .system
        # manifest's [citations] section (CLAUDE.md Sec.8.3/Sec.8.5).
        # The upfront OutputWriter receives these once assembly completes.
        cite_manifest_rows: list[dict[str, Any]] = []
        if citations:
            try:
                _db = load_default_database()
                # FFT-Poisson backend uses FFTW3 (fftw_plan_* /
                # fftw_execute in cpp/src/fft_poisson.cpp). The native
                # EWALD_3D ``fft_poisson`` method now defaults to the
                # analytical AO-pair-FT Hartree J path; it only exercises
                # FFTW3 when the legacy grid backend is explicitly
                # restored for diagnostics.
                _j_ewald_backend = os.environ.get(
                    "VIBEQC_J_EWALD3D_BACKEND", "analytic_ft"
                ).lower()
                _uses_fftw = resolved_jk in (
                    PeriodicJKMethod.GPW,
                    PeriodicJKMethod.GAPW,
                ) or (
                    resolved_jk == PeriodicJKMethod.FFT_POISSON
                    and _j_ewald_backend == "grid"
                )
                _uses_ewald_ao_ft = (
                    resolved_jk == PeriodicJKMethod.FFT_POISSON
                    and _j_ewald_backend != "grid"
                )
                _uses_gpw = resolved_jk == PeriodicJKMethod.GPW
                _uses_gapw = resolved_jk == PeriodicJKMethod.GAPW
                # GDF citation gating. The Sun-Berkelbach 2017 + MD78
                # + HJO00 stack fires whenever the GDF path runs. The
                # rsgdf sub-route in the database is reserved for the
                # GDF chat to extend with Ye-Berkelbach 2021 -- until
                # then it carries the same row as ``gdf`` and the runner
                # always uses the ``gdf`` flag (the periodic-runner
                # surface doesn't currently expose which gdf_method ran).
                _aiccm_backend_name = aiccm_backend.strip().lower().replace("-", "_")
                _uses_gdf = resolved_jk == PeriodicJKMethod.GDF or (
                    resolved_jk == PeriodicJKMethod.AICCM2026DEV_B
                    and _aiccm_backend_name != "four_center"
                )
                _uses_bipole = resolved_jk == PeriodicJKMethod.BIPOLE or (
                    resolved_jk == PeriodicJKMethod.AICCM2026DEV_B
                    and _aiccm_backend_name == "four_center"
                )
                _aiccm_acceleration = (
                    ("rijcosx",)
                    if resolved_jk == PeriodicJKMethod.AICCM2026DEV_B
                    and _aiccm_backend_name == "rijcosx"
                    else ()
                )
                # SCF initial guess. Periodic runs default to SAD
                # (superposition of atomic densities, van Lenthe 2006);
                # map the validated guess name to its citation route. AUTO /
                # HCORE carry no defining-paper route.
                _periodic_guess_routes = {
                    "sad": "sad",
                    "sap": "sap",
                    "hueckel": "huckel",
                    "huckel": "huckel",
                }
                _scf_guess = _periodic_guess_routes.get(initial_guess.strip().lower())
                # Smearing flavour. Only Fermi-Dirac ships at present (a
                # non-FD flavour raises NotImplementedError above, before any
                # SCF runs), so this pass-through fires nothing today -- but it
                # routes the Methfessel-Paxton / Marzari-Vanderbilt defining
                # papers automatically once those broadenings land. Fermi-Dirac
                # is skipped by assemble(); its Mermin citation comes from
                # uses_smearing. Hyphenated flavour names normalise to the
                # underscored route keys (e.g. "marzari-vanderbilt").
                _smearing_method = None
                if opts.smearing_temperature > 0.0:
                    _smearing_method = (
                        str(smearing_method_label).strip().lower().replace("-", "_")
                    )
                # COOP/COHP bonding analysis (Hughbanks-Hoffmann COOP,
                # Dronskowski-Bloechl COHP, LOBSTER AO projection). The
                # analysis itself runs in the DOS/QVF block further down
                # (it needs the lattice matrices assembled there), so gate
                # the citation on the exact conditions that gate that
                # block: coop_cohp requested, QVF output on, SCF
                # converged. Both routes fire together -- the runner
                # always passes H_terms, so COHP is computed alongside
                # COOP whenever the analysis runs.
                _props: list[str] = []
                if coop_cohp and output_qvf and result.converged:
                    _props = ["coop", "cohp"]
                _refs = _db.assemble(
                    method=(
                        "aiccm2026dev-b"
                        if resolved_jk == PeriodicJKMethod.AICCM2026DEV_B
                        else method_upper
                    ),
                    basis=basis.name,
                    functional=functional,
                    periodic=True,  # => spglib fires
                    uses_ecp=bool(_ecp_blocks),
                    uses_fftw_poisson=_uses_fftw,
                    uses_smearing=opts.smearing_temperature > 0.0,
                    # Saunders-Hillier level shift (constant, warm-up, or an
                    # explicit level_shift_schedule) is a cited convergence
                    # technique.
                    uses_level_shift=(
                        float(getattr(opts, "level_shift", 0.0) or 0.0) != 0.0
                        or any(
                            float(s) != 0.0
                            for s in (getattr(opts, "level_shift_schedule", None) or [])
                        )
                    ),
                    smearing_method=_smearing_method,
                    scf_guess=_scf_guess,
                    dft_plus_u=bool(dft_plus_u),
                    uses_slab_ewald_2d=(
                        resolved_jk == PeriodicJKMethod.SLAB_EWALD_2D
                    ),
                    uses_gpw=_uses_gpw,
                    uses_gapw=_uses_gapw,
                    uses_gdf=_uses_gdf,
                    uses_bipole=_uses_bipole,
                    # M4b QQR-style separation-aware SR screening. M5
                    # enables it with the default padded image domain.
                    uses_bipole_sr_range=(
                        _uses_bipole
                        and bool(
                            getattr(
                                getattr(opts, "lattice_opts", None),
                                "sr_range_screening",
                                False,
                            )
                        )
                    ),
                    uses_ewald_ao_ft=_uses_ewald_ao_ft,
                    uses_ml_kpredictor=_kpts_uses_ml,
                    acceleration=_aiccm_acceleration,
                    properties=_props,
                    numerics=_numerics,
                )
                cite_manifest_rows = citation_manifest_rows(_refs)
                _output_writer.dispatch_role(
                    "citations",
                    citations=_refs,
                    raise_on_error=True,
                )
                cite_block_text = format_references_block(_refs)
                bibtex_path = output_stem.with_suffix(".bibtex")
                write(
                    f"  Citations written to {bibtex_path.name} + "
                    f"{output_stem.with_suffix('.references').name}\n"
                )
            except Exception as exc:
                warn(
                    f"citation emission failed: {type(exc).__name__}: {exc}",
                    role="citations",
                )
                warn_writer_failure(
                    exc,
                    output_stem.with_suffix(".bibtex"),
                    role="citations",
                    category=OutputFailureKind.optional_artifact,
                    writer=_output_writer,
                )

        if write_xsf_structure_file:
            xsf_struct_path = output_stem.with_suffix(".xsf")
            try:
                # When write_density=True we'd collide with the
                # density XSF below; route the structure-only XSF to
                # a different suffix in that case.
                if write_density:
                    xsf_struct_path = output_stem.with_suffix(
                        ".structure.xsf",
                    )
                _output_writer.dispatch_role(
                    "geometry",
                    only_format="xsf",
                    system=system,
                    raise_on_error=True,
                )
                write(
                    f"  Structure written to {xsf_struct_path.name} "
                    f"(XSF crystal block)\n"
                )
            except Exception as exc:
                warn(
                    f"XSF structure write failed: {type(exc).__name__}: {exc}",
                    role="xsf_structure",
                )
                warn_writer_failure(
                    exc,
                    xsf_struct_path,
                    role="xsf_structure",
                    category=OutputFailureKind.optional_artifact,
                    writer=_output_writer,
                )

        # --- Density XSF + grid for QVF -------------------------------
        # Compute the primitive-cell density grid once; reuse it for
        # both the XSF writer and the QVF archive (when requested).
        # The grid origin is (0,0,0) and the span is the lattice-
        # vector transpose (rows = spanning vectors), matching the XSF
        # convention.  All units are bohr -- write_qvf stores them
        # as-is in the grid descriptor (see _grid_descriptor).
        _qvf_volume_data = None
        _qvf_system = system
        _qvf_mo_data = None
        _qvf_extensions = None
        _qvf_vendor_json_sections = None
        if output_qvf and resolved_jk == PeriodicJKMethod.AICCM2026DEV_B:
            _qvf_vendor_json_sections = _aiccm_b_qvf_vendor_sections(result)
            if _qvf_vendor_json_sections:
                _qvf_extensions = _qvf_extensions_with(
                    _qvf_extensions,
                    "x_vibeqc",
                )
            if qvf_wannier_centers:
                try:
                    if hasattr(result, "mo_coeffs_alpha"):
                        from .periodic_aiccm2026dev_b_localization import (
                            localize_aiccm2026dev_b_unrestricted_occupied,
                        )

                        localization = localize_aiccm2026dev_b_unrestricted_occupied(
                            result,
                            system,
                            basis,
                        )
                    else:
                        from .periodic_aiccm2026dev_b_localization import (
                            localize_aiccm2026dev_b_occupied,
                        )

                        localization = localize_aiccm2026dev_b_occupied(
                            result,
                            system,
                            basis,
                        )
                    center_sections = _aiccm_b_qvf_wannier_center_sections(
                        localization,
                    )
                    if center_sections:
                        _qvf_vendor_json_sections = list(
                            _qvf_vendor_json_sections or []
                        )
                        _qvf_vendor_json_sections.extend(center_sections)
                        _qvf_extensions = _qvf_extensions_with(
                            _qvf_extensions,
                            "x_ccm",
                        )
                except Exception as _qvf_wannier_exc:
                    warn(
                        "χ-CCM-B QVF Wannier-centre overlay failed: "
                        f"{type(_qvf_wannier_exc).__name__}: {_qvf_wannier_exc}",
                        role="qvf_wannier_overlay",
                    )
                    warn_output_failure(
                        _qvf_wannier_exc,
                        output_stem.with_suffix(".qvf"),
                        role="qvf_wannier_centers",
                        category=OutputFailureKind.optional_artifact,
                    )
        _needs_density_grid = write_density or output_qvf
        if _needs_density_grid:
            try:
                from ._vibeqc_core import (
                    CoulombMethod,
                    LatticeSumOptions,
                )
                from .periodic_density import evaluate_periodic_density_on_grid
                lat_opts = LatticeSumOptions()
                lat_opts.coulomb_method = CoulombMethod.EWALD_3D
                lat_opts.cutoff_bohr = 18.0
                L_bohr = np.asarray(system.lattice, dtype=float)
                density_write_system = system
                shape = tuple(
                    max(
                        1,
                        int(
                            np.ceil(np.linalg.norm(L_bohr[:, i]) / density_spacing_bohr)
                        ),
                    )
                    for i in range(3)
                )
                if resolved_jk == PeriodicJKMethod.AICCM2026DEV_B:
                    diag = getattr(result, "aiccm2026dev_b", None)
                    mesh = tuple(
                        int(value)
                        for value in getattr(diag, "mesh", (1, 1, 1))
                    )
                    super_system = _aiccm_b_qvf_supercell_system(system, mesh)
                    super_basis = BasisSet(
                        super_system.unit_cell_molecule(),
                        basis.name,
                    )
                    L_bohr = np.asarray(super_system.lattice, dtype=float)
                    shape = tuple(
                        max(
                            1,
                            int(
                                np.ceil(
                                    np.linalg.norm(L_bohr[:, i])
                                    / density_spacing_bohr
                                )
                            ),
                        )
                        for i in range(3)
                    )
                    density_matrix = _aiccm_b_full_density_matrix_for_qvf(
                        result,
                        system,
                        mesh,
                    )
                    if density_matrix.shape != (
                        int(super_basis.nbasis),
                        int(super_basis.nbasis),
                    ):
                        raise ValueError(
                            "aiccm2026dev-b QVF supercell density shape "
                            f"{density_matrix.shape!r} does not match "
                            f"supercell basis size {int(super_basis.nbasis)}"
                        )
                    rho, per_voxel = _evaluate_density_matrix_on_lattice_grid(
                        density_matrix,
                        super_basis,
                        L_bohr,
                        shape,
                        system=super_system,
                    )
                    _qvf_system = super_system
                    density_write_system = super_system
                    if output_qvf:
                        try:
                            _qvf_mo_data = _aiccm_b_qvf_gamma_orbital_grids(
                                result,
                                super_basis,
                                super_system,
                                L_bohr,
                                shape,
                                mesh,
                            )
                        except Exception as _qvf_mo_exc:
                            warn(
                                "periodic QVF orbital grid preparation "
                                f"failed: {type(_qvf_mo_exc).__name__}: "
                                f"{_qvf_mo_exc}",
                                role="qvf_orbital_grid_prep",
                            )
                            warn_output_failure(
                                _qvf_mo_exc,
                                output_stem.with_suffix(".qvf"),
                                role="qvf_orbital_grid_prep",
                                category=OutputFailureKind.compatibility_fallback,
                            )
                else:
                    if hasattr(result, "density_alpha") and hasattr(
                        result, "density_beta"
                    ):
                        D_alpha_set = _density_lattice_set_for_output(
                            basis,
                            system,
                            _density_proxy_with_k_metadata(
                                result,
                                result.density_alpha,
                            ),
                            lat_opts,
                        )
                        D_beta_set = _density_lattice_set_for_output(
                            basis,
                            system,
                            _density_proxy_with_k_metadata(
                                result,
                                result.density_beta,
                            ),
                            lat_opts,
                        )
                        D_set = _sum_lattice_density_sets_for_output(
                            basis,
                            system,
                            lat_opts,
                            D_alpha_set,
                            D_beta_set,
                            label="periodic spin density",
                        )
                    else:
                        D_set = _density_lattice_set_for_output(
                            basis,
                            system,
                            result,
                            lat_opts,
                        )

                    rho, _ = evaluate_periodic_density_on_grid(
                        basis,
                        system,
                        D_set,
                        grid_shape=shape,
                        spacing_bohr=density_spacing_bohr,
                    )
                    per_voxel = L_bohr.T / np.array(shape, dtype=float)

                # --- XSF output ---
                if write_density:
                    _output_writer.dispatch_role(
                        "density",
                        only_format="xsf",
                        system=density_write_system,
                        data=rho,
                        name=f"{output_stem.name}_density",
                        raise_on_error=True,
                    )
                    write(f"  Density written to {xsf_path.name}\n")

                # --- Package for QVF ----------------------------------
                # QVF writes grid in bohr.  Origin is (0,0,0) for a
                # primitive-cell grid.  voxel_vectors are per-voxel
                # step vectors: lattice column / shape_i.
                if output_qvf:
                    _qvf_volume_data = {
                        "Electron density": (
                            rho,
                            np.zeros(3, dtype=float),
                            per_voxel,
                        ),
                    }
            except Exception as exc:
                warn(
                    f"density grid evaluation failed: {type(exc).__name__}: {exc}",
                    role="density_grid",
                )
                if write_density:
                    warn_writer_failure(
                        exc,
                        xsf_path,
                        role="density_grid",
                        category=OutputFailureKind.optional_artifact,
                        writer=_output_writer,
                    )
                else:
                    warn_output_failure(
                        exc,
                        output_stem.with_suffix(".qvf"),
                        role="density_grid",
                        category=OutputFailureKind.compatibility_fallback,
                    )
        # --- Timing ---------------------------------------------------
        t_total = time.perf_counter() - t_job_start
        n_iter = int(getattr(result, "n_iter", 0))
        iter_avg = (t_scf / n_iter) if n_iter > 0 else float("nan")
        write("\n  Timings (wall clock, seconds)\n")
        write("  " + "-" * 56 + "\n")
        write(f"    {'SCF total':<28s} {render_duration(t_scf)}\n")
        write(f"    {'SCF avg per iter':<28s} {render_duration(iter_avg)}  ({n_iter} iters)\n")
        write(f"    {'Job total':<28s} {render_duration(t_total)}\n")
        flush()

        # --- TD-DFT excited states (Gamma-point, TDA) -------------------
        if tddft and _has_valid_mo_coeffs(result):
            try:
                from vibeqc.tddft import (
                    run_tddft_tda_periodic as _run_td_periodic,
                )

                # Extract n_occ from the result.
                _n_occ = None
                _occ = getattr(result, "occupations", None)
                if _occ is not None:
                    if isinstance(_occ, (list, tuple)):
                        _occ = _occ[0]
                    _occ = np.asarray(_occ, dtype=float)
                    _n_occ = int(np.sum(_occ > 1.0 - 1e-8))
                if _n_occ is None:
                    _n_el = getattr(result, "n_electrons", 0)
                    _n_occ = _n_el // 2 if _n_el else 1

                _td = _run_td_periodic(
                    result,
                    basis,
                    n_occ=_n_occ,
                    n_states=tddft_n_states,
                    functional=functional,
                )
                _td_func = f" ({functional})" if functional else ""
                write(
                    f"\n  ## TD-DFT excited states"
                    f" (TDA{_td_func}, Gamma-point)\n"
                    f"  {'─' * 50}\n"
                )
                write(
                    f"  {'State':>6s}  {'E (eV)':>10s}"
                    f"  {'λ (nm)':>10s}  {'f_osc':>10s}"
                    f"  {'Dominant transition'}\n"
                )
                for _st in _td.states:
                    _dom_str = ", ".join(
                        f"{occ}→{virt} ({abs(amp):.3f})"
                        for occ, virt, amp in _st.dominant_amplitudes[:3]
                    )
                    write(
                        f"  {_st.index:>6d}  "
                        f"{_st.excitation_energy_ev:>10.4f}  "
                        f"{_st.wavelength_nm:>10.1f}  "
                        f"{_st.oscillator_strength:>10.4f}  "
                        f"{_dom_str}\n"
                    )
                write("\n")
                flush()
            except Exception as _td_per_exc:
                write(
                    f"\n  ## TD-DFT excited states\n"
                    f"  {'─' * 50}\n"
                    f"  FAILED: {type(_td_per_exc).__name__}:"
                    f" {_td_per_exc}\n"
                )
                flush()

        # --- References block (Phase O5b) -----------------------------
        # Embed the assembled citation list in the .out so the text
        # log is self-contained -- a user reading the .out doesn't have
        # to chase the .bibtex / .references siblings to know what to
        # cite. Same content that lives in the .references file, hard-
        # wrapped to match the SCF-trace layout. Emitted through the
        # citation printer's channel writer (a first-class logger op).
        if cite_block_text:
            write_references_block(block=cite_block_text)
            flush()

        # --- Harmonic vibrational analysis (finite-difference Hessian) ---
        if hessian:
            _scf_converged = bool(getattr(result, "converged", False))
            if not _scf_converged:
                write(
                    "\n  ## Vibrational Frequencies\n"
                    "  " + "-" * 52 + "\n"
                    "  SKIPPED -- SCF did not converge.\n"
                )
                flush()
            else:
                try:
                    from ._vibeqc_core import (
                        RHFOptions,
                        RKSOptions,
                        UHFOptions,
                        UKSOptions,
                    )
                    from .hessian import (
                        HessianFDOptions,
                        compute_hessian_fd,
                        ir_intensities,
                    )

                    if method_upper == "RHF":
                        _hess_scf_opts = RHFOptions()
                    elif method_upper == "UHF":
                        _hess_scf_opts = UHFOptions()
                    elif method_upper == "RKS":
                        _hess_scf_opts = RKSOptions()
                        _hess_scf_opts.functional = str(functional)
                    elif method_upper == "UKS":
                        _hess_scf_opts = UKSOptions()
                        _hess_scf_opts.functional = str(functional)
                    else:
                        raise ValueError(
                            f"Hessian not available for method={method_upper}"
                        )

                    _hess_opts = HessianFDOptions(
                        include_dipole_derivatives=True,
                        frozen_indices=hessian_frozen_indices,
                    )
                    _uc_mol = system.unit_cell_molecule()
                    hessian_result = compute_hessian_fd(
                        _uc_mol,
                        basis.name,
                        method=method_upper,
                        scf_options=_hess_scf_opts,
                        hessian_options=_hess_opts,
                    )
                    write(
                        f"\n  ## Vibrational Frequencies\n"
                        f"  " + "-" * 52 + "\n"
                        f"  Finite-difference Hessian (unit cell)"
                        f"  (step = {_hess_opts.step_bohr:.3f} bohr,"
                        f"  {hessian_result.n_displacements} displacements)\n"
                        f"  Imaginary modes: {hessian_result.imaginary_count}\n"
                        f"  Linear molecule: {hessian_result.is_linear}\n\n"
                    )

                    _n_skip = 5 if hessian_result.is_linear else 6
                    _freqs = hessian_result.frequencies_cm1
                    _ir = None
                    try:
                        _ir = ir_intensities(hessian_result)
                    except Exception as _ir_exc:
                        write(
                            f"  (warning: IR intensities not available: "
                            f"{type(_ir_exc).__name__}: {_ir_exc})\n"
                        )

                    _header = "  {:<6s}  {:>10s}" + (
                        "  {:>12s}" if _ir is not None else ""
                    )
                    if _ir is not None:
                        write(
                            _header.format("Mode", "Freq/cm\u207b\u00b9", "IR/(km/mol)")
                            + "\n"
                        )
                        write("  " + "-" * 52 + "\n")
                        for k in range(_n_skip, len(_freqs)):
                            _label = f"{k - _n_skip + 1}"
                            _freq = _freqs[k]
                            _iri = _ir[k]
                            if _freq < 0:
                                _freq_str = f"{abs(_freq):.1f}i"
                            else:
                                _freq_str = f"{_freq:.1f}"
                            write(
                                f"  {_label:<6s}  {_freq_str:>10s}  {_iri:>12.2f}\n"
                            )
                    else:
                        write(_header.format("Mode", "Freq/cm\u207b\u00b9") + "\n")
                        write("  " + "-" * 52 + "\n")
                        for k in range(_n_skip, len(_freqs)):
                            _label = f"{k - _n_skip + 1}"
                            _freq = _freqs[k]
                            if _freq < 0:
                                _freq_str = f"{abs(_freq):.1f}i"
                            else:
                                _freq_str = f"{_freq:.1f}"
                            write(f"  {_label:<6s}  {_freq_str:>10s}\n")
                    write("\n")
                    flush()
                except Exception as _hess_exc:
                    write(
                        f"\n  ## Vibrational Frequencies\n"
                        f"  " + "-" * 52 + "\n"
                        f"  FAILED: {type(_hess_exc).__name__}: {_hess_exc}\n"
                    )
                    flush()
                    hessian_result = None

    _output_writer.record(out_path, wall_time_s=t_total)

    # Fill late-bound execution and finite-torus fields while preserving the
    # upfront ``status = running`` manifest and all dispatch outcomes.
    _manifest_run_fields = {
        "exchange_q0": exchange_q0_label(exchange_exxdiv),
        "jk_method_requested": _jk_requested_label,
        "jk_method_resolved": resolved_jk.value,
        "jk_method_executed": _executed_jk_method,
        "dft_plus_u": bool(dft_plus_u),
        "dft_plus_u_route": _dft_plus_u_route,
    }
    if resolved_jk == PeriodicJKMethod.AICCM2026DEV_B:
        _b_convention = getattr(result, "finite_torus_convention", None)
        _b_exchange_q0 = getattr(_b_convention, "exchange_q0", None)
        _manifest_run_fields["exchange_q0"] = (
            "not-recorded"
            if _b_exchange_q0 is None
            else exchange_q0_label(_b_exchange_q0)
        )
        _manifest_run_fields["ccm_approach"] = getattr(
            _b_convention,
            "ccm_approach",
            "not-recorded",
        )
        _manifest_run_fields["ccm_construction"] = getattr(
            _b_convention,
            "ccm_construction",
            "not-recorded",
        )
        _manifest_run_fields["evaluation_representation"] = getattr(
            _b_convention,
            "evaluation_representation",
            "not-recorded",
        )
        _manifest_run_fields["exchange_q0_applicability"] = getattr(
            _b_convention,
            "exchange_q0_applicability",
            "not-recorded",
        )
    _output_writer.update_run_fields(_manifest_run_fields)
    if cite_manifest_rows:
        # Feed the full assembled provenance into the .system manifest's
        # [citations] section, mirroring the molecular runner. Without
        # this the periodic .system ended with `count = 0` even though
        # the references were assembled and printed to .out / .bibtex.
        _output_writer.set_citations(cite_manifest_rows)

    # --- DOS (total + projected) for QVF embedding -------------------
    # Compute Fock lattice terms from the converged SCF density, then
    # Gaussian-broaden eigenvalues on a dense Monkhorst-Pack mesh.
    # Both total and projected DOS are serialized into the QVF archive
    # so vibe-view can render interactive side-by-side bands+DOS panels.
    _qvf_dos_data: Optional[dict[str, Any]] = None
    _qvf_pdos_data: Optional[dict[str, Any]] = None
    _qvf_coop_data: Optional[dict[str, Any]] = None
    _qvf_cohp_data: Optional[dict[str, Any]] = None
    if output_qvf and result.converged:
        try:
            system = _system_with_valid_unit_cell_multiplicity(system)

            from vibeqc._vibeqc_core import (
                LatticeSumOptions,
                build_jk_2e_real_space,
                build_fock_2e_real_space,
                compute_kinetic_lattice,
                compute_nuclear_lattice,
                compute_nuclear_lattice_ewald,
                compute_overlap_lattice,
                monkhorst_pack,
            )
            from vibeqc.bands import (
                _HARTREE_TO_EV,
                _density_of_states_from_terms,
                _projected_dos_from_terms,
                ao_groups_per_atom_l,
            )

            lat_opts = LatticeSumOptions()
            from vibeqc._vibeqc_core import CoulombMethod

            lat_opts.coulomb_method = CoulombMethod.EWALD_3D
            lat_opts.cutoff_bohr = 18.0

            # Overlap lattice (same as the density-grid path).
            S_real = compute_overlap_lattice(basis, system, lat_opts)

            _is_unrestricted = hasattr(result, "density_alpha") and hasattr(
                result, "density_beta"
            )
            if _is_unrestricted:
                D_real_alpha = _density_lattice_set_for_output(
                    basis,
                    system,
                    _density_proxy_with_k_metadata(result, result.density_alpha),
                    lat_opts,
                )
                D_real_beta = _density_lattice_set_for_output(
                    basis,
                    system,
                    _density_proxy_with_k_metadata(result, result.density_beta),
                    lat_opts,
                )
                D_real = _sum_lattice_density_sets_for_output(
                    basis,
                    system,
                    lat_opts,
                    D_real_alpha,
                    D_real_beta,
                    label="periodic spin density",
                )
            else:
                D_real = _density_lattice_set_for_output(
                    basis,
                    system,
                    result,
                    lat_opts,
                )

            # Real-space Fock terms: T(g), V(g), and 2e J-K(g).
            # For 3D bulk, use Ewald-summed nuclear attraction to avoid
            # the conditionally convergent point-charge sum.
            T_real = compute_kinetic_lattice(basis, system, lat_opts)
            if system.dim == 3:
                from vibeqc._vibeqc_core import (
                    EwaldOptions,
                    GridOptions,
                    build_grid,
                )

                # Reuse the SCF Ewald a if available (BIPOLE result
                # stores it); otherwise compute a reasonable default.
                ewald_alpha = getattr(result, "ewald_alpha_bohr_inv", None)
                if ewald_alpha is None:
                    from vibeqc.bipole_ext_el_pole import (
                        crystal_default_ewald_alpha,
                    )

                    V_cell_au = float(
                        abs(np.linalg.det(np.asarray(system.lattice, dtype=float)))
                    )
                    ewald_alpha = crystal_default_ewald_alpha(V_cell_au)

                # Molecular grid for Ewald V_ne integration.
                mol = system.unit_cell_molecule()
                grid_opts = GridOptions()
                grid_opts.n_radial = 75
                grid_opts.angular = "lebedev"
                grid_opts.lebedev_order = 29
                grid_opts.partition = "becke"
                grid = build_grid(mol, grid_opts)
                ewald_opts = EwaldOptions()
                ewald_opts.alpha = float(ewald_alpha)
                ewald_opts.real_cutoff_bohr = lat_opts.cutoff_bohr
                ewald_opts.tolerance = 1e-8
                V_real = compute_nuclear_lattice_ewald(
                    basis,
                    system,
                    grid,
                    lat_opts,
                    ewald_opts,
                )
            else:
                V_real = compute_nuclear_lattice(basis, system, lat_opts)
            F2e_real = build_fock_2e_real_space(
                basis,
                system,
                lat_opts,
                D_real,
                1.0,
                0.0,
            )
            fock_terms = [T_real, V_real, F2e_real]

            # For unrestricted, build a separate beta Fock channel.
            if _is_unrestricted:
                F_J_real = build_fock_2e_real_space(
                    basis,
                    system,
                    lat_opts,
                    D_real,
                    0.0,
                    0.0,
                )
                if method_upper == "UHF":
                    output_hf_exchange_fraction = 1.0
                elif method_upper == "UKS" and functional:
                    from vibeqc._vibeqc_core import Functional

                    output_hf_exchange_fraction = float(
                        Functional(functional, 2).hf_exchange_fraction
                    )
                else:
                    output_hf_exchange_fraction = 0.0

                if output_hf_exchange_fraction != 0.0:
                    from vibeqc.pbc_bipole_common import _copy_lattice_with_blocks

                    jk_alpha = build_jk_2e_real_space(
                        basis,
                        system,
                        lat_opts,
                        D_real_alpha,
                        0.0,
                    )
                    jk_beta = build_jk_2e_real_space(
                        basis,
                        system,
                        lat_opts,
                        D_real_beta,
                        0.0,
                    )
                    j_blocks = [
                        np.asarray(block, dtype=float).copy()
                        for block in F_J_real.blocks
                    ]
                    alpha_blocks = [
                        j - output_hf_exchange_fraction
                        * np.asarray(k, dtype=float)
                        for j, k in zip(j_blocks, jk_alpha.K.blocks)
                    ]
                    beta_blocks = [
                        j - output_hf_exchange_fraction
                        * np.asarray(k, dtype=float)
                        for j, k in zip(j_blocks, jk_beta.K.blocks)
                    ]
                    F2e_real_alpha = _copy_lattice_with_blocks(
                        basis,
                        system,
                        lat_opts,
                        F_J_real.cells,
                        alpha_blocks,
                    )
                    F2e_real_beta = _copy_lattice_with_blocks(
                        basis,
                        system,
                        lat_opts,
                        F_J_real.cells,
                        beta_blocks,
                    )
                else:
                    F2e_real_alpha = F_J_real
                    F2e_real_beta = F_J_real
                fock_terms_alpha = [T_real, V_real, F2e_real_alpha]
                fock_terms_beta = [T_real, V_real, F2e_real_beta]
                _n_spin = 2
            else:
                fock_terms_alpha = fock_terms
                fock_terms_beta = None
                _n_spin = 1

            # DOS k-mesh -- denser than the SCF mesh for smooth curves.
            _dos_mesh_ints = list(dos_kmesh) if dos_kmesh is not None else [8, 8, 8]
            dos_kmesh_obj = monkhorst_pack(system, _dos_mesh_ints)

            n_elec = system.n_electrons()
            sigma_ev = 0.05  # eV, Gaussian broadening
            sigma_ha = sigma_ev / _HARTREE_TO_EV

            if _is_unrestricted:
                mult = int(system.multiplicity)
                n_alpha = (n_elec + mult - 1) // 2
                n_beta = (n_elec - mult + 1) // 2
            else:
                n_alpha = n_elec // 2
                n_beta = 0

            # --- Total DOS -------------------------------------------------
            dos_result = _density_of_states_from_terms(
                fock_terms_alpha,
                S_real,
                dos_kmesh_obj,
                sigma=sigma_ha,
                n_grid=500,
                n_electrons_per_cell=(
                    2 * n_alpha if _is_unrestricted else n_elec
                ),
            )
            dos_result_beta = None
            if _is_unrestricted and fock_terms_beta is not None:
                dos_result_beta = _density_of_states_from_terms(
                    fock_terms_beta,
                    S_real,
                    dos_kmesh_obj,
                    sigma=sigma_ha,
                    energy_grid=dos_result.energies,
                    n_grid=500,
                    n_electrons_per_cell=2 * n_beta,
                )

            # Energies in eV, shifted to Fermi = 0 eV.
            e_fermi_values = [
                e for e in (
                    dos_result.e_fermi,
                    getattr(dos_result_beta, "e_fermi", None),
                )
                if e is not None
            ]
            e_fermi_ha = max(e_fermi_values) if e_fermi_values else 0.0
            energies_ev = (dos_result.energies - e_fermi_ha) * _HARTREE_TO_EV
            if dos_result_beta is not None:
                dos_arr = np.stack(
                    [
                        np.asarray(dos_result.dos, dtype=np.float64),
                        np.asarray(dos_result_beta.dos, dtype=np.float64),
                    ],
                    axis=0,
                )
            else:
                dos_arr = np.asarray(dos_result.dos, dtype=np.float64)
            # Convert DOS units: states / Hartree / cell -> states / eV / cell
            dos_arr = dos_arr / _HARTREE_TO_EV

            _qvf_dos_data = {
                "energies": energies_ev,
                "dos": dos_arr,
                "smearing": sigma_ev,
                "smearing_type": "gaussian",
                "fermi_energy_ev": float(e_fermi_ha * _HARTREE_TO_EV),
                "n_electrons": float(n_elec),
                "n_spin": _n_spin,
            }

            # --- Projected DOS --------------------------------------------
            groups = ao_groups_per_atom_l(system, basis)
            if groups:
                pdos_result = _projected_dos_from_terms(
                    fock_terms_alpha,
                    S_real,
                    dos_kmesh_obj,
                    groups,
                    sigma_ha,
                    dos_result.energies,
                    n_grid=500,
                    pad=5.0,
                    n_electrons_per_cell=(
                        2 * n_alpha if _is_unrestricted else n_elec
                    ),
                )
                pdos_result_beta = None
                if _is_unrestricted and fock_terms_beta is not None:
                    pdos_result_beta = _projected_dos_from_terms(
                        fock_terms_beta,
                        S_real,
                        dos_kmesh_obj,
                        groups,
                        sigma_ha,
                        dos_result.energies,
                        n_grid=500,
                        pad=5.0,
                        n_electrons_per_cell=2 * n_beta,
                    )
                # Build channel metadata: (atom_index, symbol, l, label)
                channels: list[dict[str, Any]] = []
                projections_list: list[np.ndarray] = []
                projections_beta_list: list[np.ndarray] = []
                for label in pdos_result.group_labels:
                    contrib = np.asarray(
                        pdos_result.contributions[label],
                        dtype=np.float64,
                    )
                    # Convert from states/Hartree -> states/eV
                    contrib = contrib / _HARTREE_TO_EV
                    projections_list.append(contrib)
                    if pdos_result_beta is not None:
                        contrib_beta = np.asarray(
                            pdos_result_beta.contributions[label],
                            dtype=np.float64,
                        )
                        projections_beta_list.append(
                            contrib_beta / _HARTREE_TO_EV
                        )
                    # Parse label like "Mg1-s" -> atom_index, symbol, l
                    # Use the same element-symbol table as bands.py.
                    from vibeqc.basis_crystal import _ELEMENT_SYMBOLS as _SYMS

                    parts = label.rsplit("-", 1)
                    atom_label = parts[0] if len(parts) == 2 else label
                    l_letter = parts[1] if len(parts) == 2 else "?"
                    # Extract atom index from "H1", "Mg2", etc.
                    atom_idx = 0
                    atom_symbol = label
                    for a_idx, atom in enumerate(system.unit_cell):
                        z = int(atom.Z)
                        sym = _SYMS[z] if 0 < z < len(_SYMS) else f"Z{z}"
                        expected = sym + str(a_idx + 1)
                        if expected == atom_label:
                            atom_idx = a_idx
                            atom_symbol = sym
                            break
                    l_map = {"s": 0, "p": 1, "d": 2, "f": 3, "g": 4, "h": 5}
                    l_val = l_map.get(l_letter.lower(), -1)
                    channels.append(
                        {
                            "atom_index": atom_idx,
                            "symbol": atom_symbol,
                            "l": l_val,
                            "label": label,
                        }
                    )

                projections_alpha = np.stack(projections_list, axis=0)
                if projections_beta_list:
                    projections = np.stack(
                        [
                            projections_alpha,
                            np.stack(projections_beta_list, axis=0),
                        ],
                        axis=0,
                    )
                else:
                    projections = projections_alpha

                _qvf_pdos_data = {
                    "energies": energies_ev,
                    "projections": projections,
                    "energies_units": "eV",
                    "n_spin": _n_spin,
                    "fermi_energy_ev": float(e_fermi_ha * _HARTREE_TO_EV),
                    "channels": channels,
                }

            # --- COOP/COHP bonding analysis --------------------------
            if coop_cohp:
                try:
                    from vibeqc.coop_cohp import compute_coop_cohp as _compute_coop_cohp

                    _cc_result = _compute_coop_cohp(
                        fock_terms_alpha,
                        S_real,
                        system,
                        basis,
                        dos_kmesh_obj,
                        H_terms=[T_real, V_real],  # Hcore = T + V
                        F_terms_beta=fock_terms_beta,
                        sigma=sigma_ha,
                        n_grid=500,
                        n_electrons_per_cell=n_elec,
                    )
                    # Energies in eV, shifted to Fermi = 0.
                    _cc_fermi_ha = _cc_result.fermi_energy
                    _cc_energies_ev = (
                        _cc_result.energies - _cc_fermi_ha
                    ) * _HARTREE_TO_EV

                    _qvf_coop_data = {
                        "energies": _cc_energies_ev,
                        "projections": np.asarray(_cc_result.coop, dtype=np.float64),
                        "integrated": np.asarray(
                            _cc_result.integrated_coop, dtype=np.float64
                        ),
                        "energies_units": "eV",
                        "n_spin": _n_spin,
                        "fermi_energy_ev": float(_cc_fermi_ha * _HARTREE_TO_EV),
                        "sigma_ev": float(sigma_ha * _HARTREE_TO_EV),
                        "pairs": _cc_result.pairs,
                    }

                    if _cc_result.cohp is not None:
                        _qvf_cohp_data = {
                            "energies": _cc_energies_ev,
                            "projections": np.asarray(
                                _cc_result.cohp, dtype=np.float64
                            ),
                            "integrated": np.asarray(
                                _cc_result.integrated_cohp, dtype=np.float64
                            ),
                            "energies_units": "eV",
                            "n_spin": _n_spin,
                            "fermi_energy_ev": float(_cc_fermi_ha * _HARTREE_TO_EV),
                            "sigma_ev": float(sigma_ha * _HARTREE_TO_EV),
                            "pairs": _cc_result.pairs,
                        }
                except Exception as _cc_exc:
                    warn_output_failure(
                        _cc_exc,
                        output_stem.with_suffix(".qvf"),
                        role="coop_cohp_computation",
                        category=OutputFailureKind.optional_artifact,
                    )

            # --- Periodic Mayer bond orders (k-space) -----------
            try:
                from vibeqc.basis_crystal import _ELEMENT_SYMBOLS as _SYMS
                from vibeqc.coop_cohp import periodic_mayer_bond_orders as _pmbo

                _bo_matrix = _pmbo(
                    fock_terms_alpha,
                    S_real,
                    system,
                    basis,
                    dos_kmesh_obj,
                    n_electrons_per_cell=n_elec,
                    F_terms_beta=fock_terms_beta,
                )
                _bo_pairs: list[dict[str, Any]] = []
                _atoms = list(system.unit_cell)
                _coords = np.asarray([np.asarray(a.xyz) for a in _atoms])
                for i in range(_bo_matrix.shape[0]):
                    for j in range(i + 1, _bo_matrix.shape[1]):
                        order = float(_bo_matrix[i, j])
                        if order < 0.05:
                            continue
                        zi = int(_atoms[i].Z)
                        zj = int(_atoms[j].Z)
                        sym_i = _SYMS[zi] if 0 < zi < len(_SYMS) else f"Z{zi}"
                        sym_j = _SYMS[zj] if 0 < zj < len(_SYMS) else f"Z{zj}"
                        dist = float(np.linalg.norm(_coords[i] - _coords[j]))
                        _bo_pairs.append(
                            {
                                "i": i,
                                "j": j,
                                "symbol_i": sym_i,
                                "symbol_j": sym_j,
                                "order": order,
                                "distance_ang": float(dist * 0.529177210903),
                            }
                        )
                if _bo_pairs:
                    _qvf_bond_orders = {
                        "method": "mayer",
                        "pairs": _bo_pairs,
                    }
            except Exception as _bo_exc:
                warn_output_failure(
                    _bo_exc,
                    output_stem.with_suffix(".qvf"),
                    role="bond_orders_computation",
                    category=OutputFailureKind.optional_artifact,
                )

        except Exception as _dos_exc:
            warn_output_failure(
                _dos_exc,
                output_stem.with_suffix(".qvf"),
                role="dos_computation",
                category=OutputFailureKind.optional_artifact,
            )

    # --- QVF visualisation archive (v1) ----------------------------------
    if output_qvf:
        try:
            from vibeqc.output.formats.qvf import (
                scf_history_from_result as _scf_history_from_result,
            )

            _qvf_scf_history = _scf_history_from_result(result)
            _output_writer.dispatch_role(
                "qvf",
                system=_qvf_system,
                result=result,
                method=method_upper,
                basis=basis.name,
                functional=functional,
                jk_method=_jk_requested_label,
                jk_method_resolved=resolved_jk.value,
                jk_method_executed=_executed_jk_method,
                dft_plus_u=bool(dft_plus_u),
                dft_plus_u_route=_dft_plus_u_route,
                wall_seconds=t_total,
                volume_data=_qvf_volume_data,
                mo_data=_qvf_mo_data,
                wf_data=_qvf_wf,
                bloch_wf_data=_qvf_bloch_wf,
                hessian_result=hessian_result,
                band_structure=band_structure,
                population_summary=_qvf_pop,
                bond_orders_data=_qvf_bond_orders,
                dipole_moment_data=_qvf_dipole,
                dos_data=_qvf_dos_data,
                pdos_data=_qvf_pdos_data,
                coop_data=_qvf_coop_data,
                cohp_data=_qvf_cohp_data,
                scf_history_data=_qvf_scf_history,
                extensions=_qvf_extensions,
                vendor_json_sections=_qvf_vendor_json_sections,
                raise_on_error=True,
            )
        except Exception as _qvf_exc:
            warn_writer_failure(
                _qvf_exc,
                output_stem.with_suffix(".qvf"),
                role="qvf_archive",
                category=OutputFailureKind.optional_artifact,
                writer=_output_writer,
            )

    # ---- Geometry optimization ---------------------------------------
    opt_result = None
    if optimize and result.converged:
        # The optimize+U guard from prior commits is lifted: all four
        # BIPOLE drivers now accept dft_plus_u= and the per-k Pulay
        # overlap-derivative term is plumbed through
        # bipole_optimize._compute_gradient.
        from .bipole_optimize import relax_atoms, relax_full

        plog.info("")
        plog.banner("Geometry optimization")
        basis_name = basis.name
        # Reconstruct kmesh from the options used during SCF.
        km = _runner_bloch_kmesh(system, kpoints)

        # One append-channel spans the whole relaxation: the per-step
        # relaxation trace (relax_*'s summary lines, routed via progress=)
        # and the final "Optimized geometry" block both land in the .out,
        # matching the molecular path (which passes progress=_out_channel to
        # run_geomopt). Without this the relaxation trace only reached stdout.
        with OutputChannel.to_file(out_path, mode="a") as _opt_channel:
            write(section_header("Geometry optimization", width=56))
            flush()
            if optimize_cell and system.dim == 3:
                from .bipole_optimize import relax_cell_gradient

                # First relax atoms, then cell (gradient-based), then atoms again
                opt_result = relax_atoms(
                    system,
                    basis_name,
                    km,
                    method.upper(),
                    functional=functional,
                    max_iter=optimize_max_iter,
                    conv_tol_grad=optimize_conv_tol_grad,
                    scf_options=opts,
                    cutoff_bohr=float(getattr(opts.lattice_opts, "cutoff_bohr", 8.0)),
                    ewald_precision=ewald_precision,
                    use_multipole_far_field=use_multipole_far_field,
                    multipole_l_max=multipole_l_max,
                    dft_plus_u=dft_plus_u,
                    progress=_opt_channel,
                )
                opt_result = relax_cell_gradient(
                    opt_result.system,
                    basis_name,
                    km,
                    method.upper(),
                    functional=functional,
                    max_iter=10,
                    scf_options=opts,
                    cutoff_bohr=float(getattr(opts.lattice_opts, "cutoff_bohr", 8.0)),
                    ewald_precision=ewald_precision,
                    use_multipole_far_field=use_multipole_far_field,
                    multipole_l_max=multipole_l_max,
                    dft_plus_u=dft_plus_u,
                    progress=_opt_channel,
                )
                opt_result = relax_atoms(
                    opt_result.system,
                    basis_name,
                    km,
                    method.upper(),
                    functional=functional,
                    max_iter=optimize_max_iter,
                    conv_tol_grad=optimize_conv_tol_grad,
                    scf_options=opts,
                    cutoff_bohr=float(getattr(opts.lattice_opts, "cutoff_bohr", 8.0)),
                    ewald_precision=ewald_precision,
                    use_multipole_far_field=use_multipole_far_field,
                    multipole_l_max=multipole_l_max,
                    dft_plus_u=dft_plus_u,
                    progress=_opt_channel,
                )
            else:
                opt_result = relax_atoms(
                    system,
                    basis_name,
                    km,
                    method.upper(),
                    functional=functional,
                    max_iter=optimize_max_iter,
                    conv_tol_grad=optimize_conv_tol_grad,
                    scf_options=opts,
                    cutoff_bohr=float(getattr(opts.lattice_opts, "cutoff_bohr", 8.0)),
                    ewald_precision=ewald_precision,
                    use_multipole_far_field=use_multipole_far_field,
                    multipole_l_max=multipole_l_max,
                    dft_plus_u=dft_plus_u,
                    progress=_opt_channel,
                )
            if opt_result is not None:
                opt_sys = opt_result.system
                write("\n  Optimized geometry\n")
                write("  " + "-" * 56 + "\n")
                write(f"    E_final = {render_energy_labeled(opt_result.energy, width=0, precision=10, sign=True)}\n")
                write(f"    n_iter  = {opt_result.n_iter}\n")
                write(f"    converged = {opt_result.converged}\n")
                write("\n")
                write(_system_summary(opt_sys))
                flush()
                # Write optimized geometry files as SEPARATE artefacts
                # ({stem}.opt.xyz / {stem}.opt.POSCAR), not overwriting the
                # SCF geometry. The writers do Path(stem).with_suffix(EXT),
                # which *replaces* the final suffix -- so the stem is passed
                # already ending in ".opt.<EXT>" to make that replace a no-op
                # and yield the intended compound name. Each is recorded in
                # the manifest (runtime-conditional, so not in real_plan;
                # mark_written appends it to [[outputs.files]]).
                if write_xyz_file:
                    try:
                        _output_writer.dispatch_role(
                            "geometry",
                            runtime_path=output_stem.with_suffix(".opt.xyz"),
                            runtime_format="extended-xyz",
                            runtime_description="Optimized periodic geometry.",
                            system=opt_sys,
                            energy_ha=opt_result.energy,
                            raise_on_error=True,
                        )
                    except Exception as _opt_xyz_exc:
                        warn_writer_failure(
                            _opt_xyz_exc,
                            output_stem.with_suffix(".opt.xyz"),
                            role="optimized_extended_xyz",
                            category=OutputFailureKind.optional_artifact,
                            writer=_output_writer,
                        )
                if write_poscar_file:
                    try:
                        _output_writer.dispatch_role(
                            "geometry",
                            runtime_path=output_stem.with_suffix(".opt.POSCAR"),
                            runtime_format="poscar",
                            runtime_description="Optimized periodic POSCAR geometry.",
                            system=opt_sys,
                            comment=f"vibe-qc optimized {label}",
                            raise_on_error=True,
                        )
                    except Exception as _opt_poscar_exc:
                        warn_writer_failure(
                            _opt_poscar_exc,
                            output_stem.with_suffix(".opt.POSCAR"),
                            role="optimized_poscar",
                            category=OutputFailureKind.optional_artifact,
                            writer=_output_writer,
                        )
                result = opt_result  # Return optimization result

    # Terminal checkpoint frame: label the checkpoint QVF "converged" so a
    # live viewer knows the job settled and can stop watching. Uses the
    # final (optimized, if any) system + result. The OutputWriter lifecycle
    # wrapper also marks the manifest crashed if terminal work raises.
    if _checkpointer.enabled:
        _final_system = getattr(result, "system", None) or system
        _checkpointer.finalize(
            "converged",
            system=_final_system,
            result=result,
            method=method_upper,
            basis=basis.name,
            functional=functional,
        )

    # Optimization appends to the log after the SCF-side record above.
    # Re-recording updates the existing outcome row and final hash. The public
    # lifecycle wrapper finishes the manifest after runtime ``.opt.*`` files.
    _output_writer.record(out_path)

    return result