Make MaterialsProjectDFTMixingScheme reproducible; match structures across cell sizes

[kavanase]

Fixes #3113. Fixes #4261.

MaterialsProjectDFTMixingScheme.process_entries can return a different phase diagram for the same input entries from one interpreter run to the next. The variable is PYTHONHASHSEED and nothing else (same entries, order, code etc.). This was originally reported in #3113, and noted by our doped tests now that we are following the new mp-api default of using the mixed GGA/R2SCAN thermo-types (as on the MP website).

Separately, the same grouping code splits a primitive and a conventional cell of one material before StructureMatcher ever sees them, silently discarding the r2SCAN entry in several cases (#4261) — also fixed here.

MWE (#3113)

import hashlib, os
from mp_api.client import MPRester
from pymatgen.entries.mixing_scheme import MaterialsProjectDFTMixingScheme

with MPRester() as mpr:
    entries = mpr.get_entries_in_chemsys(
        "Si-Se", additional_criteria={"thermo_types": ["GGA_GGA+U", "R2SCAN"]}
    )

entries.sort(key=lambda e: str(e.entry_id))   # fixed input order: only PYTHONHASHSEED differs
processed = MaterialsProjectDFTMixingScheme().process_entries(entries)

gs = min((e for e in processed if e.reduced_formula == "SiSe2"), key=lambda e: e.energy_per_atom)
fp = hashlib.md5(
    "".join(f"{e.entry_id}:{e.energy_per_atom:.8f}" for e in sorted(processed, key=lambda e: str(e.entry_id))).encode()
).hexdigest()[:8]
print(f"HASHSEED={os.environ.get('PYTHONHASHSEED', '<unset>')} {len(processed)} entries | "
      f"SiSe2 GS {gs.entry_id} {gs.energy_per_atom:.4f} eV/atom | fp={fp}")

On main — 3 of 8 seeds give a different hull:

$ for s in 0 1 2 3 4 5 6 7; do PYTHONHASHSEED=$s python mwe.py; done
HASHSEED=0  59 entries | SiSe2 GS mp-1250007-r2SCAN  -12.9842 eV/atom | fp=206456a5
HASHSEED=1  59 entries | SiSe2 GS mp-1250007-r2SCAN  -12.9842 eV/atom | fp=206456a5
HASHSEED=2  59 entries | SiSe2 GS mp-1250007-r2SCAN  -12.9842 eV/atom | fp=206456a5
HASHSEED=3  59 entries | SiSe2 GS mp-568264-GGA       -4.9047 eV/atom | fp=66cd2d47   <--
HASHSEED=4  59 entries | SiSe2 GS mp-1250007-r2SCAN  -12.9842 eV/atom | fp=206456a5
HASHSEED=5  59 entries | SiSe2 GS mp-1250007-r2SCAN  -12.9842 eV/atom | fp=206456a5
HASHSEED=6  59 entries | SiSe2 GS mp-568264-GGA       -4.9047 eV/atom | fp=66cd2d47   <--
HASHSEED=7  59 entries | SiSe2 GS mp-568264-GGA       -4.9047 eV/atom | fp=66cd2d47   <--

On this updated branch all eight print fp=206456a5.
Repeats at a fixed seed are identical, so this is hash-iteration order, not (other) randomness; matching the original #3113 report of the count alternating 637/638 "roughly every third run". It is also not degenerate tie-breaking. Over the 58 entries common to both outcomes, formation energies differ by up to 0.254 eV/atom and e_above_hull by 0.007, and the stable set flips (mp-1250007-r2SCAN vs mp-568264-GGA), one entry being silently dropped in each direction.

Origin

1. _filter_and_sort_entries returns list(EntrySet(...)), i.e. list(set(...)).
   ComputedEntry.__hash__ hashes an f-string, and CPython randomises str hashing per process for security (PEP 456), so the set's iteration order and the returned list order changes every run. So the method sorts, but non-deterministically.
2. The groupby(sorted(...)) calls in get_mixing_state_data are stable, so that order survives StructureMatcher.group_structures, which is then greedy (unmatched.pop(0) picks a reference).
3. StructureMatcher.fit is directional (symmetric=False by default) and group_structures always passes its as struct1. Among the five SiSe2 spacegroup-72 structures, 3 of the 10 cross pairs match one way only — e.g. fit(mp-1179443-GGA, mp-568264-r2SCAN) is True, the reverse False. Over all 120 permutations of that pre-group: 3 distinct groupings (60/30/30), vs 1 with symmetric=True (added in this PR). Changing tolerances (e.g. stol, ltol, etc) can't fix this.
4. Whether the GGA ground state keeps its r2SCAN partner flips all(mixing_state_data[mixing_state_data["is_stable_1"]]["entry_id_2"].notna()) in get_adjustments, re-anchoring the whole hull from the r2SCAN scale onto GGA.

Changes/Fixes

• Sorting: New _entry_sort_key, applied at the run_type split in get_mixing_state_data and on _filter_and_sort_entries return. Both are needed: the former is public and callable without the latter, and the latter is what process_entries builds its own return order from. The key must be a robust total order; now using the energy per atom (NaN guarded to avoid silent degradations) and str() on Entry.entry_id.

• symmetric=True in the group_structures call, so the grouping no longer depends on which structure is popped as reference. This requires the concurrent PR I just raised to pymatgen-core here: Add symmetric option to group_structures, pass through to sm.fit() and test pymatgen-core#133 (and hence the pymatgen-core requirement bump).

• entry.structure.copy() before stamping entry_id onto it. .structure returns the live structure, so two entries built from one Structure both end up labelled with the second's entry_id and the first vanishes from the mixing state; in a 3-entry test case that discards every r2SCAN entry. MP's own path is safe (from_dict builds a fresh Structure each time); user code comparing two functionals on one geometry is not.

• entries/mixing_scheme.py -> re-export stub. It was missed by 1c7086f, which stubbed its five siblings, so two identical copies have existed as separate classes. fix: correctly clear energy_adjustments in clean mode (avoid skipped removals) #4615 had to "apply fix in both mixing_scheme modules", and because the two stamp different @module into adjustment cls, re-processing across the two paths with clean=False appends a duplicate adjustment.

• Pre-group on reduced_composition (fixes Mixing scheme fails to mix entries from the same material id even when their structures match #4261). The first groupby in get_mixing_state_data keyed on the full composition, so a 1-atom and a 2-atom cell of one material ('Fe1' != 'Fe2') landed in different pre-groups and never reached group_structures — even though StructureMatcher().fit() would match them correctly. mp-13, the entry reported in Mixing scheme fails to mix entries from the same material id even when their structures match #4261, is exactly this case. Safe to pass on downstream: get_hull_energy_per_atom is scale-invariant (verified), the formula column already used reduced_formula, and the fuzzy-diatomic branch still splits by site count. MWE:

  from pymatgen.analysis.compatibility.mixing_scheme import MaterialsProjectDFTMixingScheme
  from pymatgen.core import Lattice, Structure
  from pymatgen.entries.computed_entries import ComputedStructureEntry
  
  prim = Structure(Lattice.cubic(2.87), ["Fe"], [[0, 0, 0]])  # 1-atom cell
  entries = [
      ComputedStructureEntry(prim, -8.0, parameters={"run_type": "GGA"}, entry_id="Fe-GGA"),
      ComputedStructureEntry(prim * (1, 1, 2), -18.0, parameters={"run_type": "R2SCAN"}, entry_id="Fe-r2SCAN"),
  ]
  df = MaterialsProjectDFTMixingScheme(compat_1=None).get_mixing_state_data(entries)
  print(df[["formula", "spacegroup", "num_sites", "entry_id_1", "entry_id_2"]].to_string(index=False))

  # main -- two rows, each missing its counterpart; the r2SCAN energy is discarded
  formula  spacegroup  num_sites entry_id_1 entry_id_2
       Fe         221          1     Fe-GGA        NaN
       Fe         221          2        NaN  Fe-r2SCAN
  # this PR -- one paired row; process_entries keeps the r2SCAN energy
  formula  spacegroup  num_sites entry_id_1 entry_id_2
       Fe         221          1     Fe-GGA  Fe-r2SCAN

  Unlike the reproducibility changes, this one is a deliberate behaviour change on real data. On 5 test MP chemical systems, 3 are identical and 2 move; in each changed case a GGA entry is newly paired with its r2SCAN partner and replaced by it, as the scheme intends:

  system	entries	GGA(+U) → r2SCAN pairings gained	formation E/atom	e_above_hull
  Li-Fe-O	374 → 374	mp-13 (the Mixing scheme fails to mix entries from the same material id even when their structures match #4261 entry), mp-18905, mp-757614	5/371 changed, max 0.5623 eV	140/371 changed
  Ba-Ti-O	248 → 245	mp-122, mp-12205, mp-2657, mp-2998, mp-882 — the last three's r2SCAN twins were already in the output unpaired, i.e. the same phase counted twice at inconsistent energies	92/243 changed, max 0.0805 eV	10/243 changed
  Si-Se, Zn-S, Na-Cl-O	unchanged	—	—	—

Verification

5 chemical systems (Si-Se, Li-Fe-O, Ba-Ti-O, Zn-S, Na-Cl-O), 972 processed entries, each under PYTHONHASHSEED 0/3/7, unseeded, and with a permuted input list -> exactly one fingerprint per system. New tests, each mutation-tested: removing any one of the source changes here fails exactly one of them. No behaviour change to the existing suite.

Not fixed here

• The all-or-nothing frame choice in get_adjustments: removing a single r2SCAN entry unrelated to the composition of interest re-anchors the whole diagram (from r2SCAN to GGA).

Note for downstream (namely emmet)

emmet-core imports MaterialsProjectDFTMixingScheme from entries.mixing_scheme (emmet/core/io/pymatgen.py:135), which is now deprecated (new import path is analysis.compatibility.mixing_scheme -> should be updated to avoidDeprecationWarnings. cc @esoteric-ephemera @tsmathis

Also relevant to MP directly: emmet-builders' thermo builder calls the same compatibility.process_entries(...) over mixed GGA/R2SCAN entries, so published mixed hulls are one arbitrary sample from this distribution, and a rebuild could shift reported formation energies (in the desired direction, favouring r2SCAN hulls when possible, but worth noting).

cc @rkingsbury as author of the mixing schemce

Checklist

• Google format doc strings added. Check with ruff.
• Type annotations included. Check with mypy.
• Tests added for new features/fixes.
• If applicable, new classes/functions/modules have duecredit @due.dcite decorators to reference relevant papers by DOI (example)

[kavanase]

Tests/linting currently fail as this requires a new pymatgen-core release (with materialsproject/pymatgen-core#133 included)

[shyuep]

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

Solid fix for both #3113 and #4261. The approach is right: canonical (energy, str(entry_id)) sort with NaN handling before greedy grouping, group_structures(symmetric=True) for order-independence, and grouping on reduced_composition so primitive/conventional cells land in one row. The struct.copy() before stamping entry_id also fixes a real caller-mutation bug, and the tests are well-designed (the RecordingStructureMatcher + non-vacuous-assertion guard is a nice touch).

Comments:

1. CI is red purely at the install step — all failures (lint + 3 test matrices) die at "Install dependencies via uv" because pymatgen-core>=2026.8.24 isn't resolvable on PyPI yet. This PR is blocked on pymatgen-core#133 (symmetric kwarg) being merged and released. Re-run CI after that release; no code issues are visible in the failed runs.
2. The conversion of pymatgen/entries/mixing_scheme.py to a import * deprecation stub fires a module-level DeprecationWarning on import — fine, but note downstream packages with -W error filters (atomate2, emmet) will break until they migrate. Consider coordinating.
3. symmetric=True grouping is quadratic-ish vs. the greedy path — likely acceptable, but worth a quick benchmark on a large chemsys (e.g. Li-Fe-P-O with full GGA+R2SCAN thermo types) before merge, since process_entries sits in doped/emmet hot paths.

Recommendation: approve pending pymatgen-core release + green CI.

[rkingsbury]

Wow, thanks for such detailed analysis @kavanase ! I'm no longer actively using the mixing scheme, but I reviewed your changes and they all make sense to me.

I never would have imagined that the sort order of EntrySet would be non-deterministic!

[kavanase]

Thanks @rkingsbury! Yes it was a tricky issue to nail down!

[shyuep]

Automated PR review generated by Claude (posted on behalf of @shyuep). Updated review following new commits since the last automated review (2026-08-25).

Strong PR that fixes two real bugs. The reproducibility fix (#3113) is cleanly done: _entry_sort_key establishes a total order — the NaN guard and str(entry_id) tiebreak are exactly right, since a NaN energy or unorderable entry_id would silently degrade the sort back to arrival (hash) order. Sorting at both get_mixing_state_data and _filter_and_sort_entries, plus group_structures(symmetric=True), closes the remaining order-dependence. Grouping by reduced_composition fixes the primitive/conventional cell split (#4261), and entry.structure.copy() stops mutating caller structures (a nice drive-by fix — the old code attached entry_id to the caller's Structure). CI is green.

Comments:

1. symmetric=True grouping does full pairwise matching rather than greedy first-match, so cost grows for large chemsys pulls. Probably acceptable given groups are pre-partitioned by composition + spacegroup, but worth a quick timing sanity check on a large chemsys (e.g. Li-Fe-P-O with both thermo types).
2. The module move to pymatgen.analysis.compatibility.mixing_scheme with a star-import stub + DeprecationWarning (removal v2027.1) is handled correctly, and the pymatgen-core pin bump to >=2026.8.30 for the symmetric kwarg is properly annotated. Confirm that pymatgen-core 2026.8.30 is actually released before merging, or CI on downstream consumers will break.
3. Behavior note for the changelog readers: deterministic energy-sorted ordering changes which duplicate entry is retained relative to previous (hash-ordered) runs, so processed results may shift vs. previously computed hulls. That is the point of the PR, but downstream users (e.g. doped) should be aware.

Overall: approve. Only item (1) is worth verifying pre-merge.