"""One-particle properties derivable from an AICCM2026DEV-B SCF state.
Only gauge-invariant quantities obtainable from the converged one-particle
density and finite-torus orbital spectrum live here. A cell dipole is not a
periodic observable, polarization requires a Berry phase, and correlated
response properties require relaxed correlated density matrices; those are
rejected rather than approximated with a cluster position operator.
"""
from __future__ import annotations
from dataclasses import dataclass
from typing import Sequence
import numpy as np
from ._vibeqc_core import BasisSet, PeriodicSystem
__all__ = [
"AICCM2026DevBSCFProperties",
"aiccm2026dev_b_one_electron_expectation",
"derive_aiccm2026dev_b_scf_properties",
]
@dataclass(frozen=True)
class AICCM2026DevBSCFProperties:
"""Gauge-invariant finite-torus SCF properties per primitive cell."""
electronic_method: str
lattice_extension: tuple[int, int, int]
n_electrons: float
n_alpha: float | None
n_beta: float | None
mulliken_charges: np.ndarray
mulliken_spin_populations: np.ndarray | None
band_gap_hartree: float
alpha_gap_hartree: float | None
beta_gap_hartree: float | None
is_metallic_at_finite_net: bool
density_idempotency_error: float
s_squared: float | None
s_squared_ideal: float | None
unavailable_periodic_observables: tuple[str, ...] = (
"cell dipole (origin/gauge dependent)",
"Berry-phase polarization (not implemented)",
"SCF response properties and analytic derivatives (not implemented)",
"correlated properties without a relaxed correlated 1-RDM",
)
def _blocks(value: object) -> list[np.ndarray]:
array = np.asarray(value)
if array.ndim == 2:
return [array]
return [np.asarray(block) for block in value] # type: ignore[union-attr]
def _vector_blocks(value: object) -> list[np.ndarray]:
array = np.asarray(value)
if array.ndim == 1:
return [array]
return [np.asarray(block) for block in value] # type: ignore[union-attr]
def _ao_atom_indices(basis: BasisSet) -> np.ndarray:
indices: list[int] = []
for shell in basis.shells():
angular = int(shell.l)
n_functions = (
2 * angular + 1
if bool(shell.pure)
else (angular + 1) * (angular + 2) // 2
)
indices.extend([int(shell.atom_index)] * n_functions)
if len(indices) != int(basis.nbasis):
raise RuntimeError("basis shell metadata does not match the AO dimension")
return np.asarray(indices, dtype=int)
def aiccm2026dev_b_one_electron_expectation(
result: object,
operator_k: Sequence[np.ndarray] | np.ndarray,
*,
spin: str = "total",
) -> float:
"""Return ``sum_k w_k Tr[P(k) O(k)]`` per primitive cell.
``spin`` may be ``total``, ``alpha``, ``beta``, or ``difference`` for an
unrestricted reference. Restricted references only accept ``total``.
The physical units are the units of the supplied operator matrix.
"""
if not hasattr(result, "aiccm2026dev_b"):
raise TypeError("result is not an AICCM2026DEV-B SCF result")
operators = _blocks(operator_k)
weights = np.asarray(result.kpoint_weights, dtype=float)
if len(operators) != len(weights):
raise ValueError("operator block count does not match the finite character net")
if hasattr(result, "density_alpha"):
alpha = _blocks(result.density_alpha)
beta = _blocks(result.density_beta)
if spin == "alpha":
densities = alpha
elif spin == "beta":
densities = beta
elif spin == "total":
densities = [a + b for a, b in zip(alpha, beta)]
elif spin == "difference":
densities = [a - b for a, b in zip(alpha, beta)]
else:
raise ValueError("spin must be total, alpha, beta, or difference")
else:
if spin != "total":
raise ValueError("restricted references only define spin='total'")
densities = _blocks(result.density)
value = sum(
weight * np.einsum("mn,nm->", density, operator, optimize=True)
for weight, density, operator in zip(weights, densities, operators)
)
if abs(complex(value).imag) > 1e-9:
raise RuntimeError("one-electron expectation has a non-negligible imaginary part")
return float(complex(value).real)
def _gap(energies: object, n_occ: int) -> float:
blocks = _vector_blocks(energies)
if n_occ <= 0 or n_occ >= blocks[0].size:
return float("inf")
homo = max(float(np.asarray(block)[n_occ - 1]) for block in blocks)
lumo = min(float(np.asarray(block)[n_occ]) for block in blocks)
return lumo - homo
[docs]
def derive_aiccm2026dev_b_scf_properties(
result: object,
system: PeriodicSystem,
basis: BasisSet,
*,
metallic_gap_tolerance_hartree: float = 1e-6,
) -> AICCM2026DevBSCFProperties:
"""Derive all currently supported gauge-invariant SCF properties."""
if not bool(getattr(result, "converged", False)):
raise ValueError("properties require a converged AICCM2026DEV-B SCF")
diagnostics = getattr(result, "aiccm2026dev_b", None)
if diagnostics is None:
raise TypeError("result is not an AICCM2026DEV-B SCF result")
overlap = _blocks(result.overlap)
weights = np.asarray(result.kpoint_weights, dtype=float)
ao_atoms = _ao_atom_indices(basis)
n_atoms = len(system.unit_cell)
def populations(densities: list[np.ndarray]) -> np.ndarray:
population = np.zeros(n_atoms)
for weight, density, metric in zip(weights, densities, overlap):
gross = np.real(np.diag(density @ metric))
population += weight * np.bincount(
ao_atoms,
weights=gross,
minlength=n_atoms,
)
return population
ne = int(system.n_electrons())
if hasattr(result, "density_alpha"):
alpha_density = _blocks(result.density_alpha)
beta_density = _blocks(result.density_beta)
alpha_pop = populations(alpha_density)
beta_pop = populations(beta_density)
total_pop = alpha_pop + beta_pop
spin_pop = alpha_pop - beta_pop
two_s = int(system.multiplicity) - 1
n_alpha_occ = (ne + two_s) // 2
n_beta_occ = (ne - two_s) // 2
alpha_gap = _gap(result.mo_energies_alpha, n_alpha_occ)
beta_gap = _gap(result.mo_energies_beta, n_beta_occ)
gap = min(alpha_gap, beta_gap)
n_alpha = float(alpha_pop.sum())
n_beta = float(beta_pop.sum())
else:
total_pop = populations(_blocks(result.density))
spin_pop = None
gap = _gap(result.mo_energies, ne // 2)
alpha_gap = None
beta_gap = None
n_alpha = None
n_beta = None
charges = np.asarray([float(atom.Z) for atom in system.unit_cell]) - total_pop
return AICCM2026DevBSCFProperties(
electronic_method=str(diagnostics.electronic_method),
lattice_extension=tuple(int(value) for value in diagnostics.mesh),
n_electrons=float(total_pop.sum()),
n_alpha=n_alpha,
n_beta=n_beta,
mulliken_charges=charges,
mulliken_spin_populations=spin_pop,
band_gap_hartree=float(gap),
alpha_gap_hartree=(None if alpha_gap is None else float(alpha_gap)),
beta_gap_hartree=(None if beta_gap is None else float(beta_gap)),
is_metallic_at_finite_net=bool(gap <= metallic_gap_tolerance_hartree),
density_idempotency_error=float(diagnostics.density_idempotency_error),
s_squared=getattr(diagnostics, "s_squared", None),
s_squared_ideal=getattr(diagnostics, "s_squared_ideal", None),
)