"""Exact rational Wyckoff-position algebra.
A Wyckoff position is a family of symmetry-equivalent sites parameterised by a few free
coordinates: SG 15 letter ``e`` is ``0,y,1/4``, so one free parameter ``y`` places four
symmetry-equivalent atoms. This module goes both ways between free parameters and
coordinates, entirely over the rationals.
Two properties of the vendored tables make this much simpler than it looks, both asserted
in ``tests/test_symmetry_data.py``:
* Each Wyckoff position's ``orbit`` is already the **complete, deduplicated** list of
affine maps, of length exactly ``multiplicity``, with centering translations folded in.
Generating an orbit is therefore a plain loop with no coincidence testing and no
tolerance.
* ``hasfreedom`` marks which of ``x``, ``y``, ``z`` are free, and the columns of every
orbit matrix for the non-free variables are identically zero, with
``sum(hasfreedom) == rank``. So the free parameters are read straight off, and the
``first_orbit`` strings (``"1/8,y,-y+1/4"``) never need parsing — the same information
is already present as an exact affine map.
Everything is expressed in the coordinates of whichever setting the record came from.
"""
import fractions
from collections.abc import Mapping, Sequence
from functools import cached_property
from typing import Any, Self, cast
from httk.core import FracVector
from httk.atomistic.symmetry.affine_operation import AffineOperation
__all__ = ["WyckoffBranch", "WyckoffPosition", "wyckoff_positions"]
[docs]
class WyckoffBranch:
"""Represent one member of a Wyckoff orbit as an affine parameter map.
A position of multiplicity *m* has *m* branches. The representative branch is the
first, but a coordinate may lie on any of them, which is why
:meth:`WyckoffPosition.parameters_of` tries them all.
:param operation: The affine map for this orbit branch.
:param free: The indices of the free parameters in ``(x, y, z)``.
"""
_operation: AffineOperation
_free: tuple[int, ...]
_unimodular: FracVector
_echelon: FracVector
_pivot_inverse: FracVector
def __init__(self, operation: AffineOperation, free: Sequence[int]) -> None:
self._operation = operation
self._free = tuple(free)
self._unimodular, self._echelon = _row_hermite(operation.matrix, self._free)
rank = len(self._free)
if rank:
pivot = [[self._echelon[row][column].to_fraction() for column in range(rank)] for row in range(rank)]
self._pivot_inverse = FracVector(_invert_small(pivot))
else:
self._pivot_inverse = FracVector(())
@property
[docs]
def operation(self) -> AffineOperation:
"""Return the affine map from ``(x, y, z)`` to this branch's coordinate.
:return: The branch's affine operation.
"""
return self._operation
@property
[docs]
def free(self) -> tuple[int, ...]:
"""Return the ascending indices of the free parameters in ``(x, y, z)``.
:return: The free-parameter indices.
"""
return self._free
[docs]
def coordinate(self, parameters: Any) -> FracVector:
"""Compute the coordinate this branch places at given free-parameter values.
``parameters`` has one entry per free parameter. Because the non-free columns of
the matrix are zero, the values placed at the non-free positions are irrelevant.
:param parameters: The free-parameter values for this branch.
:return: The exact reduced coordinate generated by the branch.
:raises ValueError: If the number of parameters does not match the branch.
"""
placed = [fractions.Fraction(0)] * 3
values = FracVector(parameters)
if not self._free:
# A fixed position has no parameters. FracVector has no zero-length shape, so
# "no parameters" is the empty vector, whose dim reads as () or (0,).
if values.dim not in ((), (0,)):
raise ValueError(f"a fixed Wyckoff position takes no free parameters, got dim {values.dim}")
return self._operation.apply(placed)
if values.dim != (len(self._free),):
raise ValueError(f"expected {len(self._free)} free parameter(s), got dim {values.dim}")
for slot, index in enumerate(self._free):
placed[index] = values[slot].to_fraction()
return self._operation.apply(placed)
[docs]
def parameters_of(self, coordinate: Any) -> FracVector | None:
"""Recover the free parameters putting this branch on ``coordinate``, if possible.
``None`` means the coordinate does not lie on this branch — for *any* lattice
translation, not merely the one given. That completeness is what the row-Hermite
form buys: with ``U`` unimodular over the integers, ``A t = d (mod Z^3)`` holds iff
``U A t = U d (mod Z^3)``, and the zero rows of ``U A`` turn the lattice-membership
question into "are these components integers?" with no search over translations.
The returned parameters are reduced into ``[0, 1)``, which is canonical: the pivot
block has determinant ``±1`` throughout the vendored tables, so the solution is
unique modulo one. The result is verified by re-evaluating the branch, so a table
that ever violated that assumption would yield a clean miss rather than a wrong
answer.
:param coordinate: The exact reduced coordinate to match.
:return: The normalized free parameters, or ``None`` when the coordinate is not on
this branch.
"""
point = FracVector(coordinate)
rank = len(self._free)
mapped = (point - self._operation.vector) * self._unimodular.T()
# Rows below the rank carry no free parameter, so they must already be lattice
# translations for the coordinate to lie on this branch at all. Checking that first
# is a cheap rejection for the common case of a coordinate that is simply elsewhere.
for row in range(rank, 3):
if mapped[row].to_fraction().denominator != 1:
return None
parameters = self._parameters_from_mapped(mapped)
if self.coordinate(parameters).normalize() != point.normalize():
return None
return parameters
[docs]
def nearest_parameters(self, coordinate: Any) -> FracVector:
"""Compute free parameters that put this branch as close to ``coordinate`` as possible.
Unlike :meth:`parameters_of` this always returns a value: the free directions are
solved exactly and any discrepancy is left in the *fixed* directions, where the
branch's own coordinates win. It is the projection used when recognizing a measured
structure, whose coordinates carry rounding and do not lie exactly on any position.
The projection is taken along the branch's own lattice basis rather than being
minimised in the cell metric, so for a strongly oblique cell it is a near-optimal
rather than provably optimal choice. That is safe because the caller measures the
resulting Cartesian displacement and rejects anything beyond its tolerance — the
method can cost a match, never grant a wrong one.
:param coordinate: The reduced coordinate to approximate.
:return: The normalized free parameters for the nearest branch point.
"""
point = FracVector(coordinate)
mapped = (point - self._operation.vector) * self._unimodular.T()
return self._parameters_from_mapped(mapped)
[docs]
def nearest_parameters_float(self, coordinate: Sequence[float]) -> tuple[float, ...]:
"""Project a coordinate onto this branch in floating point for screening.
This deliberately mirrors :meth:`nearest_parameters` without its exact verification.
It is only a candidate-screening aid; callers must calculate and compare the final
distance with the exact methods before accepting a match.
:param coordinate: A reduced coordinate as three floating-point values.
:return: The normalized floating-point free parameters.
"""
rank = len(self._free)
if not rank:
return ()
vector, unimodular, inverse, _ = self._float_coefficients
mapped = [
sum((coordinate[column] - vector[column]) * unimodular[row][column] for column in range(3))
for row in range(rank)
]
return tuple(sum(mapped[row] * inverse[column][row] for row in range(rank)) % 1.0 for column in range(rank))
[docs]
def coordinate_float(self, parameters: Sequence[float]) -> tuple[float, float, float]:
"""Evaluate this branch in floating point for candidate screening.
:param parameters: One floating-point value for every free parameter.
:return: The unwrapped floating-point reduced coordinate.
:raises ValueError: If the parameter count is wrong.
"""
if len(parameters) != len(self._free):
raise ValueError(f"expected {len(self._free)} free parameter(s), got {len(parameters)}")
placed = [0.0, 0.0, 0.0]
for slot, index in enumerate(self._free):
placed[index] = parameters[slot]
vector, _, _, matrix = self._float_coefficients
return cast(
tuple[float, float, float],
tuple(sum(matrix[row][column] * placed[column] for column in range(3)) + vector[row] for row in range(3)),
)
@cached_property
def _float_coefficients(
self,
) -> tuple[
tuple[float, ...],
tuple[tuple[float, ...], ...],
tuple[tuple[float, ...], ...],
tuple[tuple[float, ...], ...],
]:
"""Cache immutable coefficients shared by every floating-point screen."""
return (
tuple(self._operation.vector.to_floats()),
tuple(tuple(row) for row in self._unimodular.to_floats()),
tuple(tuple(row) for row in self._pivot_inverse.to_floats()),
tuple(tuple(row) for row in self._operation.matrix.to_floats()),
)
def _parameters_from_mapped(self, mapped: FracVector) -> FracVector:
"""Back-substitute the free parameters out of ``U * (point - b)``."""
rank = len(self._free)
if rank == 0:
return FracVector(())
head = FracVector([mapped[row].to_fraction() for row in range(rank)])
return (head * self._pivot_inverse.T()).normalize()
[docs]
class WyckoffPosition:
"""Represent a Wyckoff position of one space-group setting.
:param record: The vendored record describing the Wyckoff position.
"""
_letter: str
_multiplicity: int
_site_symmetry: str
_free: tuple[int, ...]
_branches: tuple[WyckoffBranch, ...]
def __init__(self, record: Mapping[str, Any]) -> None:
self._letter = record["letter"]
self._multiplicity = record["multiplicity"]
self._site_symmetry = record["sitesym"]
self._free = tuple(index for index, free in enumerate(record["hasfreedom"]) if free)
self._branches = tuple(WyckoffBranch(AffineOperation.from_record(item), self._free) for item in record["orbit"])
if len(self._branches) != self._multiplicity:
raise ValueError(
f"Wyckoff position {self._letter!r} has {len(self._branches)} orbit entries "
f"but multiplicity {self._multiplicity}"
)
@classmethod
[docs]
def from_record(cls, record: Mapping[str, Any]) -> Self:
"""Build a Wyckoff position from a vendored record.
:param record: The vendored Wyckoff-position record.
:return: The corresponding Wyckoff position.
"""
return cls(record)
@property
[docs]
def letter(self) -> str:
"""Return the bare Wyckoff letter, such as ``"e"``.
:return: The Wyckoff letter without a multiplicity prefix.
"""
return self._letter
@property
[docs]
def multiplicity(self) -> int:
"""Return the number of sites generated by one parameter set.
:return: The position multiplicity in the unit cell.
"""
return self._multiplicity
@property
[docs]
def site_symmetry(self) -> str:
"""Return the site-symmetry group in Hermann-Mauguin notation.
:return: The site-symmetry symbol.
"""
return self._site_symmetry
@property
[docs]
def free(self) -> tuple[int, ...]:
"""Return the indices of the free parameters in ``(x, y, z)``.
:return: The free-parameter indices.
"""
return self._free
@property
[docs]
def free_count(self) -> int:
"""Return the number of degrees of freedom of the position.
:return: The number of free parameters, from zero through three.
"""
return len(self._free)
@property
[docs]
def branches(self) -> tuple[WyckoffBranch, ...]:
"""Return the complete, deduplicated orbit branches.
:return: One branch for each equivalent site.
"""
return self._branches
@property
[docs]
def representative(self) -> WyckoffBranch:
"""Return the first orbit member printed as ``first_orbit`` in the tables.
:return: The representative orbit branch.
"""
return self._branches[0]
[docs]
def coordinates(self, parameters: Any) -> FracVector:
"""Compute every coordinate of the orbit as an exact ``(multiplicity, 3)`` block.
Not wrapped and not deduplicated: within one setting the tabulated orbit is
already distinct, so wrapping is the caller's business (and matters only once a
setting transform enters).
:param parameters: The free-parameter values for the position.
:return: The unwrapped coordinates of all orbit branches.
:raises ValueError: If the number of parameters does not match a branch.
"""
return FracVector([branch.coordinate(parameters).to_fractions() for branch in self._branches])
[docs]
def parameters_of(self, coordinate: Any) -> FracVector | None:
"""Recover free parameters placing some branch on ``coordinate``, if possible.
Tries every branch, not only the representative. That matters: across the vendored
tables, 11673 of the 20639 non-representative orbit members lie on a different
branch than the representative, so a matcher that only tested ``first_orbit``
would reject a majority of legitimate orbit points.
:param coordinate: The exact reduced coordinate to match.
:return: The normalized free parameters, or ``None`` when no branch matches.
"""
for branch in self._branches:
parameters = branch.parameters_of(coordinate)
if parameters is not None:
return parameters
return None
def __repr__(self) -> str:
return f"WyckoffPosition({self._multiplicity}{self._letter}, sitesym={self._site_symmetry!r})"
_WYCKOFF_CACHE: dict[str, tuple[WyckoffPosition, ...]] = {}
[docs]
def wyckoff_positions(record: Mapping[str, Any]) -> tuple[WyckoffPosition, ...]:
"""Build the Wyckoff positions of a setting record, most specific first.
Ordered by ``(free_count, multiplicity, letter)`` so that the first match found when
identifying a coordinate is the most specific position it lies on. Ties do not arise:
positions are affine subspaces, so a coordinate on two distinct positions of the same
dimension also lies on their lower-dimensional intersection, which is covered by an
earlier entry.
:param record: The vendored space-group setting record.
:return: The setting's ordered Wyckoff positions.
"""
key = record["hall_entry"]
cached = _WYCKOFF_CACHE.get(key)
if cached is not None:
return cached
positions = [WyckoffPosition(entry) for entry in record["wyckoff"]]
positions.sort(key=lambda position: (position.free_count, position.multiplicity, position.letter))
return _WYCKOFF_CACHE.setdefault(key, tuple(positions))
def _invert_small(rows: list[list[fractions.Fraction]]) -> list[list[fractions.Fraction]]:
"""The exact inverse of a small square rational matrix, by Gauss-Jordan elimination.
:meth:`~httk.core.FracVector.inv` covers only scalars and 3x3 matrices, but the pivot
blocks here are 1x1 and 2x2 as well (a Wyckoff position with one or two degrees of
freedom). The sizes are at most 3, so the plain algorithm is entirely adequate.
"""
size = len(rows)
work = [list(row) + [fractions.Fraction(1 if i == j else 0) for j in range(size)] for i, row in enumerate(rows)]
for column in range(size):
pivot = next((r for r in range(column, size) if work[r][column] != 0), None)
if pivot is None:
raise ValueError("singular matrix: a Wyckoff pivot block must be invertible")
work[column], work[pivot] = work[pivot], work[column]
scale = work[column][column]
work[column] = [value / scale for value in work[column]]
for row in range(size):
if row != column and work[row][column] != 0:
factor = work[row][column]
work[row] = [a - factor * b for a, b in zip(work[row], work[column])]
return [row[size:] for row in work]
def _row_hermite(matrix: FracVector, free: Sequence[int]) -> tuple[FracVector, FracVector]:
"""Row-reduce the free columns of an orbit matrix to echelon form over the integers.
Returns ``(U, H)`` with ``U`` unimodular (integer, determinant ``±1``) and
``H = U * A`` in row echelon form, where ``A`` is ``matrix`` restricted to the free
columns. Unimodularity is the point: it preserves the integer lattice, so a
congruence modulo ``Z^3`` may be applied to ``H`` and ``U d`` instead of ``A`` and
``d`` without changing which coordinates satisfy it.
Uses gcd-based elimination (integer row reduction), so no division is ever inexact.
Orbit matrices are integer with entries in ``{-2, ..., 2}``, so the values stay tiny.
"""
rows = [[int(matrix[row][column].to_fraction()) for column in free] for row in range(3)]
unimodular = [[1 if row == column else 0 for column in range(3)] for row in range(3)]
pivot_row = 0
for column in range(len(free)):
# Drive all entries below the pivot to zero by repeated gcd steps.
for row in range(pivot_row + 1, 3):
while rows[row][column] != 0:
if abs(rows[pivot_row][column]) > abs(rows[row][column]) or rows[pivot_row][column] == 0:
rows[pivot_row], rows[row] = rows[row], rows[pivot_row]
unimodular[pivot_row], unimodular[row] = unimodular[row], unimodular[pivot_row]
if rows[pivot_row][column] == 0:
break
factor = rows[row][column] // rows[pivot_row][column]
rows[row] = [a - factor * b for a, b in zip(rows[row], rows[pivot_row])]
unimodular[row] = [a - factor * b for a, b in zip(unimodular[row], unimodular[pivot_row])]
if rows[pivot_row][column] != 0:
pivot_row += 1
if pivot_row == 3:
break
return FracVector(unimodular), FracVector(rows)