[breaking] Support mixed-supercell orderings in HeisenbergMapper via a shared parent-cell sublattice definition

[Luguza]

Summary

Refactors HeisenbergMapper so magnetic sublattices are defined once on a common parent cell, and every DFT-relaxed ordering is mapped onto it by structure matching rather than by direct index comparison. Orderings relaxed in different-sized supercells of the same parent can now be combined in one fit.

Closes #4667.

Motivation

The old implementation identified "unique sites" from the geometry of ordered_structures[0] and assumed every other ordering shared that exact indexing and cell. Orderings from a MagneticOrderingsWF-style workflow often don't: they relax into different, commensurate supercells of the same parent, and the mapper had no way to reconcile them. This PR makes an explicit parent cell the single source of truth for sublattice identity and matches each ordering onto it with StructureMatcher(attempt_supercell=True).

How it works

HeisenbergMapper owns one ParentOrdering, whose symmetry orbits are the sublattices, and two or more RelaxedOrderings, each labelled by the parent. Both inherit their magnetic-only structure and coupling graph from MagneticOrdering, which builds the graph with the new SublatticeMinimumDistanceNN.

Steps 1–6 run in the constructor, step 7 in get_exchange() and step 8 in the export methods. Step 6 writes out the linear system with the column names of ex_mat and the keys of ex_params.

Changes

New behaviour

• Sublattices are the parent's symmetry orbits of the magnetic species, found with the nonmagnetic ions present so that site symmetry is correct. The new symprec argument (default 0.01) sets the tolerance.
• Every ordering is matched onto the parent. A match that fails, or that doesn't line up species by species, raises a ValueError instead of mislabelling sites.
• The new matcher argument takes a StructureMatcher with primitive_cell=False and attempt_supercell=True; any other setting raises a ValueError. Loosen its ltol, stol or angle_tol when relaxation distorts an ordering beyond the defaults.
• Coupling graphs are built on the parent geometry (each ordering's magnetic sites at their parent positions), so relaxation no longer changes which pairs couple or which shell they belong to. A bond the parent doesn't have raises a ValueError.
• Shells: for each sublattice pair, sorted bond lengths are split wherever two consecutive lengths differ by more than tol (default 0.02 Å).
• SublatticeMinimumDistanceNN keeps every bond within cutoff, plus the nearest shell (up to the first gap larger than tol) of each sublattice pair with no bond within it. Every sublattice pair with a neighbor within 10 Å therefore couples, also with a cutoff, and cutoff=0 gives one shell per pair; long intra-sublattice bonds no longer drop out behind shorter inter-sublattice ones. A pair with no neighbor that close emits a UserWarning naming it.
• get_exchange() fits by least squares over all orderings instead of solving a square system. It warns when the fit is ill-conditioned or rank deficient and reports an RMS residual in meV per magnetic ion. The residual is only meaningful with substantially more orderings than fitted parameters; otherwise the extra parameters overfit and lower it without improving the fit.
• With n orderings, at most n − 1 couplings are fitted; the longest-range ones are dropped with a UserWarning that names them.
• Magnetic species are pooled over all orderings, so an ion whose moment relaxes to zero keeps its site. Induced moments on nonmagnetic ions are ignored.

Breaking API changes (also in the migration guide in the module docstring)

• HeisenbergMapper(ordered_structures, energies, parent=None, cutoff=0, tol=0.02, symprec=0.01, matcher=None): parent is now third, and a non-Structure passed there raises a TypeError that points to the change.
• get_exchange() returns (ex_params, residual).
• sgraphs, unique_site_ids, wyckoff_ids and nn_interactions are replaced by orderings, coupling_graphs, sublattice_ids, sublattice_wyckoff_symbols and interactions. J labels are now '<i>-<j>-nn', '<i>-<j>-nnn', ... per sublattice pair.
• energies now holds total energies as passed in; per-ion energies are in energies_per_magnetic_ion. The same applies to HeisenbergModel.
• HeisenbergScreener takes RelaxedOrdering objects and exposes screened_orderings.
• HeisenbergModel.from_dict rejects 0.1 dicts with a pointer to recompute, since the old site labels don't map onto parent sublattices.

Deprecations (deadline 2027-08-01)

• estimate_exchange, get_low_energy_orderings and get_mft_temperature: use the shell-resolved couplings from get_exchange(), and a Monte Carlo code such as VAMPIRE for critical temperatures.

Units

• The fitted J are in meV/µB² because they multiply the raw moments, and E0 is in eV per magnetic ion. get_interaction_graph() converts the couplings to meV for unit spins, the convention of VAMPIRE, UppASD and TB2J.

Testing

All 35 tests in tests/analysis/magnetism/test_heisenberg.py pass:

• TestHeisenbergMapperKnownHamiltonian: fits energies generated from a known Hamiltonian and recovers the couplings, including mixed 1×1, 1×2 and 2×2 supercells and relaxed orderings; also covers shell splitting by tol, the nearest-shell floor with and without a cutoff, SublatticeMinimumDistanceNN directly, symprec, the matcher, the residual, the warnings and the errors.
• TestHeisenbergMeanFieldTemperature: the deprecated mean-field path, single- and multi-sublattice.
• TestHeisenbergMapperFullCellSymmetry: nonmagnetic ions splitting sublattices, and equivalent sites sharing one.
• TestHeisenbergMapperZeroMomentIon: an ion whose moment relaxed to zero stays on the magnetic lattice.

Reviewing

heisenberg.py reads top to bottom in the order of the figures: SublatticeMinimumDistanceNN, the three ordering classes (steps 1–3), then HeisenbergMapper (steps 1–8), and finally HeisenbergScreener and HeisenbergModel. The part that most needs a careful look is the alignment check in RelaxedOrdering._set_sublattice_ids.

[shyuep]

Automated PR review generated by Claude (on behalf of @shyuep)

Substantial, well-motivated refactor (parent-cell sublattice definition + StructureMatcher mapping + least-squares fit) with unusually thorough docstrings and new test coverage. Main comments:

1. Attribution: don't comment out the original __author__ block (ncfrey). Either keep both authors or drop the dunder metadata entirely — git history carries attribution.
2. raise SystemExit(...) was carried over from the old code — please use ValueError instead; SystemExit kills interactive sessions and notebooks.
3. This is a significant API break for HeisenbergMapper (attributes like sgraphs, unique_site_ids, wyckoff_ids are gone; estimate_exchange/get_mft_temperature deprecated). A short migration note in the PR description / docs would help downstream users.
4. CI has not run on this branch (first-time contributor — needs a maintainer to approve workflows). Please rebase on master and a maintainer should trigger CI before merge.
5. Default-parent inference from the lowest-energy ordering is a documented footgun; consider emitting a UserWarning when parent=None so users can't miss it.

Nice work overall — the shell-midpoint interaction labeling and pooled magnetic species are clear improvements.

[Luguza]

@shyuep @mkhorton As per the discussion in materialsproject/atomate2#1528, should I restore the vampire caller into pymatgen, because it is written to be specific to this workflow and relevant to local execution?

[shyuep]

Automated PR review generated by Claude (scheduled routine, posted on behalf of @shyuep)

Substantial and well-motivated refactor — defining sublattices once on a parent cell and mapping orderings via StructureMatcher(attempt_supercell=True) is the right architecture for mixed-supercell fits, and the migration guide in the module docstring is excellent. Points to address:

1. Breaking changes without deprecation shims. Constructor arg order changes (parent in third position), get_exchange now returns a tuple, and sgraphs/unique_site_ids/wyckoff_ids/ordered_structures are removed outright. Consider deprecated property aliases for one release cycle where cheap, and mark the PR title [breaking] so the compatibility log picks it up.
2. pymatgen-core submodule pointer bump rides along in the diff — drop it unless this PR actually requires it.
3. Attribution: __author__/__maintainer__ replace @ncfrey entirely; please append rather than overwrite the original authorship.
4. SublatticeMinimumDistanceNN is a nice generalization; the cached nn_graph guard raising before sublattice_ids are set is good defensive design.
5. No CI results available yet (checks likely awaiting approval) — review based on code analysis only.

[Luguza]

Agreed on the title, I'll mark it [breaking] so it shows up in the compatibility log.

On the deprecation shims, I'd rather not add them here, and it's not only about cost. Taking them one at a time:

sgraphs → nn_graphs would be an actively harmful alias. The old graphs were built on the full ordered structures, with every ion as a node; the new ones are built on the magnetic-only structure. The node sets differ, so any survivingsgraph.structure[i]or index-based loop would keep running and silently address a different site. Today that code fails loudly with an AttributeError on the first line that touches sgraphs. An alias would trade that for wrong numbers.

unique_site_ids → sublattice_ids isn't an alias at all. The old attribute was a dict[tuple[int, ...], int] mapping tuples of equivalent supercell indices to an id; the new one is a list[list[int]] giving the parent-sublattice label of each magnetic site in each ordering. Bridging them means reconstructing the per-supercell orbit map — which is precisely the single-supercell assumption this PR exists to remove. wyckoff_ids → sublattice_wyckoff_symbols does happen to keep its dict[int, str] shape, but the ids underneath now index parent orbits rather than supercell orbits, so a bare alias would hand back keys that don't mean what the caller thinks.

get_exchange is breaking regardless of the return type. The {"<J>": ...} fallback for under-determined systems is gone and now raises ValueError, and the solve became a least-squares fit over all orderings rather than an exact inverse — so numbers move even for callers who never touch residual. Keeping the dict return would hide the signature change without making any of that safe.

There's also one break no shim can cover: energies keeps its name but now holds energy per magnetic ion instead of total energy. That's a silent numeric change in an attribute that still exists, which to me is the argument that a partial soft landing here is worse than a clean one — it lets callers get most of the way through a migration while quietly wrong.

Given that, I think the honest version is a hard break with a good migration path rather than a half-compatible surface. The module docstring already carries a migration guide covering every renamed and removed attribute, the constructor order, the get_exchange change, and HeisenbergScreener — worth a look if it didn't come through in the diff view. Happy to expand it, add a CHANGES.md entry, or move it somewhere more discoverable if you'd prefer.

[shyuep]

Automated PR review generated by Claude (posted on behalf of @shyuep; scheduled review routine).

This is a well-executed refactor. The test suite is a highlight: round-tripping against a synthetic Hamiltonian with known J's, ground-truth energies computed independently of the mapper's own graph, and explicit mixed-supercell (1x1/1x2/2x2) recovery. The migration notes in the module docstring are thorough.

Comments:

1. Legacy HeisenbergModel deserialization is broken. from_dict now requires magnetic_structures, nn_graphs, sublattice_ids, sublattice_wyckoff_symbols, and residual. Any previously stored MSON dict (with sgraphs, unique_site_ids, wyckoff_ids, javg) will KeyError. You handle old ex_mat string forms, so partial back-compat is clearly intended — either extend from_dict to accept the old key set (even if lossy) or raise a clear "serialized with HeisenbergModel <0.2" error.
2. tol default silently changed from 0.02 to 0.05 (DEFAULT_TOL). Not mentioned in the migration section — please add it, since it changes shell-merging behavior for existing callers.
3. Positional-arg shift: HeisenbergMapper(structs, energies, 5.0) previously meant cutoff=5.0; now 5.0 lands in parent. It fails loudly in _nonmagnetic, but a friendlier isinstance(parent, Structure) check with a pointed error message would help migrating users.
4. self.residual is only set in get_exchange — initialize to None in __init__ alongside ex_params so attribute access before fitting doesn't raise AttributeError (the class docstring implies it's an attribute).
5. Explicit parent is silently primitivized via get_primitive_structure(). Worth a docstring note, since a user supplying a deliberate supercell parent may not expect orbit analysis on the primitive cell.
6. The pymatgen-core submodule bump is included — confirm it's intentional and points at merged main (it matches merged #130, so likely fine, just flagging).
7. CI has not run on the latest commits ("no checks reported" on the branch). A maintainer needs to approve/trigger workflows before this can be evaluated for merge.

The deprecations of estimate_exchange/get_mft_temperature with concrete alternatives (shell-resolved fit; MC solvers) are the right call, and pinning the flawed multi-sublattice MFT formula in a test rather than silently "fixing" it is good hygiene.

[shyuep]

@mkhorton I suggest you review this since I am not familiar with this code.

[shyuep]

Automated PR review (generated by Claude on behalf of @shyuep).

This is a well-motivated redesign — defining sublattices once on an explicit parent cell and mapping each relaxed ordering via StructureMatcher is the right fix for #4667, and the move from a square solve to least squares over all orderings (with rank-deficiency warning + RMS residual) is a clear improvement. The test coverage (known-Hamiltonian recovery, Wyckoff-orbit sublattices, non-supercell rejection) is thorough, and deprecating estimate_exchange/get_mft_temperature rather than removing them is the right call.

Requested changes:

1. Revert the pymatgen-core submodule pointer bump — it's unrelated to this PR and appears accidental (also present in Add whole motif flip strategy to MagneticStructureEnumerator #4699).
2. CI has not run on this branch (needs maintainer approval for first-time contributor workflows); once the submodule change is reverted, we'll approve the run.
3. HeisenbergScreener signature is now (orderings, screen=False) rather than (structures, energies, screen) — that's a second breaking public API beyond HeisenbergMapper itself. Please call it out explicitly in the PR description/CHANGES so downstream users (e.g. atomate2 magnetism workflows, which also touch sgraphs/unique_site_ids) can migrate.
4. Minor: SystemExit → ValueError is appreciated.

[shyuep]

🤖 Automated PR review generated by Claude on behalf of @shyuep. Not a human review — treat as a first pass.

The physics motivation is right, and the core insight — define sublattices once on a paramagnetic parent cell so orderings in different supercells can be fitted together — is a real improvement over deriving them from the orderings themselves. SublatticeMinimumDistanceNN fixes a genuine defect (sublattice-blind MinimumDistanceNN drops long intra-sublattice bonds entirely), and replacing the exactly-determined solve with a least-squares fit plus a reported residual is clearly better practice.

That said, at +1427/−699 in one module with a breaking public surface and no CI results reported (0 checks), this is not mergeable as-is.

Blocking

1. No CI. Please push/rebase so the suite runs. Nothing below can be confirmed without it.

2. Positional-argument insertion is a silent breakage. HeisenbergMapper(structures, energies, parent, cutoff, tol) puts parent where cutoff used to be. Existing positional calls will pass a float as parent and get a confusing downstream failure rather than a TypeError. Make parent keyword-only, or append it last, or add a type check on the third argument that raises a clear deprecation error.

3. energies changes meaning under the same name (total → per magnetic ion). Silently returning different numbers from an unchanged attribute is the hardest kind of change to detect downstream. Rename to energies_per_magnetic_ion and either drop energies or keep it as the total with a deprecation warning.

4. get_exchange return type changed from ex_params to (ex_params, residual). Every existing caller breaks. Consider attaching residual only to self.residual / HeisenbergModel and leaving the return type alone — the information is already stored on both.

5. Removals with no deprecation path: sgraphs, unique_site_ids, wyckoff_ids, ordered_structures, HeisenbergScreener.screened_structures / screened_energies, and the {"<J>": ...} under-determined fallback. monty.dev.deprecated is already imported in this file — please use it for at least one release cycle. atomate2 and emmet both touch this class.

6. HeisenbergModel.from_dict round-trip. The diff rewrites from_dict (−64/+20). Serialized HeisenbergModel documents exist in production MongoDB stores without a residual field and with the old ex_params shape. Please add an explicit test deserializing a pre-0.2 dict.

7. __maintainer__ / __email__. Adding yourself as module maintainer is a decision for the pymatgen maintainers, not the PR. __author__ is the appropriate field for a contribution of this size; please revert the maintainer change and let @shyuep decide.

Should fix

8. Default tol 0.02 → 0.05 Å silently changes which distances merge into a shell, and therefore the fitted J values, for every existing user. If the new default is better justified, say why in the PR body; otherwise keep 0.02 and let users opt in.

9. Split the PR. As one change this is very hard to review and impossible to bisect. Natural seams: (a) SublatticeMinimumDistanceNN + tests, (b) parent-cell sublattice refactor, (c) least-squares fit + residual reporting, (d) deprecations. (a) and (c) are independently valuable and could land first.

10. DISTANCE_ROUND_DECIMALS = 2 interacts with tol. Distances are rounded to 0.01 Å and then compared against tol = 0.05 Å. Two constants controlling the same shell-merging decision invites inconsistency; consider deriving the rounding from tol or documenting the relationship.

11. _interaction_label falls back to min(labels, key=|dist - self.dists[label]|) — an off-parent-lattice bond is assigned to its nearest shell with only a logger.warning. For a relaxed ordering that has genuinely drifted, this silently contaminates the fit. Given the PR already routes ill-conditioning through warnings.warn, this one should be a UserWarning too.

Positive notes

• The "Migrating from HeisenbergMapper 0.1" section is unusually thorough — please keep it and promote it to the release notes.
• Reporting rank deficiency / ill-conditioning as a UserWarning rather than a logger call is the right call.
• Deprecating get_mft_temperature in favour of a proper MC solver is correct; mean-field T_c overestimates badly for low-dimensional magnets.
• Test file grows from 95 to 477 lines — good.

[shyuep]

Automated PR review generated by Claude (sanity-check before acting on it)

The design direction is good: defining sublattices once on a shared parent cell and mapping each relaxed ordering onto it via StructureMatcher(attempt_supercell=True) is the right fix for #4667, and the least-squares get_exchange (replacing the square-system solve, adding a rank-deficiency warning and RMS residual) is a real improvement over silently discarding surplus orderings. Deprecating estimate_exchange/get_mft_temperature with a dated sunset is handled the right way.

However, this cannot merge as-is: CI is failing because the file doesn't parse. src/pymatgen/analysis/magnetism/heisenberg.py around line 856 has a stray fragment appended directly onto the closing } of the ex_params dict comprehension, with no line break:

        }if matched_parent[i].specie.symbol != self.structure[i].specie.symbol: raise ValueError(...

Ruff reports this as invalid-syntax (Expected else, found :`` / Expected ,`, found name` / `Expected `)`, found newline`) — this is a genuine syntax error, not a style nit, and it's why lint fails and every test matrix job (`macos`, `ubuntu`, `windows`) fails too (the module can't be imported). `matched_parent` and `i` also aren't otherwise bound in this scope, and the `raise ValueError(...)` looks like an unfinished edit (literal `...`) rather than intended code — looks like a merge/edit artifact that needs to be either completed properly or removed before anything else here can be evaluated.

Given the size of this diff (1427/-699 across 2 files) and that it's marked [breaking], I'd hold off on a deeper design/correctness pass until the file at least imports and CI is green, since right now none of the new test classes described in the PR body (TestHeisenbergMapperKnownHamiltonian, etc.) can actually run.

[Luguza]

Thanks for the detailed reviews. I pushed a round of fixes; here's where each point stands.

CI failure. The syntax error at heisenberg.py:856 was a stray fragment left over from an edit and is removed. Ruff and the full test_heisenberg.py now pass locally.

energies (correction to my earlier comment). I claimed energies changed from total energy to energy per magnetic ion. That was wrong: in 0.1, HeisenbergScreener._do_cleanup already divided by the number of magnetic ions, so mapper.energies and HeisenbergModel.energies were always per ion, despite the constructor taking total energies. I've now made the name mean what it says:

• mapper.energies / HeisenbergModel.energies hold the total energies of the screened orderings.
• mapper.energies_per_magnetic_ion / HeisenbergModel.energies_per_magnetic_ion hold the per-ion energies the fit uses, counted over the magnetic species pooled across orderings.
• The migration guide states the switch explicitly.

Positional parent. A non-Structure third argument raises a TypeError that points to the signature change and the migration guide, so HeisenbergMapper(structs, energies, 5.0) fails immediately with a clear message. I think this is what we want users to experience, so that they are made aware of the code changes under the hood.

get_exchange return type. I'd like to keep (ex_params, residual). The result changes for existing callers anyway (least squares over all orderings instead of an exact solve, no <J> fallback), so keeping the old return type wouldn't make old call sites safe. The residual is the measure of whether the J's can be trusted, so it should come back with them rather than sit only on an attribute.

Deprecation shims / removed attributes. I'd still rather not add them, for the reasons in my earlier comment: sgraphs, unique_site_ids and wyckoff_ids have no faithful equivalent in the new model, and aliases would turn today's loud AttributeErrors into silently wrong indices. The migration guide in the module docstring covers every renamed or removed attribute.

from_dict on pre-0.2 dicts. Covered by test_from_dict_rejects_legacy_serialization: old dicts raise a ValueError that names the version and explains that the model must be recomputed, instead of a bare KeyError.

tol default and DISTANCE_ROUND_DECIMALS. The reasoning for 0.05 Å is now a comment at DEFAULT_TOL: relaxation spreads symmetry-equivalent bond lengths by a few hundredths of an Å, and at 0.02 one physical shell was split into several J columns. How tol relates to the 0.01 Å distance rounding is documented next to DISTANCE_ROUND_DECIMALS, and HeisenbergMapper now raises for tol < 0.01, where it would have no effect.

_interaction_label warning. A bond whose sublattice pair isn't in the parent graph now emits a UserWarning instead of a log line, consistent with the rank-deficiency and ill-conditioning warnings. Assigning each bond to the nearest shell of its pair is intentional (shell boundaries sit at the midpoints between parent distances, see the docstring); the warning covers the case where no shell exists for that pair at all.

Maintainer field. @ncfrey hasn't been active on this module for a while, and I'm happy to take over maintaining it, which is why I'm listed. If you'd rather decide that separately, I'll revert the field.

Splitting the PR. I'd prefer to keep it as one. SublatticeMinimumDistanceNN, the shell labelling and the per-ion least-squares fit all use the parent-cell sublattice labels, so landing them separately would mean temporarily adapting them to the old unique_site_ids, which is exactly the single-supercell assumption this PR removes. If it helps review, I can walk through the diff in commit-sized pieces here.

[shyuep]

Automated PR review generated by Claude (posted under @shyuep's account) — please sanity-check before acting.

Good: Defining sublattices once on a parent cell and mapping each ordering onto it via StructureMatcher is the right decoupling, and the pooled magnetic-species set neatly keeps zero-moment ions on the lattice. The alignment check after get_s2_like_s1 (length + species-by-species) turns a silent mislabeling into a clear error. Module docstring migration guide is thorough, PR is marked [breaking], and CI is green on lint and all test matrices at e610008.

Substantive

• _exchange_columns truncates to len(orderings) + 1 columns, silently dropping the longest-range J's (sorted by distance) when there are fewer orderings than interactions. Yet nn_interactions/_interaction_label still contain those dropped labels, so _interaction_label can return a label with no column in ex_mat. Please confirm _build_exchange_mat/get_interaction_graph handle labels absent from j_columns, and emit a UserWarning when interactions are dropped.
• With truncation the system is exactly determined whenever interactions ≥ orderings−1, so the "surplus orderings average out noise" claim only holds when the lattice has few interactions. Worth stating in the get_exchange docstring.
• parent=None inference from orderings[0].structure is documented as unsafe after relaxation; since this is now the default path, consider whether a warning (currently emitted) should be an error when the inferred primitive cell has lower symmetry than expected (e.g. compare Wyckoff multiplicity vs. the ordering's).
• Default tol changed 0.02 → 0.05 Å and the constructor now rejects tol < 0.01. This is a silent numeric change for existing callers passing nothing; the docstring covers it.

Minor

• RelaxedOrdering.set_sublattice_ids: zip(..., strict=False) after the explicit length check is fine, but strict=True documents intent.
• self.orderings = self.parent = None chained assignment is easy to misread; two lines are clearer.
• Author/maintainer email is a student address; confirm this is intended in __email__.
• Only the first ~1000 lines of the diff (heisenberg.py) were reviewed in detail; the 491-line test additions were skimmed, not analyzed.

[Luguza]

@shyuep
Thanks for the last review. This round has six commits. They respond to it and, beyond that, try to make the PR easier to review: the module is shorter, and its methods appear in the order they run. I've also included two figures below, showing the classes and the mapping step by step, to make a human review easier. After all these rounds of improvements, I hope the PR is now ready for one.

Where the last review's points stand

• Dropped interactions: the mapper now warns with a UserWarning that lists the labels it leaves out (c7ba4a5). Labels without a column in ex_mat are handled: get_interaction_graph looks them up with ex_params.get(label, 0) and skips them, so they don't appear as edges.
• "Surplus orderings average out noise": correct, this applies only when the parent has fewer than n − 1 shells. With the default cutoff=0 that is common, because only the nearest shell of each sublattice pair is kept. With a large cutoff the system becomes square.
• tol default 0.05 and the minimum-tol check: both are gone. Coupling graphs are now built on the parent geometry (4ce9b4d), so every bond length is an exact parent distance and needs no tolerance for relaxation. DEFAULT_TOL is back to 0.02, and tol now only means "a gap larger than this between consecutive bond lengths of a sublattice pair starts a new shell" (3404426).
• parent=None: I kept a warning rather than an error. Checking Wyckoff multiplicities can't tell a relaxation-lowered symmetry from a genuinely low-symmetry material, so an error would reject valid input. Instead there is a new symprec argument for parents distorted by an antiferromagnetic relaxation. Happy to change this if you'd prefer an error.
• strict=True and the chained assignment: done.

What changed in this round

Commit	Change
e04b973	Rename nn_graph to coupling_graph and nn_interactions to interactions
c7ba4a5	Warn when interactions are left out of the fit
4ce9b4d	Build coupling graphs on the parent geometry, so relaxation no longer decides which pairs couple; a bond the parent lacks raises a ValueError
3404426	Split shells at gaps larger than tol; add symprec
8657bfd	Build the parent before the orderings, so a RelaxedOrdering is fully labelled from its constructor (removes set_parent and the half-initialized state)
cb4154c	Simplify for review: shorter docstrings, one-line helpers inlined, methods ordered as they run; deprecate get_low_energy_orderings, which only estimate_exchange uses

A guide to reviewing heisenberg.py

The module now reads top to bottom in the order things happen, so the two figures below double as a map of the file.

It starts with the classes. HeisenbergMapper owns one ParentOrdering, which defines the sublattices, and two or more RelaxedOrderings, each labelled by the parent. Both ordering classes inherit their magnetic-only structure and coupling graph from MagneticOrdering, and build that graph with the new SublatticeMinimumDistanceNN:

Reading SublatticeMinimumDistanceNN first, then the three ordering classes, covers everything below the mapper. The part I'd most like a second pair of eyes on is the alignment check in RelaxedOrdering.set_sublattice_ids, after get_s2_like_s1, because a wrong match there would shift every label silently.

The mapper itself then runs these eight steps. Each step names the attributes it produces on the right, and step 6 writes out the linear system with the column names of ex_mat and the keys of ex_params:

In the code, steps 1–6 are __init__ → _initialize_orderings → _set_interactions → _interaction_label → _build_exchange_mat, step 7 is get_exchange and step 8 is get_interaction_graph / get_heisenberg_model. After these come the deprecated methods, then HeisenbergScreener (step 4) and HeisenbergModel with its serialization, which now rejects 0.1 dicts.

For the tests, in tests/analysis/magnetism/test_heisenberg.py, TestHeisenbergMapperKnownHamiltonian is the core: it fits energies generated from a known Hamiltonian, including mixed 1×1, 1×2 and 2×2 supercells and relaxed orderings, and checks that the J's come back. All 29 tests pass locally.

Breaking changes

These are unchanged from earlier rounds and are listed in the module docstring's migration guide: parent is the third positional argument, get_exchange() returns (ex_params, residual), sgraphs/unique_site_ids/wyckoff_ids are replaced by per-ordering attributes, and 0.1 HeisenbergModel dicts are rejected with a pointer to recompute.

[JaGeo]

@shyuep This PR is part of a larger workflow that will be integrated in the atomate2 infrastructure. Thus, it would be of high interest that it is integrated into pymatgen.

[shyuep]

Automated PR review generated by Claude (posted under @shyuep's account) — sanity-check before acting. Based on code and CI only; nothing was run. This is a 117k-char diff; I read the first ~900 lines of heisenberg.py (module docstring, MagneticOrdering/ParentOrdering/RelaxedOrdering, SublatticeMinimumDistanceNN, HeisenbergMapper.__init__, _set_interactions, _build_exchange_mat, start of get_exchange) and did not review the remainder or the 594 added test lines. CI: 4 pass, 1 skipped.

Good: Defining sublattices once on a parent cell and mapping each ordering by structure matching fixes a real limitation (index/cell assumptions across orderings). The migration notes in the module docstring are excellent, parent type-check gives a useful error for the shifted positional args, and the unaligned-match guard in RelaxedOrdering.set_sublattice_ids avoids silently wrong labels.

Substantive

• _build_exchange_mat truncates to n_orderings - 1 J columns plus E0, so whenever truncation applies the system is exactly determined and the reported residual is ~0 by construction. It is then not a fit-quality metric; only the rank-deficient/overdetermined cases make it informative. The docstring/migration text ("least-squares… more orderings than parameters is now useful") overstates this. Consider documenting, or keeping all columns and warning on rank deficiency.
• SublatticeMinimumDistanceNN(self.sublattice_ids) is constructed with its default relative tol=0.1, while the mapper's tol is an absolute Å gap used only for shell labelling. With cutoff=0, shells beyond 10% of the nearest per sublattice pair are never in the graph, so "nnn" labels cannot appear, and the user-facing tol doesn't control neighbour selection as the docstring suggests. Pass through or document.
• StructureMatcher(primitive_cell=False, attempt_supercell=True) with default tolerances picks one match; for strongly relaxed orderings the match may be wrong or fail. The size/species alignment check helps, but failure mode is a hard ValueError for the whole set. Expose ltol/stol/angle_tol (or accept per-ordering failures with a warning).
• site.specie.symbol is used throughout; this raises on disordered sites. Fine if unsupported, but say so.
• Breaking API (parent inserted positionally, get_exchange return signature, energies semantics changed, sgraphs/nn_interactions removed) with no deprecation shim for the attribute renames. Also confirm as_dict/from_dict of HeisenbergModel handles previously serialized models (I did not reach that part).

Nits

• ~1.4k changed source lines mixing new classes, screener rewrite and deprecations; consider splitting for reviewability.
• __author__/__maintainer__ email edits are unrelated to the change.

[Luguza]

@shyuep Thanks for the review. Two commits address it; point by point:

Residual when the system is square. Agreed, the docs oversold it. The get_exchange docstring and the migration guide now say the residual is only meaningful with substantially more orderings than fitted parameters; otherwise extra parameters overfit and lower it without improving the fit. I kept the truncation to n − 1 interactions, since fitting more columns than orderings would only trade a ~0 residual for a rank-deficient fit.

SublatticeMinimumDistanceNN and tol. You were right that the strategy's own relative 10% window was independent of the mapper's tol, and undocumented. (One correction: "nnn" could still appear without a cutoff, when a gap larger than tol fell inside that window, which was also not intended.) The neighbour selection is now one rule, using the mapper's cutoff and tol:

• every bond within cutoff couples;
• a sublattice pair with no bond within cutoff still couples through its nearest shell, which ends at the first gap larger than tol (Å), the same rule that splits shells in _set_interactions.

So cutoff=0 gives exactly one interaction per sublattice pair, and a cutoff can no longer leave a pair uncoupled, which the old MinimumDistanceNN(cutoff=...) path could. Both paths now use the same strategy, and the 10% window is gone. The migration guide notes the change.

StructureMatcher tolerances. I'd rather not expose them in this PR. The defaults (ltol=0.2, stol=0.3, angle_tol=5) are already loose, and a wrong match fails loudly at the site-by-site alignment check rather than mislabelling sites. Relaxation-distorted parents are handled by symprec or by passing the unrelaxed parent. Adding matcher kwargs later would be non-breaking if someone needs them.

Disordered sites. Not reachable: every structure goes through CollinearMagneticStructureAnalyzer first, which raises NotImplementedError for disordered structures, so specie.symbol only ever sees ordered sites.

Breaking API / HeisenbergModel deserialization. Unchanged from earlier rounds: the breaking changes are listed in the module docstring's migration guide (see my earlier comment on why I'd rather not add shims), and from_dict rejects 0.1 dicts with a clear error that points to recomputing.

Splitting the PR / author block. The pieces depend on each other (the screener and deprecations change because the mapper's data model changes), so I'd prefer to keep it together; the figures below and the review guide in my previous comment walk through the module in execution order. The author block appends to @ncfrey's entry instead of replacing it, as requested in an earlier round.

Updated figures

SublatticeMinimumDistanceNN now derives from NearNeighbors and takes the mapper's cutoff and tol:

Step 5 now states which bonds are taken, and step 7 the caveat on the residual:

The PR description is updated to match. All 31 tests pass locally.

[shyuep]

Automated PR review generated by Claude (posted via @shyuep's account); sanity-check before acting.

Good: Defining sublattices once on a parent cell and mapping each ordering onto it is a sound design, and lets orderings in different supercells share a single J set. The StructureMatcher mapping has an explicit alignment guard, the migration notes in the module docstring are thorough, and the TypeError guard on the shifted positional parent argument is a nice touch. Lint and tests are green.

Substantive:

• Cost of SublatticeMinimumDistanceNN.get_nn_info: it calls structure.get_neighbors(site, max(cutoff, 10 Å)) for every site and sorts the result, even when the nearest shell is within about 3 Å. For large supercells this will dominate runtime. Consider growing the search radius adaptively until every neighbor sublattice has a shell.
• Silent uncoupling: a sublattice pair with no neighbor within MAX_SEARCH_DIST is dropped without notice, which yields a model missing that J. Emit a warnings.warn naming the pair.
• Parent inference from a relaxed ordering: RelaxedOrdering._set_sublattice_ids runs StructureMatcher with default tolerances against the (possibly relaxed) parent. Strongly relaxed orderings will hit the "not a supercell of the parent" ValueError with no way to loosen ltol/stol/angle_tol. Expose them (or a matcher argument).
• The PR is about 1,800 lines of source diff with several breaking changes (constructor order, get_exchange return type, removed attributes). Splitting the API changes from the fit changes would make review and bisecting easier.

Nits: the ex_params docstring mixes units (J in meV/μB², E0 in eV per ion), so it is worth a unit-table in get_exchange.

[Luguza]

@shyuep Thanks for the review. Two commits address it; point by point:

Cost of SublatticeMinimumDistanceNN.get_nn_info. The graph is built on the magnetic-only cell, so the fixed 10 Å search stays cheap: building the whole coupling graph takes 0.04 s for 32 magnetic sites and 0.32 s for 256 (fcc lattice, once per ordering). I'd rather keep the fixed radius than add an adaptive search.

Silent uncoupling. Agreed. A sublattice pair with no neighbor within the search distance now raises a UserWarning that names the pair.

StructureMatcher tolerances. Agreed on reflection. In atomate2's MagneticOrderingsMaker, which feeds this module, the parent is whatever structure the workflow starts from. In practice that is usually a Materials Project structure, relaxed with ferromagnetic initial moments, and each ordering is then relaxed separately from it under its own moments. That can distort an ordering beyond the default tolerances. HeisenbergMapper now takes a matcher. It must have primitive_cell=False and attempt_supercell=True (otherwise it raises a ValueError), and the mapping errors point to it.

Splitting the PR. Unchanged from my previous reply: the pieces depend on each other.

Units of ex_params. The mix is inherited: 0.1 also returned J in meV and E0 in eV, but documented both as meV/atom. Both units are now stated in the class docstring and in the get_exchange return value.

[shyuep]

Automated PR review generated by Claude (posted from @shyuep's account). Based on code and CI status only, nothing was run. Please sanity-check before acting.

Good: Making a paramagnetic parent cell the single source of truth for sublattices (orbits of SpacegroupAnalyzer on the nonmagnetic parent) and mapping each ordering onto it with StructureMatcher(attempt_supercell=True) is the right design. Couplings come from the parent geometry, so relaxation no longer decides shell membership. The aligned check after matching guards against silently shifted labels. The least-squares fit with residual, ill-conditioning and rank-deficiency warnings is a clear improvement over the exactly-determined inv(H). The migration docstring and the legacy from_dict error are helpful. CI: lint + 3 test jobs green.

Substantive:

1. get_exchange raises for a single J. With one sublattice and one shell (the textbook FM/AFM case) col_names is [E0, J], so len(col_names) < 3 raises ValueError, and get_heisenberg_model() fails too. [E0, J] is a well-posed 2-parameter fit with ≥2 orderings; the old <J> fallback is gone but nothing replaces it. Consider allowing 1 J (or at least say in the error how to get a usable model).
2. Deprecated get_mft_temperature silently changes meaning. It now reads ex_params, whose J_ij are in meV/μB² (not meV per ion), and the k_boltzmann and omega algebra is unchanged. T_c scales with the moment²: wrong by ~m² for m ≠ 1 μB, with no error. Either convert with the mean |m|² before use, or raise rather than return a number with the wrong units.
3. Mixed units in one dict. ex_params holds J_ij in meV/μB² and E0 in eV per ion. It is documented, but fragile for downstream consumers (e.g. anything that sums or plots the values). Consider moving E0 out of ex_params, or at least keeping units uniform.
4. Cost of the neighbor search. SublatticeMinimumDistanceNN.get_nn_info calls structure.get_neighbors(site, max(cutoff, 10 Å)) for every site of every ordering's supercell, and RelaxedOrdering runs get_s2_like_s1 with attempt_supercell per ordering. For the 100+ magnetic-site supercells that motivate this PR this may be slow. I don't see a test or note on scaling. Consider a cell-list/get_all_neighbors batch call, and a note on expected cost.
5. Truncation of J columns. When n_orderings − 1 < n_shells, the longest-range shells are dropped with a warning. Reasonable, but interactions/dists still list the dropped shells, and get_interaction_graph silently omits those bonds. Ensure downstream users can tell which shells are in ex_params.

Nits:

• HeisenbergModel.as_dict calls self.ex_mat.to_dict(), which raises on the default ex_mat=None; guard for None.
• DEFAULT_TOL / DIST_EPSILON interplay: tol splits shells but _interaction_label matches with DIST_EPSILON = 1e-6 against the parent's shell starts. For relaxed supercells with a small numerical difference between dist and the parent bond length, a bond just under a shell start would hit ValueError("No interaction ..."). Confirm the ideal_magnetic_structure path makes this exact.
• The PR is large (+~1.3k/−0.5k in heisenberg.py) and breaking. A short deprecation shim for the old positional cutoff, tol order would ease migration, though the TypeError on a non-Structure parent is a good start.

[Luguza]

@shyuep Thanks for the review. One commit addresses it.

However, the last few automated rounds have mostly repeated points that were already answered, with no change in between. I take that as a sign the PR is ready for a human review.

get_exchange raises for a single J. Agreed. E0 and one J are fixed by any two orderings that give J different coefficients, such as FM and AFM. get_exchange now needs one J instead of two. The test for this case now checks that E0 and J are recovered exactly on the square-lattice FM/AFM case.

get_mft_temperature units. This is not a regression. On master the exchange matrix already used the raw moments (s_i * s_j), so J was already in meV/μB², and the multi-sublattice branch read the same ex_params. The method is deprecated and its behavior is unchanged.

Mixed units in ex_params. Same answer as last round. The J's multiply the raw moments so that they don't absorb the differences in moment size between orderings. E0 has the same unit as the energies it is fitted to. Both units are stated in the class docstring and in the get_exchange return value.

Cost of the neighbor search. Same answer as last round. Building the coupling graph takes 0.04 s for 32 magnetic sites and 0.32 s for 256. The StructureMatcher call runs once per ordering, and it is the mapping step itself.

Truncation. dists and interactions must keep the dropped shells. Without them, _interaction_label would put long bonds into a shorter shell. The warning names the dropped shells, and the keys of ex_params show which ones were fitted.

Nits. as_dict with ex_mat=None behaves as on master; the other attributes that default to None would fail the same way. Bonds are measured on ideal_magnetic_structure, so their lengths match the parent's shells up to float noise. The positional-argument change already raises a TypeError that names the cause.

[shyuep]

Automated PR review generated by Claude (posted from @shyuep's account). Based on code and CI status only, nothing was run. Please sanity-check before acting.

Follow-up on 5e6a6c4 ("Let get_exchange fit a single J"). This resolves item 1 of my previous review: with one shell, [E0, J] is now fitted, and the ValueError is limited to the no-interaction case with an accurate message. test_single_interaction_is_fitted recovers E0 and J exactly. CI: lint + 3 test jobs green, trigger_atomate2_ci skipped.

Caveat: with one J and two orderings the fit is exactly determined, so residual is identically 0 and says nothing about model quality. The docstring warns about overfitting in general; consider a UserWarning when n_orderings == n_params.

Still open from my previous review (not touched by this commit): item 2 (deprecated get_mft_temperature now reads ex_params in meV/μB², so T_c is off by ~m² for m ≠ 1 μB), item 3 (E0 in eV per ion and J in meV/μB² in one dict), item 4 (neighbor-search cost for large supercells), item 5 (dropped shells still listed in interactions/dists).