diff --git a/other/materials_designer/workflows/analyze_convex_hull.ipynb b/other/materials_designer/workflows/analyze_convex_hull.ipynb index d50ebf299..58f4c8cd4 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" ] @@ -335,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. They apply there only if its elemental references are the same calculations: compare their energies per atom with E/atom in the results table." ] }, { @@ -344,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", @@ -354,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 7279e3f22..28bf4ae64 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), @@ -72,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 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(): + 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..286d76cd9 100644 --- a/tests/py/unit/core/entity/test_property_analysis.py +++ b/tests/py/unit/core/entity/test_property_analysis.py @@ -2,9 +2,11 @@ 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, ) @@ -29,6 +31,18 @@ }, ] +# 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", 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), +] + def test_build_convex_hull(): phase_diagram = build_convex_hull(ENTRIES_DATA) @@ -48,6 +62,15 @@ 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) + + def test_elemental_entries_on_hull(): phase_diagram = build_convex_hull(ENTRIES_DATA) df = get_results_table(phase_diagram, ENTRIES_DATA)