feat: preserve constituent data in magnetic CIF exports

CHANGELOG.md

@@ -16,6 +16,20 @@ 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.
+- 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
 
 - Initial public release

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,14 +86,22 @@ 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
 
 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.

docs/file-formats.md

@@ -31,9 +31,40 @@ 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.
+  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

src/spinforge/_cif.py

@@ -1,5 +1,13 @@
 from __future__ import annotations
 
+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, Structure
+from pymatgen.electronic_structure.core import Magmom
+
 _ATOM_SITE_LOOP_HEADER = (
     "",
     "loop_",
@@ -10,3 +18,185 @@
     "_atom_site_fract_z",
     "_atom_site_occupancy",
 )
+
+
+@dataclass(frozen=True)
+class _SiteConstituent:
+    species: Element | Species
+    occupancy: float
+    magnetic: bool
+
+
+@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:
+            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(
+                sorted(
+                    rows, key=lambda row: (_species_key(row.species), row.occupancy, row.magnetic)
+                )
+            )
+        )
+    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 = {
+        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"{_species_token(species)} {_number(_oxidation_state(species))}"
+            for species in sorted(ionic_species, key=_species_key)
+        ],
+    ]