"""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 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 (
ManifestUpdater,
OutputChannel,
OutputPlan,
dry_run_manifest,
flush,
is_dry_run_estimate_requested,
is_dry_run_requested,
section_header,
warn,
write,
write_extended_xyz,
)
from .output import (
write_cif as _write_cif,
)
from .output import (
write_population as _write_population,
)
from .output import (
write_poscar as _write_poscar,
)
from .output.formats.scf_log import write_scf_trace
from .output._errors import (
OutputFailureKind,
warn_output_failure,
)
from .output.citations import (
citation_manifest_rows,
format_references_block,
load_default_database,
write_bibtex,
write_references,
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 (
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 (
run_krhf_periodic_gdf,
run_krks_periodic_gdf,
run_kuhf_periodic_gdf,
run_kuks_periodic_gdf,
)
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)))
out = []
out.append(" Periodicity")
out.append(" " + "-" * 56)
out.append(f" dimensionality = {dim}D")
active_axes = ", ".join(f"a{i + 1}" for i in range(dim))
out.append(f" active axes = {active_axes}")
out.append(f" {measure_label:<14s} = {measure:.4f} {measure_unit}")
out.append("")
out.append(" Lattice (bohr)")
out.append(" " + "-" * 56)
# ``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):
out.append(f" a{i + 1} = {L[0, i]:14.8f} {L[1, i]:14.8f} {L[2, i]:14.8f}")
if dim < 3:
supercell_volume = float(abs(np.linalg.det(L)))
out.append(f" embedding volume = {supercell_volume:.4f} bohr^3")
out.append("")
out.append(f" Atoms (bohr) -- {len(system.unit_cell)} in unit cell")
out.append(" " + "-" * 56)
for i, atom in enumerate(system.unit_cell, start=1):
x, y, z = atom.xyz
out.append(f" {i:4d} Z={atom.Z:3d} {x:14.8f} {y:14.8f} {z:14.8f}")
out.append(
f" n_electrons = {system.n_electrons()} "
f"multiplicity = {system.multiplicity}"
)
out.append("")
return "\n".join(out)
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):
if any(np.iscomplexobj(np.asarray(block)) for block in D_attr.blocks):
from ._vibeqc_core import compute_overlap_lattice
D_set = compute_overlap_lattice(basis, system, lat_opts)
if len(D_set.blocks) != len(D_attr.blocks):
raise ValueError(
"density lattice-set output received incompatible cell lists"
)
for i, block in enumerate(D_attr.blocks):
D_set.set_block(
i,
_real_block_for_periodic_output(
block,
label=f"periodic density block g={tuple(D_attr.cells[i].index)}",
),
)
return D_set
return D_attr
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]
for other in density_sets[1:]:
if len(other.blocks) != len(ref.blocks):
raise ValueError("spin density lattice sets have different sizes")
from vibeqc.pbc_bipole_common import _copy_lattice_with_blocks
blocks = []
for i, parts in enumerate(zip(*(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",
"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)
[docs]
def run_periodic_job(
system: PeriodicSystem,
basis: 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,
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.
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."
)
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 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_kmesh_size = (
_bloch_kmesh_size(_runner_bloch_kmesh(system, kpoints))
if kpoints is not None
else 1
)
_requested_true_multik = kpoints is not None and _requested_kmesh_size > 1
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:
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"
)
# 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_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_cif=write_cif_file,
output_qvf=output_qvf,
write_density_xsf=write_density,
write_population=write_population_file,
citations=citations,
crash_dump=False,
job_kind="periodic_scf",
)
dry_run_manifest(
dry_plan,
record_hostname=record_hostname,
estimate_bytes=_estimate_bytes,
)
return None
plog = resolve_progress(progress, verbose=verbose)
# --- Opt-in live QVF checkpointing (vibe-view hot-reload) -----------
# Build the checkpointer + its plan up front (the real end-of-run
# OutputPlan isn't built until after the SCF, but write_qvf uses the
# plan only as a type-gate -- the checkpoint sections come from the
# 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 _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; GDF honours the resolved value
# verbatim, so no floor there.
ks_driver_fock_mixing_floor=(resolved_jk == PeriodicJKMethod.BIPOLE),
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
),
},
)
_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")
_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:
write(f" WARNING: {_insulator_smearing_note}\n")
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 = {float(data.get('energy', 0)):.6f} Ha)"
)
write(f" restart_from = {restart_from}\n")
write(
f" restart_energy = {float(data.get('energy', 0)):20.10f} Ha\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,
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,
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:
# 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")
_resolved_level_shift = convergence_strategy.value("level_shift")
if _resolved_level_shift != 0.0:
opts.level_shift = float(_resolved_level_shift)
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,
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,
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,
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,
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(
" 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(
" 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} {disp.energy:14.8f} Ha"
f" ({disp.energy * 627.5094740631:+.4f} kcal/mol)\n"
f" {'E_SCF':>10s} {e_scf:14.8f} Ha\n"
f" {'E_total':>10s} {e_total:14.8f} Ha\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 ._vibeqc_core import Molecule as _Mol
from .io.molden import write_molden
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)
write_molden(
molden_path, mol, basis, molden_result, title=output_stem.name
)
write(f" Molecular orbitals written to {molden_path.name}\n")
except Exception as exc:
write(
f" WARNING: molden write failed: {type(exc).__name__}: {exc}\n"
)
warn_output_failure(
exc,
molden_path,
role="molden",
category=OutputFailureKind.optional_artifact,
)
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,
write_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"
)
write_population_summary(output_stem, _qvf_pop)
elif resolved_jk == PeriodicJKMethod.AICCM2026DEV_B:
from vibeqc.output.formats.population import (
compute_aiccm2026dev_b_population_summary,
unsupported_population_summary,
write_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"
)
write_population_summary(output_stem, _qvf_pop)
else:
if isinstance(result.mo_coeffs, (list, tuple)):
pop_result = _gamma_proxy_for_multi_k(result)
_write_population(
output_stem,
pop_result,
basis,
mol_p,
)
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:
write(
f" WARNING: population dump failed: {type(exc).__name__}: {exc}\n"
)
warn_output_failure(
exc,
output_stem.with_suffix(".population.txt"),
role="population_summary",
category=OutputFailureKind.optional_artifact,
)
# --- 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
write_extended_xyz(
output_stem,
system,
energy_ha=energy_ha,
)
write(
f" Final geometry written to {xyz_path.name} "
f"(extended XYZ with lattice)\n"
)
except Exception as exc:
write(
f" WARNING: extended-xyz write failed: "
f"{type(exc).__name__}: {exc}\n"
)
warn_output_failure(
exc,
xyz_path,
role="extended_xyz",
category=OutputFailureKind.optional_artifact,
)
if write_poscar_file:
poscar_path = output_stem.with_suffix(".POSCAR")
try:
_write_poscar(
output_stem,
system,
comment=f"vibe-qc periodic {label} basis={basis.name}",
)
write(f" Structure written to {poscar_path.name} (VASP-5 POSCAR)\n")
except Exception as exc:
write(
f" WARNING: POSCAR write failed: {type(exc).__name__}: {exc}\n"
)
warn_output_failure(
exc,
poscar_path,
role="poscar",
category=OutputFailureKind.optional_artifact,
)
if write_cif_file:
cif_path = output_stem.with_suffix(".cif")
try:
_write_cif(
output_stem,
system,
comment=f"vibe-qc periodic {label} basis={basis.name}",
)
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_output_failure(
exc,
cif_path,
role="cif",
category=OutputFailureKind.optional_artifact,
)
# --- 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 ManifestUpdater is constructed after this block, so the
# rows are carried forward and fed to set_citations there.
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,
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)
write_bibtex(output_stem, _refs)
write_references(output_stem, _refs)
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:
write(
f" WARNING: citation emission failed: "
f"{type(exc).__name__}: {exc}\n"
)
warn_output_failure(
exc,
output_stem.with_suffix(".bibtex"),
role="citations",
category=OutputFailureKind.optional_artifact,
)
if write_xsf_structure_file:
xsf_struct_path = output_stem.with_suffix(".xsf")
try:
from .xsf import write_xsf_structure
# 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",
)
write_xsf_structure(xsf_struct_path, system)
write(
f" Structure written to {xsf_struct_path.name} "
f"(XSF crystal block)\n"
)
except Exception as exc:
write(
f" WARNING: XSF structure write failed: "
f"{type(exc).__name__}: {exc}\n"
)
warn_output_failure(
exc,
xsf_struct_path,
role="xsf_structure",
category=OutputFailureKind.optional_artifact,
)
# --- 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:
write(
" WARNING: χ-CCM-B QVF Wannier-centre overlay "
f"failed: {type(_qvf_wannier_exc).__name__}: "
f"{_qvf_wannier_exc}\n"
)
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
from .xsf import write_xsf_volume
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:
write(
" WARNING: periodic QVF orbital grid "
f"preparation failed: "
f"{type(_qvf_mo_exc).__name__}: "
f"{_qvf_mo_exc}\n"
)
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:
write_xsf_volume(
xsf_path,
density_write_system,
data=rho,
name=f"{output_stem.name}_density",
)
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:
write(
f" WARNING: density grid evaluation failed: "
f"{type(exc).__name__}: {exc}\n"
)
warn_output_failure(
exc,
xsf_path,
role="density_grid",
category=OutputFailureKind.optional_artifact,
)
# --- 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} {t_scf:12.3f}\n")
write(f" {'SCF avg per iter':<28s} {iter_avg:12.3f} ({n_iter} iters)\n")
write(f" {'Job total':<28s} {t_total:12.3f}\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
# --- System manifest with [plan] + [outputs] (Phase O5) -------------
# Replaces the bare write_system_manifest call with the
# OutputPlan / ManifestUpdater pipeline so the periodic .system
# carries the same declarative-plan + status-flag surface as the
# molecular path. vq's submit pre-flight + status polling reads
# the same fields whether the job is molecular or periodic.
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_population=write_population_file,
citations=citations,
crash_dump=False,
output_qvf=output_qvf,
job_kind="periodic_scf",
)
# Record the exchange-q=0 finite-size convention (§0.5 item 1 reproducibility
# gate): the J/Coulomb kernel is forced, but the exchange q=0 seam is an
# additional convention (BvK-ewald vs strict-zero-mode) that shifts the HF gap
# and hence every post-HF denominator — so A/B/KMP2 comparisons must match on it.
from .periodic.exchange_convention import exchange_q0_label
manifest = ManifestUpdater(
real_plan,
record_hostname=record_hostname,
wall_seconds=t_total,
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,
},
)
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.
manifest.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,
)
from vibeqc.output.formats.qvf import write_qvf as _write_qvf
_qvf_scf_history = _scf_history_from_result(result)
_qvf_path = _write_qvf(
output_stem,
real_plan,
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,
)
try:
manifest.mark_written(_qvf_path)
except Exception as _qvf_rec_exc:
warn_output_failure(
_qvf_rec_exc,
_qvf_path,
role="qvf_manifest_record",
category=OutputFailureKind.manifest_recording,
manifest=manifest,
)
except Exception as _qvf_exc:
warn_output_failure(
_qvf_exc,
output_stem.with_suffix(".qvf"),
role="qvf_archive",
category=OutputFailureKind.optional_artifact,
manifest=manifest,
)
# Record actually-written artefacts so [[outputs.files]] reflects
# what landed on disk (best-effort: a writer that failed silently
# leaves its row with written=false, which is the right signal
# for vq + downstream tooling).
for declared in real_plan.files:
if declared.path.is_file():
try:
manifest.mark_written(declared.path)
except Exception as _rec_exc:
# mark_written hashes the file; a brittle write
# shouldn't trip the wrap-up step -- but we still log
# it so a downstream ``vq fetch`` doesn't see a row
# quietly stuck at written=false.
warn_output_failure(
_rec_exc,
declared.path,
role=f"manifest_record_{declared.role}",
category=OutputFailureKind.manifest_recording,
manifest=manifest,
)
manifest.finish(wall_seconds=t_total)
# ---- 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)
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,
)
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,
)
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,
)
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,
)
if opt_result is not None:
opt_sys = opt_result.system
with OutputChannel.to_file(out_path, mode="a"):
write("\n Optimized geometry\n")
write(" " + "-" * 56 + "\n")
write(f" E_final = {opt_result.energy:+.10f} Ha\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
if write_xyz_file:
try:
write_extended_xyz(
output_stem.with_suffix(".opt"),
opt_sys,
energy_ha=opt_result.energy,
)
except Exception as _opt_xyz_exc:
warn_output_failure(
_opt_xyz_exc,
output_stem.with_suffix(".opt.xyz"),
role="optimized_extended_xyz",
category=OutputFailureKind.optional_artifact,
)
if write_poscar_file:
try:
_write_poscar(
output_stem.with_suffix(".opt"),
opt_sys,
comment=f"vibe-qc optimized {label}",
)
except Exception as _opt_poscar_exc:
warn_output_failure(
_opt_poscar_exc,
output_stem.with_suffix(".opt.POSCAR"),
role="optimized_poscar",
category=OutputFailureKind.optional_artifact,
)
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. Failure labeling for the
# periodic path rides with future periodic crash-status plumbing --
# run_periodic_job does not yet flip its own ``.system`` manifest to
# "crashed" either (a pre-existing gap), so we only label success here.
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,
)
return result