"""Site-permutation symmetry operations on the magnetic sublattice."""
from __future__ import annotations
import logging
from typing import List, Optional, Tuple
import numpy as np
log = logging.getLogger(__name__)
def _spglib_rotations(lat_vecs: np.ndarray,
frac: np.ndarray,
numbers: List[int],
symprec: float) -> Tuple[np.ndarray, np.ndarray]:
"""(rotations, translations) of the cell from spglib, fractional basis.
Falls back to the identity operation (with a warning) if spglib cannot
analyse the cell — the caller then still gets time-reversal-only grouping
rather than a crash.
"""
import spglib # noqa: WPS433 — runtime dep via pymatgen
cell = (np.asarray(lat_vecs, dtype=float),
np.asarray(frac, dtype=float) % 1.0,
np.asarray(numbers, dtype=int))
sym = spglib.get_symmetry(cell, symprec=symprec)
if sym is None:
log.warning("spglib could not analyse the cell — using identity only "
"(time-reversal-only config grouping).")
return (np.eye(3, dtype=int)[None, :, :], np.zeros((1, 3)))
return sym["rotations"], sym["translations"]
[docs]
def find_symmetry_ops(cart_coords: np.ndarray,
lat_vecs: np.ndarray,
tol: float = 0.05,
*,
crystal=None) -> List[Tuple[int, ...]]:
"""Find site-permutation symmetry operations of the magnetic sublattice.
Candidate operations (rotation + translation pairs) come from **spglib**,
so any axis orientation and any space group is handled — the result no
longer depends on which lattice vector is the stacking axis. Each
candidate is then verified site-by-site (tolerance ``tol`` Å) and returned
as a permutation tuple; ``perm[i]`` is the index of the site to which site
``i`` is mapped.
Parameters
----------
cart_coords : Cartesian coordinates of the magnetic sites.
lat_vecs : 3×3 lattice matrix, **rows** = a, b, c.
tol : site-matching tolerance in Å (also spglib's ``symprec``).
crystal : optional :class:`mag4.constants.CrystalData` of the **full**
structure. When given, symmetry is computed on the full decorated
crystal (magnetic + non-magnetic atoms, distinguished by element), so
every returned operation is a true crystal symmetry — configurations
it merges are guaranteed degenerate. When omitted, only the bare
magnetic sublattice is analysed (its symmetry can be higher than the
crystal's; combine with
:func:`mag4.geometry.filter_symops_by_fingerprint` if ligand
environments matter).
"""
cart_coords = np.asarray(cart_coords, dtype=float)
lat_vecs = np.asarray(lat_vecs, dtype=float)
n = len(cart_coords)
frac_mag = cart_coords @ np.linalg.inv(lat_vecs)
if crystal is not None:
# Number atoms by ELEMENT (not by label group) so crystallographically
# equivalent sites that the input file lists under different labels
# (e.g. Cu001 / Cu002) may still map onto each other.
el_num = {el: i + 1 for i, el in enumerate(dict.fromkeys(crystal.elements))}
frac_all: List[List[float]] = []
numbers: List[int] = []
for grp_idx, grp in enumerate(crystal.all_species):
num = el_num[crystal.elements[grp_idx]]
for a in grp:
frac_all.append([float(a[1]), float(a[2]), float(a[3])])
numbers.append(num)
rots, trans = _spglib_rotations(lat_vecs, np.asarray(frac_all),
numbers, tol)
else:
rots, trans = _spglib_rotations(lat_vecs, frac_mag, [1] * n, tol)
ops: List[Tuple[int, ...]] = []
seen = set()
for W, t in zip(rots, trans):
new_frac = frac_mag @ np.asarray(W, dtype=float).T + t
# Vectorised site matching: minimum-image distance of every mapped
# site to every original site in ONE (n,n,3) operation. The previous
# per-site Python loop (find_site × n) made high-symmetry supercells
# pathological — a cubic 2×2×2 cell yields 1536 candidate ops, all
# valid (so no early exit), at O(n²) small numpy calls each ≈ tens of
# minutes; this form does the identical math in milliseconds.
d = new_frac[:, None, :] - frac_mag[None, :, :]
d -= np.round(d) # minimum image, fractional
dist = np.linalg.norm(d @ lat_vecs, axis=2)
perm = np.argmin(dist, axis=1)
if not bool(np.all(dist[np.arange(n), perm] < tol)):
continue
if len(set(perm.tolist())) != n:
continue
p = tuple(int(j) for j in perm)
if p not in seen:
seen.add(p)
ops.append(p)
log.info("Found %d symmetry operations (%s, via spglib)",
len(ops),
"full crystal" if crystal is not None else "magnetic sublattice")
return ops