"""Exact one-hop backward lifts through Bärnighausen tables.
The public functions in this module invert one tabulated subgroup descent. Coordinates,
Wyckoff parameters, affine maps, modular solves, and returned shifts are rational. A
Cartesian tolerance is used only when accepting a measured structure that is not an exact
solution of the assembled equations. Cell-metric validation covers monoclinic, orthorhombic,
tetragonal, trigonal, hexagonal, and cubic systems; every tabulated trigonal and hexagonal parent
is in a hexagonal-axes standard setting, so their metric constraint is a=b with alpha=beta=90 and
gamma=120. Normalizer retry applies tabulated cosets to child fractional coordinates, maps
successful results back with the exact inverse, and follows tabulated coset order.
On top of the one-hop lift, :func:`highest_symmetry` / :func:`canonicalize` search upward for the
highest-symmetry description of a crystal and return one deterministic, normalizer-canonical
representative. Every search state is reduced to a normal form that collapses same-group
descriptions of the same crystal: mislabeled special sites are demoted, and the state is quotiented
by the group's continuous- and discrete-Euclidean-normalizer translations and its affine-normalizer
cosets. A triclinic (SG 1 or 2) entry is first Niggli-reduced so the result is independent of the
input basis choice, and the returned cell is put in the standard orientation of its metric. The
result is therefore invariant under origin shift, cell-basis choice (relabeling/shear), and site
order for the same crystal, and agrees with direct entry at the crystal's own space group. For a
P1 / unit-cell start, build the ASU in SG 1 and canonicalize it::
cell = Cell(((4, 0, 0), (0, 4, 0), (0, 0, 4)))
sites = [WyckoffSite("a", FracVector((0, 0, 0)), "Cs"),
WyckoffSite("a", FracVector((Fraction(1, 2),) * 3), "Cl")]
p1 = ASUStructure(cell, 1, sites, [Species(...), Species(...)])
result = canonicalize(p1) # result.spacegroup.it_number == 221
The upward search lifts each state through three fail-only tiers, tried in order and only when the
earlier ones return nothing for a state: (1) the direct tabulated lift; (2) tabulated
affine-normalizer-coset retries; (3) a conventional-cell re-choice. The third tier exists because
a centred-lattice parent can be presented, in the reduced cell the search carries, in an axis choice
that misses the parent's exact metric class even though the lattice admits a conforming cell -- an
F-centred cubic (NaCl from its Niggli primitive) is the motivating case. It searches the candidate
parent lattice for a conventional basis meeting the parent metric exactly, derives the implied child
re-expression, and applies it through two arms: an integer lattice-normalizer re-choice on the
conventional lattice (the F-centred cubic case), and -- when the parent conventional cell is an
intrinsic supercell of the child's conventional lattice, so the conventional-lattice search cannot
produce it (the R-centred trigonal case, Bi-166) -- a re-choice enumerated on the child's PRIMITIVE
lattice (centring vectors included; the det-3 hexagonal basis lives there) and crossed with the
child's discrete normalizer translations, which is what places the atoms on the tabulated
splitting's origin coset. Centred-arm re-expressions are verified to conjugate the child's op set
onto itself exactly and may in general be half-integer on the conventional cell (a supported
generalization; every map accepted so far has been integer). The descent round-trip gate stays
authoritative in both arms. The tier is inert whenever an earlier tier succeeds, so it never runs
on a normally-climbing (e.g. P-lattice) input.
"""
import itertools
import logging
import math
from collections import Counter, defaultdict
from dataclasses import dataclass
from fractions import Fraction
from functools import lru_cache
from typing import Any
from httk.core import FracVector, SurdScalar, SurdVector
from httk.atomistic import data
from httk.atomistic.models.cell.cell import Cell
from httk.atomistic.models.cell.params import CellParams
from httk.atomistic.models.structure.asu import ASUStructure, WyckoffSite
from httk.atomistic.models.structure.comparison import same_crystal
from httk.atomistic.symmetry._lattice import finite_translation_cosets
from httk.atomistic.symmetry._periodicity_guard import require_full_periodicity
from httk.atomistic.symmetry._standardization_common import (
_matrix_column_sum_factor,
_matrix_row_sum_factor,
_scaled_precision,
)
from httk.atomistic.symmetry.affine_operation import AffineOperation
from httk.atomistic.symmetry.recognition import (
_cartesian_distance_squared,
structure_tolerance,
)
from httk.atomistic.symmetry.setting_transform import SettingTransform
from httk.atomistic.symmetry.spacegroup import Spacegroup
from httk.atomistic.symmetry.subgroups import (
SubgroupTransform,
_child_sites,
_standard_input,
minimal_supergroups,
subgroup_closure,
subgroup_representation,
subgroup_transforms,
supergroup_closure,
)
__all__ = [
"COMPATIBLE_CRYSTAL_SYSTEMS",
"LiftResult",
"backward_lift",
"canonicalize",
"highest_symmetry",
"lift_candidates",
"rerepresent",
]
[docs]
COMPATIBLE_CRYSTAL_SYSTEMS: dict[str, frozenset[str]] = {
"triclinic": frozenset({"triclinic"}),
"monoclinic": frozenset({"triclinic", "monoclinic"}),
"orthorhombic": frozenset({"triclinic", "monoclinic", "orthorhombic"}),
"tetragonal": frozenset({"triclinic", "monoclinic", "orthorhombic", "tetragonal"}),
"trigonal": frozenset({"triclinic", "trigonal"}),
"hexagonal": frozenset({"triclinic", "trigonal", "hexagonal"}),
"cubic": frozenset({"triclinic", "monoclinic", "orthorhombic", "tetragonal", "trigonal", "cubic"}),
}
_MAX_SOLVER_BRANCHES = 200_000
# The noisy least-squares fallback only helps genuine recognition noise: the wrap box that carried an
# accepted noisy solution was <= 16 in the reviewed cases, while the full sweep on the same Bi
# 2-atom P1 -> P-1 hop enters boxes up to 1296 before this cap (and 46656 without it, ~90 s that only
# ever yields a rejected approximation). The ceiling sits ~3.2x above that measured 1296. It is a
# heuristic, not a proof: the box grows multiplicatively in the number of rows, so a genuinely noisy
# lift whose wrap box exceeds it would be silently missed -- hence the warning log at the skip. The
# full cap still guards the exact path.
_MAX_NOISY_SWEEP_BRANCHES = 4_096
_MAX_FOURIER_MOTZKIN_INEQUALITIES = 20_000
@dataclass(frozen=True, slots=True)
[docs]
class LiftResult:
"""One exact or tolerance-accepted parent representation.
:param asu: The parent-standard-setting asymmetric unit.
:param spacegroup: The parent space group in standard setting.
:param path: Child-first tabulated parent-to-child subgroup transforms used.
:param shift: The continuous-normalizer origin shift from the final hop, expressed in that
hop's parent standard frame. ``path`` and ``shift`` document the lift route; they do not by
themselves reconstruct ``asu``, since :func:`highest_symmetry` additionally passes each state
through an unrecorded normal form and canonical orientation. ``asu`` is authoritative.
When a default (chirality-normalizing) canonicalization flips an enantiomorphic terminal,
``path`` and ``shift`` describe the route to the pre-flip higher-member terminal, while
``asu``/``spacegroup`` carry the flipped lower-member result -- ``asu`` remains authoritative.
:param residual: The largest wrapped fractional residual accepted.
"""
[docs]
path: tuple[SubgroupTransform, ...]
@dataclass(frozen=True, slots=True)
class _Orbit:
index: int
site: WyckoffSite
position: Any
coordinates: tuple[FracVector, ...]
@dataclass(frozen=True, slots=True)
class _Equation:
matrix: tuple[tuple[Fraction, ...], ...]
constant: tuple[Fraction, ...]
@dataclass(frozen=True, slots=True)
class _Candidate:
parent_letter: str
species: str
pieces: tuple[tuple[int, int], ...]
equations: tuple[_Equation, ...]
covered: frozenset[int]
@dataclass(frozen=True, slots=True)
class _MetricCell:
cell: Cell
cartesian_deviation: float
fractional_deviation: Fraction
def _fractions(value: Any) -> Any:
"""Return nested Fraction values from a vector-like value."""
return value.to_fractions() if hasattr(value, "to_fractions") else value
def _matrix(value: Any) -> tuple[tuple[Fraction, ...], ...]:
return tuple(tuple(Fraction(item) for item in row) for row in _fractions(value))
def _vector(value: Any) -> tuple[Fraction, ...]:
return tuple(Fraction(item) for item in _fractions(value))
def _matvec(matrix: tuple[tuple[Fraction, ...], ...], vector: tuple[Fraction, ...]) -> tuple[Fraction, ...]:
return tuple(sum((a * b for a, b in zip(row, vector)), Fraction(0)) for row in matrix)
def _matmul(
left: tuple[tuple[Fraction, ...], ...], right: tuple[tuple[Fraction, ...], ...]
) -> tuple[tuple[Fraction, ...], ...]:
return tuple(
tuple(sum((left[i][k] * right[k][j] for k in range(len(right))), Fraction(0)) for j in range(len(right[0])))
for i in range(len(left))
)
def _transpose(matrix: tuple[tuple[Fraction, ...], ...]) -> tuple[tuple[Fraction, ...], ...]:
return tuple(tuple(row[column] for row in matrix) for column in range(len(matrix[0])))
def _wrapped(value: Fraction) -> Fraction:
"""Wrap one exact fraction into ``[-1/2, 1/2)``."""
return value - math.floor(value + Fraction(1, 2))
def _wrapped_tuple(value: Any) -> tuple[Fraction, ...]:
return tuple(_wrapped(item) for item in _vector(value))
def _ceil_fraction(value: Fraction) -> int:
return -((-value.numerator) // value.denominator)
def _integer_options(row: tuple[Fraction, ...], constant: Fraction) -> range:
low = constant + sum(min(Fraction(0), value) for value in row)
high = constant + sum(max(Fraction(0), value) for value in row)
# x is in [0, 1); including the upper endpoint only adds a harmless branch at a
# boundary. This is the finite per-row bound; no guessed global wrap range is used.
return range(_ceil_fraction(low), math.floor(high) + 1)
@dataclass(frozen=True, slots=True)
class _Inequality:
coefficients: tuple[Fraction, ...]
bound: Fraction
strict: bool
def _fourier_motzkin(inequalities: tuple[_Inequality, ...]) -> tuple[_Inequality, ...]:
"""Eliminate the first variable from a rational inequality system exactly."""
positive = [item for item in inequalities if item.coefficients[0] > 0]
negative = [item for item in inequalities if item.coefficients[0] < 0]
zero = [
_Inequality(item.coefficients[1:], item.bound, item.strict) for item in inequalities if not item.coefficients[0]
]
result = list(zero)
for upper in positive:
for lower in negative:
upper_coefficient = upper.coefficients[0]
lower_coefficient = lower.coefficients[0]
result.append(
_Inequality(
tuple(
-lower_coefficient * a + upper_coefficient * b
for a, b in zip(upper.coefficients[1:], lower.coefficients[1:], strict=True)
),
-lower_coefficient * upper.bound + upper_coefficient * lower.bound,
upper.strict or lower.strict,
)
)
if len(result) > _MAX_FOURIER_MOTZKIN_INEQUALITIES:
# ponytail: capped FM; replace with a polyhedral package only if table dimensions grow materially.
raise ValueError("Fourier-Motzkin inequality cap exceeded")
return tuple(result)
def _inequality_feasible(inequalities: tuple[_Inequality, ...]) -> bool:
return all((0 < item.bound if item.strict else 0 <= item.bound) for item in inequalities)
def _choose_feasible_free_values(inequalities: tuple[_Inequality, ...], free_count: int) -> tuple[Fraction, ...] | None:
stages = [inequalities]
for _ in range(free_count):
stages.append(_fourier_motzkin(stages[-1]))
if not _inequality_feasible(stages[-1]):
return None
chosen: list[Fraction] = []
for stage in reversed(stages[:-1]):
lower: tuple[Fraction, bool] | None = None
upper: tuple[Fraction, bool] | None = None
for inequality in stage:
coefficient = inequality.coefficients[0]
remainder = sum(value * chosen[index] for index, value in enumerate(inequality.coefficients[1:]))
bound = inequality.bound - remainder
if coefficient > 0:
candidate = (bound / coefficient, inequality.strict)
if upper is None or candidate[0] < upper[0] or (candidate[0] == upper[0] and candidate[1]):
upper = candidate
elif coefficient < 0:
candidate = (bound / coefficient, inequality.strict)
if lower is None or candidate[0] > lower[0] or (candidate[0] == lower[0] and candidate[1]):
lower = candidate
elif (
inequality.strict
and not remainder < inequality.bound
or not inequality.strict
and not remainder <= inequality.bound
):
return None
if lower is None or upper is None:
return None
if lower[0] > upper[0] or (lower[0] == upper[0] and (lower[1] or upper[1])):
return None
chosen.append(lower[0] if lower[0] == upper[0] else (lower[0] + upper[0]) / 2)
chosen.reverse()
return tuple(chosen)
def _linear_solve(matrix: tuple[tuple[Fraction, ...], ...], rhs: tuple[Fraction, ...]) -> tuple[Fraction, ...] | None:
"""Return a deterministic exact solution of a boxed rational linear system.
Free variables are selected at the midpoint of their exact Fourier--Motzkin-feasible interval.
"""
rows = [list(row) + [value] for row, value in zip(matrix, rhs)]
if not rows:
return ()
width = len(matrix[0])
pivot_columns: list[int] = []
pivot_row = 0
for column in range(width):
pivot = next((row for row in range(pivot_row, len(rows)) if rows[row][column]), None)
if pivot is None:
continue
rows[pivot_row], rows[pivot] = rows[pivot], rows[pivot_row]
divisor = rows[pivot_row][column]
rows[pivot_row] = [item / divisor for item in rows[pivot_row]]
for row in range(len(rows)):
if row == pivot_row or not rows[row][column]:
continue
factor = rows[row][column]
rows[row] = [a - factor * b for a, b in zip(rows[row], rows[pivot_row])]
pivot_columns.append(column)
pivot_row += 1
if pivot_row == len(rows):
break
if any(not any(row[:width]) and row[width] for row in rows):
return None
free_columns = tuple(column for column in range(width) if column not in pivot_columns)
if len(free_columns) > 12:
raise ValueError("Fourier-Motzkin free-variable cap exceeded")
inequalities: list[_Inequality] = []
for index in range(len(free_columns)):
coefficients = [Fraction(0)] * len(free_columns)
coefficients[index] = Fraction(-1)
inequalities.append(_Inequality(tuple(coefficients), Fraction(0), False))
coefficients[index] = Fraction(1)
inequalities.append(_Inequality(tuple(coefficients), Fraction(1), True))
for row, column in enumerate(pivot_columns):
pivot_coefficients = tuple(-rows[row][free_column] for free_column in free_columns)
constant = rows[row][width]
inequalities.append(_Inequality(tuple(-value for value in pivot_coefficients), constant, False))
inequalities.append(_Inequality(pivot_coefficients, Fraction(1) - constant, True))
free_values = _choose_feasible_free_values(tuple(inequalities), len(free_columns))
if free_values is None:
return None
solution = [Fraction(0)] * width
for column, value in zip(free_columns, free_values, strict=True):
solution[column] = value
for row, column in enumerate(pivot_columns):
solution[column] = rows[row][width] - sum(
rows[row][free_column] * solution[free_column] for free_column in free_columns
)
if any(value < 0 or value >= 1 for value in solution):
return None
return tuple(solution)
def _least_squares(matrix: tuple[tuple[Fraction, ...], ...], rhs: tuple[Fraction, ...]) -> tuple[Fraction, ...] | None:
transpose = _transpose(matrix)
normal = _matmul(transpose, matrix)
target = _matvec(transpose, rhs)
return _linear_solve(normal, target)
def _rational_null_space(rows: tuple[tuple[Fraction, ...], ...], width: int) -> tuple[tuple[Fraction, ...], ...]:
"""Return an exact basis of ``{x : rows @ x = 0}`` by reduced row echelon elimination."""
matrix = [list(row) for row in rows]
pivot_columns: list[int] = []
pivot_row = 0
for column in range(width):
if pivot_row == len(matrix):
break
pivot = next((row for row in range(pivot_row, len(matrix)) if matrix[row][column]), None)
if pivot is None:
continue
matrix[pivot_row], matrix[pivot] = matrix[pivot], matrix[pivot_row]
divisor = matrix[pivot_row][column]
matrix[pivot_row] = [item / divisor for item in matrix[pivot_row]]
for row in range(len(matrix)):
if row != pivot_row and matrix[row][column]:
factor = matrix[row][column]
matrix[row] = [a - factor * b for a, b in zip(matrix[row], matrix[pivot_row])]
pivot_columns.append(column)
pivot_row += 1
free_columns = [column for column in range(width) if column not in pivot_columns]
basis: list[tuple[Fraction, ...]] = []
for free_column in free_columns:
vector = [Fraction(0)] * width
vector[free_column] = Fraction(1)
for row, column in enumerate(pivot_columns):
vector[column] = -matrix[row][free_column]
basis.append(tuple(vector))
return tuple(basis)
def _rational_inverse(matrix: tuple[tuple[Fraction, ...], ...]) -> tuple[tuple[Fraction, ...], ...]:
"""Return the exact inverse of a nonsingular rational square matrix by Gauss--Jordan."""
size = len(matrix)
augmented = [
list(matrix[row]) + [Fraction(1) if row == column else Fraction(0) for column in range(size)]
for row in range(size)
]
for column in range(size):
pivot = next(row for row in range(column, size) if augmented[row][column])
augmented[column], augmented[pivot] = augmented[pivot], augmented[column]
divisor = augmented[column][column]
augmented[column] = [item / divisor for item in augmented[column]]
for row in range(size):
if row != column and augmented[row][column]:
factor = augmented[row][column]
augmented[row] = [a - factor * b for a, b in zip(augmented[row], augmented[column])]
return tuple(tuple(row[size:]) for row in augmented)
def _integer_consistency_system(
matrix: tuple[tuple[Fraction, ...], ...], constants: tuple[Fraction, ...]
) -> tuple[tuple[tuple[int, ...], ...], tuple[int, ...]]:
"""Return the integer Diophantine system ``A n = d`` that ``n - c in col(M)`` requires.
Each left-null-space vector ``L`` of ``M`` gives one exact constraint ``L n = L c``; clearing
denominators turns the rational rows into integers.
"""
height = len(matrix)
unknowns = len(matrix[0]) if matrix else 0
if unknowns:
left_null = _rational_null_space(_transpose(matrix), height)
else:
# Zero unknowns: ``col(M) = {0}``, so every row must equal its constant and the full
# identity is the constraint set.
left_null = tuple(
tuple(Fraction(1) if index == row else Fraction(0) for index in range(height)) for row in range(height)
)
rows: list[tuple[int, ...]] = []
targets: list[int] = []
for vector in left_null:
rhs = sum((coefficient * constant for coefficient, constant in zip(vector, constants)), Fraction(0))
denominator = math.lcm(*(value.denominator for value in (*vector, rhs)))
rows.append(tuple(int(value * denominator) for value in vector))
targets.append(int(rhs * denominator))
return tuple(rows), tuple(targets)
def _integer_diophantine(
rows: tuple[tuple[int, ...], ...], targets: tuple[int, ...], width: int
) -> tuple[tuple[int, ...], tuple[tuple[int, ...], ...]] | None:
"""Solve ``A n = d`` exactly over the integers by unimodular row reduction of ``[A^T | I]``.
Returns one particular integer solution and an integer basis of ``{v : A v = 0}``, or ``None``
when the system has no integer solution.
"""
height = len(rows)
# Each work row is ``(left, right)`` with ``left`` an integer combination of ``A``'s columns and
# ``right`` the coefficients producing it; unimodular row operations preserve that invariant.
work: list[tuple[list[int], list[int]]] = [
([rows[constraint][column] for constraint in range(height)], [1 if k == column else 0 for k in range(width)])
for column in range(width)
]
pivots: list[int] = []
pivot_row = 0
for column in range(height):
if pivot_row == len(work):
break
while True:
candidates = [row for row in range(pivot_row, len(work)) if work[row][0][column]]
if not candidates:
break
leader = min(candidates, key=lambda row: abs(work[row][0][column]))
work[pivot_row], work[leader] = work[leader], work[pivot_row]
settled = True
for row in range(pivot_row + 1, len(work)):
if work[row][0][column]:
factor = work[row][0][column] // work[pivot_row][0][column]
if factor:
work[row] = (
[a - factor * b for a, b in zip(work[row][0], work[pivot_row][0])],
[a - factor * b for a, b in zip(work[row][1], work[pivot_row][1])],
)
if work[row][0][column]:
settled = False
if settled:
break
if work[pivot_row][0][column]:
if work[pivot_row][0][column] < 0:
work[pivot_row] = ([-value for value in work[pivot_row][0]], [-value for value in work[pivot_row][1]])
pivots.append(pivot_row)
pivot_row += 1
null_basis = tuple(tuple(row[1]) for row in work if not any(row[0]))
particular = [0] * width
remainder = list(targets)
for pivot in pivots:
left = work[pivot][0]
column = next(index for index, value in enumerate(left) if value)
if remainder[column] % left[column]:
return None
factor = remainder[column] // left[column]
if factor:
remainder = [a - factor * b for a, b in zip(remainder, left)]
particular = [a + factor * b for a, b in zip(particular, work[pivot][1])]
if any(remainder):
return None
return tuple(particular), null_basis
def _lattice_box_points(
particular: tuple[int, ...],
null_basis: tuple[tuple[int, ...], ...],
options: tuple[range, ...],
) -> tuple[tuple[int, ...], ...] | None:
"""Return every ``n = particular + sum(t_i * null_basis_i)`` lying in the per-row integer box.
Returns ``None`` when the bounded coset enumeration would exceed the solver cap.
"""
if any(not option for option in options):
return ()
lows = tuple(option[0] for option in options)
highs = tuple(option[-1] for option in options)
if not null_basis:
point = particular
return (point,) if all(low <= value <= high for value, low, high in zip(point, lows, highs)) else ()
dimension = len(null_basis)
columns = tuple(tuple(vector[row] for vector in null_basis) for row in range(len(particular)))
gram = tuple(
tuple(sum(a * b for a, b in zip(null_basis[i], null_basis[j])) for j in range(dimension))
for i in range(dimension)
)
pseudo_inverse = _rational_inverse(tuple(tuple(Fraction(value) for value in row) for row in gram))
# ``left = (K^T K)^-1 K^T`` maps a box point back to exact coset coefficients, so each
# coefficient's finite range follows from the box corners componentwise.
ranges: list[range] = []
for index in range(dimension):
projector = tuple(
sum(pseudo_inverse[index][j] * null_basis[j][row] for j in range(dimension))
for row in range(len(particular))
)
low_sum = Fraction(0)
high_sum = Fraction(0)
for row, coefficient in enumerate(projector):
first = coefficient * (lows[row] - particular[row])
second = coefficient * (highs[row] - particular[row])
low_sum += min(first, second)
high_sum += max(first, second)
ranges.append(range(math.floor(low_sum), math.floor(high_sum) + 1))
if math.prod(len(item) for item in ranges) > _MAX_SOLVER_BRANCHES:
return None
points: list[tuple[int, ...]] = []
for coefficients in itertools.product(*ranges):
point = tuple(
particular[row] + sum(coefficients[index] * columns[row][index] for index in range(dimension))
for row in range(len(particular))
)
if all(low <= value <= high for value, low, high in zip(point, lows, highs)):
points.append(point)
return tuple(sorted(points))
def _exact_modular_solution(
matrix: tuple[tuple[Fraction, ...], ...],
constants: tuple[Fraction, ...],
options: tuple[range, ...],
) -> tuple[tuple[Fraction, ...], Fraction, bool] | None:
"""Return the lexicographically first exact boxed solution, or ``None`` to use the fallback.
Only the integer wraps ``n`` that keep ``n - c`` in ``col(M)`` can yield an exact solution;
they form a lattice coset enumerated inside the finite per-row box, in the same deterministic
order the product enumeration visits, so the first-hit result is preserved.
"""
diophantine_rows, diophantine_targets = _integer_consistency_system(matrix, constants)
solved = _integer_diophantine(diophantine_rows, diophantine_targets, len(options))
if solved is None:
return None
particular, null_basis = solved
points = _lattice_box_points(particular, null_basis, options)
if points is None:
return None
for integers in points:
rhs = tuple(Fraction(integer) - constant for integer, constant in zip(integers, constants))
solution = _linear_solve(matrix, rhs)
if solution is not None:
return solution, Fraction(0), True
return None
def _solve_modular(equations: tuple[_Equation, ...]) -> tuple[tuple[Fraction, ...], Fraction, bool] | None:
"""Solve exact modular rows, then return the best exact rational least square fit."""
if not equations:
return (), Fraction(0), True
matrix = tuple(row for equation in equations for row in equation.matrix)
constants = tuple(value for equation in equations for value in equation.constant)
options = tuple(_integer_options(row, constant) for row, constant in zip(matrix, constants))
exact = _exact_modular_solution(matrix, constants, options)
if exact is not None:
return exact
# Noisy inputs have no exact wrap, so fall back to the capped least-squares sweep over the full
# per-row product. The exact path above already ran unconditionally, so only this sweep is capped.
branches = math.prod(len(item) for item in options)
if branches > _MAX_SOLVER_BRANCHES:
raise ValueError("exact modular lift solver branch cap exceeded")
if branches > _MAX_NOISY_SWEEP_BRANCHES:
# No exact wrap exists (the exact path already ran); a noisy match this deep in the wrap box
# is not a real lift, so fail the candidate cheaply rather than grinding the full sweep. This
# can in principle silently drop a genuinely noisy lift in a large box, so record it on the
# same warning channel the branch-cap skip uses.
logging.getLogger(__name__).warning(
"skipping noisy least-squares sweep over %d wrap boxes (cap %d)",
branches,
_MAX_NOISY_SWEEP_BRANCHES,
extra={"context": "symmetry"},
)
return None
best: tuple[tuple[Fraction, ...], Fraction, bool] | None = None
for integers in itertools.product(*(tuple(item) for item in options)):
rhs = tuple(Fraction(integer) - constant for integer, constant in zip(integers, constants))
solution = _linear_solve(matrix, rhs)
if solution is not None:
return solution, Fraction(0), True
approximation = _least_squares(matrix, rhs)
if approximation is None:
continue
residual = max(
(
abs(
sum(row[column] * approximation[column] for column in range(len(approximation)))
+ constant
- integer
)
for row, constant, integer in zip(matrix, constants, integers)
),
default=Fraction(0),
)
if best is None or residual < best[1]:
best = (approximation, residual, False)
return best
def _solve_for_transform(
equations: tuple[_Equation, ...], transform: SubgroupTransform
) -> tuple[tuple[Fraction, ...], Fraction, bool] | None:
try:
return _solve_modular(equations)
except ValueError as error:
if "branch cap exceeded" not in str(error):
raise
raise ValueError(
f"exact modular lift solver branch cap exceeded for "
f"{transform.parent.setting} -> {transform.subgroup.setting}"
) from error
def _shift_basis(structure: ASUStructure) -> tuple[tuple[Fraction, ...], ...]:
record = structure.spacegroup.it_number
from httk.atomistic import data
return tuple(
tuple(Fraction(value) for value in vector)
for vector in data.spacegroup_subgroup_record(record)["continuous_normalizer"]["basis_vectors"]
)
def _equation(
parent_position: Any,
piece: Any,
child_orbit: _Orbit,
child_branch: int,
shift_basis: tuple[tuple[Fraction, ...], ...],
) -> _Equation:
parent_branch = parent_position.representative
parent_matrix = _matrix(parent_branch.operation.matrix)
parent_vector = _vector(parent_branch.operation.vector)
piece_matrix = _matrix(piece.operation.matrix)
piece_vector = _vector(piece.operation.vector)
image_matrix = _matmul(piece_matrix, parent_matrix)
free = parent_branch.free
coefficients = tuple(tuple(image_matrix[row][column] for column in free) for row in range(3))
shift_columns = tuple(tuple(-vector[row] for vector in shift_basis) for row in range(3))
matrix = tuple(coefficients[row] + shift_columns[row] for row in range(3))
image_constant = tuple(value + piece_vector[row] for row, value in enumerate(_matvec(piece_matrix, parent_vector)))
observed = _vector(child_orbit.position.branches[child_branch].coordinate(child_orbit.site.free_params))
return _Equation(matrix, tuple(a - b for a, b in zip(image_constant, observed)))
def _orbit_distance(predicted: FracVector, observed: tuple[FracVector, ...], cell: Cell) -> tuple[float, Fraction]:
best_distance: float | None = None
best_fraction = Fraction(0)
for target in observed:
difference = FracVector(_wrapped_tuple(predicted - target))
distance = math.sqrt(_cartesian_distance_squared(difference, cell))
if best_distance is None or distance < best_distance:
best_distance = distance
best_fraction = max((abs(value) for value in difference.to_fractions()), default=Fraction(0))
assert best_distance is not None
return best_distance, best_fraction
def _validate_candidate(
candidate: _Candidate,
solution: tuple[Fraction, ...],
parent_position: Any,
pieces: tuple[Any, ...],
orbits: tuple[_Orbit, ...],
shift_basis: tuple[tuple[Fraction, ...], ...],
translation_cosets: tuple[FracVector, ...],
cell: Cell,
tolerance: float,
parent_count: int,
) -> tuple[float, Fraction] | None:
shift = _matvec(_transpose(shift_basis), solution[parent_count:]) if shift_basis else (Fraction(0),) * 3
predicted_points: list[FracVector] = []
observed_points = [point for _, orbit_index in candidate.pieces for point in orbits[orbit_index].coordinates]
for piece_index, _ in candidate.pieces:
piece = pieces[piece_index]
for branch in parent_position.branches:
parent_point = branch.coordinate(solution[:parent_count])
image = piece.operation.apply(parent_point) - FracVector(shift)
predicted_points.extend((image + coset).normalize() for coset in translation_cosets)
worst_distance = 0.0
worst_fraction = Fraction(0)
for predicted in predicted_points:
distance, residual = _orbit_distance(predicted, tuple(observed_points), cell)
if distance > tolerance:
return None
worst_distance = max(worst_distance, distance)
worst_fraction = max(worst_fraction, residual)
for observed in observed_points:
distance, residual = _orbit_distance(observed, tuple(predicted_points), cell)
if distance > tolerance:
return None
worst_distance = max(worst_distance, distance)
worst_fraction = max(worst_fraction, residual)
return worst_distance, worst_fraction
def _multiplicity_possible(structure: ASUStructure, transform: SubgroupTransform) -> bool:
targets: dict[str, dict[str, int]] = defaultdict(lambda: defaultdict(int))
for site in structure.wyckoff_sites:
targets[site.species][site.wyckoff] += 1
labels = tuple(sorted({label for target in targets.values() for label in target}))
for target in targets.values():
need = tuple(target[label] for label in labels)
vectors = []
for pieces in transform.splittings.values():
counts: defaultdict[str, int] = defaultdict(int)
for piece in pieces:
counts[piece.letter] += 1
vectors.append(tuple(counts[label] for label in labels))
vectors = [vector for vector in vectors if any(vector)]
seen: set[tuple[int, ...]] = set()
def reachable(remaining: tuple[int, ...], *, _seen=seen, _vectors=vectors) -> bool:
if not any(remaining):
return True
if remaining in _seen:
return False
_seen.add(remaining)
for vector in _vectors:
if all(value <= residual for value, residual in zip(vector, remaining)) and reachable(
tuple(residual - value for residual, value in zip(remaining, vector))
):
return True
return False
if not reachable(need):
return False
return True
def _metric_requirements(system: str) -> tuple[tuple[tuple[int, int], ...], tuple[tuple[int, Fraction], ...]]:
right_angle = Fraction(90)
if system == "monoclinic":
return (), ((0, right_angle), (2, right_angle))
if system == "orthorhombic":
return (), ((0, right_angle), (1, right_angle), (2, right_angle))
if system == "tetragonal":
return ((0, 1),), ((0, right_angle), (1, right_angle), (2, right_angle))
if system == "cubic":
return ((0, 1), (0, 2)), ((0, right_angle), (1, right_angle), (2, right_angle))
if system in ("trigonal", "hexagonal"):
# Every tabulated trigonal/hexagonal parent is reached in its hexagonal-axes standard
# setting (all trigonal R groups are the ":H" settings), so the metric constraint is
# a=b with alpha=beta=90, gamma=120. No rhombohedral-axes (":R") parent occurs.
return ((0, 1),), ((0, right_angle), (1, right_angle), (2, Fraction(120)))
return (), ()
def _cell_for_transform(structure: ASUStructure, transform: SubgroupTransform, tolerance: float) -> _MetricCell | None:
child_matrix = transform.operation.matrix.T()
parent_basis = SurdVector(child_matrix.inv()) * structure.cell.basis
measured = Cell(parent_basis, precision=structure.cell.precision, periodicity=structure.cell.periodicity)
equal_lengths, fixed_angles = _metric_requirements(transform.parent.crystal_system)
if not equal_lengths and not fixed_angles:
return _MetricCell(measured, 0.0, Fraction(0))
metric = measured.metric()
exact = all(metric._element((pair[0], pair[0])) == metric._element((pair[1], pair[1])) for pair in equal_lengths)
exact = exact and all(measured.angles[index] == angle for index, angle in fixed_angles)
if exact:
return _MetricCell(measured, 0.0, Fraction(0))
lengths = measured.lengths
cartesian = 0.0
fractional = Fraction(0)
for first, second in equal_lengths:
difference = abs(float(lengths[first]) - float(lengths[second]))
cartesian = max(cartesian, difference)
first_value = (
lengths[first]._rational_fraction() if lengths[first].is_rational else lengths[first].to_fractions_approx()
)
second_value = (
lengths[second]._rational_fraction()
if lengths[second].is_rational
else lengths[second].to_fractions_approx()
)
fractional = max(fractional, abs(first_value - second_value) / max(first_value, second_value, Fraction(1)))
for index, angle in fixed_angles:
measured_angle = measured.angles[index]
if measured_angle == angle:
continue
other = ((1, 2), (0, 2), (0, 1))[index]
cartesian = max(
cartesian,
min(float(lengths[other[0]]), float(lengths[other[1]]))
* abs(math.sin(math.radians(float(measured_angle - angle)))),
)
fractional = max(fractional, abs(measured_angle - angle) / Fraction(180))
if cartesian > tolerance:
return None
if any(not length.is_rational for length in lengths):
return None
params = [length._rational_fraction() for length in lengths] + list(measured.angles)
for first, second in equal_lengths:
params[second] = params[first]
for index, angle in fixed_angles:
params[3 + index] = angle
snapped = Cell(
CellParams(params).basis,
precision=structure.cell.precision,
periodicity=structure.cell.periodicity,
)
return _MetricCell(snapped, cartesian, fractional)
def _parent_charge(structure: ASUStructure, transform: SubgroupTransform) -> Fraction | None:
"""Undo descent's exact child/parent content scaling for one hop."""
if structure.charge is None:
return None
multiplier = abs(transform.operation.determinant())
if not multiplier:
raise ValueError(f"singular subgroup transform {transform.parent.setting} -> {transform.subgroup.setting}")
return structure.charge / multiplier
def _translation_cosets(transform: SubgroupTransform) -> tuple[FracVector, ...]:
"""Return child-coordinate representatives of parent-lattice translations."""
matrix = transform.operation.matrix.T().inv()
return finite_translation_cosets(tuple(FracVector(row) for row in matrix.to_fractions()))
def _candidate_list(
structure: ASUStructure,
transform: SubgroupTransform,
tolerance: float,
parent_cell: Cell,
) -> tuple[_Candidate, ...]:
orbits = tuple(
_Orbit(
index,
site,
structure.spacegroup.wyckoff_position(site.wyckoff),
tuple(
FracVector(point).normalize()
for point in structure.spacegroup.wyckoff_position(site.wyckoff).coordinates(site.free_params)
),
)
for index, site in enumerate(structure.wyckoff_sites)
)
shift_vectors = _shift_basis(structure)
translation_cosets = _translation_cosets(transform)
parent_labels = tuple(sorted(transform.splittings))
target_counts: dict[str, int] = defaultdict(int)
for orbit in orbits:
target_counts[orbit.site.wyckoff] += 1
candidates: list[_Candidate] = []
for parent_label in parent_labels:
parent_position = transform.parent.wyckoff_position(parent_label)
pieces = transform.splittings[parent_label]
piece_counts: dict[str, int] = defaultdict(int)
for piece in pieces:
piece_counts[piece.letter] += 1
if any(label not in target_counts for label in piece_counts):
continue
piece_labels = {piece.letter for piece in pieces}
species_values = sorted({orbit.site.species for orbit in orbits})
for species in species_values:
choices = [
tuple(
orbit.index
for orbit in orbits
if orbit.site.species == species and orbit.site.wyckoff in piece_labels
)
for piece in pieces
]
if any(not choice for choice in choices):
continue
order = sorted(range(len(pieces)), key=lambda index: (len(choices[index]), index))
selected: list[tuple[int, int]] = []
def visit(
depth: int,
*,
_order=order,
_selected=selected,
_pieces=pieces,
_parent_label=parent_label,
_species=species,
_parent_position=parent_position,
_choices=choices,
) -> None:
if depth == len(_order):
selected_labels = [orbits[index].site.wyckoff for _, index in _selected]
if sorted(selected_labels) != sorted(piece.letter for piece in _pieces):
return
# The tabulated split maps each child orbit into the parent orbit, but the anchor
# need not close on the child's branch-0 representative: a correspondence that only
# solves on another branch was silently lost when this was hard-coded to branch 0.
# Enumerate every child anchor branch (deterministic product order); any spurious
# modular match this admits is removed by the exact descent round trip in
# _lift_transform, which is the authoritative correctness gate.
# ponytail: full branch product costs O(branches^pieces) solver calls -- negligible
# for the one-branch P-lattice majority, but it lengthens the failed-lift search on
# high-multiplicity trigonal parents; dedup combos by their solved placement if that
# class of input matters.
branch_counts = [len(orbits[index].position.branches) for _, index in _selected]
for combo in itertools.product(*(range(count) for count in branch_counts)):
equations = tuple(
_equation(
_parent_position,
_pieces[piece_index],
orbits[orbit_index],
combo[position],
shift_vectors,
)
for position, (piece_index, orbit_index) in enumerate(_selected)
)
candidate = _Candidate(
_parent_label,
_species,
tuple(sorted(_selected)),
equations,
frozenset(index for _, index in _selected),
)
solved = _solve_for_transform(candidate.equations, transform)
if solved is None:
continue
check = _validate_candidate(
candidate,
solved[0],
_parent_position,
_pieces,
orbits,
shift_vectors,
translation_cosets,
structure.cell,
tolerance,
_parent_position.free_count,
)
if check is not None:
candidates.append(candidate)
return
piece_index = _order[depth]
for orbit_index in _choices[piece_index]:
_selected.append((piece_index, orbit_index))
visit(depth + 1)
_selected.pop()
visit(0)
return tuple(candidates)
def _embed_equation(
equation: _Equation, parent_offset: int, parent_count: int, total_parent: int, shift_dim: int
) -> _Equation:
"""Place one per-parent equation into a complete-cover variable vector."""
width = total_parent + shift_dim
rows = []
for row in equation.matrix:
embedded = [Fraction(0)] * width
embedded[parent_offset : parent_offset + parent_count] = row[:parent_count]
embedded[total_parent:] = row[parent_count:]
rows.append(tuple(embedded))
return _Equation(tuple(rows), equation.constant)
def _solve_with_fixed_shift(
equations: tuple[_Equation, ...], total_parent: int, shift_dim: int, transform: SubgroupTransform
) -> tuple[Fraction, ...] | None:
"""Solve a second deterministic gauge when continuous shift freedom exists."""
if not shift_dim:
return None
coefficients = tuple(Fraction(1, 2) if index == 0 else Fraction(0) for index in range(shift_dim))
fixed: list[_Equation] = []
for equation in equations:
constants = tuple(
value + sum(row[total_parent + index] * coefficients[index] for index in range(shift_dim))
for row, value in zip(equation.matrix, equation.constant)
)
fixed.append(_Equation(tuple(row[:total_parent] for row in equation.matrix), constants))
solved = _solve_for_transform(tuple(fixed), transform)
return None if solved is None else solved[0] + coefficients
@lru_cache(maxsize=100_000)
def _canonical_representative_cached(
spacegroup: Spacegroup, wyckoff: str, free_params: tuple[Fraction, ...]
) -> tuple[str, tuple[Fraction, ...]]:
"""Cached core of :func:`_canonical_representative`, keyed on hashable exact inputs.
``Spacegroup`` hashes and compares by setting, so states repeating a group and site share entries.
"""
position = spacegroup.wyckoff_position(wyckoff)
points = [FracVector(point).normalize() for point in position.coordinates(FracVector(free_params))]
least = min(points, key=lambda point: tuple(point.to_fractions()))
# A non-degenerate orbit has the declared position's full multiplicity. Its least member is
# therefore still on that same position, so recovering parameters there avoids scanning the
# complete Wyckoff table. A coordinate on a more-special position would collapse equivalent
# members; retain the full most-specific search for precisely that degenerate case.
if len({tuple(point.to_fractions()) for point in points}) == len(points):
parameters = position.parameters_of(least)
if parameters is not None:
return wyckoff, tuple(parameters.to_fractions())
match = spacegroup.identify_wyckoff(least)
if match is None:
return wyckoff, free_params
matched_position, parameters = match
return matched_position.letter, tuple(parameters.to_fractions())
def _canonical_representative(spacegroup: Spacegroup, site: WyckoffSite) -> tuple[str, Any]:
"""Return ``(letter, free_params)`` derived from the site's lexicographically least orbit point.
Two same-crystal descriptions can store different points of the same orbit (e.g. ``x`` versus
``-x``); re-identifying the least wrapped orbit point makes the reported letter and parameters a
function of the orbit alone.
"""
letter, parameters = _canonical_representative_cached(
spacegroup, site.wyckoff, tuple(site.free_params.to_fractions())
)
return letter, FracVector(parameters)
def _site_key(structure: ASUStructure) -> tuple[tuple[str, str, tuple[Fraction, ...]], ...]:
"""Return the exact translation- and orbit-representative-invariant sorted-site key of one ASU.
The key uses each site's canonical orbit representative, which may differ from the params stored
on a structure that has not passed :func:`_canonical_sites`, so ``backward_lift``-family ordering
follows this canonical key rather than the stored params.
"""
entries: list[tuple[str, str, tuple[Fraction, ...]]] = []
for site in structure.wyckoff_sites:
letter, parameters = _canonical_representative(structure.spacegroup, site)
entries.append((site.species, letter, tuple(parameters.to_fractions())))
return tuple(sorted(entries))
def _canonical_sites(structure: ASUStructure) -> ASUStructure:
"""Rewrite each site with the deterministic orbit representative from :func:`_canonical_representative`."""
sites: list[WyckoffSite] = []
changed = False
for site in structure.wyckoff_sites:
letter, parameters = _canonical_representative(structure.spacegroup, site)
if letter != site.wyckoff or tuple(parameters.to_fractions()) != tuple(site.free_params.to_fractions()):
changed = True
sites.append(WyckoffSite(letter, parameters, site.species))
ordered = sorted(
sites,
key=lambda site: (site.species, site.wyckoff, tuple(site.free_params.to_fractions())),
)
changed = changed or ordered != sites
if not changed:
return structure
try:
return ASUStructure(
structure.cell,
structure.spacegroup,
ordered,
structure.species,
transform=SettingTransform.identity(),
coordinate_precision=structure.coordinate_precision,
charge=structure.charge,
)
except ValueError:
return structure
def _canonical_result_key(result: LiftResult) -> tuple[Any, ...]:
sites = _site_key(result.asu)
path_key = tuple((transform.index, transform.subgroup_type) for transform in result.path)
metric = result.asu.cell.metric()
gram = tuple(metric._element((row, column)) for row in range(3) for column in range(3))
return (sites, tuple(result.shift.to_fractions()), gram, path_key)
def _crystals_match_within(left: ASUStructure, right: ASUStructure, tolerance: float) -> bool:
"""Return whether two same-cell crystals agree atom-for-atom within ``tolerance`` (Cartesian).
A bijective, species-labelled, minimum-image match -- the tolerant counterpart of
:func:`same_crystal` for the recognition-snapped path, where an accepted lift reproduces the
child only up to the hop's tolerance rather than exactly.
The basis-equality guard below is a defensive precondition, not a working branch: the round-trip
caller (:func:`_round_trip_reproduces`) always passes two structures built on the child's own
cell, so it never fires there. It is kept so the tolerant match stays sound if reused with
mismatched cells, where a Cartesian comparison would be meaningless.
"""
from httk.atomistic.models.structure.unitcell_view import UnitcellStructureView
left_view = UnitcellStructureView(left)
right_view = UnitcellStructureView(right)
if left_view.cell.basis != right_view.cell.basis:
return False
left_points = list(zip(left_view.species_at_sites, left_view.sites.reduced_coords.to_fractions()))
remaining = list(zip(right_view.species_at_sites, right_view.sites.reduced_coords.to_fractions()))
if len(left_points) != len(remaining):
return False
for name, coordinate in left_points:
match_index: int | None = None
for index, (other_name, other_coordinate) in enumerate(remaining):
if other_name != name:
continue
difference = FracVector(_wrapped_tuple(FracVector(coordinate) - FracVector(other_coordinate)))
# ponytail: greedy first-match within tolerance; sites sit far more than one tolerance
# apart, so the first admissible partner is the only one -- no assignment search needed.
if math.sqrt(_cartesian_distance_squared(difference, left_view.cell)) <= tolerance:
match_index = index
break
if match_index is None:
return False
remaining.pop(match_index)
return True
def _round_trip_reproduces(
child: ASUStructure, parent: ASUStructure, transform: SubgroupTransform, tolerance: float
) -> bool:
"""Return whether descending ``parent`` through ``transform`` reproduces ``child``.
A backward-lift candidate is correct exactly when the descent it claims to invert recovers the
child crystal. Descent (:func:`_child_sites`) is exact and already trusted, so this is the
authoritative correctness gate: the per-orbit distance check in :func:`_validate_candidate` is
only a cheap pre-filter, because a modular anchor match on any child branch can pass it while
placing the other sites on a different crystal. Exact input reproduces the child exactly, so the
fast path is :func:`same_crystal`; a recognition-snapped input is reproduced only up to the hop's
Cartesian ``tolerance``, so the fallback matches within it -- one rule, exact when it can be.
The lift's continuous-normalizer origin freedom (:attr:`LiftResult.shift`, and the residual
origin a polar or P1 child carries) leaves the descended child at a canonical origin the raw
input need not share. Both sides are reduced through :func:`_translation_normal_form` first, so
the comparison is invariant to exactly that freedom -- a no-op for the non-polar majority, where
the continuous normalizer is trivial. It shifts only along continuous directions, so a genuinely
different crystal from a spurious branch match still cannot be reconciled.
"""
try:
rebuilt = ASUStructure(
child.cell,
transform.subgroup,
_child_sites(parent, transform),
child.species,
transform=SettingTransform.identity(),
coordinate_precision=child.coordinate_precision,
charge=child.charge,
)
except ValueError:
return False
rebuilt = _translation_normal_form(rebuilt)
target = _translation_normal_form(child)
return same_crystal(rebuilt, target) or _crystals_match_within(rebuilt, target, tolerance)
def _lift_transform(structure: ASUStructure, transform: SubgroupTransform, tolerance: float) -> tuple[LiftResult, ...]:
if not _multiplicity_possible(structure, transform):
return ()
metric_cell = _cell_for_transform(structure, transform, tolerance)
if metric_cell is None:
return ()
parent_cell = metric_cell.cell
candidates = _candidate_list(structure, transform, tolerance, parent_cell)
if not candidates:
return ()
by_orbit: dict[int, list[int]] = defaultdict(list)
for index, candidate in enumerate(candidates):
for orbit in candidate.covered:
by_orbit[orbit].append(index)
all_orbits = frozenset(range(len(structure.wyckoff_sites)))
shift_vectors = _shift_basis(structure)
translation_cosets = _translation_cosets(transform)
shift_dim = len(shift_vectors)
orbits = tuple(
_Orbit(
index,
site,
structure.spacegroup.wyckoff_position(site.wyckoff),
tuple(
FracVector(point).normalize()
for point in structure.spacegroup.wyckoff_position(site.wyckoff).coordinates(site.free_params)
),
)
for index, site in enumerate(structure.wyckoff_sites)
)
results: list[LiftResult] = []
used: set[int] = set()
chosen: list[_Candidate] = []
def search(remaining: frozenset[int]) -> None:
if not remaining:
offsets: list[tuple[int, int]] = []
total_parent = 0
for candidate in chosen:
count = transform.parent.wyckoff_position(candidate.parent_letter).free_count
offsets.append((total_parent, count))
total_parent += count
equations = tuple(
_embed_equation(equation, offsets[index][0], offsets[index][1], total_parent, shift_dim)
for index, candidate in enumerate(chosen)
for equation in candidate.equations
)
solved = _solve_for_transform(equations, transform)
if solved is None:
return
solution, _, _ = solved
parent_parameters: list[tuple[str, FracVector, str]] = []
for candidate, (offset, count) in zip(chosen, offsets, strict=True):
parent_parameters.append(
(candidate.parent_letter, FracVector(solution[offset : offset + count]), candidate.species)
)
sites = tuple(WyckoffSite(letter, params, species) for letter, params, species in parent_parameters)
try:
asu = ASUStructure(
parent_cell,
transform.parent,
sites,
structure.species,
transform=SettingTransform.identity(),
coordinate_precision=structure.coordinate_precision,
charge=_parent_charge(structure, transform),
)
except ValueError:
return
shift = _matvec(_transpose(shift_vectors), solution[total_parent:]) if shift_vectors else (Fraction(0),) * 3
residual = metric_cell.fractional_deviation
for candidate, (offset, count) in zip(chosen, offsets, strict=True):
position = transform.parent.wyckoff_position(candidate.parent_letter)
check = _validate_candidate(
candidate,
solution[offset : offset + count] + solution[total_parent:],
position,
transform.splittings[candidate.parent_letter],
orbits,
shift_vectors,
translation_cosets,
structure.cell,
tolerance,
count,
)
if check is None:
return
residual = max(residual, check[1])
results.append(LiftResult(asu, transform.parent, (transform,), FracVector(shift), residual))
alternate = _solve_with_fixed_shift(equations, total_parent, shift_dim, transform)
if alternate is not None and alternate != solution:
alternate_sites = []
for candidate, (offset, count) in zip(chosen, offsets, strict=True):
alternate_sites.append(
WyckoffSite(
candidate.parent_letter,
FracVector(alternate[offset : offset + count]),
candidate.species,
)
)
try:
alternate_asu = ASUStructure(
parent_cell,
transform.parent,
alternate_sites,
structure.species,
transform=SettingTransform.identity(),
coordinate_precision=structure.coordinate_precision,
charge=_parent_charge(structure, transform),
)
except ValueError:
alternate_asu = None
if alternate_asu is not None:
alternate_shift = _matvec(_transpose(shift_vectors), alternate[total_parent:])
alternate_residual = metric_cell.fractional_deviation
valid = True
for candidate, (offset, count) in zip(chosen, offsets, strict=True):
check = _validate_candidate(
candidate,
alternate[offset : offset + count] + alternate[total_parent:],
transform.parent.wyckoff_position(candidate.parent_letter),
transform.splittings[candidate.parent_letter],
orbits,
shift_vectors,
translation_cosets,
structure.cell,
tolerance,
count,
)
if check is None:
valid = False
break
alternate_residual = max(alternate_residual, check[1])
if valid:
results.append(
LiftResult(
alternate_asu,
transform.parent,
(transform,),
FracVector(alternate_shift),
alternate_residual,
)
)
return
pivot = min(remaining, key=lambda item: sum(index not in used for index in by_orbit[item]))
for candidate_index in by_orbit[pivot]:
if candidate_index in used:
continue
candidate = candidates[candidate_index]
if not candidate.covered <= remaining:
continue
used.add(candidate_index)
chosen.append(candidate)
search(remaining - candidate.covered)
chosen.pop()
used.remove(candidate_index)
search(all_orbits)
# Authoritative correctness gate: keep only candidates whose exact descent reproduces the child.
# Branch-free anchoring above admits more modular matches; the ones that do not round-trip -- the
# unsound lifts branch-0 anchoring used to mask by never proposing them -- are removed here.
results = [result for result in results if _round_trip_reproduces(structure, result.asu, transform, tolerance)]
deduplicated: dict[tuple[Any, ...], LiftResult] = {}
for result in results:
deduplicated.setdefault(_canonical_result_key(result), result)
return tuple(deduplicated[key] for key in sorted(deduplicated))
def _apply_normalizer(structure: ASUStructure, record: dict[str, Any]) -> ASUStructure | None:
"""Apply one exact child-setting normalizer and rematch its Wyckoff sites."""
return _apply_normalizer_operation(structure, AffineOperation.from_record(record))
def _apply_normalizer_operation(
structure: ASUStructure,
operation: AffineOperation,
*,
trusted: bool = False,
) -> ASUStructure | None:
"""Re-express a standard-setting ASU through a point map, or ``None`` if the image is invalid.
Sites and basis transform JOINTLY (points by the operation, basis by its inverse-transpose), so
the Cartesian crystal is invariant under ANY invertible map -- a wrong or inapplicable operation
can only fail the per-orbit re-identification below and strand, never corrupt the emitted
structure; the descent round-trip gate downstream stays authoritative regardless.
"""
sites: list[WyckoffSite] = []
for site in structure.wyckoff_sites:
position = structure.spacegroup.wyckoff_position(site.wyckoff)
if trusted:
representative = FracVector(position.representative.coordinate(site.free_params)).normalize()
transformed_representative = operation.apply_wrapped(representative)
# A trusted affine normalizer maps a complete orbit bijectively, hence preserves its
# multiplicity. Search only that multiplicity tier while retaining the table's
# most-specific ordering within the tier.
first_match = None
for candidate_position in structure.spacegroup.wyckoff:
if candidate_position.multiplicity != position.multiplicity:
continue
candidate_parameters = candidate_position.parameters_of(transformed_representative)
if candidate_parameters is not None:
first_match = candidate_position, candidate_parameters
break
if first_match is None:
return None
first_position, first_parameters = first_match
# Callers use this only for tabulated or algebraically derived members of the affine
# normalizer. Such an operation maps a complete orbit bijectively; identifying one
# representative and preserving its multiplicity is therefore sufficient. The generic
# path below retains exact full-orbit set validation for arbitrary proposed operations.
sites.append(WyckoffSite(first_position.letter, first_parameters, site.species))
continue
original = tuple(FracVector(point).normalize() for point in position.coordinates(site.free_params))
transformed = tuple(operation.apply_wrapped(point) for point in original)
if not transformed:
return None
# One identified member proposes the target orbit; exact set equality below is already the
# authoritative proof that every transformed member lands in that complete orbit. Calling
# ``identify_wyckoff`` separately for all remaining members repeated a full Wyckoff-table
# scan up to hundreds of times per site and dominated high-symmetry normal forms.
first_match = structure.spacegroup.identify_wyckoff(transformed[0])
if first_match is None:
return None
first_position, first_parameters = first_match
expected = {
tuple(FracVector(point).normalize().to_fractions())
for point in first_position.coordinates(first_parameters)
}
actual = {tuple(point.to_fractions()) for point in transformed}
if expected != actual:
return None
sites.append(WyckoffSite(first_position.letter, first_parameters, site.species))
basis = SurdVector(operation.matrix.T().inv()) * structure.cell.basis
try:
return ASUStructure(
Cell(basis, precision=structure.cell.precision, periodicity=structure.cell.periodicity),
structure.spacegroup,
sites,
structure.species,
transform=SettingTransform.identity(),
coordinate_precision=structure.coordinate_precision,
charge=structure.charge,
)
except ValueError:
return None
def _enantiomorph(structure: ASUStructure) -> ASUStructure | None:
"""Flip a standard-setting structure in the HIGHER enantiomorph to its LOWER-numbered partner.
The genuine chirality flip negates-and-wraps every fractional coordinate (``f -> (-f) mod 1``,
i.e. Cartesian ``r -> -r``, an improper map) with the cell basis LEFT UNCHANGED -- the Gram matrix
and basis handedness are untouched -- and swaps the group to its enantiomorphic partner
(``it_number_enantiomorphic``). This differs from :func:`_apply_normalizer_operation`, which
transforms sites and basis JOINTLY and re-identifies against the SAME group and so cannot express
this map: here the mirrored sites are re-identified against the PARTNER group and orbit
completeness is verified there.
Returns ``None`` (no flip) unless the structure is in the higher member of a pair, so it only ever
normalizes higher to lower. Also returns ``None`` -- keeping the group -- when any site carries a
moment, since transforming an axial vector under an improper map is out of scope; magnetic
enantiomorphs are therefore left in their own group. Precondition: ``structure`` is in the IT
standard setting (every enantiomorphic group is single-setting, and both call sites guarantee it).
"""
partner_number = structure.spacegroup.record["it_number_enantiomorphic"]
if partner_number is None or partner_number >= structure.spacegroup.it_number:
return None
if any(site.moment is not None for site in structure.wyckoff_sites):
return None
partner = Spacegroup.standard(partner_number)
negate = AffineOperation(FracVector(((-1, 0, 0), (0, -1, 0), (0, 0, -1))), (0, 0, 0))
sites: list[WyckoffSite] = []
for site in structure.wyckoff_sites:
position = structure.spacegroup.wyckoff_position(site.wyckoff)
original = tuple(FracVector(point).normalize() for point in position.coordinates(site.free_params))
transformed = tuple(negate.apply_wrapped(point) for point in original)
matches = [partner.identify_wyckoff(point) for point in transformed]
if not matches or any(match is None for match in matches):
return None
first_match = matches[0]
if first_match is None:
return None
first_position, first_parameters = first_match
if any(match[0].letter != first_position.letter for match in matches[1:] if match is not None):
return None
expected = {
tuple(FracVector(point).normalize().to_fractions())
for point in first_position.coordinates(first_parameters)
}
actual = {tuple(point.to_fractions()) for point in transformed}
if expected != actual:
return None
sites.append(WyckoffSite(first_position.letter, first_parameters, site.species))
try:
return ASUStructure(
structure.cell,
partner,
sites,
structure.species,
transform=SettingTransform.identity(),
coordinate_precision=structure.coordinate_precision,
charge=structure.charge,
)
except ValueError:
return None
def _integer_diagonalize(rows: tuple[tuple[int, ...], ...]) -> tuple[tuple[int, ...], list[list[int]]]:
"""Diagonalize an integer ``n x 3`` matrix by unimodular row and column operations.
Returns the positive diagonal entries ``(d0, d1, d2)`` and the accumulated column transform
``V`` (``3 x 3`` unimodular), so that ``M t`` is integral iff ``d_i * (V^-1 t)_i`` is integral;
a zero diagonal marks a continuous direction (column of ``V``). Euclidean row/column swaps
drive each pivot to the gcd, so this terminates on any integer input.
"""
work = [list(row) for row in rows]
height = len(work)
transform = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
for pivot in range(3):
while True:
if pivot >= height or work[pivot][pivot] == 0:
swapped = False
for row in range(pivot, height):
for column in range(pivot, 3):
if work[row][column]:
if row != pivot:
work[pivot], work[row] = work[row], work[pivot]
if column != pivot:
for line in work:
line[pivot], line[column] = line[column], line[pivot]
for line in transform:
line[pivot], line[column] = line[column], line[pivot]
swapped = True
break
if swapped:
break
if not swapped:
break
leader = work[pivot][pivot]
settled = True
for row in range(pivot + 1, height):
if work[row][pivot]:
factor = work[row][pivot] // leader
for column in range(3):
work[row][column] -= factor * work[pivot][column]
if work[row][pivot]:
work[pivot], work[row] = work[row], work[pivot]
settled = False
break
if not settled:
continue
for column in range(pivot + 1, 3):
if work[pivot][column]:
factor = work[pivot][column] // leader
for line in work:
line[column] -= factor * line[pivot]
for line in transform:
line[column] -= factor * line[pivot]
if work[pivot][column]:
for line in work:
line[pivot], line[column] = line[column], line[pivot]
for line in transform:
line[pivot], line[column] = line[column], line[pivot]
settled = False
break
if settled:
break
diagonal = tuple(abs(work[index][index]) if index < height else 0 for index in range(3))
return diagonal, transform
def _ext_gcd(first: int, second: int) -> tuple[int, int, int]:
"""Return ``(g, x, y)`` with ``g = gcd`` positive and ``x*first + y*second = g``."""
old_r, r = first, second
old_s, s = 1, 0
old_t, t = 0, 1
while r:
quotient = old_r // r
old_r, r = r, old_r - quotient * r
old_s, s = s, old_s - quotient * s
old_t, t = t, old_t - quotient * t
if old_r < 0:
return -old_r, -old_s, -old_t
return old_r, old_s, old_t
def _lattice_basis(generators: list[tuple[int, ...]]) -> list[list[int]]:
"""Return a 3x3 integer basis (rows) of the rank-3 sublattice of ``Z^3`` spanned by generators."""
pivots: list[list[int]] = []
leads: list[int] = []
for generator in generators:
vector = list(generator)
for pivot, lead in zip(pivots, leads):
if vector[lead]:
gcd, left, right = _ext_gcd(pivot[lead], vector[lead])
combined = [left * pivot[k] + right * vector[k] for k in range(3)]
reduced = [(vector[lead] // gcd) * pivot[k] - (pivot[lead] // gcd) * vector[k] for k in range(3)]
pivot[:] = combined
vector = reduced
head = next((k for k in range(3) if vector[k]), None)
if head is not None:
if vector[head] < 0:
vector = [-value for value in vector]
pivots.append(vector)
leads.append(head)
order = sorted(range(len(leads)), key=lambda index: leads[index])
pivots = [pivots[index] for index in order]
leads = [leads[index] for index in order]
basis = [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
for pivot, lead in zip(pivots, leads):
basis[lead] = pivot
return basis
def _translation_lattice(spacegroup: Spacegroup) -> tuple[tuple[Fraction, ...], ...]:
"""Return a 3x3 rational matrix whose columns are a basis of the group's translation lattice.
The lattice is ``Z^3`` plus the centring vectors -- the translation parts of the ``W = I``
symmetry operations -- so for a P lattice this is the identity and for a centred lattice it is
finer than ``Z^3``.
"""
identity = FracVector.eye((3, 3))
centrings = [
tuple(FracVector(operation.vector).normalize().to_fractions())
for operation in spacegroup.symmetry_operations
if operation.matrix == identity
]
denominators = [value.denominator for vector in centrings for value in vector]
scale = math.lcm(*denominators) if denominators else 1
generators = [tuple(int(value * scale) for value in vector) for vector in centrings]
generators += [tuple(scale if i == j else 0 for j in range(3)) for i in range(3)]
rows = _lattice_basis(generators)
vectors = [tuple(Fraction(value, scale) for value in row) for row in rows]
return tuple(tuple(vectors[column][row] for column in range(3)) for row in range(3))
_DISCRETE_NORMALIZER_CACHE: dict[str, tuple[tuple[Fraction, ...], ...]] = {}
def _discrete_normalizer_translations(spacegroup: Spacegroup) -> tuple[tuple[Fraction, ...], ...]:
"""Return the finite group of discrete Euclidean-normalizer translations, reps in ``[0, 1)^3``.
A translation ``t`` normalizes ``G`` iff ``(I - W) t`` lies in the group's translation lattice
``T_G`` for every distinct linear part ``W`` (conjugation sends ``(W, w)`` to
``(W, w + (I - W) t)``, and ops sharing ``W`` differ only by translations in ``T_G``). For a
centred lattice ``T_G`` is finer than ``Z^3``, so the criterion uses ``T_G`` rather than ``Z^3``
-- otherwise quarter-translations such as ``(1/4, 1/4, 1/4)`` for F-43m are missed. With ``C``
a column basis of ``T_G``, the substitution ``t = C u`` turns ``(I - W) t in T_G`` into the
integral ``M' = C^-1 (I - W) C`` acting on ``u in Z^3``; diagonalizing ``M'`` enumerates the
finite solution set exactly, mapped back by ``C``. Continuous null-space directions are excluded
here -- the continuous-translation quotient owns them. Always includes the identity.
"""
key = spacegroup.hall_entry
cached = _DISCRETE_NORMALIZER_CACHE.get(key)
if cached is not None:
return cached
lattice = _translation_lattice(spacegroup)
lattice_inverse = _rational_inverse(lattice)
seen: set[tuple[tuple[int, ...], ...]] = set()
rows: list[tuple[int, ...]] = []
for operation in spacegroup.symmetry_operations:
linear = tuple(tuple(int(value) for value in row) for row in operation.matrix.to_fractions())
if linear in seen:
continue
seen.add(linear)
i_minus_w = tuple(tuple(Fraction((1 if i == j else 0) - linear[i][j]) for j in range(3)) for i in range(3))
conjugated = _matmul(_matmul(lattice_inverse, i_minus_w), lattice)
for row in conjugated:
assert all(value.denominator == 1 for value in row)
rows.append(tuple(value.numerator for value in row))
diagonal, transform = _integer_diagonalize(tuple(rows))
axes = [[Fraction(step, divisor) for step in range(divisor)] if divisor else [Fraction(0)] for divisor in diagonal]
representatives: set[tuple[Fraction, ...]] = set()
for first in axes[0]:
for second in axes[1]:
for third in axes[2]:
coefficients = (first, second, third)
unit = tuple(
sum((transform[row][column] * coefficients[column] for column in range(3)), Fraction(0))
for row in range(3)
)
representatives.add(
tuple(
sum((lattice[row][column] * unit[column] for column in range(3)), Fraction(0)) % 1
for row in range(3)
)
)
result = tuple(sorted(representatives))
_DISCRETE_NORMALIZER_CACHE[key] = result
return result
def _translation_normal_form(structure: ASUStructure) -> ASUStructure:
"""Return the origin-canonical image under the group's continuous-normalizer translations.
The continuous normalizer basis spans the directions along which the whole structure may be
translated while staying a valid same-group description. Candidate origins are, for every
expanded orbit point of every site, the pure translation that cancels that point's components
along the continuous directions; the least exact sorted-site key wins, and the identity
translation is always a candidate so the result never regresses. Pure translations leave the
cell unchanged, so the site key alone is a sound comparison.
"""
shift_basis = _shift_basis(structure)
if not shift_basis:
return structure
# Every tabulated continuous-normalizer vector is an axis-aligned unit vector, so the continuous
# directions are simply the axes those vectors point along.
axes = sorted({index for vector in shift_basis for index in range(3) if vector[index]})
identity = FracVector.eye((3, 3))
candidates: set[tuple[Fraction, ...]] = {(Fraction(0), Fraction(0), Fraction(0))}
for site in structure.wyckoff_sites:
position = structure.spacegroup.wyckoff_position(site.wyckoff)
for point in position.coordinates(site.free_params):
values = FracVector(point).normalize().to_fractions()
translation = [Fraction(0), Fraction(0), Fraction(0)]
for index in axes:
translation[index] = (-values[index]) % 1
candidates.add(tuple(translation))
if structure.spacegroup.it_number == 1:
# In P1 every site is its own general-position orbit and all three origin directions are
# free. The generic path below needlessly re-identifies all N sites through the Wyckoff
# tables for each of the N candidate origins (and this is repeated for every boundary-metric
# automorphism). Compare the translated rows directly, then build only the winning image.
# This is the same exact _site_key: P1's canonical orbit representative is the wrapped point.
def p1_key(candidate: tuple[Fraction, ...]) -> tuple[tuple[str, str, tuple[Fraction, ...]], ...]:
shift = FracVector(candidate)
return tuple(
sorted(
(
site.species,
"a",
tuple((site.free_params + shift).normalize().to_fractions()),
)
for site in structure.wyckoff_sites
)
)
winner = min(candidates, key=p1_key)
if not any(winner):
return structure
shift = FracVector(winner)
return ASUStructure(
structure.cell,
structure.spacegroup,
[
WyckoffSite("a", (site.free_params + shift).normalize(), site.species, moment=site.moment)
for site in structure.wyckoff_sites
],
structure.species,
transform=SettingTransform.identity(),
coordinate_precision=structure.coordinate_precision,
charge=structure.charge,
)
best = structure
best_key = _site_key(structure)
for candidate in sorted(candidates):
if not any(candidate):
continue
image = _apply_normalizer_operation(
structure,
AffineOperation(identity, FracVector(candidate)),
trusted=True,
)
if image is None:
continue
key = _site_key(image)
if key < best_key:
best, best_key = image, key
return best
@lru_cache(maxsize=50_000)
def _metric_automorphism_operations(
gram: tuple[tuple[Fraction, ...], ...],
) -> tuple[AffineOperation, ...]:
"""Return every integral automorphism of one exact rational three-dimensional Gram matrix.
A Niggli-reduced Gram can still sit on a boundary with finitely many equally reduced bases.
For a row basis change ``B' = U B`` those are precisely the integral ``U`` satisfying
``U G U.T = G``. Each row ``u`` has fixed squared length ``u G u.T = G_ii`` and Cauchy in
the ``G`` metric gives ``u_k**2 <= G_ii * (G^-1)_kk``. Those exact bounds make the search
finite without a heuristic radius. The affine point map corresponding to ``U`` is
``U.T().inv()``, because :func:`_apply_normalizer_operation` constructs
``B' = M.T().inv() B`` for point matrix ``M``.
This is deliberately a metric stabilizer, not a search over adjacent Niggli domains: it closes
the finite ambiguity of a fixed reduced Gram (cubic has 48 signed permutations) while keeping
generic P1 input finite and exact.
"""
# Rational LDL^T decomposition. Enumerating the ellipsoid from the last coordinate backwards
# lets each assigned coordinate tighten the remaining radius. The final coordinate is solved as
# an exact rational square, so an elongated cell costs linearly rather than as the volume of its
# axis-aligned bounding box.
lower = [[Fraction(int(row == column)) for column in range(3)] for row in range(3)]
diagonal = [Fraction(0), Fraction(0), Fraction(0)]
for column in range(3):
diagonal[column] = gram[column][column] - sum(
(lower[column][prior] * lower[column][prior] * diagonal[prior] for prior in range(column)),
Fraction(0),
)
for row in range(column + 1, 3):
lower[row][column] = (
gram[row][column]
- sum(
(lower[row][prior] * lower[column][prior] * diagonal[prior] for prior in range(column)),
Fraction(0),
)
) / diagonal[column]
def exact_square_roots(value: Fraction) -> tuple[Fraction, ...]:
if value < 0:
return ()
numerator = math.isqrt(value.numerator)
denominator = math.isqrt(value.denominator)
if numerator * numerator != value.numerator or denominator * denominator != value.denominator:
return ()
root = Fraction(numerator, denominator)
return (root,) if root == 0 else (-root, root)
def candidates(target: Fraction) -> tuple[tuple[int, ...], ...]:
vector = [0, 0, 0]
found: list[tuple[int, ...]] = []
def search(index: int, used: Fraction) -> None:
remaining = target - used
if remaining < 0:
return
offset = sum((lower[row][index] * vector[row] for row in range(index + 1, 3)), Fraction(0))
bound = remaining / diagonal[index]
if index == 0:
for root in exact_square_roots(bound):
exact_coordinate = root - offset
if exact_coordinate.denominator == 1:
vector[0] = exact_coordinate.numerator
found.append(tuple(vector))
return
radius = math.isqrt(bound.numerator // bound.denominator)
# ``sqrt(bound) < radius + 1``. Widening by two around ``-offset`` gives a compact
# conservative integer interval; the exact inequality below removes its fringe.
centre_floor = (-offset).numerator // (-offset).denominator
for integer_coordinate in range(centre_floor - radius - 2, centre_floor + radius + 3):
term = diagonal[index] * (Fraction(integer_coordinate) + offset) ** 2
if term > remaining:
continue
vector[index] = integer_coordinate
search(index - 1, used + term)
search(2, Fraction(0))
return tuple(found)
row_candidates = tuple(candidates(gram[index][index]) for index in range(3))
operations: list[AffineOperation] = []
for first in row_candidates[0]:
for second in row_candidates[1]:
if _bilinear(gram, first, second) != gram[0][1]:
continue
for third in row_candidates[2]:
if _bilinear(gram, first, third) != gram[0][2] or _bilinear(gram, second, third) != gram[1][2]:
continue
rows = (first, second, third)
if abs(_integer_determinant(rows)) != 1:
continue
operations.append(AffineOperation(FracVector(rows).T().inv(), (0, 0, 0)))
return tuple(operations)
def _p1_metric_automorphism_operations(structure: ASUStructure) -> tuple[AffineOperation, ...]:
"""Return exact finite Gram stabilizers for a P1 cell, or none for a non-rational metric."""
metric = structure.cell.metric()
if not metric.is_rational:
return ()
gram = tuple(tuple(value for value in row) for row in metric.coefficient(1).to_fractions())
return _metric_automorphism_operations(gram)
def _normalizer_metric_key(structure: ASUStructure, operation: AffineOperation) -> tuple[Any, ...]:
"""Return the exact Gram key after one joint affine re-expression, without rematching sites."""
basis_change = SurdVector(operation.matrix.T().inv())
metric = basis_change * structure.cell.metric() * basis_change.T()
return tuple(metric._element((row, column)) for row in range(3) for column in range(3))
def _basis_key(basis: SurdVector) -> tuple[Any, ...]:
"""Return Cartesian basis components in deterministic row-major order."""
return tuple(basis._element((row, column)) for row in range(3) for column in range(3))
def _demote_sites(structure: ASUStructure) -> ASUStructure:
"""Re-label any site whose expanded orbit degenerates onto a more-special Wyckoff position.
A lift can leave an atom on a special coordinate while still carrying a general (or less
special) Wyckoff letter; its expanded orbit then contains coincident points. Re-identifying the
representative coordinate through the exact Wyckoff machinery demotes it to its true most-special
letter and free params, so equal crystals carry equal site keys.
"""
demoted: list[WyckoffSite] = []
changed = False
for site in structure.wyckoff_sites:
position = structure.spacegroup.wyckoff_position(site.wyckoff)
points = [FracVector(point).normalize() for point in position.coordinates(site.free_params)]
if len({tuple(point.to_fractions()) for point in points}) == len(points):
demoted.append(site)
continue
match = structure.spacegroup.identify_wyckoff(points[0])
if match is None:
demoted.append(site)
continue
letter_position, parameters = match
demoted.append(WyckoffSite(letter_position.letter, parameters, site.species))
changed = True
if not changed:
return structure
try:
return ASUStructure(
structure.cell,
structure.spacegroup,
demoted,
structure.species,
transform=SettingTransform.identity(),
coordinate_precision=structure.coordinate_precision,
charge=structure.charge,
)
except ValueError:
return structure
def _normal_form(structure: ASUStructure) -> ASUStructure:
"""Return the compact normal-form state used while traversing upward lifts.
Search states retain the established, setting-compatible affine-coset quotient: it removes
origin and ordinary normalizer duplicates without changing a cell into an alternative valid
setting that could strand a later tabulated lift. The broad coset-times-group metric quotient
is intentionally deferred to :func:`_terminal_normal_form`, where no further lift needs the
current setting. P1 has only the identity space-group operation, so its finite exact Gram
stabilizer is safe and closes boundary-equivalent Niggli bases already at entry.
"""
structure = _demote_sites(structure)
identity_matrix = FracVector.eye((3, 3))
inversion_matrix = FracVector(((-1, 0, 0), (0, -1, 0), (0, 0, -1)))
trusted_inversion = _resetting_preserves_group(structure.spacegroup, inversion_matrix)
operations = [AffineOperation.identity()]
if structure.spacegroup.it_number == 1:
operations.extend(_p1_metric_automorphism_operations(structure))
else:
try:
record = data.affine_normalizer_coset_record(structure.spacegroup.hall_entry)
except KeyError:
record = None
if record is not None:
system = structure.spacegroup.crystal_system
operations.extend(
AffineOperation.from_record(coset)
for coset in record.get("affine_normalizer_cosets", ())
if system in coset["compatible_systems"]
)
operations = list(dict.fromkeys(operations))
translations = _discrete_normalizer_translations(structure.spacegroup)
best: ASUStructure | None = None
best_key: tuple[Any, ...] | None = None
for operation in operations:
image = _apply_normalizer_operation(structure, operation)
if image is None:
continue
for translation in translations:
if any(translation):
shifted = _apply_normalizer_operation(image, AffineOperation(identity_matrix, FracVector(translation)))
else:
shifted = image
if shifted is None:
continue
reduced = _translation_normal_form(shifted)
if reduced.cell.basis.det().sign() < 0:
flipped = _apply_normalizer_operation(
reduced,
AffineOperation(inversion_matrix, (0, 0, 0)),
trusted=trusted_inversion,
)
if flipped is not None:
reduced = _translation_normal_form(flipped)
key = _site_key(reduced)
if best_key is None or key < best_key:
best, best_key = reduced, key
return _canonical_sites(best) if best is not None else _canonical_sites(structure)
def _terminal_normal_form(structure: ASUStructure) -> ASUStructure:
"""Return a deterministic canonical representative of one state within its own space group.
Sites mislabeled onto a special coordinate are demoted first, then the normalizer quotients are
collapsed: the candidate images are every tabulated affine-normalizer coset expanded by every
member of the space group, crossed with every discrete Euclidean-normalizer translation (all
including the identity), each followed by the continuous-translation quotient. The data records
are representatives modulo the space group, but applying one member jointly transforms the basis
and the coordinates, so applying only the representative is not equivalent to applying its whole
coset here. Nor may the fixed-axes ``compatible_systems`` annotation prune this enumeration:
a valid normalizer can re-choose the axes of an equal-metric lattice (notably at Niggli boundary
metrics) even when it does not preserve the generic metric written in those axes.
Every accepted image is the same Cartesian crystal re-expressed in the same group --
``_apply_normalizer_operation`` re-identifies the Wyckoff sites and rejects anything that is not
-- so collapsing them to the least exact joint metric/site key cannot lose a reachable terminal
and yields a description-invariant representative. The least key wins in deterministic operation
and translation order; each site is finally stored at its canonical orbit representative.
Basis-choice invariance for a generic same lattice comes separately from the Niggli reduction of
triclinic (SG 1/2) entries in :func:`highest_symmetry`; this finite enumeration additionally
resolves its boundary-equivalent reduced bases.
"""
structure = _demote_sites(structure)
identity_matrix = FracVector.eye((3, 3))
inversion_matrix = FracVector(((-1, 0, 0), (0, -1, 0), (0, 0, -1)))
trusted_inversion = _resetting_preserves_group(structure.spacegroup, inversion_matrix)
representatives = [AffineOperation.identity()]
if structure.spacegroup.it_number == 1:
# A P1 Niggli entry has no point-group operations to supply the finite stabilizer of a
# boundary Gram. Enumerate that exact GL(3, Z) stabilizer directly instead of relying on
# the bounded, system-labelled affine table.
representatives.extend(_p1_metric_automorphism_operations(structure))
else:
try:
record = data.affine_normalizer_coset_record(structure.spacegroup.hall_entry)
except KeyError:
record = None
if record is not None:
representatives.extend(
AffineOperation.from_record(coset) for coset in record.get("affine_normalizer_cosets", ())
)
representatives = list(dict.fromkeys(representatives))
# This quotient factors each space-group member as a basis-only point operation below. Members
# that share a matrix and differ only by a centring translation therefore produce identical
# candidates; keep one per matrix (preferring the true identity for the final tie-break).
point_operations: dict[tuple[tuple[Fraction, ...], ...], AffineOperation] = {}
for operation in structure.spacegroup.symmetry_operations:
matrix_key = tuple(tuple(row) for row in operation.matrix.to_fractions())
if matrix_key not in point_operations or operation.is_identity():
point_operations[matrix_key] = operation
# Metric is the primary normal-form key and neither a discrete nor a continuous translation
# changes it. Select its least exact tier before doing expensive Wyckoff rematching. Every
# candidate is a tabulated normalizer (or a checked P1 Gram automorphism), so an operation in a
# higher metric tier cannot beat a successful member of this one.
by_metric: dict[tuple[Any, ...], list[tuple[AffineOperation, AffineOperation]]] = {}
for representative in representatives:
for group_operation in point_operations.values():
operation = representative * group_operation
metric_key = _normalizer_metric_key(structure, operation)
by_metric.setdefault(metric_key, []).append((representative, group_operation))
translations = _discrete_normalizer_translations(structure.spacegroup)
best: ASUStructure | None = None
best_key: tuple[Any, ...] | None = None
translation_normal_forms: dict[tuple[Any, ...], ASUStructure] = {}
def translated_normal_form(candidate: ASUStructure) -> ASUStructure:
"""Memoize the pure translation quotient within this terminal orbit."""
# Every candidate below derives from the same immutable input and therefore shares group,
# species definitions, precision, charge, and periodicity. The exact basis and immutable
# Wyckoff-site tuple fully distinguish the inputs on which the translation quotient depends.
cache_key = (
candidate.cell.basis,
tuple(
(site.wyckoff, tuple(site.free_params.to_fractions()), site.species) for site in candidate.wyckoff_sites
),
)
cached = translation_normal_forms.get(cache_key)
if cached is None:
cached = _translation_normal_form(candidate)
translation_normal_forms[cache_key] = cached
return cached
for metric_key in sorted(by_metric):
grouped: dict[AffineOperation, list[AffineOperation]] = {}
for representative, group_operation in by_metric[metric_key]:
grouped.setdefault(representative, []).append(group_operation)
for representative, group_operations in grouped.items():
image = _apply_normalizer_operation(structure, representative, trusted=True)
if image is None:
continue
for translation in translations:
if any(translation):
shifted = _apply_normalizer_operation(
image,
AffineOperation(identity_matrix, FracVector(translation)),
trusted=True,
)
else:
shifted = image
if shifted is None:
continue
# A coset or discrete shift can change which origin is canonical, so re-run the
# continuous-translation quotient on each candidate.
reduced = translated_normal_form(shifted)
# The table is quotient-ed by G. Because C normalizes G, the left-coset elements
# g*C cover the same quotient class as C*g. Applying g after C leaves every complete
# transformed Wyckoff orbit unchanged as a set, but selects the jointly transformed
# basis g.T^-1 C.T^-1 B. Factor that basis-only step here so Wyckoff rematching and
# the continuous origin gauge run once per C/translation rather than once per g*C.
# The sites are likewise shared by every factored group operation. Compute their
# inversion lazily, at most once, only if a det=-1 basis actually needs it: the full-
# orbit rematch is expensive for groups whose operations are all orientation-preserving.
site_keys = {id(reduced): _site_key(reduced)}
inverted_reduced: ASUStructure | None = None
inversion_attempted = False
for group_operation in group_operations:
basis = SurdVector(group_operation.matrix.T().inv()) * reduced.cell.basis
site_source = reduced
if basis.det().sign() < 0:
if not inversion_attempted:
inverted_reduced = _apply_normalizer_operation(
reduced,
AffineOperation(inversion_matrix, (0, 0, 0)),
trusted=trusted_inversion,
)
inversion_attempted = True
if inverted_reduced is not None:
inverted_reduced = translated_normal_form(inverted_reduced)
site_keys[id(inverted_reduced)] = _site_key(inverted_reduced)
if inverted_reduced is not None:
basis = -basis
site_source = inverted_reduced
# The site key and handedness can be decided without rebuilding an ASUStructure;
# group operations in this factored loop change only the basis. Construct just a
# candidate that can improve the current best, rather than hundreds of identical-
# site temporary structures per metric tier.
handedness_key = 0 if basis.det().sign() > 0 else 1
basis_key = _basis_key(basis)
identity_key = 0 if representative.is_identity() and group_operation.is_identity() else 1
key = (metric_key, site_keys[id(site_source)], handedness_key, basis_key, identity_key)
if best_key is not None and key >= best_key:
continue
try:
candidate = ASUStructure(
Cell(
basis,
precision=site_source.cell.precision,
periodicity=site_source.cell.periodicity,
),
site_source.spacegroup,
site_source.wyckoff_sites,
site_source.species,
transform=SettingTransform.identity(),
coordinate_precision=site_source.coordinate_precision,
charge=site_source.charge,
)
except ValueError:
continue
# Prefer a right-handed representative when the orbit supplies one. The raw
# Cartesian components are only the final tie-break after metric, sites, and
# handedness. Normally the orientation rebuild replaces them with a function of
# the Gram matrix. When that rebuild would require unsupported nested radicals,
# this tie-break still collapses the finite point-group frame ambiguity exposed
# by equivalent basis scrambles in the input's shared Cartesian frame.
if best_key is None or key < best_key:
best, best_key = candidate, key
if best is not None:
break
# Store each site at its canonical orbit representative so the returned free params, not just the
# comparison key, are independent of which orbit point the input happened to carry.
return _canonical_sites(best) if best is not None else _canonical_sites(structure)
def _normalizer_retries(structure: ASUStructure, target: Spacegroup, tolerance: float) -> tuple[LiftResult, ...]:
try:
record = data.affine_normalizer_coset_record(structure.spacegroup.hall_entry)
except KeyError:
return ()
transforms = subgroup_transforms(target, structure.spacegroup)
results: list[LiftResult] = []
images: dict[tuple[Any, ...], ASUStructure | None] = {}
for transform in transforms:
if not _multiplicity_possible(structure, transform):
continue
for coset in record.get("affine_normalizer_cosets", ()):
if target.crystal_system not in coset["compatible_systems"]:
continue
affine = coset["affine_transformation"]
image_key = (tuple(tuple(row) for row in affine["matrix"]), tuple(affine["vector"]))
if image_key not in images:
images[image_key] = _apply_normalizer(structure, coset)
image = images[image_key]
if image is None:
continue
matches = tuple(
result
for result in _lift_transform(image, transform, tolerance)
if result.path and result.path[0] == transform
)
operation = AffineOperation.from_record(coset)
correction = transform.operation * operation.inverse() * transform.operation.inverse()
for match in matches:
restored = _apply_normalizer_operation(match.asu, correction)
if restored is None:
continue
results.append(LiftResult(restored, match.spacegroup, match.path, match.shift, match.residual))
deduplicated: dict[tuple[Any, ...], LiftResult] = {}
for result in results:
deduplicated.setdefault(_canonical_result_key(result), result)
return tuple(deduplicated[key] for key in sorted(deduplicated))
def _bilinear(gram: tuple[tuple[Fraction, ...], ...], left: tuple[int, ...], right: tuple[int, ...]) -> Fraction:
return sum((left[i] * gram[i][j] * right[j] for i in range(3) for j in range(3)), Fraction(0))
def _integer_determinant(rows: tuple[tuple[int, ...], ...]) -> int:
a, b, c = rows
return a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) + a[2] * (b[0] * c[1] - b[1] * c[0])
@lru_cache(maxsize=50_000)
def _search_conventional_basis(
gram: tuple[tuple[Fraction, ...], ...], system: str
) -> tuple[tuple[int, ...], ...] | None:
"""Return a canonical integer basis of ``gram``'s lattice meeting the system's exact metric.
The lattice is the candidate parent lattice ``inv(M^T)*B_child`` (as its exact rational gram);
the returned rows are integer coordinates on it whose cell satisfies the parent crystal system's
:func:`_metric_requirements` exactly. Enumeration is a bounded exact short-vector search --
``bound`` 5 was what the trigonal Bi case needed in the parameterization these grams use, kept as
a documented bounded search rather than a proven radius. The choice is deterministic: least
``(abs(det), row-order squared lengths, rows)``, then flipped to a right-handed (positive
determinant) basis so the applied re-expression is orientation-preserving. ``None`` when no such
basis is found in the bound.
"""
equal_lengths, fixed_angles = _metric_requirements(system)
if not equal_lengths and not fixed_angles:
return None
everything = [n for n in itertools.product(range(-5, 6), repeat=3) if any(n)]
everything.sort(key=lambda n: (_bilinear(gram, n, n), n))
norms = sorted({_bilinear(gram, n, n) for n in everything})
cap = norms[min(len(norms) - 1, 9)]
short = [n for n in everything if _bilinear(gram, n, n) <= cap]
right = Fraction(90)
# Every target system that reaches here has alpha = beta = 90 (parents are orthorhombic or more
# symmetric), so the c axis is perpendicular to whichever earlier axes those angles name -- a
# precomputed perpendicular set keeps the innermost loop to genuine candidates instead of all of
# ``short``, which is what makes the exact search tractable inside the breadth-first search.
perpendicular = {vector: {other for other in short if _bilinear(gram, vector, other) == 0} for vector in short}
c_perpendicular_to_a = any(index == 1 and degrees == right for index, degrees in fixed_angles)
c_perpendicular_to_b = any(index == 0 and degrees == right for index, degrees in fixed_angles)
basal_angles = tuple((index, degrees) for index, degrees in fixed_angles if index == 2)
def angle_ok(rows: tuple[tuple[int, ...], ...], index: int, degrees: Fraction) -> bool:
first, second = ((1, 2), (0, 2), (0, 1))[index]
product = _bilinear(gram, rows[first], rows[second])
if degrees == right:
return product == 0
length = _bilinear(gram, rows[first], rows[first])
return 2 * product == -length and length == _bilinear(gram, rows[second], rows[second])
best: tuple[tuple[Any, ...], tuple[tuple[int, ...], ...]] | None = None
for a in short:
length_a = _bilinear(gram, a, a)
for b in short:
# Prune on the constraints that involve only the first two axes before the c loop.
if (0, 1) in equal_lengths and _bilinear(gram, b, b) != length_a:
continue
if any(not angle_ok((a, b, a), index, degrees) for index, degrees in basal_angles):
continue
candidates: Any = short
if c_perpendicular_to_a:
candidates = perpendicular[a]
if c_perpendicular_to_b:
candidates = candidates & perpendicular[b] if c_perpendicular_to_a else perpendicular[b]
for c in candidates:
rows: tuple[tuple[int, ...], ...] = (a, b, c)
if _integer_determinant(rows) == 0:
continue
if not all(
_bilinear(gram, rows[i], rows[i]) == _bilinear(gram, rows[j], rows[j]) for i, j in equal_lengths
):
continue
if not all(angle_ok(rows, index, degrees) for index, degrees in fixed_angles):
continue
key = (
abs(_integer_determinant(rows)),
tuple(_bilinear(gram, rows[k], rows[k]) for k in range(3)),
rows,
)
if best is None or key < best[0]:
best = (key, rows)
if best is None:
return None
rows = best[1]
if _integer_determinant(rows) < 0:
# The abs(det) key leaves the handedness of the winner undetermined; negate one axis so the
# returned basis is right-handed and every re-expression built from it preserves orientation.
rows = (rows[0], rows[1], tuple(-value for value in rows[2]))
return rows
def _exact_rational_matrix(matrix: Any) -> tuple[tuple[Fraction, ...], ...]:
return tuple(tuple(value for value in row) for row in matrix.to_fractions())
def _rational_determinant(matrix: tuple[tuple[Fraction, ...], ...]) -> Fraction:
a, b, c = matrix
return a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) + a[2] * (b[0] * c[1] - b[1] * c[0])
@lru_cache(maxsize=50_000)
def _primitive_conventional_bases(
gram: tuple[tuple[Fraction, ...], ...], system: str
) -> tuple[tuple[tuple[int, ...], ...], ...]:
"""Return ALL integer bases on ``gram``'s lattice meeting the system's exact metric, in canonical order.
The centred-supercell recell arm enumerates candidate parent conventional cells over the child's
PRIMITIVE lattice, so unlike :func:`_search_conventional_basis` it must return every metric-valid
basis -- which candidate is usable depends on the transform (the implied child re-expression must
be a lattice automorphism), a filter only the caller can apply. The ``c`` axis is additionally
allowed past the short-vector norm cap when it is fully perpendicular-constrained: a hexagonal
conventional ``c`` (e.g. Bi-166's c = 11.8 vs a = 4.5) is far longer than the primitive vectors,
and perpendicularity to both basal axes is itself a strong exact filter. Bounded search
(entries in [-5, 5]), same heritage as :func:`_search_conventional_basis`.
"""
equal_lengths, fixed_angles = _metric_requirements(system)
if not equal_lengths and not fixed_angles:
return ()
everything = [n for n in itertools.product(range(-5, 6), repeat=3) if any(n)]
everything.sort(key=lambda n: (_bilinear(gram, n, n), n))
norms = sorted({_bilinear(gram, n, n) for n in everything})
cap = norms[min(len(norms) - 1, 9)]
short = [n for n in everything if _bilinear(gram, n, n) <= cap]
right = Fraction(90)
c_perpendicular_to_a = any(index == 1 and degrees == right for index, degrees in fixed_angles)
c_perpendicular_to_b = any(index == 0 and degrees == right for index, degrees in fixed_angles)
basal_angles = tuple((index, degrees) for index, degrees in fixed_angles if index == 2)
def angle_ok(rows: tuple[tuple[int, ...], ...], index: int, degrees: Fraction) -> bool:
first, second = ((1, 2), (0, 2), (0, 1))[index]
product = _bilinear(gram, rows[first], rows[second])
if degrees == right:
return product == 0
length = _bilinear(gram, rows[first], rows[first])
return 2 * product == -length and length == _bilinear(gram, rows[second], rows[second])
found: list[tuple[tuple[int, ...], ...]] = []
for a in short:
length_a = _bilinear(gram, a, a)
for b in short:
if (0, 1) in equal_lengths and _bilinear(gram, b, b) != length_a:
continue
if any(not angle_ok((a, b, a), index, degrees) for index, degrees in basal_angles):
continue
if c_perpendicular_to_a and c_perpendicular_to_b:
candidates = [n for n in everything if _bilinear(gram, a, n) == 0 and _bilinear(gram, b, n) == 0]
else:
candidates = short
for c in candidates:
rows = (a, b, c)
if _integer_determinant(rows) == 0:
continue
if not all(
_bilinear(gram, rows[i], rows[i]) == _bilinear(gram, rows[j], rows[j]) for i, j in equal_lengths
):
continue
if not all(angle_ok(rows, index, degrees) for index, degrees in fixed_angles):
continue
found.append(rows)
found.sort(
key=lambda rows: (
abs(_integer_determinant(rows)),
tuple(_bilinear(gram, rows[k], rows[k]) for k in range(3)),
rows,
)
)
return tuple(found)
@lru_cache(maxsize=10_000)
def _resetting_preserves_group(spacegroup: Spacegroup, point_matrix: FracVector) -> bool:
"""Whether the pure-linear point map ``x -> A x`` conjugates the group's op set onto itself exactly.
The centred recell arm's re-expressions are lattice automorphisms unimodular over the centred
translation lattice; expressed on the conventional cell they MAY be half-integer (a supported
generalization -- no probed centred group has yet produced a non-integer map that passes this
check, so in practice the accepted maps have been integer), so op-set preservation is verified
directly rather than assumed from integrality: the conjugate image set
``{(A W A^-1, wrap(A w))}`` must EQUAL the tabulated op set of the standard setting -- set
equality, not membership, since a non-injective conjugation could collapse two operations onto
one and still pass a subset check. This is the exact-op-set contract the SettingTransform
machinery pins across all 527 settings, applied here to a candidate re-choice.
"""
inverse = point_matrix.inv()
operations = {
(
tuple(tuple(value for value in row) for row in operation.matrix.to_fractions()),
tuple(FracVector(operation.vector).normalize().to_fractions()),
)
for operation in spacegroup.symmetry_operations
}
conjugates = set()
for operation in spacegroup.symmetry_operations:
conjugated = point_matrix * operation.matrix * inverse
translated = point_matrix * FracVector(operation.vector)
conjugates.add(
(
tuple(tuple(value for value in row) for row in conjugated.to_fractions()),
tuple(FracVector(translated).normalize().to_fractions()),
)
)
return conjugates == operations
def _integer_recell_lifts(
structure: ASUStructure,
transform: SubgroupTransform,
child_gram: tuple[tuple[Fraction, ...], ...],
tolerance: float,
) -> tuple[LiftResult, ...]:
"""First recell arm: integer lattice-normalizer re-choice (see :func:`_recell_lifts`)."""
matrix_transpose = _exact_rational_matrix(transform.operation.matrix.T())
matrix_transpose_inverse = _exact_rational_matrix(transform.operation.matrix.T().inv())
parent_gram = _matmul3(
_matmul3(matrix_transpose_inverse, child_gram),
tuple(tuple(matrix_transpose_inverse[j][i] for j in range(3)) for i in range(3)),
)
basis = _search_conventional_basis(parent_gram, transform.parent.crystal_system)
if basis is None:
return ()
rechoice = _matmul3(
_matmul3(matrix_transpose, tuple(tuple(Fraction(v) for v in row) for row in basis)),
matrix_transpose_inverse,
)
if any(value.denominator != 1 for row in rechoice for value in row):
return ()
# A re-expression is a same-crystal re-choice only if it preserves the cell volume: an integer
# but volume-changing matrix (real runs produced det -2 and -6) is a supercell, not a lattice
# normalizer, and must not be applied even though the downstream re-identification would also
# reject it. The right-handed search basis makes the surviving determinant exactly +1.
if abs(_rational_determinant(rechoice)) != 1:
return ()
# Pre-check the rebased parent metric before paying for the re-expression and lift: applying
# the operation sets the child basis to inv(U)*B, so the parent cell it derives is
# inv(M^T)*inv(U)*B. Most searched re-choices do not clear the parent metric here (8 of 13
# firings in review), so this exact check -- the same one _cell_for_transform applies, on the
# exact rational grams this tier is guarded to -- skips them cheaply.
rebased_parent = SurdVector(transform.operation.matrix.T().inv()) * (
SurdVector(FracVector(rechoice).inv()) * structure.cell.basis
)
measured = Cell(rebased_parent)
equal_lengths, fixed_angles = _metric_requirements(transform.parent.crystal_system)
metric = measured.metric()
if not all(metric._element((i, i)) == metric._element((j, j)) for i, j in equal_lengths):
return ()
if not all(measured.angles[index] == angle for index, angle in fixed_angles):
return ()
transpose = FracVector(tuple(tuple(rechoice[j][i] for j in range(3)) for i in range(3)))
image = _apply_normalizer_operation(structure, AffineOperation(transpose, (0, 0, 0)))
if image is None:
return ()
return _lift_transform(image, transform, tolerance)
def _centred_recell_lifts(
structure: ASUStructure,
transform: SubgroupTransform,
child_gram: tuple[tuple[Fraction, ...], ...],
tolerance: float,
) -> tuple[LiftResult, ...]:
"""Second recell arm: centred-lattice re-choice of the child's own conventional cell.
An R-centred parent's conventional cell is an intrinsic supercell of the child's CONVENTIONAL
lattice, so the first arm's conventional-lattice search can never produce it. The two
load-bearing ingredients here are (a) candidate parent cells enumerated over the child's
PRIMITIVE lattice (conventional basis + centring vectors -- the det-3 hexagonal basis exists
only there), and (b) the discrete-normalizer-translation crossing of the applied image, which
places the atoms on the tabulated splitting's origin coset (observed: Bi's C2/m 'i' at
x = 1/2 vs the split's x = 0 family). A candidate ``M_int`` is usable for this transform iff
``X = M^T * M_int`` is integer (the new child cell expressed in primitive coordinates), giving
the child re-expression ``V = X * C^T`` with exact ``det V = +1`` demanded. ``V`` may in
general be half-integer on the conventional cell (unimodular over the centred lattice only) --
a supported generalization, though every map accepted so far, including Bi's landing one, has
been integer. ``V`` is verified to conjugate the child's op set onto itself exactly
(:func:`_resetting_preserves_group`), applied as the point map ``(V^T)^-1`` through the
validity-gated :func:`_apply_normalizer_operation`, then crossed with the discrete normalizer
translations before lifting. The committed descent round-trip gate in :func:`_lift_transform`
stays authoritative -- and since a joint site+basis transformation is Cartesian-invariant for
ANY invertible map, a wrong re-expression can only strand (fail identification or the gate),
never corrupt the emitted crystal.
Deterministic: candidate bases in canonical order, translations in sorted order, first landing
basis per transform wins. # ponytail: first-landing stop bounds cost; enumerate-and-dedup all
landing bases if a centred case ever needs a non-first representative.
"""
lattice = _translation_lattice(structure.spacegroup)
identity = tuple(tuple(Fraction(1) if row == column else Fraction(0) for column in range(3)) for row in range(3))
if lattice == identity:
return ()
lattice_transpose = tuple(tuple(lattice[row][column] for row in range(3)) for column in range(3))
primitive_gram = _matmul3(
_matmul3(lattice_transpose, child_gram),
tuple(tuple(lattice_transpose[j][i] for j in range(3)) for i in range(3)),
)
matrix_transpose = _exact_rational_matrix(transform.operation.matrix.T())
identity_matrix = FracVector.eye((3, 3))
results: list[LiftResult] = []
for basis in _primitive_conventional_bases(primitive_gram, transform.parent.crystal_system):
integer_candidate = tuple(tuple(Fraction(value) for value in row) for row in basis)
new_cell_primitive = _matmul3(matrix_transpose, integer_candidate)
if any(value.denominator != 1 for row in new_cell_primitive for value in row):
continue
rechoice = _matmul3(new_cell_primitive, lattice_transpose)
if _rational_determinant(rechoice) != 1:
continue
point_matrix = FracVector(rechoice).T().inv()
if not _resetting_preserves_group(structure.spacegroup, point_matrix):
continue
image = _apply_normalizer_operation(structure, AffineOperation(point_matrix, (0, 0, 0)))
if image is None:
continue
if _cell_for_transform(image, transform, tolerance) is None:
continue
for translation in _discrete_normalizer_translations(structure.spacegroup):
shifted = (
image
if not any(translation)
else _apply_normalizer_operation(image, AffineOperation(identity_matrix, FracVector(translation)))
)
if shifted is None:
continue
results.extend(_lift_transform(shifted, transform, tolerance))
if results:
break
return tuple(results)
def _matmul3(
left: tuple[tuple[Fraction, ...], ...], right: tuple[tuple[Fraction, ...], ...]
) -> tuple[tuple[Fraction, ...], ...]:
return tuple(
tuple(sum((left[i][k] * right[k][j] for k in range(3)), Fraction(0)) for j in range(3)) for i in range(3)
)
def _recell_lifts(structure: ASUStructure, target: Spacegroup, tolerance: float) -> tuple[LiftResult, ...]:
"""Third, fail-only lift tier: re-choose the child's conventional cell so a metric-rejected hop fits.
A tabulated transform is rejected by :func:`_cell_for_transform` when the child, carried through
the breadth-first search in a reduced cell, presents the parent lattice in an axis choice that
misses the parent's exact metric class -- even though the lattice does admit a conforming cell.
For each such transform this searches the candidate parent lattice for a conventional basis
meeting the parent metric exactly, derives the implied child re-expression ``U = M^T*N*inv(M^T)``,
and -- when ``U`` is a unimodular (``abs(det) == 1``) integer lattice normalizer of the child --
applies ``U^T`` through the validity-gated :func:`_apply_normalizer_operation` and lifts the
re-expressed child. A ``U`` that is non-integer or volume-changing is not an integer lattice
normalizer and is skipped by this first arm; for a CENTRED child (the R-centred trigonal case,
where the parent conventional cell is an intrinsic supercell of the child's conventional lattice)
the second arm :func:`_centred_recell_lifts` then searches the child's primitive lattice for a
half-integer, centred-lattice-unimodular re-choice instead. The committed descent round-trip
gate in :func:`_lift_transform` remains authoritative against the original child.
``_lift_transform`` is called directly, deliberately bypassing the same-setting
:data:`COMPATIBLE_CRYSTAL_SYSTEMS` filter that cell re-choice is designed to escape.
"""
child_gram_surd = SurdVector(structure.cell.basis) * SurdVector(structure.cell.basis).T()
if not child_gram_surd.is_rational:
return ()
child_gram = tuple(tuple(value for value in row) for row in child_gram_surd.to_fractions_approx())
results: list[LiftResult] = []
for transform in subgroup_transforms(target, structure.spacegroup):
if not _multiplicity_possible(structure, transform):
continue
if _cell_for_transform(structure, transform, tolerance) is not None:
continue
lifted = _integer_recell_lifts(structure, transform, child_gram, tolerance)
if not lifted:
# The integer arm cannot re-choose across a centring (the parent conventional cell may be
# an intrinsic supercell of the child's conventional lattice); the centred arm searches
# the child's primitive lattice for a half-integer lattice-unimodular re-choice instead.
lifted = _centred_recell_lifts(structure, transform, child_gram, tolerance)
results.extend(lifted)
deduplicated: dict[tuple[Any, ...], LiftResult] = {}
for result in results:
deduplicated.setdefault(_canonical_result_key(result), result)
return tuple(deduplicated[key] for key in sorted(deduplicated))
def _raw_lifts(structure: ASUStructure, target: Spacegroup, tolerance: float) -> tuple[LiftResult, ...]:
if structure.spacegroup.crystal_system not in COMPATIBLE_CRYSTAL_SYSTEMS[target.crystal_system]:
return ()
results: list[LiftResult] = []
for transform in subgroup_transforms(target, structure.spacegroup):
results.extend(_lift_transform(structure, transform, tolerance))
deduplicated: dict[tuple[Any, ...], LiftResult] = {}
for result in results:
deduplicated.setdefault(_canonical_result_key(result), result)
return tuple(deduplicated[key] for key in sorted(deduplicated))
[docs]
def backward_lift(
structure: ASUStructure,
supergroup: Spacegroup | int,
*,
tolerance: float | None = None,
) -> tuple[LiftResult, ...]:
"""Return all exact or tolerance-accepted lifts into one minimal supergroup.
:param structure: The child-group asymmetric unit to lift.
:param supergroup: The one-hop parent space group or IT number.
:param tolerance: Cartesian acceptance tolerance, or the recognition-derived default.
:return: Distinct parent representations in table order and canonical order.
:raises ValueError: If the input is unsupported or the target is not one hop above it.
A bounded normalizer retry applies tabulated cosets to child fractional coordinates and maps
successful results back with the exact inverse, in tabulated coset order.
"""
if not isinstance(structure, ASUStructure):
raise TypeError(f"expected ASUStructure, got {type(structure).__name__}")
require_full_periodicity(structure.cell, "backward_lift")
if any(site.moment is not None for site in structure.wyckoff_sites):
raise ValueError("backward_lift does not support structures with site moments")
if structure.assemblies is not None:
raise ValueError("backward_lift does not support structures with assemblies")
if structure.molecular:
raise ValueError("backward_lift does not support molecular structures")
current = _standard_input(structure)
target = supergroup if isinstance(supergroup, Spacegroup) else Spacegroup.standard(supergroup)
if target.it_number not in minimal_supergroups(current.spacegroup):
raise ValueError(f"space group {target.setting} is not a minimal supergroup of {current.spacegroup.setting}")
accepted_tolerance = structure_tolerance(current) if tolerance is None else float(tolerance)
results = list(_raw_lifts(current, target, accepted_tolerance))
results.extend(_normalizer_retries(current, target, accepted_tolerance))
deduplicated: dict[tuple[Any, ...], LiftResult] = {}
for result in results:
deduplicated.setdefault(_canonical_result_key(result), result)
return tuple(deduplicated[key] for key in sorted(deduplicated))
[docs]
def lift_candidates(structure: ASUStructure, *, tolerance: float | None = None) -> tuple[LiftResult, ...]:
"""Return all one-hop parent lifts in deterministic order.
:param structure: The child-group asymmetric unit to lift.
:param tolerance: Cartesian acceptance tolerance, or the recognition-derived default.
:return: Results ordered by parent IT number, table order, and exact key.
"""
parents = minimal_supergroups(structure.spacegroup)
results: list[tuple[int, int, LiftResult]] = []
for parent_number in parents:
transforms = subgroup_transforms(parent_number, structure.spacegroup)
table_indices = {transform: index for index, transform in enumerate(transforms)}
for result in backward_lift(structure, parent_number, tolerance=tolerance):
if result.path and result.path[0] in table_indices:
results.append((parent_number, table_indices[result.path[0]], result))
return tuple(
result for _, _, result in sorted(results, key=lambda item: (item[0], item[1], _canonical_result_key(item[2])))
)
def _highest_lifts(structure: ASUStructure, tolerance: float) -> tuple[LiftResult, ...]:
results: list[tuple[int, int, LiftResult]] = []
parents = minimal_supergroups(structure.spacegroup)
indices_by_parent: dict[int, dict[SubgroupTransform, int]] = {}
for parent_number in parents:
transforms = subgroup_transforms(parent_number, structure.spacegroup)
indices = {transform: index for index, transform in enumerate(transforms)}
indices_by_parent[parent_number] = indices
target = Spacegroup.standard(parent_number)
try:
raw = _raw_lifts(structure, target, tolerance)
# The state normal form now collapses the origin/normalizer-equivalent variants that
# unconditional retries used to contribute, so only fall back to retries when the direct
# lift found nothing for this parent (backward_lift keeps unconditional retries).
candidates = raw or _normalizer_retries(structure, target, tolerance)
except ValueError as error:
if "branch cap exceeded" not in str(error):
raise
# One hopeless parent target may saturate the modular solver's per-row product; skip it
# so the breadth-first search still explores the remaining minimal supergroups.
logging.getLogger(__name__).warning(
"skipping lift into %s: %s", target.setting, error, extra={"context": "symmetry"}
)
continue
for result in candidates:
if result.path and result.path[0] in indices:
results.append((parent_number, indices[result.path[0]], result))
if not results:
# Only a genuinely stranded state -- no direct or retry lift into any parent -- pays for the
# exact conventional-cell re-choice search, so its cost never touches states that climb
# normally (the search is dormant on every P-lattice battery).
for parent_number in parents:
for result in _recell_lifts(structure, Spacegroup.standard(parent_number), tolerance):
if result.path and result.path[0] in indices_by_parent[parent_number]:
results.append((parent_number, indices_by_parent[parent_number][result.path[0]], result))
deduplicated: dict[tuple[Any, ...], tuple[int, int, LiftResult]] = {}
for parent_number, table_index, result in results:
key = (parent_number, table_index, _canonical_result_key(result))
deduplicated.setdefault(key, (parent_number, table_index, result))
return tuple(
result
for _, _, result in sorted(
deduplicated.values(), key=lambda item: (item[0], item[1], _canonical_result_key(item[2]))
)
)
def _structure_signature(structure: ASUStructure) -> tuple[Any, ...]:
sites = []
for site in structure.wyckoff_sites:
position = structure.spacegroup.wyckoff_position(site.wyckoff)
orbit = tuple(sorted(_wrapped_tuple(point) for point in position.coordinates(site.free_params)))
sites.append((site.species, site.wyckoff, orbit))
metric = structure.cell.metric()
gram = tuple(metric._element((row, column)) for row in range(3) for column in range(3))
return tuple(sorted(sites)), gram
def _terminal_result(
structure: ASUStructure, path: tuple[SubgroupTransform, ...], shift: FracVector, residual: Fraction
) -> LiftResult:
return LiftResult(structure, structure.spacegroup, path, shift, residual)
def _exact_triangular_basis(metric: SurdVector, *, left_handed: bool) -> SurdVector | None:
"""Return an exact canonical Cartesian basis for a rational Gram matrix.
This is the three-dimensional Cholesky construction written in the surd field. Every square
root taken is of a positive rational Schur complement, so the result is exact even when no
rational-angle :class:`CellParams` round-trip can reproduce the Gram. A fixed Cartesian
reflection gives the corresponding left-handed orientation when required.
"""
if not metric.is_rational:
return None
gram = metric.coefficient(1).to_fractions()
# Keep the Cholesky arithmetic in Q until the three final square roots. Expressing, for
# example, g10/sqrt(g00) through generic SurdScalar division is algebraically sound but can
# build a very large multiradical field before cancellation discovers that the square is
# rational again. Strongly sheared (yet perfectly ordinary) cells made that route take
# minutes. The rationalized forms below are identical:
#
# g10/sqrt(g00) = (g10/g00) sqrt(g00)
#
# and likewise for the third row. Consequently every Schur complement is computed as one
# Fraction and exact work stays bounded independently of the input basis shear.
first_square = gram[0][0]
if first_square <= 0:
return None
second_x_factor = gram[1][0] / first_square
second_square = gram[1][1] - gram[1][0] * second_x_factor
if second_square <= 0:
return None
third_x_factor = gram[2][0] / first_square
third_y_numerator = gram[2][1] - gram[2][0] * second_x_factor
third_y_factor = third_y_numerator / second_square
third_square = gram[2][2] - gram[2][0] * third_x_factor - third_y_numerator * third_y_factor
if third_square <= 0:
return None
# SurdScalar canonicalizes a radical by extracting the square part of numerator*denominator.
# That is cheap for the small crystallographic radicands the exact path is intended for, but an
# approximate CIF angle can make a rational Gram with hundred-digit, essentially arbitrary
# radicands whose exact factorization is not a bounded operation. Perfect rational squares are
# always cheap; otherwise keep this optional orientation step to radicands small enough for the
# square-part implementation's trial division. The caller retains the already exact basis when
# this returns None.
def bounded_sqrt(value: Fraction) -> SurdScalar | None:
numerator_root = math.isqrt(value.numerator)
denominator_root = math.isqrt(value.denominator)
if (
numerator_root * numerator_root == value.numerator
and denominator_root * denominator_root == value.denominator
):
return SurdScalar(Fraction(numerator_root, denominator_root))
# Decimal-scaled conventional cells commonly have q = r*s^2 with a tiny crystallographic
# radicand (notably r=2 or 3) but a large raw numerator/denominator. Recognize those forms
# without factoring the large square coefficient at all.
for radicand in range(2, 49):
scaled = value / radicand
scaled_numerator_root = math.isqrt(scaled.numerator)
scaled_denominator_root = math.isqrt(scaled.denominator)
if (
scaled_numerator_root * scaled_numerator_root == scaled.numerator
and scaled_denominator_root * scaled_denominator_root == scaled.denominator
):
return SurdScalar.from_radicand_map(
{radicand: Fraction(scaled_numerator_root, scaled_denominator_root)}
)._as_scalar()
if (value.numerator * value.denominator).bit_length() <= 48:
return SurdScalar.sqrt_of(value)
return None
first = bounded_sqrt(first_square)
second = bounded_sqrt(second_square)
third = bounded_sqrt(third_square)
if first is None or second is None or third is None:
return None
second_x = first * SurdScalar(second_x_factor)
third_x = first * SurdScalar(third_x_factor)
third_y = second * SurdScalar(third_y_factor)
basis = SurdVector._from_scalar_grid(
[[first, SurdScalar(0), SurdScalar(0)], [second_x, second, SurdScalar(0)], [third_x, third_y, third]],
(3, 3),
)
if left_handed:
basis = basis * SurdVector(((1, 0, 0), (0, 1, 0), (0, 0, -1)))
return basis
def _canonical_orientation(structure: ASUStructure) -> ASUStructure:
"""Put the cell in the standard crystallographic orientation for its exact Gram matrix.
Equal crystals reached through different lifts can carry the same metric in different rotated
frames; fractional coordinates are orientation-independent, so re-expressing the cell in the
unique standard orientation of its parameters makes the returned basis a function of the metric
alone. A det<0 (left-handed) basis -- which a det=-1 normalizer coset can leave behind -- is
first re-expressed through inversion so the standard-orientation rebuild does not silently
produce the enantiomorph. In an enantiomorphic group inversion conjugates to the other group and
is rejected; the rebuild then uses a fixed left-handed Cartesian orientation instead. That cell
is related to any other left-handed realization of the same Gram by a proper global rotation, so
physical chirality is preserved while raw Cartesian-frame dependence is removed. The rebuild is
accepted only when it reproduces the Gram matrix exactly (the stored angles are approximate
fractions, so it need not). A non-rational Gram whose canonical triangular factor would require
nested radicals outside :class:`SurdVector`'s multiquadratic field cannot be reoriented exactly;
that uncommon fallback is idempotent but retains the input's global Cartesian rotation.
"""
left_handed = False
if structure.cell.basis.det().sign() < 0:
# Inversion as a basis change: basis -> -B (right-handed), coords -> -x, Cartesian identical.
inversion_matrix = FracVector(((-1, 0, 0), (0, -1, 0), (0, 0, -1)))
inverted = _apply_normalizer_operation(
structure,
AffineOperation(inversion_matrix, (0, 0, 0)),
trusted=_resetting_preserves_group(structure.spacegroup, inversion_matrix),
)
if inverted is None:
# Enantiomorphic group: inversion would change the crystal. Keep the handedness and use
# a fixed left-handed Cartesian orientation below.
left_handed = True
else:
structure = inverted
cell = structure.cell
metric = cell.metric()
oriented_basis = _exact_triangular_basis(metric, left_handed=left_handed)
if oriented_basis is None:
# A rational metric reaches this branch only when exact squarefree-radical normalization was
# deliberately bounded above. Asking for ``cell.lengths`` would immediately attempt the same
# unbounded factorization again, so retain the exact input orientation. Standard recognition
# and the lattice normal forms have already removed basis shear; this fallback only leaves a
# global Cartesian rotation uncanonicalized for unwieldy approximate-angle Gram matrices.
if metric.is_rational:
return structure
lengths = cell.lengths
if any(not length.is_rational for length in lengths):
return structure
params = [length._rational_fraction() for length in lengths] + list(cell.angles)
try:
oriented_basis = CellParams(params).basis
if left_handed:
oriented_basis = oriented_basis * SurdVector(((1, 0, 0), (0, 1, 0), (0, 0, -1)))
except ValueError:
return structure
try:
oriented = Cell(oriented_basis, precision=cell.precision, periodicity=cell.periodicity)
except ValueError:
return structure
if oriented.basis == cell.basis or oriented.metric() != metric:
return structure
return ASUStructure(
oriented,
structure.spacegroup,
structure.wyckoff_sites,
structure.species,
transform=SettingTransform.identity(),
coordinate_precision=structure.coordinate_precision,
charge=structure.charge,
)
def _exact_translations(atoms: list[tuple[str, tuple[Fraction, ...]]]) -> list[tuple[Fraction, ...]]:
"""Return every exact non-lattice translation that maps a P1 atom set onto itself.
Candidates are the differences within the smallest species class (any self-translation must send
that class's first atom to another of its atoms), and each is verified against the full per-species
wrapped-coordinate **multisets** by exact rational equality -- no tolerance. The multiset (not a
set) matters: a duplicated ``(species, coordinate)`` site finds no self-translation and passes
through untouched instead of reaching the reduction's internal consistency error.
"""
by_species: dict[str, Counter[tuple[Fraction, ...]]] = defaultdict(Counter)
for species, coordinate in atoms:
by_species[species][coordinate] += 1
smallest = min(by_species.values(), key=lambda counts: counts.total())
anchor = FracVector(next(iter(smallest)))
translations: list[tuple[Fraction, ...]] = []
for point in smallest:
candidate = (FracVector(point) - anchor).normalize()
if not any(candidate.to_fractions()):
continue
if all(
Counter(
{
tuple((FracVector(coordinate) + candidate).normalize().to_fractions()): count
for coordinate, count in counts.items()
}
)
== counts
for counts in by_species.values()
):
translations.append(tuple(candidate.to_fractions()))
return translations
def _primitive_reduced_entry(structure: ASUStructure) -> ASUStructure:
"""Collapse an exact P1 supercell to its unique primitive description before the search.
Any exact supercell -- any multiplicity, any (diagonal or sheared) sublattice orientation -- has
a translation group finer than the cell lattice. The finer lattice ``L'`` is generated by ``Z^3``
together with **all** exact self-translations. ``L'`` itself is unique, but the basis emitted here
(``_lattice_basis`` echelon output) is order-dependent, so it is *not* a canonical cell: this step
only reaches *a* primitive description. Description-invariance of the final result comes from the
downstream Niggli reduction and normal form, not from this helper, which must not be reused
standalone as a canonicalizer. With ``S`` the new lattice basis in old fractional units (rows),
the new cell is ``B' = S B`` and coordinates remap as ``c' = c S^-1``; each primitive atom then
appears exactly ``n = 1/|det S|`` times and is deduplicated. Fractional coordinate precision is
rescaled through the subdivision (conservatively, as one hop does) so ``cartesian_precision`` and
the derived tolerance are preserved.
Only exact rational invariance fires this: a noisy (near-coincident) supercell is left untouched
here -- snapping such input is :func:`~httk.atomistic.canonical_asu`'s tolerant job.
"""
atoms = [
(
site.species,
tuple(
FracVector(structure.spacegroup.wyckoff_position(site.wyckoff).coordinates(site.free_params)[0])
.normalize()
.to_fractions()
),
)
for site in structure.wyckoff_sites
]
if len(atoms) < 2:
return structure
translations = _exact_translations(atoms)
if not translations:
return structure
denominator = math.lcm(*(value.denominator for translation in translations for value in translation))
generators = [tuple(denominator if i == j else 0 for j in range(3)) for i in range(3)]
generators += [tuple(int(value * denominator) for value in translation) for translation in translations]
hermite = _lattice_basis(generators)
lattice = tuple(tuple(Fraction(hermite[i][j], denominator) for j in range(3)) for i in range(3))
determinant = (
hermite[0][0] * (hermite[1][1] * hermite[2][2] - hermite[1][2] * hermite[2][1])
- hermite[0][1] * (hermite[1][0] * hermite[2][2] - hermite[1][2] * hermite[2][0])
+ hermite[0][2] * (hermite[1][0] * hermite[2][1] - hermite[1][1] * hermite[2][0])
)
multiplicity, remainder = divmod(denominator**3, abs(determinant))
assert remainder == 0
inverse = _rational_inverse(lattice)
collapsed: dict[tuple[str, tuple[Fraction, ...]], int] = defaultdict(int)
for species, coordinate in atoms:
reduced = tuple((FracVector(coordinate) * FracVector(inverse)).normalize().to_fractions())
collapsed[(species, reduced)] += 1
if any(count != multiplicity for count in collapsed.values()):
raise ValueError("exact supercell reduction produced uneven site copies; this is a logic error")
sites = [WyckoffSite("a", FracVector(coordinate), species) for species, coordinate in collapsed]
charge = None if structure.charge is None else structure.charge / multiplicity
# B' = S * B, so scale precisions exactly as one subgroup hop does for a basis change: the cell's
# absolute precision by the row-sum factor of S, and the fractional coordinate precision by the
# column-sum factor of S^-1 -- the conservative bound that keeps cartesian_precision equal-or-wider.
matrix = FracVector(lattice)
return ASUStructure(
Cell(
SurdVector(lattice) * structure.cell.basis,
precision=_scaled_precision(structure.cell.precision, _matrix_row_sum_factor(matrix)),
periodicity=structure.cell.periodicity,
),
structure.spacegroup,
sites,
structure.species,
transform=SettingTransform.identity(),
coordinate_precision=_scaled_precision(
structure.coordinate_precision, _matrix_column_sum_factor(FracVector(inverse))
),
charge=charge,
)
def _pseudo_translations(structure: ASUStructure) -> list[tuple[Fraction, ...]]:
"""Full-cell exact self-translations lying OUTSIDE the group's own translation lattice.
For a declared-symmetry structure the full-cell atom set is always invariant under the group's
centring translations, so those are filtered out through ``_translation_lattice``; anything that
survives is genuine translation pseudosymmetry -- the declared description is an exact supercell.
Exact rational check only (a noisy near-supercell is untouched; that is ``canonical_asu``'s job).
"""
atoms: list[tuple[str, tuple[Fraction, ...]]] = []
for site in structure.wyckoff_sites:
position = structure.spacegroup.wyckoff_position(site.wyckoff)
for point in position.coordinates(site.free_params):
atoms.append((site.species, tuple(FracVector(point).normalize().to_fractions())))
if len(atoms) < 2:
return []
translations = _exact_translations(atoms)
if not translations:
return []
lattice_inverse = _rational_inverse(_translation_lattice(structure.spacegroup))
return [
translation
for translation in translations
if any(
sum(lattice_inverse[row][column] * translation[column] for column in range(3)).denominator != 1
for row in range(3)
)
]
def _isomorphic_reduced_entry(structure: ASUStructure) -> ASUStructure:
"""Collapse a declared-symmetry (SG >= 2) exact supercell via tabulated isomorphic self-lifts.
The cheap guard is :func:`_pseudo_translations`: one full-cell expansion plus an exact O(N^2)
translation check, empty for the overwhelmingly common non-supercell input -- measured
millisecond-scale to ~16 full-cell atoms and ~250 ms at ~100 atoms (about 9% of a bulk
``canonical_asu`` call at that size). When genuine pseudo-translations exist, their count
determines the total collapse index exactly (the pseudo-translation cosets number
``len(pseudo) + 1`` including identity), so backward lifts through same-IT-number entries --
the vendored Bärnighausen self-entries (tabulated only where isomorphic subgroups are maximal)
plus the dedicated ``isomorphic_subgroups_std`` transforms -- are FILTERED to indices dividing
that implied index (a self-lift of index k collapses a k-subgroup of the pseudo-translation
group, so any other index cannot land and only burns solver time) and tried ascending, so
composite supercells collapse stepwise through the cheapest applicable steps (tabulated
indices reach 9; prime factors 2, 3, 5, 7 compose by iteration). Exact tolerance (0) only;
the descent round-trip gate inside :func:`_lift_transform` keeps every landing sound. If no
self-lift lands, the loop warns and proceeds with whatever was achieved -- the upward search
may still collapse the supercell through a cross-group chain (the cubic
221(2a) -> 225(2a) -> 221(a) route works through the plain search, since cubic isomorphic
subgroups are never maximal).
"""
from httk.atomistic.symmetry.subgroups import isomorphic_subgroup_transforms
while True:
pseudo = _pseudo_translations(structure)
if not pseudo:
return structure
implied_index = len(pseudo) + 1
number = structure.spacegroup.it_number
transforms: list[SubgroupTransform] = []
seen: set[SubgroupTransform] = set()
for transform in subgroup_transforms(number, number) + isomorphic_subgroup_transforms(number):
if transform.index <= 1 or implied_index % transform.index != 0:
continue
if transform in seen:
continue
seen.add(transform)
transforms.append(transform)
transforms.sort(key=lambda transform: transform.index) # stable: table order within an index
landed: ASUStructure | None = None
for transform in transforms:
results = _lift_transform(structure, transform, 0.0)
if results:
landed = min(results, key=_canonical_result_key).asu
break
if landed is None:
logging.getLogger(__name__).warning(
"declared supercell in %s (implied index %d) not collapsed by tabulated isomorphic "
"self-lifts (indices up to 9); proceeding with the supercell description",
structure.spacegroup.setting,
implied_index,
)
return structure
structure = landed
def _niggli_reduced_entry(structure: ASUStructure) -> ASUStructure:
"""Re-express a triclinic-entry (SG 1 or 2) structure in its exact Niggli-reduced basis.
Different bases of the same lattice reduce to the same Niggli cell, so this removes the
basis-choice freedom before the search. A non-rational (surd) Gram matrix that
``niggli_reduce`` cannot handle degrades gracefully to the unreduced structure.
"""
# Local import: reduction is a sibling top-level module; importing it at module load would run
# before httk.atomistic finishes wiring symmetry, so keep it lazy.
from httk.atomistic.reduction import niggli_reduce
try:
reduced = niggli_reduce(structure.cell)
except (ValueError, ArithmeticError):
return structure
if reduced.transform == FracVector.eye((3, 3)):
return structure
rebased = _apply_normalizer_operation(structure, AffineOperation(reduced.transform.T().inv(), (0, 0, 0)))
return rebased if rebased is not None else structure
def _canonical_entry(structure: ASUStructure) -> ASUStructure:
"""Return the normal-form state ``highest_symmetry`` starts its search from.
Standardizes the input, collapses an exact supercell (a P1 one to its primitive cell; a
declared-symmetry one via tabulated isomorphic self-lifts), Niggli-reduces a triclinic
(SG 1/2) cell, and reduces to the normalizer-canonical normal form.
"""
current = _standard_input(structure)
if current.spacegroup.it_number == 1:
current = _primitive_reduced_entry(current)
else:
current = _isomorphic_reduced_entry(current)
if current.spacegroup.it_number in (1, 2):
current = _niggli_reduced_entry(current)
return _normal_form(current)
def _terminal_entry(structure: ASUStructure) -> ASUStructure:
"""Return the no-BFS entry reduced by the terminal-only joint normal form.
This shares standardization, exact supercell collapse, and triclinic Niggli reduction with
:func:`_canonical_entry`, but deliberately does not pass through the search-state normal form.
A compatible-system state quotient can select a setting that is useful while lifting yet hides
an alternative terminal metric minimum; direct no-BFS canonicalization must compare the full
terminal orbit from the pre-state-normal-form entry.
"""
current = _standard_input(structure)
if current.spacegroup.it_number == 1:
current = _primitive_reduced_entry(current)
else:
current = _isomorphic_reduced_entry(current)
if current.spacegroup.it_number in (1, 2):
current = _niggli_reduced_entry(current)
return _terminal_normal_form(current)
def _canonical_without_bfs(
structure: ASUStructure,
*,
preserve_chirality: bool = True,
) -> ASUStructure:
"""Return the canonical representative of ``structure`` within its own group -- no upward search.
It applies the same terminal normal form ``highest_symmetry`` uses at a no-lift terminal
(:func:`_terminal_normal_form` placed in the standard orientation of its metric), but reaches it
directly from the terminal-only entry (:func:`_terminal_entry`). Unlike the search, which reduces
each state through :func:`_normal_form` first, this does not pass through the search-state
quotient -- whose setting choice can hide an alternative terminal metric minimum -- so it compares
the full terminal orbit from the pre-state-normal-form entry. It is the deterministic, fully
invariant representation of the *recognized* symmetry, without hunting for pseudosymmetry above it.
``preserve_chirality`` defaults to ``True``: the recognized group is kept. When it is ``False``
and the entry lands in the higher member of an enantiomorphic pair, an exact chirality flip
(:func:`_enantiomorph`) normalizes it to the lower-numbered partner, re-taking the partner group's
normal form before the orientation step; a moment-carrying structure is left in its own group.
"""
entry = _canonical_orientation(_terminal_entry(structure))
# The normal-form pipeline strips moments before the flip would see them, so gate on the ORIGINAL
# input's moments -- a magnetic structure is left in its own group (see :func:`_enantiomorph`).
if not preserve_chirality and not any(site.moment is not None for site in structure.wyckoff_sites):
return normalize_chirality(entry)
return entry
[docs]
def normalize_chirality(structure: ASUStructure) -> ASUStructure:
"""Collapse a canonical ASU's higher enantiomorph onto its lower-numbered partner.
Structures outside the higher member of an enantiomorphic pair, including magnetic structures,
are returned by identity. A flipped partner is put back into its canonical normal form.
"""
if not isinstance(structure, ASUStructure):
raise TypeError(f"expected ASUStructure, got {type(structure).__name__}")
flipped = _enantiomorph(structure)
if flipped is None:
return structure
return _canonical_orientation(_terminal_normal_form(flipped))
[docs]
def highest_symmetry(
structure: ASUStructure,
*,
tolerance: float | None = None,
all_paths: bool = False,
preserve_chirality: bool = True,
) -> tuple[LiftResult, ...]:
"""Return all terminal upward lifts reached by breadth-first search.
:param structure: The starting asymmetric-unit structure.
:param tolerance: Cartesian acceptance tolerance, or the recognition-derived default.
:param all_paths: When ``False`` (default) the visited set collapses alternate Bärnighausen
routes to one entry per state, so each terminal appears once. When ``True`` the visited set
also keys on the accumulated path, so every distinct ``(terminal, path)`` pair is returned;
the ``.asu`` representatives of the extra results are identical, only ``path`` differs. The
state cap therefore binds sooner under the flag.
:param preserve_chirality: When ``True`` (the default), the recognized group is kept. When
``False``, a terminal landing in the higher member of an enantiomorphic pair is flipped by an
exact chirality transformation to its lower-numbered partner and re-normal-formed in that
partner group. The flip is applied only at terminal emission, never to mid-search states (the
Bärnighausen tables are per-group).
:return: Deterministically ordered highest-symmetry representations.
:raises ValueError: If the input is unsupported, or if the breadth-first search exceeds its
visited-state cap. A per-parent modular-solver branch-cap failure is not raised: that parent
target is skipped and reported through the ``"symmetry"`` warning channel, so in that rare
case the returned symmetry may be lower than the true maximum. The noisy least-squares
fallback is capped the same way -- a candidate whose integer-wrap box exceeds the noisy cap is
skipped (also on the ``"symmetry"`` warning channel), which could likewise lower the returned
symmetry for a genuinely noisy large-box lift.
An exact P1 supercell entry is first collapsed to its unique primitive description (any
multiplicity or sublattice orientation), and a triclinic (SG 1 or 2) entry is then Niggli-reduced
so the search is independent of the input basis choice. Each search state is then reduced to its
normalizer-canonical normal form --
special-site demotion plus the continuous- and discrete-Euclidean-normalizer translation
quotients and the affine-normalizer coset quotient -- collapsing origin-, basis- and
normalizer-equivalent representations to one visited entry so the search terminates from a raw
P1 input. The returned ``asu`` is that normalizer-canonical representative, its cell placed in
the standard orientation of its metric, so the result is invariant under origin shift, cell-basis
choice, and site order for the same crystal. ``path`` records the tabulated hops of the route
that reached it, and a bounded normalizer retry along that route applies tabulated cosets to
child coordinates and maps results back with the exact inverse, in tabulated coset order.
"""
if not isinstance(structure, ASUStructure):
raise TypeError(f"expected ASUStructure, got {type(structure).__name__}")
require_full_periodicity(structure.cell, "highest_symmetry")
if any(site.moment is not None for site in structure.wyckoff_sites):
raise ValueError("highest_symmetry does not support structures with site moments")
if structure.assemblies is not None:
raise ValueError("highest_symmetry does not support structures with assemblies")
if structure.molecular:
raise ValueError("highest_symmetry does not support molecular structures")
# Standardize, collapse exact P1 supercells, Niggli-reduce triclinic cells, and take the normal
# form -- the state the breadth-first search starts from.
current = _canonical_entry(structure)
accepted_tolerance = structure_tolerance(current) if tolerance is None else float(tolerance)
queue: list[tuple[ASUStructure, tuple[SubgroupTransform, ...], FracVector, Fraction]] = [
(current, (), FracVector((0, 0, 0)), Fraction(0))
]
start_signature = (current.spacegroup.it_number, _structure_signature(current))
visited: set[tuple[Any, ...]] = {(*start_signature, ()) if all_paths else start_signature}
terminals: list[LiftResult] = []
while queue:
state, path, shift, residual = queue.pop(0)
lifts = _highest_lifts(state, accepted_tolerance)
if not lifts:
emitted = _canonical_orientation(_terminal_normal_form(state))
if not preserve_chirality:
emitted = normalize_chirality(emitted)
terminals.append(_terminal_result(emitted, path, shift, residual))
continue
for result in lifts:
next_state = _normal_form(result.asu)
next_path = path + result.path
signature = (next_state.spacegroup.it_number, _structure_signature(next_state))
# all_paths keys the visited set on the route too, so alternate routes to a state are all
# explored rather than collapsed to the first one reached.
key = (*signature, next_path) if all_paths else signature
if key in visited:
continue
if len(visited) >= 10_000:
# ponytail: bounded state scan; add quotient-graph memoization if tables grow materially.
raise ValueError(f"highest_symmetry search cap exceeded for {current.spacegroup.setting}")
visited.add(key)
queue.append((next_state, next_path, result.shift, max(residual, result.residual)))
deduplicated: dict[tuple[Any, ...], LiftResult] = {}
for result in terminals:
deduplicated.setdefault((result.spacegroup.it_number, _structure_signature(result.asu), result.path), result)
return tuple(
sorted(
deduplicated.values(),
key=lambda result: (
-len(result.spacegroup.symmetry_operations),
result.spacegroup.it_number,
result.residual,
_canonical_result_key(result),
),
)
)
[docs]
def canonicalize(
structure: ASUStructure, *, tolerance: float | None = None, preserve_chirality: bool = True
) -> LiftResult:
"""Return the first deterministic highest-symmetry representation.
The result is the normalizer-canonical representative of the input's crystal: the same exact
``(it_number, sorted (species, wyckoff, free_params), cell basis)`` for any origin shift,
cell-basis choice (relabeling/shear), or site ordering of that crystal, and coherent with direct
entry at its own space group. See :func:`highest_symmetry` for the full contract; for a
P1/unit-cell start build the ASU in SG 1 and pass it here.
:param structure: The structure to canonicalize.
:param tolerance: Cartesian acceptance tolerance, or the recognition-derived default.
:param preserve_chirality: When ``True`` (the default), the recognized group is kept, so a
genuinely chiral crystal retains its handedness. When ``False``, a result in the higher
member of an enantiomorphic pair is normalized to its lower-numbered partner by an exact
chirality flip (magnetic structures excepted; :func:`~httk.atomistic.normalize_chirality` is the standalone
bridge). See :func:`highest_symmetry`.
:return: The canonical terminal lift.
"""
return highest_symmetry(structure, tolerance=tolerance, preserve_chirality=preserve_chirality)[0]
def _supergroup_path(start: int, target: int) -> tuple[int, ...] | None:
queue: list[tuple[int, tuple[int, ...]]] = [(start, ())]
visited = {start}
while queue:
current, path = queue.pop(0)
for parent in minimal_supergroups(current):
if parent == target:
return path + (parent,)
if parent not in visited:
visited.add(parent)
queue.append((parent, path + (parent,)))
return None
[docs]
def rerepresent(
structure: ASUStructure,
target: Spacegroup | int,
*,
tolerance: float | None = None,
) -> ASUStructure:
"""Express a structure in a reachable subgroup or supergroup setting.
:param structure: The input asymmetric-unit structure.
:param target: The target space group or IT number.
:param tolerance: Cartesian acceptance tolerance for upward lifts.
:return: The target-group asymmetric unit.
:raises ValueError: If the target is unrelated, an upward hop has no lift, or a cross-group
rerepresentation requires descending or lifting a structure with site moments, assemblies,
or molecular semantics.
"""
if not isinstance(structure, ASUStructure):
raise TypeError(f"expected ASUStructure, got {type(structure).__name__}")
target_group = (target if isinstance(target, Spacegroup) else Spacegroup.standard(target)).standard_setting()
structure = _standard_input(structure)
current_number = structure.spacegroup.it_number
if target_group.it_number == current_number:
return structure
if target_group.it_number in subgroup_closure(current_number):
return subgroup_representation(structure, target_group).asu
if target_group.it_number in supergroup_closure(current_number):
path = _supergroup_path(current_number, target_group.it_number)
if path is None:
raise ValueError(f"no supergroup path from {structure.spacegroup.setting} to {target_group.setting}")
current = structure
for parent_number in path:
results = backward_lift(current, parent_number, tolerance=tolerance)
if not results:
raise ValueError(
f"no lift from {current.spacegroup.setting} to {Spacegroup.standard(parent_number).setting}"
)
current = min(results, key=lambda result: (result.residual, _canonical_result_key(result))).asu
return current
raise ValueError(f"space groups {structure.spacegroup.setting} and {target_group.setting} are unrelated")