"""Periodic PAO, pair-orbit, and PNO primitives for AICCM2026DEV-B.
The functions here are deliberately independent of the molecular DLPNO
driver. They define the finite Born-von Karman (BvK) objects consumed by
that driver and expose invariants that can be tested before a correlation
energy is trusted.
For an AO overlap ``S`` and an ``S``-orthonormal occupied coefficient matrix
``C_o``, the coefficient-space projector is
``Q = I - C_o C_o^H S``.
Columns of ``Q`` selected by a pair domain are redundant PAOs. Canonical
orthogonalization of their Gram matrix removes the null space without
changing their span. For a pair amplitude matrix ``T_ij`` the Hermitian
model pair density is
``D_ij = (T T^H + T^H T) / (1 + delta_ij)``.
Its eigenvectors are the PNOs and its eigenvalues are non-negative pair
occupations. Tightening the occupation threshold therefore gives nested
ranks and the zero-threshold limit retains the complete PAO pair space.
The MP2 energy itself is not variational, so strict monotonicity of the
energy error remains a numerical validation criterion rather than a theorem.
This is the historical bare-amplitude density, not the spin-adapted
density of Riplinger and Neese, J. Chem. Phys. 138, 034106 (2013), Eq. (23),
or Pinski et al., J. Chem. Phys. 143, 034108 (2015), the paragraph following
Eq. (18). Those use the contravariant amplitude in the pair density; see
``vibeqc.dlpno.pno_density`` for the explicit normalization conventions.
The finite-BvK pair-orbit bookkeeping follows Nejad et al.,
J. Chem. Phys. 163, 214107 (2025),
doi:10.1063/5.0290816: simultaneous translations act on both occupied
indices, and optional space-group permutations enlarge those orbits.
"""
from __future__ import annotations
from dataclasses import dataclass, field
from hashlib import sha256
from itertools import product
from math import prod
from typing import Iterable, Sequence
from ._transform_inputs import integer_array, checked_mesh_size
import numpy as np
__all__ = [
"AICCM2026DevBPAOSpace",
"AICCM2026DevBPNOSpace",
"AICCM2026DevBPairOrbit",
"AICCM2026DevBPointPairPlan",
"AICCM2026DevBTranslationPairPlan",
"build_point_pair_plan",
"build_translation_pair_plan",
"build_pair_natural_orbitals",
"enumerate_pair_orbits",
"full_space_noncanonical_mp2_energy",
"projected_atomic_orbitals",
"translation_permutations",
]
@dataclass(frozen=True)
class AICCM2026DevBPAOSpace:
"""An orthonormal PAO basis for one finite-torus domain."""
coefficients: np.ndarray
orbital_energies: np.ndarray
projector: np.ndarray
gram_eigenvalues: np.ndarray
source_ao_indices: np.ndarray
orthonormality_error: float
occupied_leakage_error: float
discarded_rank: int
@property
def rank(self) -> int:
return int(self.coefficients.shape[1])
@dataclass(frozen=True)
class AICCM2026DevBPNOSpace:
"""Pair natural orbitals obtained from one PAO pair space."""
coefficients: np.ndarray
occupations: np.ndarray
retained_occupations: np.ndarray
orbital_energies: np.ndarray | None
retained_rank: int
discarded_occupation: float
hermiticity_error: float
minimum_occupation: float
@dataclass(frozen=True)
class AICCM2026DevBPairOrbit:
"""One orbit of unordered occupied pairs under finite permutations."""
representative: tuple[int, int]
members: tuple[tuple[int, int], ...]
n_cells: int
@property
def multiplicity(self) -> int:
return len(self.members)
@property
def per_cell_weight(self) -> float:
"""Coefficient multiplying a representative pair energy per cell."""
return self.multiplicity / self.n_cells
@dataclass(frozen=True)
class AICCM2026DevBTranslationPairPlan:
"""Compact exact plan for translation-unique occupied-pair work.
Occupieds are labelled ``(band, cell)`` on the finite group
``Z_n0 x Z_n1 x Z_n2``. Every row represents the unordered pair
``{(band_i, 0), (band_j, displacement)}``; simultaneous translation
supplies its placed members. Same-band displacements are canonical
under ``L ~ -L``. A nonzero self-inverse displacement on an even mesh
has a two-element stabilizer and therefore multiplicity ``n_cells / 2``;
all other rows have multiplicity ``n_cells``.
The structure-of-arrays payload deliberately contains no orbit members
and no ``(n_cells, n_occupied)`` permutation matrix. It is therefore the
workload description that a later streamed MPI/GPU executor can partition
without first materialising the quadratic placed-pair space. This class
is a plan only: constructing it does not enable representative execution
or change a correlation energy.
"""
mesh: tuple[int, int, int]
n_bands: int
band_pairs: np.ndarray
displacements: np.ndarray
multiplicities: np.ndarray
schema: str = field(
init=False,
default="vibeqc.aiccm2026dev-b.translation-pair-plan/v1",
)
fingerprint: str = field(init=False)
def __post_init__(self) -> None:
mesh = _validated_pair_mesh(self.mesh)
n_bands = _validated_positive_integer(self.n_bands, "n_bands")
n_cells, n_occupied, n_placed_pairs = _checked_pair_counts(
mesh,
n_bands,
)
raw_band_pairs = np.asarray(self.band_pairs)
raw_displacements = np.asarray(self.displacements)
raw_multiplicities = np.asarray(self.multiplicities)
if any(
array.dtype.kind not in "iu"
for array in (raw_band_pairs, raw_displacements, raw_multiplicities)
):
raise ValueError("translation-pair plan arrays must contain integers")
if raw_band_pairs.ndim != 2 or raw_band_pairs.shape[1:] != (2,):
raise ValueError("translation-pair plan band_pairs must have shape (n, 2)")
n_representatives = int(raw_band_pairs.shape[0])
if raw_displacements.shape != (n_representatives, 3):
raise ValueError(
"translation-pair plan displacements must have shape (n, 3)"
)
if raw_multiplicities.shape != (n_representatives,):
raise ValueError(
"translation-pair plan multiplicities must have shape (n,)"
)
maximum = np.iinfo(np.int64).max
for name, array in (
("band_pairs", raw_band_pairs),
("displacements", raw_displacements),
("multiplicities", raw_multiplicities),
):
if array.size and (
int(np.min(array)) < 0 or int(np.max(array)) > maximum
):
raise ValueError(f"translation-pair plan {name} is out of range")
band_pairs = np.ascontiguousarray(
raw_band_pairs,
dtype=np.int64,
)
displacements = np.ascontiguousarray(
raw_displacements,
dtype=np.int64,
)
multiplicities = np.ascontiguousarray(
raw_multiplicities,
dtype=np.int64,
)
if band_pairs.size and int(np.max(band_pairs)) >= n_bands:
raise ValueError("translation-pair plan band index is out of range")
for axis, extent in enumerate(mesh):
if displacements.size and int(np.max(displacements[:, axis])) >= extent:
raise ValueError(
"translation-pair plan displacement is outside its finite mesh"
)
if multiplicities.size and int(np.min(multiplicities)) <= 0:
raise ValueError("translation-pair plan multiplicities must be positive")
rows = [
(
int(band_pairs[index, 0]),
int(band_pairs[index, 1]),
*(int(value) for value in displacements[index]),
)
for index in range(n_representatives)
]
if rows != sorted(rows) or len(set(rows)) != len(rows):
raise ValueError(
"translation-pair plan representatives must be unique canonical order"
)
for index, row in enumerate(rows):
band_i, band_j, *raw_displacement = row
displacement = tuple(raw_displacement)
inverse = tuple(
(-displacement[axis]) % mesh[axis]
for axis in range(3)
)
if band_i > band_j or (
band_i == band_j and displacement > inverse
):
raise ValueError(
"translation-pair plan representative is not canonical"
)
expected_multiplicity = n_cells
if (
band_i == band_j
and displacement != (0, 0, 0)
and displacement == inverse
):
expected_multiplicity //= 2
if int(multiplicities[index]) != expected_multiplicity:
raise ValueError(
"translation-pair plan stabilizer multiplicity is inconsistent"
)
if sum(int(value) for value in multiplicities) != n_placed_pairs:
raise ValueError("translation-pair plan placed-pair census is incomplete")
expected_representatives = _translation_pair_representative_count(
mesh,
n_bands,
)
if n_representatives != expected_representatives:
raise ValueError(
"translation-pair plan representative census is inconsistent"
)
del n_cells, n_occupied
# A merely read-only owning ndarray can be made writable again with
# ``setflags(write=True)``. Rebuild each canonical payload from a
# Python ``bytes`` buffer so neither a caller alias nor a flag change
# can invalidate the fingerprint after construction.
band_pairs = np.frombuffer(
band_pairs.tobytes(order="C"),
dtype=np.int64,
).reshape(n_representatives, 2)
displacements = np.frombuffer(
displacements.tobytes(order="C"),
dtype=np.int64,
).reshape(n_representatives, 3)
multiplicities = np.frombuffer(
multiplicities.tobytes(order="C"),
dtype=np.int64,
)
object.__setattr__(self, "mesh", mesh)
object.__setattr__(self, "n_bands", n_bands)
object.__setattr__(self, "band_pairs", band_pairs)
object.__setattr__(self, "displacements", displacements)
object.__setattr__(self, "multiplicities", multiplicities)
object.__setattr__(self, "fingerprint", self._build_fingerprint())
def _build_fingerprint(self) -> str:
digest = sha256()
digest.update(self.schema.encode("ascii"))
digest.update(b"\0")
for value in (*self.mesh, self.n_bands):
digest.update(int(value).to_bytes(8, "big", signed=False))
digest.update(np.asarray(self.band_pairs, dtype=">i8").tobytes())
digest.update(np.asarray(self.displacements, dtype=">i8").tobytes())
digest.update(np.asarray(self.multiplicities, dtype=">i8").tobytes())
return digest.hexdigest()
@property
def n_cells(self) -> int:
return prod(self.mesh)
@property
def n_occupied(self) -> int:
return self.n_cells * self.n_bands
@property
def n_representatives(self) -> int:
return int(self.band_pairs.shape[0])
@property
def n_placed_pairs(self) -> int:
return self.n_occupied * (self.n_occupied + 1) // 2
@property
def per_cell_weights(self) -> np.ndarray:
return np.asarray(self.multiplicities, dtype=float) / self.n_cells
@property
def payload_nbytes(self) -> int:
return sum(
int(array.nbytes)
for array in (
self.band_pairs,
self.displacements,
self.multiplicities,
)
)
@dataclass(frozen=True, init=False)
class AICCM2026DevBPointPairPlan:
"""Compact point-symmetry quotient of a translation-pair plan.
``member_indices`` partitions the rows of ``translation_plan`` into
orbits generated by caller-supplied occupied-index support permutations.
The permutations are not retained: ``action_fingerprint`` binds their
exact canonical support-generator set, while ``parent_fingerprint`` binds
the D123 plan on which they acted. ``combined_multiplicities`` therefore
preserves the full placed-pair census without materialising a single
placed pair.
This factory-only immutable object is a workload plan only. It carries
no orbital signs or phases, does not attest an orbital gauge or a
covariant PAO/PNO domain, and is not consumed by an energy path.
"""
translation_plan: AICCM2026DevBTranslationPairPlan
representative_indices: np.ndarray
orbit_offsets: np.ndarray
member_indices: np.ndarray
combined_multiplicities: np.ndarray
n_action_generators: int
action_fingerprint: str
schema: str = field(
init=False,
default="vibeqc.aiccm2026dev-b.point-pair-plan/v1",
)
parent_fingerprint: str = field(init=False)
fingerprint: str = field(init=False)
def __init__(self, *args: object, **kwargs: object) -> None:
del args, kwargs
raise TypeError(
"AICCM2026DevBPointPairPlan is factory-only; "
"use build_point_pair_plan(...)"
)
@classmethod
def _from_payload(
cls,
*,
translation_plan: AICCM2026DevBTranslationPairPlan,
representative_indices: np.ndarray,
orbit_offsets: np.ndarray,
member_indices: np.ndarray,
combined_multiplicities: np.ndarray,
n_action_generators: int,
action_fingerprint: str,
) -> AICCM2026DevBPointPairPlan:
instance = object.__new__(cls)
object.__setattr__(instance, "translation_plan", translation_plan)
object.__setattr__(
instance,
"representative_indices",
representative_indices,
)
object.__setattr__(instance, "orbit_offsets", orbit_offsets)
object.__setattr__(instance, "member_indices", member_indices)
object.__setattr__(
instance,
"combined_multiplicities",
combined_multiplicities,
)
object.__setattr__(
instance,
"n_action_generators",
n_action_generators,
)
object.__setattr__(instance, "action_fingerprint", action_fingerprint)
instance._validate_freeze_and_fingerprint()
return instance
def _validate_freeze_and_fingerprint(self) -> None:
parent = self.translation_plan
if not isinstance(parent, AICCM2026DevBTranslationPairPlan):
raise ValueError(
"point-pair plan requires an AICCM2026DevBTranslationPairPlan"
)
generator_count = _validated_positive_integer(
self.n_action_generators,
"n_action_generators",
)
if generator_count > int(np.iinfo(np.int64).max):
raise ValueError("point-pair plan generator count exceeds int64")
action_fingerprint = _validated_sha256(
self.action_fingerprint,
"point-pair action fingerprint",
)
raw_representatives = np.asarray(self.representative_indices)
raw_offsets = np.asarray(self.orbit_offsets)
raw_members = np.asarray(self.member_indices)
raw_multiplicities = np.asarray(self.combined_multiplicities)
arrays = (
("representative_indices", raw_representatives),
("orbit_offsets", raw_offsets),
("member_indices", raw_members),
("combined_multiplicities", raw_multiplicities),
)
if any(array.dtype.kind not in "iu" for _, array in arrays):
raise ValueError("point-pair plan arrays must contain integers")
if raw_representatives.ndim != 1:
raise ValueError(
"point-pair plan representative_indices must be one-dimensional"
)
n_orbits = int(raw_representatives.size)
if n_orbits < 1:
raise ValueError("point-pair plan must contain at least one orbit")
if raw_offsets.shape != (n_orbits + 1,):
raise ValueError(
"point-pair plan orbit_offsets must have shape (n_orbits + 1,)"
)
if raw_members.shape != (parent.n_representatives,):
raise ValueError(
"point-pair plan member_indices must partition every parent row"
)
if raw_multiplicities.shape != (n_orbits,):
raise ValueError(
"point-pair plan combined_multiplicities must have shape "
"(n_orbits,)"
)
maximum = int(np.iinfo(np.int64).max)
for name, array in arrays:
if array.size and (
int(np.min(array)) < 0 or int(np.max(array)) > maximum
):
raise ValueError(f"point-pair plan {name} is out of range")
representatives = np.ascontiguousarray(raw_representatives, dtype=np.int64)
offsets = np.ascontiguousarray(raw_offsets, dtype=np.int64)
members = np.ascontiguousarray(raw_members, dtype=np.int64)
multiplicities = np.ascontiguousarray(raw_multiplicities, dtype=np.int64)
if int(offsets[0]) != 0 or int(offsets[-1]) != parent.n_representatives:
raise ValueError("point-pair plan orbit offsets do not cover parent rows")
if np.any(offsets[1:] <= offsets[:-1]):
raise ValueError("point-pair plan orbits must be nonempty")
if not np.array_equal(
np.sort(members),
np.arange(parent.n_representatives, dtype=np.int64),
):
raise ValueError(
"point-pair plan members must be an exact partition of parent rows"
)
if np.any(representatives[1:] <= representatives[:-1]):
raise ValueError(
"point-pair plan representatives must be in canonical order"
)
if multiplicities.size and int(np.min(multiplicities)) <= 0:
raise ValueError(
"point-pair plan combined multiplicities must be positive"
)
for orbit_index in range(n_orbits):
start = int(offsets[orbit_index])
stop = int(offsets[orbit_index + 1])
orbit_members = members[start:stop]
if np.any(orbit_members[1:] <= orbit_members[:-1]):
raise ValueError(
"point-pair plan members must be sorted within each orbit"
)
if int(representatives[orbit_index]) != int(orbit_members[0]):
raise ValueError(
"point-pair plan representative must be its minimum member"
)
parent_multiplicities = parent.multiplicities[orbit_members]
if np.any(parent_multiplicities != parent_multiplicities[0]):
raise ValueError(
"point-pair plan orbit changes a translation stabilizer"
)
expected = sum(int(value) for value in parent_multiplicities)
if int(multiplicities[orbit_index]) != expected:
raise ValueError(
"point-pair plan combined multiplicity is inconsistent"
)
if sum(int(value) for value in multiplicities) != parent.n_placed_pairs:
raise ValueError("point-pair plan placed-pair census is incomplete")
# Rebuild the arrays from immutable bytes-backed storage. This
# detaches caller aliases and prevents changing WRITEABLE after the
# plan fingerprint has been established.
representatives = np.frombuffer(
representatives.tobytes(order="C"),
dtype=np.int64,
)
offsets = np.frombuffer(offsets.tobytes(order="C"), dtype=np.int64)
members = np.frombuffer(members.tobytes(order="C"), dtype=np.int64)
multiplicities = np.frombuffer(
multiplicities.tobytes(order="C"),
dtype=np.int64,
)
object.__setattr__(self, "representative_indices", representatives)
object.__setattr__(self, "orbit_offsets", offsets)
object.__setattr__(self, "member_indices", members)
object.__setattr__(self, "combined_multiplicities", multiplicities)
object.__setattr__(self, "n_action_generators", generator_count)
object.__setattr__(self, "action_fingerprint", action_fingerprint)
object.__setattr__(self, "parent_fingerprint", parent.fingerprint)
object.__setattr__(self, "fingerprint", self._build_fingerprint())
def _build_fingerprint(self) -> str:
digest = sha256()
digest.update(self.schema.encode("ascii"))
digest.update(b"\0")
digest.update(bytes.fromhex(self.parent_fingerprint))
digest.update(bytes.fromhex(self.action_fingerprint))
digest.update(int(self.n_action_generators).to_bytes(8, "big"))
for array in (
self.representative_indices,
self.orbit_offsets,
self.member_indices,
self.combined_multiplicities,
):
digest.update(np.asarray(array, dtype=">i8").tobytes())
return digest.hexdigest()
@property
def n_representatives(self) -> int:
return int(self.representative_indices.size)
@property
def n_translation_representatives(self) -> int:
return self.translation_plan.n_representatives
@property
def n_cells(self) -> int:
return self.translation_plan.n_cells
@property
def n_occupied(self) -> int:
return self.translation_plan.n_occupied
@property
def n_placed_pairs(self) -> int:
return self.translation_plan.n_placed_pairs
@property
def per_cell_weights(self) -> np.ndarray:
return np.asarray(self.combined_multiplicities, dtype=float) / self.n_cells
@property
def payload_nbytes(self) -> int:
return sum(
int(array.nbytes)
for array in (
self.representative_indices,
self.orbit_offsets,
self.member_indices,
self.combined_multiplicities,
)
)
def _validated_positive_integer(value: object, label: str) -> int:
if isinstance(value, (bool, np.bool_)) or not isinstance(
value,
(int, np.integer),
):
raise ValueError(f"translation-pair plan {label} must be an integer")
result = int(value)
if result <= 0:
raise ValueError(f"translation-pair plan {label} must be positive")
return result
def _validated_pair_mesh(mesh: Sequence[int]) -> tuple[int, int, int]:
try:
values = tuple(mesh)
except TypeError as exc:
raise ValueError(
"translation-pair plan mesh must contain three positive integers"
) from exc
if len(values) != 3:
raise ValueError(
"translation-pair plan mesh must contain three positive integers"
)
return tuple(
_validated_positive_integer(value, f"mesh[{axis}]")
for axis, value in enumerate(values)
)
def _checked_pair_counts(
mesh: tuple[int, int, int],
n_bands: int,
) -> tuple[int, int, int]:
maximum = int(np.iinfo(np.int64).max)
n_cells = prod(mesh)
if n_cells > maximum:
raise OverflowError("translation-pair plan cell count exceeds int64")
if n_cells > maximum // n_bands:
raise OverflowError("translation-pair plan occupied count exceeds int64")
n_occupied = n_cells * n_bands
if n_occupied > (2 * maximum) // (n_occupied + 1):
raise OverflowError("translation-pair plan pair count exceeds int64")
n_placed_pairs = n_occupied * (n_occupied + 1) // 2
return n_cells, n_occupied, n_placed_pairs
def _two_torsion_count(mesh: tuple[int, int, int]) -> int:
"""Number of finite translations satisfying ``2 L = 0``."""
return prod(2 if extent % 2 == 0 else 1 for extent in mesh)
def _translation_pair_representative_count(
mesh: tuple[int, int, int],
n_bands: int,
) -> int:
n_cells = prod(mesh)
two_torsion = _two_torsion_count(mesh)
same_band = n_bands * (n_cells + two_torsion) // 2
cross_band = n_bands * (n_bands - 1) // 2 * n_cells
return same_band + cross_band
def build_translation_pair_plan(
n_bands: int,
mesh: Sequence[int],
) -> AICCM2026DevBTranslationPairPlan:
"""Build the compact exact finite-translation pair-orbit plan.
The returned rows are canonical under simultaneous translation and
unordered-pair exchange. Their multiplicities sum to the full placed
pair count ``N_occ (N_occ + 1) / 2`` without allocating that space.
"""
mesh_tuple = _validated_pair_mesh(mesh)
bands = _validated_positive_integer(n_bands, "n_bands")
n_cells, _, _ = _checked_pair_counts(mesh_tuple, bands)
n_representatives = _translation_pair_representative_count(
mesh_tuple,
bands,
)
maximum_index = int(np.iinfo(np.intp).max)
payload_bytes_per_row = 6 * np.dtype(np.int64).itemsize
if n_representatives > maximum_index // payload_bytes_per_row:
raise OverflowError("translation-pair plan payload exceeds addressable memory")
band_pairs = np.empty((n_representatives, 2), dtype=np.int64)
displacements = np.empty((n_representatives, 3), dtype=np.int64)
multiplicities = np.empty(n_representatives, dtype=np.int64)
cursor = 0
for band_i in range(bands):
for band_j in range(band_i, bands):
for displacement in product(*(range(value) for value in mesh_tuple)):
inverse = tuple(
(-displacement[axis]) % mesh_tuple[axis]
for axis in range(3)
)
if band_i == band_j and displacement > inverse:
continue
multiplicity = n_cells
if (
band_i == band_j
and displacement != (0, 0, 0)
and displacement == inverse
):
multiplicity //= 2
band_pairs[cursor] = (band_i, band_j)
displacements[cursor] = displacement
multiplicities[cursor] = multiplicity
cursor += 1
if cursor != n_representatives:
raise RuntimeError("translation-pair plan construction broke its census")
return AICCM2026DevBTranslationPairPlan(
mesh=mesh_tuple,
n_bands=bands,
band_pairs=band_pairs,
displacements=displacements,
multiplicities=multiplicities,
)
def _validated_sha256(value: object, label: str) -> str:
if not isinstance(value, str) or len(value) != 64:
raise ValueError(f"{label} must be a lowercase SHA-256 digest")
if any(character not in "0123456789abcdef" for character in value):
raise ValueError(f"{label} must be a lowercase SHA-256 digest")
return value
def _validated_point_action_generators(
permutations: Iterable[Sequence[int] | np.ndarray],
*,
n_occupied: int,
) -> tuple[np.ndarray, ...]:
"""Return unique occupied permutations in deterministic byte order."""
try:
supplied = tuple(permutations)
except TypeError as exc:
raise ValueError(
"point-pair plan occupied_permutations must be an iterable"
) from exc
if not supplied:
raise ValueError(
"point-pair plan requires at least one occupied permutation"
)
expected = np.arange(n_occupied, dtype=np.int64)
canonical: dict[bytes, np.ndarray] = {}
for supplied_permutation in supplied:
raw = np.asarray(supplied_permutation)
if raw.dtype.kind not in "iu":
raise ValueError("point-pair occupied mappings must contain integers")
if raw.shape != (n_occupied,):
raise ValueError("point-pair occupied permutation shape mismatch")
maximum = int(np.iinfo(np.int64).max)
if raw.size and (
int(np.min(raw)) < 0 or int(np.max(raw)) > maximum
):
raise ValueError("point-pair occupied mapping is out of range")
permutation = np.ascontiguousarray(raw, dtype=np.int64)
if not np.array_equal(np.sort(permutation), expected):
raise ValueError("point-pair occupied mapping is not a permutation")
key = permutation.astype(">i8", copy=False).tobytes()
canonical[key] = permutation
return tuple(canonical[key] for key in sorted(canonical))
def _cell_coordinates(
cell_index: int,
mesh: tuple[int, int, int],
) -> tuple[int, int, int]:
plane = mesh[1] * mesh[2]
coordinate_0, remainder = divmod(cell_index, plane)
coordinate_1, coordinate_2 = divmod(remainder, mesh[2])
return coordinate_0, coordinate_1, coordinate_2
def _cell_index(
coordinates: tuple[int, int, int],
mesh: tuple[int, int, int],
) -> int:
return (
(coordinates[0] * mesh[1] + coordinates[1]) * mesh[2]
+ coordinates[2]
)
def _occupied_translation_generator(
plan: AICCM2026DevBTranslationPairPlan,
axis: int,
) -> np.ndarray:
"""Build one axis generator, rather than the full translation table."""
permutation = np.empty(plan.n_occupied, dtype=np.int64)
for cell_index in range(plan.n_cells):
cell = _cell_coordinates(cell_index, plan.mesh)
target = tuple(
(cell[component] + (1 if component == axis else 0))
% plan.mesh[component]
for component in range(3)
)
target_start = _cell_index(target, plan.mesh) * plan.n_bands
source_start = cell_index * plan.n_bands
permutation[source_start : source_start + plan.n_bands] = (
target_start + np.arange(plan.n_bands, dtype=np.int64)
)
return permutation
def _validate_point_action_normalizer(
plan: AICCM2026DevBTranslationPairPlan,
permutation: np.ndarray,
) -> None:
"""Require conjugation of every axis translation to be a translation."""
inverse = np.empty(plan.n_occupied, dtype=np.int64)
inverse[permutation] = np.arange(plan.n_occupied, dtype=np.int64)
occupied_indices = np.arange(plan.n_occupied, dtype=np.int64)
source_cells = occupied_indices // plan.n_bands
source_bands = occupied_indices % plan.n_bands
source_coordinates = np.asarray(
[_cell_coordinates(int(index), plan.mesh) for index in source_cells],
dtype=np.int64,
)
mesh_array = np.asarray(plan.mesh, dtype=np.int64)
for axis in range(3):
translation = _occupied_translation_generator(plan, axis)
conjugated = permutation[translation[inverse]]
target_cells = conjugated // plan.n_bands
target_bands = conjugated % plan.n_bands
target_coordinates = np.asarray(
[_cell_coordinates(int(index), plan.mesh) for index in target_cells],
dtype=np.int64,
)
displacement = (target_coordinates[0] - source_coordinates[0]) % mesh_array
if not np.array_equal(target_bands, source_bands) or not np.all(
(target_coordinates - source_coordinates) % mesh_array == displacement
):
raise ValueError(
"point-pair occupied action does not normalize the finite "
"translation subgroup"
)
def _canonical_translation_pair_row(
plan: AICCM2026DevBTranslationPairPlan,
occupied_i: int,
occupied_j: int,
) -> tuple[int, int, int, int, int]:
cell_i, band_i = divmod(occupied_i, plan.n_bands)
cell_j, band_j = divmod(occupied_j, plan.n_bands)
coordinates_i = _cell_coordinates(cell_i, plan.mesh)
coordinates_j = _cell_coordinates(cell_j, plan.mesh)
displacement = tuple(
(coordinates_j[axis] - coordinates_i[axis]) % plan.mesh[axis]
for axis in range(3)
)
if band_i > band_j:
band_i, band_j = band_j, band_i
displacement = tuple(
(-displacement[axis]) % plan.mesh[axis]
for axis in range(3)
)
if band_i == band_j:
inverse = tuple(
(-displacement[axis]) % plan.mesh[axis]
for axis in range(3)
)
displacement = min(displacement, inverse)
return band_i, band_j, *displacement
def _induced_translation_row_action(
plan: AICCM2026DevBTranslationPairPlan,
permutation: np.ndarray,
row_lookup: dict[tuple[int, int, int, int, int], int],
) -> np.ndarray:
action = np.empty(plan.n_representatives, dtype=np.int64)
for row_index in range(plan.n_representatives):
band_i, band_j = (int(value) for value in plan.band_pairs[row_index])
displacement = tuple(int(value) for value in plan.displacements[row_index])
occupied_i = band_i
occupied_j = _cell_index(displacement, plan.mesh) * plan.n_bands + band_j
key = _canonical_translation_pair_row(
plan,
int(permutation[occupied_i]),
int(permutation[occupied_j]),
)
try:
action[row_index] = row_lookup[key]
except KeyError as exc:
raise RuntimeError(
"point-pair action did not induce a parent-plan row"
) from exc
if not np.array_equal(
np.sort(action),
np.arange(plan.n_representatives, dtype=np.int64),
):
raise RuntimeError("point-pair action on parent rows is not a permutation")
if np.any(plan.multiplicities[action] != plan.multiplicities):
raise RuntimeError("point-pair action changed a translation stabilizer")
return action
def _point_action_fingerprint(
plan: AICCM2026DevBTranslationPairPlan,
generators: tuple[np.ndarray, ...],
) -> str:
digest = sha256()
digest.update(b"vibeqc.aiccm2026dev-b.occupied-point-action/v1\0")
digest.update(bytes.fromhex(plan.fingerprint))
digest.update(len(generators).to_bytes(8, "big"))
for generator in generators:
digest.update(np.asarray(generator, dtype=">i8").tobytes())
return digest.hexdigest()
def build_point_pair_plan(
translation_plan: AICCM2026DevBTranslationPairPlan,
occupied_permutations: Iterable[Sequence[int] | np.ndarray],
) -> AICCM2026DevBPointPairPlan:
"""Build an exact compact point quotient of D123 pair-plan rows.
Each supplied support permutation maps a source occupied index to its
target and must normalize the finite translation subgroup. The induced
support action is evaluated on D123 row anchors, never on their
placed-pair members. Orbit closure over those exact row permutations
produces a second compact plan whose multiplicities retain the complete
placed-pair census.
The API is deliberately permutation-support-only. It neither accepts nor
fingerprints the signs/phases of a physical monomial representation.
General unitary occupied actions cannot map an individual pair row to
another row and are rejected by construction rather than approximated
geometrically.
"""
if not isinstance(translation_plan, AICCM2026DevBTranslationPairPlan):
raise ValueError(
"build_point_pair_plan requires an "
"AICCM2026DevBTranslationPairPlan"
)
generators = _validated_point_action_generators(
occupied_permutations,
n_occupied=translation_plan.n_occupied,
)
row_lookup = {
(
int(translation_plan.band_pairs[row, 0]),
int(translation_plan.band_pairs[row, 1]),
*(int(value) for value in translation_plan.displacements[row]),
): row
for row in range(translation_plan.n_representatives)
}
row_actions: list[np.ndarray] = []
for generator in generators:
_validate_point_action_normalizer(translation_plan, generator)
row_actions.append(
_induced_translation_row_action(
translation_plan,
generator,
row_lookup,
)
)
visited = np.zeros(translation_plan.n_representatives, dtype=bool)
orbits: list[tuple[int, ...]] = []
for seed in range(translation_plan.n_representatives):
if visited[seed]:
continue
members = {seed}
queue = [seed]
cursor = 0
while cursor < len(queue):
member = queue[cursor]
cursor += 1
for action in row_actions:
mapped = int(action[member])
if mapped not in members:
members.add(mapped)
queue.append(mapped)
orbit = tuple(sorted(members))
visited[np.asarray(orbit, dtype=np.int64)] = True
orbits.append(orbit)
if not bool(np.all(visited)):
raise RuntimeError("point-pair orbit closure broke its parent-row census")
representative_indices = np.asarray(
[orbit[0] for orbit in orbits],
dtype=np.int64,
)
orbit_offsets = np.empty(len(orbits) + 1, dtype=np.int64)
orbit_offsets[0] = 0
for index, orbit in enumerate(orbits, start=1):
orbit_offsets[index] = orbit_offsets[index - 1] + len(orbit)
member_indices = np.asarray(
[member for orbit in orbits for member in orbit],
dtype=np.int64,
)
combined_multiplicities = np.asarray(
[
sum(int(translation_plan.multiplicities[member]) for member in orbit)
for orbit in orbits
],
dtype=np.int64,
)
return AICCM2026DevBPointPairPlan._from_payload(
translation_plan=translation_plan,
representative_indices=representative_indices,
orbit_offsets=orbit_offsets,
member_indices=member_indices,
combined_multiplicities=combined_multiplicities,
n_action_generators=len(generators),
action_fingerprint=_point_action_fingerprint(
translation_plan,
generators,
),
)
def _validate_occupied_metric(
occupied: np.ndarray,
overlap: np.ndarray,
*,
tolerance: float,
) -> None:
metric = occupied.conj().T @ overlap @ occupied
error = float(np.max(np.abs(metric - np.eye(metric.shape[0]))))
if error > tolerance:
raise ValueError(
"AICCM2026DEV-B PAOs require S-orthonormal occupied orbitals; "
f"maximum error is {error:.3e}"
)
[docs]
def projected_atomic_orbitals(
occupied_coefficients: np.ndarray,
overlap: np.ndarray,
fock: np.ndarray,
*,
ao_indices: Sequence[int] | np.ndarray | None = None,
linear_dependence_threshold: float = 1e-8,
occupied_tolerance: float = 1e-8,
) -> AICCM2026DevBPAOSpace:
"""Project occupieds from selected AOs and canonically orthogonalize.
``ao_indices`` is the union of the atom-centred AOs in a pair domain.
Passing ``None`` selects the full AO space and provides the exact virtual
limit. Eigenvalues are retained relative to the largest PAO Gram
eigenvalue, which makes the threshold invariant to a uniform AO scaling.
"""
occupied = np.asarray(occupied_coefficients)
S = np.asarray(overlap)
F = np.asarray(fock)
if S.ndim != 2 or S.shape[0] != S.shape[1]:
raise ValueError("AICCM2026DEV-B PAO overlap must be square")
if occupied.ndim != 2 or occupied.shape[0] != S.shape[0]:
raise ValueError("AICCM2026DEV-B PAO occupied/overlap shape mismatch")
if F.shape != S.shape:
raise ValueError("AICCM2026DEV-B PAO Fock/overlap shape mismatch")
if linear_dependence_threshold <= 0.0:
raise ValueError("PAO linear-dependence threshold must be positive")
_validate_occupied_metric(
occupied,
S,
tolerance=occupied_tolerance,
)
n_ao = S.shape[0]
if ao_indices is None:
selected = np.arange(n_ao, dtype=int)
else:
selected = np.unique(np.asarray(ao_indices, dtype=int))
if selected.ndim != 1 or selected.size == 0:
raise ValueError("AICCM2026DEV-B PAO domain must select at least one AO")
if selected[0] < 0 or selected[-1] >= n_ao:
raise ValueError("AICCM2026DEV-B PAO domain AO index out of range")
projector = np.eye(n_ao, dtype=np.result_type(occupied, S)) - (
occupied @ (occupied.conj().T @ S)
)
raw = projector[:, selected]
gram = raw.conj().T @ S @ raw
gram = 0.5 * (gram + gram.conj().T)
eigenvalues, eigenvectors = np.linalg.eigh(gram)
order = np.argsort(eigenvalues)[::-1]
eigenvalues = np.asarray(eigenvalues[order], dtype=float)
eigenvectors = eigenvectors[:, order]
scale = max(float(eigenvalues[0]), np.finfo(float).eps)
keep = eigenvalues > linear_dependence_threshold * scale
if not np.any(keep):
raise RuntimeError("AICCM2026DEV-B PAO domain has zero virtual rank")
coefficients = raw @ (
eigenvectors[:, keep] / np.sqrt(eigenvalues[keep])[None, :]
)
pao_metric = coefficients.conj().T @ S @ coefficients
orthogonality_error = float(
np.max(np.abs(pao_metric - np.eye(pao_metric.shape[0])))
)
leakage = float(np.max(np.abs(occupied.conj().T @ S @ coefficients)))
if max(orthogonality_error, leakage) > 10.0 * occupied_tolerance:
raise RuntimeError(
"AICCM2026DEV-B PAO invariant failure: "
f"orthonormality={orthogonality_error:.3e}, "
f"occupied leakage={leakage:.3e}"
)
fock_pao = coefficients.conj().T @ F @ coefficients
fock_pao = 0.5 * (fock_pao + fock_pao.conj().T)
orbital_energies, rotation = np.linalg.eigh(fock_pao)
coefficients = coefficients @ rotation
return AICCM2026DevBPAOSpace(
coefficients=np.asarray(np.real_if_close(coefficients, tol=1000)),
orbital_energies=np.asarray(orbital_energies, dtype=float),
projector=projector,
gram_eigenvalues=eigenvalues,
source_ao_indices=selected,
orthonormality_error=orthogonality_error,
occupied_leakage_error=leakage,
discarded_rank=int(np.count_nonzero(~keep)),
)
[docs]
def build_pair_natural_orbitals(
pair_amplitudes: np.ndarray,
virtual_coefficients: np.ndarray,
*,
diagonal_pair: bool,
occupation_threshold: float,
fock: np.ndarray | None = None,
minimum_rank: int = 1,
negative_occupation_tolerance: float = 1e-10,
) -> AICCM2026DevBPNOSpace:
"""Diagonalize the Hermitian MP2 model pair density.
A zero threshold retains the complete supplied virtual space, including
exactly zero-occupation directions, and is therefore the algebraic
no-truncation limit. Positive thresholds retain occupations strictly
larger than the threshold, subject to ``minimum_rank``.
"""
amplitudes = np.asarray(pair_amplitudes)
virtuals = np.asarray(virtual_coefficients)
if amplitudes.ndim != 2 or amplitudes.shape[0] != amplitudes.shape[1]:
raise ValueError("AICCM2026DEV-B pair amplitudes must be square")
if virtuals.ndim != 2 or virtuals.shape[1] != amplitudes.shape[0]:
raise ValueError("AICCM2026DEV-B PNO virtual/amplitude shape mismatch")
if occupation_threshold < 0.0:
raise ValueError("PNO occupation threshold must be non-negative")
if minimum_rank < 0 or minimum_rank > amplitudes.shape[0]:
raise ValueError("PNO minimum rank is out of range")
adjoint = amplitudes.conj().T
divisor = 2.0 if diagonal_pair else 1.0
pair_density = (amplitudes @ adjoint + adjoint @ amplitudes) / divisor
hermiticity_error = float(
np.max(np.abs(pair_density - pair_density.conj().T))
)
pair_density = 0.5 * (pair_density + pair_density.conj().T)
occupations, eigenvectors = np.linalg.eigh(pair_density)
order = np.argsort(occupations)[::-1]
occupations = np.asarray(occupations[order], dtype=float)
eigenvectors = eigenvectors[:, order]
minimum_occupation = float(occupations[-1])
if minimum_occupation < -negative_occupation_tolerance:
raise RuntimeError(
"AICCM2026DEV-B pair density is not positive semidefinite; "
f"minimum occupation is {minimum_occupation:.3e}"
)
occupations = np.maximum(occupations, 0.0)
if occupation_threshold == 0.0:
keep = np.ones_like(occupations, dtype=bool)
else:
keep = occupations > occupation_threshold
if int(np.count_nonzero(keep)) < minimum_rank:
keep[:minimum_rank] = True
retained_vectors = eigenvectors[:, keep]
coefficients = virtuals @ retained_vectors
orbital_energies: np.ndarray | None = None
if fock is not None:
F = np.asarray(fock)
if F.shape != (virtuals.shape[0], virtuals.shape[0]):
raise ValueError("AICCM2026DEV-B PNO Fock matrix shape mismatch")
fock_pno = coefficients.conj().T @ F @ coefficients
fock_pno = 0.5 * (fock_pno + fock_pno.conj().T)
orbital_energies, rotation = np.linalg.eigh(fock_pno)
coefficients = coefficients @ rotation
return AICCM2026DevBPNOSpace(
coefficients=np.asarray(np.real_if_close(coefficients, tol=1000)),
occupations=occupations,
retained_occupations=occupations[keep],
orbital_energies=(
None
if orbital_energies is None
else np.asarray(orbital_energies, dtype=float)
),
retained_rank=int(np.count_nonzero(keep)),
discarded_occupation=float(np.sum(occupations[~keep])),
hermiticity_error=hermiticity_error,
minimum_occupation=minimum_occupation,
)
def _full_space_mp2_energy_python(
factors_ov: np.ndarray,
eps_occ: np.ndarray,
eps_vir: np.ndarray,
denominator_tolerance: float,
) -> float:
"""Reference tensor expression for the complete-domain MP2 audit."""
integrals = np.einsum(
"Pia,Pjb->ijab",
factors_ov,
factors_ov,
optimize=True,
)
denominator = (
eps_occ[:, None, None, None]
+ eps_occ[None, :, None, None]
- eps_vir[None, None, :, None]
- eps_vir[None, None, None, :]
)
if np.any(denominator >= -denominator_tolerance):
raise RuntimeError(
"AICCM2026DEV-B MP2 audit encountered a non-negative denominator"
)
amplitudes = integrals / denominator
return float(
np.einsum(
"ijab,ijab->",
amplitudes,
2.0 * integrals - integrals.swapaxes(2, 3),
optimize=True,
)
)
def _full_space_mp2_energy_native(
factors: np.ndarray,
occupied_canonical: np.ndarray,
virtual_canonical: np.ndarray,
eps_occ: np.ndarray,
eps_vir: np.ndarray,
denominator_tolerance: float,
) -> float | None:
"""Return native streaming MP2 audit energy, or ``None`` if unavailable."""
try:
from ... import _vibeqc_core
except ImportError:
return None
transform = getattr(_vibeqc_core, "aiccm2026dev_b_3index_mo_transform", None)
kernel = getattr(_vibeqc_core, "aiccm2026dev_b_real_mp2_energy_from_lov", None)
if transform is None or kernel is None:
return None
factors_ov = transform(
np.ascontiguousarray(factors, dtype=float),
np.ascontiguousarray(occupied_canonical, dtype=float),
np.ascontiguousarray(virtual_canonical, dtype=float),
)
try:
return float(
kernel(
np.ascontiguousarray(factors_ov, dtype=float),
np.ascontiguousarray(eps_occ, dtype=float),
np.ascontiguousarray(eps_vir, dtype=float),
float(denominator_tolerance),
)
)
finally:
del factors_ov
def full_space_noncanonical_mp2_energy(
occupied_coefficients: np.ndarray,
virtual_coefficients: np.ndarray,
fock: np.ndarray,
three_center_factors: np.ndarray,
*,
denominator_tolerance: float = 1e-12,
) -> float:
"""Return the closed-shell DF-MP2 energy in rotation-invariant form.
The occupied and virtual inputs may be any real orthonormal gauges. We
diagonalize the two Fock blocks, rotate the fitted three-center factors,
and evaluate the ordinary ordered-pair MP2 expression. This is the
complete-domain, zero-PNO-threshold limit used to audit the iterative
local-pair equations; it is not used for a truncated calculation.
"""
occupied = np.asarray(occupied_coefficients)
virtual = np.asarray(virtual_coefficients)
F = np.asarray(fock)
factors = np.asarray(three_center_factors)
if any(np.iscomplexobj(value) for value in (occupied, virtual, F, factors)):
maximum_imaginary = max(
float(np.max(np.abs(np.asarray(value).imag)))
for value in (occupied, virtual, F, factors)
)
if maximum_imaginary > 1e-10:
raise ValueError(
"AICCM2026DEV-B real-torus MP2 audit requires real tensors"
)
occupied = occupied.real
virtual = virtual.real
F = F.real
factors = factors.real
if occupied.ndim != 2 or virtual.ndim != 2:
raise ValueError("AICCM2026DEV-B MP2 audit coefficients must be matrices")
if occupied.shape[0] != virtual.shape[0] or F.shape != (
occupied.shape[0],
occupied.shape[0],
):
raise ValueError("AICCM2026DEV-B MP2 audit one-particle shape mismatch")
if factors.ndim != 3 or factors.shape[1:] != F.shape:
raise ValueError("AICCM2026DEV-B MP2 audit factor shape mismatch")
f_occ = occupied.T @ F @ occupied
f_vir = virtual.T @ F @ virtual
eps_occ, rotation_occ = np.linalg.eigh(0.5 * (f_occ + f_occ.T))
eps_vir, rotation_vir = np.linalg.eigh(0.5 * (f_vir + f_vir.T))
occupied_canonical = occupied @ rotation_occ
virtual_canonical = virtual @ rotation_vir
native = _full_space_mp2_energy_native(
factors,
occupied_canonical,
virtual_canonical,
eps_occ,
eps_vir,
denominator_tolerance,
)
if native is not None:
return native
factors_ov = np.einsum(
"mi,Pmn,na->Pia",
occupied_canonical,
factors,
virtual_canonical,
optimize=True,
)
return _full_space_mp2_energy_python(
factors_ov,
eps_occ,
eps_vir,
denominator_tolerance,
)
def _translation_input_contract(n_bands, mesh):
bands = int(integer_array(n_bands, name="n_bands", shape=()))
shape = tuple(integer_array(mesh, name="mesh", shape=(3,)).tolist())
cells = checked_mesh_size(shape)
if bands <= 0 or cells * cells * bands > np.iinfo(np.intp).max // 8:
raise ValueError("translation permutations exceed addressable array bounds")
return bands, shape, cells
def _pair_input_contract(permutations, n_cells):
cells = int(integer_array(n_cells, name="n_cells", shape=()))
if cells <= 0:
raise ValueError("pair orbit n_cells must be positive")
return _validate_permutations(permutations), cells
def translation_permutations(
n_bands: int,
mesh: Sequence[int],
) -> tuple[np.ndarray, ...]:
"""Return occupied-index permutations for every finite translation.
Occupied indices use ``index = cell_index * n_bands + band`` with cells
in lexicographic product order. Each permutation acts simultaneously on
both indices of an occupied pair.
"""
n_bands, mesh_tuple, _ = _translation_input_contract(n_bands, mesh)
native = _translation_permutations_native(n_bands, mesh_tuple)
if native is not None:
return native
return _translation_permutations_python(n_bands, mesh_tuple)
def _translation_permutations_native(
n_bands: int,
mesh: tuple[int, int, int],
) -> tuple[np.ndarray, ...] | None:
"""Return native translation permutations, or ``None`` when unavailable."""
n_bands, mesh, expected = _translation_input_contract(n_bands, mesh)
try:
from ... import _vibeqc_core
except ImportError:
return None
kernel = getattr(_vibeqc_core, "aiccm2026dev_b_translation_permutations", None)
if kernel is None:
return None
matrix = np.asarray(
kernel(
int(n_bands),
np.ascontiguousarray(mesh, dtype=np.int64),
),
dtype=int,
)
if matrix.shape != (expected, expected * int(n_bands)):
raise RuntimeError(
"AICCM2026DEV-B native translation-permutation kernel broke shape"
)
return tuple(np.array(row, dtype=int, copy=True) for row in matrix)
def _translation_permutations_python(
n_bands: int,
mesh_tuple: tuple[int, int, int],
) -> tuple[np.ndarray, ...]:
"""Python oracle for finite-translation occupied-index permutations."""
n_bands, mesh_tuple, _ = _translation_input_contract(n_bands, mesh_tuple)
cells = tuple(product(*(range(value) for value in mesh_tuple)))
cell_index = {cell: index for index, cell in enumerate(cells)}
permutations: list[np.ndarray] = []
for shift in cells:
permutation = np.empty(len(cells) * n_bands, dtype=int)
for source_index, cell in enumerate(cells):
target = tuple(
(cell[axis] + shift[axis]) % mesh_tuple[axis]
for axis in range(3)
)
target_index = cell_index[target]
for band in range(n_bands):
permutation[source_index * n_bands + band] = (
target_index * n_bands + band
)
permutations.append(permutation)
return tuple(permutations)
def _canonical_pair(i: int, j: int) -> tuple[int, int]:
return (i, j) if i <= j else (j, i)
def _validate_permutations(
permutations: Iterable[Sequence[int] | np.ndarray],
) -> tuple[np.ndarray, ...]:
output = tuple(integer_array(value, name="permutations") for value in permutations)
if not output:
raise ValueError("AICCM2026DEV-B pair orbits require at least one permutation")
n_orbitals = int(output[0].size)
if n_orbitals <= 0 or n_orbitals * (n_orbitals + 1) > np.iinfo(np.intp).max // 8:
raise ValueError("pair orbit space must be nonempty and addressable")
expected = np.arange(n_orbitals)
for permutation in output:
if permutation.shape != (n_orbitals,):
raise ValueError("AICCM2026DEV-B pair permutation shape mismatch")
if not np.array_equal(np.sort(permutation), expected):
raise ValueError("AICCM2026DEV-B pair mapping is not a permutation")
return output
def _pair_orbits_from_native(
permutations: tuple[np.ndarray, ...],
*,
n_cells: int,
) -> tuple[AICCM2026DevBPairOrbit, ...] | None:
"""Return native pair orbits, or ``None`` when the extension is older."""
permutations, n_cells = _pair_input_contract(permutations, n_cells)
try:
from ... import _vibeqc_core
except ImportError:
return None
kernel = getattr(_vibeqc_core, "aiccm2026dev_b_pair_orbits", None)
if kernel is None:
return None
permutation_matrix = np.ascontiguousarray(np.vstack(permutations), dtype=np.int64)
representatives, offsets, members = kernel(permutation_matrix, int(n_cells))
reps = np.asarray(representatives, dtype=int)
starts = np.asarray(offsets, dtype=int)
flat_members = np.asarray(members, dtype=int)
if reps.ndim != 2 or reps.shape[1] != 2:
raise RuntimeError("AICCM2026DEV-B native pair-orbit representatives broke shape")
if starts.shape != (reps.shape[0] + 1,):
raise RuntimeError("AICCM2026DEV-B native pair-orbit offsets broke shape")
if flat_members.ndim != 2 or flat_members.shape[1] != 2:
raise RuntimeError("AICCM2026DEV-B native pair-orbit members broke shape")
orbits: list[AICCM2026DevBPairOrbit] = []
for orbit_index, representative in enumerate(reps):
start = int(starts[orbit_index])
stop = int(starts[orbit_index + 1])
if start < 0 or stop < start or stop > flat_members.shape[0]:
raise RuntimeError("AICCM2026DEV-B native pair-orbit offsets are invalid")
orbit_members = tuple(
(int(i), int(j))
for i, j in flat_members[start:stop]
)
if not orbit_members:
raise RuntimeError("AICCM2026DEV-B native pair orbit is empty")
orbits.append(
AICCM2026DevBPairOrbit(
representative=(int(representative[0]), int(representative[1])),
members=orbit_members,
n_cells=n_cells,
)
)
return tuple(orbits)
def _enumerate_pair_orbits_python(
permutations: tuple[np.ndarray, ...],
*,
n_cells: int,
) -> tuple[AICCM2026DevBPairOrbit, ...]:
permutations, n_cells = _pair_input_contract(permutations, n_cells)
n_orbitals = int(permutations[0].size)
all_pairs = {
(i, j)
for i in range(n_orbitals)
for j in range(i, n_orbitals)
}
unassigned = set(all_pairs)
orbits: list[AICCM2026DevBPairOrbit] = []
while unassigned:
seed = min(unassigned)
members = {seed}
frontier = [seed]
while frontier:
i, j = frontier.pop()
for permutation in permutations:
mapped = _canonical_pair(
int(permutation[i]),
int(permutation[j]),
)
if mapped not in members:
members.add(mapped)
frontier.append(mapped)
if not members <= all_pairs:
raise RuntimeError("AICCM2026DEV-B pair orbit escaped pair space")
unassigned -= members
ordered = tuple(sorted(members))
orbits.append(
AICCM2026DevBPairOrbit(
representative=ordered[0],
members=ordered,
n_cells=n_cells,
)
)
flattened = [member for orbit in orbits for member in orbit.members]
if len(flattened) != len(all_pairs) or set(flattened) != all_pairs:
raise RuntimeError("AICCM2026DEV-B pair orbits do not partition pair space")
return tuple(sorted(orbits, key=lambda orbit: orbit.representative))
[docs]
def enumerate_pair_orbits(
permutations: Iterable[Sequence[int] | np.ndarray],
*,
n_cells: int,
) -> tuple[AICCM2026DevBPairOrbit, ...]:
"""Partition all unordered occupied pairs under supplied generators.
The breadth-first closure makes the result correct whether callers pass
the full finite symmetry group or only a generating set. The hard gate
verifies that every unordered pair occurs in exactly one orbit. Production
uses the private native χ-CCM kernel when it is present; the Python closure
remains as the reference path for source checkouts with an older extension.
"""
group_generators, n_cells = _pair_input_contract(permutations, n_cells)
native = _pair_orbits_from_native(group_generators, n_cells=n_cells)
if native is not None:
return native
return _enumerate_pair_orbits_python(group_generators, n_cells=n_cells)