Deep, code-cited reference. See the readable summary for the public version. This page is the full provenance; a Verification addendum (run wf_7fbcd118-848, 2026-06-08) at the bottom lists every point where the code differs from this text.

Relaxed tracks replaced 2026-06-15 (hydrostatic + full relax). The site now shows three ordering modes: static (sp, unchanged) plus two relaxed modes from the dedicated relax-cell campaign — hydro (volume-only, cell shape frozen) and full (cell shape + volume, ASE FrechetCellFilter). These replace the old fixed-cell relax track (Track B) described below. The MLIP relaxed energies come from relax_export/relax_structures.parquet (59 potentials × 834 structures × 2 modes), regenerated into ordering_rollup.csv by tools/regen_ordering_from_relax_export.py. The metric definitions (§3) and the energy_above_hull ranking argument are unchanged — only the MLIP geometry/energy source differs. Sections describing mlip_relax/Track B and build_rollups document the superseded fixed-cell pipeline; m3gnet_MP-2021.2.8-EFS was not re-run and is absent from the relaxed tracks (still present in sp).

Polymorph Stability-Ordering Study — Full Method & Provenance

Benchmark: mlip-elemental-benchmark (60 ASSYST MLIP potentials × 834 Materials-Project elemental crystals) Question answered: For each chemical element, does an MLIP reproduce the DFT energy ranking of that element's polymorphs, and does it pick the correct ground-state polymorph? DB: dft_surrogate_mlip_elemental on db-host:db-port Repo / dev tree: github.com/ligerzero-ai/mlip-elemental-benchmark (private) / .../mlip-elemental-benchmark Date documented: 2026-06-08 · Author: Han

0. TL;DR — what "the stability study" is

The stability study is the ordering rollup. It has two evaluation modes ("tracks"):

For each (potential, element) we collect all that element's polymorphs, rank them by MLIP energy/atom, rank them by DFT, and score three metrics: Spearman ρ, ground-state-correct (does the MLIP's lowest-energy polymorph match DFT's?), and dE-MAE (relative-energy error).

The output is results/ordering_rollup.csv — one row per (potential, element, track).

1. Provenance chain (end-to-end, file by file)

  Materials Project API
        │  ingest_mp_elementals.py        (login node, needs internet + MP_API_KEY)
        ▼
  structure  table  ───────────────────────────────────────────────┐
   834 rows: atoms_pickle, element, mp_id, energy_above_hull,        │  DFT reference
   is_ground_state, spacegroup, (+ dft_elastic_* for ~310)          │  (ranking key)
        │                                                            │
        │  seed_track.py  --track sp / --track relax                 │
        ▼                                                            │
  mlip_claim_queue / mlip_relax_claim_queue                          │
   (potential, structure_id) rows, one per structure × potential     │
        │                                                            │
        │  launch_all.sh ab → slurm_elem_{sp,relax}.sh               │
        │  → {worker_v2,relax_worker}.py  (A100 GPU, DELETE-on-claim)│
        ▼                                                            │
  mlip_calc  (Track A)   /   mlip_relax  (Track B)                   │  MLIP prediction
   energy, energy_per_atom, forces, stress, status='done'/'error'    │  (ranking subject)
        │                                                            │
        │  build_rollups.py  (joins track ⨝ structure on element)    │
        ▼              ◄─────────────────────────────────────────────┘
  results/ordering_rollup.csv
   per (potential, element, track): n, spearman, ground_state_correct, dE_mae

1.1 Ingest — benchmark/ingest_mp_elementals.py

• Where the structures come from. MPRester.materials.summary.search(num_elements=1, …) pulls every single-element crystal in Materials Project (→ 834 structures), with fields material_id, formula_pretty, elements, structure, symmetry, energy_above_hull.
• A second call materials.elasticity.search(material_ids=…) pulls DFT elastic refs for the ~310 that have them (this feeds the elastic track, not stability, but is ingested in the same pass).
• Each MP structure → plain ase.Atoms via plain_atoms() (rebuilds positions/cell/numbers/pbc so the pickle depends only on ase+numpy — model envs without pymatgen can unpickle it). Stored as pickled atoms_pickle BYTEA.
• Ground-state flag: is_ground_state = (energy_above_hull < 1e-6).
• Idempotent: dedup on atoms_hash (sha256 of rounded positions/cell/numbers/pbc), ON CONFLICT (atoms_hash) DO NOTHING.
• ⚠️ dft_energy_per_atom is ingested as None (line 84: dft_energy_per_atom=None). MP's summary search does not return a per-atom formation/total energy in this call, so this column is empty for all 834 rows. This is the single most important provenance fact for the ranking — see §3.

1.2 Seed — benchmark/seed_track.py

• For each potential, inserts (potential, structure_id) for all 834 structures into the track's claim queue (mlip_claim_queue for sp, mlip_relax_claim_queue for relax). ON CONFLICT DO NOTHING → re-seedable.
• --only-ground-states / --elements filters exist but were not used for the stability study (full 834 corpus per potential).

1.3 Compute — workers (one row per structure, on A100 GPU)

Track A — benchmark/worker_v2.py (single-point):

1. DELETE … FOR UPDATE SKIP LOCKED RETURNING structure_id — atomically claims a batch off the queue (DELETE-on-claim; the queue only ever holds pending work).
2. Fetch atoms_pickle, attach the calculator (factory dispatch by potential name, _resolve_factory).
3. e = atoms.get_potential_energy(), f = atoms.get_forces(), s = atoms.get_stress(voigt=True). No relaxation — energy at the MP geometry.
4. UPSERT … status='done' into mlip_calc with energy, energy_per_atom = e/len(atoms), forces, stress, runtime, gpu mem, host, slurm job id. Failures → status='error' + error_msg.
5. CUDA-context-poison guard: device-side asserts/illegal-access exit code 75 and re-queue the un-computed remainder so a fresh worker retries them (prevents one bad structure faking a whole batch of errors).

Track B — benchmark/relax_worker.py (relaxation):

• Same claim/fetch/upsert skeleton, but per structure runs an ASE optimizer: opt = BFGS(atoms, logfile=None); opt.run(fmax=0.01, steps=10000).
• Defaults burned into the row: algo=BFGS, fmax_target=0.01, max_steps=10000. Live DB confirms only this one config exists (BFGS / 0.01 / 10000).
• ⚠️ Fixed-cell ion relaxation. The optimizer wraps the bare atoms object — there is no UnitCellFilter / ExpCellFilter / FrechetCellFilter. So only ionic positions relax; the cell is held fixed at the MP cell. (The repo README's "cell+ion" wording is inaccurate for what the code does — the code relaxes ions only. Documented here because it changes how Track B should be interpreted: it is not a full geometry optimization.)
• Writes energy_final, energy_per_atom, converged = (max_force_final ≤ fmax), n_steps, the relaxed atoms_final pickle, forces, stress into mlip_relax.
• Natural key UNIQUE (structure_id, potential, algo, fmax_target) — so the same potential could be relaxed at multiple fmax targets, but for this study there is exactly one.

Energy convention (both tracks): the energy is whatever atoms.get_potential_energy() returns from the MLIP calculator — there is no DFT-style smearing/sigma split here (that energy vs energy_zero distinction applies to the VASP DFT training labels, not to MLIP inference). For the stability ranking only relative energies/atom within one element matter.

1.4 Launch — launch_all.sh ab → slurm_elem_{sp,relax}.sh

• One SLURM job per (potential, track); 4-GPU flux fan-out per node; auto-resubmit chains drain the queue. Tracks A/B reuse the existing model envs via env_map.sh (sp/relax need only model+ASE, no elastic stack). Driver runs in the executorlib launcher env; workers run in the model env.

1.5 Rollup — benchmark/eval/build_rollups.py + ordering_metrics.py

This is where "the study" is actually computed. See §2–§3.

2. The metric definitions — benchmark/eval/ordering_metrics.py

per_element_ordering(dft, mlip) takes two dicts {mp_id: energy_per_atom} for one element and returns:

keys = [k for k in dft if k in mlip and both not None]         # shared polymorphs
if len(keys) < 2:                                              # need ≥2 to rank
    return {n, spearman=NaN, ground_state_correct=None, dE_mae=NaN}
d = [dft[k] for k in keys];  m = [mlip[k] for k in keys]
rho      = spearmanr(d, m).correlation                          # rank correlation of the two orderings
gs_correct = (argmin(m) structure == argmin(d) structure)       # does MLIP's GS == DFT's GS?
dE_mae   = mean(| (m - m.min()) - (d - d.min()) |)              # error in energy-above-own-GS

• spearman ∈ [−1, 1] — how well the MLIP reproduces the rank order of polymorph energies. Higher = better. NaN when an element has <2 polymorphs (nothing to rank).
• ground_state_correct ∈ {True, False, None} — does the MLIP's lowest-energy polymorph match DFT's lowest-energy polymorph? None when <2 polymorphs.
• dE_mae (eV/atom) — referencing each set to its own ground state, the MAE of the relative energies across polymorphs. Lower = better. Penalizes wrong energy spacing even when the rank order is right.
• n — number of shared polymorphs that actually entered the metric.

3. ⚠️ The DFT ranking key — what we ranked against, and why it's valid

The MLIP side of every metric uses energy_per_atom from mlip_calc / mlip_relax — the actual MLIP-predicted energy. No issue.

The DFT side does NOT use a DFT energy. Because structure.dft_energy_per_atom is unpopulated (§1.1), build_rollups._ordering_for() ranks the DFT side by energy_above_hull instead (line 20):

sql = "SELECT s.element, s.mp_id, s.energy_above_hull, t.potential, t.energy_per_atom
       FROM {table} t JOIN structure s ON s.id=t.structure_id WHERE t.status='done'"
...
dft  = {mp_id: energy_above_hull}     # ← DFT ranking key
mlip = {mp_id: energy_per_atom}       # ← MLIP ranking key
per_element_ordering(dft, mlip)

Why this is correct for ordering (the code comment, expanded): within a single element, energy_above_hull = (total energy per atom) − (the element's hull/ground-state energy per atom). The subtracted term is a single constant for that element. Therefore, for one element:

• Spearman ρ is invariant — rank order is unchanged by subtracting a constant. ✓ identical.
• ground_state_correct is invariant — argmin is unchanged by an additive constant. ✓ identical.
• dE_mae is invariant — it already references each set to its own minimum, so the constant cancels in (d − d.min()). ✓ identical.

So all three stability metrics are mathematically identical whether the DFT side is ranked by the (missing) dft_energy_per_atom or by energy_above_hull. The substitution is exact, not an approximation — but only because every comparison is within one element. energy_above_hull must never be used to compare energies across elements (the per-element constant differs); the rollup never does this (it groups by element first).

One caveat to flag for any re-analysis: energy_above_hull is clamped ≥ 0 and is 0 for every ground state. If an element has two polymorphs both exactly on the hull (E_above_hull = 0), DFT sees a tie and argmin picks the first by dict order — a (rare) tie-break artifact. With true total energies this wouldn't tie. Worth checking if a given element shows a suspicious gs result.

4. The rollup builder — benchmark/eval/build_rollups.py

ordering_rollup(conn) = concat(
    _ordering_for(conn, "mlip_calc",  "sp"),      # Track A
    _ordering_for(conn, "mlip_relax", "relax"))   # Track B

• Each _ordering_for selects status='done' rows only (errors are excluded — they never enter the CSV), joins to structure for element + energy_above_hull, groups by (potential, element), and applies per_element_ordering.
• Tags each row with track ∈ {sp, relax}.
• Run command: python -m benchmark.eval.build_rollups --out results/ (python = python).
• Output: results/ordering_rollup.csv (707 KB). Columns: n, spearman, ground_state_correct, dE_mae, potential, element, track.

Naming hazard for downstream consumers: the CSV's track column means sp vs relax, which is a facet within the stability study, not the site-level track concept (stability / elastic / vacancy / phonon). Consider renaming it eval_mode in any frontend ingest.

Per-potential aggregation (CSV → one leaderboard number)

The CSV is per-element; to get one score per potential, aggregate across that potential's elements:

df = pd.read_csv("results/ordering_rollup.csv")
g  = df[df.track=="sp"].groupby("potential")
summary = g.agg(mean_spearman=("spearman","mean"),
                gs_hit_rate =("ground_state_correct","mean"),   # bool→0/1, NaN/None dropped
                mean_dE_mae =("dE_mae","mean"),
                n_elements  =("element","nunique"))

Sort directions: mean_spearman ↓, gs_hit_rate ↓, mean_dE_mae ↑. Denominator caveat: gs_hit_rate from .mean() drops elements with <2 polymorphs (their value is None/NaN), so its denominator < the displayed n_elements. Flag in any tooltip.

5. Current state of the study (live DB, 2026-06-08)

• Both tracks fully attempted (queues drained to 0); all 60 potentials present. 60 × 834 = 50,040 cells per track; the gaps from 50,040 are status='error' rows (real model failures + the known transient/element-coverage classes — see project memory mlip_elemental_benchmark.md).
• Relax config is uniform: BFGS / fmax 0.01 / max_steps 10000, fixed cell.
• Known data-quality flags that bias the energy side (carry into stability too): UMA/equiformer magnetic-TM blowups, m3gnet stress/EFS fragility, dpa transient CUDA casts. These are energy outliers; for ordering they mostly hurt dE_mae more than rank ρ.

6. How to reproduce / refresh the whole study

PY=python
cd .../mlip-elemental-benchmark

# 0. (one-time) ingest — login node, needs internet + key
source .../.mp_api_key
$PY benchmark/ingest_mp_elementals.py

# 1. seed both tracks for a potential (repeat per potential, or via launch_all)
$PY -m benchmark.seed_track --track sp    --potential <POT>
$PY -m benchmark.seed_track --track relax --potential <POT>

# 2. compute on GPUs (all 60, both tracks)
bash benchmark/launch_all.sh ab

# 3. rebuild the rollup once queues drain
$PY -m benchmark.eval.build_rollups --out results/
#   → results/ordering_rollup.csv

Inputs that fully determine the result: the 834 MP structures + their energy_above_hull (ingest), each potential's calculator factory (_resolve_factory), the relax config (BFGS/0.01/10000, fixed cell), and per_element_ordering. Everything else is bookkeeping.

7. Files (authoritative source, in dependency order)

Related: project memory mlip_elemental_benchmark.md; frontend brief FRONTEND_HANDOFF.md; full DB export results/db_export/{structure,mlip_calc,mlip_relax,…}.parquet.

Verification addendum (run wf_7fbcd118-848, 2026-06-08)

• Superseded 2026-06-15 (see top note): the old ordering_rollup.csv carried sp (5239) + fixed-cell relax (5322) from mlip_relax; the variable-cell relax_cell was defined but never run. The site now carries sp (5239, unchanged) + hydro (5233) + full (5233) — the two relaxed modes rebuilt from the standalone relax-cell campaign (relax_export, tools/regen_ordering_from_relax_export.py), not the in-DB mlip_relax Track B.
• Corpus = every num_elements=1 MP entry (ingest_mp_elementals.py:58-62), no convex-hull cutoff. 834 structures / 89 elements / 90 GS rows: Br, I, O have two on-hull GS each; Hg, Pr have none.
• DFT energy — only energy_above_hull is stored (dft_energy_per_atom is NULL for all 834); the smeared-energy/ISMEAR text is an upstream MP/VASP convention, not a choice this track makes.
• CUDA requeue/restart guard exists only in the Track A single-point worker (worker_v2.py, exit 75 + requeue); the relax worker writes an error row and continues.
• Convergence flag = opt.converged() (primary), max|Fᵢ|≤fmax (fallback) (relax_worker.py:236-239).
• Rollup column is literally track (values sp/relax); the rename to eval_mode happens site-side only. schema_relax.sql:1 header names the wrong DB (gbsegbench, copy-paste; cosmetic).