"""Geometric fingerprints for J-shell discrimination via M–X–M bond angles.
The four-state method needs each magnetic coupling to correspond to a unique
physical exchange path. Grouping pairs by **distance alone** lumps together
pairs that share the same M–M distance but route their super-exchange through
**different** bridging ligands at **different** M–L–M angles — physically
distinct couplings the four-state method would conflate into a single J.
This module adds a *secondary discriminator* derived from the bridging-ligand
geometry of each pair: a quantised, sorted tuple of M–L–M angles. Pairs are
clustered first by distance, then sub-clustered by fingerprint. Resulting
sub-shells inherit their parent's index plus a letter suffix (``J1a``,
``J1b``, …).
The fingerprint computation reuses
:func:`mag4.analysis.bond_angles.find_bridging_angles`; no new physics is
introduced.
"""
from __future__ import annotations
import logging
from typing import Dict, List, Sequence, Tuple
import numpy as np
from mag4.analysis.bond_angles import find_bridging_angles, find_dihedral_paths
from mag4.constants import CrystalData, DIST_EQUIV_TOL
log = logging.getLogger(__name__)
#: A geometric fingerprint for an M–M pair, packed as a 4-tuple:
#:
#: ``(n_bridges, sorted bridge-angle bins,
#: n_dihedrals, sorted signed dihedral bins)``
#:
#: * ``n_bridges`` / ``bridge bins`` describe the direct M–L–M paths.
#: * ``n_dihedrals`` / ``dihedral bins`` describe the M–L–L–M
#: super-super-exchange paths through two **distinct** ligands joined
#: by an L–L contact within ``ll_cutoff``.
#:
#: The dihedral component is populated **only when** ``n_bridges == 0`` and
#: ``ll_cutoff > 0`` — direct-bridge pairs keep ``(0, ())`` for the
#: dihedral slot so they are not over-split by what is, for them, a
#: redundant secondary signal.
AngleFingerprint = Tuple[int, Tuple[int, ...], int, Tuple[int, ...]]
[docs]
def quantise_angles(angles: Sequence[float], angle_tol: float) -> Tuple[int, ...]:
"""Bin a sequence of angles to integer multiples of ``angle_tol``.
Returns a sorted tuple so equal fingerprints compare equal regardless of
the input order. Empty input yields the empty tuple.
"""
if not angles:
return ()
if angle_tol <= 0:
raise ValueError(f"angle_tol must be > 0, got {angle_tol}")
return tuple(sorted(int(round(a / angle_tol)) for a in angles))
[docs]
def compute_pair_fingerprint(
structure,
site_i: int,
site_j: int,
image: Tuple[int, int, int],
ligand_elements: Sequence[str],
bond_cutoff: float,
angle_tol: float,
*,
ll_cutoff: float = 0.0,
signed_dihedral: bool = False,
) -> AngleFingerprint:
"""Geometric fingerprint of one M–M pair.
Parameters
----------
structure : pymatgen.core.Structure
The cell containing the pair (primitive or supercell — both work).
site_i, site_j : int
Indices of the two magnetic atoms in ``structure``.
image : tuple of 3 int
Periodic image of ``site_j`` relative to ``site_i``. Use
``(0, 0, 0)`` for direct intra-supercell pairs.
ligand_elements : sequence of str
Element symbols treated as bridging ligands (e.g. ``["O"]``).
bond_cutoff : float
Maximum M–L distance (Å) for a ligand to be considered bonded.
angle_tol : float
Tolerance (°) for quantising both M–L–M angles and M–L–L–M
dihedrals into integer bins.
ll_cutoff : float
Maximum L–L distance (Å) for a pair of distinct ligands to be
counted as a super-super-exchange bridge. Used **only** when the
pair has no direct M–L–M bridge (``n_bridges == 0``). Set to
``0`` to disable dihedral analysis.
signed_dihedral : bool
Keep the IUPAC sign of each M–L–L–M dihedral when ``True`` —
this distinguishes mirror-image (chiral) paths and is
appropriate when Dzyaloshinskii–Moriya / antisymmetric exchange
matters. Default ``False`` folds every dihedral to ``|θ|`` in
``[0°, 180°]`` so chiral counterparts are merged, which is the
right choice for plain Heisenberg J.
Returns
-------
AngleFingerprint
4-tuple ``(n_bridges, bridge bins, n_dihedrals, dihedral bins)``.
For pairs with a direct bridge the dihedral slot is
``(0, ())`` — selective fallback so bridged shells aren't
over-split.
"""
raw = find_bridging_angles(
structure, site_i, site_j, image, ligand_elements, bond_cutoff
)
bridge_angles = [a["angle"] for a in raw]
n_bridges = len(bridge_angles)
bridge_bins = quantise_angles(bridge_angles, angle_tol)
if n_bridges == 0 and ll_cutoff > 0:
dihedrals = find_dihedral_paths(
structure, site_i, site_j, image,
ligand_elements, bond_cutoff, ll_cutoff,
)
dihedral_angles = [p["dihedral"] for p in dihedrals]
if not signed_dihedral:
# Fold to [0, 180°]: chiral mirror-image paths merge.
dihedral_angles = [abs(d) for d in dihedral_angles]
n_dihedrals = len(dihedral_angles)
dihedral_bins = quantise_angles(dihedral_angles, angle_tol)
else:
n_dihedrals = 0
dihedral_bins: Tuple[int, ...] = ()
return (n_bridges, bridge_bins, n_dihedrals, dihedral_bins)
def _format_fingerprint(fp: AngleFingerprint, angle_tol: float,
*, signed_dihedral: bool = False) -> str:
"""Human-readable fingerprint summary.
* With a direct bridge: ``"2×96.0°"`` (M–L–M angles).
* No direct bridge but dihedrals (unsigned): ``"dih: 2×150.0°"``.
* No direct bridge but dihedrals (signed): ``"dih: 2×-150.0°"``.
* No bridge and no dihedral: ``"no bridge"``.
"""
from collections import Counter
n_bridges, bridge_bins, n_dihedrals, dihedral_bins = fp
if n_bridges > 0:
degs = [b * angle_tol for b in bridge_bins]
cnt = Counter(degs)
parts = []
for deg, k in sorted(cnt.items()):
parts.append(f"{deg:.1f}°" if k == 1 else f"{k}×{deg:.1f}°")
return ", ".join(parts)
if n_dihedrals > 0:
degs = [b * angle_tol for b in dihedral_bins]
cnt = Counter(degs)
parts = []
fmt = "{:+.1f}°" if signed_dihedral else "{:.1f}°"
for deg, k in sorted(cnt.items()):
txt = fmt.format(deg)
parts.append(txt if k == 1 else f"{k}×{txt}")
return "dih: " + ", ".join(parts)
return "no bridge"
[docs]
def group_inequivalent_geom(
results: list,
structure,
*,
ligand_elements: Sequence[str],
bond_cutoff: float = 2.5,
dist_tol: float = DIST_EQUIV_TOL,
angle_tol: float = 1.0,
ll_cutoff: float = 0.0,
signed_dihedral: bool = False,
) -> list:
"""Cluster M–M pairs by ``(distance, angle fingerprint)``.
First-level clustering by distance is identical to
:func:`mag4.io_cif.group_inequivalent`. Within each distance shell, pairs
are then sub-grouped by their bridging-ligand fingerprint. Sub-shells
take suffix letters (``J1a``, ``J1b``, …) when a split occurs; singletons
keep the plain ``J1`` label.
Each output dict carries an extra ``fingerprint`` field (``AngleFingerprint``)
and a ``fingerprint_str`` field (human-readable summary).
"""
if not results:
return []
# First-level: by distance.
by_dist: List[List[dict]] = []
current = [results[0]]
for pair in results[1:]:
if abs(pair["distance"] - current[-1]["distance"]) <= dist_tol:
current.append(pair)
else:
by_dist.append(current)
current = [pair]
by_dist.append(current)
# Annotate every pair with its fingerprint.
for shell in by_dist:
for pair in shell:
pair["fingerprint"] = compute_pair_fingerprint(
structure, pair["site_i"], pair["site_j"], pair["image"],
ligand_elements, bond_cutoff, angle_tol,
ll_cutoff=ll_cutoff,
signed_dihedral=signed_dihedral,
)
# Second-level: by fingerprint within each distance shell.
suffixes = "abcdefghijklmnopqrstuvwxyz"
grouped: list = []
shell_idx = 1
for shell in by_dist:
buckets: Dict[AngleFingerprint, List[dict]] = {}
for pair in shell:
buckets.setdefault(pair["fingerprint"], []).append(pair)
keys = sorted(buckets.keys()) # deterministic suffix order
for k, fp in enumerate(keys):
sub = buckets[fp]
dists = [p["distance"] for p in sub]
label = f"J{shell_idx}" if len(keys) == 1 else f"J{shell_idx}{suffixes[k]}"
grouped.append(dict(
label=label,
distance=sum(dists) / len(dists),
dist_min=min(dists),
dist_max=max(dists),
multiplicity=len(sub),
pairs=sub,
fingerprint=fp,
fingerprint_str=_format_fingerprint(
fp, angle_tol, signed_dihedral=signed_dihedral),
))
shell_idx += 1
return grouped
# ---------------------------------------------------------------------------
# Supercell ↔ pymatgen.Structure adapter
# ---------------------------------------------------------------------------
[docs]
def filter_symops_by_fingerprint(
sym_ops: Sequence[Tuple[int, ...]],
structure,
n_mag: int,
*,
ligand_elements: Sequence[str],
bond_cutoff: float,
angle_tol: float,
ll_cutoff: float = 0.0,
signed_dihedral: bool = False,
) -> List[Tuple[int, ...]]:
"""Drop magnetic-sublattice symmetry ops that disagree with the bridging
environment.
The default :func:`mag4.symmetry.find_symmetry_ops` considers **only**
the magnetic sites: two pairs related by such an op are considered
equivalent even if their bridging ligands sit very differently in space.
When that happens, :func:`mag4.spinconfig.find_unique_configs` merges
configurations from physically distinct sub-shells (e.g. ``J1a`` vs
``J1b``) and the extraction formulas degenerate.
This filter keeps an op only when the permutation maps every M–M pair
onto another pair with the **same** bridging fingerprint — i.e. the op
is a symmetry of the *decorated* magnetic structure, not just of the
bare sublattice. Fingerprints are computed once per pair and cached.
"""
from itertools import combinations
# Pre-compute fingerprints for every unordered pair of magnetic atoms.
fps: Dict[Tuple[int, int], AngleFingerprint] = {}
for i, j in combinations(range(n_mag), 2):
fps[(i, j)] = compute_pair_fingerprint(
structure, i, j, (0, 0, 0),
ligand_elements, bond_cutoff, angle_tol,
ll_cutoff=ll_cutoff,
signed_dihedral=signed_dihedral,
)
valid: List[Tuple[int, ...]] = []
for op in sym_ops:
keep = True
for (i, j), fp_ij in fps.items():
pi, pj = op[i], op[j]
key = (pi, pj) if pi < pj else (pj, pi)
if fps[key] != fp_ij:
keep = False
break
if keep:
valid.append(tuple(op))
return valid
[docs]
def reduce_to_wien2k_cell(
crystal: CrystalData,
group: List,
*,
spin_magnitude: float = 1.0,
symprec: float = 1e-3,
) -> CrystalData:
"""Reduce a :class:`CrystalData` to the **standardized conventional
cell** of its spin-decorated (colored) structure — the setting
WIEN2k's ``init_lapw`` pipeline expects.
Uses ``spglib.standardize_cell(no_idealize=False, to_primitive=False)``
on the same colored ``(lattice, frac_coords, atom_types)`` tuple
that :func:`analyze_config_symmetry` already builds.
* ``to_primitive=False`` — return the **conventional** cell, not
the primitive cell. WIEN2k's lattice-type letters (F, I, R, …)
are defined relative to the conventional cell of the parent
space group; for F-centered cubic groups like Fm-3m the
primitive cell is rhombohedral with 60° angles, which WIEN2k's
``sgroup`` immediately rewrites to the cubic conventional
setting — yielding a sym-op mismatch with whatever mag4 had
written. Emitting the conventional cell directly avoids that
round trip.
* ``no_idealize=False`` — apply spglib's **origin shift** so the
cell sits at the canonical origin of its space group (4/m m m
at the corner for P4/mmm, etc.). Without this, ``sgroup``
complains "Wrong multiplicity for atom 1" / "Struct file is
not consistent" because the orbit-grouping mag4 reports
(computed in spglib's canonical setting) does not match the
atom positions (kept in mag4's input setting). Idealizing the
cell aligns positions and orbits in one go.
For primitive groups (Pm-3m, P4/mmm, …) the conventional and
primitive cells coincide; only the origin shift is non-trivial.
Returns a brand-new ``CrystalData`` whose ``cell``,
``lattice_vectors``, ``all_species``, ``target_species``,
``nb_species``, ``total_atoms`` and ``elements`` describe the
reduced cell. Atom labels mirror the ``crystal.elements`` order
(``Ni1``, ``Ni2``, …) re-numbered inside the smaller cell.
Falls back to the input ``crystal`` when spglib refuses to reduce
*and* would not shift the origin either (typically: cell already
in canonical setting and primitive). The fallback path emits a
single debug log entry and never raises.
"""
try:
import spglib # noqa: WPS433 — runtime dep via pymatgen
except ImportError: # pragma: no cover - spglib ships with pymatgen
return crystal
# The helpers live in mag4.wien2k (they were factored there first).
from mag4.wien2k import (_elements_per_atom, _frac_coords_array,
_spglib_atom_types)
lat = np.asarray(crystal.lattice_vectors, dtype=float)
frac = _frac_coords_array(crystal)
types = _spglib_atom_types(crystal, group, spin_magnitude)
if len(types) != len(frac): # defensive
return crystal
# First inspect the symmetry dataset of the input cell to find out
# (a) whether spglib wants a true cell reduction (T ≠ I), and
# (b) the origin shift needed to put atoms at the canonical origin
# of the colored space group.
try:
ds = spglib.get_symmetry_dataset((lat, frac, types), symprec=symprec)
except Exception as exc: # pragma: no cover - spglib internal failure
log.warning("get_symmetry_dataset failed (%s); keeping working cell",
exc)
return crystal
if ds is None:
return crystal
from mag4._spglib_compat import ds_get
T = np.asarray(ds_get(ds, "transformation_matrix"), dtype=float)
p = np.asarray(ds_get(ds, "origin_shift"), dtype=float)
T_is_identity = np.allclose(T, np.eye(3), atol=1e-3)
p_is_zero = np.allclose(p, 0.0, atol=1e-6)
# No reduction (T ≈ I) → apply origin shift manually, preserving
# input atom order so the spin-vector indexing in
# :func:`_spin_per_atom` keeps matching the labels. We pick this
# path even when ``p == 0`` (a no-op) so the function has a single
# behaviour for primitive groups.
if T_is_identity:
if p_is_zero:
return crystal
new_frac = np.mod(frac + p[np.newaxis, :], 1.0)
# Build a copy of the input crystal with shifted positions but
# original labels / element ordering preserved.
elements_order = list(crystal.elements)
new_all_species: List[list] = [[] for _ in elements_order]
global_idx = 1
flat_in = [a for grp in crystal.all_species for a in grp]
flat_groups = [crystal.elements.index(_elements_per_atom(crystal)[i])
for i in range(len(flat_in))]
for in_i, atom in enumerate(flat_in):
x, y, z = (float(new_frac[in_i, 0]),
float(new_frac[in_i, 1]),
float(new_frac[in_i, 2]))
# Inherit the input atom's label so _spin_per_atom keeps
# finding it in target_species at the same index.
new_all_species[flat_groups[in_i]].append(
[atom[0], x, y, z, 0, global_idx]
)
global_idx += 1
new_cell = list(crystal.cell)
new_lat_mat = np.asarray(crystal.lattice_matrix, dtype=float)
new_lat_vecs = np.asarray(crystal.lattice_vectors, dtype=float)
target_group_idx = crystal.all_species.index(crystal.target_species)
return CrystalData(
cell=new_cell,
all_species=new_all_species,
target_species=new_all_species[target_group_idx],
nb_species=[len(g) for g in new_all_species],
total_atoms=global_idx,
elements=elements_order,
lattice_matrix=new_lat_mat,
lattice_vectors=new_lat_vecs,
)
# T ≠ I → genuine cell reduction (e.g. 1×1×2 supercell → F-cubic
# conv cell). Let spglib standardize, then re-label the smaller
# set of atoms. spglib *may* reorder atoms here; for label/spin
# consistency we inherit each output atom's label from the **first**
# input atom that lands on it after the (T, p) transformation.
try:
result = spglib.standardize_cell(
(lat, frac, types),
no_idealize=False,
to_primitive=False,
symprec=symprec,
)
except Exception as exc: # pragma: no cover - spglib internal failure
log.warning("standardize_cell failed (%s); keeping working cell", exc)
return crystal
if result is None:
return crystal
new_lat, new_frac, new_types = result
new_lat = np.asarray(new_lat, dtype=float)
new_frac = np.asarray(new_frac, dtype=float)
new_types = list(int(t) for t in new_types)
# Recover lattice parameters (a, b, c, α, β, γ).
from pymatgen.core import Lattice
pmg_lat = Lattice(new_lat)
new_cell = [pmg_lat.a, pmg_lat.b, pmg_lat.c,
pmg_lat.alpha, pmg_lat.beta, pmg_lat.gamma]
# Map each output atom to the first matching input atom.
# Empirically (cf. NiO config1 with T = diag(1, 1, 2)), spglib's
# standardized fractional position of an input atom is
# ``x_std = T @ x_input + origin_shift`` where T is the
# ``transformation_matrix``. Atoms in the same centring coset of
# the input cell collapse onto the same output position.
input_elements = _elements_per_atom(crystal)
input_labels = [a[0] for grp in crystal.all_species for a in grp]
output_to_input: Dict[int, int] = {}
for out_i in range(new_frac.shape[0]):
for in_j in range(frac.shape[0]):
if types[in_j] != new_types[out_i]:
continue
pred = np.mod(T @ frac[in_j] + p, 1.0)
diff = np.mod(new_frac[out_i] - pred + 0.5, 1.0) - 0.5
if np.max(np.abs(diff)) < 1e-3:
output_to_input[out_i] = in_j
break
# Reverse-map: each input atom owns at most one output atom (the
# first output that picked it). We keep input atoms in input order
# so target_species stays aligned with the original spin_vector.
input_to_output: Dict[int, int] = {}
for out_i, in_j in output_to_input.items():
input_to_output.setdefault(in_j, out_i)
elements_order = list(crystal.elements)
new_all_species: List[list] = [[] for _ in elements_order]
global_idx = 1
for in_j in sorted(input_to_output.keys()):
out_i = input_to_output[in_j]
el = input_elements[in_j]
x, y, z = (float(new_frac[out_i, 0]) % 1.0,
float(new_frac[out_i, 1]) % 1.0,
float(new_frac[out_i, 2]) % 1.0)
# Inherit input label so spin_vector indexing keeps working.
label = input_labels[in_j]
new_all_species[elements_order.index(el)].append(
[label, x, y, z, 0, global_idx]
)
global_idx += 1
target_group_idx = crystal.all_species.index(crystal.target_species)
new_target_species = new_all_species[target_group_idx]
from mag4.lattice import lattice_matrix
new_lat_mat = lattice_matrix(*new_cell)
new_lat_vecs = np.array([new_lat_mat[:, 0],
new_lat_mat[:, 1],
new_lat_mat[:, 2]])
return CrystalData(
cell=new_cell,
all_species=new_all_species,
target_species=new_target_species,
nb_species=[len(g) for g in new_all_species],
total_atoms=global_idx,
elements=elements_order,
lattice_matrix=new_lat_mat,
lattice_vectors=new_lat_vecs,
)
#: Backwards-compatible alias. The function used to return the
#: rhombohedral primitive for F-centered groups (which WIEN2k's
#: ``init_lapw`` rejected), it now returns the conventional cell.
#: The old name is preserved so existing imports keep working.
reduce_to_primitive = reduce_to_wien2k_cell
[docs]
def crystal_to_pymatgen(crystal: CrystalData):
"""Build a ``pymatgen.core.Structure`` from a :class:`CrystalData`.
Magnetic atoms (``crystal.target_species``) are placed **first** in the
structure so that their mag4 supercell indices coincide with their
pymatgen indices. This makes it trivial for :func:`find_dimers` to look
up bridging-angle fingerprints by passing the same integer index it
already uses internally.
"""
from pymatgen.core import Lattice, Structure
lattice = Lattice.from_parameters(*crystal.cell)
target_group_idx = crystal.all_species.index(crystal.target_species)
species: List[str] = []
fracs: List[List[float]] = []
target_el = crystal.elements[target_group_idx]
for atom in crystal.target_species:
species.append(target_el)
fracs.append([float(atom[1]), float(atom[2]), float(atom[3])])
for grp_idx, grp in enumerate(crystal.all_species):
if grp_idx == target_group_idx:
continue
el = crystal.elements[grp_idx]
for atom in grp:
species.append(el)
fracs.append([float(atom[1]), float(atom[2]), float(atom[3])])
return Structure(lattice, species, fracs)