"""Space-group analysis for the finite AICCM2026DEV-B torus.
The primitive-cell space group does not automatically act on an arbitrary
Born--von Karman cluster. For a diagonal cluster matrix
``N = diag(N1, N2, N3)``, a point operation ``W`` descends to the finite
translation quotient exactly when ``N^-1 W N`` is an integer matrix. This
module applies that test, builds the resulting atom/cell permutations, and
partitions the exact Gamma-centred reciprocal net into symmetry orbits.
This first increment is intentionally diagnostic. It does not replace the
full k net used by SCF and does not skip libint shell pairs or quartets. AO
group averaging is exposed only at Gamma, where nonsymmorphic Bloch phases
are unity. General-k AO sewing matrices and petite-list integral scattering
must be validated before either optimization can be enabled.
"""
from __future__ import annotations
from dataclasses import dataclass
from enum import Enum
from typing import Sequence
import numpy as np
from ._vibeqc_core import BasisSet, PeriodicSystem
from .symmetry_ao import AtomPermutation, build_ao_permutation_matrix
from .symmetry_lattice import lattice_to_cartesian_rotation
__all__ = [
"AICCM2026DevBSymmetryMode",
"AICCM2026DevBSymmetryDiagnostics",
"AICCM2026DevBSymmetryOperation",
"AICCM2026DevBSymmetryPlan",
"build_aiccm2026dev_b_symmetry_plan",
"gamma_matrix_symmetry_residual",
"shell_pair_orbits",
"shell_quartet_orbits",
"symmetrize_gamma_ao_matrix",
]
class AICCM2026DevBSymmetryMode(str, Enum):
"""Supported B-stream symmetry behavior.
``DIAGNOSTIC`` constructs and verifies the symmetry plan but deliberately
leaves the SCF's full k net and integral build unchanged. ``INTEGRALS``
names the requested future petite-list route and fails closed today.
"""
OFF = "off"
DIAGNOSTIC = "diagnostic"
INTEGRALS = "integrals"
@dataclass(frozen=True)
class AICCM2026DevBSymmetryOperation:
"""One cluster-compatible ``{W|w}`` operation and its atom mapping.
``atom_lattice_shifts[a]`` is the integer vector ``q_a`` defined by
``W f_a + w = f_perm[a] + q_a``. Thus a basis center on cell ``r``
maps to cell ``W r + q_a`` modulo the cyclic mesh.
"""
full_group_index: int
rotation: np.ndarray
translation: np.ndarray
atom_permutation: np.ndarray
atom_lattice_shifts: np.ndarray
max_atom_mapping_residual_bohr: float
has_fractional_translation: bool
[docs]
@dataclass(frozen=True)
class AICCM2026DevBSymmetryPlan:
"""Verified symmetry metadata for one finite cyclic cluster."""
mesh: tuple[int, int, int]
space_group_number: int
international_symbol: str
hall_number: int
point_group: str
wyckoff_letters: tuple[str, ...]
site_symmetry_symbols: tuple[str, ...]
equivalent_atoms: tuple[int, ...]
n_operations_full: int
operations: tuple[AICCM2026DevBSymmetryOperation, ...]
incompatible_operation_indices: tuple[int, ...]
full_kpoints_frac: np.ndarray
full_to_irreducible: np.ndarray
irreducible_representative_indices: np.ndarray
irreducible_weights: np.ndarray
time_reversal: bool
acceleration_applied: bool = False
ao_acceleration_status: str = (
"diagnostic only; general-k AO sewing and libint petite lists disabled"
)
@property
def n_operations_compatible(self) -> int:
return len(self.operations)
@property
def n_kpoints_full(self) -> int:
return int(self.full_kpoints_frac.shape[0])
@property
def n_kpoints_irreducible(self) -> int:
return int(self.irreducible_representative_indices.size)
[docs]
def map_cell(
self,
operation_index: int,
atom_index: int,
cell: Sequence[int],
) -> tuple[int, tuple[int, int, int]]:
"""Map ``(atom, cell)`` through one compatible operation."""
operation = self.operations[operation_index]
r = np.asarray(cell, dtype=np.int64).reshape(3)
image = operation.rotation @ r + operation.atom_lattice_shifts[atom_index]
residue = np.mod(image, np.asarray(self.mesh, dtype=np.int64))
return (
int(operation.atom_permutation[atom_index]),
tuple(int(value) for value in residue),
)
@dataclass(frozen=True)
class AICCM2026DevBSymmetryDiagnostics:
"""Post-SCF witness for the non-mutating diagnostic route."""
plan: AICCM2026DevBSymmetryPlan
gamma_fock_residual: float | None
gamma_density_residual: float | None
n_shell_pairs: int
n_unique_shell_pairs: int
n_shell_quartets: int | None
n_unique_shell_quartets: int | None
energy_change_hartree: float = 0.0
def _dataset_value(dataset: object, name: str):
"""Read a spglib dataset field across old/new Python APIs."""
if hasattr(dataset, name):
return getattr(dataset, name)
return dataset[name] # type: ignore[index]
def _normalise_mesh(
system: PeriodicSystem,
mesh: int | Sequence[int],
) -> tuple[int, int, int]:
dim = int(system.dim)
if isinstance(mesh, (int, np.integer)):
values = [int(mesh)] * dim
else:
values = [int(value) for value in mesh]
if len(values) == dim:
values.extend([1] * (3 - dim))
if len(values) != 3 or any(value < 1 for value in values):
raise ValueError("AICCM2026DEV-B symmetry requires a positive 3-axis mesh")
if any(values[axis] != 1 for axis in range(dim, 3)):
raise ValueError("inactive periodic directions must have mesh size one")
return tuple(values) # type: ignore[return-value]
def _cluster_compatible(rotation: np.ndarray, mesh: tuple[int, int, int]) -> bool:
"""Return whether ``N^-1 W N`` is integral, using integer arithmetic."""
W = np.asarray(rotation, dtype=np.int64)
n = np.asarray(mesh, dtype=np.int64)
# (N^-1 W N)_ij = W_ij N_j / N_i. Avoid floating-point tests.
return all(
int(W[i, j] * n[j]) % int(n[i]) == 0
for i in range(3)
for j in range(3)
)
def _fractional_positions(system: PeriodicSystem) -> tuple[np.ndarray, np.ndarray]:
lattice = np.asarray(system.lattice, dtype=float)
inv_lattice = np.linalg.inv(lattice)
positions = np.asarray(
[np.asarray(atom.xyz, dtype=float) @ inv_lattice for atom in system.unit_cell]
)
species = np.asarray([int(atom.Z) for atom in system.unit_cell], dtype=np.int64)
return positions % 1.0, species
def _system_to_spglib_cell(
system: PeriodicSystem,
) -> tuple[np.ndarray, np.ndarray, list[int]]:
"""Convert the repository's row-lattice convention for spglib."""
positions, species = _fractional_positions(system)
return (
np.asarray(system.lattice, dtype=float),
positions,
[int(value) for value in species],
)
def _map_atoms(
system: PeriodicSystem,
rotation: np.ndarray,
translation: np.ndarray,
*,
symprec: float,
) -> tuple[np.ndarray, np.ndarray, float]:
"""Derive ``W f_a + w = f_b + q_a`` with a Cartesian residual gate."""
lattice = np.asarray(system.lattice, dtype=float)
positions, species = _fractional_positions(system)
n_atoms = len(positions)
permutation = np.full(n_atoms, -1, dtype=np.int64)
shifts = np.zeros((n_atoms, 3), dtype=np.int64)
max_residual = 0.0
for atom_index, (position, atomic_number) in enumerate(zip(positions, species)):
image = rotation @ position + translation
matches: list[tuple[int, np.ndarray, float]] = []
for candidate in np.flatnonzero(species == atomic_number):
delta = image - positions[candidate]
shift = np.rint(delta).astype(np.int64)
residual = float(np.linalg.norm((delta - shift) @ lattice))
if residual <= symprec:
matches.append((int(candidate), shift, residual))
if len(matches) != 1:
raise ValueError(
"AICCM2026DEV-B symmetry atom mapping is not unique for "
f"atom {atom_index}: found {len(matches)} matches at "
f"symprec={symprec:.3e} bohr"
)
candidate, shift, residual = matches[0]
permutation[atom_index] = candidate
shifts[atom_index] = shift
max_residual = max(max_residual, residual)
if len(set(int(value) for value in permutation)) != n_atoms:
raise RuntimeError("space-group operation did not induce an atom permutation")
return permutation, shifts, max_residual
def _gamma_mesh_fractional(mesh: tuple[int, int, int]) -> np.ndarray:
"""Full Gamma-centred character mesh in deterministic C order."""
return np.asarray(
[
(i / mesh[0], j / mesh[1], k / mesh[2])
for i in range(mesh[0])
for j in range(mesh[1])
for k in range(mesh[2])
],
dtype=float,
)
def _k_grid_index(kpoint: np.ndarray, mesh: tuple[int, int, int]) -> int:
scaled = np.asarray(kpoint, dtype=float) * np.asarray(mesh, dtype=float)
rounded = np.rint(scaled).astype(np.int64)
if np.max(np.abs(scaled - rounded)) > 1.0e-8:
raise RuntimeError("compatible symmetry operation left the cyclic k net")
residue = np.mod(rounded, np.asarray(mesh, dtype=np.int64))
return int((residue[0] * mesh[1] + residue[1]) * mesh[2] + residue[2])
def _k_orbits(
mesh: tuple[int, int, int],
operations: Sequence[AICCM2026DevBSymmetryOperation],
*,
time_reversal: bool,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
points = _gamma_mesh_fractional(mesh)
parent = np.arange(len(points), dtype=np.int64)
def find(index: int) -> int:
while parent[index] != index:
parent[index] = parent[parent[index]]
index = int(parent[index])
return index
def union(left: int, right: int) -> None:
left_root, right_root = find(left), find(right)
if left_root != right_root:
parent[max(left_root, right_root)] = min(left_root, right_root)
for source, kpoint in enumerate(points):
for operation in operations:
reciprocal_rotation = np.linalg.inv(operation.rotation).T
target = reciprocal_rotation @ kpoint
union(source, _k_grid_index(target, mesh))
if time_reversal:
union(source, _k_grid_index(-kpoint, mesh))
roots = np.asarray([find(index) for index in range(len(points))], dtype=np.int64)
representatives = np.asarray(
sorted(set(int(root) for root in roots)), dtype=np.int64
)
root_to_ir = {int(root): ir for ir, root in enumerate(representatives)}
mapping = np.asarray([root_to_ir[int(root)] for root in roots], dtype=np.int64)
weights = np.bincount(mapping, minlength=len(representatives)).astype(float)
weights /= len(points)
return points, mapping, representatives, weights
[docs]
def build_aiccm2026dev_b_symmetry_plan(
system: PeriodicSystem,
mesh: int | Sequence[int],
*,
symprec: float = 1.0e-5,
time_reversal: bool = True,
require_full_space_group: bool = False,
) -> AICCM2026DevBSymmetryPlan:
"""Analyze the space group and its exact action on a finite BvK torus.
The current implementation is restricted to 3D. Applying ordinary
3D spglib to a slab or wire cell would treat the artificial vacuum as a
genuine lattice direction and can create nonphysical operations; that
case therefore fails closed pending a layer/rod-group implementation.
"""
if int(system.dim) != 3:
raise NotImplementedError(
"AICCM2026DEV-B space-group diagnostics currently require dim=3; "
"layer and rod groups for 2D/1D are not implemented"
)
if not np.isfinite(symprec) or symprec <= 0.0:
raise ValueError("symprec must be a finite positive distance in bohr")
mesh_tuple = _normalise_mesh(system, mesh)
try:
import spglib
except ImportError as exc: # pragma: no cover - core dependency
raise ImportError("AICCM2026DEV-B symmetry diagnostics require spglib") from exc
dataset = spglib.get_symmetry_dataset(
_system_to_spglib_cell(system), symprec=float(symprec)
)
if dataset is None:
raise RuntimeError("spglib failed to determine the crystal space group")
rotations = np.asarray(_dataset_value(dataset, "rotations"), dtype=np.int64)
translations = np.asarray(_dataset_value(dataset, "translations"), dtype=float)
operations: list[AICCM2026DevBSymmetryOperation] = []
incompatible: list[int] = []
for operation_index, (rotation, translation) in enumerate(
zip(rotations, translations)
):
if not _cluster_compatible(rotation, mesh_tuple):
incompatible.append(operation_index)
continue
permutation, shifts, residual = _map_atoms(
system,
rotation,
translation,
symprec=float(symprec),
)
operations.append(
AICCM2026DevBSymmetryOperation(
full_group_index=operation_index,
rotation=rotation.copy(),
translation=translation.copy(),
atom_permutation=permutation,
atom_lattice_shifts=shifts,
max_atom_mapping_residual_bohr=residual,
has_fractional_translation=bool(
np.max(np.abs(translation - np.rint(translation))) > 1.0e-10
),
)
)
if not operations:
raise RuntimeError("cyclic-cluster-compatible subgroup is empty")
if require_full_space_group and incompatible:
raise ValueError(
"cyclic mesh is incompatible with the full crystal space group: "
f"{len(incompatible)} of {len(rotations)} operations fail "
"N^-1 W N integral; use an isotropic/symmetry-compatible mesh or "
"set require_full_space_group=False to use the exact subgroup"
)
points, mapping, representatives, weights = _k_orbits(
mesh_tuple, operations, time_reversal=bool(time_reversal)
)
return AICCM2026DevBSymmetryPlan(
mesh=mesh_tuple,
space_group_number=int(_dataset_value(dataset, "number")),
international_symbol=str(_dataset_value(dataset, "international")),
hall_number=int(_dataset_value(dataset, "hall_number")),
point_group=str(_dataset_value(dataset, "pointgroup")),
wyckoff_letters=tuple(
str(value) for value in _dataset_value(dataset, "wyckoffs")
),
site_symmetry_symbols=tuple(
str(value) for value in _dataset_value(dataset, "site_symmetry_symbols")
),
equivalent_atoms=tuple(
int(value) for value in _dataset_value(dataset, "equivalent_atoms")
),
n_operations_full=len(rotations),
operations=tuple(operations),
incompatible_operation_indices=tuple(incompatible),
full_kpoints_frac=points,
full_to_irreducible=mapping,
irreducible_representative_indices=representatives,
irreducible_weights=weights,
time_reversal=bool(time_reversal),
)
def _gamma_ao_actions(
system: PeriodicSystem,
basis: BasisSet,
plan: AICCM2026DevBSymmetryPlan,
) -> tuple[np.ndarray, ...]:
"""Build Gamma AO actions; all nonsymmorphic Bloch phases equal one."""
lattice = np.asarray(system.lattice, dtype=float)
actions: list[np.ndarray] = []
for operation in plan.operations:
rotation_cart = lattice_to_cartesian_rotation(
operation.rotation,
lattice.T,
)
orthogonality_error = np.linalg.norm(
rotation_cart.T @ rotation_cart - np.eye(3)
)
if orthogonality_error > 1.0e-8:
raise RuntimeError(
"space-group Cartesian rotation is not orthogonal; check lattice "
f"conventions (residual {orthogonality_error:.3e})"
)
atom_mapping = AtomPermutation(
operation.atom_permutation,
operation.atom_lattice_shifts,
)
actions.append(
build_ao_permutation_matrix(basis, rotation_cart, atom_mapping)
)
return tuple(actions)
def _shell_permutations(
basis: BasisSet,
plan: AICCM2026DevBSymmetryPlan,
) -> tuple[np.ndarray, ...]:
"""Map primitive-cell shell indices under compatible operations."""
shells = list(basis.shells())
by_atom: dict[int, list[int]] = {}
positions: dict[int, int] = {}
for shell_index, shell in enumerate(shells):
atom = int(shell.atom_index)
positions[shell_index] = len(by_atom.setdefault(atom, []))
by_atom[atom].append(shell_index)
permutations: list[np.ndarray] = []
for operation in plan.operations:
permutation = np.empty(len(shells), dtype=int)
for shell_index, shell in enumerate(shells):
destination_atom = int(
operation.atom_permutation[int(shell.atom_index)]
)
destination_shells = by_atom.get(destination_atom, [])
slot = positions[shell_index]
if slot >= len(destination_shells):
raise RuntimeError(
"space-group atom permutation changed shell structure"
)
destination = destination_shells[slot]
if int(shells[destination].l) != int(shell.l):
raise RuntimeError(
"space-group shell permutation changed angular momentum"
)
permutation[shell_index] = destination
if not np.array_equal(np.sort(permutation), np.arange(len(shells))):
raise RuntimeError("space-group shell mapping is not a permutation")
permutations.append(permutation)
return tuple(permutations)
def _partition_discrete_orbits(
objects: set[tuple],
actions,
) -> tuple[tuple[tuple, ...], ...]:
unassigned = set(objects)
orbits: list[tuple[tuple, ...]] = []
while unassigned:
seed = min(unassigned)
members = {seed}
frontier = [seed]
while frontier:
current = frontier.pop()
for action in actions:
image = action(current)
if image not in members:
if image not in objects:
raise RuntimeError("shell symmetry orbit escaped its space")
members.add(image)
frontier.append(image)
unassigned -= members
orbits.append(tuple(sorted(members)))
flattened = [member for orbit in orbits for member in orbit]
if len(flattened) != len(objects) or set(flattened) != objects:
raise RuntimeError("shell symmetry orbits do not form a partition")
return tuple(sorted(orbits, key=lambda orbit: orbit[0]))
def _canonical_shell_pair(left: int, right: int) -> tuple[int, int]:
return (left, right) if left <= right else (right, left)
[docs]
def shell_pair_orbits(
basis: BasisSet,
plan: AICCM2026DevBSymmetryPlan,
) -> tuple[tuple[tuple[int, int], ...], ...]:
"""Partition intrinsic-unique primitive shell pairs by point symmetry."""
n_shells = len(list(basis.shells()))
pairs = {
(left, right)
for left in range(n_shells)
for right in range(left, n_shells)
}
permutations = _shell_permutations(basis, plan)
actions = tuple(
lambda pair, permutation=permutation: _canonical_shell_pair(
int(permutation[pair[0]]),
int(permutation[pair[1]]),
)
for permutation in permutations
)
return _partition_discrete_orbits(pairs, actions) # type: ignore[return-value]
def _canonical_shell_quartet(
first: tuple[int, int],
second: tuple[int, int],
) -> tuple[tuple[int, int], tuple[int, int]]:
left = _canonical_shell_pair(*first)
right = _canonical_shell_pair(*second)
return (left, right) if left <= right else (right, left)
[docs]
def shell_quartet_orbits(
basis: BasisSet,
plan: AICCM2026DevBSymmetryPlan,
) -> tuple[tuple[tuple[tuple[int, int], tuple[int, int]], ...], ...]:
"""Partition intrinsic 8-fold-unique primitive shell quartets."""
n_shells = len(list(basis.shells()))
pairs = tuple(
(left, right)
for left in range(n_shells)
for right in range(left, n_shells)
)
quartets = {
(left, right)
for left_index, left in enumerate(pairs)
for right in pairs[left_index:]
}
permutations = _shell_permutations(basis, plan)
actions = tuple(
lambda quartet, permutation=permutation: _canonical_shell_quartet(
(
int(permutation[quartet[0][0]]),
int(permutation[quartet[0][1]]),
),
(
int(permutation[quartet[1][0]]),
int(permutation[quartet[1][1]]),
),
)
for permutation in permutations
)
return _partition_discrete_orbits( # type: ignore[return-value]
quartets,
actions,
)
def symmetrize_gamma_ao_matrix(
matrix: np.ndarray,
system: PeriodicSystem,
basis: BasisSet,
plan: AICCM2026DevBSymmetryPlan,
) -> np.ndarray:
"""Project a Gamma AO matrix onto the compatible space-group invariant.
This is the Reynolds operator ``|G|^-1 sum_g U_g M U_g^T``. It is
rigorous only at Gamma in this increment. The routine is a diagnostic
utility and is not inserted into DIIS or the B-stream SCF build.
"""
array = np.asarray(matrix)
if array.shape != (basis.nbasis, basis.nbasis):
raise ValueError(
f"Gamma AO matrix must have shape {(basis.nbasis, basis.nbasis)}; "
f"got {array.shape}"
)
actions = _gamma_ao_actions(system, basis, plan)
projected = np.zeros_like(array, dtype=np.result_type(array, np.float64))
for action in actions:
projected += action @ array @ action.T
return projected / len(actions)
def gamma_matrix_symmetry_residual(
matrix: np.ndarray,
system: PeriodicSystem,
basis: BasisSet,
plan: AICCM2026DevBSymmetryPlan,
) -> float:
"""Maximum Frobenius violation ``||U_g M U_g^T - M||`` at Gamma."""
array = np.asarray(matrix)
actions = _gamma_ao_actions(system, basis, plan)
return max(
float(np.linalg.norm(action @ array @ action.T - array))
for action in actions
)