From 9133942a57d56c47aff7d0e608958034db19eb71 Mon Sep 17 00:00:00 2001 From: lan496 Date: Mon, 12 Oct 2026 20:37:24 +0900 Subject: [PATCH 1/3] feat: preserve mixed compositions and explicit spin frames in CIF exports Co-Authored-By: Codex --- CHANGELOG.md | 8 + docs/file-formats.md | 18 ++- src/spinforge/_cif.py | 92 +++++++++++ src/spinforge/mcif.py | 102 ++++++------ src/spinforge/scif/__init__.py | 8 + src/spinforge/scif/_reader.py | 188 ++++++++++++++++++++-- src/spinforge/scif/_writer.py | 137 +++++++++++----- tests/scif_tests/test_compositions.py | 210 ++++++++++++++++++++++++ tests/scif_tests/test_reader.py | 49 ++++++ tests/scif_tests/test_writer.py | 10 +- tests/test_mcif.py | 221 ++++++++++++++++++++++++-- 11 files changed, 916 insertions(+), 127 deletions(-) create mode 100644 tests/scif_tests/test_compositions.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 6848d64..aae60a7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,14 @@ The initial public release includes: - Public project, contribution, security, citation, and data-attribution documentation. +## Unreleased + +- Preserve mixed-site compositions, partial occupancies, oxidation states, and + selected magnetic constituents in MCIF and spinCIF exports. +- Accept explicit orthogonal spin frames in `SpinCifWriter`. +- Restore occupancies when reading spinCIF and retain constituent moments when + merging disordered sites. + ## v1.0.0 - 2026/09/06 - Initial public release diff --git a/docs/file-formats.md b/docs/file-formats.md index 5157888..171b09c 100644 --- a/docs/file-formats.md +++ b/docs/file-formats.md @@ -31,9 +31,21 @@ you need both representations. revision, [`SPINCIF_REVISION`][spinforge.scif.SPINCIF_REVISION], into each file. The reader supports files written by SpinForge and the FINDSPINGROUP reference files; it is not a general-purpose spinCIF parser. -- MCIF output supports partial occupancy of a single species at a site. It does - not support mixed solid-solution sites because those sites cannot pass - through the symmetry identification used by the writer. +- Both writers preserve mixed species, oxidation states, and partial occupancy. + Use `magnetic_species` to select constituents that share each site's + occupied-atom moment. Other constituents have zero moment. Occupancy is written + separately and is not multiplied into the moment. +- For spinCIF, pass full `site_compositions` when atomic numbers do not describe + the chemistry. `spin_rotation` accepts the explicit orthogonal matrix that + maps intrinsic spins to Cartesian spins, including an improper frame. Supply + moments already rotated by that matrix. +- `SpinCifReader.get_structure()` merges coincident constituents and retains + `species_magmoms` and `magnetic_species` site properties. Its `magmom` remains + the shared selected occupied-atom vector. It rejects unequal selected moments + at one position rather than averaging them. +- Pymatgen's structure parser does not support disordered magnetic CIFs. + Mixed MCIF output therefore requires a consumer that preserves constituent + rows and their individual occupancies and moments. To write and check either format, continue to [Export result files](export-files.md). For the scientific relation between spin and magnetic symmetry, see the diff --git a/src/spinforge/_cif.py b/src/spinforge/_cif.py index 2217454..7297378 100644 --- a/src/spinforge/_cif.py +++ b/src/spinforge/_cif.py @@ -1,5 +1,13 @@ from __future__ import annotations +from collections.abc import Sequence +from dataclasses import dataclass + +import numpy as np +from pymatgen.core import Composition, DummySpecies, Element, Species +from pymatgen.core.periodic_table import get_el_sp +from spgrep.utils import NDArrayFloat + _ATOM_SITE_LOOP_HEADER = ( "", "loop_", @@ -10,3 +18,87 @@ "_atom_site_fract_z", "_atom_site_occupancy", ) + + +@dataclass(frozen=True) +class _SiteConstituent: + species: Element | Species + occupancy: float + magnetic: bool + + +def _normalize_site_species( + site_compositions: Sequence[Composition], + magnetic_species: Sequence[Sequence[str]] | None, + magnetic_moments: NDArrayFloat, +) -> tuple[tuple[_SiteConstituent, ...], ...]: + """Validate chemistry and shared occupied-atom moments for CIF writers.""" + num_sites = len(site_compositions) + moments = np.asarray(magnetic_moments, dtype=float) + if moments.shape != (num_sites, 3) or not np.all(np.isfinite(moments)): + raise ValueError(f"magnetic_moments must be a finite ({num_sites}, 3) array") + if magnetic_species is not None and len(magnetic_species) != num_sites: + raise ValueError("magnetic_species must have one selection per site") + + normalized = [] + for index, composition in enumerate(site_compositions): + items = sorted(composition.items(), key=lambda item: str(item[0])) + if not items: + raise ValueError(f"Site {index} must contain at least one species") + for species, occupancy in items: + if not isinstance(species, (Element, Species)) or isinstance(species, DummySpecies): + raise TypeError(f"Unsupported species {species!r} at site {index}") + if isinstance(species, Species) and ( + species.spin is not None + or species.oxi_state is None + or not np.isfinite(species.oxi_state) + ): + raise ValueError(f"Unsupported species {species!r} at site {index}") + if not np.isfinite(occupancy) or occupancy <= 0: + raise ValueError(f"Site {index} occupancies must be finite and positive") + if sum(occupancy for _, occupancy in items) > 1 + 1e-8: + raise ValueError(f"Site {index} total occupancy exceeds 1") + + identities = {str(species) for species, _ in items} + if magnetic_species is None: + selected = identities + else: + selection = magnetic_species[index] + if isinstance(selection, str) or any(not isinstance(name, str) for name in selection): + raise ValueError(f"Site {index} magnetic_species must be a sequence of species") + selected = {str(get_el_sp(name)) for name in selection} + if len(selected) != len(selection): + raise ValueError(f"Site {index} magnetic_species contains duplicate species") + if not selected <= identities: + raise ValueError(f"Site {index} magnetic_species contains an absent species") + if not selected and np.any(moments[index] != 0): + raise ValueError(f"Site {index} has a nonzero moment but no selected magnetic species") + normalized.append( + tuple( + _SiteConstituent(species, float(occupancy), str(species) in selected) + for species, occupancy in items + ) + ) + return tuple(normalized) + + +def _atom_type_lines(site_species: Sequence[Sequence[_SiteConstituent]]) -> list[str]: + """Declare ionic atom types while retaining uncharged Element identities.""" + ionic_species = { + str(row.species): row.species + for site in site_species + for row in site + if isinstance(row.species, Species) + } + if not ionic_species: + return [] + return [ + "", + "loop_", + "_atom_type_symbol", + "_atom_type_oxidation_number", + *[ + f"{symbol} {species.oxi_state:.12g}" + for symbol, species in sorted(ionic_species.items()) + ], + ] diff --git a/src/spinforge/mcif.py b/src/spinforge/mcif.py index d13413b..9c5b5e6 100644 --- a/src/spinforge/mcif.py +++ b/src/spinforge/mcif.py @@ -1,12 +1,14 @@ """Write magnetic structures as mcif files (BNS/MAGNDATA convention). -Reading needs no local code: ``pymatgen.core.Structure.from_file`` handles -mcif, and its symmetry expansion of ``_atom_site_moment.crystalaxis_*`` +For single-species sites, ``pymatgen.core.Structure.from_file`` handles mcif, +and its symmetry expansion of ``_atom_site_moment.crystalaxis_*`` components is correct as of pymatgen-core v2026.7.16 (materialsproject/pymatgen-core#76): moments are converted to Cartesian before applying an operation, so |m| is preserved even when an operation mixes unequal or non-orthogonal axes (e.g. hexagonal settings; see materialsproject/pymatgen-core#76 for the old failure mode). +Mixed-species files use constituent atom rows, which pymatgen's magnetic +structure reader does not currently support. :class:`MCifWriter` has no upstream equivalent (``CifWriter(write_magmoms=True)`` dumps every atom in P1 with no magnetic operations). @@ -14,6 +16,8 @@ from __future__ import annotations +from collections.abc import Sequence + import numpy as np from moyopy import ( MagneticSpaceGroupType, @@ -24,7 +28,12 @@ from pymatgen.core.operations import MagSymmOp from pymatgen.electronic_structure.core import Magmom -from spinforge._cif import _ATOM_SITE_LOOP_HEADER +from spinforge._cif import ( + _ATOM_SITE_LOOP_HEADER, + _atom_type_lines, + _normalize_site_species, + _SiteConstituent, +) __all__ = [ "MCifWriter", @@ -43,8 +52,7 @@ class MCifWriter: magnetic operations -- this finds the magnetic space group of ``struct`` with moyopy and writes the asymmetric unit plus the ``_space_group_symop_magn_operation.xyz`` / - ``_space_group_symop_magn_centering.xyz`` loops, so the file round-trips - through ``pymatgen.core.Structure.from_file``. + ``_space_group_symop_magn_centering.xyz`` loops. The structure is symmetrized in its own input setting (the basis is not standardized). ``struct`` must carry a ``magmom`` site property (pymatgen @@ -52,12 +60,11 @@ class MCifWriter: ``_atom_site_moment.crystalaxis_*`` in the unit-axis convention read back by ``CifParser`` (``m = mx*a_hat + my*b_hat + mz*c_hat``). - Partial occupancy is preserved: each site's fractional occupancy enters the - moyopy equivalence label (so a 0.79-occupied site is never folded into a fully - occupied orbit) and is written as an ``_atom_site_occupancy`` column. Only - single-species sites are supported; a mixed solid-solution site (two species - sharing one position) raises ``NotImplementedError``, since moyopy cannot - symmetrize it. + Each constituent is written as an atom row with its occupancy and full species + identity, including oxidation state. Compositions and magnetic selections enter + the moyopy equivalence labels, so symmetry cannot exchange chemically distinct + sites. Selected constituents share the site's occupied-atom moment; unselected + constituents have zero moment. Occupancy never scales the written moment. Parameters ---------- @@ -65,6 +72,10 @@ class MCifWriter: Magnetic structure with a ``magmom`` site property. symprec, mag_symprec moyopy symmetry-search tolerances (basis units / moment units). + magnetic_species + Species identities selected at each site, such as ``[("Fe2+",), ()]``. + Defaults to all constituents. An empty selection requires a zero moment. + Selected constituents retain moment rows even when their moment is zero. """ def __init__( @@ -73,15 +84,19 @@ def __init__( symprec: float = 1e-4, *, mag_symprec: float | None = None, + magnetic_species: Sequence[Sequence[str]] | None = None, ) -> None: if "magmom" not in struct.site_properties: raise ValueError("structure must have a 'magmom' site property to write an mcif") basis = struct.lattice.matrix moments = np.array([Magmom(m).global_moment for m in struct.site_properties["magmom"]]) + site_species = _normalize_site_species( + [site.species for site in struct], magnetic_species, moments + ) cell = NonCollinearMagneticCell( basis=basis.tolist(), positions=struct.frac_coords.tolist(), - numbers=self._composition_numbers(struct), + numbers=self._composition_numbers(site_species), magnetic_moments=moments.tolist(), ) dataset = MoyoNonCollinearMagneticDataset( @@ -104,7 +119,9 @@ def __init__( inv_unit_axes = np.linalg.inv(unit_axes) # crystalaxis comps = m_cart @ inv(unit_axes) bns = MagneticSpaceGroupType(dataset.uni_number).bns_number - self._mcif_str = self._format_mcif(struct, reps, op_strings, moments, inv_unit_axes, bns) + self._mcif_str = self._format_mcif( + struct, site_species, reps, op_strings, moments, inv_unit_axes, bns + ) def __str__(self) -> str: """The symmetrized mcif as a string.""" @@ -118,13 +135,12 @@ def write_file(self, filename: str, mode: str = "wt") -> None: @staticmethod def _format_mcif( structure: Structure, + site_species: tuple[tuple[_SiteConstituent, ...], ...], reps: list[int], op_strings: list[str], moments: np.ndarray, inv_unit_axes: np.ndarray, bns: str, - *, - atol: float = 1e-6, ) -> str: """Assemble symmetrized-mcif text from the asymmetric unit and magnetic operations.""" a, b, c = structure.lattice.abc @@ -132,6 +148,7 @@ def _format_mcif( lines = [ "# Symmetrized magnetic CIF written by spinforge.mcif.MCifWriter", "data_magnetic_structure", + f"_chemical_formula_sum '{structure.composition.element_composition.formula}'", f"_space_group_magn.number_BNS {bns}", f"_cell_length_a {a:.6f}", f"_cell_length_b {b:.6f}", @@ -139,6 +156,7 @@ def _format_mcif( f"_cell_angle_alpha {alpha:.6f}", f"_cell_angle_beta {beta:.6f}", f"_cell_angle_gamma {gamma:.6f}", + *_atom_type_lines(site_species), "", "loop_", "_space_group_symop_magn_operation.id", @@ -153,22 +171,21 @@ def _format_mcif( "1 x,y,z,+1", *_ATOM_SITE_LOOP_HEADER, ] - labels: dict[int, str] = {} counts: dict[str, int] = {} - for idx in reps: - symbol, occ = MCifWriter._single_species(structure[idx]) - counts[symbol] = counts.get(symbol, 0) + 1 - label = f"{symbol}{counts[symbol]}" - labels[idx] = label - fx, fy, fz = structure[idx].frac_coords - lines.append(f"{label} {symbol} {fx:.6f} {fy:.6f} {fz:.6f} {occ:.6f}") - moment_lines = [] for idx in reps: + fx, fy, fz = structure[idx].frac_coords comps = moments[idx] @ inv_unit_axes - if np.linalg.norm(comps) <= atol: - continue - moment_lines.append(f"{labels[idx]} {comps[0]:.6f} {comps[1]:.6f} {comps[2]:.6f}") + for constituent in site_species[idx]: + symbol = constituent.species.symbol + counts[symbol] = counts.get(symbol, 0) + 1 + label = f"{symbol}{counts[symbol]}" + lines.append( + f"{label} {constituent.species} {fx:.6f} {fy:.6f} {fz:.6f} " + f"{constituent.occupancy:.12g}" + ) + if constituent.magnetic: + moment_lines.append(f"{label} {comps[0]:.6f} {comps[1]:.6f} {comps[2]:.6f}") if moment_lines: lines += [ "", @@ -183,38 +200,15 @@ def _format_mcif( return "\n".join(lines) @staticmethod - def _composition_numbers(structure: Structure) -> list[int]: - """Equivalence labels for moyopy that fold occupancy into the site identity. - - moyopy treats ``numbers`` as opaque per-site equivalence classes. Encoding the - full ``(Z, occupancy)`` signature means a magnetic operation may only map a site - onto another with the *same* species and occupancy, so partially occupied sites - are never merged into a fully occupied orbit. - """ + def _composition_numbers(site_species: tuple[tuple[_SiteConstituent, ...], ...]) -> list[int]: + """moyopy equivalence labels for complete compositions and magnetic selections.""" signatures: dict[tuple, int] = {} numbers = [] - for site in structure: - signature = tuple(sorted((sp.Z, round(occ, 6)) for sp, occ in site.species.items())) + for site in site_species: + signature = tuple((str(row.species), row.occupancy, row.magnetic) for row in site) numbers.append(signatures.setdefault(signature, len(signatures) + 1)) return numbers - @staticmethod - def _single_species(site) -> tuple[str, float]: - """Return the ``(symbol, occupancy)`` of a single-species site. - - A mixed solid-solution site (two species sharing one position) cannot be - symmetrized by moyopy nor written as one mcif ``_atom_site_*`` row, so it is - rejected explicitly rather than silently mangled. - """ - items = list(site.species.items()) - if len(items) != 1: - raise NotImplementedError( - f"Mixed-occupancy site {site.frac_coords} has species {site.species}; " - "only single-species (optionally partial) occupancy is supported" - ) - element, occ = items[0] - return element.symbol, float(occ) - @staticmethod def _snap_translation( trans: list[float], *, denom: int = 24, atol: float = 1e-3 diff --git a/src/spinforge/scif/__init__.py b/src/spinforge/scif/__init__.py index 1661230..366b740 100644 --- a/src/spinforge/scif/__init__.py +++ b/src/spinforge/scif/__init__.py @@ -22,6 +22,14 @@ it enters through ``_space_group_spin.collinear_direction_xyz`` / ``_space_group_spin.coplanar_perp_uvw``, filled from the ``spinspg.spin.SpinOnlyGroup`` the structure was enumerated with. +- ``spin_rotation`` accepts the explicit orthogonal frame used for the sampled + moments, independently of an MSG object. ``site_compositions`` preserves + partial and mixed occupancies and oxidation states; ``magnetic_species`` + assigns the shared occupied-atom moment only to selected constituents. +- The reader merges constituent rows into disordered sites. Its ``magmom`` + property remains the shared occupied-atom vector, while ``species_magmoms`` + and ``magnetic_species`` retain the constituent values and selections. Compute + a mean-site moment by weighting constituent moments by occupancy once. The dictionary is explicitly preliminary upstream ("all dictionary aspects and file formats are subject to change"): the targeted revision is pinned as diff --git a/src/spinforge/scif/_reader.py b/src/spinforge/scif/_reader.py index c27f733..0105f58 100644 --- a/src/spinforge/scif/_reader.py +++ b/src/spinforge/scif/_reader.py @@ -16,7 +16,7 @@ from os import PathLike import numpy as np -from pymatgen.core import Lattice, Structure +from pymatgen.core import Element, Lattice, Species, Structure, get_el_sp from pymatgen.core.operations import MagSymmOp from pymatgen.io.cif import CifBlock, CifFile from spgrep.utils import NDArrayFloat, is_integer_array @@ -45,6 +45,12 @@ class SpinCifReader: ``(label, symbol, frac_position, moment_uvw, multiplicity)`` per asymmetric-unit site (zero moment when the site has no ``_atom_site_spin_moment`` row). + site_occupancies + Occupancy keyed by atom-site label, defaulting to one when absent. + site_moment_subspaces + Coefficient matrices from ``symmform_uvw``, keyed by magnetic atom-site + label. Each matrix's columns span the permitted moment subspace in the + declared spin basis. Missing or unknown forms have no entry. lattice :class:`pymatgen.core.Lattice` from the cell tags. spin_basis @@ -75,15 +81,32 @@ def __init__(self, text: str, *, symprec: float = 1e-4) -> None: else: # the loop is optional in the dictionary self.spin_lattice = [(np.eye(3), np.zeros(3), 1, np.eye(3))] + moment_columns = _read_loop_columns( + block, + tuple( + f"_atom_site_spin_moment.{column}" + for column in ("label", "axis_u", "axis_v", "axis_w") + ), + optional=True, + ) + moment_labels, _, _, _ = moment_columns + if len(set(moment_labels)) != len(moment_labels): + raise ValueError("Spin-moment labels must be unique") moments = { label: np.array([float(u), float(v), float(w)]) - for label, u, v, w in zip( - *( - block.data.get(f"_atom_site_spin_moment.{column}", []) - for column in ("label", "axis_u", "axis_v", "axis_w") - ) - ) + for label, u, v, w in zip(*moment_columns, strict=True) } + self._magnetic_labels = frozenset(moments) + self.site_moment_subspaces = {} + if "_atom_site_spin_moment.symmform_uvw" in block.data: + subspace_columns = _read_loop_columns( + block, ("_atom_site_spin_moment.label", "_atom_site_spin_moment.symmform_uvw") + ) + self.site_moment_subspaces = { + label: _linear_matrix(expression, "uvw") + for label, expression in zip(*subspace_columns, strict=True) + if expression not in (".", "?") + } self.sites = [ ( label, @@ -101,6 +124,34 @@ def __init__(self, text: str, *, symprec: float = 1e-4) -> None: block.data["_atom_site_symmetry_multiplicity"], ) ] + labels = block.data["_atom_site_label"] + if len(set(labels)) != len(labels): + raise ValueError("Atom-site labels must be unique") + if not self._magnetic_labels.issubset(labels): + raise ValueError("A spin moment refers to an unknown atom-site label") + self.site_occupancies = { + label: float(occupancy) + for label, occupancy in zip( + labels, block.data.get("_atom_site_occupancy", ["1"] * len(labels)), strict=True + ) + } + if any( + not np.isfinite(occupancy) or not 0 < occupancy <= 1 + for occupancy in self.site_occupancies.values() + ): + raise ValueError("Atom-site occupancies must be finite and lie in (0, 1]") + oxidation_states = _read_oxidation_states(block) + self._site_species = {} + for label, symbol, _, moment, _ in self.sites: + species = get_el_sp(symbol) + oxidation = oxidation_states.get(symbol) + if oxidation is not None: + species = Species(species.symbol, oxidation) + if not isinstance(species, Element | Species): + raise ValueError(f"Unsupported atom-site species: {symbol}") + if not np.all(np.isfinite(moment)): + raise ValueError(f"Spin moment for {label} must be finite") + self._site_species[label] = species self.lattice = Lattice.from_parameters( *(float(block.data[f"_cell_length_{name}"]) for name in "abc"), @@ -126,14 +177,29 @@ def expanded_operations(self) -> list[SymmetryOperation]: def get_structure(self) -> Structure: """Expand the asymmetric unit into the full magnetic structure. - Moments are returned as Cartesian ``magmom`` site properties in the - lattice's own Cartesian frame. Raises if an orbit disagrees with its - ``_atom_site_symmetry_multiplicity`` or transforms inconsistently. + Coincident constituent rows are merged into one disordered site with + their occupancies and oxidation identities. ``species_magmoms`` maps + every constituent identity to its Cartesian occupied-atom moment; + ``magnetic_species`` lists the selected constituents, identified by + the presence of a moment row (including explicit zeros). + + The existing ``magmom`` property is the shared selected occupied-atom + vector, or zero when none are selected. Unequal selected constituent + vectors are outside this reader's shared-moment model and raise an + error. No occupancy averaging is applied. All vectors use the read + lattice's Cartesian frame. + + Raises if an orbit disagrees with its declared multiplicity, transforms + inconsistently, or repeats the same constituent at a position. """ - species, coords, magmoms = [], [], [] + compositions: list[dict[Element | Species, float]] = [] + coords: list[np.ndarray] = [] + species_magmoms: list[dict[str, np.ndarray]] = [] + selections: list[set[str]] = [] + operations = self.expanded_operations() for label, symbol, frac, moment_uvw, multiplicity in self.sites: orbit: list[tuple[np.ndarray, np.ndarray]] = [] - for w, t, _tr, u in self.expanded_operations(): + for w, t, _tr, u in operations: position = (w @ frac + t) % 1.0 moment = u @ moment_uvw for known_position, known_moment in orbit: @@ -151,10 +217,100 @@ def get_structure(self) -> Structure: f"Orbit of {label} has {len(orbit)} sites, multiplicity says {multiplicity}" ) for position, moment in orbit: - species.append(symbol) - coords.append(position) - magmoms.append(moment @ self.spin_basis) - return Structure(self.lattice, species, coords, site_properties={"magmom": magmoms}) + index = next( + ( + index + for index, known_position in enumerate(coords) + if is_integer_array(position - known_position, atol=self._symprec) + ), + len(coords), + ) + if index == len(coords): + coords.append(position) + compositions.append({}) + species_magmoms.append({}) + selections.append(set()) + species = self._site_species[label] + identity = str(species) + if species in compositions[index]: + raise ValueError(f"Repeated constituent {identity} at {position}") + compositions[index][species] = self.site_occupancies[label] + species_magmoms[index][identity] = moment @ self.spin_basis + if label in self._magnetic_labels: + selections[index].add(identity) + + magmoms = [] + for composition, moments, selected in zip( + compositions, species_magmoms, selections, strict=True + ): + if sum(composition.values()) > 1 + 1e-6: + raise ValueError("Constituent occupancies at one site sum to more than one") + shared = moments[min(selected)] if selected else np.zeros(3) + if any( + not np.allclose(moments[identity], shared, atol=self._symprec, rtol=0) + for identity in selected + ): + raise ValueError( + "Unequal selected constituent moments cannot be represented by the shared " + "occupied-atom magmom property" + ) + magmoms.append(shared) + return Structure( + self.lattice, + compositions, + coords, + site_properties={ + "magmom": magmoms, + "species_magmoms": species_magmoms, + "magnetic_species": [tuple(sorted(selected)) for selected in selections], + }, + ) + + +def _read_loop_columns( + block: CifBlock, tags: tuple[str, ...], *, optional: bool = False +) -> list[list[str]]: + """Read named loop columns without silently truncating malformed categories.""" + if optional and not any(tag in block.data for tag in tags): + return [[] for _ in tags] + columns = [] + for tag in tags: + column = block.data.get(tag) + if not isinstance(column, list): + raise ValueError(f"Expected equal-length loop columns for {', '.join(tags)}") + columns.append(column) + if len({len(column) for column in columns}) != 1: + raise ValueError(f"Expected equal-length loop columns for {', '.join(tags)}") + return columns + + +def _read_oxidation_states(block: CifBlock) -> dict[str, float]: + """Read unambiguous atom types, checking ionic tokens against their charges.""" + tags = ("_atom_type_symbol",) + has_oxidation = "_atom_type_oxidation_number" in block.data + if has_oxidation: + tags += ("_atom_type_oxidation_number",) + columns = _read_loop_columns(block, tags, optional=True) + symbols = columns[0] + if len(set(symbols)) != len(symbols): + raise ValueError("Atom-type symbols must be unique") + if not has_oxidation: + return {} + _, values = columns + oxidation_states = {} + for symbol, value in zip(symbols, values, strict=True): + if value in (".", "?"): + continue + oxidation = float(value) + if not np.isfinite(oxidation): + raise ValueError("Oxidation states must be finite") + species = get_el_sp(symbol) + if isinstance(species, Species) and species.oxi_state != oxidation: + raise ValueError( + f"Oxidation state for {symbol} contradicts its explicit species identity" + ) + oxidation_states[symbol] = oxidation + return oxidation_states def _read_symop_loop(block: CifBlock, category: str) -> list[SymmetryOperation]: diff --git a/src/spinforge/scif/_writer.py b/src/spinforge/scif/_writer.py index 480749c..4c0cfd0 100644 --- a/src/spinforge/scif/_writer.py +++ b/src/spinforge/scif/_writer.py @@ -3,15 +3,21 @@ from __future__ import annotations -from collections.abc import Sequence +from collections.abc import Mapping, Sequence import numpy as np from moyopy import SpaceGroupType -from pymatgen.core import Element, Lattice, Structure +from pymatgen.core import Composition, Element, Lattice, Structure +from pymatgen.electronic_structure.core import Magmom from spgrep.utils import NDArrayFloat from spinspg.spin import SpinOnlyGroup, SpinOnlyGroupType -from spinforge._cif import _ATOM_SITE_LOOP_HEADER +from spinforge._cif import ( + _ATOM_SITE_LOOP_HEADER, + _atom_type_lines, + _normalize_site_species, + _SiteConstituent, +) from spinforge.configuration import SpinSymmetryAdaptedStructure from spinforge.msg import MagneticSpaceSubgroup from spinforge.ssg import NontrivialSpinSpaceGroup @@ -60,11 +66,26 @@ class SpinCifWriter: type and axis fill ``_space_group_spin.collinear_direction_xyz`` / ``coplanar_perp_uvw``. magnetic_moments - Cartesian magnetic moments per supercell site, shape + Cartesian moments per occupied selected magnetic atom, shape ``(num_supercell_sites, 3)``, exactly as passed to :meth:`SpinSymmetryAdaptedStructure.generate_with_magnetic_moments`. For structures from ``generate_oriented`` use :meth:`from_oriented`, which also conjugates the spin operations by the MSG rotation. + spin_rotation + Orthogonal 3 by 3 matrix taking intrinsic spins to spatial Cartesian + spins. Moments must already include this rotation, as ``intrinsic @ + spin_rotation.T``. Both proper and improper frames are supported; + ``None`` uses the identity. The operations and SSA restrictions are + rotated with the same matrix. + site_compositions + Full composition of each supercell site, including fractional + occupancies and oxidation states. Defaults to the supercell's atomic + elements. Cannot be combined with ``site_occupancies``. + magnetic_species + Species identities selected at each supercell site, as strings such as + ``"Fe2+"``. Selected constituents share the site's occupied-atom moment; + unselected constituents carry zero. ``None`` selects every constituent. + Compositions and selections must be invariant under the SSG operations. site_occupancies Fractional occupancy per supercell site (``None`` = fully occupied). Written as ``_atom_site_occupancy``; must be invariant under the SSG @@ -87,6 +108,9 @@ def __init__( magnetic_moments: NDArrayFloat, *, spin_only_group: SpinOnlyGroup, + spin_rotation: NDArrayFloat | None = None, + site_compositions: Sequence[Composition | Mapping[str, float]] | None = None, + magnetic_species: Sequence[Sequence[str]] | None = None, data_name: str = "spinforge", parent_space_group_type: SpaceGroupType | None = None, site_occupancies: Sequence[float] | NDArrayFloat | None = None, @@ -97,10 +121,12 @@ def __init__( nssg, np.asarray(magnetic_moments, dtype=float), spin_only_group=spin_only_group, - q_mat=np.eye(3), + q_mat=np.eye(3) if spin_rotation is None else np.asarray(spin_rotation, dtype=float), data_name=data_name, parent_space_group_type=parent_space_group_type, site_occupancies=site_occupancies, + site_compositions=site_compositions, + magnetic_species=magnetic_species, atol=atol, ) @@ -113,6 +139,8 @@ def from_oriented( msg: MagneticSpaceSubgroup, *, spin_only_group: SpinOnlyGroup, + site_compositions: Sequence[Composition | Mapping[str, float]] | None = None, + magnetic_species: Sequence[Sequence[str]] | None = None, data_name: str = "spinforge", parent_space_group_type: SpaceGroupType | None = None, site_occupancies: Sequence[float] | NDArrayFloat | None = None, @@ -124,22 +152,22 @@ def from_oriented( conjugates the spin parts of the SSG operations by the MSG rotation ``msg.Q`` that ``generate_oriented`` applied to them. """ - writer = cls.__new__(cls) moments = np.array( - [np.asarray(m) for m in structure.site_properties["magmom"]], dtype=float + [Magmom(m).global_moment for m in structure.site_properties["magmom"]], dtype=float ) - writer._build( + return cls( sas, nssg, moments, spin_only_group=spin_only_group, - q_mat=np.asarray(msg.Q, dtype=float), + spin_rotation=np.asarray(msg.Q, dtype=float), data_name=data_name, parent_space_group_type=parent_space_group_type, site_occupancies=site_occupancies, + site_compositions=site_compositions, + magnetic_species=magnetic_species, atol=atol, ) - return writer def _build( self, @@ -152,6 +180,8 @@ def _build( data_name: str, parent_space_group_type: SpaceGroupType | None, site_occupancies: Sequence[float] | NDArrayFloat | None = None, + site_compositions: Sequence[Composition | Mapping[str, float]] | None = None, + magnetic_species: Sequence[Sequence[str]] | None = None, atol: float, ) -> None: positions = np.asarray(sas.supercell.positions) @@ -160,8 +190,16 @@ def _build( f"magnetic_moments shape {moments.shape} does not match " f"({len(positions)}, 3) supercell sites" ) - if not np.allclose(q_mat @ q_mat.T, np.eye(3), rtol=0, atol=atol): - raise ValueError(f"MSG rotation Q must be orthogonal, got {q_mat}") + if not np.all(np.isfinite(moments)): + raise ValueError("magnetic_moments must be finite") + if ( + q_mat.shape != (3, 3) + or not np.all(np.isfinite(q_mat)) + or not np.allclose(q_mat @ q_mat.T, np.eye(3), rtol=0, atol=atol) + ): + raise ValueError("spin_rotation must be a finite orthogonal 3 by 3 matrix") + if site_compositions is not None and site_occupancies is not None: + raise ValueError("Use either site_compositions or site_occupancies, not both") lattice = np.asarray(sas.supercell.basis) frame = _crystal_cartesian_frame(lattice) # columns = spin basis vectors (Cartesian) @@ -191,6 +229,26 @@ def _build( occupancies = _validated_occupancies( site_occupancies, len(positions), permutations, atol=atol ) + compositions = ( + [Composition(composition) for composition in site_compositions] + if site_compositions is not None + else [ + Composition({Element.from_Z(number): occupancy}) + for number, occupancy in zip(sas.supercell.numbers, occupancies, strict=True) + ] + ) + constituents = _normalize_site_species(compositions, magnetic_species, moments) + signatures = [ + tuple((str(c.species), c.occupancy, c.magnetic) for c in site) for site in constituents + ] + for permutation in permutations: + if any( + signatures[index] != signatures[int(image)] + for index, image in enumerate(permutation) + ): + raise ValueError( + "site_compositions and magnetic_species must be invariant under the SSG operations" + ) collinear_xyz = _collinear_direction_xyz(spin_only_group, q_mat, lattice) coplanar_uvw = _coplanar_perp_uvw(spin_only_group, q_mat, frame) @@ -200,10 +258,9 @@ def _build( coplanar_uvw=coplanar_uvw, lattice=lattice, frame=frame, - numbers=sas.supercell.numbers, + constituents=constituents, positions=positions, moments=moments, - occupancies=occupancies, orbits=orbits, operations=operations, spin_lattice=spin_lattice, @@ -293,7 +350,11 @@ def _validated_occupancies( f"site_occupancies shape {occupancies.shape} does not match " f"({num_sites},) supercell sites" ) - if np.any(occupancies <= 0.0) or np.any(occupancies > 1.0 + atol): + if ( + not np.all(np.isfinite(occupancies)) + or np.any(occupancies <= 0.0) + or np.any(occupancies > 1.0 + atol) + ): raise ValueError(f"site_occupancies must lie in (0, 1], got {occupancies.tolist()}") for permutation in permutations: permuted = occupancies[np.asarray(permutation, dtype=int)] @@ -312,10 +373,9 @@ def _format_scif( coplanar_uvw: str, lattice: NDArrayFloat, frame: NDArrayFloat, - numbers: list[int], + constituents: tuple[tuple[_SiteConstituent, ...], ...], positions: NDArrayFloat, moments: NDArrayFloat, - occupancies: NDArrayFloat, orbits: list[tuple[int, int]], operations: list[tuple[np.ndarray, np.ndarray, np.ndarray]], spin_lattice: list[tuple[np.ndarray, np.ndarray]], @@ -370,38 +430,37 @@ def _format_scif( for i, (t, u) in enumerate(spin_lattice): lines.append(f"{i + 1} {_xyzt(np.eye(3, dtype=int), t, u)} {_uvw_expr(u, atol=atol)}") - symbols = [Element.from_Z(z).symbol for z in numbers] - lines += ["", "loop_", "_atom_type_symbol"] - lines += list(dict.fromkeys(symbols)) + lines += _atom_type_lines(constituents) lines += (*_ATOM_SITE_LOOP_HEADER, "_atom_site_symmetry_multiplicity") - labels = [] + magnetic_labels: list[list[str]] = [] counts: dict[str, int] = {} for rep, multiplicity in orbits: - symbol = symbols[rep] - counts[symbol] = counts.get(symbol, 0) + 1 - label = f"{symbol}{counts[symbol]}" - labels.append(label) + labels = [] fx, fy, fz = positions[rep] - lines.append( - f"{label} {symbol} {_fmt(fx)} {_fmt(fy)} {_fmt(fz)} " - f"{_fmt(float(occupancies[rep]))} {multiplicity}" - ) + for constituent in constituents[rep]: + symbol = constituent.species.symbol + counts[symbol] = counts.get(symbol, 0) + 1 + label = f"{symbol}{counts[symbol]}" + if constituent.magnetic: + labels.append(label) + lines.append( + f"{label} {constituent.species} {_fmt(fx)} {_fmt(fy)} {_fmt(fz)} " + f"{constituent.occupancy:.12g} {multiplicity}" + ) + magnetic_labels.append(labels) moment_lines = [] - for (rep, _), label, subspace in zip(orbits, labels, site_moment_subspaces): + for (rep, _), labels, subspace in zip( + orbits, magnetic_labels, site_moment_subspaces, strict=True + ): m_uvw = frame.T @ moments[rep] components = [_fmt(x) for x in m_uvw] - # The reader treats an omitted row as zero, so omit one only when the - # written values are all zero anyway. The decision is read off the - # formatted text rather than a cutoff, so weak moments that still print - # a nonzero digit are never silently dropped. - if all(float(component) == 0.0 for component in components): - continue symmform = _symmform(subspace, atol=atol) - moment_lines.append( - f"{label} {' '.join(components)} {symmform} {_fmt(float(np.linalg.norm(m_uvw)))}" - ) + for label in labels: + moment_lines.append( + f"{label} {' '.join(components)} {symmform} {_fmt(float(np.linalg.norm(m_uvw)))}" + ) if moment_lines: lines += [ "", diff --git a/tests/scif_tests/test_compositions.py b/tests/scif_tests/test_compositions.py new file mode 100644 index 0000000..e08a5ac --- /dev/null +++ b/tests/scif_tests/test_compositions.py @@ -0,0 +1,210 @@ +"""Constituent chemistry and explicit spin frames survive spinCIF expansion.""" + +from __future__ import annotations + +import numpy as np +import pytest +from pymatgen.core import Composition +from pymatgen.io.cif import CifFile +from spinspg.spin import SpinOnlyGroupType + +from spinforge.scif import SpinCifReader, SpinCifWriter + + +@pytest.fixture +def antitranslation_model(enumerated): + return next( + (sog, group, ssa) + for sog, group, ssa in enumerated("MnTe", SpinOnlyGroupType.COLLINEAR, 2) + if group.spin_translation_coset.size > 1 and ssa.dim == 1 + ) + + +def _chemistry(ssa, magnetic_composition, selected): + compositions = [ + Composition(magnetic_composition if number == 25 else {"Te2-": 0.9}) + for number in ssa.supercell.numbers + ] + selections = [selected if number == 25 else () for number in ssa.supercell.numbers] + return compositions, selections + + +@pytest.mark.parametrize( + ("magnetic_composition", "selected"), + [ + ({"Mn2+": 0.6, "Zn2+": 0.2}, ("Mn2+",)), + ({"Mn2+": 0.6, "Fe2+": 0.2}, ("Mn2+", "Fe2+")), + ({"Mn2+": 0.6, "Mn3+": 0.2}, ("Mn2+",)), + ({"Mn": 0.6, "Mn0+": 0.2}, ("Mn0+",)), + ({"Mn2+": 0.6, "Zn2+": 1e-8}, ("Mn2+",)), + ], +) +@pytest.mark.parametrize("improper", [False, True]) +def test_mixed_antitranslation_with_explicit_frame( + antitranslation_model, magnetic_composition, selected, improper +): + sog, group, ssa = antitranslation_model + rotation = np.array([[0.0, 0.0, 1.0], [0.0, 1.0, 0.0], [-1.0, 0.0, 0.0]]) + if improper: + rotation = rotation @ np.diag([1.0, 1.0, -1.0]) + (intrinsic,) = ssa.magnetic_moments_basis + moments = intrinsic @ rotation.T + compositions, selections = _chemistry(ssa, magnetic_composition, selected) + text = str( + SpinCifWriter( + ssa, + group, + moments, + spin_only_group=sog, + spin_rotation=rotation, + site_compositions=compositions, + magnetic_species=selections, + ) + ) + reader = SpinCifReader(text) + structure = reader.get_structure() + assert len(structure) == ssa.supercell.num_supercell_sites + assert any(reversal == -1 for _, _, reversal, _ in reader.spin_lattice) + embedding = np.linalg.inv(structure.lattice.matrix) @ ssa.supercell.basis + for site in structure: + indices = [ + index + for index, position in enumerate(ssa.supercell.positions) + if np.allclose( + site.frac_coords - position - np.rint(site.frac_coords - position), + 0, + atol=1e-5, + ) + ] + (index,) = indices + assert site.species.as_dict() == compositions[index].as_dict() + assert set(site.properties["magnetic_species"]) == set(selections[index]) + np.testing.assert_allclose( + site.properties["magmom"] @ embedding, moments[index], atol=1e-5 + ) + assert set(site.properties["species_magmoms"]) == { + str(species) for species in compositions[index] + } + mean = np.zeros(3) + for species, occupancy in compositions[index].items(): + actual = site.properties["species_magmoms"][str(species)] @ embedding + expected = moments[index] if str(species) in selections[index] else np.zeros(3) + np.testing.assert_allclose(actual, expected, atol=1e-5) + mean += occupancy * actual + expected_mean = moments[index] * sum( + occupancy + for species, occupancy in compositions[index].items() + if str(species) in selections[index] + ) + np.testing.assert_allclose(mean, expected_mean, atol=1e-5) + + direction = np.array([float(value) for value in reader.collinear_direction_xyz.split(",")]) + actual_axis = direction @ ssa.supercell.basis + expected_axis = rotation @ sog.axis + np.testing.assert_allclose(np.cross(actual_axis, expected_axis), 0, atol=1e-5) + + original_spin_basis = reader.spin_basis @ embedding + operations, _ = group.get_full_operations_and_table() + transformation = np.asarray(ssa.supercell.sublattice.transformation) + inverse = np.linalg.inv(transformation) + expanded = reader.expanded_operations() + assert len(expanded) == len(operations.rotations) + for w, t, u in zip( + operations.rotations, operations.translations, operations.spin_rotations, strict=True + ): + expected_w = inverse @ w @ transformation + expected_t = inverse @ t + expected_u = original_spin_basis @ rotation @ u @ rotation.T @ original_spin_basis.T + assert any( + np.allclose(read_w, expected_w, atol=1e-6) + and np.allclose(read_t - expected_t - np.rint(read_t - expected_t), 0, atol=1e-6) + and np.allclose(read_u, expected_u, atol=1e-5) + for read_w, read_t, _, read_u in expanded + ) + + +def test_selected_zero_moments_remain_selected(antitranslation_model): + sog, group, ssa = antitranslation_model + compositions, selections = _chemistry(ssa, {"Mn2+": 0.6, "Zn2+": 0.2}, ("Mn2+",)) + text = str( + SpinCifWriter( + ssa, + group, + np.zeros((ssa.supercell.num_supercell_sites, 3)), + spin_only_group=sog, + site_compositions=compositions, + magnetic_species=selections, + ) + ) + (block,) = CifFile.from_str(text).data.values() + moment_labels = set(block.data["_atom_site_spin_moment.label"]) + assert moment_labels == { + label + for label, symbol in zip( + block.data["_atom_site_label"], block.data["_atom_site_type_symbol"], strict=True + ) + if symbol == "Mn2+" + } + structure = SpinCifReader(text).get_structure() + assert all( + site.properties["magnetic_species"] == (("Mn2+",) if "Mn2+" in site.species else ()) + for site in structure + ) + np.testing.assert_array_equal(structure.site_properties["magmom"], 0) + + +@pytest.mark.parametrize("change_selection", [False, True]) +def test_orbit_inconsistent_constituents_are_rejected(antitranslation_model, change_selection): + sog, group, ssa = antitranslation_model + (moments,) = ssa.magnetic_moments_basis + compositions, selections = _chemistry(ssa, {"Mn2+": 0.6, "Fe2+": 0.2}, ("Mn2+",)) + index = next(index for index, number in enumerate(ssa.supercell.numbers) if number == 25) + if change_selection: + selections[index] = ("Fe2+",) + else: + compositions[index] = Composition({"Mn2+": 0.5, "Fe2+": 0.2}) + with pytest.raises(ValueError, match="invariant under the SSG operations"): + SpinCifWriter( + ssa, + group, + moments, + spin_only_group=sog, + site_compositions=compositions, + magnetic_species=selections, + ) + + +@pytest.mark.parametrize("rotation", [np.eye(2), np.ones((3, 3)), np.eye(3) * np.nan]) +def test_invalid_explicit_frames_are_rejected(antitranslation_model, rotation): + sog, group, ssa = antitranslation_model + (moments,) = ssa.magnetic_moments_basis + with pytest.raises(ValueError, match="finite orthogonal 3 by 3"): + SpinCifWriter(ssa, group, moments, spin_only_group=sog, spin_rotation=rotation) + + +def test_compositions_and_legacy_occupancies_are_exclusive(antitranslation_model): + sog, group, ssa = antitranslation_model + (moments,) = ssa.magnetic_moments_basis + compositions, _ = _chemistry(ssa, {"Mn": 0.5}, ("Mn",)) + with pytest.raises(ValueError, match="either site_compositions or site_occupancies"): + SpinCifWriter( + ssa, + group, + moments, + spin_only_group=sog, + site_compositions=compositions, + site_occupancies=[0.5] * len(compositions), + ) + + +def test_reader_rejects_different_selected_constituent_moments(external_scif): + text = external_scif.replace( + "_atom_site_symmetry_multiplicity\nMn1 Mn 0.25 0.0 0.0 2", + "_atom_site_symmetry_multiplicity\n_atom_site_occupancy\n" + "Mn1 Mn 0.25 0.0 0.0 2 0.5\nFe1 Fe 0.25 0.0 0.0 2 0.5", + ).replace( + "Mn1 0.0 0.0 3.0 0,0,w 0,0,w 3.0", + "Mn1 0.0 0.0 3.0 0,0,w 0,0,w 3.0\nFe1 0.0 0.0 2.0 0,0,w 0,0,w 2.0", + ) + with pytest.raises(ValueError, match="Unequal selected constituent moments"): + SpinCifReader(text).get_structure() diff --git a/tests/scif_tests/test_reader.py b/tests/scif_tests/test_reader.py index 656dfd1..7f00e74 100644 --- a/tests/scif_tests/test_reader.py +++ b/tests/scif_tests/test_reader.py @@ -21,6 +21,55 @@ def test_reader_external_style_columns(external_scif): assert np.allclose(structure.site_properties["magmom"], [[0, 0, 3.0], [0, 0, 3.0]]) +def test_reader_rejects_duplicate_moment_labels(external_scif): + text = external_scif + "Mn1 0.0 0.0 3.0 0,0,w 0,0,w 3.0\n" + with pytest.raises(ValueError, match="Spin-moment labels must be unique"): + SpinCifReader(text) + + +@pytest.mark.parametrize( + "columns", + [ + "loop_\n_atom_site_spin_moment.label\nMn1 Mn2\n" + "loop_\n_atom_site_spin_moment.axis_u\n_atom_site_spin_moment.axis_v\n" + "_atom_site_spin_moment.axis_w\n0 0 3\n", + "loop_\n_atom_site_spin_moment.label\n_atom_site_spin_moment.axis_u\n" + "_atom_site_spin_moment.axis_v\n_atom_site_spin_moment.axis_w\nMn1 0 0 3\n" + "loop_\n_atom_site_spin_moment.symmform_uvw\n0,0,w 0,0,w\n", + "_atom_site_spin_moment.label Mn1\n_atom_site_spin_moment.axis_u 0\n" + "_atom_site_spin_moment.axis_v 0\n_atom_site_spin_moment.axis_w 3\n", + ], +) +def test_reader_rejects_malformed_moment_columns(external_scif, columns): + text = external_scif.split("loop_\n_atom_site_spin_moment.label")[0] + columns + with pytest.raises(ValueError, match="equal-length loop columns"): + SpinCifReader(text) + + +@pytest.mark.parametrize( + ("type_rows", "message"), + [ + ("Mn2+ 3", "contradicts its explicit species identity"), + ("Mn2+ 2\nMn2+ 2", "Atom-type symbols must be unique"), + ("Mn2+ nan", "Oxidation states must be finite"), + ("Mn2+ inf", "Oxidation states must be finite"), + ], +) +def test_reader_rejects_invalid_oxidation_metadata(external_scif, type_rows, message): + text = external_scif.replace("Mn1 Mn ", "Mn1 Mn2+ ") + text += "loop_\n_atom_type_symbol\n_atom_type_oxidation_number\n" + type_rows + "\n" + with pytest.raises(ValueError, match=message): + SpinCifReader(text) + + +def test_reader_rejects_mismatched_atom_type_columns(external_scif): + text = external_scif + ( + "loop_\n_atom_type_symbol\nMn2+ Fe2+\nloop_\n_atom_type_oxidation_number\n2\n" + ) + with pytest.raises(ValueError, match="equal-length loop columns"): + SpinCifReader(text) + + def test_reader_rejects_non_grammar_expressions(external_scif): # uvw values from the file must never reach code execution; anything # outside the linear-expression grammar is rejected while parsing. diff --git a/tests/scif_tests/test_writer.py b/tests/scif_tests/test_writer.py index c7a676e..1f033ef 100644 --- a/tests/scif_tests/test_writer.py +++ b/tests/scif_tests/test_writer.py @@ -52,10 +52,12 @@ def test_site_occupancies_are_written_and_checked(enumerated): text = str( SpinCifWriter(sas, nssg, moments, spin_only_group=sog, site_occupancies=occupancies) ) - mn_rows = [line for line in text.splitlines() if line.startswith("Mn1 Mn ")] - assert len(mn_rows) == 1 and "0.500000" in mn_rows[0] - te_rows = [line for line in text.splitlines() if line.startswith("Te1 Te ")] - assert len(te_rows) == 1 and "1.000000" in te_rows[0] + structure = SpinCifReader(text).get_structure() + assert len(structure) == len(sas.supercell.positions) + assert all( + site.species.get_el_amt_dict() == ({"Mn": 0.5} if "Mn" in site.species else {"Te": 1.0}) + for site in structure + ) with pytest.raises(ValueError, match="does not match"): SpinCifWriter(sas, nssg, moments, spin_only_group=sog, site_occupancies=[0.5]) diff --git a/tests/test_mcif.py b/tests/test_mcif.py index 7532fcb..bb9c99f 100644 --- a/tests/test_mcif.py +++ b/tests/test_mcif.py @@ -10,7 +10,9 @@ import numpy as np import pytest -from pymatgen.core import Structure +from pymatgen.core import Composition, DummySpecies, Lattice, Species, Structure +from pymatgen.electronic_structure.core import Magmom +from pymatgen.io.cif import CifFile, CifParser from spinforge.mcif import MCifWriter from spinforge.testing import get_magnetic_space_group_type @@ -147,7 +149,8 @@ def test_partial_occupancy_roundtrips(tmp_path) -> None: text = out.read_text() assert "_atom_site_occupancy" in text - assert "0.800000" in text # the partial Mn occupancy survived into the asym unit + data = next(iter(CifFile.from_str(text).data.values())).data + assert 0.8 in map(float, data["_atom_site_occupancy"]) reloaded = Structure.from_file(str(out)) mn_occupancies = sorted( @@ -160,15 +163,211 @@ def test_partial_occupancy_roundtrips(tmp_path) -> None: assert sum(abs(occ - 0.8) < 1e-6 for occ in mn_occupancies) == 1 -def test_save_rejects_mixed_solid_solution() -> None: - from pymatgen.core import Lattice, Structure - from pymatgen.electronic_structure.core import Magmom - +def test_mixed_site_retains_minority_species_and_occupied_atom_moment() -> None: mixed = Structure( - Lattice.cubic(4.0), - [{"Mn": 0.5, "Ge": 0.5}], + Lattice.hexagonal(4.0, 6.0), + [{"Fe2+": 0.2, "Zn2+": 0.7}, {"Zn2+": 0.7, "Fe2+": 0.2}], + [[0, 0, 0], [0, 0, 0.5]], + site_properties={"magmom": [[1, 2, 3], [-1, -2, -3]]}, + ) + text = str(MCifWriter(mixed, magnetic_species=[("Fe2+",), ("Fe2+",)])) + images = _expand_constituent_rows(text) + assert len(images) == 4 + for symbol, occupancy, position, moment, selected in images: + is_iron = symbol == "Fe2+" + assert symbol in {"Fe2+", "Zn2+"} + assert occupancy == (0.2 if is_iron else 0.7) + assert selected == is_iron + expected = np.array([1, 2, 3]) * (1 if abs(position[2]) < 1e-5 else -1) + assert np.allclose(moment, expected if is_iron else 0, atol=2e-6) + # Reordering the composition mapping cannot alter the exported constituent model. + block = next(iter(CifFile.from_str(text).data.values())).data + assert Composition(block["_chemical_formula_sum"]).as_dict() == {"Fe": 0.4, "Zn": 1.4} + assert len(block["_atom_site_label"]) == 2 + assert len(block["_atom_site_moment.label"]) == 1 + + +@pytest.mark.parametrize("moment", [[0, 0, 2], [0, 0, 0]]) +def test_symmetry_preserves_constituent_selection_even_for_zero_moments(moment) -> None: + structure = Structure( + Lattice.cubic(4), + [{"Fe": 0.5, "Mn": 0.5}] * 2, + [[0, 0, 0], [0, 0, 0.5]], + site_properties={"magmom": [moment, moment]}, + ) + text = str(MCifWriter(structure, magnetic_species=[("Fe",), ("Mn",)])) + block = next(iter(CifFile.from_str(text).data.values())).data + assert len(block["_atom_site_label"]) == 4 + assert len(block["_atom_site_moment.label"]) == 2 + images = _expand_constituent_rows(text) + assert len(images) == 4 + for symbol, _, position, vector, selected in images: + expected_selected = symbol == ("Fe" if abs(position[2]) < 1e-5 else "Mn") + assert selected == expected_selected + assert np.allclose(vector, moment if selected else 0) + + +def test_oxidation_states_and_neutral_elements_remain_distinct() -> None: + structure = Structure( + Lattice.cubic(4), + [{"Fe": 0.25, "Fe0+": 0.25, "Fe2+": 0.5}], + [[0, 0, 0]], + site_properties={"magmom": [[0, 0, 2]]}, + ) + text = str(MCifWriter(structure, magnetic_species=[("Fe2+",)])) + data = next(iter(CifFile.from_str(text).data.values())).data + assert set(data["_atom_site_type_symbol"]) == {"Fe", "Fe0+", "Fe2+"} + assert dict(zip(data["_atom_type_symbol"], data["_atom_type_oxidation_number"])) == { + "Fe0+": "0", + "Fe2+": "2", + } + assert len(data["_atom_site_moment.label"]) == 1 + + +def test_ionic_single_species_roundtrips_with_pymatgen() -> None: + structure = Structure( + Lattice.cubic(4), + [{"Fe2+": 0.7}], [[0, 0, 0]], - site_properties={"magmom": [Magmom([0, 0, 2.0])]}, + site_properties={"magmom": [[0, 0, 3]]}, + ) + parser = CifParser.from_str(str(MCifWriter(structure))) + (data,) = parser.as_dict().values() + assert Composition(data["_chemical_formula_sum"]) == Composition({"Fe": 0.7}) + (reloaded,) = parser.parse_structures(primitive=False) + (site,) = reloaded + assert site.species == Composition({"Fe2+": 0.7}) + assert np.allclose(Magmom(site.properties["magmom"]).global_moment, [0, 0, 3]) + + +def test_oxidation_state_decorations_constrain_symmetry() -> None: + structure = Structure( + Lattice.cubic(4), + ["Fe2+", "Fe3+"], + [[0, 0, 0], [0, 0, 0.5]], + site_properties={"magmom": [[0, 0, 2], [0, 0, 2]]}, + ) + images = _expand_constituent_rows(str(MCifWriter(structure))) + assert len(images) == 2 + assert {symbol for symbol, *_ in images} == {"Fe2+", "Fe3+"} + + +def test_small_minority_occupancy_is_not_rounded_to_zero() -> None: + structure = Structure( + Lattice.cubic(4), + [{"Fe": 1e-7, "Zn": 0.9}], + [[0, 0, 0]], + site_properties={"magmom": [[0, 0, 2]]}, + ) + images = _expand_constituent_rows(str(MCifWriter(structure, magnetic_species=[("Fe",)]))) + (iron,) = (row for row in images if row[0] == "Fe") + assert iron[1] == 1e-7 + assert np.allclose(iron[3], [0, 0, 2]) + + +@pytest.mark.parametrize( + ("selection", "message"), + [ + ([], "one selection per site"), + ([("Mn",)], "absent species"), + ([("Fe", "Fe")], "duplicate species"), + ([()], "nonzero moment"), + (["Fe"], "sequence of species"), + ], +) +def test_rejects_invalid_magnetic_selection(selection, message) -> None: + structure = Structure( + Lattice.cubic(4), ["Fe"], [[0, 0, 0]], site_properties={"magmom": [[0, 0, 2]]} + ) + with pytest.raises(ValueError, match=message): + MCifWriter(structure, magnetic_species=selection) + + +@pytest.mark.parametrize( + ("composition", "exception", "message"), + [ + (Composition({"Fe": -0.2}, allow_negative=True), ValueError, "positive"), + (Composition({"Fe": np.nan}), ValueError, "at least one species"), + (Composition({DummySpecies("X"): 1}), TypeError, "Unsupported species"), + (Composition({Species("Fe", None): 1}), ValueError, "Unsupported species"), + (Composition({Species("Fe", 2, spin=1): 1}), ValueError, "Unsupported species"), + ], +) +def test_rejects_invalid_constituents(composition, exception, message) -> None: + structure = Structure( + Lattice.cubic(4), [composition], [[0, 0, 0]], site_properties={"magmom": [[0, 0, 2]]} + ) + with pytest.raises(exception, match=message): + MCifWriter(structure) + + +@pytest.mark.parametrize("value", [np.nan, np.inf]) +def test_rejects_nonfinite_moments(value) -> None: + structure = Structure( + Lattice.cubic(4), ["Fe"], [[0, 0, 0]], site_properties={"magmom": [[0, 0, value]]} + ) + with pytest.raises(ValueError, match="finite"): + MCifWriter(structure) + + +def _expand_constituent_rows(text: str) -> list[tuple[str, float, np.ndarray, np.ndarray, bool]]: + """Read CIF tokens and expand each constituent independently of the writer.""" + parser = CifParser.from_str(text) + block = next(iter(CifFile.from_str(text).data.values())) + data = block.data + lattice = Lattice.from_parameters( + *[float(data[f"_cell_length_{axis}"]) for axis in ("a", "b", "c")], + *[float(data[f"_cell_angle_{axis}"]) for axis in ("alpha", "beta", "gamma")], + ) + moments = dict( + zip( + data.get("_atom_site_moment.label", []), + zip(*[data.get(f"_atom_site_moment.crystalaxis_{axis}", []) for axis in "xyz"]), + ) ) - with pytest.raises(NotImplementedError, match="single-species"): - MCifWriter(mixed) + images = [] + for label, symbol, x, y, z, occupancy in zip( + *[ + data[name] + for name in ( + "_atom_site_label", + "_atom_site_type_symbol", + "_atom_site_fract_x", + "_atom_site_fract_y", + "_atom_site_fract_z", + "_atom_site_occupancy", + ) + ] + ): + position = np.array([x, y, z], dtype=float) + mx, my, mz = moments.get(label, (0, 0, 0)) + moment = Magmom.from_moment_relative_to_crystal_axes( + (float(mx), float(my), float(mz)), lattice + ).global_moment + for operation in parser.get_magsymops(block): + new_position = operation.operate(position) % 1 + cartesian_rotation = ( + lattice.matrix.T @ operation.rotation_matrix @ np.linalg.inv(lattice.matrix.T) + ) + new_moment = ( + np.linalg.det(cartesian_rotation) + * operation.time_reversal + * cartesian_rotation + @ moment + ) + matches = [ + row + for row in images + if row[0] == symbol + and np.allclose(row[2] - new_position, np.rint(row[2] - new_position), atol=1e-5) + ] + if matches: + (match,) = matches + assert match[1] == float(occupancy) + assert np.allclose(match[3], new_moment, atol=2e-6) + assert match[4] == (label in moments) + else: + images.append( + (symbol, float(occupancy), new_position, new_moment, label in moments) + ) + return images From 95b822f7dd305586f1f784c1c914e7c565e5f655 Mon Sep 17 00:00:00 2001 From: lan496 Date: Tue, 13 Oct 2026 07:29:12 +0900 Subject: [PATCH 2/3] test: cover caller-controlled CIF tolerances Co-Authored-By: Codex --- tests/scif_tests/test_compositions.py | 66 ++++++++++++++++++++++ tests/scif_tests/test_reader.py | 80 +++++++++++++++++++++++++++ tests/test_mcif.py | 18 ++++++ 3 files changed, 164 insertions(+) diff --git a/tests/scif_tests/test_compositions.py b/tests/scif_tests/test_compositions.py index e08a5ac..c133e7c 100644 --- a/tests/scif_tests/test_compositions.py +++ b/tests/scif_tests/test_compositions.py @@ -197,6 +197,72 @@ def test_compositions_and_legacy_occupancies_are_exclusive(antitranslation_model ) +@pytest.mark.parametrize("legacy_occupancies", [False, True]) +def test_writer_occupancy_bound_uses_caller_tolerance(antitranslation_model, legacy_occupancies): + sog, group, ssa = antitranslation_model + (moments,) = ssa.magnetic_moments_basis + occupancy = 1 + 5e-9 + if legacy_occupancies: + chemistry = {"site_occupancies": [occupancy] * ssa.supercell.num_supercell_sites} + else: + compositions, selections = _chemistry(ssa, {"Mn": occupancy}, ("Mn",)) + chemistry = {"site_compositions": compositions, "magnetic_species": selections} + + with pytest.raises(ValueError, match="occupanc"): + SpinCifWriter(ssa, group, moments, spin_only_group=sog, atol=1e-12, **chemistry) + + text = str(SpinCifWriter(ssa, group, moments, spin_only_group=sog, atol=1e-6, **chemistry)) + (block,) = CifFile.from_str(text).data.values() + magnetic_occupancies = [ + float(amount) + for species, amount in zip( + block.data["_atom_site_type_symbol"], block.data["_atom_site_occupancy"], strict=True + ) + if species == "Mn" + ] + assert magnetic_occupancies + assert all(amount == occupancy for amount in magnetic_occupancies) + if legacy_occupancies: + chemistry = {"site_occupancies": [1 + 5e-7] * ssa.supercell.num_supercell_sites} + else: + compositions, selections = _chemistry(ssa, {"Mn": 1 + 5e-7}, ("Mn",)) + chemistry = {"site_compositions": compositions, "magnetic_species": selections} + with pytest.raises(ValueError, match="occupanc"): + SpinCifWriter(ssa, group, moments, spin_only_group=sog, atol=1e-6, **chemistry) + + +@pytest.mark.parametrize("legacy_occupancies", [False, True]) +def test_writer_occupancy_invariance_uses_caller_tolerance( + antitranslation_model, legacy_occupancies +): + sog, group, ssa = antitranslation_model + (moments,) = ssa.magnetic_moments_basis + index = next(index for index, number in enumerate(ssa.supercell.numbers) if number == 25) + if legacy_occupancies: + occupancies = [0.8] * ssa.supercell.num_supercell_sites + occupancies[index] += 5e-7 + chemistry = {"site_occupancies": occupancies} + else: + compositions, selections = _chemistry(ssa, {"Mn": 0.8}, ("Mn",)) + compositions[index] = Composition({"Mn": 0.8 + 5e-7}) + chemistry = {"site_compositions": compositions, "magnetic_species": selections} + + with pytest.raises(ValueError, match="invariant under the SSG operations"): + SpinCifWriter(ssa, group, moments, spin_only_group=sog, atol=1e-8, **chemistry) + + text = str(SpinCifWriter(ssa, group, moments, spin_only_group=sog, atol=1e-6, **chemistry)) + (block,) = CifFile.from_str(text).data.values() + magnetic_occupancies = [ + float(amount) + for species, amount in zip( + block.data["_atom_site_type_symbol"], block.data["_atom_site_occupancy"], strict=True + ) + if species == "Mn" + ] + assert magnetic_occupancies + np.testing.assert_allclose(magnetic_occupancies, 0.8, atol=1e-6, rtol=0) + + def test_reader_rejects_different_selected_constituent_moments(external_scif): text = external_scif.replace( "_atom_site_symmetry_multiplicity\nMn1 Mn 0.25 0.0 0.0 2", diff --git a/tests/scif_tests/test_reader.py b/tests/scif_tests/test_reader.py index 7f00e74..ae9a9e8 100644 --- a/tests/scif_tests/test_reader.py +++ b/tests/scif_tests/test_reader.py @@ -6,6 +6,7 @@ import numpy as np import pytest +from pymatgen.core import Composition from spinspg.spin import SpinOnlyGroupType from spinforge.scif import SpinCifReader, SpinCifWriter @@ -111,6 +112,85 @@ def test_reader_rejects_degenerate_spin_frame(external_scif): SpinCifReader(external_scif.replace("'a,b,c'", "'a,a,c'")) +@pytest.mark.parametrize("from_file", [False, True]) +def test_reader_spin_frame_uses_symprec(external_scif, tmp_path, from_file): + text = external_scif.replace("'a,b,c'", "'a,b+0.0005*a,c'") + source = text + read = SpinCifReader + if from_file: + path = tmp_path / "frame.scif" + path.write_text(text) + source = str(path) + read = SpinCifReader.from_file + with pytest.raises(ValueError, match="not orthogonal"): + read(source, symprec=1e-6) + structure = read(source, symprec=1e-2).get_structure() + assert len(structure) == 2 + assert np.allclose(structure.site_properties["magmom"], [[0, 0, 3], [0, 0, 3]]) + + +@pytest.mark.parametrize("scale", [1e-200, 1e-10, 1e200]) +def test_reader_normalizes_nonzero_spin_axes(external_scif, scale): + text = external_scif.replace("'a,b,c'", f"'{scale:.200f}*a,b,c'") + reader = SpinCifReader(text, symprec=1e-6) + assert np.allclose(reader.spin_basis, np.eye(3), atol=1e-6, rtol=0) + assert np.allclose(reader.get_structure().site_properties["magmom"], [[0, 0, 3], [0, 0, 3]]) + + +@pytest.mark.parametrize("mixed", [False, True]) +def test_reader_occupancy_uses_symprec_without_normalizing(external_scif, mixed): + text = external_scif.replace( + "_atom_site_symmetry_multiplicity\n", + "_atom_site_symmetry_multiplicity\n_atom_site_occupancy\n", + ).replace( + "Mn1 Mn 0.25 0.0 0.0 2", + "Mn1 Mn 0.25 0.0 0.0 2 0.5\nFe1 Fe 0.25 0.0 0.0 2 0.500000005" + if mixed + else "Mn1 Mn 0.25 0.0 0.0 2 1.000000005", + ) + with pytest.raises(ValueError, match="occupanc"): + SpinCifReader(text, symprec=1e-10).get_structure() + structure = SpinCifReader(text, symprec=1e-6).get_structure() + expected = {"Mn": 0.5, "Fe": 0.500000005} if mixed else {"Mn": 1.000000005} + assert len(structure) == 2 + assert all(site.species.as_dict() == expected for site in structure) + + +def test_reader_occupancy_respects_structure_limit(external_scif): + occupancy = 0.5 + 2 * Composition.amount_tolerance + text = external_scif.replace( + "_atom_site_symmetry_multiplicity\n", + "_atom_site_symmetry_multiplicity\n_atom_site_occupancy\n", + ).replace( + "Mn1 Mn 0.25 0.0 0.0 2", + f"Mn1 Mn 0.25 0.0 0.0 2 0.5\nFe1 Fe 0.25 0.0 0.0 2 {occupancy}", + ) + with pytest.raises(ValueError, match="occupanc"): + SpinCifReader(text, symprec=1e-2).get_structure() + + +def test_reader_orbit_positions_use_absolute_symprec(external_scif): + text = external_scif.replace("Mn1 Mn 0.25", "Mn1 Mn 0.000002") + structure = SpinCifReader(text, symprec=1e-7).get_structure() + assert sorted(site.frac_coords[0] for site in structure) == [0.000002, 0.999998] + with pytest.raises(ValueError, match="multiplicity"): + SpinCifReader(text, symprec=1e-3).get_structure() + + +def test_reader_constituent_positions_use_absolute_symprec(external_scif): + text = external_scif.replace("2 -x,-y,z,+1 u,v,w\n", "").replace( + "Mn1 Mn 0.25 0.0 0.0 2", + "Mn1 Mn 0.000002 0.0 0.0 1\nFe1 Fe 0.999998 0.0 0.0 1", + ) + structure = SpinCifReader(text, symprec=1e-7).get_structure() + assert {site.species_string: site.frac_coords[0] for site in structure} == { + "Mn": 0.000002, + "Fe": 0.999998, + } + with pytest.raises(ValueError, match="occupanc"): + SpinCifReader(text, symprec=1e-3).get_structure() + + def test_reader_legacy_cartn_frame_fallback(enumerated): # Files written before the P_abc switch declare the frame as # spinframe_orientation_cartn [0 0 0]; both routes must decode to the diff --git a/tests/test_mcif.py b/tests/test_mcif.py index bb9c99f..747adcd 100644 --- a/tests/test_mcif.py +++ b/tests/test_mcif.py @@ -121,6 +121,24 @@ def test_mcif_writer_str_matches_write_file(mcif_path: str, tmp_path) -> None: assert out.read_text() == str(writer) +@pytest.mark.parametrize("symprec", [1e-8, 3e-3, 1e-2]) +def test_translation_snapping_respects_cartesian_symprec(symprec) -> None: + structure = Structure( + Lattice.from_parameters(4, 5, 6, 70, 80, 75), + ["Fe"], + [[0.0002, 0.0002, 0.0002]], + site_properties={"magmom": [[0, 0, 2]]}, + ) + text = str(MCifWriter(structure, symprec=symprec)) + (block,) = CifFile.from_str(text).data.values() + operations = CifParser.from_str(text).get_magsymops(block) + (site,) = structure + for operation in operations: + displacement = operation.operate(site.frac_coords) - site.frac_coords + displacement -= np.rint(displacement) + assert np.linalg.norm(displacement @ structure.lattice.matrix) <= symprec + + def test_save_requires_magmom_property() -> None: from pymatgen.core import Lattice, Structure From 7d1e0bce981c0fc0f97bf3f923094803e65baf95 Mon Sep 17 00:00:00 2001 From: lan496 Date: Tue, 13 Oct 2026 07:56:30 +0900 Subject: [PATCH 3/3] fix: use constituent structures for lossless CIF exports Replace parallel chemistry arguments with Structure metadata, preserve reader rows, and honor caller tolerances through symmetry and serialization. Co-Authored-By: Codex --- CHANGELOG.md | 12 +- docs/export-files.md | 43 ++++- docs/file-formats.md | 41 +++-- src/spinforge/_cif.py | 214 ++++++++++++++++------ src/spinforge/mcif.py | 78 ++++---- src/spinforge/scif/__init__.py | 15 +- src/spinforge/scif/_format.py | 28 +-- src/spinforge/scif/_reader.py | 171 +++++++++--------- src/spinforge/scif/_writer.py | 222 +++++++++++++++-------- tests/scif_tests/test_compositions.py | 249 ++++++++++++-------------- tests/scif_tests/test_reader.py | 122 +++++++++++-- tests/scif_tests/test_structure.py | 183 +++++++++++++++++++ tests/scif_tests/test_writer.py | 36 +++- tests/test_mcif.py | 152 ++++++++++++++-- 14 files changed, 1098 insertions(+), 468 deletions(-) create mode 100644 tests/scif_tests/test_structure.py diff --git a/CHANGELOG.md b/CHANGELOG.md index aae60a7..cae7cff 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,9 +20,15 @@ The initial public release includes: - Preserve mixed-site compositions, partial occupancies, oxidation states, and selected magnetic constituents in MCIF and spinCIF exports. -- Accept explicit orthogonal spin frames in `SpinCifWriter`. -- Restore occupancies when reading spinCIF and retain constituent moments when - merging disordered sites. +- Add `SpinCifWriter.from_structure` for complete constituent chemistry and + explicit orthogonal spin frames; preserve chemistry in `from_oriented`. +- Preserve separate spinCIF constituent rows with occupancies, individual + moments, and explicit magnetic declarations. +- Preserve fractional oxidation identities and weak moments without display + rounding. Reject loss during pymatgen conversion in the spinCIF reader and + Structure-based exports. +- Use caller tolerances for position, moment, and frame checks and MCIF + translation snapping. ## v1.0.0 - 2026/09/06 diff --git a/docs/export-files.md b/docs/export-files.md index f15e9aa..eddcbab 100644 --- a/docs/export-files.md +++ b/docs/export-files.md @@ -31,9 +31,11 @@ spin_cif = SpinCifWriter.from_oriented( spin_cif.write_file("candidate.scif") ``` -Do not replace this with the plain `SpinCifWriter(...)` constructor for an -oriented result; that constructor expects magnetic moments in the original -enumeration frame. +For an orientation obtained without an MSG object, use +`SpinCifWriter.from_structure(..., spin_rotation=rotation, symprec=symprec)`. +The structure must already contain the rotated moments. The plain array +constructor assumes the original enumeration frame unless `spin_rotation` is +supplied. ## Write MCIF @@ -49,6 +51,31 @@ MCifWriter(magnetic_structure).write_file("candidate.mcif") Pass `symprec` or `mag_symprec` to `MCifWriter` when the default symmetry-search tolerances are unsuitable for the structure. +## Represent mixed magnetic sites + +Use separate constituent rows when only some species at a position carry a +magnetic declaration: + +```python +from pymatgen.core import Lattice, Species, Structure + +mixed_structure = Structure( + Lattice.cubic(4), + [{Species("Fe", 2): 0.6}, {Species("Zn", 2): 0.2}], + [[0, 0, 0], [0, 0, 0]], + site_properties={ + "magmom": [[0, 0, 2], [0, 0, 0]], + "magnetic": [True, False], + }, +) +MCifWriter(mixed_structure).write_file("mixed.mcif") +``` + +The vacancy fraction remains 0.2. Fe's written moment is 2, and its contribution +to the mean-site moment is 1.2. Set `magnetic=True` with a zero vector to preserve +an explicitly magnetic zero. For spinCIF, pass the same structure to +`SpinCifWriter.from_structure` together with its SSA model and orientation. + ## Check the files Parse the spinCIF with SpinForge and expand its asymmetric unit: @@ -59,7 +86,12 @@ from spinforge.scif import SpinCifReader spin_structure = SpinCifReader.from_file("candidate.scif").get_structure() ``` -Parse the MCIF with pymatgen: +Each returned row retains its constituent occupancy, Cartesian `magmom`, and +`magnetic` declaration. Distinct constituents can occupy the same position and +have different moments. The current writers require a shared magnetic vector +when grouping such rows for their symmetry backends. + +Parse single-species-site MCIFs with pymatgen: ```python from pymatgen.core import Structure @@ -67,6 +99,9 @@ from pymatgen.core import Structure magnetic_structure_from_file = Structure.from_file("candidate.mcif") ``` +Pymatgen's magnetic structure parser does not support mixed-site MCIFs. Those +files need a consumer that preserves their constituent rows and occupancies. + A successful parse catches syntax and expansion errors. Add application-level checks for composition, lattice, positions, and moment vectors before using a round-tripped structure in a calculation. diff --git a/docs/file-formats.md b/docs/file-formats.md index 171b09c..eeb5434 100644 --- a/docs/file-formats.md +++ b/docs/file-formats.md @@ -32,21 +32,40 @@ you need both representations. file. The reader supports files written by SpinForge and the FINDSPINGROUP reference files; it is not a general-purpose spinCIF parser. - Both writers preserve mixed species, oxidation states, and partial occupancy. - Use `magnetic_species` to select constituents that share each site's - occupied-atom moment. Other constituents have zero moment. Occupancy is written - separately and is not multiplied into the moment. -- For spinCIF, pass full `site_compositions` when atomic numbers do not describe - the chemistry. `spin_rotation` accepts the explicit orthogonal matrix that - maps intrinsic spins to Cartesian spins, including an improper frame. Supply - moments already rotated by that matrix. -- `SpinCifReader.get_structure()` merges coincident constituents and retains - `species_magmoms` and `magnetic_species` site properties. Its `magmom` remains - the shared selected occupied-atom vector. It rejects unequal selected moments - at one position rather than averaging them. + Their `Structure` input carries Cartesian occupied-atom `magmom` values and + an optional boolean `magnetic` property, which defaults to true. False omits a + magnetic declaration; true preserves it even at zero moment. To assign a + moment to only part of a mixed site, supply separate constituent rows at the + same position. Occupancy is written separately from the moment. +- Use `SpinCifWriter.from_structure` for complete chemistry and an explicit + `spin_rotation`. The matrix maps intrinsic spins to Cartesian spins, including + improper frames. Supply moments already rotated by that matrix. +- The writers group coincident rows for their symmetry backends. Magnetic + constituents at one physical site must share a vector. The spinCIF reader + preserves separate constituent rows and permits unequal coincident vectors. + Its output therefore needs application-specific checks before use by routines + that require one site per position. +- Occupancy must survive pymatgen's representation. The reader reports an error + if conversion would discard a constituent; it does not normalize occupancies. - Pymatgen's structure parser does not support disordered magnetic CIFs. Mixed MCIF output therefore requires a consumer that preserves constituent rows and their individual occupancies and moments. +## Numerical tolerances + +`MCifWriter` and `SpinCifWriter.from_structure` use `symprec` in Cartesian length +units. MCIF also uses it to bound displacement from rational translation snapping. +An explicit `mag_symprec` controls moment comparisons in moment units, while +spinCIF's `atol` controls dimensionless frame and coefficient checks. + +The reader's `symprec` retains its fractional-position units. Its optional +`mag_symprec` and `atol` override moment and normalized-frame tolerances; +omitting them retains the existing `symprec` fallback. Absolute comparisons use +`rtol=0`. For the reader and Structure-based exports, chemical identities and +occupancies are independent of these geometric and magnetic tolerances. The +legacy `site_occupancies` argument retains its `atol`-based occupancy bounds and +orbit comparisons; it writes the representative occupancy for each orbit. + To write and check either format, continue to [Export result files](export-files.md). For the scientific relation between spin and magnetic symmetry, see the [SpinForge article](https://doi.org/10.1103/8n3w-h2t1). diff --git a/src/spinforge/_cif.py b/src/spinforge/_cif.py index 7297378..ed023a9 100644 --- a/src/spinforge/_cif.py +++ b/src/spinforge/_cif.py @@ -2,11 +2,11 @@ from collections.abc import Sequence from dataclasses import dataclass +from math import fsum import numpy as np -from pymatgen.core import Composition, DummySpecies, Element, Species -from pymatgen.core.periodic_table import get_el_sp -from spgrep.utils import NDArrayFloat +from pymatgen.core import Composition, DummySpecies, Element, Species, Structure +from pymatgen.electronic_structure.core import Magmom _ATOM_SITE_LOOP_HEADER = ( "", @@ -27,68 +27,166 @@ class _SiteConstituent: magnetic: bool -def _normalize_site_species( - site_compositions: Sequence[Composition], - magnetic_species: Sequence[Sequence[str]] | None, - magnetic_moments: NDArrayFloat, -) -> tuple[tuple[_SiteConstituent, ...], ...]: - """Validate chemistry and shared occupied-atom moments for CIF writers.""" - num_sites = len(site_compositions) - moments = np.asarray(magnetic_moments, dtype=float) - if moments.shape != (num_sites, 3) or not np.all(np.isfinite(moments)): - raise ValueError(f"magnetic_moments must be a finite ({num_sites}, 3) array") - if magnetic_species is not None and len(magnetic_species) != num_sites: - raise ValueError("magnetic_species must have one selection per site") - - normalized = [] - for index, composition in enumerate(site_compositions): - items = sorted(composition.items(), key=lambda item: str(item[0])) - if not items: - raise ValueError(f"Site {index} must contain at least one species") - for species, occupancy in items: - if not isinstance(species, (Element, Species)) or isinstance(species, DummySpecies): - raise TypeError(f"Unsupported species {species!r} at site {index}") - if isinstance(species, Species) and ( - species.spin is not None - or species.oxi_state is None - or not np.isfinite(species.oxi_state) - ): - raise ValueError(f"Unsupported species {species!r} at site {index}") - if not np.isfinite(occupancy) or occupancy <= 0: - raise ValueError(f"Site {index} occupancies must be finite and positive") - if sum(occupancy for _, occupancy in items) > 1 + 1e-8: - raise ValueError(f"Site {index} total occupancy exceeds 1") - - identities = {str(species) for species, _ in items} - if magnetic_species is None: - selected = identities +@dataclass(frozen=True) +class _PhysicalSites: + structure: Structure + constituents: tuple[tuple[_SiteConstituent, ...], ...] + + +def _number(value: float) -> str: + """A positional decimal that round-trips to the same float.""" + return np.format_float_positional(float(value), unique=True, trim="-") + + +def _xyzt_string(rotation: np.ndarray, translation: np.ndarray, time_reversal: int) -> str: + """Serialize integral rotations and float translations without rational approximation.""" + if not np.array_equal(rotation, np.rint(rotation)): + raise ValueError("Spatial rotation must have integer coefficients") + expressions = [] + for row, offset in zip(rotation, translation, strict=True): + expression = "" + for coefficient, variable in zip(row, "xyz", strict=True): + if coefficient == 0: + continue + sign = "-" if coefficient < 0 else "+" if expression else "" + magnitude = abs(int(coefficient)) + expression += sign + (str(magnitude) if magnitude != 1 else "") + variable + if offset != 0: + expression += ("+" if expression and offset > 0 else "") + _number(offset) + expressions.append(expression or "0") + return ",".join(expressions) + f",{time_reversal:+d}" + + +def _species_token(species: Element | Species) -> str: + if isinstance(species, Species): + charge = _oxidation_state(species) + return f"{species.symbol}{_number(abs(charge))}{'+' if charge >= 0 else '-'}" + return species.symbol + + +def _species_key(species: Element | Species) -> tuple[str, bool, float]: + return ( + species.symbol, + isinstance(species, Species), + _oxidation_state(species) if isinstance(species, Species) else 0.0, + ) + + +def _oxidation_state(species: Species) -> float: + charge = species.oxi_state + if charge is None or not np.isfinite(charge): + raise ValueError(f"Unsupported species {species!r}") + return float(charge) + + +def _cartesian_moment(moment) -> np.ndarray: + """Convert Magmom's spin axis without its small-component zero flushing.""" + if not isinstance(moment, Magmom): + value = np.asarray(moment, dtype=float) + return np.array([0.0, 0.0, float(value)]) if value.ndim == 0 else value + x, y, z = moment.saxis + alpha = np.arctan2(y, x) + beta = np.arctan2(np.hypot(x, y), z) + ca, sa, cb, sb = np.cos(alpha), np.sin(alpha), np.cos(beta), np.sin(beta) + return np.asarray(moment.moment) @ np.array( + [[cb * ca, cb * sa, -sb], [-sa, ca, 0.0], [sb * ca, sb * sa, cb]] + ) + + +def _physical_sites(structure: Structure, *, symprec: float, mag_symprec: float) -> _PhysicalSites: + """Adapt constituent rows to the one-vector-per-site writer backends.""" + if "magmom" not in structure.site_properties: + raise ValueError("structure must have a 'magmom' site property") + moments = np.asarray( + [_cartesian_moment(moment) for moment in structure.site_properties["magmom"]] + ) + if moments.shape != (len(structure), 3) or not np.all(np.isfinite(moments)): + raise ValueError("magmom must contain finite Cartesian 3-vectors") + magnetic = structure.site_properties.get("magnetic", [True] * len(structure)) + if any(not isinstance(flag, (bool, np.bool_)) for flag in magnetic): + raise ValueError("magnetic must contain boolean flags") + if any(not flag and np.any(moment != 0) for flag, moment in zip(magnetic, moments)): + raise ValueError("A nonmagnetic row must have a zero moment") + + distances = structure.lattice.get_all_distances(structure.frac_coords, structure.frac_coords) + groups: list[list[int]] = [] + for index in range(len(structure)): + matches = [group for group in groups if np.any(distances[index, group] <= symprec)] + if len(matches) > 1 or (matches and np.any(distances[index, matches[0]] > symprec)): + raise ValueError("Ambiguous coincident-site grouping within symprec") + if matches: + matches[0].append(index) else: - selection = magnetic_species[index] - if isinstance(selection, str) or any(not isinstance(name, str) for name in selection): - raise ValueError(f"Site {index} magnetic_species must be a sequence of species") - selected = {str(get_el_sp(name)) for name in selection} - if len(selected) != len(selection): - raise ValueError(f"Site {index} magnetic_species contains duplicate species") - if not selected <= identities: - raise ValueError(f"Site {index} magnetic_species contains an absent species") - if not selected and np.any(moments[index] != 0): - raise ValueError(f"Site {index} has a nonzero moment but no selected magnetic species") - normalized.append( + groups.append([index]) + + compositions = [] + positions = [] + shared_moments = [] + constituents = [] + for group in groups: + selected = [moments[index] for index in group if magnetic[index]] + moment = selected[0] if selected else np.zeros(3) + if any( + np.linalg.norm(first - second) > mag_symprec + for first in selected + for second in selected + ): + raise ValueError("Unequal magnetic moments at one physical site are unsupported") + row_amounts: dict[tuple[Element | Species, bool], list[float]] = {} + amounts: dict[Element | Species, list[float]] = {} + for index in group: + composition = structure[index].species + if not composition: + raise ValueError("Every row must contain at least one species") + for species, occupancy in composition.items(): + if not isinstance(species, (Element, Species)) or isinstance( + species, DummySpecies + ): + raise TypeError(f"Unsupported species {species!r}") + if isinstance(species, Species): + if species.spin is not None: + raise ValueError(f"Unsupported species {species!r}") + _oxidation_state(species) + if not np.isfinite(occupancy) or occupancy <= 0: + raise ValueError("Occupancies must be finite and positive") + row_amounts.setdefault((species, bool(magnetic[index])), []).append(occupancy) + amounts.setdefault(species, []).append(occupancy) + expected = {species: fsum(values) for species, values in amounts.items()} + composition = Composition(expected) + if dict(composition.items()) != expected: + raise ValueError("The site composition is not representable by pymatgen") + compositions.append(composition) + positions.append(structure[group[0]].frac_coords) + shared_moments.append(moment) + rows = [ + _SiteConstituent(species, fsum(values), flag) + for (species, flag), values in row_amounts.items() + ] + constituents.append( tuple( - _SiteConstituent(species, float(occupancy), str(species) in selected) - for species, occupancy in items + sorted( + rows, key=lambda row: (_species_key(row.species), row.occupancy, row.magnetic) + ) ) ) - return tuple(normalized) + physical = Structure( + structure.lattice, + compositions, + positions, + site_properties={"magmom": shared_moments}, + ) + if any( + dict(site.species.items()) != dict(composition.items()) + for site, composition in zip(physical, compositions, strict=True) + ): + raise ValueError("The site composition is not representable by pymatgen") + return _PhysicalSites(physical, tuple(constituents)) def _atom_type_lines(site_species: Sequence[Sequence[_SiteConstituent]]) -> list[str]: """Declare ionic atom types while retaining uncharged Element identities.""" ionic_species = { - str(row.species): row.species - for site in site_species - for row in site - if isinstance(row.species, Species) + row.species for site in site_species for row in site if isinstance(row.species, Species) } if not ionic_species: return [] @@ -98,7 +196,7 @@ def _atom_type_lines(site_species: Sequence[Sequence[_SiteConstituent]]) -> list "_atom_type_symbol", "_atom_type_oxidation_number", *[ - f"{symbol} {species.oxi_state:.12g}" - for symbol, species in sorted(ionic_species.items()) + f"{_species_token(species)} {_number(_oxidation_state(species))}" + for species in sorted(ionic_species, key=_species_key) ], ] diff --git a/src/spinforge/mcif.py b/src/spinforge/mcif.py index 9c5b5e6..2804c88 100644 --- a/src/spinforge/mcif.py +++ b/src/spinforge/mcif.py @@ -16,8 +16,6 @@ from __future__ import annotations -from collections.abc import Sequence - import numpy as np from moyopy import ( MagneticSpaceGroupType, @@ -25,14 +23,15 @@ NonCollinearMagneticCell, ) from pymatgen.core import Structure -from pymatgen.core.operations import MagSymmOp -from pymatgen.electronic_structure.core import Magmom from spinforge._cif import ( _ATOM_SITE_LOOP_HEADER, _atom_type_lines, - _normalize_site_species, + _number, + _physical_sites, _SiteConstituent, + _species_token, + _xyzt_string, ) __all__ = [ @@ -61,10 +60,15 @@ class MCifWriter: ``CifParser`` (``m = mx*a_hat + my*b_hat + mz*c_hat``). Each constituent is written as an atom row with its occupancy and full species - identity, including oxidation state. Compositions and magnetic selections enter + identity, including oxidation state. Compositions and magnetic flags enter the moyopy equivalence labels, so symmetry cannot exchange chemically distinct sites. Selected constituents share the site's occupied-atom moment; unselected constituents have zero moment. Occupancy never scales the written moment. + The optional boolean ``magnetic`` site property defaults to true. A false + flag requires a zero moment; a true flag retains a moment row even for zero. + A mixed site applies its moment and flag to all its species. Use coincident + constituent rows for partial magnetic selection. Their selected moments must + agree within ``mag_symprec`` (which defaults to ``symprec`` as in moyopy). Parameters ---------- @@ -72,10 +76,8 @@ class MCifWriter: Magnetic structure with a ``magmom`` site property. symprec, mag_symprec moyopy symmetry-search tolerances (basis units / moment units). - magnetic_species - Species identities selected at each site, such as ``[("Fe2+",), ()]``. - Defaults to all constituents. An empty selection requires a zero moment. - Selected constituents retain moment rows even when their moment is zero. + ``symprec`` also bounds the Cartesian displacement introduced when + snapping translations to the crystallographic fraction grid. """ def __init__( @@ -84,15 +86,14 @@ def __init__( symprec: float = 1e-4, *, mag_symprec: float | None = None, - magnetic_species: Sequence[Sequence[str]] | None = None, ) -> None: - if "magmom" not in struct.site_properties: - raise ValueError("structure must have a 'magmom' site property to write an mcif") - basis = struct.lattice.matrix - moments = np.array([Magmom(m).global_moment for m in struct.site_properties["magmom"]]) - site_species = _normalize_site_species( - [site.species for site in struct], magnetic_species, moments + physical = _physical_sites( + struct, symprec=symprec, mag_symprec=symprec if mag_symprec is None else mag_symprec ) + struct = physical.structure + basis = struct.lattice.matrix + moments = np.asarray(struct.site_properties["magmom"]) + site_species = physical.constituents cell = NonCollinearMagneticCell( basis=basis.tolist(), positions=struct.frac_coords.tolist(), @@ -104,13 +105,11 @@ def __init__( ) ops = dataset.magnetic_operations op_strings = [ - MagSymmOp.from_rotation_and_translation_and_time_reversal( - rotation_matrix=np.array(rot), - translation_vec=self._snap_translation(trans), - time_reversal=-1 if rev else 1, + _xyzt_string( + np.array(rot), + self._snap_translation(trans, basis, symprec=symprec), + -1 if rev else 1, ) - .as_xyzt_str() - .replace(" ", "") # CIF loop columns are whitespace-delimited; xyzt must be one token for rot, trans, rev in zip(ops.rotations, ops.translations, ops.time_reversals) ] @@ -150,12 +149,12 @@ def _format_mcif( "data_magnetic_structure", f"_chemical_formula_sum '{structure.composition.element_composition.formula}'", f"_space_group_magn.number_BNS {bns}", - f"_cell_length_a {a:.6f}", - f"_cell_length_b {b:.6f}", - f"_cell_length_c {c:.6f}", - f"_cell_angle_alpha {alpha:.6f}", - f"_cell_angle_beta {beta:.6f}", - f"_cell_angle_gamma {gamma:.6f}", + f"_cell_length_a {_number(a)}", + f"_cell_length_b {_number(b)}", + f"_cell_length_c {_number(c)}", + f"_cell_angle_alpha {_number(alpha)}", + f"_cell_angle_beta {_number(beta)}", + f"_cell_angle_gamma {_number(gamma)}", *_atom_type_lines(site_species), "", "loop_", @@ -181,11 +180,11 @@ def _format_mcif( counts[symbol] = counts.get(symbol, 0) + 1 label = f"{symbol}{counts[symbol]}" lines.append( - f"{label} {constituent.species} {fx:.6f} {fy:.6f} {fz:.6f} " - f"{constituent.occupancy:.12g}" + f"{label} {_species_token(constituent.species)} " + f"{_number(fx)} {_number(fy)} {_number(fz)} {_number(constituent.occupancy)}" ) if constituent.magnetic: - moment_lines.append(f"{label} {comps[0]:.6f} {comps[1]:.6f} {comps[2]:.6f}") + moment_lines.append(f"{label} {' '.join(_number(value) for value in comps)}") if moment_lines: lines += [ "", @@ -205,21 +204,20 @@ def _composition_numbers(site_species: tuple[tuple[_SiteConstituent, ...], ...]) signatures: dict[tuple, int] = {} numbers = [] for site in site_species: - signature = tuple((str(row.species), row.occupancy, row.magnetic) for row in site) + signature = tuple((row.species, row.occupancy, row.magnetic) for row in site) numbers.append(signatures.setdefault(signature, len(signatures) + 1)) return numbers @staticmethod def _snap_translation( - trans: list[float], *, denom: int = 24, atol: float = 1e-3 + trans: list[float], basis: np.ndarray, *, symprec: float, denom: int = 24 ) -> np.ndarray: - """Snap a fractional translation to the nearest k/denom when within atol. + """Snap a fractional translation when its Cartesian change is within symprec. - moyopy's tolerant symmetry search can return translations carrying ~1e-6 noise - on a fitted (numerically imperfect) structure. The crystallographic translations - here are multiples of 1/2, 1/3, 1/4, 1/6 (all on the 1/24 grid), so snapping - removes the noise without altering a genuine translation. + The 1/24 grid includes usual crystallographic fractions. General-origin + translations need not lie on that grid, so preserve the complete vector + unless snapping it respects the caller's Cartesian symmetry tolerance. """ t = np.array(trans, dtype=float) % 1.0 snapped = np.round(t * denom) / denom - return np.where(np.abs(snapped - t) <= atol, snapped, t) + return snapped if np.linalg.norm((snapped - t) @ basis) <= symprec else t diff --git a/src/spinforge/scif/__init__.py b/src/spinforge/scif/__init__.py index 366b740..6511162 100644 --- a/src/spinforge/scif/__init__.py +++ b/src/spinforge/scif/__init__.py @@ -22,14 +22,13 @@ it enters through ``_space_group_spin.collinear_direction_xyz`` / ``_space_group_spin.coplanar_perp_uvw``, filled from the ``spinspg.spin.SpinOnlyGroup`` the structure was enumerated with. -- ``spin_rotation`` accepts the explicit orthogonal frame used for the sampled - moments, independently of an MSG object. ``site_compositions`` preserves - partial and mixed occupancies and oxidation states; ``magnetic_species`` - assigns the shared occupied-atom moment only to selected constituents. -- The reader merges constituent rows into disordered sites. Its ``magmom`` - property remains the shared occupied-atom vector, while ``species_magmoms`` - and ``magnetic_species`` retain the constituent values and selections. Compute - a mean-site moment by weighting constituent moments by occupancy once. +- ``SpinCifWriter.from_structure`` preserves constituent chemistry and accepts + the explicit ``spin_rotation`` used for sampled moments, independently of + an MSG object. The Structure's ``magnetic`` flag distinguishes an absent + magnetic declaration from an explicitly zero ``magmom``. +- The reader preserves separate constituent rows, including coincident rows + with different moments. Occupancies and occupied-atom moments remain separate; + a mean-site moment requires weighting constituent moments by occupancy once. The dictionary is explicitly preliminary upstream ("all dictionary aspects and file formats are subject to change"): the targeted revision is pinned as diff --git a/src/spinforge/scif/_format.py b/src/spinforge/scif/_format.py index 3b59e32..6f07594 100644 --- a/src/spinforge/scif/_format.py +++ b/src/spinforge/scif/_format.py @@ -4,20 +4,10 @@ from __future__ import annotations import numpy as np -from pymatgen.core.operations import MagSymmOp from spgrep.utils import NDArrayFloat from spinspg.spin import SpinOnlyGroup, SpinOnlyGroupType - -def _fmt(value: float, *, decimals: int = 6) -> str: - """Format with fixed decimals, avoiding a ``-0.000000`` output. - - The zero-flush threshold is half a unit in the last printed place, i.e. - exactly the rounding boundary of the format itself, so no tolerance beyond - the output precision is introduced. - """ - v = 0.0 if abs(value) < 0.5 * 10.0**-decimals else float(value) - return f"{v:.{decimals}f}" +from spinforge._cif import _number, _xyzt_string def _spinframe_p_abc(frame: NDArrayFloat, lattice: NDArrayFloat, *, atol: float) -> str: @@ -52,7 +42,7 @@ def _collinear_direction_xyz( axis = q_mat @ np.asarray(spin_only_group.axis, dtype=float) direction = axis @ np.linalg.inv(lattice) direction = direction / np.max(np.abs(direction)) - return "'" + ",".join(_fmt(x) for x in direction) + "'" + return "'" + ",".join(_number(x) for x in direction) + "'" def _coplanar_perp_uvw( @@ -66,22 +56,16 @@ def _coplanar_perp_uvw( if spin_only_group.spin_only_group_type != SpinOnlyGroupType.COPLANAR: return "." normal = frame.T @ q_mat @ np.asarray(spin_only_group.axis, dtype=float) - return "'" + ",".join(_fmt(x) for x in normal) + "'" + return "'" + ",".join(_number(x) for x in normal) + "'" def _xyzt(w: np.ndarray, t: np.ndarray, u: np.ndarray) -> str: time_reversal = -1 if np.linalg.det(u) < 0 else 1 - op = MagSymmOp.from_rotation_and_translation_and_time_reversal( - rotation_matrix=np.array(w), - translation_vec=np.remainder(t, 1.0), - time_reversal=time_reversal, - ) - # CIF loop columns are whitespace-delimited; xyzt must be one token - return op.as_xyzt_str().replace(" ", "") + return _xyzt_string(w, np.remainder(t, 1.0), time_reversal) # Exact coefficient strings for common rational values; anything else -# (e.g. sqrt(3)/2) is written as a plain 15-significant-digit number to +# (e.g. sqrt(3)/2) is written as a lossless decimal to # keep the expression grammar free of math functions. _EXACT_COEFFICIENTS = [ (1.0, "1"), @@ -101,7 +85,7 @@ def _coefficient_str(value: float, *, atol: float) -> str: for ref, text in _EXACT_COEFFICIENTS: if np.isclose(magnitude, ref, rtol=0, atol=atol): return sign + text - return f"{value:.15g}" + return _number(value) def _linear_expr(coefficients: NDArrayFloat, variables: tuple[str, ...], *, atol: float) -> str: diff --git a/src/spinforge/scif/_reader.py b/src/spinforge/scif/_reader.py index 0105f58..43c5853 100644 --- a/src/spinforge/scif/_reader.py +++ b/src/spinforge/scif/_reader.py @@ -64,19 +64,35 @@ class SpinCifReader: text The spinCIF content. symprec - Tolerance (fractional units) for merging symmetry-equivalent - positions and for moment consistency when expanding the orbits. + Absolute fractional-position tolerance for expanding each constituent's + orbit. This also supplies the default for ``mag_symprec`` and ``atol``. + mag_symprec + Absolute occupied-atom moment tolerance for consistency within an orbit. + Defaults to ``symprec``. + atol + Absolute tolerance for orthogonality of the normalized spin frame. + Defaults to ``symprec``. All comparisons use zero relative tolerance. """ - def __init__(self, text: str, *, symprec: float = 1e-4) -> None: + def __init__( + self, + text: str, + *, + symprec: float = 1e-4, + mag_symprec: float | None = None, + atol: float | None = None, + ) -> None: self._symprec = symprec + self._mag_symprec = symprec if mag_symprec is None else mag_symprec blocks = CifFile.from_str(text).data if len(blocks) != 1: raise ValueError(f"Expected a single data block, found {len(blocks)}") (block,) = blocks.values() self.operations = _read_symop_loop(block, "_space_group_symop_spin_operation") - if "_space_group_symop_spin_lattice.xyzt" in block.data: + if any( + f"_space_group_symop_spin_lattice.{column}" in block.data for column in ("xyzt", "uvw") + ): self.spin_lattice = _read_symop_loop(block, "_space_group_symop_spin_lattice") else: # the loop is optional in the dictionary self.spin_lattice = [(np.eye(3), np.zeros(3), 1, np.eye(3))] @@ -107,6 +123,20 @@ def __init__(self, text: str, *, symprec: float = 1e-4) -> None: for label, expression in zip(*subspace_columns, strict=True) if expression not in (".", "?") } + site_columns = _read_loop_columns( + block, + tuple( + f"_atom_site_{column}" + for column in ( + "label", + "type_symbol", + "fract_x", + "fract_y", + "fract_z", + "symmetry_multiplicity", + ) + ), + ) self.sites = [ ( label, @@ -115,16 +145,9 @@ def __init__(self, text: str, *, symprec: float = 1e-4) -> None: moments.get(label, np.zeros(3)), int(multiplicity), ) - for label, symbol, fx, fy, fz, multiplicity in zip( - block.data["_atom_site_label"], - block.data["_atom_site_type_symbol"], - block.data["_atom_site_fract_x"], - block.data["_atom_site_fract_y"], - block.data["_atom_site_fract_z"], - block.data["_atom_site_symmetry_multiplicity"], - ) + for label, symbol, fx, fy, fz, multiplicity in zip(*site_columns, strict=True) ] - labels = block.data["_atom_site_label"] + labels = site_columns[0] if len(set(labels)) != len(labels): raise ValueError("Atom-site labels must be unique") if not self._magnetic_labels.issubset(labels): @@ -136,10 +159,10 @@ def __init__(self, text: str, *, symprec: float = 1e-4) -> None: ) } if any( - not np.isfinite(occupancy) or not 0 < occupancy <= 1 + not np.isfinite(occupancy) or occupancy <= 0 for occupancy in self.site_occupancies.values() ): - raise ValueError("Atom-site occupancies must be finite and lie in (0, 1]") + raise ValueError("Atom-site occupancies must be finite and positive") oxidation_states = _read_oxidation_states(block) self._site_species = {} for label, symbol, _, moment, _ in self.sites: @@ -157,14 +180,23 @@ def __init__(self, text: str, *, symprec: float = 1e-4) -> None: *(float(block.data[f"_cell_length_{name}"]) for name in "abc"), *(float(block.data[f"_cell_angle_{name}"]) for name in ("alpha", "beta", "gamma")), ) - self.spin_basis = _read_spin_basis(block, text, np.array(self.lattice.matrix)) + self.spin_basis = _read_spin_basis( + block, text, np.array(self.lattice.matrix), atol=symprec if atol is None else atol + ) self.collinear_direction_xyz = block.data["_space_group_spin.collinear_direction_xyz"] self.coplanar_perp_uvw = block.data["_space_group_spin.coplanar_perp_uvw"] @classmethod - def from_file(cls, filename: str | PathLike[str], *, symprec: float = 1e-4) -> SpinCifReader: + def from_file( + cls, + filename: str | PathLike[str], + *, + symprec: float = 1e-4, + mag_symprec: float | None = None, + atol: float | None = None, + ) -> SpinCifReader: with open(filename) as f: - return cls(f.read(), symprec=symprec) + return cls(f.read(), symprec=symprec, mag_symprec=mag_symprec, atol=atol) def expanded_operations(self) -> list[SymmetryOperation]: """Products of the spin-lattice and operation loops (the full listed group).""" @@ -177,25 +209,20 @@ def expanded_operations(self) -> list[SymmetryOperation]: def get_structure(self) -> Structure: """Expand the asymmetric unit into the full magnetic structure. - Coincident constituent rows are merged into one disordered site with - their occupancies and oxidation identities. ``species_magmoms`` maps - every constituent identity to its Cartesian occupied-atom moment; - ``magnetic_species`` lists the selected constituents, identified by - the presence of a moment row (including explicit zeros). - - The existing ``magmom`` property is the shared selected occupied-atom - vector, or zero when none are selected. Unequal selected constituent - vectors are outside this reader's shared-moment model and raise an - error. No occupancy averaging is applied. All vectors use the read - lattice's Cartesian frame. + Each expanded constituent remains a separate row with its original + species and occupancy, including coincident rows with different moments. + ``magmom`` is the Cartesian occupied-atom vector in the read lattice's + frame. ``magnetic`` records the presence of a moment row, including an + explicit zero. No occupancy averaging or normalization is applied. Raises if an orbit disagrees with its declared multiplicity, transforms - inconsistently, or repeats the same constituent at a position. + inconsistently, or pymatgen cannot preserve a constituent's species and + occupancy exactly. """ compositions: list[dict[Element | Species, float]] = [] coords: list[np.ndarray] = [] - species_magmoms: list[dict[str, np.ndarray]] = [] - selections: list[set[str]] = [] + magmoms: list[np.ndarray] = [] + magnetic: list[bool] = [] operations = self.expanded_operations() for label, symbol, frac, moment_uvw, multiplicity in self.sites: orbit: list[tuple[np.ndarray, np.ndarray]] = [] @@ -203,8 +230,8 @@ def get_structure(self) -> Structure: position = (w @ frac + t) % 1.0 moment = u @ moment_uvw for known_position, known_moment in orbit: - if is_integer_array(position - known_position, atol=self._symprec): - if not np.allclose(moment, known_moment, rtol=0, atol=self._symprec): + if is_integer_array(position - known_position, atol=self._symprec, rtol=0): + if not np.allclose(moment, known_moment, rtol=0, atol=self._mag_symprec): raise ValueError( f"Inconsistent moments for {label} at {position}: " f"{moment} vs {known_moment}" @@ -217,54 +244,26 @@ def get_structure(self) -> Structure: f"Orbit of {label} has {len(orbit)} sites, multiplicity says {multiplicity}" ) for position, moment in orbit: - index = next( - ( - index - for index, known_position in enumerate(coords) - if is_integer_array(position - known_position, atol=self._symprec) - ), - len(coords), - ) - if index == len(coords): - coords.append(position) - compositions.append({}) - species_magmoms.append({}) - selections.append(set()) - species = self._site_species[label] - identity = str(species) - if species in compositions[index]: - raise ValueError(f"Repeated constituent {identity} at {position}") - compositions[index][species] = self.site_occupancies[label] - species_magmoms[index][identity] = moment @ self.spin_basis - if label in self._magnetic_labels: - selections[index].add(identity) - - magmoms = [] - for composition, moments, selected in zip( - compositions, species_magmoms, selections, strict=True - ): - if sum(composition.values()) > 1 + 1e-6: - raise ValueError("Constituent occupancies at one site sum to more than one") - shared = moments[min(selected)] if selected else np.zeros(3) - if any( - not np.allclose(moments[identity], shared, atol=self._symprec, rtol=0) - for identity in selected - ): - raise ValueError( - "Unequal selected constituent moments cannot be represented by the shared " - "occupied-atom magmom property" - ) - magmoms.append(shared) - return Structure( + coords.append(position) + compositions.append({self._site_species[label]: self.site_occupancies[label]}) + magmoms.append(moment @ self.spin_basis) + magnetic.append(label in self._magnetic_labels) + + structure = Structure( self.lattice, compositions, coords, site_properties={ "magmom": magmoms, - "species_magmoms": species_magmoms, - "magnetic_species": [tuple(sorted(selected)) for selected in selections], + "magnetic": magnetic, }, ) + for composition, site in zip(compositions, structure, strict=True): + if dict(site.species.items()) != composition: + raise ValueError( + "Pymatgen cannot preserve atom-site species and occupancy exactly" + ) + return structure def _read_loop_columns( @@ -315,7 +314,8 @@ def _read_oxidation_states(block: CifBlock) -> dict[str, float]: def _read_symop_loop(block: CifBlock, category: str) -> list[SymmetryOperation]: ops = [] - for xyzt, uvw in zip(block.data[f"{category}.xyzt"], block.data[f"{category}.uvw"]): + columns = _read_loop_columns(block, (f"{category}.xyzt", f"{category}.uvw")) + for xyzt, uvw in zip(*columns, strict=True): op = MagSymmOp.from_xyzt_str(xyzt) ops.append( ( @@ -407,7 +407,9 @@ def _linear_form(node: ast.expr, variables: str) -> tuple[float, NDArrayFloat]: raise ValueError(f"Unsupported token in spinCIF expression: {ast.unparse(node)}") -def _read_spin_basis(block: CifBlock, text: str, lattice: np.ndarray) -> np.ndarray: +def _read_spin_basis( + block: CifBlock, text: str, lattice: np.ndarray, *, atol: float +) -> np.ndarray: """Unit spin-basis vectors (rows, Cartesian) from the frame declaration. Supports the ``transform_spinframe_P_abc`` route the writer emits (basis @@ -419,12 +421,17 @@ def _read_spin_basis(block: CifBlock, text: str, lattice: np.ndarray) -> np.ndar """ expression = block.data.get("_space_group_spin.transform_spinframe_P_abc") if expression is not None: - basis = _linear_matrix(expression, "abc") @ lattice - norms = np.linalg.norm(basis, axis=1) - if not np.all(np.isfinite(basis)) or np.any(norms < 1e-8): + coefficients = _linear_matrix(expression, "abc") + scales = np.max(np.abs(coefficients), axis=1) + if not np.all(np.isfinite(coefficients)) or np.any(scales == 0): + raise ValueError(f"Degenerate spin-frame axis in declaration: {expression}") + basis = (coefficients / scales[:, None]) @ lattice + scales = np.max(np.abs(basis), axis=1) + if not np.all(np.isfinite(basis)) or np.any(scales == 0): raise ValueError(f"Degenerate spin-frame axis in declaration: {expression}") - basis = basis / norms[:, None] - if not np.allclose(basis @ basis.T, np.eye(3), rtol=0, atol=1e-3): + basis = basis / scales[:, None] + basis = basis / np.linalg.norm(basis, axis=1)[:, None] + if not np.allclose(basis @ basis.T, np.eye(3), rtol=0, atol=atol): raise ValueError(f"Spin-frame axes are not orthogonal: {expression}") return basis match = re.search( diff --git a/src/spinforge/scif/_writer.py b/src/spinforge/scif/_writer.py index 4c0cfd0..22029f8 100644 --- a/src/spinforge/scif/_writer.py +++ b/src/spinforge/scif/_writer.py @@ -3,20 +3,21 @@ from __future__ import annotations -from collections.abc import Mapping, Sequence +from collections.abc import Sequence import numpy as np from moyopy import SpaceGroupType -from pymatgen.core import Composition, Element, Lattice, Structure -from pymatgen.electronic_structure.core import Magmom +from pymatgen.core import Element, Lattice, Structure from spgrep.utils import NDArrayFloat from spinspg.spin import SpinOnlyGroup, SpinOnlyGroupType from spinforge._cif import ( _ATOM_SITE_LOOP_HEADER, _atom_type_lines, - _normalize_site_species, + _number, + _physical_sites, _SiteConstituent, + _species_token, ) from spinforge.configuration import SpinSymmetryAdaptedStructure from spinforge.msg import MagneticSpaceSubgroup @@ -25,7 +26,6 @@ from ._format import ( _collinear_direction_xyz, _coplanar_perp_uvw, - _fmt, _spinframe_p_abc, _symmform, _uvw_expr, @@ -77,15 +77,6 @@ class SpinCifWriter: spin_rotation.T``. Both proper and improper frames are supported; ``None`` uses the identity. The operations and SSA restrictions are rotated with the same matrix. - site_compositions - Full composition of each supercell site, including fractional - occupancies and oxidation states. Defaults to the supercell's atomic - elements. Cannot be combined with ``site_occupancies``. - magnetic_species - Species identities selected at each supercell site, as strings such as - ``"Fe2+"``. Selected constituents share the site's occupied-atom moment; - unselected constituents carry zero. ``None`` selects every constituent. - Compositions and selections must be invariant under the SSG operations. site_occupancies Fractional occupancy per supercell site (``None`` = fully occupied). Written as ``_atom_site_occupancy``; must be invariant under the SSG @@ -97,8 +88,9 @@ class SpinCifWriter: ``_parent_space_group`` tags when given. atol Tolerance for spin-space and translation checks (coefficient - snapping, orthogonality). Site matching uses the supercell's own - symprec via :meth:`Supercell.site_permutations`. + snapping, orthogonality) and absolute legacy occupancy comparisons. + Site matching uses the supercell's own symprec via + :meth:`Supercell.site_permutations`. """ def __init__( @@ -109,8 +101,6 @@ def __init__( *, spin_only_group: SpinOnlyGroup, spin_rotation: NDArrayFloat | None = None, - site_compositions: Sequence[Composition | Mapping[str, float]] | None = None, - magnetic_species: Sequence[Sequence[str]] | None = None, data_name: str = "spinforge", parent_space_group_type: SpaceGroupType | None = None, site_occupancies: Sequence[float] | NDArrayFloat | None = None, @@ -125,11 +115,59 @@ def __init__( data_name=data_name, parent_space_group_type=parent_space_group_type, site_occupancies=site_occupancies, - site_compositions=site_compositions, - magnetic_species=magnetic_species, atol=atol, ) + @classmethod + def from_structure( + cls, + sas: SpinSymmetryAdaptedStructure, + nssg: NontrivialSpinSpaceGroup, + structure: Structure, + *, + spin_only_group: SpinOnlyGroup, + symprec: float, + spin_rotation: NDArrayFloat | None = None, + mag_symprec: float | None = None, + data_name: str = "spinforge", + parent_space_group_type: SpaceGroupType | None = None, + atol: float = 1e-5, + ) -> SpinCifWriter: + """Write full chemistry from constituent rows matched to the SSA sites. + + ``symprec`` bounds Cartesian position and lattice matching. Coincident + constituent rows share a physical site, and their selected moments must + agree within ``mag_symprec`` (defaulting to ``atol``). The optional + boolean ``magnetic`` property defaults to true; false requires zero. + A mixed row applies its moment and flag to all constituents. Chemistry + and magnetic flags must agree exactly across each SSG orbit. + + An explicit ``mag_symprec`` is an absolute moment tolerance for both + the SSA residual and spin-only-group checks. When omitted, these checks + retain the array constructor's moment-scale-relative ``atol`` policy. + """ + physical = _physical_sites( + structure, + symprec=symprec, + mag_symprec=atol if mag_symprec is None else mag_symprec, + ) + order = _match_ssa_sites(sas, physical.structure, symprec=symprec) + moments = np.asarray(physical.structure.site_properties["magmom"])[order] + writer = cls.__new__(cls) + writer._build( + sas, + nssg, + moments, + spin_only_group=spin_only_group, + q_mat=np.eye(3) if spin_rotation is None else np.asarray(spin_rotation, dtype=float), + data_name=data_name, + parent_space_group_type=parent_space_group_type, + constituents=tuple(physical.constituents[index] for index in order), + mag_symprec=mag_symprec, + atol=atol, + ) + return writer + @classmethod def from_oriented( cls, @@ -139,11 +177,11 @@ def from_oriented( msg: MagneticSpaceSubgroup, *, spin_only_group: SpinOnlyGroup, - site_compositions: Sequence[Composition | Mapping[str, float]] | None = None, - magnetic_species: Sequence[Sequence[str]] | None = None, data_name: str = "spinforge", parent_space_group_type: SpaceGroupType | None = None, site_occupancies: Sequence[float] | NDArrayFloat | None = None, + symprec: float = 1e-4, + mag_symprec: float | None = None, atol: float = 1e-5, ) -> SpinCifWriter: """Writer for one ``(structure, msg)`` pair from ``generate_oriented``. @@ -152,20 +190,41 @@ def from_oriented( conjugates the spin parts of the SSG operations by the MSG rotation ``msg.Q`` that ``generate_oriented`` applied to them. """ - moments = np.array( - [Magmom(m).global_moment for m in structure.site_properties["magmom"]], dtype=float - ) - return cls( + if site_occupancies is not None: + if any(not site.is_ordered for site in structure): + raise ValueError("site_occupancies can only override ordered input") + occupancies = _validated_occupancies(site_occupancies, len(structure), [], atol=atol) + physical = _physical_sites( + structure, + symprec=symprec, + mag_symprec=atol if mag_symprec is None else mag_symprec, + ) + order = _match_ssa_sites(sas, physical.structure, symprec=symprec) + writer = cls.__new__(cls) + writer._build( + sas, + nssg, + np.asarray(physical.structure.site_properties["magmom"])[order], + spin_only_group=spin_only_group, + q_mat=np.asarray(msg.Q, dtype=float), + data_name=data_name, + parent_space_group_type=parent_space_group_type, + constituents=tuple(physical.constituents[index] for index in order), + site_occupancies=occupancies[order], + mag_symprec=mag_symprec, + atol=atol, + ) + return writer + return cls.from_structure( sas, nssg, - moments, + structure, spin_only_group=spin_only_group, spin_rotation=np.asarray(msg.Q, dtype=float), + symprec=symprec, + mag_symprec=mag_symprec, data_name=data_name, parent_space_group_type=parent_space_group_type, - site_occupancies=site_occupancies, - site_compositions=site_compositions, - magnetic_species=magnetic_species, atol=atol, ) @@ -180,8 +239,8 @@ def _build( data_name: str, parent_space_group_type: SpaceGroupType | None, site_occupancies: Sequence[float] | NDArrayFloat | None = None, - site_compositions: Sequence[Composition | Mapping[str, float]] | None = None, - magnetic_species: Sequence[Sequence[str]] | None = None, + constituents: tuple[tuple[_SiteConstituent, ...], ...] | None = None, + mag_symprec: float | None = None, atol: float, ) -> None: positions = np.asarray(sas.supercell.positions) @@ -198,14 +257,13 @@ def _build( or not np.allclose(q_mat @ q_mat.T, np.eye(3), rtol=0, atol=atol) ): raise ValueError("spin_rotation must be a finite orthogonal 3 by 3 matrix") - if site_compositions is not None and site_occupancies is not None: - raise ValueError("Use either site_compositions or site_occupancies, not both") - lattice = np.asarray(sas.supercell.basis) frame = _crystal_cartesian_frame(lattice) # columns = spin basis vectors (Cartesian) oriented_basis = [q_mat @ np.asarray(b).T for b in sas.magnetic_moments_basis] - _validate_moments(moments, oriented_basis, spin_only_group, q_mat, atol=atol) + _validate_moments( + moments, oriented_basis, spin_only_group, q_mat, atol=atol, mag_symprec=mag_symprec + ) prim_operations = _nontrivial_operations( nssg.nontrivial_coset, nssg.invariant_rotations, nssg.invariant_translations @@ -226,29 +284,36 @@ def _build( [t + t2 for _, t, _ in prim_operations for t2 in spin_translations], ) orbits = _orbit_partition(len(positions), permutations) - occupancies = _validated_occupancies( - site_occupancies, len(positions), permutations, atol=atol - ) - compositions = ( - [Composition(composition) for composition in site_compositions] - if site_compositions is not None - else [ - Composition({Element.from_Z(number): occupancy}) + if constituents is None: + occupancies = _validated_occupancies( + site_occupancies, len(positions), permutations, atol=atol + ) + constituents = tuple( + (_SiteConstituent(Element.from_Z(number), occupancy, True),) for number, occupancy in zip(sas.supercell.numbers, occupancies, strict=True) - ] - ) - constituents = _normalize_site_species(compositions, magnetic_species, moments) - signatures = [ - tuple((str(c.species), c.occupancy, c.magnetic) for c in site) for site in constituents - ] - for permutation in permutations: - if any( - signatures[index] != signatures[int(image)] - for index, image in enumerate(permutation) - ): - raise ValueError( - "site_compositions and magnetic_species must be invariant under the SSG operations" + ) + else: + if site_occupancies is not None: + occupancies = _validated_occupancies( + site_occupancies, len(positions), permutations, atol=atol + ) + constituents = tuple( + tuple(_SiteConstituent(row.species, amount, row.magnetic) for row in site) + for site, amount in zip(constituents, occupancies, strict=True) ) + signatures = [ + tuple((row.species, row.magnetic) for row in site) for site in constituents + ] + else: + signatures = constituents + for permutation in permutations: + if any( + signatures[index] != signatures[int(image)] + for index, image in enumerate(permutation) + ): + raise ValueError( + "Chemistry and magnetic flags must be invariant under the SSG operations" + ) collinear_xyz = _collinear_direction_xyz(spin_only_group, q_mat, lattice) coplanar_uvw = _coplanar_perp_uvw(spin_only_group, q_mat, frame) @@ -281,6 +346,21 @@ def write_file(self, filename: str, mode: str = "wt") -> None: f.write(self._scif_str) +def _match_ssa_sites( + sas: SpinSymmetryAdaptedStructure, structure: Structure, *, symprec: float +) -> np.ndarray: + lattice = np.asarray(sas.supercell.basis) + if np.any(np.linalg.norm(structure.lattice.matrix - lattice, axis=1) > symprec): + raise ValueError("structure lattice does not match the SSA supercell") + matches = ( + Lattice(lattice).get_all_distances(structure.frac_coords, sas.supercell.positions) + <= symprec + ) + if not np.all(matches.sum(axis=0) == 1) or not np.all(matches.sum(axis=1) == 1): + raise ValueError("structure sites must match the SSA supercell uniquely") + return np.argmax(matches, axis=0) + + def _validate_moments( moments: NDArrayFloat, oriented_basis: list[np.ndarray], @@ -288,6 +368,7 @@ def _validate_moments( q_mat: NDArrayFloat, *, atol: float, + mag_symprec: float | None = None, ) -> None: """Fail fast on caller mismatches, in O(n). @@ -306,7 +387,8 @@ def _validate_moments( basis_matrix = np.stack([b.T.reshape(-1) for b in oriented_basis], axis=1) coefficients, *_ = np.linalg.lstsq(basis_matrix, target, rcond=None) residual = np.linalg.norm(basis_matrix @ coefficients - target) - if residual > atol * scale: + tolerance = atol * scale if mag_symprec is None else mag_symprec + if residual > tolerance: raise ValueError( f"magnetic_moments are not a combination of the SSA basis " f"(residual {residual:.3e}); check that they come from this " @@ -323,7 +405,7 @@ def _validate_moments( deviation = np.abs(moments @ axis) # in the plane normal to axis else: return - if np.any(deviation > atol * scale): + if np.any(deviation > tolerance): raise ValueError( f"magnetic_moments are incompatible with the " f"{spin_only_group.spin_only_group_type} spin-only axis {axis}" @@ -358,7 +440,7 @@ def _validated_occupancies( raise ValueError(f"site_occupancies must lie in (0, 1], got {occupancies.tolist()}") for permutation in permutations: permuted = occupancies[np.asarray(permutation, dtype=int)] - if not np.allclose(permuted, occupancies, atol=atol): + if not np.allclose(permuted, occupancies, atol=atol, rtol=0): raise ValueError( "site_occupancies are not invariant under the SSG operations; " "symmetry-related sites must share one occupancy." @@ -405,12 +487,12 @@ def _format_scif( a, b, c, alpha, beta, gamma = Lattice(lattice).parameters lines += [ "", - f"_cell_length_a {a:.6f}", - f"_cell_length_b {b:.6f}", - f"_cell_length_c {c:.6f}", - f"_cell_angle_alpha {alpha:.6f}", - f"_cell_angle_beta {beta:.6f}", - f"_cell_angle_gamma {gamma:.6f}", + f"_cell_length_a {_number(a)}", + f"_cell_length_b {_number(b)}", + f"_cell_length_c {_number(c)}", + f"_cell_angle_alpha {_number(alpha)}", + f"_cell_angle_beta {_number(beta)}", + f"_cell_angle_gamma {_number(gamma)}", "", "loop_", "_space_group_symop_spin_operation.id", @@ -445,8 +527,8 @@ def _format_scif( if constituent.magnetic: labels.append(label) lines.append( - f"{label} {constituent.species} {_fmt(fx)} {_fmt(fy)} {_fmt(fz)} " - f"{constituent.occupancy:.12g} {multiplicity}" + f"{label} {_species_token(constituent.species)} {_number(fx)} {_number(fy)} {_number(fz)} " + f"{_number(constituent.occupancy)} {multiplicity}" ) magnetic_labels.append(labels) @@ -455,11 +537,11 @@ def _format_scif( orbits, magnetic_labels, site_moment_subspaces, strict=True ): m_uvw = frame.T @ moments[rep] - components = [_fmt(x) for x in m_uvw] + components = [_number(x) for x in m_uvw] symmform = _symmform(subspace, atol=atol) for label in labels: moment_lines.append( - f"{label} {' '.join(components)} {symmform} {_fmt(float(np.linalg.norm(m_uvw)))}" + f"{label} {' '.join(components)} {symmform} {_number(float(np.linalg.norm(m_uvw)))}" ) if moment_lines: lines += [ diff --git a/tests/scif_tests/test_compositions.py b/tests/scif_tests/test_compositions.py index c133e7c..bde5890 100644 --- a/tests/scif_tests/test_compositions.py +++ b/tests/scif_tests/test_compositions.py @@ -4,10 +4,11 @@ import numpy as np import pytest -from pymatgen.core import Composition +from pymatgen.core import Composition, Structure, get_el_sp from pymatgen.io.cif import CifFile from spinspg.spin import SpinOnlyGroupType +from spinforge.mcif import MCifWriter from spinforge.scif import SpinCifReader, SpinCifWriter @@ -20,13 +21,25 @@ def antitranslation_model(enumerated): ) -def _chemistry(ssa, magnetic_composition, selected): - compositions = [ - Composition(magnetic_composition if number == 25 else {"Te2-": 0.9}) - for number in ssa.supercell.numbers - ] - selections = [selected if number == 25 else () for number in ssa.supercell.numbers] - return compositions, selections +def _structure(ssa, magnetic_composition, selected, moments): + selected = {get_el_sp(species) for species in selected} + compositions, positions, row_moments, magnetic = [], [], [], [] + for number, position, moment in zip( + ssa.supercell.numbers, ssa.supercell.positions, moments, strict=True + ): + composition = Composition(magnetic_composition if number == 25 else {"Te2-": 0.9}) + for species, occupancy in composition.items(): + is_magnetic = species in selected + compositions.append({species: occupancy}) + positions.append(position) + row_moments.append(moment if is_magnetic else np.zeros(3)) + magnetic.append(is_magnetic) + return Structure( + ssa.supercell.basis, + compositions, + positions, + site_properties={"magmom": row_moments, "magnetic": magnetic}, + ) @pytest.mark.parametrize( @@ -37,6 +50,10 @@ def _chemistry(ssa, magnetic_composition, selected): ({"Mn2+": 0.6, "Mn3+": 0.2}, ("Mn2+",)), ({"Mn": 0.6, "Mn0+": 0.2}, ("Mn0+",)), ({"Mn2+": 0.6, "Zn2+": 1e-8}, ("Mn2+",)), + ( + {"Mn0.3333333333333333+": 0.6, "Mn0.3333333333333334+": 0.2}, + ("Mn0.3333333333333333+",), + ), ], ) @pytest.mark.parametrize("improper", [False, True]) @@ -49,54 +66,45 @@ def test_mixed_antitranslation_with_explicit_frame( rotation = rotation @ np.diag([1.0, 1.0, -1.0]) (intrinsic,) = ssa.magnetic_moments_basis moments = intrinsic @ rotation.T - compositions, selections = _chemistry(ssa, magnetic_composition, selected) + original = _structure(ssa, magnetic_composition, selected, moments) text = str( - SpinCifWriter( + SpinCifWriter.from_structure( ssa, group, - moments, + Structure.from_sites(list(reversed(original))), spin_only_group=sog, spin_rotation=rotation, - site_compositions=compositions, - magnetic_species=selections, + symprec=1e-5, ) ) reader = SpinCifReader(text) structure = reader.get_structure() - assert len(structure) == ssa.supercell.num_supercell_sites + assert len(structure) == len(original) assert any(reversal == -1 for _, _, reversal, _ in reader.spin_lattice) embedding = np.linalg.inv(structure.lattice.matrix) @ ssa.supercell.basis for site in structure: - indices = [ - index - for index, position in enumerate(ssa.supercell.positions) + matches = [ + expected + for expected in original if np.allclose( - site.frac_coords - position - np.rint(site.frac_coords - position), + site.frac_coords + - expected.frac_coords + - np.rint(site.frac_coords - expected.frac_coords), 0, atol=1e-5, + rtol=0, ) + and set(site.species) == set(expected.species) ] - (index,) = indices - assert site.species.as_dict() == compositions[index].as_dict() - assert set(site.properties["magnetic_species"]) == set(selections[index]) + (expected,) = matches + assert dict(site.species.items()) == dict(expected.species.items()) + assert site.properties["magnetic"] == expected.properties["magnetic"] np.testing.assert_allclose( - site.properties["magmom"] @ embedding, moments[index], atol=1e-5 + site.properties["magmom"] @ embedding, + expected.properties["magmom"], + atol=1e-5, + rtol=0, ) - assert set(site.properties["species_magmoms"]) == { - str(species) for species in compositions[index] - } - mean = np.zeros(3) - for species, occupancy in compositions[index].items(): - actual = site.properties["species_magmoms"][str(species)] @ embedding - expected = moments[index] if str(species) in selections[index] else np.zeros(3) - np.testing.assert_allclose(actual, expected, atol=1e-5) - mean += occupancy * actual - expected_mean = moments[index] * sum( - occupancy - for species, occupancy in compositions[index].items() - if str(species) in selections[index] - ) - np.testing.assert_allclose(mean, expected_mean, atol=1e-5) direction = np.array([float(value) for value in reader.collinear_direction_xyz.split(",")]) actual_axis = direction @ ssa.supercell.basis @@ -125,15 +133,19 @@ def test_mixed_antitranslation_with_explicit_frame( def test_selected_zero_moments_remain_selected(antitranslation_model): sog, group, ssa = antitranslation_model - compositions, selections = _chemistry(ssa, {"Mn2+": 0.6, "Zn2+": 0.2}, ("Mn2+",)) + original = _structure( + ssa, + {"Mn2+": 0.6, "Zn2+": 0.2}, + ("Mn2+",), + np.zeros((ssa.supercell.num_supercell_sites, 3)), + ) text = str( - SpinCifWriter( + SpinCifWriter.from_structure( ssa, group, - np.zeros((ssa.supercell.num_supercell_sites, 3)), + original, spin_only_group=sog, - site_compositions=compositions, - magnetic_species=selections, + symprec=1e-5, ) ) (block,) = CifFile.from_str(text).data.values() @@ -147,8 +159,7 @@ def test_selected_zero_moments_remain_selected(antitranslation_model): } structure = SpinCifReader(text).get_structure() assert all( - site.properties["magnetic_species"] == (("Mn2+",) if "Mn2+" in site.species else ()) - for site in structure + site.properties["magnetic"] == (get_el_sp("Mn2+") in site.species) for site in structure ) np.testing.assert_array_equal(structure.site_properties["magmom"], 0) @@ -157,20 +168,31 @@ def test_selected_zero_moments_remain_selected(antitranslation_model): def test_orbit_inconsistent_constituents_are_rejected(antitranslation_model, change_selection): sog, group, ssa = antitranslation_model (moments,) = ssa.magnetic_moments_basis - compositions, selections = _chemistry(ssa, {"Mn2+": 0.6, "Fe2+": 0.2}, ("Mn2+",)) - index = next(index for index, number in enumerate(ssa.supercell.numbers) if number == 25) + structure = _structure(ssa, {"Mn2+": 0.6, "Fe2+": 0.2}, ("Mn2+",), moments) + index = next( + index for index, site in enumerate(structure) if get_el_sp("Mn2+") in site.species + ) + site = structure[index] if change_selection: - selections[index] = ("Fe2+",) + other = next( + other + for other in structure + if get_el_sp("Fe2+") in other.species + and np.array_equal(other.frac_coords, site.frac_coords) + ) + other.properties["magnetic"] = True + other.properties["magmom"] = site.properties["magmom"] + site.properties["magnetic"] = False + site.properties["magmom"] = np.zeros(3) else: - compositions[index] = Composition({"Mn2+": 0.5, "Fe2+": 0.2}) + structure.replace(index, {"Mn2+": 0.5}, properties=site.properties) with pytest.raises(ValueError, match="invariant under the SSG operations"): - SpinCifWriter( + SpinCifWriter.from_structure( ssa, group, - moments, + structure, spin_only_group=sog, - site_compositions=compositions, - magnetic_species=selections, + symprec=1e-5, ) @@ -182,88 +204,7 @@ def test_invalid_explicit_frames_are_rejected(antitranslation_model, rotation): SpinCifWriter(ssa, group, moments, spin_only_group=sog, spin_rotation=rotation) -def test_compositions_and_legacy_occupancies_are_exclusive(antitranslation_model): - sog, group, ssa = antitranslation_model - (moments,) = ssa.magnetic_moments_basis - compositions, _ = _chemistry(ssa, {"Mn": 0.5}, ("Mn",)) - with pytest.raises(ValueError, match="either site_compositions or site_occupancies"): - SpinCifWriter( - ssa, - group, - moments, - spin_only_group=sog, - site_compositions=compositions, - site_occupancies=[0.5] * len(compositions), - ) - - -@pytest.mark.parametrize("legacy_occupancies", [False, True]) -def test_writer_occupancy_bound_uses_caller_tolerance(antitranslation_model, legacy_occupancies): - sog, group, ssa = antitranslation_model - (moments,) = ssa.magnetic_moments_basis - occupancy = 1 + 5e-9 - if legacy_occupancies: - chemistry = {"site_occupancies": [occupancy] * ssa.supercell.num_supercell_sites} - else: - compositions, selections = _chemistry(ssa, {"Mn": occupancy}, ("Mn",)) - chemistry = {"site_compositions": compositions, "magnetic_species": selections} - - with pytest.raises(ValueError, match="occupanc"): - SpinCifWriter(ssa, group, moments, spin_only_group=sog, atol=1e-12, **chemistry) - - text = str(SpinCifWriter(ssa, group, moments, spin_only_group=sog, atol=1e-6, **chemistry)) - (block,) = CifFile.from_str(text).data.values() - magnetic_occupancies = [ - float(amount) - for species, amount in zip( - block.data["_atom_site_type_symbol"], block.data["_atom_site_occupancy"], strict=True - ) - if species == "Mn" - ] - assert magnetic_occupancies - assert all(amount == occupancy for amount in magnetic_occupancies) - if legacy_occupancies: - chemistry = {"site_occupancies": [1 + 5e-7] * ssa.supercell.num_supercell_sites} - else: - compositions, selections = _chemistry(ssa, {"Mn": 1 + 5e-7}, ("Mn",)) - chemistry = {"site_compositions": compositions, "magnetic_species": selections} - with pytest.raises(ValueError, match="occupanc"): - SpinCifWriter(ssa, group, moments, spin_only_group=sog, atol=1e-6, **chemistry) - - -@pytest.mark.parametrize("legacy_occupancies", [False, True]) -def test_writer_occupancy_invariance_uses_caller_tolerance( - antitranslation_model, legacy_occupancies -): - sog, group, ssa = antitranslation_model - (moments,) = ssa.magnetic_moments_basis - index = next(index for index, number in enumerate(ssa.supercell.numbers) if number == 25) - if legacy_occupancies: - occupancies = [0.8] * ssa.supercell.num_supercell_sites - occupancies[index] += 5e-7 - chemistry = {"site_occupancies": occupancies} - else: - compositions, selections = _chemistry(ssa, {"Mn": 0.8}, ("Mn",)) - compositions[index] = Composition({"Mn": 0.8 + 5e-7}) - chemistry = {"site_compositions": compositions, "magnetic_species": selections} - - with pytest.raises(ValueError, match="invariant under the SSG operations"): - SpinCifWriter(ssa, group, moments, spin_only_group=sog, atol=1e-8, **chemistry) - - text = str(SpinCifWriter(ssa, group, moments, spin_only_group=sog, atol=1e-6, **chemistry)) - (block,) = CifFile.from_str(text).data.values() - magnetic_occupancies = [ - float(amount) - for species, amount in zip( - block.data["_atom_site_type_symbol"], block.data["_atom_site_occupancy"], strict=True - ) - if species == "Mn" - ] - assert magnetic_occupancies - np.testing.assert_allclose(magnetic_occupancies, 0.8, atol=1e-6, rtol=0) - - -def test_reader_rejects_different_selected_constituent_moments(external_scif): +def test_reader_preserves_different_coincident_constituent_moments(external_scif): text = external_scif.replace( "_atom_site_symmetry_multiplicity\nMn1 Mn 0.25 0.0 0.0 2", "_atom_site_symmetry_multiplicity\n_atom_site_occupancy\n" @@ -272,5 +213,41 @@ def test_reader_rejects_different_selected_constituent_moments(external_scif): "Mn1 0.0 0.0 3.0 0,0,w 0,0,w 3.0", "Mn1 0.0 0.0 3.0 0,0,w 0,0,w 3.0\nFe1 0.0 0.0 2.0 0,0,w 0,0,w 2.0", ) - with pytest.raises(ValueError, match="Unequal selected constituent moments"): - SpinCifReader(text).get_structure() + structure = SpinCifReader(text).get_structure() + assert len(structure) == 4 + assert all(site.properties["magnetic"] for site in structure) + for site in structure: + (species,) = site.species + np.testing.assert_allclose( + site.properties["magmom"], [0, 0, 3 if species == get_el_sp("Mn") else 2] + ) + assert dict(site.species.items()) == {species: 0.5} + + +def test_reader_constituent_rows_can_be_exported_to_mcif(external_scif): + text = external_scif.replace( + "_atom_site_symmetry_multiplicity\nMn1 Mn 0.25 0.0 0.0 2", + "_atom_site_symmetry_multiplicity\n_atom_site_occupancy\n" + "Mn1 Mn0.3333333333333333+ 0.25 0.0 0.0 2 0.5\n" + "Mn2 Mn0.3333333333333334+ 0.25 0.0 0.0 2 0.25", + ).replace("Mn1 0.0 0.0 3.0 0,0,w 0,0,w 3.0", "Mn1 0.0 0.0 0.0 0,0,w 0,0,w 0.0") + output = str(MCifWriter(SpinCifReader(text).get_structure(), symprec=1e-5)) + (block,) = CifFile.from_str(output).data.values() + assert set( + zip(block.data["_atom_site_type_symbol"], block.data["_atom_site_occupancy"], strict=True) + ) == { + ("Mn0.3333333333333333+", "0.5"), + ("Mn0.3333333333333334+", "0.25"), + } + magnetic_labels = set(block.data["_atom_site_moment.label"]) + assert magnetic_labels == { + label + for label, symbol in zip( + block.data["_atom_site_label"], block.data["_atom_site_type_symbol"], strict=True + ) + if symbol == "Mn0.3333333333333333+" + } + for axis in "xyz": + assert all( + float(value) == 0 for value in block.data[f"_atom_site_moment.crystalaxis_{axis}"] + ) diff --git a/tests/scif_tests/test_reader.py b/tests/scif_tests/test_reader.py index ae9a9e8..2154781 100644 --- a/tests/scif_tests/test_reader.py +++ b/tests/scif_tests/test_reader.py @@ -6,7 +6,7 @@ import numpy as np import pytest -from pymatgen.core import Composition +from pymatgen.core import Element, Species from spinspg.spin import SpinOnlyGroupType from spinforge.scif import SpinCifReader, SpinCifWriter @@ -47,6 +47,54 @@ def test_reader_rejects_malformed_moment_columns(external_scif, columns): SpinCifReader(text) +@pytest.mark.parametrize( + "short_column", + ["label", "type_symbol", "fract_x", "fract_y", "fract_z", "symmetry_multiplicity"], +) +def test_reader_rejects_mismatched_atom_site_columns(external_scif, short_column): + columns = { + "label": ["Mn1", "Mn2"], + "type_symbol": ["Mn", "Mn"], + "fract_x": ["0.25", "0.75"], + "fract_y": ["0", "0"], + "fract_z": ["0", "0"], + "symmetry_multiplicity": ["2", "2"], + } + columns[short_column].pop() + before, rest = external_scif.split("loop_\n_atom_site_label") + _, after = rest.split("loop_\n_atom_site_spin_moment.label") + text = ( + before + + "".join( + f"loop_\n_atom_site_{column}\n{' '.join(values)}\n" + for column, values in columns.items() + ) + + "loop_\n_atom_site_spin_moment.label" + + after + ) + with pytest.raises(ValueError, match="equal-length loop columns"): + SpinCifReader(text) + + +@pytest.mark.parametrize("category", ["operation", "lattice"]) +@pytest.mark.parametrize("missing_position", [False, True]) +def test_reader_rejects_mismatched_symmetry_columns(external_scif, category, missing_position): + columns = ( + f"loop_\n_space_group_symop_spin_{category}.xyzt\nx,y,z,+1 -x,-y,z,+1\n" + f"loop_\n_space_group_symop_spin_{category}.uvw\nu,v,w\n" + ) + if missing_position: + columns = f"loop_\n_space_group_symop_spin_{category}.uvw\nu,v,w\n" + if category == "operation": + before, rest = external_scif.split("loop_\n_space_group_symop_spin_operation.id") + _, after = rest.split("loop_\n_atom_site_label") + text = before + columns + "loop_\n_atom_site_label" + after + else: + text = external_scif + columns + with pytest.raises(ValueError, match="equal-length loop columns"): + SpinCifReader(text) + + @pytest.mark.parametrize( ("type_rows", "message"), [ @@ -129,7 +177,7 @@ def test_reader_spin_frame_uses_symprec(external_scif, tmp_path, from_file): assert np.allclose(structure.site_properties["magmom"], [[0, 0, 3], [0, 0, 3]]) -@pytest.mark.parametrize("scale", [1e-200, 1e-10, 1e200]) +@pytest.mark.parametrize("scale", [1e-200, 1e-10, 1e200, 1e308]) def test_reader_normalizes_nonzero_spin_axes(external_scif, scale): text = external_scif.replace("'a,b,c'", f"'{scale:.200f}*a,b,c'") reader = SpinCifReader(text, symprec=1e-6) @@ -137,33 +185,28 @@ def test_reader_normalizes_nonzero_spin_axes(external_scif, scale): assert np.allclose(reader.get_structure().site_properties["magmom"], [[0, 0, 3], [0, 0, 3]]) -@pytest.mark.parametrize("mixed", [False, True]) -def test_reader_occupancy_uses_symprec_without_normalizing(external_scif, mixed): +@pytest.mark.parametrize("symprec", [1e-10, 1e-2]) +def test_reader_occupancy_does_not_depend_on_geometric_tolerance(external_scif, symprec): text = external_scif.replace( "_atom_site_symmetry_multiplicity\n", "_atom_site_symmetry_multiplicity\n_atom_site_occupancy\n", ).replace( "Mn1 Mn 0.25 0.0 0.0 2", - "Mn1 Mn 0.25 0.0 0.0 2 0.5\nFe1 Fe 0.25 0.0 0.0 2 0.500000005" - if mixed - else "Mn1 Mn 0.25 0.0 0.0 2 1.000000005", + "Mn1 Mn 0.25 0.0 0.0 2 1.000000005", ) - with pytest.raises(ValueError, match="occupanc"): - SpinCifReader(text, symprec=1e-10).get_structure() - structure = SpinCifReader(text, symprec=1e-6).get_structure() - expected = {"Mn": 0.5, "Fe": 0.500000005} if mixed else {"Mn": 1.000000005} + structure = SpinCifReader(text, symprec=symprec).get_structure() assert len(structure) == 2 - assert all(site.species.as_dict() == expected for site in structure) + assert all(dict(site.species.items()) == {Element("Mn"): 1.000000005} for site in structure) -def test_reader_occupancy_respects_structure_limit(external_scif): - occupancy = 0.5 + 2 * Composition.amount_tolerance +@pytest.mark.parametrize("occupancy", [1.01, 1e-10]) +def test_reader_rejects_occupancy_that_pymatgen_cannot_preserve(external_scif, occupancy): text = external_scif.replace( "_atom_site_symmetry_multiplicity\n", "_atom_site_symmetry_multiplicity\n_atom_site_occupancy\n", ).replace( "Mn1 Mn 0.25 0.0 0.0 2", - f"Mn1 Mn 0.25 0.0 0.0 2 0.5\nFe1 Fe 0.25 0.0 0.0 2 {occupancy}", + f"Mn1 Mn 0.25 0.0 0.0 2 {occupancy}", ) with pytest.raises(ValueError, match="occupanc"): SpinCifReader(text, symprec=1e-2).get_structure() @@ -177,18 +220,59 @@ def test_reader_orbit_positions_use_absolute_symprec(external_scif): SpinCifReader(text, symprec=1e-3).get_structure() -def test_reader_constituent_positions_use_absolute_symprec(external_scif): +@pytest.mark.parametrize("symprec", [1e-7, 1e-3]) +def test_reader_preserves_distinct_constituent_rows(external_scif, symprec): text = external_scif.replace("2 -x,-y,z,+1 u,v,w\n", "").replace( "Mn1 Mn 0.25 0.0 0.0 2", "Mn1 Mn 0.000002 0.0 0.0 1\nFe1 Fe 0.999998 0.0 0.0 1", ) - structure = SpinCifReader(text, symprec=1e-7).get_structure() + structure = SpinCifReader(text, symprec=symprec).get_structure() assert {site.species_string: site.frac_coords[0] for site in structure} == { "Mn": 0.000002, "Fe": 0.999998, } - with pytest.raises(ValueError, match="occupanc"): - SpinCifReader(text, symprec=1e-3).get_structure() + assert structure.site_properties["magnetic"] == [True, False] + + +@pytest.mark.parametrize("from_file", [False, True]) +def test_reader_frame_and_moment_tolerances_are_independent(external_scif, tmp_path, from_file): + text = ( + external_scif.replace("'a,b,c'", "'a,b+0.0005*a,c'") + .replace("Mn1 Mn 0.25 0.0 0.0 2", "Mn1 Mn 0.0 0.0 0.0 1") + .replace("2 -x,-y,z,+1 u,v,w", "2 -x,-y,z,+1 u,v,1.0001*w") + ) + source = text + read = SpinCifReader + if from_file: + path = tmp_path / "tolerances.scif" + path.write_text(text) + source = str(path) + read = SpinCifReader.from_file + with pytest.raises(ValueError, match="not orthogonal"): + read(source, symprec=1e-2, mag_symprec=1e-2, atol=1e-6) + with pytest.raises(ValueError, match="Inconsistent moments"): + read(source, symprec=1e-2, mag_symprec=1e-6, atol=1e-2).get_structure() + with pytest.raises(ValueError, match="Inconsistent moments"): + read(source, symprec=1e-6, atol=1e-2).get_structure() + assert len(read(source, symprec=1e-2).get_structure()) == 1 + structure = read(source, symprec=1e-6, mag_symprec=1e-2, atol=1e-2).get_structure() + assert len(structure) == 1 + np.testing.assert_allclose(structure.site_properties["magmom"], [[0, 0, 3]]) + + +def test_reader_preserves_fractional_oxidation_identities(external_scif): + text = external_scif.replace( + "_atom_site_symmetry_multiplicity\nMn1 Mn 0.25 0.0 0.0 2", + "_atom_site_symmetry_multiplicity\n_atom_site_occupancy\n" + "Mn1 Mn0.3333333333333333+ 0.25 0.0 0.0 2 0.5\n" + "Mn2 Mn0.3333333333333334+ 0.25 0.0 0.0 2 0.5", + ) + structure = SpinCifReader(text).get_structure() + assert len(structure) == 4 + assert {tuple(site.species.items()) for site in structure} == { + ((Species("Mn", 0.3333333333333333), 0.5),), + ((Species("Mn", 0.3333333333333334), 0.5),), + } def test_reader_legacy_cartn_frame_fallback(enumerated): diff --git a/tests/scif_tests/test_structure.py b/tests/scif_tests/test_structure.py new file mode 100644 index 0000000..333eb15 --- /dev/null +++ b/tests/scif_tests/test_structure.py @@ -0,0 +1,183 @@ +"""Structure-based adaptation to the SSA writer.""" + +from __future__ import annotations + +import numpy as np +import pytest +from pymatgen.core import Species, Structure +from pymatgen.io.cif import CifFile +from spinspg.spin import SpinOnlyGroup, SpinOnlyGroupType + +from spinforge.scif import SpinCifReader, SpinCifWriter + + +@pytest.fixture +def collinear_model(enumerated): + return next( + entry for entry in enumerated("MnTe", SpinOnlyGroupType.COLLINEAR, 1) if entry[2].dim == 1 + ) + + +def test_reordered_structure_matches_ssa_sites(collinear_model): + sog, group, ssa = collinear_model + (moments,) = ssa.magnetic_moments_basis + structure = ssa.generate_with_magnetic_moments(moments) + reversed_structure = Structure.from_sites(list(reversed(structure))) + reference = SpinCifWriter.from_structure( + ssa, group, structure, spin_only_group=sog, symprec=1e-8 + ) + reordered = SpinCifWriter.from_structure( + ssa, group, reversed_structure, spin_only_group=sog, symprec=1e-8 + ) + assert str(reordered) == str(reference) + reader = SpinCifReader(str(reordered)) + for rotation, translation, _, _ in reader.expanded_operations(): + transformed = structure.frac_coords @ rotation.T + translation + distances = structure.lattice.get_all_distances(transformed, structure.frac_coords) + assert np.all(np.min(distances, axis=1) <= 1e-8) + + +def test_constituent_partition_does_not_break_ssg_orbits(collinear_model): + sog, group, ssa = collinear_model + (moments,) = ssa.magnetic_moments_basis + structure = ssa.generate_with_magnetic_moments(moments) + species, positions, vectors = [], [], [] + split_done = False + for site in structure: + amounts = [0.5, 0.5] if site.specie.symbol == "Mn" and not split_done else [1.0] + split_done = split_done or site.specie.symbol == "Mn" + for amount in amounts: + species.append({site.specie: amount}) + positions.append(site.frac_coords) + vectors.append(site.properties["magmom"]) + split = Structure(structure.lattice, species, positions, site_properties={"magmom": vectors}) + original = SpinCifWriter.from_structure( + ssa, group, structure, spin_only_group=sog, symprec=1e-8 + ) + partitioned = SpinCifWriter.from_structure( + ssa, group, split, spin_only_group=sog, symprec=1e-8 + ) + assert str(partitioned) == str(original) + + +@pytest.mark.parametrize("mismatch", ["position", "lattice", "missing"]) +def test_structure_must_match_ssa_supercell(collinear_model, mismatch): + sog, group, ssa = collinear_model + (moments,) = ssa.magnetic_moments_basis + structure = ssa.generate_with_magnetic_moments(moments) + index = next(index for index, site in enumerate(structure) if site.specie.symbol == "Mn") + if mismatch == "position": + structure.translate_sites(index, [0.01, 0, 0]) + elif mismatch == "lattice": + structure.scale_lattice(structure.volume * 1.01) + else: + structure.remove_sites([index]) + with pytest.raises(ValueError, match="match.*SSA supercell"): + SpinCifWriter.from_structure(ssa, group, structure, spin_only_group=sog, symprec=1e-8) + + +def test_explicit_moment_tolerance_controls_ssa_residual(collinear_model): + sog, group, ssa = collinear_model + (basis,) = ssa.magnetic_moments_basis + moments = basis * 1000 + index = next(index for index, number in enumerate(ssa.supercell.numbers) if number == 25) + moments[index, 0] += 1e-4 + structure = ssa.generate_with_magnetic_moments(moments) + with pytest.raises(ValueError, match="SSA basis"): + SpinCifWriter.from_structure( + ssa, group, structure, spin_only_group=sog, symprec=1e-8, mag_symprec=1e-6 + ) + SpinCifWriter.from_structure( + ssa, group, structure, spin_only_group=sog, symprec=1e-8, mag_symprec=1e-3 + ) + + +def test_explicit_moment_tolerance_controls_spin_only_axis(collinear_model): + _, group, ssa = collinear_model + (moments,) = ssa.magnetic_moments_basis + structure = ssa.generate_with_magnetic_moments(moments * 1000) + axis = np.array([1e-4, 0, 1]) + sog = SpinOnlyGroup(SpinOnlyGroupType.COLLINEAR, axis / np.linalg.norm(axis)) + with pytest.raises(ValueError, match="spin-only axis"): + SpinCifWriter.from_structure( + ssa, group, structure, spin_only_group=sog, symprec=1e-8, mag_symprec=1e-6 + ) + SpinCifWriter.from_structure( + ssa, group, structure, spin_only_group=sog, symprec=1e-8, mag_symprec=1 + ) + + +def test_from_oriented_preserves_structure_chemistry(collinear_model): + sog, group, ssa = collinear_model + for structure, msg in ssa.generate_oriented(sog, group, rng=np.random.default_rng(4)): + species = [ + {Species("Mn", 0.331): 0.6, Species("Zn", 2): 0.2} + if site.specie.symbol == "Mn" + else {Species("Te", -2): 0.9} + for site in structure + ] + mixed = Structure( + structure.lattice, + species, + structure.frac_coords, + site_properties=structure.site_properties, + ) + text = str(SpinCifWriter.from_oriented(ssa, group, mixed, msg, spin_only_group=sog)) + read = SpinCifReader(text).get_structure() + actual = {(next(iter(site.species)), next(iter(site.species.values()))) for site in read} + assert actual == { + (Species("Mn", 0.331), 0.6), + (Species("Zn", 2), 0.2), + (Species("Te", -2), 0.9), + } + with pytest.raises(ValueError, match="only override ordered input"): + SpinCifWriter.from_oriented( + ssa, group, mixed, msg, spin_only_group=sog, site_occupancies=[1] * len(mixed) + ) + + +def test_from_oriented_legacy_occupancy_override_respects_atol(collinear_model): + sog, group, ssa = collinear_model + for structure, msg in ssa.generate_oriented(sog, group, rng=np.random.default_rng(4)): + occupancies = [0.8] * len(structure) + index = next(index for index, site in enumerate(structure) if site.specie.symbol == "Mn") + occupancies[index] += 5e-7 + SpinCifWriter.from_oriented( + ssa, + group, + structure, + msg, + spin_only_group=sog, + site_occupancies=occupancies, + atol=1e-6, + ) + with pytest.raises(ValueError, match="not invariant"): + SpinCifWriter.from_oriented( + ssa, + group, + structure, + msg, + spin_only_group=sog, + site_occupancies=occupancies, + atol=1e-8, + ) + + +def test_structure_occupancy_is_independent_of_geometric_tolerance(collinear_model): + sog, group, ssa = collinear_model + (moments,) = ssa.magnetic_moments_basis + source = ssa.generate_with_magnetic_moments(moments) + occupancy = 1 + 5e-9 + structure = Structure( + source.lattice, + [{site.specie: occupancy} for site in source], + source.frac_coords, + site_properties=source.site_properties, + ) + text = str( + SpinCifWriter.from_structure( + ssa, group, structure, spin_only_group=sog, symprec=1e-10, atol=1e-12 + ) + ) + (block,) = CifFile.from_str(text).data.values() + assert all(float(amount) == occupancy for amount in block.data["_atom_site_occupancy"]) diff --git a/tests/scif_tests/test_writer.py b/tests/scif_tests/test_writer.py index 1f033ef..99dea0d 100644 --- a/tests/scif_tests/test_writer.py +++ b/tests/scif_tests/test_writer.py @@ -5,6 +5,7 @@ import numpy as np import pytest from moyopy import SpaceGroupType +from pymatgen.io.cif import CifFile from spinspg.spin import SpinOnlyGroupType from spinforge.scif import SpinCifReader, SpinCifWriter @@ -97,19 +98,50 @@ def test_snap_translation_respects_supercell_grid(): assert np.array_equal(_snap_translation(noisy, denom=denom, atol=1e-5), [0.0, 0.0, 0.2]) -def test_weak_moments_are_not_dropped(enumerated): +@pytest.mark.parametrize("magnitude", [4e-6, 1e-10]) +def test_weak_moments_are_not_dropped(enumerated, magnitude): # Moments below the old norm-atol cutoff but above the print resolution # must still be written; otherwise weak magnetic structures silently # round-trip as nonmagnetic (the reader treats omitted rows as zero). sog, nssg, sas = enumerated("MnTe", SpinOnlyGroupType.COLLINEAR, 1)[0] moments = sas._sample_magnetic_moments(rng=np.random.default_rng(0)) - weak = moments * (4e-6 / np.max(np.abs(moments))) + weak = moments * (magnitude / np.max(np.abs(moments))) text = str(SpinCifWriter(sas, nssg, weak, spin_only_group=sog)) assert "_atom_site_spin_moment.label" in text structure = SpinCifReader(text).get_structure() assert np.max(np.abs(structure.site_properties["magmom"])) > 0 +def test_legacy_occupancy_tolerance_is_preserved(enumerated): + sog, group, ssa = next( + entry + for entry in enumerated("MnTe", SpinOnlyGroupType.COLLINEAR, 2) + if entry[1].spin_translation_coset.size > 1 and entry[2].dim == 1 + ) + (moments,) = ssa.magnetic_moments_basis + amounts = [0.8] * ssa.supercell.num_supercell_sites + index = next(index for index, number in enumerate(ssa.supercell.numbers) if number == 25) + amounts[index] += 5e-7 + with pytest.raises(ValueError, match="not invariant"): + SpinCifWriter( + ssa, group, moments, spin_only_group=sog, site_occupancies=amounts, atol=1e-8 + ) + SpinCifWriter(ssa, group, moments, spin_only_group=sog, site_occupancies=amounts, atol=1e-6) + + amounts = [1 + 5e-7] * ssa.supercell.num_supercell_sites + with pytest.raises(ValueError, match="occupancies"): + SpinCifWriter( + ssa, group, moments, spin_only_group=sog, site_occupancies=amounts, atol=1e-8 + ) + text = str( + SpinCifWriter( + ssa, group, moments, spin_only_group=sog, site_occupancies=amounts, atol=1e-6 + ) + ) + (block,) = CifFile.from_str(text).data.values() + assert all(float(amount) == 1 + 5e-7 for amount in block.data["_atom_site_occupancy"]) + + def test_writer_declares_p_abc_frame(enumerated): # The frame must go out via transform_spinframe_P_abc (the tag the # COMCIFS reference files and VESTA understand), also on a diff --git a/tests/test_mcif.py b/tests/test_mcif.py index 747adcd..1bae232 100644 --- a/tests/test_mcif.py +++ b/tests/test_mcif.py @@ -139,6 +139,23 @@ def test_translation_snapping_respects_cartesian_symprec(symprec) -> None: assert np.linalg.norm(displacement @ structure.lattice.matrix) <= symprec +def test_general_origin_translation_retains_requested_precision() -> None: + structure = Structure( + Lattice.from_parameters(4, 5, 6, 70, 80, 75), + ["Fe"], + [[0.123456789, 0.234567891, 0.345678912]], + site_properties={"magmom": [[0, 0, 2]]}, + ) + text = str(MCifWriter(structure, symprec=1e-12)) + (block,) = CifFile.from_str(text).data.values() + (site,) = structure + for operation in CifParser.from_str(text).get_magsymops(block): + distance, _ = structure.lattice.get_distance_and_image( + site.frac_coords, operation.operate(site.frac_coords) + ) + assert distance <= 1e-12 + + def test_save_requires_magmom_property() -> None: from pymatgen.core import Lattice, Structure @@ -188,7 +205,7 @@ def test_mixed_site_retains_minority_species_and_occupied_atom_moment() -> None: [[0, 0, 0], [0, 0, 0.5]], site_properties={"magmom": [[1, 2, 3], [-1, -2, -3]]}, ) - text = str(MCifWriter(mixed, magnetic_species=[("Fe2+",), ("Fe2+",)])) + text = str(MCifWriter(_constituent_structure(mixed, [("Fe2+",), ("Fe2+",)]))) images = _expand_constituent_rows(text) assert len(images) == 4 for symbol, occupancy, position, moment, selected in images: @@ -213,7 +230,7 @@ def test_symmetry_preserves_constituent_selection_even_for_zero_moments(moment) [[0, 0, 0], [0, 0, 0.5]], site_properties={"magmom": [moment, moment]}, ) - text = str(MCifWriter(structure, magnetic_species=[("Fe",), ("Mn",)])) + text = str(MCifWriter(_constituent_structure(structure, [("Fe",), ("Mn",)]))) block = next(iter(CifFile.from_str(text).data.values())).data assert len(block["_atom_site_label"]) == 4 assert len(block["_atom_site_moment.label"]) == 2 @@ -232,7 +249,7 @@ def test_oxidation_states_and_neutral_elements_remain_distinct() -> None: [[0, 0, 0]], site_properties={"magmom": [[0, 0, 2]]}, ) - text = str(MCifWriter(structure, magnetic_species=[("Fe2+",)])) + text = str(MCifWriter(_constituent_structure(structure, [("Fe2+",)]))) data = next(iter(CifFile.from_str(text).data.values())).data assert set(data["_atom_site_type_symbol"]) == {"Fe", "Fe0+", "Fe2+"} assert dict(zip(data["_atom_type_symbol"], data["_atom_type_oxidation_number"])) == { @@ -277,28 +294,120 @@ def test_small_minority_occupancy_is_not_rounded_to_zero() -> None: [[0, 0, 0]], site_properties={"magmom": [[0, 0, 2]]}, ) - images = _expand_constituent_rows(str(MCifWriter(structure, magnetic_species=[("Fe",)]))) + images = _expand_constituent_rows( + str(MCifWriter(_constituent_structure(structure, [("Fe",)]))) + ) (iron,) = (row for row in images if row[0] == "Fe") assert iron[1] == 1e-7 assert np.allclose(iron[3], [0, 0, 2]) +def test_fractional_oxidation_states_and_weak_moments_are_lossless() -> None: + structure = Structure( + Lattice.cubic(4), + [{Species("Fe", 0.331): 0.4, Species("Fe", 0.334): 0.5}], + [[0, 0, 0]], + site_properties={"magmom": [[0, 0, 1e-10]]}, + ) + (block,) = CifFile.from_str(str(MCifWriter(structure))).data.values() + assert block.data["_atom_site_type_symbol"] == ["Fe0.331+", "Fe0.334+"] + assert dict( + zip( + block.data["_atom_type_symbol"], + map(float, block.data["_atom_type_oxidation_number"]), + strict=True, + ) + ) == {"Fe0.331+": 0.331, "Fe0.334+": 0.334} + assert list(map(float, block.data["_atom_site_moment.crystalaxis_z"])) == [1e-10, 1e-10] + + +def test_coincident_rows_retain_magnetic_flags() -> None: + structure = Structure( + Lattice.cubic(4), + [{"Fe2+": 0.6}, {"Zn2+": 0.2}], + [[0, 0, 0], [0, 0, 0]], + site_properties={"magmom": [[0, 0, 0], [0, 0, 0]], "magnetic": [True, False]}, + ) + (block,) = CifFile.from_str(str(MCifWriter(structure))).data.values() + labels = dict(zip(block.data["_atom_site_type_symbol"], block.data["_atom_site_label"])) + assert block.data["_atom_site_moment.label"] == [labels["Fe2+"]] + assert list(map(float, block.data["_atom_site_occupancy"])) == [0.6, 0.2] + + +def test_equivalent_constituent_partitions_keep_magnetic_symmetry() -> None: + combined = Structure( + Lattice.cubic(4), + ["Fe", "Fe"], + [[0, 0, 0], [0.5, 0.5, 0.5]], + site_properties={"magmom": [[0, 0, 2]] * 2}, + ) + split = Structure( + combined.lattice, + [{"Fe": 0.5}, {"Fe": 0.5}, "Fe"], + [[0, 0, 0], [0, 0, 0], [0.5, 0.5, 0.5]], + site_properties={"magmom": [[0, 0, 2]] * 3}, + ) + assert str(MCifWriter(split)) == str(MCifWriter(combined)) + + +def test_magmom_spin_axis_preserves_weak_components() -> None: + structure = Structure( + Lattice.cubic(4), + ["Fe"], + [[0, 0, 0]], + site_properties={"magmom": [Magmom([1e-10, 2e-10, 3e-10], saxis=[1, 0, 0])]}, + ) + (block,) = CifFile.from_str(str(MCifWriter(structure))).data.values() + (moment,) = zip( + *[list(map(float, block.data[f"_atom_site_moment.crystalaxis_{axis}"])) for axis in "xyz"] + ) + # SAXIS=x sends local x to global -z, local y to global y, and local z to global x. + np.testing.assert_allclose(moment, [3e-10, 2e-10, -1e-10], atol=1e-25, rtol=0) + + +@pytest.mark.parametrize("order", [(0, 1, 2), (2, 0, 1)]) +def test_ambiguous_coincident_grouping_is_rejected(order) -> None: + positions = [[0, 0, 0], [0.0002, 0, 0], [0.0004, 0, 0]] + structure = Structure( + Lattice.cubic(4), + [{"Fe": 0.2}] * 3, + [positions[index] for index in order], + site_properties={"magmom": [[0, 0, 2]] * 3}, + ) + with pytest.raises(ValueError, match="Ambiguous coincident-site grouping"): + MCifWriter(structure, symprec=1e-3) + + +@pytest.mark.parametrize("difference", [1e-7, 1e-3]) +def test_shared_site_moment_comparison_uses_mag_symprec(difference) -> None: + structure = Structure( + Lattice.cubic(4), + [{"Fe": 0.6}, {"Mn": 0.3}], + [[0, 0, 0]] * 2, + site_properties={"magmom": [[0, 0, 2], [0, 0, 2 + difference]]}, + ) + with pytest.raises(ValueError, match="Unequal magnetic moments"): + MCifWriter(structure, mag_symprec=difference / 10) + assert "_atom_site_moment.label" in str(MCifWriter(structure, mag_symprec=difference * 10)) + + @pytest.mark.parametrize( - ("selection", "message"), + ("flag", "message"), [ - ([], "one selection per site"), - ([("Mn",)], "absent species"), - ([("Fe", "Fe")], "duplicate species"), - ([()], "nonzero moment"), - (["Fe"], "sequence of species"), + (False, "nonmagnetic row must have a zero moment"), + ("true", "boolean flags"), + (None, "boolean flags"), ], ) -def test_rejects_invalid_magnetic_selection(selection, message) -> None: +def test_rejects_invalid_magnetic_flags(flag, message) -> None: structure = Structure( - Lattice.cubic(4), ["Fe"], [[0, 0, 0]], site_properties={"magmom": [[0, 0, 2]]} + Lattice.cubic(4), + ["Fe"], + [[0, 0, 0]], + site_properties={"magmom": [[0, 0, 2]], "magnetic": [flag]}, ) with pytest.raises(ValueError, match=message): - MCifWriter(structure, magnetic_species=selection) + MCifWriter(structure) @pytest.mark.parametrize( @@ -328,6 +437,23 @@ def test_rejects_nonfinite_moments(value) -> None: MCifWriter(structure) +def _constituent_structure(structure: Structure, selections: list[tuple[str, ...]]) -> Structure: + compositions, positions, moments, magnetic = [], [], [], [] + for site, selected in zip(structure, selections, strict=True): + for species, occupancy in site.species.items(): + flag = str(species) in selected + compositions.append({species: occupancy}) + positions.append(site.frac_coords) + moments.append(site.properties["magmom"] if flag else np.zeros(3)) + magnetic.append(flag) + return Structure( + structure.lattice, + compositions, + positions, + site_properties={"magmom": moments, "magnetic": magnetic}, + ) + + def _expand_constituent_rows(text: str) -> list[tuple[str, float, np.ndarray, np.ndarray, bool]]: """Read CIF tokens and expand each constituent independently of the writer.""" parser = CifParser.from_str(text)