From 9b8543c797265e9748fec2ccb1f9af3ba72e9c1e Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Tue, 29 Sep 2026 23:46:51 -0700 Subject: [PATCH 1/3] fix(SOF-8067): analyze_convex_hull takes each material's total energy from its job on SCF_KGRID when set, and the results table reads Material ID from the phase-diagram entry: `precision.value` is not a k-grid measure (monoclinic HfO2 carries a Gamma-only energy at precision 2000 and a 4x4x4 one at 768), so selecting by highest precision mixed grids; `find_job_for_material_with_property(..., kgrid=SCF_KGRID)` picks the job whose pw_scf unit ran on that grid, searched over every account material with the same exabyteId because the job may sit on a saved copy ("final structure") rather than the material found, and its total_energy is read with `properties.get_for_job`, while SCF_KGRID = None keeps the highest-precision rule; `get_results_table` indexed `entries_data` by position in `phase_diagram.all_entries`, which pymatgen sorts by reduced composition, so IDs were shifted against formulas (HfO2 showed the ZrO2 id, ZrO2 the O2 id); each entry already carries its `material_id` as `entry_id`. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../workflows/analyze_convex_hull.ipynb | 14 +++++++++++--- .../core/entity/property/analysis.py | 5 ++--- 2 files changed, 13 insertions(+), 6 deletions(-) diff --git a/other/materials_designer/workflows/analyze_convex_hull.ipynb b/other/materials_designer/workflows/analyze_convex_hull.ipynb index d50ebf299..13c941661 100644 --- a/other/materials_designer/workflows/analyze_convex_hull.ipynb +++ b/other/materials_designer/workflows/analyze_convex_hull.ipynb @@ -75,6 +75,8 @@ "# Precision value. None = use the best (highest) available. Specific value = use exactly that.\n", "PRECISION = None\n", "\n", + "SCF_KGRID = None # e.g. [4, 4, 4]; None = highest precision\n", + "\n", "# Organization to select.\n", "ORGANIZATION_NAME = None\n" ] @@ -223,7 +225,7 @@ "source": [ "## 4. Retrieve total energies\n", "\n", - "For each material, find completed jobs and extract total energy properties." + "For each material, find completed jobs and extract total energy properties. With `SCF_KGRID` set, the energy comes from the material's finished job whose `pw_scf` unit ran on that k-grid; otherwise from the highest-precision property (or exactly `PRECISION`)." ] }, { @@ -234,6 +236,7 @@ "outputs": [], "source": [ "from collections import Counter\n", + "from mat3ra.notebooks_utils.core.entity.job.api import find_job_for_material_with_property\n", "\n", "exabyte_ids = list(set(material[\"exabyteId\"] for material in all_materials if material.get(\"exabyteId\")))\n", "\n", @@ -267,6 +270,11 @@ " seen_exabyte_ids.add(exabyte_id)\n", "\n", " property_holder = properties_by_exabyte_id.get(exabyte_id)\n", + " if SCF_KGRID:\n", + " materials_with_exabyte_id = client.materials.list({\"exabyteId\": exabyte_id, \"owner._id\": ACCOUNT_ID})\n", + " jobs = (find_job_for_material_with_property(client, candidate[\"_id\"], \"total_energy\", ACCOUNT_ID, kgrid=SCF_KGRID) for candidate in materials_with_exabyte_id)\n", + " job = next(filter(None, jobs), None)\n", + " property_holder = {\"data\": client.properties.get_for_job(job[\"_id\"], property_name=\"total_energy\")[0]} if job else None\n", " if not property_holder:\n", " print(f\"\\u26a0\\ufe0f {material['formula']}: no total energy\")\n", " continue\n", @@ -275,10 +283,10 @@ " number_of_atoms = len(material[\"basis\"][\"coordinates\"])\n", " composition = dict(Counter(element[\"value\"] for element in elements))\n", " energy = property_holder[\"data\"][\"value\"]\n", - " precision = property_holder.get(\"precision\", {}).get(\"value\", \"?\")\n", + " source = f\"job {job['_id']}\" if SCF_KGRID else f\"precision={property_holder.get('precision', {}).get('value', '?')}\"\n", "\n", " entries_data.append({\"material_id\": material[\"_id\"], \"formula\": material[\"formula\"], \"composition\": composition, \"n_atoms\": number_of_atoms, \"total_energy\": energy})\n", - " print(f\"\\u2705 {material['formula']} (precision={precision}): {energy:.4f} eV ({energy / number_of_atoms:.4f} eV/atom)\")\n", + " print(f\"\\u2705 {material['formula']} ({source}): {energy:.4f} eV ({energy / number_of_atoms:.4f} eV/atom)\")\n", "\n", "print(f\"\\n{len(entries_data)} entries.\")\n" ] diff --git a/src/py/mat3ra/notebooks_utils/core/entity/property/analysis.py b/src/py/mat3ra/notebooks_utils/core/entity/property/analysis.py index 7279e3f22..9cffb6b4d 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/property/analysis.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/property/analysis.py @@ -54,15 +54,14 @@ def get_results_table(phase_diagram: PhaseDiagram, entries_data: List[PhaseStabi DataFrame with formula, material ID, energies, stability, and decomposition. """ results = [] - for i, entry in enumerate(phase_diagram.all_entries): + for entry in phase_diagram.all_entries: energy_above_hull = phase_diagram.get_e_above_hull(entry) decomposition = phase_diagram.get_decomposition(entry.composition) decomposition_str = " + ".join([e.composition.reduced_formula for e in decomposition]) - data = entries_data[i] results.append( { "Formula": entry.composition.reduced_formula, - "Material ID": data.get("material_id", ""), + "Material ID": entry.entry_id, "E/atom (eV)": round(entry.energy_per_atom, 4), "Eform/atom (eV)": round(phase_diagram.get_form_energy_per_atom(entry), 4), "Above hull (eV)": round(energy_above_hull, 4), From 8f26dbc6a71e044340695edeb08609d99880a3bf Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Wed, 7 Oct 2026 14:22:41 -0700 Subject: [PATCH 2/3] feat(SOF-8067): analyze_convex_hull prints the chemical potentials the defect formation energy notebooks take MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `get_chemical_potentials_table` lists, for every stable compound, the corners of its stability region as pymatgen's `get_all_chempots` returns them (named by the phases in equilibrium there, e.g. "ZrO2-O2-HfO2"), with each element's chemical potential minus its elemental reference energy per atom: Δμ ≤ 0 in eV, the flat dict `CHEMICAL_POTENTIALS` takes in the defect formation energy notebooks. Until now these were worked out by hand from the formation energies. The notebook shows the table as section 6.2 after the results table. The test checks the five corners of the Hf-Zr-O fixture against values computed by hand, and that Σ nᵢΔμᵢ equals the compound's formation energy per formula unit on every row. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../workflows/analyze_convex_hull.ipynb | 26 +++++++++++++++-- .../core/entity/property/analysis.py | 25 ++++++++++++++++ .../core/entity/test_property_analysis.py | 29 +++++++++++++++++++ 3 files changed, 78 insertions(+), 2 deletions(-) diff --git a/other/materials_designer/workflows/analyze_convex_hull.ipynb b/other/materials_designer/workflows/analyze_convex_hull.ipynb index 13c941661..6c9c2638a 100644 --- a/other/materials_designer/workflows/analyze_convex_hull.ipynb +++ b/other/materials_designer/workflows/analyze_convex_hull.ipynb @@ -343,7 +343,9 @@ "id": "21", "metadata": {}, "source": [ - "## 7. Plot" + "### 6.2. Chemical potentials\n", + "\n", + "Δμ per element at each corner of every stable compound's stability region, relative to the elemental references — the values `CHEMICAL_POTENTIALS` takes in the defect formation energy notebooks." ] }, { @@ -352,6 +354,26 @@ "id": "22", "metadata": {}, "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.core.entity.property.analysis import get_chemical_potentials_table\n", + "\n", + "get_chemical_potentials_table(phase_diagram)\n" + ] + }, + { + "cell_type": "markdown", + "id": "23", + "metadata": {}, + "source": [ + "## 7. Plot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24", + "metadata": {}, + "outputs": [], "source": [ "from mat3ra.notebooks_utils.ipython.plot._plotly import render_figure\n", "\n", @@ -362,7 +384,7 @@ { "cell_type": "code", "execution_count": null, - "id": "23", + "id": "25", "metadata": {}, "outputs": [], "source": [] diff --git a/src/py/mat3ra/notebooks_utils/core/entity/property/analysis.py b/src/py/mat3ra/notebooks_utils/core/entity/property/analysis.py index 9cffb6b4d..a90cd3225 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/property/analysis.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/property/analysis.py @@ -71,3 +71,28 @@ def get_results_table(phase_diagram: PhaseDiagram, entries_data: List[PhaseStabi ) return pd.DataFrame(results).sort_values("Above hull (eV)") + + +def get_chemical_potentials_table(phase_diagram: PhaseDiagram) -> pd.DataFrame: + """Build a chemical potentials DataFrame from phase diagram analysis. + + Args: + phase_diagram: pymatgen PhaseDiagram object. + + Returns: + DataFrame with Δμ per element (eV, relative to the elemental reference) + at each corner of every stable compound's stability region. + """ + results = [] + for entry in phase_diagram.stable_entries: + if entry.is_element: + continue + for facet_name, chemical_potentials in phase_diagram.get_all_chempots(entry.composition).items(): + row = {"Formula": entry.composition.reduced_formula, "Phases in equilibrium": facet_name} + for element, chemical_potential in chemical_potentials.items(): + row[f"Δμ_{element} (eV)"] = round( + chemical_potential - phase_diagram.el_refs[element].energy_per_atom, 4 + ) + results.append(row) + + return pd.DataFrame(results) diff --git a/tests/py/unit/core/entity/test_property_analysis.py b/tests/py/unit/core/entity/test_property_analysis.py index 35eaeba8f..50df3db2c 100644 --- a/tests/py/unit/core/entity/test_property_analysis.py +++ b/tests/py/unit/core/entity/test_property_analysis.py @@ -2,11 +2,14 @@ from typing import List +import pytest from mat3ra.notebooks_utils.core.entity.property.analysis import ( PhaseStabilityEntry, build_convex_hull, + get_chemical_potentials_table, get_results_table, ) +from pymatgen.core import Composition # Realistic total energies (eV) for the Hf-Zr-O system (PBE-level) ENTRIES_DATA: List[PhaseStabilityEntry] = [ @@ -29,6 +32,19 @@ }, ] +# eV per formula unit from ENTRIES_DATA against Hf -11.93, Zr -9.28, O -4.93 eV/atom. +FORMATION_ENERGY_PER_FORMULA_UNIT = {"HfO2": -11.12, "ZrO2": -9.045} +METAL_RICH = {"Hf": 0.0, "Zr": 0.0, "O": -5.56} +ZIRCONIUM_RICH = {"Hf": -2.075, "Zr": 0.0, "O": -4.5225} +OXYGEN_RICH = {"Hf": -11.12, "Zr": -9.045, "O": 0.0} +STABILITY_REGION_CORNERS = [ + ("HfO2", "Hf-Zr-HfO2", METAL_RICH), + ("HfO2", "ZrO2-Zr-HfO2", ZIRCONIUM_RICH), + ("HfO2", "ZrO2-O2-HfO2", OXYGEN_RICH), + ("ZrO2", "ZrO2-Zr-HfO2", ZIRCONIUM_RICH), + ("ZrO2", "ZrO2-O2-HfO2", OXYGEN_RICH), +] + def test_build_convex_hull(): phase_diagram = build_convex_hull(ENTRIES_DATA) @@ -48,6 +64,19 @@ def test_get_results_table(): assert "Eform/atom (eV)" in df.columns +@pytest.mark.parametrize("formula, phases_in_equilibrium, expected_delta_mu", STABILITY_REGION_CORNERS) +def test_get_chemical_potentials_table(formula, phases_in_equilibrium, expected_delta_mu): + df = get_chemical_potentials_table(build_convex_hull(ENTRIES_DATA)) + assert len(df) == len(STABILITY_REGION_CORNERS) + row = df[(df["Formula"] == formula) & (df["Phases in equilibrium"] == phases_in_equilibrium)].iloc[0] + delta_mu = {element: row[f"Δμ_{element} (eV)"] for element in expected_delta_mu} + assert delta_mu == pytest.approx(expected_delta_mu, abs=1e-4) + amounts = Composition(formula).get_el_amt_dict() + assert sum(amounts[element] * delta_mu[element] for element in amounts) == pytest.approx( + FORMATION_ENERGY_PER_FORMULA_UNIT[formula], abs=1e-3 + ) + + def test_elemental_entries_on_hull(): phase_diagram = build_convex_hull(ENTRIES_DATA) df = get_results_table(phase_diagram, ENTRIES_DATA) From 3458082ca72308964212d8c195732ff89d7bca21 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Wed, 7 Oct 2026 14:32:35 -0700 Subject: [PATCH 3/3] fix(SOF-8067): chemical potentials table in a fixed order; test drops the identity check that held by construction MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `stable_entries` is a set, so the compounds came out in a different order per session; they are now sorted by formula. The Σ nᵢΔμᵢ = ΔH_f assertion added in 8f26dbc6 could not fail on its own — the expected Δμ satisfy it by construction — so the test keeps only the comparison with the hand-computed corners. The notebook says the Δμ apply in the defect notebook only when its elemental references are the same calculations as the hull's. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../workflows/analyze_convex_hull.ipynb | 2 +- .../core/entity/property/analysis.py | 2 +- .../core/entity/test_property_analysis.py | 24 +++++++------------ 3 files changed, 11 insertions(+), 17 deletions(-) diff --git a/other/materials_designer/workflows/analyze_convex_hull.ipynb b/other/materials_designer/workflows/analyze_convex_hull.ipynb index 6c9c2638a..58f4c8cd4 100644 --- a/other/materials_designer/workflows/analyze_convex_hull.ipynb +++ b/other/materials_designer/workflows/analyze_convex_hull.ipynb @@ -345,7 +345,7 @@ "source": [ "### 6.2. Chemical potentials\n", "\n", - "Δμ per element at each corner of every stable compound's stability region, relative to the elemental references — the values `CHEMICAL_POTENTIALS` takes in the defect formation energy notebooks." + "Δμ per element at each corner of every stable compound's stability region, relative to the elemental references — the values `CHEMICAL_POTENTIALS` takes in the defect formation energy notebooks. They apply there only if its elemental references are the same calculations: compare their energies per atom with E/atom in the results table." ] }, { diff --git a/src/py/mat3ra/notebooks_utils/core/entity/property/analysis.py b/src/py/mat3ra/notebooks_utils/core/entity/property/analysis.py index a90cd3225..28bf4ae64 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/property/analysis.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/property/analysis.py @@ -84,7 +84,7 @@ def get_chemical_potentials_table(phase_diagram: PhaseDiagram) -> pd.DataFrame: at each corner of every stable compound's stability region. """ results = [] - for entry in phase_diagram.stable_entries: + for entry in sorted(phase_diagram.stable_entries, key=lambda entry: entry.composition.reduced_formula): if entry.is_element: continue for facet_name, chemical_potentials in phase_diagram.get_all_chempots(entry.composition).items(): diff --git a/tests/py/unit/core/entity/test_property_analysis.py b/tests/py/unit/core/entity/test_property_analysis.py index 50df3db2c..286d76cd9 100644 --- a/tests/py/unit/core/entity/test_property_analysis.py +++ b/tests/py/unit/core/entity/test_property_analysis.py @@ -9,7 +9,6 @@ get_chemical_potentials_table, get_results_table, ) -from pymatgen.core import Composition # Realistic total energies (eV) for the Hf-Zr-O system (PBE-level) ENTRIES_DATA: List[PhaseStabilityEntry] = [ @@ -32,17 +31,16 @@ }, ] -# eV per formula unit from ENTRIES_DATA against Hf -11.93, Zr -9.28, O -4.93 eV/atom. -FORMATION_ENERGY_PER_FORMULA_UNIT = {"HfO2": -11.12, "ZrO2": -9.045} -METAL_RICH = {"Hf": 0.0, "Zr": 0.0, "O": -5.56} -ZIRCONIUM_RICH = {"Hf": -2.075, "Zr": 0.0, "O": -4.5225} -OXYGEN_RICH = {"Hf": -11.12, "Zr": -9.045, "O": 0.0} +# From ENTRIES_DATA: Hf -11.93, Zr -9.28, O -4.93 eV/atom; formation energy HfO2 -11.12, ZrO2 -9.045 eV/formula unit. +DELTA_MU_METAL_RICH = {"Hf": 0.0, "Zr": 0.0, "O": -5.56} +DELTA_MU_ZIRCONIUM_RICH = {"Hf": -2.075, "Zr": 0.0, "O": -4.5225} +DELTA_MU_OXYGEN_RICH = {"Hf": -11.12, "Zr": -9.045, "O": 0.0} STABILITY_REGION_CORNERS = [ - ("HfO2", "Hf-Zr-HfO2", METAL_RICH), - ("HfO2", "ZrO2-Zr-HfO2", ZIRCONIUM_RICH), - ("HfO2", "ZrO2-O2-HfO2", OXYGEN_RICH), - ("ZrO2", "ZrO2-Zr-HfO2", ZIRCONIUM_RICH), - ("ZrO2", "ZrO2-O2-HfO2", OXYGEN_RICH), + ("HfO2", "Hf-Zr-HfO2", DELTA_MU_METAL_RICH), + ("HfO2", "ZrO2-Zr-HfO2", DELTA_MU_ZIRCONIUM_RICH), + ("HfO2", "ZrO2-O2-HfO2", DELTA_MU_OXYGEN_RICH), + ("ZrO2", "ZrO2-Zr-HfO2", DELTA_MU_ZIRCONIUM_RICH), + ("ZrO2", "ZrO2-O2-HfO2", DELTA_MU_OXYGEN_RICH), ] @@ -71,10 +69,6 @@ def test_get_chemical_potentials_table(formula, phases_in_equilibrium, expected_ row = df[(df["Formula"] == formula) & (df["Phases in equilibrium"] == phases_in_equilibrium)].iloc[0] delta_mu = {element: row[f"Δμ_{element} (eV)"] for element in expected_delta_mu} assert delta_mu == pytest.approx(expected_delta_mu, abs=1e-4) - amounts = Composition(formula).get_el_amt_dict() - assert sum(amounts[element] * delta_mu[element] for element in amounts) == pytest.approx( - FORMATION_ENERGY_PER_FORMULA_UNIT[formula], abs=1e-3 - ) def test_elemental_entries_on_hull():