From dab746531489e7a3ee24f5366052d42f5262db8d Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 28 Sep 2026 15:22:16 -0700 Subject: [PATCH 01/18] feat(SOF-7975): optionally relax each defective supercell in its charge state before its SCF With RELAX_DEFECTIVE_MATERIAL, the charged defect notebook runs the Standata Fixed-cell Relaxation on every defective supercell for every charge, with tot_charge patched into pw_relax and the supercell's own k-grid, reuses a finished relaxation of the same cell and charge by its workflow name, submits the relaxations one at a time, and runs each Defect Formation Energy job on its relaxed structure. The relaxations stay untagged, so their total energy can never be taken for a pristine reference. A live run on the unrelaxed cells left residual forces of 0.7 eV/A at q=0 and 8.5 eV/A at q=+2, which makes the charged energies meaningless. The header now also states that no potential alignment is applied. The flag defaults to False, so the notebook behaves as before unless it is set. Co-Authored-By: Claude Fable 5.1 --- .../defect_formation_energy_charged.ipynb | 75 +++++++++++++++---- 1 file changed, 59 insertions(+), 16 deletions(-) diff --git a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb index ff38759ef..236c2c934 100644 --- a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb @@ -24,9 +24,9 @@ "\n", "$E_{\\text{tot}}[\\text{bulk}]$, $E_{\\text{VBM}}$ and $E_{\\text{gap}}$ come from a Total Energy and a Band Gap job on the pristine supercell; $\\mu_i$ from the Standata elemental reference materials, as in [Formation Energy](formation_energy.ipynb). The charge $q$ enters as `tot_charge` in the QE `&SYSTEM` namelist and is compensated by a uniform jellium background.\n", "\n", - "Every Total Energy, Band Gap and defect job is tagged with its charge, `charge:`. The pristine references run neutral, tagged `charge:0`, and the notebook names those jobs to the workflow, which reads its references from them alone, so a charged calculation on the same pristine cell is never taken for the neutral reference. The pristine cell can optionally be relaxed first (variable cell), and every supercell is then built from the relaxed cell.\n", + "Every Total Energy, Band Gap and defect job is tagged with its charge, `charge:`. The pristine references run neutral, tagged `charge:0`, and the notebook names those jobs to the workflow, which reads its references from them alone, so a charged calculation on the same pristine cell is never taken for the neutral reference. The pristine cell can optionally be relaxed first (variable cell), and every supercell is then built from the relaxed cell. With `RELAX_DEFECTIVE_MATERIAL`, each defective supercell is also relaxed at fixed cell in its charge state before its SCF.\n", "\n", - "No analytical finite-size correction (Freysoldt-Neugebauer-Van de Walle, Makov-Payne) is applied. Running several `SUPERCELL_SCALINGS` and extrapolating $E_f(L\\to\\infty)$ from the fitted image-charge terms takes its place; a single size gives the uncorrected value for that cell.\n", + "No analytical finite-size correction (Freysoldt-Neugebauer-Van de Walle, Makov-Payne) is applied. Running several `SUPERCELL_SCALINGS` and extrapolating $E_f(L\\to\\infty)$ from the fitted image-charge terms takes its place; a single size gives the uncorrected value for that cell. No potential alignment ($q\\,\\Delta V$) is applied either.\n", "\n", "

Usage

\n", "\n", @@ -133,7 +133,8 @@ "\n", "SCF_KGRID = None # e.g. [8, 8, 8]\n", "\n", - "RELAX_PRISTINE_MATERIAL = False" + "RELAX_PRISTINE_MATERIAL = False\n", + "RELAX_DEFECTIVE_MATERIAL = False" ] }, { @@ -647,7 +648,8 @@ "id": "41", "metadata": {}, "source": [ - "### 6.2. Create the Defect Formation Energy jobs, one per size and charge" + "### 6.2. Relax the defective supercells (optional)\n", + "With `RELAX_DEFECTIVE_MATERIAL`, the atomic positions of each defective supercell are relaxed in each charge state with the Fixed-cell Relaxation workflow, one relaxation after another, or taken from an earlier finished relaxation of the same cell and charge, and the Defect Formation Energy jobs below run on the relaxed structures. The lattice stays that of the pristine supercell." ] }, { @@ -656,6 +658,46 @@ "id": "42", "metadata": {}, "outputs": [], + "source": [ + "relaxed_defective_supercells = {}\n", + "if RELAX_DEFECTIVE_MATERIAL:\n", + " defect_relax_workflow_config = WorkflowStandata.filter_by_application(APPLICATION_NAME).get_by_name_first_match(\n", + " \"fixed_cell_relaxation.json\"\n", + " )\n", + " for scaling, (_, defective_supercell) in saved_supercell_pairs.items():\n", + " for charge in CHARGES:\n", + " workflow = Workflow.create(defect_relax_workflow_config)\n", + " workflow.name = f\"{workflow.name} {defective_supercell.name} q={charge:+d}\"\n", + " apply_scf_kgrid(\n", + " workflow, get_scf_kgrid_for_supercell(scaling), material=defective_supercell, unit_name=\"pw_relax\"\n", + " )\n", + " if charge:\n", + " patch_workflow_qe_input(workflow, {\"system\": {\"tot_charge\": charge}}, unit_names=[\"pw_relax\"])\n", + " job = find_job_for_material(client, defective_supercell.id, workflow.name, ACCOUNT_ID)\n", + " if job is not None:\n", + " print(f\"♻️ n={scaling} q={charge:+d}: reusing the relaxation from job {job['_id']}\")\n", + " else:\n", + " job = create_job_for_materials([defective_supercell], workflow)\n", + " print(f\"✅ n={scaling} q={charge:+d}: created a relaxation job {job['_id']}\")\n", + " submit_jobs(client.jobs, [job[\"_id\"]])\n", + " await wait_for_jobs_to_finish_async(client.jobs, [job[\"_id\"]], poll_interval=POLL_INTERVAL)\n", + " relaxed_defective_supercells[(scaling, charge)] = get_final_structure_for_job(client, job[\"_id\"])" + ] + }, + { + "cell_type": "markdown", + "id": "43", + "metadata": {}, + "source": [ + "### 6.3. Create the Defect Formation Energy jobs, one per size and charge" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "44", + "metadata": {}, + "outputs": [], "source": [ "import pandas as pd\n", "\n", @@ -663,8 +705,9 @@ "for scaling, (pristine_supercell, defective_supercell) in saved_supercell_pairs.items():\n", " for charge in CHARGES:\n", " workflow = create_defect_workflow(scaling, charge, pristine_reference_jobs[scaling])\n", + " defective_material = relaxed_defective_supercells.get((scaling, charge), defective_supercell)\n", " # Order matters: [0] defective (computed), [1] pristine (reference).\n", - " job = create_job_for_materials([defective_supercell, pristine_supercell], workflow, [f\"charge:{charge}\"])\n", + " job = create_job_for_materials([defective_material, pristine_supercell], workflow, [f\"charge:{charge}\"])\n", " print(f\"n={scaling} q={charge:+d}: k-grid {get_scf_kgrid_for_supercell(scaling) or 'from KPPRA'}, \"\n", " f\"job {job['_id']}\")\n", " job_records.append({\n", @@ -681,16 +724,16 @@ }, { "cell_type": "markdown", - "id": "43", + "id": "45", "metadata": {}, "source": [ - "### 6.3. Submit the jobs and monitor the statuses" + "### 6.4. Submit the jobs and monitor the statuses" ] }, { "cell_type": "code", "execution_count": null, - "id": "44", + "id": "46", "metadata": {}, "outputs": [], "source": [ @@ -702,7 +745,7 @@ { "cell_type": "code", "execution_count": null, - "id": "45", + "id": "47", "metadata": {}, "outputs": [], "source": [ @@ -712,7 +755,7 @@ }, { "cell_type": "markdown", - "id": "46", + "id": "48", "metadata": {}, "source": [ "## 7. Retrieve results\n", @@ -723,7 +766,7 @@ { "cell_type": "code", "execution_count": null, - "id": "47", + "id": "49", "metadata": {}, "outputs": [], "source": [ @@ -747,7 +790,7 @@ }, { "cell_type": "markdown", - "id": "48", + "id": "50", "metadata": {}, "source": [ "### 7.2. Formation energy versus the electron chemical potential\n", @@ -757,7 +800,7 @@ { "cell_type": "code", "execution_count": null, - "id": "49", + "id": "51", "metadata": {}, "outputs": [], "source": [ @@ -781,7 +824,7 @@ { "cell_type": "code", "execution_count": null, - "id": "50", + "id": "52", "metadata": {}, "outputs": [], "source": [ @@ -805,7 +848,7 @@ }, { "cell_type": "markdown", - "id": "51", + "id": "53", "metadata": {}, "source": [ "### 7.3. Finite-size extrapolation\n", @@ -815,7 +858,7 @@ { "cell_type": "code", "execution_count": null, - "id": "52", + "id": "54", "metadata": {}, "outputs": [], "source": [ From 489ba30b37609a26b3878e08a988ab417945f89d Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 28 Sep 2026 17:53:37 -0700 Subject: [PATCH 02/18] fix: stop re-exporting load_material from the top-level material module MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit mat3ra.notebooks_utils.material is imported by the made notebooks, which JupyterLite runs with only the default + made packages from config.yml. Since fdb62d27 (SOF-8044) it re-exported load_material from core/entity/material/api.py, which imports mat3ra.api_client and mat3ra.prode — neither is installed in that environment, so every made notebook that imports set_materials or load_material_from_folder failed at import. Nothing imports load_material through this module (the defect notebooks and test_material_api.py import it from core.entity.material.api), so the two lines go and nothing else changes. Co-Authored-By: Claude Fable 5.1 --- src/py/mat3ra/notebooks_utils/material.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/py/mat3ra/notebooks_utils/material.py b/src/py/mat3ra/notebooks_utils/material.py index b1f04ee1c..69cb41cfd 100644 --- a/src/py/mat3ra/notebooks_utils/material.py +++ b/src/py/mat3ra/notebooks_utils/material.py @@ -1,4 +1,3 @@ -from .core.entity.material.api import load_material from .core.entity.material.io import get_materials, load_material_from_folder, load_materials_from_folder, set_materials __all__ = [ @@ -6,5 +5,4 @@ "set_materials", "load_materials_from_folder", "load_material_from_folder", - "load_material", ] From e4cebc228c7f1fc7fd8e49e6b4da4bc8b4628ee1 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 28 Sep 2026 18:11:10 -0700 Subject: [PATCH 03/18] feat(SOF-7975): reuse a pristine reference only when it ran on the same k-grid A rerun of the charged defect notebook with SCF_KGRID=[1,1,1] silently reused the charge:0 Total Energy and Band Gap jobs computed at 4x4x4, because find_job_for_material_with_property matched only material, owner, status and tag, so the formation energy mixed two grids. The finder now takes an optional kgrid and asks the platform only for jobs whose pw_scf unit carries a kgrid context with those dimensions (new get_kgrid_query, matching workflow.subworkflows[].units[name].context[name="kgrid"].data.dimensions, where apply_scf_kgrid writes it; checked on production jobs rbGKefNLtim9FmGKY and R2xmsS6QEz6jjx6oa). The notebook passes the supercell's grid, prints it on the reuse line, and the section text no longer says a reference is reused whatever grid it ran with. With SCF_KGRID unset the grid is left to the platform's KPPRA default and the lookup is unfiltered as before. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../defect_formation_energy_charged.ipynb | 6 +++--- .../notebooks_utils/core/entity/job/api.py | 19 +++++++++++++++-- tests/py/unit/core/entity/test_job_api.py | 21 +++++++++++++++++++ 3 files changed, 41 insertions(+), 5 deletions(-) diff --git a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb index 236c2c934..a8974a43c 100644 --- a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb @@ -286,7 +286,7 @@ "metadata": {}, "source": [ "### 3.3. Define how jobs are created\n", - "Every job is tagged with its charge; the pristine references run neutral, tagged `charge:0`. A reference is reused when a finished `charge:0` job on that pristine cell already reported the property (`total_energy`, `band_gaps`), whatever k-grid it ran with. The Defect Formation Energy workflow reads each reference from the one job the notebook names to it, so no other calculation on the same cell can stand in; the relaxation carries no charge tag and is never taken for one." + "Every job is tagged with its charge; the pristine references run neutral, tagged `charge:0`. A reference is reused when a finished `charge:0` job on that pristine cell already reported the property (`total_energy`, `band_gaps`), on the same k-grid when `SCF_KGRID` is set. The Defect Formation Energy workflow reads each reference from the one job the notebook names to it, so no other calculation on the same cell can stand in; the relaxation carries no charge tag and is never taken for one." ] }, { @@ -614,10 +614,10 @@ " kgrid = get_scf_kgrid_for_supercell(scaling)\n", " for property_name, workflow_config in pristine_workflow_configs.items():\n", " job = find_job_for_material_with_property(\n", - " client, pristine_supercell.id, property_name, ACCOUNT_ID, tags=[\"charge:0\"]\n", + " client, pristine_supercell.id, property_name, ACCOUNT_ID, tags=[\"charge:0\"], kgrid=kgrid\n", " )\n", " if job is not None:\n", - " print(f\"♻️ n={scaling}: reusing {property_name} from job {job['name']}\")\n", + " print(f\"♻️ n={scaling}: reusing {property_name} from job {job['name']}, k-grid {kgrid or 'from KPPRA'}\")\n", " else:\n", " workflow = Workflow.create(workflow_config)\n", " workflow.name = f\"{workflow.name} {pristine_supercell.name}\"\n", diff --git a/src/py/mat3ra/notebooks_utils/core/entity/job/api.py b/src/py/mat3ra/notebooks_utils/core/entity/job/api.py index f2a6ff443..0df7fa1a0 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/job/api.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/job/api.py @@ -153,16 +153,30 @@ def find_job_for_material( return existing[0] if existing else None +def get_kgrid_query(kgrid: Optional[List[int]], unit_name: str = "pw_scf") -> Dict[str, Any]: + """ + `jobs.list` condition for jobs whose `unit_name` unit ran on `kgrid` (empty when `kgrid` is None), matched where + `apply_scf_kgrid` sets it: `workflow.subworkflows[].units[name].context[name="kgrid"].data.dimensions`. + A job created without an explicit k-grid has no such context and never matches. + """ + if kgrid is None: + return {} + kgrid_context = {"$elemMatch": {"name": "kgrid", "data.dimensions": list(kgrid)}} + return {"workflow.subworkflows.units": {"$elemMatch": {"name": unit_name, "context": kgrid_context}}} + + def find_job_for_material_with_property( api_client: APIClient, material_id: str, property_name: str, owner_id: str, tags: Optional[List[str]] = None, + kgrid: Optional[List[int]] = None, ) -> Optional[dict]: """ Finds a finished job on a material that reported the given property, optionally among the jobs - carrying every one of `tags` (e.g. ["charge:0"] for a reference computed in the neutral state). + carrying every one of `tags` (e.g. ["charge:0"] for a reference computed in the neutral state) + and among those whose `pw_scf` unit ran on `kgrid`. Args: api_client (APIClient): API client instance carrying the authorization context. @@ -170,6 +184,7 @@ def find_job_for_material_with_property( property_name (str): Property the job must have reported, e.g. "total_energy". owner_id (str): Account ID the job must belong to. tags (List[str], optional): Tags the job must all carry. + kgrid (List[int], optional): K-grid dimensions the job's `pw_scf` unit ran on, see `get_kgrid_query`. Returns: dict, optional: The first matching job, or None if none exists. @@ -177,7 +192,7 @@ def find_job_for_material_with_property( query: Dict[str, Any] = {"_material._id": material_id, "owner._id": owner_id, "status": "finished"} if tags: query["tags"] = {"$all": list(tags)} - jobs = api_client.jobs.list(query) + jobs = api_client.jobs.list({**query, **get_kgrid_query(kgrid)}) return next( (job for job in jobs if api_client.properties.get_for_job(job["_id"], property_name=property_name)), None ) diff --git a/tests/py/unit/core/entity/test_job_api.py b/tests/py/unit/core/entity/test_job_api.py index 8d90484b5..5d766d24d 100644 --- a/tests/py/unit/core/entity/test_job_api.py +++ b/tests/py/unit/core/entity/test_job_api.py @@ -175,3 +175,24 @@ def test_find_job_for_material_with_property_returns_none_when_no_job_reported_i assert job is None assert "tags" not in client.jobs.list.call_args.args[0] + + +KGRID_QUERY: Dict[str, Any] = { + "workflow.subworkflows.units": { + "$elemMatch": {"name": "pw_scf", "context": {"$elemMatch": {"name": "kgrid", "data.dimensions": [4, 4, 4]}}} + } +} + + +@pytest.mark.parametrize(("kgrid", "expected_kgrid_query"), [(None, {}), ([4, 4, 4], KGRID_QUERY)]) +def test_find_job_for_material_with_property_matches_the_pw_scf_kgrid(kgrid, expected_kgrid_query): + client = MagicMock() + client.jobs.list.return_value = [EXISTING_JOB] + client.properties.get_for_job.return_value = [{"name": PROPERTY_NAME}] + + job = find_job_for_material_with_property(client, MATERIAL_INITIAL["_id"], PROPERTY_NAME, OWNER_ID, kgrid=kgrid) + + assert job == EXISTING_JOB + client.jobs.list.assert_called_once_with( + {"_material._id": MATERIAL_INITIAL["_id"], "owner._id": OWNER_ID, "status": "finished", **expected_kgrid_query} + ) From 7cc92ea97c3b4ecf4fb0c4656d0d71d6beec899b Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 28 Sep 2026 18:13:25 -0700 Subject: [PATCH 04/18] feat(SOF-7975): find the relaxed pristine cell by its structure and k-grid, not by workflow name The pristine relaxation was reused through find_job_for_material on the exact workflow name "Variable-cell Relaxation ", which breaks on a rename and checks neither the kind of relaxation nor its grid. The notebook now asks find_relaxed_material for a finished relaxation of any material with the same structural hash, and find_relaxed_material takes the same optional kgrid and unit_name as the reference finder, so only a pw_vc-relax unit run on SCF_KGRID counts; a fixed-cell relaxation of the same cell, which exists on the production account for GaN, is not taken for it. The relaxation is created only when none is found and stays untagged, so its total energy can never serve as a pristine reference. The find_job_for_material import moves to the defective-cell relaxation cell, its only user. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../defect_formation_energy_charged.ipynb | 31 ++++++++++++------- .../core/entity/material/api.py | 12 +++++-- .../py/unit/core/entity/test_material_api.py | 26 ++++++++++++++++ 3 files changed, 54 insertions(+), 15 deletions(-) diff --git a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb index a8974a43c..50e21b1b4 100644 --- a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb @@ -352,7 +352,7 @@ "metadata": {}, "source": [ "### 4.2. Relax the pristine cell (optional)\n", - "With `RELAX_PRISTINE_MATERIAL`, the lattice and the atomic positions of the pristine cell are relaxed with the Variable-cell Relaxation workflow, or taken from an earlier finished relaxation of the same cell, and the supercells below are built from the relaxed cell." + "With `RELAX_PRISTINE_MATERIAL`, the lattice and the atomic positions of the pristine cell are relaxed with the Variable-cell Relaxation workflow, or taken from an earlier finished variable-cell relaxation of the same structure, on the same k-grid when `SCF_KGRID` is set, and the supercells below are built from the relaxed cell." ] }, { @@ -362,29 +362,34 @@ "metadata": {}, "outputs": [], "source": [ - "from mat3ra.notebooks_utils.core.entity.job.api import find_job_for_material\n", - "from mat3ra.notebooks_utils.core.entity.material.api import get_final_structure_for_job, get_or_create_material\n", + "from mat3ra.notebooks_utils.core.entity.material.api import (\n", + " find_relaxed_material,\n", + " get_final_structure_for_job,\n", + " get_or_create_material,\n", + ")\n", "from mat3ra.notebooks_utils.workflow import apply_scf_kgrid\n", "from mat3ra.standata.workflows import WorkflowStandata\n", "from mat3ra.wode.workflows import Workflow\n", "\n", "if RELAX_PRISTINE_MATERIAL:\n", " saved_pristine_material = Material.create(get_or_create_material(client, pristine_material, ACCOUNT_ID))\n", - " relax_workflow_config = WorkflowStandata.filter_by_application(APPLICATION_NAME).get_by_name_first_match(\n", - " \"variable_cell_relaxation.json\"\n", + " relaxed_pristine_material = find_relaxed_material(\n", + " client, saved_pristine_material, ACCOUNT_ID, kgrid=SCF_KGRID, unit_name=\"pw_vc-relax\"\n", " )\n", - " relax_workflow = Workflow.create(relax_workflow_config)\n", - " relax_workflow.name = f\"{relax_workflow.name} {saved_pristine_material.name}\"\n", - " apply_scf_kgrid(relax_workflow, SCF_KGRID, material=saved_pristine_material, unit_name=\"pw_vc-relax\")\n", - " relax_job = find_job_for_material(client, saved_pristine_material.id, relax_workflow.name, ACCOUNT_ID)\n", - " if relax_job is None:\n", + " if relaxed_pristine_material is None:\n", + " relax_workflow_config = WorkflowStandata.filter_by_application(APPLICATION_NAME).get_by_name_first_match(\n", + " \"variable_cell_relaxation.json\"\n", + " )\n", + " relax_workflow = Workflow.create(relax_workflow_config)\n", + " relax_workflow.name = f\"{relax_workflow.name} {saved_pristine_material.name}\"\n", + " apply_scf_kgrid(relax_workflow, SCF_KGRID, material=saved_pristine_material, unit_name=\"pw_vc-relax\")\n", " relax_job = create_job_for_materials([saved_pristine_material], relax_workflow)\n", " submit_jobs(client.jobs, [relax_job[\"_id\"]])\n", " await wait_for_jobs_to_finish_async(client.jobs, [relax_job[\"_id\"]], poll_interval=POLL_INTERVAL)\n", - " relaxed_pristine_material = get_final_structure_for_job(client, relax_job[\"_id\"])\n", + " relaxed_pristine_material = get_final_structure_for_job(client, relax_job[\"_id\"])\n", " relaxed_pristine_material.name = f\"{pristine_material.name} relaxed\"\n", " print(f\"Lattice constant a: {pristine_material.lattice.a:.4f} Å as given, \"\n", - " f\"{relaxed_pristine_material.lattice.a:.4f} Å relaxed (job {relax_job['_id']})\")\n", + " f\"{relaxed_pristine_material.lattice.a:.4f} Å relaxed\")\n", " pristine_material = relaxed_pristine_material" ] }, @@ -659,6 +664,8 @@ "metadata": {}, "outputs": [], "source": [ + "from mat3ra.notebooks_utils.core.entity.job.api import find_job_for_material\n", + "\n", "relaxed_defective_supercells = {}\n", "if RELAX_DEFECTIVE_MATERIAL:\n", " defect_relax_workflow_config = WorkflowStandata.filter_by_application(APPLICATION_NAME).get_by_name_first_match(\n", diff --git a/src/py/mat3ra/notebooks_utils/core/entity/material/api.py b/src/py/mat3ra/notebooks_utils/core/entity/material/api.py index 8cd3b3e03..d32c82c8e 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/material/api.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/material/api.py @@ -6,6 +6,7 @@ from mat3ra.made.material import Material from mat3ra.prode import PropertyName +from ..job.api import get_kgrid_query from ..property.api import get_properties_for_job from .analysis import get_slab_bulk_crystal, resolve_bulk_query_from_crystal from .io import load_materials_from_folder @@ -97,22 +98,27 @@ def get_final_structure_for_job(api_client: APIClient, job_id: str) -> Material: return Material.create(api_client.materials.get(properties[-1]["materialId"])) -def find_relaxed_material(api_client: APIClient, material, owner_id: str) -> Optional[Material]: +def find_relaxed_material( + api_client: APIClient, material, owner_id: str, kgrid: Optional[List[int]] = None, unit_name: str = "pw_scf" +) -> Optional[Material]: """ Finds a relaxed version of a material: the final structure of a finished job on a material - with the same structural hash, where the geometry has changed. + with the same structural hash, where the geometry has changed, optionally among the jobs whose + `unit_name` unit ran on `kgrid`. Args: api_client (APIClient): API client instance carrying the authorization context. material: mat3ra-made Material object (must have a .hash property). owner_id (str): Account ID under which to search. + kgrid (List[int], optional): K-grid dimensions the relaxation ran on, see `get_kgrid_query`. + unit_name (str): Name of the relaxation unit, e.g. "pw_vc-relax". Returns: Material, optional: The relaxed structure, or None if none exists. """ ids = [m["_id"] for m in api_client.materials.list({"hash": material.hash, "owner._id": owner_id})] query = {"_material._id": {"$in": ids}, "owner._id": owner_id, "status": "finished"} - for job in api_client.jobs.list(query): + for job in api_client.jobs.list({**query, **get_kgrid_query(kgrid, unit_name)}): properties = api_client.properties.get_for_job(job["_id"], PropertyName.non_scalar.final_structure.value) if not properties: continue diff --git a/tests/py/unit/core/entity/test_material_api.py b/tests/py/unit/core/entity/test_material_api.py index 2985ef887..c6c39d558 100644 --- a/tests/py/unit/core/entity/test_material_api.py +++ b/tests/py/unit/core/entity/test_material_api.py @@ -359,6 +359,32 @@ def test_find_relaxed_material_skips_a_final_structure_with_the_same_hash(): assert client.properties.get_for_job.call_count == 2 +def test_find_relaxed_material_matches_the_relaxation_kgrid(): + client = MagicMock() + client.materials.list.return_value = [SAVED_DEFECTIVE] + client.jobs.list.return_value = [FINISHED_JOB] + client.properties.get_for_job.return_value = [{"materialId": "m-relaxed"}] + client.materials.get.return_value = RELAXED_MATERIAL_DOC + + relaxed = find_relaxed_material(client, DEFECTIVE_MATERIAL, OWNER_ID, kgrid=[4, 4, 4], unit_name="pw_vc-relax") + + assert relaxed is not None + assert relaxed.name == "B-vacancy h-BN relaxed" + client.jobs.list.assert_called_once_with( + { + "_material._id": {"$in": [SAVED_DEFECTIVE["_id"]]}, + "owner._id": OWNER_ID, + "status": "finished", + "workflow.subworkflows.units": { + "$elemMatch": { + "name": "pw_vc-relax", + "context": {"$elemMatch": {"name": "kgrid", "data.dimensions": [4, 4, 4]}}, + } + }, + } + ) + + @pytest.mark.parametrize( ("folder_names", "account_names", "name", "expected_name"), [ From 5ca77d51b7699ea43cf79d0df41cd6a1a7b4b39c Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 28 Sep 2026 18:13:58 -0700 Subject: [PATCH 05/18] feat(SOF-7975): set the job time limit from a TIME_LIMIT parameter, 4 hours by default Compute defaults timeLimit to "01:00:00", and a live defective-cell relaxation of V_O in HfO2 was killed by it before it could close. The compute parameters now carry TIME_LIMIT = "04:00:00" next to the queue and ppn, and the notebook passes it to Compute, so every relaxation, reference and defect job gets it and the reader can raise it for larger supercells. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../workflows/defect_formation_energy_charged.ipynb | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb index 50e21b1b4..668b1397f 100644 --- a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb @@ -102,6 +102,7 @@ "CLUSTER_NAME = None # specify full or partial name i.e. \"cluster-001\" to select\n", "QUEUE_NAME = QueueName.D\n", "PPN = 1\n", + "TIME_LIMIT = \"04:00:00\"\n", "\n", "# 6. Job parameters\n", "timestamp = datetime.now().strftime(\"%Y-%m-%d %H:%M\")\n", @@ -276,7 +277,7 @@ "else:\n", " cluster = clusters[0]\n", "\n", - "compute = Compute(cluster=cluster, queue=QUEUE_NAME, ppn=PPN)\n", + "compute = Compute(cluster=cluster, queue=QUEUE_NAME, ppn=PPN, timeLimit=TIME_LIMIT)\n", "print(f\"Using cluster: {compute.cluster.hostname}, queue: {QUEUE_NAME}, ppn: {PPN}\")" ] }, From 50a5358e14add2a5e10505c7247de4aa4c3f6b4b Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 28 Sep 2026 18:18:30 -0700 Subject: [PATCH 06/18] feat(SOF-7975): relax and submit per charge state, and wait on a running relaxation instead of duplicating it With RELAX_DEFECTIVE_MATERIAL the notebook relaxed every charge state before creating any Defect Formation Energy job, so the first result waited on all relaxations, and a rerun while a relaxation was still running created a second one because reuse accepted only finished jobs. The relaxation and the defect job now share one loop over (size, charge): the relaxation is reused or created, waited on, and that charge's defect job is created and submitted at once, so q=+1's SCF runs while q=0 relaxes; the separate batch submission cell goes, which is fewer lines than keeping it. Reuse accepts finished, active, submitted and queued relaxations (not pre-submission, which a never-submitted job would keep forever) and, through the new kgrid and unit_name arguments of find_job_for_material, only those whose pw_relax unit ran on the supercell's grid. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../defect_formation_energy_charged.ipynb | 104 +++++++----------- .../notebooks_utils/core/entity/job/api.py | 8 +- tests/py/unit/core/entity/test_job_api.py | 12 ++ 3 files changed, 56 insertions(+), 68 deletions(-) diff --git a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb index 668b1397f..4c910e751 100644 --- a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb @@ -654,8 +654,8 @@ "id": "41", "metadata": {}, "source": [ - "### 6.2. Relax the defective supercells (optional)\n", - "With `RELAX_DEFECTIVE_MATERIAL`, the atomic positions of each defective supercell are relaxed in each charge state with the Fixed-cell Relaxation workflow, one relaxation after another, or taken from an earlier finished relaxation of the same cell and charge, and the Defect Formation Energy jobs below run on the relaxed structures. The lattice stays that of the pristine supercell." + "### 6.2. Create and submit the Defect Formation Energy jobs, one per size and charge\n", + "With `RELAX_DEFECTIVE_MATERIAL`, the atomic positions of each defective supercell are first relaxed in that charge state with the Fixed-cell Relaxation workflow, and the Defect Formation Energy job of that size and charge is submitted on the relaxed structure as soon as its relaxation has finished, before the next charge state is relaxed. An earlier relaxation of the same cell and charge (on the same k-grid when `SCF_KGRID` is set) is reused, and waited on while it is still running. The lattice stays that of the pristine supercell." ] }, { @@ -664,58 +664,40 @@ "id": "42", "metadata": {}, "outputs": [], - "source": [ - "from mat3ra.notebooks_utils.core.entity.job.api import find_job_for_material\n", - "\n", - "relaxed_defective_supercells = {}\n", - "if RELAX_DEFECTIVE_MATERIAL:\n", - " defect_relax_workflow_config = WorkflowStandata.filter_by_application(APPLICATION_NAME).get_by_name_first_match(\n", - " \"fixed_cell_relaxation.json\"\n", - " )\n", - " for scaling, (_, defective_supercell) in saved_supercell_pairs.items():\n", - " for charge in CHARGES:\n", - " workflow = Workflow.create(defect_relax_workflow_config)\n", - " workflow.name = f\"{workflow.name} {defective_supercell.name} q={charge:+d}\"\n", - " apply_scf_kgrid(\n", - " workflow, get_scf_kgrid_for_supercell(scaling), material=defective_supercell, unit_name=\"pw_relax\"\n", - " )\n", - " if charge:\n", - " patch_workflow_qe_input(workflow, {\"system\": {\"tot_charge\": charge}}, unit_names=[\"pw_relax\"])\n", - " job = find_job_for_material(client, defective_supercell.id, workflow.name, ACCOUNT_ID)\n", - " if job is not None:\n", - " print(f\"♻️ n={scaling} q={charge:+d}: reusing the relaxation from job {job['_id']}\")\n", - " else:\n", - " job = create_job_for_materials([defective_supercell], workflow)\n", - " print(f\"✅ n={scaling} q={charge:+d}: created a relaxation job {job['_id']}\")\n", - " submit_jobs(client.jobs, [job[\"_id\"]])\n", - " await wait_for_jobs_to_finish_async(client.jobs, [job[\"_id\"]], poll_interval=POLL_INTERVAL)\n", - " relaxed_defective_supercells[(scaling, charge)] = get_final_structure_for_job(client, job[\"_id\"])" - ] - }, - { - "cell_type": "markdown", - "id": "43", - "metadata": {}, - "source": [ - "### 6.3. Create the Defect Formation Energy jobs, one per size and charge" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "44", - "metadata": {}, - "outputs": [], "source": [ "import pandas as pd\n", + "from mat3ra.notebooks_utils.core.entity.job.api import find_job_for_material\n", "\n", + "defect_relax_workflow_config = WorkflowStandata.filter_by_application(APPLICATION_NAME).get_by_name_first_match(\n", + " \"fixed_cell_relaxation.json\"\n", + ")\n", "job_records = []\n", "for scaling, (pristine_supercell, defective_supercell) in saved_supercell_pairs.items():\n", + " kgrid = get_scf_kgrid_for_supercell(scaling)\n", " for charge in CHARGES:\n", + " defective_material = defective_supercell\n", + " if RELAX_DEFECTIVE_MATERIAL:\n", + " relax_workflow = Workflow.create(defect_relax_workflow_config)\n", + " relax_workflow.name = f\"{relax_workflow.name} {defective_supercell.name} q={charge:+d}\"\n", + " apply_scf_kgrid(relax_workflow, kgrid, material=defective_supercell, unit_name=\"pw_relax\")\n", + " if charge:\n", + " patch_workflow_qe_input(relax_workflow, {\"system\": {\"tot_charge\": charge}}, unit_names=[\"pw_relax\"])\n", + " relax_job = find_job_for_material(\n", + " client, defective_supercell.id, relax_workflow.name, ACCOUNT_ID, kgrid=kgrid, unit_name=\"pw_relax\",\n", + " statuses=(\"finished\", \"active\", \"submitted\", \"queued\"),\n", + " )\n", + " if relax_job is None:\n", + " relax_job = create_job_for_materials([defective_supercell], relax_workflow)\n", + " submit_jobs(client.jobs, [relax_job[\"_id\"]])\n", + " print(f\"✅ n={scaling} q={charge:+d}: created a relaxation job {relax_job['_id']}\")\n", + " else:\n", + " print(f\"♻️ n={scaling} q={charge:+d}: reusing the relaxation from job {relax_job['_id']}\")\n", + " await wait_for_jobs_to_finish_async(client.jobs, [relax_job[\"_id\"]], poll_interval=POLL_INTERVAL)\n", + " defective_material = get_final_structure_for_job(client, relax_job[\"_id\"])\n", " workflow = create_defect_workflow(scaling, charge, pristine_reference_jobs[scaling])\n", - " defective_material = relaxed_defective_supercells.get((scaling, charge), defective_supercell)\n", " # Order matters: [0] defective (computed), [1] pristine (reference).\n", " job = create_job_for_materials([defective_material, pristine_supercell], workflow, [f\"charge:{charge}\"])\n", + " submit_jobs(client.jobs, [job[\"_id\"]])\n", " print(f\"n={scaling} q={charge:+d}: k-grid {get_scf_kgrid_for_supercell(scaling) or 'from KPPRA'}, \"\n", " f\"job {job['_id']}\")\n", " job_records.append({\n", @@ -732,28 +714,16 @@ }, { "cell_type": "markdown", - "id": "45", - "metadata": {}, - "source": [ - "### 6.4. Submit the jobs and monitor the statuses" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "46", + "id": "43", "metadata": {}, - "outputs": [], "source": [ - "if job_ids:\n", - " submit_jobs(client.jobs, job_ids)\n", - " print(f\"✅ Submitted {len(job_ids)} Defect Formation Energy job(s).\")" + "### 6.3. Monitor the job statuses" ] }, { "cell_type": "code", "execution_count": null, - "id": "47", + "id": "44", "metadata": {}, "outputs": [], "source": [ @@ -763,7 +733,7 @@ }, { "cell_type": "markdown", - "id": "48", + "id": "45", "metadata": {}, "source": [ "## 7. Retrieve results\n", @@ -774,7 +744,7 @@ { "cell_type": "code", "execution_count": null, - "id": "49", + "id": "46", "metadata": {}, "outputs": [], "source": [ @@ -798,7 +768,7 @@ }, { "cell_type": "markdown", - "id": "50", + "id": "47", "metadata": {}, "source": [ "### 7.2. Formation energy versus the electron chemical potential\n", @@ -808,7 +778,7 @@ { "cell_type": "code", "execution_count": null, - "id": "51", + "id": "48", "metadata": {}, "outputs": [], "source": [ @@ -832,7 +802,7 @@ { "cell_type": "code", "execution_count": null, - "id": "52", + "id": "49", "metadata": {}, "outputs": [], "source": [ @@ -856,7 +826,7 @@ }, { "cell_type": "markdown", - "id": "53", + "id": "50", "metadata": {}, "source": [ "### 7.3. Finite-size extrapolation\n", @@ -866,7 +836,7 @@ { "cell_type": "code", "execution_count": null, - "id": "54", + "id": "51", "metadata": {}, "outputs": [], "source": [ diff --git a/src/py/mat3ra/notebooks_utils/core/entity/job/api.py b/src/py/mat3ra/notebooks_utils/core/entity/job/api.py index 0df7fa1a0..6a93e3567 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/job/api.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/job/api.py @@ -127,9 +127,12 @@ def find_job_for_material( workflow_name: str, owner_id: str, statuses: Iterable[str] = ("finished",), + kgrid: Optional[List[int]] = None, + unit_name: str = "pw_scf", ) -> Optional[dict]: """ - Finds a job for a material and workflow name under the given owner, filtered by status. + Finds a job for a material and workflow name under the given owner, filtered by status and, + optionally, by the k-grid its `unit_name` unit ran on. Args: api_client (APIClient): API client instance carrying the authorization context. @@ -137,6 +140,8 @@ def find_job_for_material( workflow_name (str): Exact workflow name the job was created with. owner_id (str): Account ID the job must belong to. statuses (Iterable[str]): Job statuses that count as a match. + kgrid (List[int], optional): K-grid dimensions the job's `unit_name` unit ran on, see `get_kgrid_query`. + unit_name (str): Name of the unit the k-grid was set on. Returns: dict, optional: The matching job, or None if none exists. @@ -147,6 +152,7 @@ def find_job_for_material( "owner._id": owner_id, "workflow.name": workflow_name, "status": {"$in": list(statuses)}, + **get_kgrid_query(kgrid, unit_name), }, {"limit": 1}, ) diff --git a/tests/py/unit/core/entity/test_job_api.py b/tests/py/unit/core/entity/test_job_api.py index 5d766d24d..f560d0f48 100644 --- a/tests/py/unit/core/entity/test_job_api.py +++ b/tests/py/unit/core/entity/test_job_api.py @@ -196,3 +196,15 @@ def test_find_job_for_material_with_property_matches_the_pw_scf_kgrid(kgrid, exp client.jobs.list.assert_called_once_with( {"_material._id": MATERIAL_INITIAL["_id"], "owner._id": OWNER_ID, "status": "finished", **expected_kgrid_query} ) + + +def test_find_job_for_material_matches_the_kgrid(): + client = MagicMock() + client.jobs.list.return_value = [EXISTING_JOB] + + job = find_job_for_material(client, MATERIAL_INITIAL["_id"], RELAX_WORKFLOW_NAME, OWNER_ID, kgrid=[4, 4, 4]) + + assert job == EXISTING_JOB + assert ( + client.jobs.list.call_args.args[0]["workflow.subworkflows.units"] == KGRID_QUERY["workflow.subworkflows.units"] + ) From bfe27aa3bcea7199ff4c1682943b5910c18467b3 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 28 Sep 2026 18:26:37 -0700 Subject: [PATCH 07/18] feat(SOF-7975): treat the platform default k-grid as its own grid when reusing jobs With SCF_KGRID unset the lookups were unfiltered, so a default-grid run could reuse a reference or relaxation computed on an explicit grid. A job created without an explicit grid carries no kgrid context on its unit, and the platform then applies its default, so kgrid=None now matches exactly those jobs: get_kgrid_query asks for a unit of that name whose context has no element named kgrid (context.name $ne "kgrid"; on production it returns BLmZo5WZfFXKKTb2H and oSKAgmDmQ5HPFdZyH for the unrelaxed HfO2 cell and nothing for the 4x4x4 one). Same grid now holds both ways: an explicit run never reuses a default-grid job and a default run never reuses an explicit one. The finders default to the new ANY_KGRID, so callers that pass no kgrid, the SOF-8044 notebooks among them, keep the unfiltered lookup. The reuse line prints "k-grid: platform default", and the section texts say the platform default counts as the grid when SCF_KGRID is not set. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../defect_formation_energy_charged.ipynb | 8 +++--- .../notebooks_utils/core/entity/job/api.py | 28 +++++++++++-------- .../core/entity/material/api.py | 13 ++++++--- tests/py/unit/core/entity/test_job_api.py | 19 +++++++++---- .../py/unit/core/entity/test_material_api.py | 21 ++++++++------ 5 files changed, 56 insertions(+), 33 deletions(-) diff --git a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb index 4c910e751..0bb983e8d 100644 --- a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb @@ -287,7 +287,7 @@ "metadata": {}, "source": [ "### 3.3. Define how jobs are created\n", - "Every job is tagged with its charge; the pristine references run neutral, tagged `charge:0`. A reference is reused when a finished `charge:0` job on that pristine cell already reported the property (`total_energy`, `band_gaps`), on the same k-grid when `SCF_KGRID` is set. The Defect Formation Energy workflow reads each reference from the one job the notebook names to it, so no other calculation on the same cell can stand in; the relaxation carries no charge tag and is never taken for one." + "Every job is tagged with its charge; the pristine references run neutral, tagged `charge:0`. A reference is reused when a finished `charge:0` job on that pristine cell already reported the property (`total_energy`, `band_gaps`), on the same k-grid (the platform default when `SCF_KGRID` is not set). The Defect Formation Energy workflow reads each reference from the one job the notebook names to it, so no other calculation on the same cell can stand in; the relaxation carries no charge tag and is never taken for one." ] }, { @@ -353,7 +353,7 @@ "metadata": {}, "source": [ "### 4.2. Relax the pristine cell (optional)\n", - "With `RELAX_PRISTINE_MATERIAL`, the lattice and the atomic positions of the pristine cell are relaxed with the Variable-cell Relaxation workflow, or taken from an earlier finished variable-cell relaxation of the same structure, on the same k-grid when `SCF_KGRID` is set, and the supercells below are built from the relaxed cell." + "With `RELAX_PRISTINE_MATERIAL`, the lattice and the atomic positions of the pristine cell are relaxed with the Variable-cell Relaxation workflow, or taken from an earlier finished variable-cell relaxation of the same structure, on the same k-grid (the platform default when `SCF_KGRID` is not set), and the supercells below are built from the relaxed cell." ] }, { @@ -623,7 +623,7 @@ " client, pristine_supercell.id, property_name, ACCOUNT_ID, tags=[\"charge:0\"], kgrid=kgrid\n", " )\n", " if job is not None:\n", - " print(f\"♻️ n={scaling}: reusing {property_name} from job {job['name']}, k-grid {kgrid or 'from KPPRA'}\")\n", + " print(f\"♻️ n={scaling}: reusing {property_name} from job {job['name']}, k-grid: {kgrid or 'platform default'}\")\n", " else:\n", " workflow = Workflow.create(workflow_config)\n", " workflow.name = f\"{workflow.name} {pristine_supercell.name}\"\n", @@ -655,7 +655,7 @@ "metadata": {}, "source": [ "### 6.2. Create and submit the Defect Formation Energy jobs, one per size and charge\n", - "With `RELAX_DEFECTIVE_MATERIAL`, the atomic positions of each defective supercell are first relaxed in that charge state with the Fixed-cell Relaxation workflow, and the Defect Formation Energy job of that size and charge is submitted on the relaxed structure as soon as its relaxation has finished, before the next charge state is relaxed. An earlier relaxation of the same cell and charge (on the same k-grid when `SCF_KGRID` is set) is reused, and waited on while it is still running. The lattice stays that of the pristine supercell." + "With `RELAX_DEFECTIVE_MATERIAL`, the atomic positions of each defective supercell are first relaxed in that charge state with the Fixed-cell Relaxation workflow, and the Defect Formation Energy job of that size and charge is submitted on the relaxed structure as soon as its relaxation has finished, before the next charge state is relaxed. An earlier relaxation of the same cell, charge and k-grid (the platform default when `SCF_KGRID` is not set) is reused, and waited on while it is still running. The lattice stays that of the pristine supercell." ] }, { diff --git a/src/py/mat3ra/notebooks_utils/core/entity/job/api.py b/src/py/mat3ra/notebooks_utils/core/entity/job/api.py index 6a93e3567..2a77f9a89 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/job/api.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/job/api.py @@ -4,6 +4,7 @@ from mat3ra.api_client import APIClient, JobEndpoints MATERIALS_SET_ENTITY_CLASS = "Material" +ANY_KGRID = "any" def save_files(job_id: str, job_endpoint: JobEndpoints, filename_on_cloud: str, filename_on_disk: str) -> None: @@ -127,7 +128,7 @@ def find_job_for_material( workflow_name: str, owner_id: str, statuses: Iterable[str] = ("finished",), - kgrid: Optional[List[int]] = None, + kgrid: Union[List[int], str, None] = ANY_KGRID, unit_name: str = "pw_scf", ) -> Optional[dict]: """ @@ -140,7 +141,8 @@ def find_job_for_material( workflow_name (str): Exact workflow name the job was created with. owner_id (str): Account ID the job must belong to. statuses (Iterable[str]): Job statuses that count as a match. - kgrid (List[int], optional): K-grid dimensions the job's `unit_name` unit ran on, see `get_kgrid_query`. + kgrid (List[int], optional): K-grid the job's `unit_name` unit ran on, None for the platform default, + see `get_kgrid_query`. unit_name (str): Name of the unit the k-grid was set on. Returns: @@ -159,16 +161,19 @@ def find_job_for_material( return existing[0] if existing else None -def get_kgrid_query(kgrid: Optional[List[int]], unit_name: str = "pw_scf") -> Dict[str, Any]: +def get_kgrid_query(kgrid: Union[List[int], str, None], unit_name: str = "pw_scf") -> Dict[str, Any]: """ - `jobs.list` condition for jobs whose `unit_name` unit ran on `kgrid` (empty when `kgrid` is None), matched where - `apply_scf_kgrid` sets it: `workflow.subworkflows[].units[name].context[name="kgrid"].data.dimensions`. - A job created without an explicit k-grid has no such context and never matches. + `jobs.list` condition on the k-grid of the job's `unit_name` unit, read where `apply_scf_kgrid` sets it: + `workflow.subworkflows[].units[name].context[name="kgrid"].data.dimensions`. `kgrid` is the dimensions to match, + None for the platform's default grid (a unit without a kgrid context), or ANY_KGRID for no condition. """ - if kgrid is None: + if kgrid == ANY_KGRID: return {} - kgrid_context = {"$elemMatch": {"name": "kgrid", "data.dimensions": list(kgrid)}} - return {"workflow.subworkflows.units": {"$elemMatch": {"name": unit_name, "context": kgrid_context}}} + if kgrid is None: + unit = {"name": unit_name, "context.name": {"$ne": "kgrid"}} + else: + unit = {"name": unit_name, "context": {"$elemMatch": {"name": "kgrid", "data.dimensions": list(kgrid)}}} + return {"workflow.subworkflows.units": {"$elemMatch": unit}} def find_job_for_material_with_property( @@ -177,7 +182,7 @@ def find_job_for_material_with_property( property_name: str, owner_id: str, tags: Optional[List[str]] = None, - kgrid: Optional[List[int]] = None, + kgrid: Union[List[int], str, None] = ANY_KGRID, ) -> Optional[dict]: """ Finds a finished job on a material that reported the given property, optionally among the jobs @@ -190,7 +195,8 @@ def find_job_for_material_with_property( property_name (str): Property the job must have reported, e.g. "total_energy". owner_id (str): Account ID the job must belong to. tags (List[str], optional): Tags the job must all carry. - kgrid (List[int], optional): K-grid dimensions the job's `pw_scf` unit ran on, see `get_kgrid_query`. + kgrid (List[int], optional): K-grid the job's `pw_scf` unit ran on, None for the platform default, + see `get_kgrid_query`. Returns: dict, optional: The first matching job, or None if none exists. diff --git a/src/py/mat3ra/notebooks_utils/core/entity/material/api.py b/src/py/mat3ra/notebooks_utils/core/entity/material/api.py index d32c82c8e..7c09d0f02 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/material/api.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/material/api.py @@ -1,12 +1,12 @@ import os import re -from typing import Any, Dict, List, Optional +from typing import Any, Dict, List, Optional, Union from mat3ra.api_client import APIClient from mat3ra.made.material import Material from mat3ra.prode import PropertyName -from ..job.api import get_kgrid_query +from ..job.api import ANY_KGRID, get_kgrid_query from ..property.api import get_properties_for_job from .analysis import get_slab_bulk_crystal, resolve_bulk_query_from_crystal from .io import load_materials_from_folder @@ -99,7 +99,11 @@ def get_final_structure_for_job(api_client: APIClient, job_id: str) -> Material: def find_relaxed_material( - api_client: APIClient, material, owner_id: str, kgrid: Optional[List[int]] = None, unit_name: str = "pw_scf" + api_client: APIClient, + material, + owner_id: str, + kgrid: Union[List[int], str, None] = ANY_KGRID, + unit_name: str = "pw_scf", ) -> Optional[Material]: """ Finds a relaxed version of a material: the final structure of a finished job on a material @@ -110,7 +114,8 @@ def find_relaxed_material( api_client (APIClient): API client instance carrying the authorization context. material: mat3ra-made Material object (must have a .hash property). owner_id (str): Account ID under which to search. - kgrid (List[int], optional): K-grid dimensions the relaxation ran on, see `get_kgrid_query`. + kgrid (List[int], optional): K-grid the relaxation ran on, None for the platform default, + see `get_kgrid_query`. unit_name (str): Name of the relaxation unit, e.g. "pw_vc-relax". Returns: diff --git a/tests/py/unit/core/entity/test_job_api.py b/tests/py/unit/core/entity/test_job_api.py index f560d0f48..e78490073 100644 --- a/tests/py/unit/core/entity/test_job_api.py +++ b/tests/py/unit/core/entity/test_job_api.py @@ -3,6 +3,7 @@ import pytest from mat3ra.notebooks_utils.core.entity.job.api import ( + ANY_KGRID, create_job, find_job_for_material, find_job_for_material_with_property, @@ -182,9 +183,15 @@ def test_find_job_for_material_with_property_returns_none_when_no_job_reported_i "$elemMatch": {"name": "pw_scf", "context": {"$elemMatch": {"name": "kgrid", "data.dimensions": [4, 4, 4]}}} } } +DEFAULT_KGRID_QUERY: Dict[str, Any] = { + "workflow.subworkflows.units": {"$elemMatch": {"name": "pw_scf", "context.name": {"$ne": "kgrid"}}} +} -@pytest.mark.parametrize(("kgrid", "expected_kgrid_query"), [(None, {}), ([4, 4, 4], KGRID_QUERY)]) +@pytest.mark.parametrize( + ("kgrid", "expected_kgrid_query"), + [(ANY_KGRID, {}), (None, DEFAULT_KGRID_QUERY), ([4, 4, 4], KGRID_QUERY)], +) def test_find_job_for_material_with_property_matches_the_pw_scf_kgrid(kgrid, expected_kgrid_query): client = MagicMock() client.jobs.list.return_value = [EXISTING_JOB] @@ -198,13 +205,13 @@ def test_find_job_for_material_with_property_matches_the_pw_scf_kgrid(kgrid, exp ) -def test_find_job_for_material_matches_the_kgrid(): +@pytest.mark.parametrize(("kgrid", "expected_kgrid_query"), [(None, DEFAULT_KGRID_QUERY), ([4, 4, 4], KGRID_QUERY)]) +def test_find_job_for_material_matches_the_kgrid(kgrid, expected_kgrid_query): client = MagicMock() client.jobs.list.return_value = [EXISTING_JOB] - job = find_job_for_material(client, MATERIAL_INITIAL["_id"], RELAX_WORKFLOW_NAME, OWNER_ID, kgrid=[4, 4, 4]) + job = find_job_for_material(client, MATERIAL_INITIAL["_id"], RELAX_WORKFLOW_NAME, OWNER_ID, kgrid=kgrid) assert job == EXISTING_JOB - assert ( - client.jobs.list.call_args.args[0]["workflow.subworkflows.units"] == KGRID_QUERY["workflow.subworkflows.units"] - ) + kgrid_key = "workflow.subworkflows.units" + assert client.jobs.list.call_args.args[0][kgrid_key] == expected_kgrid_query[kgrid_key] diff --git a/tests/py/unit/core/entity/test_material_api.py b/tests/py/unit/core/entity/test_material_api.py index c6c39d558..81d4b1b8e 100644 --- a/tests/py/unit/core/entity/test_material_api.py +++ b/tests/py/unit/core/entity/test_material_api.py @@ -359,14 +359,24 @@ def test_find_relaxed_material_skips_a_final_structure_with_the_same_hash(): assert client.properties.get_for_job.call_count == 2 -def test_find_relaxed_material_matches_the_relaxation_kgrid(): +@pytest.mark.parametrize( + ("kgrid", "expected_unit_condition"), + [ + (None, {"name": "pw_vc-relax", "context.name": {"$ne": "kgrid"}}), + ( + [4, 4, 4], + {"name": "pw_vc-relax", "context": {"$elemMatch": {"name": "kgrid", "data.dimensions": [4, 4, 4]}}}, + ), + ], +) +def test_find_relaxed_material_matches_the_relaxation_kgrid(kgrid, expected_unit_condition): client = MagicMock() client.materials.list.return_value = [SAVED_DEFECTIVE] client.jobs.list.return_value = [FINISHED_JOB] client.properties.get_for_job.return_value = [{"materialId": "m-relaxed"}] client.materials.get.return_value = RELAXED_MATERIAL_DOC - relaxed = find_relaxed_material(client, DEFECTIVE_MATERIAL, OWNER_ID, kgrid=[4, 4, 4], unit_name="pw_vc-relax") + relaxed = find_relaxed_material(client, DEFECTIVE_MATERIAL, OWNER_ID, kgrid=kgrid, unit_name="pw_vc-relax") assert relaxed is not None assert relaxed.name == "B-vacancy h-BN relaxed" @@ -375,12 +385,7 @@ def test_find_relaxed_material_matches_the_relaxation_kgrid(): "_material._id": {"$in": [SAVED_DEFECTIVE["_id"]]}, "owner._id": OWNER_ID, "status": "finished", - "workflow.subworkflows.units": { - "$elemMatch": { - "name": "pw_vc-relax", - "context": {"$elemMatch": {"name": "kgrid", "data.dimensions": [4, 4, 4]}}, - } - }, + "workflow.subworkflows.units": {"$elemMatch": expected_unit_condition}, } ) From 78044e01f9eb578b76ce8f9f029ca63aa2b208e5 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 28 Sep 2026 18:35:03 -0700 Subject: [PATCH 08/18] fix(SOF-7975): keep the defect formation energy notebook neutral The CHARGE parameter and its tot_charge patch were SOF-7917's first attempt at charged defects, now superseded by defect_formation_energy_charged.ipynb, and the header contradicted itself: it announced the neutral formation energy while its formula carried [q] and + q(E_VBM + E_F), and a paragraph then said that term was never computed. CHARGE, its comment, the tot_charge patch and the now-unused patch_workflow_qe_input import go; the formula reads E_defect = E_defective - E_pristine - sum_i dN_i mu_i; the CHARGE paragraph goes, and the opening paragraph points to the charged notebook the way that notebook points back here. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../workflows/defect_formation_energy.ipynb | 20 ++----------------- 1 file changed, 2 insertions(+), 18 deletions(-) diff --git a/other/materials_designer/workflows/defect_formation_energy.ipynb b/other/materials_designer/workflows/defect_formation_energy.ipynb index 10651a132..5891f40c0 100644 --- a/other/materials_designer/workflows/defect_formation_energy.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy.ipynb @@ -7,7 +7,7 @@ "source": [ "# Defect Formation Energy\n", "\n", - "Calculate the neutral defect formation energy (eV) of a defective supercell using a multi-material DFT workflow on the Mat3ra platform.\n", + "Calculate the neutral defect formation energy (eV) of a defective supercell using a multi-material DFT workflow on the Mat3ra platform. For charged defects, supercell-size series and formation energy vs Fermi level, use [Defect Formation Energy of Charged Defects](defect_formation_energy_charged.ipynb).\n", "\n", "The workflow takes **two materials, in order**:\n", "\n", @@ -18,9 +18,7 @@ "\n", "Formula:\n", "\n", - "$$E_{\\text{defect}} = E_{\\text{defective}}[q] - E_{\\text{pristine}} - \\sum_i \\Delta N_i\\, \\mu_i + q\\,(E_{\\text{VBM}} + E_F) \\quad [\\text{eV}]$$\n", - "\n", - "Set `CHARGE` (cell 1.3) to a non-zero value to charge the defective supercell (`tot_charge` in the QE `&SYSTEM` namelist, compensated by a uniform jellium background). **This notebook does not compute a charged-defect finite-size correction** (e.g. Freysoldt-Neugebauer-Van de Walle) or track the Fermi-level term $q(E_{\\text{VBM}}+E_F)$ -- for `CHARGE = 0` (the default) neither is needed and the value below is directly physical; for `CHARGE != 0` the reported value is the raw, uncorrected total-energy difference only.\n", + "$$E_{\\text{defect}} = E_{\\text{defective}} - E_{\\text{pristine}} - \\sum_i \\Delta N_i\\, \\mu_i \\quad [\\text{eV}]$$\n", "\n", "where $\\Delta N_i = \\text{count}_i(\\text{defective}) - \\text{count}_i(\\text{pristine})$ is the per-species atom-count change and $\\mu_i = E_{\\text{elemental},i} / n_{\\text{atoms},i}$.\n", "\n", @@ -125,15 +123,6 @@ "# K-grid for the defective-cell SCF (if not set, KPPRA is used by default)\n", "SCF_KGRID = None # e.g. [4, 4, 4]\n", "\n", - "# Net charge on the defective supercell (e.g. -3 for a triply negatively charged\n", - "# defect, +1 for a singly positively charged one). 0 = neutral defect (default).\n", - "# NOTE: a non-zero charge requires a charged-defect finite-size correction (e.g.\n", - "# Freysoldt-Neugebauer-Van de Walle) to remove the spurious electrostatic\n", - "# interaction between the charged defect and its periodic images -- this\n", - "# notebook does NOT compute that correction, so CHARGE != 0 results are raw,\n", - "# uncorrected values only.\n", - "CHARGE = 0 # e.g. -3, +1\n", - "\n", "# Whose total_energy properties to consider for the pristine reference:\n", "# \"public\" (any owner, highest precision wins), \"curators\" (only curators'),\n", "# or \"my_account\" (curators' or your own).\n", @@ -434,7 +423,6 @@ "from mat3ra.wode.workflows import Workflow\n", "from mat3ra.wode.context.providers import PointsGridDataProvider\n", "from mat3ra.notebooks_utils.ipython.entity.workflow.visualize import visualize_workflow\n", - "from mat3ra.notebooks_utils.workflow import patch_workflow_qe_input\n", "\n", "defect_workflow_config = WorkflowStandata.filter_by_application(app.name).get_by_name_first_match(\n", " WORKFLOW_SEARCH_TERM\n", @@ -453,10 +441,6 @@ " unit.add_context(new_context)\n", " subworkflow.set_unit(unit)\n", "\n", - "# Net charge on the defective-cell SCF: adds `tot_charge` to the &SYSTEM namelist\n", - "if CHARGE:\n", - " patch_workflow_qe_input(defect_workflow, {\"system\": {\"tot_charge\": CHARGE}}, unit_names=[\"pw_scf\"])\n", - "\n", "visualize_workflow(defect_workflow)" ] }, From 67245218416579dfccbccc8969ef43cfcbf737a4 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 28 Sep 2026 19:05:57 -0700 Subject: [PATCH 09/18] fix(SOF-7975): resolve SCF_KGRID = None to the platform default grid, keep the notebook's cell layout, tag the defect relaxation With SCF_KGRID unset, each supercell size now gets one explicit grid, the platform default of the pristine supercell, computed the way the platform does (made ReciprocalLattice.get_dimensions_from_points_count at web-app defaultKPPRA 10 per atom: 12-atom HfO2 mp-352 gives 1x1x1, the Si default gives 2x2x2 where the 1-atom defective cell alone would get 3x3x3); the pristine references, the pristine and defective relaxations and the defect jobs all run on it and every lookup matches it, so the None-as-own-grid rule and the ANY_KGRID sentinel go and the finders take kgrid=None as no condition, as the SOF-8044 callers expect. The notebook keeps origin/main's 53 cells in order, because the web-app Cypress feature asserts on cells 51 and 53: the defective-cell relaxation stays inside the defect-job cell and the existing submit cell again submits the defect jobs. The defective relaxation job is tagged charge:q (its lookup stays by name so a running one can be waited on), the pristine relaxation prints a created or reused line with its id, the compute line prints the time limit, the k-grid prints show the grid itself, and the defect relaxation workflow loads through app.name like the other workflows in section 5-6. get_kgrid_query is tested directly for pw_scf and pw_relax and each finder keeps one case showing its condition is merged, so dropping unit_name in find_job_for_material now fails a test. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../defect_formation_energy_charged.ipynb | 72 ++++++++++++------- .../notebooks_utils/core/entity/job/api.py | 28 +++----- .../core/entity/material/api.py | 13 ++-- tests/py/unit/core/entity/test_job_api.py | 36 +++++----- .../py/unit/core/entity/test_material_api.py | 24 ++----- 5 files changed, 87 insertions(+), 86 deletions(-) diff --git a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb index 0bb983e8d..2bb593f15 100644 --- a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb @@ -24,7 +24,7 @@ "\n", "$E_{\\text{tot}}[\\text{bulk}]$, $E_{\\text{VBM}}$ and $E_{\\text{gap}}$ come from a Total Energy and a Band Gap job on the pristine supercell; $\\mu_i$ from the Standata elemental reference materials, as in [Formation Energy](formation_energy.ipynb). The charge $q$ enters as `tot_charge` in the QE `&SYSTEM` namelist and is compensated by a uniform jellium background.\n", "\n", - "Every Total Energy, Band Gap and defect job is tagged with its charge, `charge:`. The pristine references run neutral, tagged `charge:0`, and the notebook names those jobs to the workflow, which reads its references from them alone, so a charged calculation on the same pristine cell is never taken for the neutral reference. The pristine cell can optionally be relaxed first (variable cell), and every supercell is then built from the relaxed cell. With `RELAX_DEFECTIVE_MATERIAL`, each defective supercell is also relaxed at fixed cell in its charge state before its SCF.\n", + "Every Total Energy, Band Gap and defect job is tagged with its charge, `charge:`. The pristine references run neutral, tagged `charge:0`, and the notebook names those jobs to the workflow, which reads its references from them alone, so a charged calculation on the same pristine cell is never taken for the neutral reference. The pristine cell can optionally be relaxed first (variable cell), and every supercell is then built from the relaxed cell. With `RELAX_DEFECTIVE_MATERIAL`, each defective supercell is also relaxed at fixed cell in its charge state before its SCF, in a job tagged with that charge.\n", "\n", "No analytical finite-size correction (Freysoldt-Neugebauer-Van de Walle, Makov-Payne) is applied. Running several `SUPERCELL_SCALINGS` and extrapolating $E_f(L\\to\\infty)$ from the fitted image-charge terms takes its place; a single size gives the uncorrected value for that cell. No potential alignment ($q\\,\\Delta V$) is applied either.\n", "\n", @@ -278,7 +278,7 @@ " cluster = clusters[0]\n", "\n", "compute = Compute(cluster=cluster, queue=QUEUE_NAME, ppn=PPN, timeLimit=TIME_LIMIT)\n", - "print(f\"Using cluster: {compute.cluster.hostname}, queue: {QUEUE_NAME}, ppn: {PPN}\")" + "print(f\"Using cluster: {compute.cluster.hostname}, queue: {QUEUE_NAME}, ppn: {PPN}, time limit: {TIME_LIMIT}\")" ] }, { @@ -297,9 +297,12 @@ "metadata": {}, "outputs": [], "source": [ + "from mat3ra.made.reciprocal.lattice_reciprocal import ReciprocalLattice\n", "from mat3ra.notebooks_utils.api.job import submit_jobs, wait_for_jobs_to_finish_async\n", "from mat3ra.notebooks_utils.job import create_job\n", "\n", + "PLATFORM_DEFAULT_KPPRA = 10 # defaultKPPRA in web-app src/application/imports/app_settings/settings.ts\n", + "\n", "\n", "def create_job_for_materials(materials, workflow, tags=None):\n", " job = create_job(\n", @@ -312,7 +315,13 @@ " compute=compute.to_dict(),\n", " tags=tags,\n", " )\n", - " return job[0] if isinstance(job, list) else job" + " return job[0] if isinstance(job, list) else job\n", + "\n", + "\n", + "def get_platform_default_kgrid(material):\n", + " \"\"\"The k-grid the platform gives `material` when none is set, from PLATFORM_DEFAULT_KPPRA.\"\"\"\n", + " number_of_kpoints = PLATFORM_DEFAULT_KPPRA / len(material.basis.elements.ids)\n", + " return ReciprocalLattice(**material.lattice.to_dict()).get_dimensions_from_points_count(number_of_kpoints)" ] }, { @@ -374,8 +383,9 @@ "\n", "if RELAX_PRISTINE_MATERIAL:\n", " saved_pristine_material = Material.create(get_or_create_material(client, pristine_material, ACCOUNT_ID))\n", + " kgrid = SCF_KGRID or get_platform_default_kgrid(saved_pristine_material)\n", " relaxed_pristine_material = find_relaxed_material(\n", - " client, saved_pristine_material, ACCOUNT_ID, kgrid=SCF_KGRID, unit_name=\"pw_vc-relax\"\n", + " client, saved_pristine_material, ACCOUNT_ID, kgrid=kgrid, unit_name=\"pw_vc-relax\"\n", " )\n", " if relaxed_pristine_material is None:\n", " relax_workflow_config = WorkflowStandata.filter_by_application(APPLICATION_NAME).get_by_name_first_match(\n", @@ -383,11 +393,14 @@ " )\n", " relax_workflow = Workflow.create(relax_workflow_config)\n", " relax_workflow.name = f\"{relax_workflow.name} {saved_pristine_material.name}\"\n", - " apply_scf_kgrid(relax_workflow, SCF_KGRID, material=saved_pristine_material, unit_name=\"pw_vc-relax\")\n", + " apply_scf_kgrid(relax_workflow, kgrid, material=saved_pristine_material, unit_name=\"pw_vc-relax\")\n", " relax_job = create_job_for_materials([saved_pristine_material], relax_workflow)\n", " submit_jobs(client.jobs, [relax_job[\"_id\"]])\n", + " print(f\"✅ created a relaxation job {relax_job['_id']}\")\n", " await wait_for_jobs_to_finish_async(client.jobs, [relax_job[\"_id\"]], poll_interval=POLL_INTERVAL)\n", " relaxed_pristine_material = get_final_structure_for_job(client, relax_job[\"_id\"])\n", + " else:\n", + " print(f\"♻️ reusing the relaxed structure {relaxed_pristine_material.id}\")\n", " relaxed_pristine_material.name = f\"{pristine_material.name} relaxed\"\n", " print(f\"Lattice constant a: {pristine_material.lattice.a:.4f} Å as given, \"\n", " f\"{relaxed_pristine_material.lattice.a:.4f} Å relaxed\")\n", @@ -553,8 +566,10 @@ "\n", "\n", "def get_scf_kgrid_for_supercell(scaling):\n", - " \"\"\"SCF_KGRID divided by the scaling, rounded, never below 1; None when SCF_KGRID is not set.\"\"\"\n", - " return None if SCF_KGRID is None else [max(1, round(dimension / scaling)) for dimension in SCF_KGRID]\n", + " \"\"\"SCF_KGRID divided by the scaling, rounded, never below 1; without SCF_KGRID, the pristine supercell's default.\"\"\"\n", + " if SCF_KGRID is None:\n", + " return get_platform_default_kgrid(saved_supercell_pairs[scaling][0])\n", + " return [max(1, round(dimension / scaling)) for dimension in SCF_KGRID]\n", "\n", "\n", "def create_defect_workflow(scaling, charge, reference_jobs):\n", @@ -623,14 +638,14 @@ " client, pristine_supercell.id, property_name, ACCOUNT_ID, tags=[\"charge:0\"], kgrid=kgrid\n", " )\n", " if job is not None:\n", - " print(f\"♻️ n={scaling}: reusing {property_name} from job {job['name']}, k-grid: {kgrid or 'platform default'}\")\n", + " print(f\"♻️ n={scaling}: reusing {property_name} from job {job['name']}, k-grid {kgrid}\")\n", " else:\n", " workflow = Workflow.create(workflow_config)\n", " workflow.name = f\"{workflow.name} {pristine_supercell.name}\"\n", " apply_scf_kgrid(workflow, kgrid, material=pristine_supercell)\n", " job = create_job_for_materials([pristine_supercell], workflow, [\"charge:0\"])\n", " pristine_job_ids.append(job[\"_id\"])\n", - " print(f\"✅ n={scaling}: created a {property_name} job {job['_id']}, k-grid {kgrid or 'from KPPRA'}\")\n", + " print(f\"✅ n={scaling}: created a {property_name} job {job['_id']}, k-grid {kgrid}\")\n", " pristine_reference_jobs[scaling][property_name] = job\n", "\n", "if pristine_job_ids:\n", @@ -654,8 +669,7 @@ "id": "41", "metadata": {}, "source": [ - "### 6.2. Create and submit the Defect Formation Energy jobs, one per size and charge\n", - "With `RELAX_DEFECTIVE_MATERIAL`, the atomic positions of each defective supercell are first relaxed in that charge state with the Fixed-cell Relaxation workflow, and the Defect Formation Energy job of that size and charge is submitted on the relaxed structure as soon as its relaxation has finished, before the next charge state is relaxed. An earlier relaxation of the same cell, charge and k-grid (the platform default when `SCF_KGRID` is not set) is reused, and waited on while it is still running. The lattice stays that of the pristine supercell." + "### 6.2. Create the Defect Formation Energy jobs, one per size and charge" ] }, { @@ -668,7 +682,7 @@ "import pandas as pd\n", "from mat3ra.notebooks_utils.core.entity.job.api import find_job_for_material\n", "\n", - "defect_relax_workflow_config = WorkflowStandata.filter_by_application(APPLICATION_NAME).get_by_name_first_match(\n", + "defect_relax_workflow_config = WorkflowStandata.filter_by_application(app.name).get_by_name_first_match(\n", " \"fixed_cell_relaxation.json\"\n", ")\n", "job_records = []\n", @@ -687,7 +701,7 @@ " statuses=(\"finished\", \"active\", \"submitted\", \"queued\"),\n", " )\n", " if relax_job is None:\n", - " relax_job = create_job_for_materials([defective_supercell], relax_workflow)\n", + " relax_job = create_job_for_materials([defective_supercell], relax_workflow, [f\"charge:{charge}\"])\n", " submit_jobs(client.jobs, [relax_job[\"_id\"]])\n", " print(f\"✅ n={scaling} q={charge:+d}: created a relaxation job {relax_job['_id']}\")\n", " else:\n", @@ -697,9 +711,7 @@ " workflow = create_defect_workflow(scaling, charge, pristine_reference_jobs[scaling])\n", " # Order matters: [0] defective (computed), [1] pristine (reference).\n", " job = create_job_for_materials([defective_material, pristine_supercell], workflow, [f\"charge:{charge}\"])\n", - " submit_jobs(client.jobs, [job[\"_id\"]])\n", - " print(f\"n={scaling} q={charge:+d}: k-grid {get_scf_kgrid_for_supercell(scaling) or 'from KPPRA'}, \"\n", - " f\"job {job['_id']}\")\n", + " print(f\"n={scaling} q={charge:+d}: k-grid {kgrid}, job {job['_id']}\")\n", " job_records.append({\n", " \"scaling\": scaling,\n", " \"charge\": charge,\n", @@ -717,7 +729,7 @@ "id": "43", "metadata": {}, "source": [ - "### 6.3. Monitor the job statuses" + "### 6.3. Submit the jobs and monitor the statuses" ] }, { @@ -726,6 +738,18 @@ "id": "44", "metadata": {}, "outputs": [], + "source": [ + "if job_ids:\n", + " submit_jobs(client.jobs, job_ids)\n", + " print(f\"✅ Submitted {len(job_ids)} Defect Formation Energy job(s).\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "45", + "metadata": {}, + "outputs": [], "source": [ "if job_ids:\n", " await wait_for_jobs_to_finish_async(client.jobs, job_ids, poll_interval=POLL_INTERVAL)" @@ -733,7 +757,7 @@ }, { "cell_type": "markdown", - "id": "45", + "id": "46", "metadata": {}, "source": [ "## 7. Retrieve results\n", @@ -744,7 +768,7 @@ { "cell_type": "code", "execution_count": null, - "id": "46", + "id": "47", "metadata": {}, "outputs": [], "source": [ @@ -768,7 +792,7 @@ }, { "cell_type": "markdown", - "id": "47", + "id": "48", "metadata": {}, "source": [ "### 7.2. Formation energy versus the electron chemical potential\n", @@ -778,7 +802,7 @@ { "cell_type": "code", "execution_count": null, - "id": "48", + "id": "49", "metadata": {}, "outputs": [], "source": [ @@ -802,7 +826,7 @@ { "cell_type": "code", "execution_count": null, - "id": "49", + "id": "50", "metadata": {}, "outputs": [], "source": [ @@ -826,7 +850,7 @@ }, { "cell_type": "markdown", - "id": "50", + "id": "51", "metadata": {}, "source": [ "### 7.3. Finite-size extrapolation\n", @@ -836,7 +860,7 @@ { "cell_type": "code", "execution_count": null, - "id": "51", + "id": "52", "metadata": {}, "outputs": [], "source": [ diff --git a/src/py/mat3ra/notebooks_utils/core/entity/job/api.py b/src/py/mat3ra/notebooks_utils/core/entity/job/api.py index 2a77f9a89..57c33dcd4 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/job/api.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/job/api.py @@ -4,7 +4,6 @@ from mat3ra.api_client import APIClient, JobEndpoints MATERIALS_SET_ENTITY_CLASS = "Material" -ANY_KGRID = "any" def save_files(job_id: str, job_endpoint: JobEndpoints, filename_on_cloud: str, filename_on_disk: str) -> None: @@ -128,7 +127,7 @@ def find_job_for_material( workflow_name: str, owner_id: str, statuses: Iterable[str] = ("finished",), - kgrid: Union[List[int], str, None] = ANY_KGRID, + kgrid: Optional[List[int]] = None, unit_name: str = "pw_scf", ) -> Optional[dict]: """ @@ -141,8 +140,7 @@ def find_job_for_material( workflow_name (str): Exact workflow name the job was created with. owner_id (str): Account ID the job must belong to. statuses (Iterable[str]): Job statuses that count as a match. - kgrid (List[int], optional): K-grid the job's `unit_name` unit ran on, None for the platform default, - see `get_kgrid_query`. + kgrid (List[int], optional): Exact k-grid dimensions the job's `unit_name` unit ran on; None for no condition. unit_name (str): Name of the unit the k-grid was set on. Returns: @@ -161,19 +159,16 @@ def find_job_for_material( return existing[0] if existing else None -def get_kgrid_query(kgrid: Union[List[int], str, None], unit_name: str = "pw_scf") -> Dict[str, Any]: +def get_kgrid_query(kgrid: Optional[List[int]], unit_name: str = "pw_scf") -> Dict[str, Any]: """ - `jobs.list` condition on the k-grid of the job's `unit_name` unit, read where `apply_scf_kgrid` sets it: - `workflow.subworkflows[].units[name].context[name="kgrid"].data.dimensions`. `kgrid` is the dimensions to match, - None for the platform's default grid (a unit without a kgrid context), or ANY_KGRID for no condition. + `jobs.list` condition for jobs whose `unit_name` unit ran on `kgrid` (empty when `kgrid` is None), matched where + `apply_scf_kgrid` sets it: `workflow.subworkflows[].units[name].context[name="kgrid"].data.dimensions`. + A job created without an explicit k-grid has no such context and never matches. """ - if kgrid == ANY_KGRID: - return {} if kgrid is None: - unit = {"name": unit_name, "context.name": {"$ne": "kgrid"}} - else: - unit = {"name": unit_name, "context": {"$elemMatch": {"name": "kgrid", "data.dimensions": list(kgrid)}}} - return {"workflow.subworkflows.units": {"$elemMatch": unit}} + return {} + kgrid_context = {"$elemMatch": {"name": "kgrid", "data.dimensions": list(kgrid)}} + return {"workflow.subworkflows.units": {"$elemMatch": {"name": unit_name, "context": kgrid_context}}} def find_job_for_material_with_property( @@ -182,7 +177,7 @@ def find_job_for_material_with_property( property_name: str, owner_id: str, tags: Optional[List[str]] = None, - kgrid: Union[List[int], str, None] = ANY_KGRID, + kgrid: Optional[List[int]] = None, ) -> Optional[dict]: """ Finds a finished job on a material that reported the given property, optionally among the jobs @@ -195,8 +190,7 @@ def find_job_for_material_with_property( property_name (str): Property the job must have reported, e.g. "total_energy". owner_id (str): Account ID the job must belong to. tags (List[str], optional): Tags the job must all carry. - kgrid (List[int], optional): K-grid the job's `pw_scf` unit ran on, None for the platform default, - see `get_kgrid_query`. + kgrid (List[int], optional): Exact k-grid dimensions the job's `pw_scf` unit ran on; None for no condition. Returns: dict, optional: The first matching job, or None if none exists. diff --git a/src/py/mat3ra/notebooks_utils/core/entity/material/api.py b/src/py/mat3ra/notebooks_utils/core/entity/material/api.py index 7c09d0f02..4848f1a65 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/material/api.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/material/api.py @@ -1,12 +1,12 @@ import os import re -from typing import Any, Dict, List, Optional, Union +from typing import Any, Dict, List, Optional from mat3ra.api_client import APIClient from mat3ra.made.material import Material from mat3ra.prode import PropertyName -from ..job.api import ANY_KGRID, get_kgrid_query +from ..job.api import get_kgrid_query from ..property.api import get_properties_for_job from .analysis import get_slab_bulk_crystal, resolve_bulk_query_from_crystal from .io import load_materials_from_folder @@ -99,11 +99,7 @@ def get_final_structure_for_job(api_client: APIClient, job_id: str) -> Material: def find_relaxed_material( - api_client: APIClient, - material, - owner_id: str, - kgrid: Union[List[int], str, None] = ANY_KGRID, - unit_name: str = "pw_scf", + api_client: APIClient, material, owner_id: str, kgrid: Optional[List[int]] = None, unit_name: str = "pw_scf" ) -> Optional[Material]: """ Finds a relaxed version of a material: the final structure of a finished job on a material @@ -114,8 +110,7 @@ def find_relaxed_material( api_client (APIClient): API client instance carrying the authorization context. material: mat3ra-made Material object (must have a .hash property). owner_id (str): Account ID under which to search. - kgrid (List[int], optional): K-grid the relaxation ran on, None for the platform default, - see `get_kgrid_query`. + kgrid (List[int], optional): Exact k-grid dimensions the relaxation ran on; None for no condition. unit_name (str): Name of the relaxation unit, e.g. "pw_vc-relax". Returns: diff --git a/tests/py/unit/core/entity/test_job_api.py b/tests/py/unit/core/entity/test_job_api.py index e78490073..dc3a87934 100644 --- a/tests/py/unit/core/entity/test_job_api.py +++ b/tests/py/unit/core/entity/test_job_api.py @@ -3,10 +3,10 @@ import pytest from mat3ra.notebooks_utils.core.entity.job.api import ( - ANY_KGRID, create_job, find_job_for_material, find_job_for_material_with_property, + get_kgrid_query, ) OWNER_ID = "account-1" @@ -178,40 +178,44 @@ def test_find_job_for_material_with_property_returns_none_when_no_job_reported_i assert "tags" not in client.jobs.list.call_args.args[0] -KGRID_QUERY: Dict[str, Any] = { +SCF_KGRID_QUERY: Dict[str, Any] = { "workflow.subworkflows.units": { "$elemMatch": {"name": "pw_scf", "context": {"$elemMatch": {"name": "kgrid", "data.dimensions": [4, 4, 4]}}} } } -DEFAULT_KGRID_QUERY: Dict[str, Any] = { - "workflow.subworkflows.units": {"$elemMatch": {"name": "pw_scf", "context.name": {"$ne": "kgrid"}}} +RELAX_KGRID_QUERY: Dict[str, Any] = { + "workflow.subworkflows.units": { + "$elemMatch": {"name": "pw_relax", "context": {"$elemMatch": {"name": "kgrid", "data.dimensions": [4, 4, 4]}}} + } } @pytest.mark.parametrize( - ("kgrid", "expected_kgrid_query"), - [(ANY_KGRID, {}), (None, DEFAULT_KGRID_QUERY), ([4, 4, 4], KGRID_QUERY)], + ("kgrid", "unit_name", "expected_query"), + [(None, "pw_relax", {}), ([4, 4, 4], "pw_scf", SCF_KGRID_QUERY), ([4, 4, 4], "pw_relax", RELAX_KGRID_QUERY)], ) -def test_find_job_for_material_with_property_matches_the_pw_scf_kgrid(kgrid, expected_kgrid_query): +def test_get_kgrid_query(kgrid, unit_name, expected_query): + assert get_kgrid_query(kgrid, unit_name) == expected_query + + +def test_find_job_for_material_with_property_matches_the_pw_scf_kgrid(): client = MagicMock() client.jobs.list.return_value = [EXISTING_JOB] client.properties.get_for_job.return_value = [{"name": PROPERTY_NAME}] - job = find_job_for_material_with_property(client, MATERIAL_INITIAL["_id"], PROPERTY_NAME, OWNER_ID, kgrid=kgrid) + job = find_job_for_material_with_property(client, MATERIAL_INITIAL["_id"], PROPERTY_NAME, OWNER_ID, kgrid=[4, 4, 4]) assert job == EXISTING_JOB - client.jobs.list.assert_called_once_with( - {"_material._id": MATERIAL_INITIAL["_id"], "owner._id": OWNER_ID, "status": "finished", **expected_kgrid_query} - ) + assert SCF_KGRID_QUERY.items() <= client.jobs.list.call_args.args[0].items() -@pytest.mark.parametrize(("kgrid", "expected_kgrid_query"), [(None, DEFAULT_KGRID_QUERY), ([4, 4, 4], KGRID_QUERY)]) -def test_find_job_for_material_matches_the_kgrid(kgrid, expected_kgrid_query): +def test_find_job_for_material_matches_the_kgrid_of_the_unit(): client = MagicMock() client.jobs.list.return_value = [EXISTING_JOB] - job = find_job_for_material(client, MATERIAL_INITIAL["_id"], RELAX_WORKFLOW_NAME, OWNER_ID, kgrid=kgrid) + job = find_job_for_material( + client, MATERIAL_INITIAL["_id"], RELAX_WORKFLOW_NAME, OWNER_ID, kgrid=[4, 4, 4], unit_name="pw_relax" + ) assert job == EXISTING_JOB - kgrid_key = "workflow.subworkflows.units" - assert client.jobs.list.call_args.args[0][kgrid_key] == expected_kgrid_query[kgrid_key] + assert RELAX_KGRID_QUERY.items() <= client.jobs.list.call_args.args[0].items() diff --git a/tests/py/unit/core/entity/test_material_api.py b/tests/py/unit/core/entity/test_material_api.py index 81d4b1b8e..31e607df9 100644 --- a/tests/py/unit/core/entity/test_material_api.py +++ b/tests/py/unit/core/entity/test_material_api.py @@ -6,6 +6,7 @@ import pytest from mat3ra.made.material import Material +from mat3ra.notebooks_utils.core.entity.job.api import get_kgrid_query from mat3ra.notebooks_utils.core.entity.material.api import ( find_material_set, find_relaxed_material, @@ -359,35 +360,18 @@ def test_find_relaxed_material_skips_a_final_structure_with_the_same_hash(): assert client.properties.get_for_job.call_count == 2 -@pytest.mark.parametrize( - ("kgrid", "expected_unit_condition"), - [ - (None, {"name": "pw_vc-relax", "context.name": {"$ne": "kgrid"}}), - ( - [4, 4, 4], - {"name": "pw_vc-relax", "context": {"$elemMatch": {"name": "kgrid", "data.dimensions": [4, 4, 4]}}}, - ), - ], -) -def test_find_relaxed_material_matches_the_relaxation_kgrid(kgrid, expected_unit_condition): +def test_find_relaxed_material_matches_the_relaxation_kgrid(): client = MagicMock() client.materials.list.return_value = [SAVED_DEFECTIVE] client.jobs.list.return_value = [FINISHED_JOB] client.properties.get_for_job.return_value = [{"materialId": "m-relaxed"}] client.materials.get.return_value = RELAXED_MATERIAL_DOC - relaxed = find_relaxed_material(client, DEFECTIVE_MATERIAL, OWNER_ID, kgrid=kgrid, unit_name="pw_vc-relax") + relaxed = find_relaxed_material(client, DEFECTIVE_MATERIAL, OWNER_ID, kgrid=[4, 4, 4], unit_name="pw_vc-relax") assert relaxed is not None assert relaxed.name == "B-vacancy h-BN relaxed" - client.jobs.list.assert_called_once_with( - { - "_material._id": {"$in": [SAVED_DEFECTIVE["_id"]]}, - "owner._id": OWNER_ID, - "status": "finished", - "workflow.subworkflows.units": {"$elemMatch": expected_unit_condition}, - } - ) + assert get_kgrid_query([4, 4, 4], "pw_vc-relax").items() <= client.jobs.list.call_args.args[0].items() @pytest.mark.parametrize( From c4408f6aa4f406f274e2bd98d0c799bc220d8d74 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 28 Sep 2026 19:27:40 -0700 Subject: [PATCH 10/18] fix(SOF-7975): take the neutral notebook's pristine total energy from a job on the defect job's k-grid The Standata workflow picks the pristine total_energy with the highest precision.value, and precision does not measure the k-grid: on production job u5EcXhCocETccDB6a a platform-default (Gamma-only) pristine job carried precision 2000 against 768 for the 4x4x4 one, so a 4x4x4 defect SCF was referenced to the Gamma-only energy (E_f -2.09 eV instead of +0.37). The notebook now resolves one grid, SCF_KGRID or the pristine's platform default (get_platform_default_kgrid, copied from the charged notebook), applies it to the defect SCF through apply_scf_kgrid in place of the inline PointsGridDataProvider block, finds the pristine's finished total_energy job on that grid with find_job_for_material_with_property, raises naming the material and grid when there is none, and names that job to the workflow's assign-reference-job-filter, as the charged notebook does. That lookup replaces the find_total_energy_for_material check, which was PRISTINE_TOTAL_ENERGY_SOURCE's only reader; the parameter described highest-precision selection by any owner, which the job filter no longer allows, so it goes. The header says the Total Energy job must be on the same k-grid. Same 42 cells in the same order as origin/main. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../workflows/defect_formation_energy.ipynb | 46 ++++++++++--------- 1 file changed, 25 insertions(+), 21 deletions(-) diff --git a/other/materials_designer/workflows/defect_formation_energy.ipynb b/other/materials_designer/workflows/defect_formation_energy.ipynb index 5891f40c0..0577a730d 100644 --- a/other/materials_designer/workflows/defect_formation_energy.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy.ipynb @@ -14,7 +14,7 @@ "- **[0] Defective supercell** — the structure containing the defect(s) (vacancy, substitution, interstitial, or a mix).\n", "- **[1] Pristine supercell** — the defect-free version of the *same* supercell.\n", "\n", - "Only the defective cell is computed here (`pw_scf`). The pristine reference energy is fetched from a previously-finished **Total Energy** job, so the pristine must already be converged/relaxed and have such a job on the platform. Elemental chemical potentials are taken from Standata elemental reference materials (each needs a total energy too), as in [Formation Energy](formation_energy.ipynb).\n", + "Only the defective cell is computed here (`pw_scf`). The pristine reference energy is fetched from a previously-finished **Total Energy** job on the same k-grid, so the pristine must already be converged/relaxed and have such a job on the platform. Elemental chemical potentials are taken from Standata elemental reference materials (each needs a total energy too), as in [Formation Energy](formation_energy.ipynb).\n", "\n", "Formula:\n", "\n", @@ -121,12 +121,7 @@ "outputs": [], "source": [ "# K-grid for the defective-cell SCF (if not set, KPPRA is used by default)\n", - "SCF_KGRID = None # e.g. [4, 4, 4]\n", - "\n", - "# Whose total_energy properties to consider for the pristine reference:\n", - "# \"public\" (any owner, highest precision wins), \"curators\" (only curators'),\n", - "# or \"my_account\" (curators' or your own).\n", - "PRISTINE_TOTAL_ENERGY_SOURCE = \"my_account\"\n" + "SCF_KGRID = None # e.g. [4, 4, 4]\n" ] }, { @@ -364,17 +359,32 @@ "metadata": {}, "outputs": [], "source": [ + "from mat3ra.made.reciprocal.lattice_reciprocal import ReciprocalLattice\n", + "from mat3ra.notebooks_utils.core.entity.job.api import find_job_for_material_with_property\n", "from mat3ra.notebooks_utils.core.entity.material.api import get_or_create_material\n", - "from mat3ra.notebooks_utils.core.entity.property.api import find_total_energy_for_material\n", + "\n", + "PLATFORM_DEFAULT_KPPRA = 10 # defaultKPPRA in web-app src/application/imports/app_settings/settings.ts\n", + "\n", + "\n", + "def get_platform_default_kgrid(material):\n", + " \"\"\"The k-grid the platform gives `material` when none is set, from PLATFORM_DEFAULT_KPPRA.\"\"\"\n", + " number_of_kpoints = PLATFORM_DEFAULT_KPPRA / len(material.basis.elements.ids)\n", + " return ReciprocalLattice(**material.lattice.to_dict()).get_dimensions_from_points_count(number_of_kpoints)\n", + "\n", "\n", "saved_defective = Material.create(get_or_create_material(client, defective_material, ACCOUNT_ID))\n", "saved_pristine = Material.create(get_or_create_material(client, pristine_material, ACCOUNT_ID))\n", "\n", - "pristine_te_property = find_total_energy_for_material(\n", - " client, saved_pristine.id, source=PRISTINE_TOTAL_ENERGY_SOURCE\n", + "kgrid = SCF_KGRID or get_platform_default_kgrid(saved_pristine)\n", + "pristine_reference_job = find_job_for_material_with_property(\n", + " client, saved_pristine.id, \"total_energy\", ACCOUNT_ID, kgrid=kgrid\n", ")\n", - "if pristine_te_property is None:\n", - " raise RuntimeError(\"Run total_energy.ipynb for the pristine material first.\")\n", + "if pristine_reference_job is None:\n", + " raise RuntimeError(\n", + " f\"No finished Total Energy job for '{saved_pristine.name}' on k-grid {kgrid}: \"\n", + " \"run Total Energy on the pristine at this k-grid.\"\n", + " )\n", + "print(f\"♻️ pristine reference: job {pristine_reference_job['_id']}, k-grid {kgrid}\")\n", "\n", "# Order matters: [0] defective (computed), [1] pristine (reference).\n", "materials = [saved_defective, saved_pristine]\n" @@ -421,8 +431,8 @@ "source": [ "from mat3ra.standata.workflows import WorkflowStandata\n", "from mat3ra.wode.workflows import Workflow\n", - "from mat3ra.wode.context.providers import PointsGridDataProvider\n", "from mat3ra.notebooks_utils.ipython.entity.workflow.visualize import visualize_workflow\n", + "from mat3ra.notebooks_utils.workflow import apply_scf_kgrid, set_assignment_value\n", "\n", "defect_workflow_config = WorkflowStandata.filter_by_application(app.name).get_by_name_first_match(\n", " WORKFLOW_SEARCH_TERM\n", @@ -432,14 +442,8 @@ "print(f\"Loaded workflow: {defect_workflow.name}\")\n", "print(f\"Multi-material: {getattr(defect_workflow, 'isMultiMaterial', False)}\")\n", "\n", - "# K-grid for the defective-cell SCF.\n", - "if SCF_KGRID is not None:\n", - " new_context = PointsGridDataProvider(material=defective_material, dimensions=SCF_KGRID, isEdited=True).get_context_item_data()\n", - " for subworkflow in defect_workflow.subworkflows:\n", - " unit = subworkflow.get_unit_by_name(name=\"pw_scf\")\n", - " if unit:\n", - " unit.add_context(new_context)\n", - " subworkflow.set_unit(unit)\n", + "apply_scf_kgrid(defect_workflow, kgrid, material=defective_material)\n", + "set_assignment_value(defect_workflow, \"assign-reference-job-filter\", str({\"$in\": [pristine_reference_job[\"_id\"]]}))\n", "\n", "visualize_workflow(defect_workflow)" ] From f1cec3445be523d6b816097e84cacb9e4c0cb6f5 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 28 Sep 2026 20:02:14 -0700 Subject: [PATCH 11/18] fix(SOF-7975): set SCF_KGRID explicitly in both defect formation energy notebooks, with no platform-default grid VB's decision: the k-grid is never left to the platform. SCF_KGRID is one explicit grid, the same for the pristine references and the defective cells (divided by the supercell scaling in the charged notebook), because a KPPRA default resolved on a vacancy cell is different from the pristine's by definition. Both notebooks now default to SCF_KGRID = [4, 4, 4] and drop PLATFORM_DEFAULT_KPPRA, get_platform_default_kgrid, the ReciprocalLattice import and every None branch; get_scf_kgrid_for_supercell keeps only the division, and sections 3.3 and 4.2 of the charged notebook name SCF_KGRID instead of the platform default. The neutral notebook's header shows the general formula again, as on main, because it is what the Standata workflow computes: the job reports it at E_F = 0 and this notebook runs q = 0. No cells are added or removed. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../workflows/defect_formation_energy.ipynb | 28 ++++++------------- .../defect_formation_energy_charged.ipynb | 26 +++++------------ 2 files changed, 16 insertions(+), 38 deletions(-) diff --git a/other/materials_designer/workflows/defect_formation_energy.ipynb b/other/materials_designer/workflows/defect_formation_energy.ipynb index 0577a730d..b6e184b74 100644 --- a/other/materials_designer/workflows/defect_formation_energy.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy.ipynb @@ -7,7 +7,7 @@ "source": [ "# Defect Formation Energy\n", "\n", - "Calculate the neutral defect formation energy (eV) of a defective supercell using a multi-material DFT workflow on the Mat3ra platform. For charged defects, supercell-size series and formation energy vs Fermi level, use [Defect Formation Energy of Charged Defects](defect_formation_energy_charged.ipynb).\n", + "Calculate the neutral defect formation energy (eV) of a defective supercell using a multi-material DFT workflow on the Mat3ra platform.\n", "\n", "The workflow takes **two materials, in order**:\n", "\n", @@ -18,7 +18,9 @@ "\n", "Formula:\n", "\n", - "$$E_{\\text{defect}} = E_{\\text{defective}} - E_{\\text{pristine}} - \\sum_i \\Delta N_i\\, \\mu_i \\quad [\\text{eV}]$$\n", + "$$E_{\\text{defect}} = E_{\\text{defective}}[q] - E_{\\text{pristine}} - \\sum_i \\Delta N_i\\, \\mu_i + q\\,(E_{\\text{VBM}} + E_F) \\quad [\\text{eV}]$$\n", + "\n", + "The job reports this value at $E_F = 0$ (Fermi level at the VBM); the charged notebook adds $q\\,E_F$ when plotting against the Fermi level. Here $q = 0$. For charged defects, supercell-size series and formation energy vs Fermi level, use [Defect Formation Energy of Charged Defects](defect_formation_energy_charged.ipynb).\n", "\n", "where $\\Delta N_i = \\text{count}_i(\\text{defective}) - \\text{count}_i(\\text{pristine})$ is the per-species atom-count change and $\\mu_i = E_{\\text{elemental},i} / n_{\\text{atoms},i}$.\n", "\n", @@ -120,8 +122,7 @@ "metadata": {}, "outputs": [], "source": [ - "# K-grid for the defective-cell SCF (if not set, KPPRA is used by default)\n", - "SCF_KGRID = None # e.g. [4, 4, 4]\n" + "SCF_KGRID = [4, 4, 4] # same grid for the pristine and defective cells\n" ] }, { @@ -359,32 +360,21 @@ "metadata": {}, "outputs": [], "source": [ - "from mat3ra.made.reciprocal.lattice_reciprocal import ReciprocalLattice\n", "from mat3ra.notebooks_utils.core.entity.job.api import find_job_for_material_with_property\n", "from mat3ra.notebooks_utils.core.entity.material.api import get_or_create_material\n", "\n", - "PLATFORM_DEFAULT_KPPRA = 10 # defaultKPPRA in web-app src/application/imports/app_settings/settings.ts\n", - "\n", - "\n", - "def get_platform_default_kgrid(material):\n", - " \"\"\"The k-grid the platform gives `material` when none is set, from PLATFORM_DEFAULT_KPPRA.\"\"\"\n", - " number_of_kpoints = PLATFORM_DEFAULT_KPPRA / len(material.basis.elements.ids)\n", - " return ReciprocalLattice(**material.lattice.to_dict()).get_dimensions_from_points_count(number_of_kpoints)\n", - "\n", - "\n", "saved_defective = Material.create(get_or_create_material(client, defective_material, ACCOUNT_ID))\n", "saved_pristine = Material.create(get_or_create_material(client, pristine_material, ACCOUNT_ID))\n", "\n", - "kgrid = SCF_KGRID or get_platform_default_kgrid(saved_pristine)\n", "pristine_reference_job = find_job_for_material_with_property(\n", - " client, saved_pristine.id, \"total_energy\", ACCOUNT_ID, kgrid=kgrid\n", + " client, saved_pristine.id, \"total_energy\", ACCOUNT_ID, kgrid=SCF_KGRID\n", ")\n", "if pristine_reference_job is None:\n", " raise RuntimeError(\n", - " f\"No finished Total Energy job for '{saved_pristine.name}' on k-grid {kgrid}: \"\n", + " f\"No finished Total Energy job for '{saved_pristine.name}' on k-grid {SCF_KGRID}: \"\n", " \"run Total Energy on the pristine at this k-grid.\"\n", " )\n", - "print(f\"♻️ pristine reference: job {pristine_reference_job['_id']}, k-grid {kgrid}\")\n", + "print(f\"♻️ pristine reference: job {pristine_reference_job['_id']}, k-grid {SCF_KGRID}\")\n", "\n", "# Order matters: [0] defective (computed), [1] pristine (reference).\n", "materials = [saved_defective, saved_pristine]\n" @@ -442,7 +432,7 @@ "print(f\"Loaded workflow: {defect_workflow.name}\")\n", "print(f\"Multi-material: {getattr(defect_workflow, 'isMultiMaterial', False)}\")\n", "\n", - "apply_scf_kgrid(defect_workflow, kgrid, material=defective_material)\n", + "apply_scf_kgrid(defect_workflow, SCF_KGRID, material=defective_material)\n", "set_assignment_value(defect_workflow, \"assign-reference-job-filter\", str({\"$in\": [pristine_reference_job[\"_id\"]]}))\n", "\n", "visualize_workflow(defect_workflow)" diff --git a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb index 2bb593f15..5eed0e213 100644 --- a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb @@ -132,7 +132,7 @@ "\n", "CHARGES = [0] # e.g. [1, 0, -1, -2, -3]\n", "\n", - "SCF_KGRID = None # e.g. [8, 8, 8]\n", + "SCF_KGRID = [4, 4, 4] # same grid for the pristine and defective cells, divided by the supercell scaling\n", "\n", "RELAX_PRISTINE_MATERIAL = False\n", "RELAX_DEFECTIVE_MATERIAL = False" @@ -287,7 +287,7 @@ "metadata": {}, "source": [ "### 3.3. Define how jobs are created\n", - "Every job is tagged with its charge; the pristine references run neutral, tagged `charge:0`. A reference is reused when a finished `charge:0` job on that pristine cell already reported the property (`total_energy`, `band_gaps`), on the same k-grid (the platform default when `SCF_KGRID` is not set). The Defect Formation Energy workflow reads each reference from the one job the notebook names to it, so no other calculation on the same cell can stand in; the relaxation carries no charge tag and is never taken for one." + "Every job is tagged with its charge; the pristine references run neutral, tagged `charge:0`. A reference is reused when a finished `charge:0` job on that pristine cell already reported the property (`total_energy`, `band_gaps`), on the same k-grid, `SCF_KGRID` divided by the supercell scaling. The Defect Formation Energy workflow reads each reference from the one job the notebook names to it, so no other calculation on the same cell can stand in; the relaxation carries no charge tag and is never taken for one." ] }, { @@ -297,12 +297,9 @@ "metadata": {}, "outputs": [], "source": [ - "from mat3ra.made.reciprocal.lattice_reciprocal import ReciprocalLattice\n", "from mat3ra.notebooks_utils.api.job import submit_jobs, wait_for_jobs_to_finish_async\n", "from mat3ra.notebooks_utils.job import create_job\n", "\n", - "PLATFORM_DEFAULT_KPPRA = 10 # defaultKPPRA in web-app src/application/imports/app_settings/settings.ts\n", - "\n", "\n", "def create_job_for_materials(materials, workflow, tags=None):\n", " job = create_job(\n", @@ -315,13 +312,7 @@ " compute=compute.to_dict(),\n", " tags=tags,\n", " )\n", - " return job[0] if isinstance(job, list) else job\n", - "\n", - "\n", - "def get_platform_default_kgrid(material):\n", - " \"\"\"The k-grid the platform gives `material` when none is set, from PLATFORM_DEFAULT_KPPRA.\"\"\"\n", - " number_of_kpoints = PLATFORM_DEFAULT_KPPRA / len(material.basis.elements.ids)\n", - " return ReciprocalLattice(**material.lattice.to_dict()).get_dimensions_from_points_count(number_of_kpoints)" + " return job[0] if isinstance(job, list) else job" ] }, { @@ -362,7 +353,7 @@ "metadata": {}, "source": [ "### 4.2. Relax the pristine cell (optional)\n", - "With `RELAX_PRISTINE_MATERIAL`, the lattice and the atomic positions of the pristine cell are relaxed with the Variable-cell Relaxation workflow, or taken from an earlier finished variable-cell relaxation of the same structure, on the same k-grid (the platform default when `SCF_KGRID` is not set), and the supercells below are built from the relaxed cell." + "With `RELAX_PRISTINE_MATERIAL`, the lattice and the atomic positions of the pristine cell are relaxed with the Variable-cell Relaxation workflow, or taken from an earlier finished variable-cell relaxation of the same structure on the k-grid `SCF_KGRID`, and the supercells below are built from the relaxed cell." ] }, { @@ -383,9 +374,8 @@ "\n", "if RELAX_PRISTINE_MATERIAL:\n", " saved_pristine_material = Material.create(get_or_create_material(client, pristine_material, ACCOUNT_ID))\n", - " kgrid = SCF_KGRID or get_platform_default_kgrid(saved_pristine_material)\n", " relaxed_pristine_material = find_relaxed_material(\n", - " client, saved_pristine_material, ACCOUNT_ID, kgrid=kgrid, unit_name=\"pw_vc-relax\"\n", + " client, saved_pristine_material, ACCOUNT_ID, kgrid=SCF_KGRID, unit_name=\"pw_vc-relax\"\n", " )\n", " if relaxed_pristine_material is None:\n", " relax_workflow_config = WorkflowStandata.filter_by_application(APPLICATION_NAME).get_by_name_first_match(\n", @@ -393,7 +383,7 @@ " )\n", " relax_workflow = Workflow.create(relax_workflow_config)\n", " relax_workflow.name = f\"{relax_workflow.name} {saved_pristine_material.name}\"\n", - " apply_scf_kgrid(relax_workflow, kgrid, material=saved_pristine_material, unit_name=\"pw_vc-relax\")\n", + " apply_scf_kgrid(relax_workflow, SCF_KGRID, material=saved_pristine_material, unit_name=\"pw_vc-relax\")\n", " relax_job = create_job_for_materials([saved_pristine_material], relax_workflow)\n", " submit_jobs(client.jobs, [relax_job[\"_id\"]])\n", " print(f\"✅ created a relaxation job {relax_job['_id']}\")\n", @@ -566,9 +556,7 @@ "\n", "\n", "def get_scf_kgrid_for_supercell(scaling):\n", - " \"\"\"SCF_KGRID divided by the scaling, rounded, never below 1; without SCF_KGRID, the pristine supercell's default.\"\"\"\n", - " if SCF_KGRID is None:\n", - " return get_platform_default_kgrid(saved_supercell_pairs[scaling][0])\n", + " \"\"\"SCF_KGRID divided by the scaling, rounded, never below 1.\"\"\"\n", " return [max(1, round(dimension / scaling)) for dimension in SCF_KGRID]\n", "\n", "\n", From f1ab0f9ef2deded04f20093de66d6dc5e2c2b6e2 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Tue, 29 Sep 2026 10:16:34 -0700 Subject: [PATCH 12/18] feat(SOF-7975): chemical-potential references in both defect formation energy notebooks CHEMICAL_POTENTIAL_REFERENCES names one material per element that the host is in equilibrium with. solve_chemical_potentials turns their total energies into mu_i, one equation E = sum_i n_i mu_i per material, and raises ValueError when the materials' elements differ from the keys or the system is singular. The results cell then shows E_f = E_f^elem + sum_i dN_i (E_i - mu_i), taking sum_i dN_i E_i and dN_i from the defect job's own scope (SUM_DELTA_N_TIMES_MU, DELTA_N_BY_SYMBOL). A reference material's total energy comes from its Total Energy job on SCF_KGRID, else from its highest-precision one. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../workflows/defect_formation_energy.ipynb | 26 +++++++++++++- .../defect_formation_energy_charged.ipynb | 26 +++++++++++++- .../entity/property/chemical_potentials.py | 22 ++++++++++++ .../core/entity/test_chemical_potentials.py | 35 +++++++++++++++++++ 4 files changed, 107 insertions(+), 2 deletions(-) create mode 100644 src/py/mat3ra/notebooks_utils/core/entity/property/chemical_potentials.py create mode 100644 tests/py/unit/core/entity/test_chemical_potentials.py diff --git a/other/materials_designer/workflows/defect_formation_energy.ipynb b/other/materials_designer/workflows/defect_formation_energy.ipynb index b6e184b74..316112ff4 100644 --- a/other/materials_designer/workflows/defect_formation_energy.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy.ipynb @@ -122,6 +122,7 @@ "metadata": {}, "outputs": [], "source": [ + "CHEMICAL_POTENTIAL_REFERENCES = None # {\"O\": \"\", ...}: one reference material per element; None = elemental\n", "SCF_KGRID = [4, 4, 4] # same grid for the pristine and defective cells\n" ] }, @@ -557,7 +558,8 @@ "metadata": {}, "source": [ "## 7. Retrieve results\n", - "### 7.1. Retrieve and visualize defect formation energy" + "### 7.1. Retrieve and visualize defect formation energy\n", + "With `CHEMICAL_POTENTIAL_REFERENCES`, one material per element that the host is in equilibrium with, each giving $E = \\sum_i n_i\\, \\mu_i$, the formation energy is also printed at those references, $E_f = E_f^{\\text{elem}} + \\sum_i \\Delta N_i\\,(E_i - \\mu_i)$, with $E_i$ the elemental energy per atom the job used. Each material's total energy comes from its Total Energy job on `SCF_KGRID`, else from its highest-precision one." ] }, { @@ -567,8 +569,30 @@ "metadata": {}, "outputs": [], "source": [ + "from collections import Counter\n", + "from mat3ra.notebooks_utils.core.entity.property.api import find_total_energy_for_material\n", + "from mat3ra.notebooks_utils.core.entity.property.chemical_potentials import solve_chemical_potentials\n", "from mat3ra.notebooks_utils.ipython.entity.property.visualize import visualize_properties\n", "\n", + "if CHEMICAL_POTENTIAL_REFERENCES:\n", + " reference_energies = {}\n", + " for element, name in CHEMICAL_POTENTIAL_REFERENCES.items():\n", + " material = load_material(client, FOLDER, name, ACCOUNT_ID)\n", + " material_id = client.materials.list({\"hash\": material.hash, \"owner._id\": ACCOUNT_ID})[0][\"_id\"]\n", + " job = find_job_for_material_with_property(client, material_id, \"total_energy\", ACCOUNT_ID, kgrid=SCF_KGRID)\n", + " energy = (client.properties.get_for_job(job[\"_id\"], property_name=\"total_energy\")[0][\"value\"] if job\n", + " else find_total_energy_for_material(client, material_id, source=\"public\")[\"data\"][\"value\"])\n", + " reference_energies[element] = (dict(Counter(material.basis.elements.values)), energy)\n", + " chemical_potentials = solve_chemical_potentials(reference_energies)\n", + " print(f\"Chemical potentials (eV/atom): {chemical_potentials}\")\n", + " job_scope = client.jobs.get(defect_job_id)[\"scopeTrack\"]\n", + " scope = {key: value for item in job_scope for key, value in item[\"scope\"][\"global\"].items()}\n", + " formation_energy_at_references = scope[\"DEFECT_FORMATION_ENERGY\"] + scope[\"SUM_DELTA_N_TIMES_MU\"] - sum(\n", + " delta_n * chemical_potentials[element] for element, delta_n in scope[\"DELTA_N_BY_SYMBOL\"].items()\n", + " )\n", + " print(f\"Defect formation energy: {scope['DEFECT_FORMATION_ENERGY']:.4f} eV with elemental chemical potentials, \"\n", + " f\"{formation_energy_at_references:.4f} eV at the references\")\n", + "\n", "defect_energy_data = client.properties.get_for_job(defect_job_id)\n", "visualize_properties(defect_energy_data, title=\"Defect Formation Energy\")" ] diff --git a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb index 5eed0e213..fb08c69ab 100644 --- a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb @@ -133,6 +133,7 @@ "CHARGES = [0] # e.g. [1, 0, -1, -2, -3]\n", "\n", "SCF_KGRID = [4, 4, 4] # same grid for the pristine and defective cells, divided by the supercell scaling\n", + "CHEMICAL_POTENTIAL_REFERENCES = None # {\"O\": \"\", ...}: one reference material per element; None = elemental\n", "\n", "RELAX_PRISTINE_MATERIAL = False\n", "RELAX_DEFECTIVE_MATERIAL = False" @@ -750,7 +751,7 @@ "source": [ "## 7. Retrieve results\n", "### 7.1. Formation energies\n", - "`formation_energy` is $E_f[X^q]$ at $\\mu_e = 0$, as the workflow reports it. `length` is $L = V^{1/3}$ of the supercell, the length scale the image-charge terms in section 7.3 are expressed in." + "`formation_energy` is $E_f[X^q]$ at $\\mu_e = 0$, as the workflow reports it. `length` is $L = V^{1/3}$ of the supercell, the length scale the image-charge terms in section 7.3 are expressed in. With `CHEMICAL_POTENTIAL_REFERENCES`, one material per element that the host is in equilibrium with, each giving $E = \\sum_i n_i\\, \\mu_i$, `formation_energy_at_references` is $E_f = E_f^{\\text{elem}} + \\sum_i \\Delta N_i\\,(E_i - \\mu_i)$, with $E_i$ the elemental energy per atom the job used. Each material's total energy comes from its Total Energy job on `SCF_KGRID`, else from its highest-precision one." ] }, { @@ -760,6 +761,23 @@ "metadata": {}, "outputs": [], "source": [ + "from collections import Counter\n", + "from mat3ra.notebooks_utils.core.entity.property.api import find_total_energy_for_material\n", + "from mat3ra.notebooks_utils.core.entity.property.chemical_potentials import solve_chemical_potentials\n", + "\n", + "if CHEMICAL_POTENTIAL_REFERENCES:\n", + " reference_energies = {}\n", + " for element, name in CHEMICAL_POTENTIAL_REFERENCES.items():\n", + " matches = [Material.create(data) for data in Materials.get_by_name(name) if data[\"name\"] == name]\n", + " material = matches[0] if matches else load_material(client, FOLDER, name, ACCOUNT_ID)\n", + " material_id = client.materials.list({\"hash\": material.hash, \"owner._id\": ACCOUNT_ID})[0][\"_id\"]\n", + " job = find_job_for_material_with_property(client, material_id, \"total_energy\", ACCOUNT_ID, kgrid=SCF_KGRID)\n", + " energy = (client.properties.get_for_job(job[\"_id\"], property_name=\"total_energy\")[0][\"value\"] if job\n", + " else find_total_energy_for_material(client, material_id, source=\"public\")[\"data\"][\"value\"])\n", + " reference_energies[element] = (dict(Counter(material.basis.elements.values)), energy)\n", + " chemical_potentials = solve_chemical_potentials(reference_energies)\n", + " print(f\"Chemical potentials (eV/atom): {chemical_potentials}\")\n", + "\n", "results_records = []\n", "for record in job_records:\n", " job = client.jobs.get(record[\"job_id\"])\n", @@ -769,6 +787,12 @@ " \"final_status\": job.get(\"status\"),\n", " \"formation_energy\": properties[0].get(\"value\") if properties else None,\n", " })\n", + " if CHEMICAL_POTENTIAL_REFERENCES and properties:\n", + " scope = {key: value for item in job[\"scopeTrack\"] for key, value in item[\"scope\"][\"global\"].items()}\n", + " results_records[-1][\"formation_energy_at_references\"] = (\n", + " properties[0][\"value\"] + scope[\"SUM_DELTA_N_TIMES_MU\"]\n", + " - sum(delta_n * chemical_potentials[element] for element, delta_n in scope[\"DELTA_N_BY_SYMBOL\"].items())\n", + " )\n", "\n", "results_df = pd.DataFrame(results_records)\n", "successful_df = results_df[(results_df[\"final_status\"] == \"finished\") & results_df[\"formation_energy\"].notna()]\n", diff --git a/src/py/mat3ra/notebooks_utils/core/entity/property/chemical_potentials.py b/src/py/mat3ra/notebooks_utils/core/entity/property/chemical_potentials.py new file mode 100644 index 000000000..b7755de41 --- /dev/null +++ b/src/py/mat3ra/notebooks_utils/core/entity/property/chemical_potentials.py @@ -0,0 +1,22 @@ +from typing import Dict, Tuple + +import numpy as np + + +def solve_chemical_potentials(reference_energies: Dict[str, Tuple[Dict[str, int], float]]) -> Dict[str, float]: + """ + Chemical potential of each element (eV/atom) from one reference material per element, given as its atom counts + and its total energy (eV) for them: each material gives E = sum_i n_i * mu_i, and the square system is solved. + + Raises: + ValueError: If the materials contain other elements than the keys, or do not determine every potential. + """ + elements = sorted(reference_energies) + compositions = [composition for composition, _ in reference_energies.values()] + if set().union(*compositions) != set(elements): + raise ValueError(f"The reference materials must contain exactly the elements {elements}.") + matrix = np.array([[composition.get(element, 0) for element in elements] for composition in compositions]) + if np.linalg.matrix_rank(matrix) < len(elements): + raise ValueError("The reference materials do not determine every chemical potential.") + energies = [energy for _, energy in reference_energies.values()] + return dict(zip(elements, np.linalg.solve(matrix, energies).tolist())) diff --git a/tests/py/unit/core/entity/test_chemical_potentials.py b/tests/py/unit/core/entity/test_chemical_potentials.py new file mode 100644 index 000000000..123ca5820 --- /dev/null +++ b/tests/py/unit/core/entity/test_chemical_potentials.py @@ -0,0 +1,35 @@ +"""Unit tests for chemical potentials from reference materials.""" + +import pytest +from mat3ra.notebooks_utils.core.entity.property.chemical_potentials import solve_chemical_potentials + +# Total energies (eV) as (atom counts, energy) of the Standata materials, PBE; ZrO2 is illustrative. +O2 = ({"O": 8}, -3499.6187) +HF = ({"Hf": 2}, -4320.4730) +ZR = ({"Zr": 2}, -2698.0925) +HFO2 = ({"Hf": 4, "O": 8}, -12183.3350) +HFO2_3X3X3 = ({"Hf": 108, "O": 216}, 27 * -12183.3350) +ZRO2 = ({"Zr": 4, "O": 8}, -9800.0) +MU_O_RICH = -3499.6187 / 8 +MU_HF_O_POOR = -4320.4730 / 2 + + +@pytest.mark.parametrize( + "reference_energies, expected", + [ + ({"O": O2, "Hf": HF, "Zr": ZR}, {"O": MU_O_RICH, "Hf": MU_HF_O_POOR, "Zr": -2698.0925 / 2}), + ( + {"O": O2, "Hf": HFO2, "Zr": ZRO2}, + {"O": MU_O_RICH, "Hf": -12183.3350 / 4 - 2 * MU_O_RICH, "Zr": -9800.0 / 4 - 2 * MU_O_RICH}, + ), + ({"Hf": HF, "O": HFO2}, {"Hf": MU_HF_O_POOR, "O": (-12183.3350 / 4 - MU_HF_O_POOR) / 2}), + ], +) +def test_solve_chemical_potentials(reference_energies, expected): + assert solve_chemical_potentials(reference_energies) == pytest.approx(expected) + + +@pytest.mark.parametrize("reference_energies", [{"Hf": HFO2, "O": HFO2}, {"Hf": HFO2, "O": HFO2_3X3X3}, {"Hf": HFO2}]) +def test_solve_chemical_potentials_raises(reference_energies): + with pytest.raises(ValueError): + solve_chemical_potentials(reference_energies) From b6e88661cd0ade0884da9c3878497e9c683ffb24 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Tue, 29 Sep 2026 11:11:18 -0700 Subject: [PATCH 13/18] fix(SOF-7975): resolve chemical-potential references before the defect job, from this account's Total Energy jobs get_reference_energies takes the platform materials the notebook loaded by name and finds each one's finished Total Energy job in this account: the one on SCF_KGRID, else the one other k-grid it was run on; several other grids or none raise, naming the jobs. The pristine takes the energy of its reference job, the one the workflow reads E_PRISTINE from, so the chemical potentials add up to the pristine energy; the charged notebook solves them per supercell size. Each reference prints its job, k-grid and energy. The highest-precision public fallback is gone. get_formation_energy_at_references reads DEFECT_FORMATION_ENERGY, SUM_DELTA_N_TIMES_MU and DELTA_N_BY_SYMBOL from the job's flattened scopeTrack (flatten_scope_track) and is tested on the Zr_Hf job's numbers: 0.3664 -> 0.0134 eV at the oxide references. get_kgrid_of_job reads a job's pw_scf k-grid where get_kgrid_query matches it. Both notebooks resolve the references before creating any defect job; the results cells only apply the shift. The charged 7.2 section states that the plots and the fit use elemental chemical potentials. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../workflows/defect_formation_energy.ipynb | 53 +++++++++-------- .../defect_formation_energy_charged.ipynb | 57 +++++++++++-------- .../notebooks_utils/core/entity/job/api.py | 7 +++ .../entity/property/chemical_potentials.py | 53 ++++++++++++++++- tests/py/unit/core/entity/test_job_api.py | 9 +++ ...y => test_property_chemical_potentials.py} | 33 +++++++++-- 6 files changed, 159 insertions(+), 53 deletions(-) rename tests/py/unit/core/entity/{test_chemical_potentials.py => test_property_chemical_potentials.py} (50%) diff --git a/other/materials_designer/workflows/defect_formation_energy.ipynb b/other/materials_designer/workflows/defect_formation_energy.ipynb index 316112ff4..fdd4ad9d2 100644 --- a/other/materials_designer/workflows/defect_formation_energy.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy.ipynb @@ -122,8 +122,8 @@ "metadata": {}, "outputs": [], "source": [ - "CHEMICAL_POTENTIAL_REFERENCES = None # {\"O\": \"\", ...}: one reference material per element; None = elemental\n", - "SCF_KGRID = [4, 4, 4] # same grid for the pristine and defective cells\n" + "SCF_KGRID = [4, 4, 4] # same grid for the pristine and defective cells\n", + "CHEMICAL_POTENTIAL_REFERENCES = None # {\"O\": \"\", ...}: one reference material per element; None = elemental\n" ] }, { @@ -363,6 +363,12 @@ "source": [ "from mat3ra.notebooks_utils.core.entity.job.api import find_job_for_material_with_property\n", "from mat3ra.notebooks_utils.core.entity.material.api import get_or_create_material\n", + "from mat3ra.notebooks_utils.core.entity.property.chemical_potentials import (\n", + " flatten_scope_track,\n", + " get_formation_energy_at_references,\n", + " get_reference_energies,\n", + " solve_chemical_potentials,\n", + ")\n", "\n", "saved_defective = Material.create(get_or_create_material(client, defective_material, ACCOUNT_ID))\n", "saved_pristine = Material.create(get_or_create_material(client, pristine_material, ACCOUNT_ID))\n", @@ -378,7 +384,22 @@ "print(f\"♻️ pristine reference: job {pristine_reference_job['_id']}, k-grid {SCF_KGRID}\")\n", "\n", "# Order matters: [0] defective (computed), [1] pristine (reference).\n", - "materials = [saved_defective, saved_pristine]\n" + "materials = [saved_defective, saved_pristine]\n", + "\n", + "if CHEMICAL_POTENTIAL_REFERENCES:\n", + " reference_materials = {\n", + " element: Material.create(\n", + " get_or_create_material(client, load_material(client, FOLDER, name, ACCOUNT_ID), ACCOUNT_ID)\n", + " )\n", + " for element, name in CHEMICAL_POTENTIAL_REFERENCES.items()\n", + " }\n", + " pristine_energy = client.properties.get_for_job(pristine_reference_job[\"_id\"], \"total_energy\")[0][\"value\"]\n", + " reference_energies = get_reference_energies(\n", + " client, reference_materials, ACCOUNT_ID, SCF_KGRID, {saved_pristine.id: pristine_energy}\n", + " )\n", + " chemical_potentials = solve_chemical_potentials(reference_energies)\n", + " print(\"Chemical potentials (eV/atom):\",\n", + " {element: round(value, 3) for element, value in chemical_potentials.items()})\n" ] }, { @@ -559,7 +580,7 @@ "source": [ "## 7. Retrieve results\n", "### 7.1. Retrieve and visualize defect formation energy\n", - "With `CHEMICAL_POTENTIAL_REFERENCES`, one material per element that the host is in equilibrium with, each giving $E = \\sum_i n_i\\, \\mu_i$, the formation energy is also printed at those references, $E_f = E_f^{\\text{elem}} + \\sum_i \\Delta N_i\\,(E_i - \\mu_i)$, with $E_i$ the elemental energy per atom the job used. Each material's total energy comes from its Total Energy job on `SCF_KGRID`, else from its highest-precision one." + "With `CHEMICAL_POTENTIAL_REFERENCES`, one material per element that the host is in equilibrium with, each giving $E = \\sum_i n_i\\, \\mu_i$, the formation energy is printed at those references, $E_f = E_f^{\\text{elem}} + \\sum_i \\Delta N_i\\,(E_i - \\mu_i)$, with $E_i$ the elemental energy per atom the job used. The chemical potentials are solved in section 3.5, before any defect job is created. Each material's total energy comes from its finished Total Energy job in this account, on `SCF_KGRID` if there is one, else on the one other k-grid it was run on. A reference named `PRISTINE_NAME` is the pristine supercell and takes the energy of its Total Energy job, so the chemical potentials add up to the pristine energy." ] }, { @@ -569,29 +590,13 @@ "metadata": {}, "outputs": [], "source": [ - "from collections import Counter\n", - "from mat3ra.notebooks_utils.core.entity.property.api import find_total_energy_for_material\n", - "from mat3ra.notebooks_utils.core.entity.property.chemical_potentials import solve_chemical_potentials\n", "from mat3ra.notebooks_utils.ipython.entity.property.visualize import visualize_properties\n", "\n", "if CHEMICAL_POTENTIAL_REFERENCES:\n", - " reference_energies = {}\n", - " for element, name in CHEMICAL_POTENTIAL_REFERENCES.items():\n", - " material = load_material(client, FOLDER, name, ACCOUNT_ID)\n", - " material_id = client.materials.list({\"hash\": material.hash, \"owner._id\": ACCOUNT_ID})[0][\"_id\"]\n", - " job = find_job_for_material_with_property(client, material_id, \"total_energy\", ACCOUNT_ID, kgrid=SCF_KGRID)\n", - " energy = (client.properties.get_for_job(job[\"_id\"], property_name=\"total_energy\")[0][\"value\"] if job\n", - " else find_total_energy_for_material(client, material_id, source=\"public\")[\"data\"][\"value\"])\n", - " reference_energies[element] = (dict(Counter(material.basis.elements.values)), energy)\n", - " chemical_potentials = solve_chemical_potentials(reference_energies)\n", - " print(f\"Chemical potentials (eV/atom): {chemical_potentials}\")\n", - " job_scope = client.jobs.get(defect_job_id)[\"scopeTrack\"]\n", - " scope = {key: value for item in job_scope for key, value in item[\"scope\"][\"global\"].items()}\n", - " formation_energy_at_references = scope[\"DEFECT_FORMATION_ENERGY\"] + scope[\"SUM_DELTA_N_TIMES_MU\"] - sum(\n", - " delta_n * chemical_potentials[element] for element, delta_n in scope[\"DELTA_N_BY_SYMBOL\"].items()\n", - " )\n", - " print(f\"Defect formation energy: {scope['DEFECT_FORMATION_ENERGY']:.4f} eV with elemental chemical potentials, \"\n", - " f\"{formation_energy_at_references:.4f} eV at the references\")\n", + " scope = flatten_scope_track(client.jobs.get(defect_job_id)[\"scopeTrack\"])\n", + " references = \", \".join(f\"{element}: {name}\" for element, name in CHEMICAL_POTENTIAL_REFERENCES.items())\n", + " print(f\"Defect formation energy at the references ({references}): \"\n", + " f\"{get_formation_energy_at_references(scope, chemical_potentials):.4f} eV\")\n", "\n", "defect_energy_data = client.properties.get_for_job(defect_job_id)\n", "visualize_properties(defect_energy_data, title=\"Defect Formation Energy\")" diff --git a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb index fb08c69ab..2d1394ef4 100644 --- a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb @@ -649,8 +649,36 @@ "metadata": {}, "outputs": [], "source": [ + "from mat3ra.notebooks_utils.core.entity.property.chemical_potentials import (\n", + " flatten_scope_track,\n", + " get_formation_energy_at_references,\n", + " get_reference_energies,\n", + " solve_chemical_potentials,\n", + ")\n", + "\n", "if pristine_job_ids:\n", - " await wait_for_jobs_to_finish_async(client.jobs, pristine_job_ids, poll_interval=POLL_INTERVAL)" + " await wait_for_jobs_to_finish_async(client.jobs, pristine_job_ids, poll_interval=POLL_INTERVAL)\n", + "\n", + "if CHEMICAL_POTENTIAL_REFERENCES:\n", + " reference_materials = {}\n", + " for element, name in CHEMICAL_POTENTIAL_REFERENCES.items():\n", + " if name != PRISTINE_NAME:\n", + " matches = [Material.create(data) for data in Materials.get_by_name(name) if data[\"name\"] == name]\n", + " material = matches[0] if matches else load_material(client, FOLDER, name, ACCOUNT_ID)\n", + " reference_materials[element] = Material.create(get_or_create_material(client, material, ACCOUNT_ID))\n", + " chemical_potentials_by_scaling = {}\n", + " for scaling, (pristine_supercell, _) in saved_supercell_pairs.items():\n", + " pristine_job_id = pristine_reference_jobs[scaling][\"total_energy\"][\"_id\"]\n", + " pristine_energy = client.properties.get_for_job(pristine_job_id, \"total_energy\")[0][\"value\"]\n", + " materials_by_element = {\n", + " element: reference_materials.get(element, pristine_supercell) for element in CHEMICAL_POTENTIAL_REFERENCES\n", + " }\n", + " reference_energies = get_reference_energies(\n", + " client, materials_by_element, ACCOUNT_ID, SCF_KGRID, {pristine_supercell.id: pristine_energy}\n", + " )\n", + " chemical_potentials_by_scaling[scaling] = solve_chemical_potentials(reference_energies)\n", + " print(f\"n={scaling}: chemical potentials (eV/atom):\",\n", + " {element: round(value, 3) for element, value in chemical_potentials_by_scaling[scaling].items()})" ] }, { @@ -751,7 +779,7 @@ "source": [ "## 7. Retrieve results\n", "### 7.1. Formation energies\n", - "`formation_energy` is $E_f[X^q]$ at $\\mu_e = 0$, as the workflow reports it. `length` is $L = V^{1/3}$ of the supercell, the length scale the image-charge terms in section 7.3 are expressed in. With `CHEMICAL_POTENTIAL_REFERENCES`, one material per element that the host is in equilibrium with, each giving $E = \\sum_i n_i\\, \\mu_i$, `formation_energy_at_references` is $E_f = E_f^{\\text{elem}} + \\sum_i \\Delta N_i\\,(E_i - \\mu_i)$, with $E_i$ the elemental energy per atom the job used. Each material's total energy comes from its Total Energy job on `SCF_KGRID`, else from its highest-precision one." + "`formation_energy` is $E_f[X^q]$ at $\\mu_e = 0$, as the workflow reports it. `length` is $L = V^{1/3}$ of the supercell, the length scale the image-charge terms in section 7.3 are expressed in. With `CHEMICAL_POTENTIAL_REFERENCES`, one material per element that the host is in equilibrium with, each giving $E = \\sum_i n_i\\, \\mu_i$, `formation_energy_at_references` is $E_f = E_f^{\\text{elem}} + \\sum_i \\Delta N_i\\,(E_i - \\mu_i)$, with $E_i$ the elemental energy per atom the job used. The chemical potentials are solved in section 6.1 for each supercell size, before any defect job is created. Each material's total energy comes from its finished Total Energy job in this account, on `SCF_KGRID` if there is one, else on the one other k-grid it was run on. A reference named `PRISTINE_NAME` is the pristine supercell of that size and takes the energy of its Total Energy job, so the chemical potentials add up to the pristine energy." ] }, { @@ -761,23 +789,6 @@ "metadata": {}, "outputs": [], "source": [ - "from collections import Counter\n", - "from mat3ra.notebooks_utils.core.entity.property.api import find_total_energy_for_material\n", - "from mat3ra.notebooks_utils.core.entity.property.chemical_potentials import solve_chemical_potentials\n", - "\n", - "if CHEMICAL_POTENTIAL_REFERENCES:\n", - " reference_energies = {}\n", - " for element, name in CHEMICAL_POTENTIAL_REFERENCES.items():\n", - " matches = [Material.create(data) for data in Materials.get_by_name(name) if data[\"name\"] == name]\n", - " material = matches[0] if matches else load_material(client, FOLDER, name, ACCOUNT_ID)\n", - " material_id = client.materials.list({\"hash\": material.hash, \"owner._id\": ACCOUNT_ID})[0][\"_id\"]\n", - " job = find_job_for_material_with_property(client, material_id, \"total_energy\", ACCOUNT_ID, kgrid=SCF_KGRID)\n", - " energy = (client.properties.get_for_job(job[\"_id\"], property_name=\"total_energy\")[0][\"value\"] if job\n", - " else find_total_energy_for_material(client, material_id, source=\"public\")[\"data\"][\"value\"])\n", - " reference_energies[element] = (dict(Counter(material.basis.elements.values)), energy)\n", - " chemical_potentials = solve_chemical_potentials(reference_energies)\n", - " print(f\"Chemical potentials (eV/atom): {chemical_potentials}\")\n", - "\n", "results_records = []\n", "for record in job_records:\n", " job = client.jobs.get(record[\"job_id\"])\n", @@ -788,10 +799,8 @@ " \"formation_energy\": properties[0].get(\"value\") if properties else None,\n", " })\n", " if CHEMICAL_POTENTIAL_REFERENCES and properties:\n", - " scope = {key: value for item in job[\"scopeTrack\"] for key, value in item[\"scope\"][\"global\"].items()}\n", - " results_records[-1][\"formation_energy_at_references\"] = (\n", - " properties[0][\"value\"] + scope[\"SUM_DELTA_N_TIMES_MU\"]\n", - " - sum(delta_n * chemical_potentials[element] for element, delta_n in scope[\"DELTA_N_BY_SYMBOL\"].items())\n", + " results_records[-1][\"formation_energy_at_references\"] = get_formation_energy_at_references(\n", + " flatten_scope_track(job[\"scopeTrack\"]), chemical_potentials_by_scaling[record[\"scaling\"]]\n", " )\n", "\n", "results_df = pd.DataFrame(results_records)\n", @@ -808,7 +817,7 @@ "metadata": {}, "source": [ "### 7.2. Formation energy versus the electron chemical potential\n", - "Each charge state gives a straight line of slope $q$ over $\\mu_e \\in [0, E_{\\text{gap}}]$; the lower envelope is the charge state the defect actually adopts, and the crossings between the lines are the charge transition levels. The largest supercell with a result for every charge state is used, with $E_{\\text{gap}}$ from its own Band Gap job. The table lists each charge state's formation energy at the VBM and the range of the Fermi level $\\mu_e$ in which it is the stable state." + "Each charge state gives a straight line of slope $q$ over $\\mu_e \\in [0, E_{\\text{gap}}]$; the lower envelope is the charge state the defect actually adopts, and the crossings between the lines are the charge transition levels. The largest supercell with a result for every charge state is used, with $E_{\\text{gap}}$ from its own Band Gap job. The table lists each charge state's formation energy at the VBM and the range of the Fermi level $\\mu_e$ in which it is the stable state. The lines, and the fit in section 7.3, use the elemental $\\mu_i$ the workflow stores; `formation_energy_at_references` in section 7.1 is shifted by the same amount for every charge state of a size, so the transition levels do not change." ] }, { diff --git a/src/py/mat3ra/notebooks_utils/core/entity/job/api.py b/src/py/mat3ra/notebooks_utils/core/entity/job/api.py index 57c33dcd4..74d774a1e 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/job/api.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/job/api.py @@ -171,6 +171,13 @@ def get_kgrid_query(kgrid: Optional[List[int]], unit_name: str = "pw_scf") -> Di return {"workflow.subworkflows.units": {"$elemMatch": {"name": unit_name, "context": kgrid_context}}} +def get_kgrid_of_job(job: dict) -> Optional[List[int]]: + """The k-grid dimensions the job's `pw_scf` unit ran on, read where `get_kgrid_query` matches; None if unset.""" + units = [unit for subworkflow in job["workflow"]["subworkflows"] for unit in subworkflow["units"]] + contexts = [item for unit in units if unit["name"] == "pw_scf" for item in unit["context"]] + return next((item["data"]["dimensions"] for item in contexts if item["name"] == "kgrid"), None) + + def find_job_for_material_with_property( api_client: APIClient, material_id: str, diff --git a/src/py/mat3ra/notebooks_utils/core/entity/property/chemical_potentials.py b/src/py/mat3ra/notebooks_utils/core/entity/property/chemical_potentials.py index b7755de41..f8293904f 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/property/chemical_potentials.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/property/chemical_potentials.py @@ -1,6 +1,11 @@ -from typing import Dict, Tuple +from collections import Counter +from typing import Any, Dict, List, Tuple import numpy as np +from mat3ra.api_client import APIClient +from mat3ra.made.material import Material + +from ..job.api import get_kgrid_of_job def solve_chemical_potentials(reference_energies: Dict[str, Tuple[Dict[str, int], float]]) -> Dict[str, float]: @@ -20,3 +25,49 @@ def solve_chemical_potentials(reference_energies: Dict[str, Tuple[Dict[str, int] raise ValueError("The reference materials do not determine every chemical potential.") energies = [energy for _, energy in reference_energies.values()] return dict(zip(elements, np.linalg.solve(matrix, energies).tolist())) + + +def get_reference_energies( + client: APIClient, + materials_by_element: Dict[str, Material], + owner_id: str, + kgrid: List[int], + host_energy: Dict[str, float], +) -> Dict[str, Tuple[Dict[str, int], float]]: + """ + Atom counts and total energy of each platform material, as `solve_chemical_potentials` takes them. The energy is + `host_energy[material.id]` when given there, else that of the material's finished Total Energy job of `owner_id`: + the one on `kgrid`, else the only one, or one of several on the same other k-grid. + + Raises: + RuntimeError: If a material has no such job, or has them only on several other k-grids. + """ + reference_energies = {} + for element, material in materials_by_element.items(): + if material.id in host_energy: + energy, source = host_energy[material.id], "pristine reference" + else: + query = {"_material._id": material.id, "owner._id": owner_id, "status": "finished"} + jobs = [job for job in client.jobs.list(query) if client.properties.get_for_job(job["_id"], "total_energy")] + jobs = [job for job in jobs if get_kgrid_of_job(job) == kgrid] or jobs + if len({str(get_kgrid_of_job(job)) for job in jobs}) != 1: + found = ", ".join(f"job {job['_id']} on {get_kgrid_of_job(job)}" for job in jobs) or "none" + raise RuntimeError(f"Run Total Energy on '{material.name}' at k-grid {kgrid}; found: {found}.") + energy = client.properties.get_for_job(jobs[0]["_id"], property_name="total_energy")[0]["value"] + source = f"job {jobs[0]['_id']}, k-grid {get_kgrid_of_job(jobs[0])}" + print(f"{element}: {material.name}, {source}, E = {energy:.4f} eV") + reference_energies[element] = (dict(Counter(material.basis.elements.values)), energy) + return reference_energies + + +def flatten_scope_track(scope_track: List[dict]) -> Dict[str, Any]: + """The global scope of a job's `scopeTrack`, later values overriding earlier ones.""" + return {key: value for item in scope_track for key, value in item["scope"]["global"].items()} + + +def get_formation_energy_at_references(scope: Dict[str, Any], chemical_potentials: Dict[str, float]) -> float: + """Defect formation energy at `chemical_potentials` from its job's flattened scope, stored at elemental E_i.""" + delta_n_times_mu = sum( + count * chemical_potentials[element] for element, count in scope["DELTA_N_BY_SYMBOL"].items() + ) + return scope["DEFECT_FORMATION_ENERGY"] + scope["SUM_DELTA_N_TIMES_MU"] - delta_n_times_mu diff --git a/tests/py/unit/core/entity/test_job_api.py b/tests/py/unit/core/entity/test_job_api.py index dc3a87934..604a1d81d 100644 --- a/tests/py/unit/core/entity/test_job_api.py +++ b/tests/py/unit/core/entity/test_job_api.py @@ -6,6 +6,7 @@ create_job, find_job_for_material, find_job_for_material_with_property, + get_kgrid_of_job, get_kgrid_query, ) @@ -198,6 +199,14 @@ def test_get_kgrid_query(kgrid, unit_name, expected_query): assert get_kgrid_query(kgrid, unit_name) == expected_query +@pytest.mark.parametrize( + ("context", "expected_kgrid"), [([{"name": "kgrid", "data": {"dimensions": [4, 4, 4]}}], [4, 4, 4]), ([], None)] +) +def test_get_kgrid_of_job(context, expected_kgrid): + job = {"workflow": {"subworkflows": [{"units": [{"name": "pw_scf", "context": context}]}]}} + assert get_kgrid_of_job(job) == expected_kgrid + + def test_find_job_for_material_with_property_matches_the_pw_scf_kgrid(): client = MagicMock() client.jobs.list.return_value = [EXISTING_JOB] diff --git a/tests/py/unit/core/entity/test_chemical_potentials.py b/tests/py/unit/core/entity/test_property_chemical_potentials.py similarity index 50% rename from tests/py/unit/core/entity/test_chemical_potentials.py rename to tests/py/unit/core/entity/test_property_chemical_potentials.py index 123ca5820..e57a817de 100644 --- a/tests/py/unit/core/entity/test_chemical_potentials.py +++ b/tests/py/unit/core/entity/test_property_chemical_potentials.py @@ -1,15 +1,19 @@ """Unit tests for chemical potentials from reference materials.""" import pytest -from mat3ra.notebooks_utils.core.entity.property.chemical_potentials import solve_chemical_potentials +from mat3ra.notebooks_utils.core.entity.property.chemical_potentials import ( + flatten_scope_track, + get_formation_energy_at_references, + solve_chemical_potentials, +) -# Total energies (eV) as (atom counts, energy) of the Standata materials, PBE; ZrO2 is illustrative. +# Total energies (eV) as (atom counts, energy) of the Standata materials, PBE. O2 = ({"O": 8}, -3499.6187) HF = ({"Hf": 2}, -4320.4730) ZR = ({"Zr": 2}, -2698.0925) HFO2 = ({"Hf": 4, "O": 8}, -12183.3350) HFO2_3X3X3 = ({"Hf": 108, "O": 216}, 27 * -12183.3350) -ZRO2 = ({"Zr": 4, "O": 8}, -9800.0) +ZRO2 = ({"Zr": 4, "O": 8}, -8937.1617) MU_O_RICH = -3499.6187 / 8 MU_HF_O_POOR = -4320.4730 / 2 @@ -20,7 +24,7 @@ ({"O": O2, "Hf": HF, "Zr": ZR}, {"O": MU_O_RICH, "Hf": MU_HF_O_POOR, "Zr": -2698.0925 / 2}), ( {"O": O2, "Hf": HFO2, "Zr": ZRO2}, - {"O": MU_O_RICH, "Hf": -12183.3350 / 4 - 2 * MU_O_RICH, "Zr": -9800.0 / 4 - 2 * MU_O_RICH}, + {"O": MU_O_RICH, "Hf": -12183.3350 / 4 - 2 * MU_O_RICH, "Zr": -8937.1617 / 4 - 2 * MU_O_RICH}, ), ({"Hf": HF, "O": HFO2}, {"Hf": MU_HF_O_POOR, "O": (-12183.3350 / 4 - MU_HF_O_POOR) / 2}), ], @@ -33,3 +37,24 @@ def test_solve_chemical_potentials(reference_energies, expected): def test_solve_chemical_potentials_raises(reference_energies): with pytest.raises(ValueError): solve_chemical_potentials(reference_energies) + + +# scopeTrack globals of the m-HfO2 jobs rxCNizLKg7hkrPgAh (Zr_Hf) and mpxeDWYNKPBr6zc8Z (V_O, E_O from O2). +ZR_HF_SCOPE_TRACK = [ + {"scope": {"global": {"SUM_DELTA_N_TIMES_MU": 0.0, "DELTA_N_BY_SYMBOL": {"Hf": -1, "O": 0, "Zr": 1}}}}, + {"scope": {"global": {"SUM_DELTA_N_TIMES_MU": 811.1902858605226, "DEFECT_FORMATION_ENERGY": 0.3664}}}, +] +V_O_SCOPE_TRACK = [ + {"scope": {"global": {"DELTA_N_BY_SYMBOL": {"Hf": 0, "O": -1}, "SUM_DELTA_N_TIMES_MU": -MU_O_RICH}}}, + {"scope": {"global": {"DEFECT_FORMATION_ENERGY": 6.356}}}, +] + + +@pytest.mark.parametrize( + "scope_track, reference_energies, expected", + [(ZR_HF_SCOPE_TRACK, {"O": O2, "Hf": HFO2, "Zr": ZRO2}, 0.0134), (V_O_SCOPE_TRACK, {"O": O2, "Hf": HFO2}, 6.356)], +) +def test_get_formation_energy_at_references(scope_track, reference_energies, expected): + chemical_potentials = solve_chemical_potentials(reference_energies) + scope = flatten_scope_track(scope_track) + assert get_formation_energy_at_references(scope, chemical_potentials) == pytest.approx(expected, abs=1e-3) From 11b5a1772edf460d08889c34fc058dd2ecc6fc46 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Tue, 29 Sep 2026 13:23:36 -0700 Subject: [PATCH 14/18] feat(SOF-7975): chemical potentials as typed inputs, formation energy printed and plotted per condition A defect formation energy is a function of the chemical potentials, not one number: the job stores it at the elemental references it used, mu_i = E_i, and at any other mu_i it is E_f(dmu) = E_f - sum_i dN_i dmu_i with dmu_i = mu_i - E_i <= 0. Which mu_i apply is the user's input, read from a phase diagram run separately, so both notebooks now take CHEMICAL_POTENTIALS, dmu per element for each named condition (None = the stored value), instead of resolving reference materials and solving for mu. get_formation_energy_at_chemical_potentials applies the shift from the job's flattened scope and raises KeyError when a condition lacks an element of the job; the results cells print E_f per condition (charged: per size, charge and condition) and plot it against -sum_i dN_i dmu_i, one line per charge for the largest size. get_reference_energies, solve_chemical_potentials and get_kgrid_of_job had no other callers and are removed with their tests. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../workflows/defect_formation_energy.ipynb | 54 ++++++-------- .../defect_formation_energy_charged.ipynb | 69 ++++++++--------- .../notebooks_utils/core/entity/job/api.py | 7 -- .../entity/property/chemical_potentials.py | 74 +++---------------- tests/py/unit/core/entity/test_job_api.py | 9 --- .../test_property_chemical_potentials.py | 62 +++++----------- 6 files changed, 84 insertions(+), 191 deletions(-) diff --git a/other/materials_designer/workflows/defect_formation_energy.ipynb b/other/materials_designer/workflows/defect_formation_energy.ipynb index fdd4ad9d2..4033a43b4 100644 --- a/other/materials_designer/workflows/defect_formation_energy.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy.ipynb @@ -123,7 +123,7 @@ "outputs": [], "source": [ "SCF_KGRID = [4, 4, 4] # same grid for the pristine and defective cells\n", - "CHEMICAL_POTENTIAL_REFERENCES = None # {\"O\": \"\", ...}: one reference material per element; None = elemental\n" + "CHEMICAL_POTENTIALS = None # Δμ per element, eV/atom, per condition, e.g. {\"O-rich\": {\"O\": 0.0, \"Hf\": -10.69, \"Zr\": -10.34}, \"O-poor\": {\"O\": -5.35, \"Hf\": 0.0, \"Zr\": 0.0}}\n" ] }, { @@ -363,12 +363,6 @@ "source": [ "from mat3ra.notebooks_utils.core.entity.job.api import find_job_for_material_with_property\n", "from mat3ra.notebooks_utils.core.entity.material.api import get_or_create_material\n", - "from mat3ra.notebooks_utils.core.entity.property.chemical_potentials import (\n", - " flatten_scope_track,\n", - " get_formation_energy_at_references,\n", - " get_reference_energies,\n", - " solve_chemical_potentials,\n", - ")\n", "\n", "saved_defective = Material.create(get_or_create_material(client, defective_material, ACCOUNT_ID))\n", "saved_pristine = Material.create(get_or_create_material(client, pristine_material, ACCOUNT_ID))\n", @@ -384,22 +378,7 @@ "print(f\"♻️ pristine reference: job {pristine_reference_job['_id']}, k-grid {SCF_KGRID}\")\n", "\n", "# Order matters: [0] defective (computed), [1] pristine (reference).\n", - "materials = [saved_defective, saved_pristine]\n", - "\n", - "if CHEMICAL_POTENTIAL_REFERENCES:\n", - " reference_materials = {\n", - " element: Material.create(\n", - " get_or_create_material(client, load_material(client, FOLDER, name, ACCOUNT_ID), ACCOUNT_ID)\n", - " )\n", - " for element, name in CHEMICAL_POTENTIAL_REFERENCES.items()\n", - " }\n", - " pristine_energy = client.properties.get_for_job(pristine_reference_job[\"_id\"], \"total_energy\")[0][\"value\"]\n", - " reference_energies = get_reference_energies(\n", - " client, reference_materials, ACCOUNT_ID, SCF_KGRID, {saved_pristine.id: pristine_energy}\n", - " )\n", - " chemical_potentials = solve_chemical_potentials(reference_energies)\n", - " print(\"Chemical potentials (eV/atom):\",\n", - " {element: round(value, 3) for element, value in chemical_potentials.items()})\n" + "materials = [saved_defective, saved_pristine]\n" ] }, { @@ -580,7 +559,7 @@ "source": [ "## 7. Retrieve results\n", "### 7.1. Retrieve and visualize defect formation energy\n", - "With `CHEMICAL_POTENTIAL_REFERENCES`, one material per element that the host is in equilibrium with, each giving $E = \\sum_i n_i\\, \\mu_i$, the formation energy is printed at those references, $E_f = E_f^{\\text{elem}} + \\sum_i \\Delta N_i\\,(E_i - \\mu_i)$, with $E_i$ the elemental energy per atom the job used. The chemical potentials are solved in section 3.5, before any defect job is created. Each material's total energy comes from its finished Total Energy job in this account, on `SCF_KGRID` if there is one, else on the one other k-grid it was run on. A reference named `PRISTINE_NAME` is the pristine supercell and takes the energy of its Total Energy job, so the chemical potentials add up to the pristine energy." + "The job reports $E_f$ at the elemental chemical potentials it used, $\\mu_i = E_i$. `CHEMICAL_POTENTIALS` sets $\\Delta\\mu_i = \\mu_i - E_i \\le 0$ for every element of the job, per named condition; take them from a phase diagram computed separately. The table lists $E_f(\\Delta\\mu) = E_f - \\sum_i \\Delta N_i\\,\\Delta\\mu_i$ for each condition, and the plot draws this line through them." ] }, { @@ -590,16 +569,29 @@ "metadata": {}, "outputs": [], "source": [ + "import pandas as pd\n", + "import plotly.graph_objects as go\n", + "from mat3ra.notebooks_utils.core.entity.property.chemical_potentials import (\n", + " flatten_scope_track,\n", + " get_formation_energy_at_chemical_potentials,\n", + ")\n", "from mat3ra.notebooks_utils.ipython.entity.property.visualize import visualize_properties\n", - "\n", - "if CHEMICAL_POTENTIAL_REFERENCES:\n", - " scope = flatten_scope_track(client.jobs.get(defect_job_id)[\"scopeTrack\"])\n", - " references = \", \".join(f\"{element}: {name}\" for element, name in CHEMICAL_POTENTIAL_REFERENCES.items())\n", - " print(f\"Defect formation energy at the references ({references}): \"\n", - " f\"{get_formation_energy_at_references(scope, chemical_potentials):.4f} eV\")\n", + "from mat3ra.notebooks_utils.ipython.plot._plotly import render_figure\n", "\n", "defect_energy_data = client.properties.get_for_job(defect_job_id)\n", - "visualize_properties(defect_energy_data, title=\"Defect Formation Energy\")" + "visualize_properties(defect_energy_data, title=\"Defect Formation Energy\")\n", + "\n", + "if CHEMICAL_POTENTIALS:\n", + " scope = flatten_scope_track(client.jobs.get(defect_job_id)[\"scopeTrack\"])\n", + " chemical_potentials_df = pd.DataFrame.from_dict(CHEMICAL_POTENTIALS, orient=\"index\")\n", + " chemical_potentials_df[\"formation_energy\"] = [\n", + " get_formation_energy_at_chemical_potentials(scope, delta_mu) for delta_mu in CHEMICAL_POTENTIALS.values()\n", + " ]\n", + " print(chemical_potentials_df)\n", + " energies = chemical_potentials_df[\"formation_energy\"]\n", + " figure = go.Figure(go.Scatter(x=energies - scope[\"DEFECT_FORMATION_ENERGY\"], y=energies, text=energies.index))\n", + " figure.update_layout(xaxis_title=\"−Σ ΔN_i Δμ_i (eV)\", yaxis_title=\"Defect formation energy (eV)\")\n", + " render_figure(figure)" ] } ], diff --git a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb index 2d1394ef4..5eec99343 100644 --- a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb @@ -133,7 +133,7 @@ "CHARGES = [0] # e.g. [1, 0, -1, -2, -3]\n", "\n", "SCF_KGRID = [4, 4, 4] # same grid for the pristine and defective cells, divided by the supercell scaling\n", - "CHEMICAL_POTENTIAL_REFERENCES = None # {\"O\": \"\", ...}: one reference material per element; None = elemental\n", + "CHEMICAL_POTENTIALS = None # Δμ per element, eV/atom, per condition, e.g. {\"O-rich\": {\"O\": 0.0, \"Hf\": -10.69, \"Zr\": -10.34}, \"O-poor\": {\"O\": -5.35, \"Hf\": 0.0, \"Zr\": 0.0}}\n", "\n", "RELAX_PRISTINE_MATERIAL = False\n", "RELAX_DEFECTIVE_MATERIAL = False" @@ -649,36 +649,8 @@ "metadata": {}, "outputs": [], "source": [ - "from mat3ra.notebooks_utils.core.entity.property.chemical_potentials import (\n", - " flatten_scope_track,\n", - " get_formation_energy_at_references,\n", - " get_reference_energies,\n", - " solve_chemical_potentials,\n", - ")\n", - "\n", "if pristine_job_ids:\n", - " await wait_for_jobs_to_finish_async(client.jobs, pristine_job_ids, poll_interval=POLL_INTERVAL)\n", - "\n", - "if CHEMICAL_POTENTIAL_REFERENCES:\n", - " reference_materials = {}\n", - " for element, name in CHEMICAL_POTENTIAL_REFERENCES.items():\n", - " if name != PRISTINE_NAME:\n", - " matches = [Material.create(data) for data in Materials.get_by_name(name) if data[\"name\"] == name]\n", - " material = matches[0] if matches else load_material(client, FOLDER, name, ACCOUNT_ID)\n", - " reference_materials[element] = Material.create(get_or_create_material(client, material, ACCOUNT_ID))\n", - " chemical_potentials_by_scaling = {}\n", - " for scaling, (pristine_supercell, _) in saved_supercell_pairs.items():\n", - " pristine_job_id = pristine_reference_jobs[scaling][\"total_energy\"][\"_id\"]\n", - " pristine_energy = client.properties.get_for_job(pristine_job_id, \"total_energy\")[0][\"value\"]\n", - " materials_by_element = {\n", - " element: reference_materials.get(element, pristine_supercell) for element in CHEMICAL_POTENTIAL_REFERENCES\n", - " }\n", - " reference_energies = get_reference_energies(\n", - " client, materials_by_element, ACCOUNT_ID, SCF_KGRID, {pristine_supercell.id: pristine_energy}\n", - " )\n", - " chemical_potentials_by_scaling[scaling] = solve_chemical_potentials(reference_energies)\n", - " print(f\"n={scaling}: chemical potentials (eV/atom):\",\n", - " {element: round(value, 3) for element, value in chemical_potentials_by_scaling[scaling].items()})" + " await wait_for_jobs_to_finish_async(client.jobs, pristine_job_ids, poll_interval=POLL_INTERVAL)" ] }, { @@ -779,7 +751,7 @@ "source": [ "## 7. Retrieve results\n", "### 7.1. Formation energies\n", - "`formation_energy` is $E_f[X^q]$ at $\\mu_e = 0$, as the workflow reports it. `length` is $L = V^{1/3}$ of the supercell, the length scale the image-charge terms in section 7.3 are expressed in. With `CHEMICAL_POTENTIAL_REFERENCES`, one material per element that the host is in equilibrium with, each giving $E = \\sum_i n_i\\, \\mu_i$, `formation_energy_at_references` is $E_f = E_f^{\\text{elem}} + \\sum_i \\Delta N_i\\,(E_i - \\mu_i)$, with $E_i$ the elemental energy per atom the job used. The chemical potentials are solved in section 6.1 for each supercell size, before any defect job is created. Each material's total energy comes from its finished Total Energy job in this account, on `SCF_KGRID` if there is one, else on the one other k-grid it was run on. A reference named `PRISTINE_NAME` is the pristine supercell of that size and takes the energy of its Total Energy job, so the chemical potentials add up to the pristine energy." + "`formation_energy` is $E_f[X^q]$ at $\\mu_e = 0$, as the workflow reports it. `length` is $L = V^{1/3}$ of the supercell, the length scale the image-charge terms in section 7.3 are expressed in. The workflow uses the elemental chemical potentials $\\mu_i = E_i$. `CHEMICAL_POTENTIALS` sets $\\Delta\\mu_i = \\mu_i - E_i \\le 0$ for every element of the job, per named condition; take them from a phase diagram computed separately. A table then lists $E_f(\\Delta\\mu) = E_f - \\sum_i \\Delta N_i\\,\\Delta\\mu_i$ per size, charge and condition, and a plot draws these lines for the largest size, one per charge state." ] }, { @@ -789,7 +761,14 @@ "metadata": {}, "outputs": [], "source": [ - "results_records = []\n", + "import plotly.graph_objects as go\n", + "from mat3ra.notebooks_utils.core.entity.property.chemical_potentials import (\n", + " flatten_scope_track,\n", + " get_formation_energy_at_chemical_potentials,\n", + ")\n", + "from mat3ra.notebooks_utils.ipython.plot._plotly import render_figure\n", + "\n", + "results_records, chemical_potential_records = [], []\n", "for record in job_records:\n", " job = client.jobs.get(record[\"job_id\"])\n", " properties = client.properties.get_for_job(record[\"job_id\"], property_name=\"defect_formation_energy\")\n", @@ -798,16 +777,32 @@ " \"final_status\": job.get(\"status\"),\n", " \"formation_energy\": properties[0].get(\"value\") if properties else None,\n", " })\n", - " if CHEMICAL_POTENTIAL_REFERENCES and properties:\n", - " results_records[-1][\"formation_energy_at_references\"] = get_formation_energy_at_references(\n", - " flatten_scope_track(job[\"scopeTrack\"]), chemical_potentials_by_scaling[record[\"scaling\"]]\n", - " )\n", + " if CHEMICAL_POTENTIALS and properties:\n", + " scope = flatten_scope_track(job[\"scopeTrack\"])\n", + " for condition, delta_mu in CHEMICAL_POTENTIALS.items():\n", + " formation_energy = get_formation_energy_at_chemical_potentials(scope, delta_mu)\n", + " chemical_potential_records.append({\n", + " \"scaling\": record[\"scaling\"], \"charge\": record[\"charge\"], \"condition\": condition, **delta_mu,\n", + " \"formation_energy\": formation_energy,\n", + " \"formation_energy_shift\": formation_energy - scope[\"DEFECT_FORMATION_ENERGY\"],\n", + " })\n", "\n", "results_df = pd.DataFrame(results_records)\n", "successful_df = results_df[(results_df[\"final_status\"] == \"finished\") & results_df[\"formation_energy\"].notna()]\n", "if len(successful_df) < len(results_df):\n", " print(f\"⚠️ {len(results_df) - len(successful_df)} of {len(results_df)} job(s) returned no formation energy; \"\n", " \"they are left out of sections 7.2 and 7.3.\")\n", + "if CHEMICAL_POTENTIALS:\n", + " chemical_potentials_df = pd.DataFrame(chemical_potential_records)\n", + " print(chemical_potentials_df.to_string(index=False))\n", + " largest_df = chemical_potentials_df[chemical_potentials_df[\"scaling\"] == chemical_potentials_df[\"scaling\"].max()]\n", + " figure = go.Figure()\n", + " for charge, group in largest_df.groupby(\"charge\"):\n", + " figure.add_scatter(x=group[\"formation_energy_shift\"], y=group[\"formation_energy\"], text=group[\"condition\"],\n", + " name=f\"q = {charge:+d}\")\n", + " figure.update_layout(title=f\"Defect formation energy (n={largest_df['scaling'].iloc[0]})\",\n", + " xaxis_title=\"−Σ ΔN_i Δμ_i (eV)\", yaxis_title=\"Formation energy at the VBM (eV)\")\n", + " render_figure(figure)\n", "results_df" ] }, @@ -817,7 +812,7 @@ "metadata": {}, "source": [ "### 7.2. Formation energy versus the electron chemical potential\n", - "Each charge state gives a straight line of slope $q$ over $\\mu_e \\in [0, E_{\\text{gap}}]$; the lower envelope is the charge state the defect actually adopts, and the crossings between the lines are the charge transition levels. The largest supercell with a result for every charge state is used, with $E_{\\text{gap}}$ from its own Band Gap job. The table lists each charge state's formation energy at the VBM and the range of the Fermi level $\\mu_e$ in which it is the stable state. The lines, and the fit in section 7.3, use the elemental $\\mu_i$ the workflow stores; `formation_energy_at_references` in section 7.1 is shifted by the same amount for every charge state of a size, so the transition levels do not change." + "Each charge state gives a straight line of slope $q$ over $\\mu_e \\in [0, E_{\\text{gap}}]$; the lower envelope is the charge state the defect actually adopts, and the crossings between the lines are the charge transition levels. The largest supercell with a result for every charge state is used, with $E_{\\text{gap}}$ from its own Band Gap job. The table lists each charge state's formation energy at the VBM and the range of the Fermi level $\\mu_e$ in which it is the stable state." ] }, { diff --git a/src/py/mat3ra/notebooks_utils/core/entity/job/api.py b/src/py/mat3ra/notebooks_utils/core/entity/job/api.py index 74d774a1e..57c33dcd4 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/job/api.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/job/api.py @@ -171,13 +171,6 @@ def get_kgrid_query(kgrid: Optional[List[int]], unit_name: str = "pw_scf") -> Di return {"workflow.subworkflows.units": {"$elemMatch": {"name": unit_name, "context": kgrid_context}}} -def get_kgrid_of_job(job: dict) -> Optional[List[int]]: - """The k-grid dimensions the job's `pw_scf` unit ran on, read where `get_kgrid_query` matches; None if unset.""" - units = [unit for subworkflow in job["workflow"]["subworkflows"] for unit in subworkflow["units"]] - contexts = [item for unit in units if unit["name"] == "pw_scf" for item in unit["context"]] - return next((item["data"]["dimensions"] for item in contexts if item["name"] == "kgrid"), None) - - def find_job_for_material_with_property( api_client: APIClient, material_id: str, diff --git a/src/py/mat3ra/notebooks_utils/core/entity/property/chemical_potentials.py b/src/py/mat3ra/notebooks_utils/core/entity/property/chemical_potentials.py index f8293904f..423185576 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/property/chemical_potentials.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/property/chemical_potentials.py @@ -1,73 +1,21 @@ -from collections import Counter -from typing import Any, Dict, List, Tuple +from typing import Any, Dict, List -import numpy as np -from mat3ra.api_client import APIClient -from mat3ra.made.material import Material - -from ..job.api import get_kgrid_of_job +def flatten_scope_track(scope_track: List[dict]) -> Dict[str, Any]: + """The global scope of a job's `scopeTrack`, later values overriding earlier ones.""" + return {key: value for item in scope_track for key, value in item["scope"]["global"].items()} -def solve_chemical_potentials(reference_energies: Dict[str, Tuple[Dict[str, int], float]]) -> Dict[str, float]: - """ - Chemical potential of each element (eV/atom) from one reference material per element, given as its atom counts - and its total energy (eV) for them: each material gives E = sum_i n_i * mu_i, and the square system is solved. - Raises: - ValueError: If the materials contain other elements than the keys, or do not determine every potential. +def get_formation_energy_at_chemical_potentials(scope: Dict[str, Any], delta_mu: Dict[str, float]) -> float: """ - elements = sorted(reference_energies) - compositions = [composition for composition, _ in reference_energies.values()] - if set().union(*compositions) != set(elements): - raise ValueError(f"The reference materials must contain exactly the elements {elements}.") - matrix = np.array([[composition.get(element, 0) for element in elements] for composition in compositions]) - if np.linalg.matrix_rank(matrix) < len(elements): - raise ValueError("The reference materials do not determine every chemical potential.") - energies = [energy for _, energy in reference_energies.values()] - return dict(zip(elements, np.linalg.solve(matrix, energies).tolist())) + Defect formation energy at the chemical potentials mu_i = E_i + delta_mu[i], from its job's flattened scope. - -def get_reference_energies( - client: APIClient, - materials_by_element: Dict[str, Material], - owner_id: str, - kgrid: List[int], - host_energy: Dict[str, float], -) -> Dict[str, Tuple[Dict[str, int], float]]: - """ - Atom counts and total energy of each platform material, as `solve_chemical_potentials` takes them. The energy is - `host_energy[material.id]` when given there, else that of the material's finished Total Energy job of `owner_id`: - the one on `kgrid`, else the only one, or one of several on the same other k-grid. + E_i is the elemental energy per atom the job used, so the stored value is the one at delta_mu = 0: + E_f(delta_mu) = E_f - sum_i dN_i * delta_mu[i]. Raises: - RuntimeError: If a material has no such job, or has them only on several other k-grids. + KeyError: If `delta_mu` has no value for an element of the job, including one with dN_i = 0. """ - reference_energies = {} - for element, material in materials_by_element.items(): - if material.id in host_energy: - energy, source = host_energy[material.id], "pristine reference" - else: - query = {"_material._id": material.id, "owner._id": owner_id, "status": "finished"} - jobs = [job for job in client.jobs.list(query) if client.properties.get_for_job(job["_id"], "total_energy")] - jobs = [job for job in jobs if get_kgrid_of_job(job) == kgrid] or jobs - if len({str(get_kgrid_of_job(job)) for job in jobs}) != 1: - found = ", ".join(f"job {job['_id']} on {get_kgrid_of_job(job)}" for job in jobs) or "none" - raise RuntimeError(f"Run Total Energy on '{material.name}' at k-grid {kgrid}; found: {found}.") - energy = client.properties.get_for_job(jobs[0]["_id"], property_name="total_energy")[0]["value"] - source = f"job {jobs[0]['_id']}, k-grid {get_kgrid_of_job(jobs[0])}" - print(f"{element}: {material.name}, {source}, E = {energy:.4f} eV") - reference_energies[element] = (dict(Counter(material.basis.elements.values)), energy) - return reference_energies - - -def flatten_scope_track(scope_track: List[dict]) -> Dict[str, Any]: - """The global scope of a job's `scopeTrack`, later values overriding earlier ones.""" - return {key: value for item in scope_track for key, value in item["scope"]["global"].items()} - - -def get_formation_energy_at_references(scope: Dict[str, Any], chemical_potentials: Dict[str, float]) -> float: - """Defect formation energy at `chemical_potentials` from its job's flattened scope, stored at elemental E_i.""" - delta_n_times_mu = sum( - count * chemical_potentials[element] for element, count in scope["DELTA_N_BY_SYMBOL"].items() + return scope["DEFECT_FORMATION_ENERGY"] - sum( + count * delta_mu[element] for element, count in scope["DELTA_N_BY_SYMBOL"].items() ) - return scope["DEFECT_FORMATION_ENERGY"] + scope["SUM_DELTA_N_TIMES_MU"] - delta_n_times_mu diff --git a/tests/py/unit/core/entity/test_job_api.py b/tests/py/unit/core/entity/test_job_api.py index 604a1d81d..dc3a87934 100644 --- a/tests/py/unit/core/entity/test_job_api.py +++ b/tests/py/unit/core/entity/test_job_api.py @@ -6,7 +6,6 @@ create_job, find_job_for_material, find_job_for_material_with_property, - get_kgrid_of_job, get_kgrid_query, ) @@ -199,14 +198,6 @@ def test_get_kgrid_query(kgrid, unit_name, expected_query): assert get_kgrid_query(kgrid, unit_name) == expected_query -@pytest.mark.parametrize( - ("context", "expected_kgrid"), [([{"name": "kgrid", "data": {"dimensions": [4, 4, 4]}}], [4, 4, 4]), ([], None)] -) -def test_get_kgrid_of_job(context, expected_kgrid): - job = {"workflow": {"subworkflows": [{"units": [{"name": "pw_scf", "context": context}]}]}} - assert get_kgrid_of_job(job) == expected_kgrid - - def test_find_job_for_material_with_property_matches_the_pw_scf_kgrid(): client = MagicMock() client.jobs.list.return_value = [EXISTING_JOB] diff --git a/tests/py/unit/core/entity/test_property_chemical_potentials.py b/tests/py/unit/core/entity/test_property_chemical_potentials.py index e57a817de..f1c2ec34d 100644 --- a/tests/py/unit/core/entity/test_property_chemical_potentials.py +++ b/tests/py/unit/core/entity/test_property_chemical_potentials.py @@ -1,60 +1,34 @@ -"""Unit tests for chemical potentials from reference materials.""" +"""Unit tests for defect formation energies at chemical potentials relative to the elemental references.""" import pytest from mat3ra.notebooks_utils.core.entity.property.chemical_potentials import ( flatten_scope_track, - get_formation_energy_at_references, - solve_chemical_potentials, + get_formation_energy_at_chemical_potentials, ) -# Total energies (eV) as (atom counts, energy) of the Standata materials, PBE. -O2 = ({"O": 8}, -3499.6187) -HF = ({"Hf": 2}, -4320.4730) -ZR = ({"Zr": 2}, -2698.0925) -HFO2 = ({"Hf": 4, "O": 8}, -12183.3350) -HFO2_3X3X3 = ({"Hf": 108, "O": 216}, 27 * -12183.3350) -ZRO2 = ({"Zr": 4, "O": 8}, -8937.1617) -MU_O_RICH = -3499.6187 / 8 -MU_HF_O_POOR = -4320.4730 / 2 - - -@pytest.mark.parametrize( - "reference_energies, expected", - [ - ({"O": O2, "Hf": HF, "Zr": ZR}, {"O": MU_O_RICH, "Hf": MU_HF_O_POOR, "Zr": -2698.0925 / 2}), - ( - {"O": O2, "Hf": HFO2, "Zr": ZRO2}, - {"O": MU_O_RICH, "Hf": -12183.3350 / 4 - 2 * MU_O_RICH, "Zr": -8937.1617 / 4 - 2 * MU_O_RICH}, - ), - ({"Hf": HF, "O": HFO2}, {"Hf": MU_HF_O_POOR, "O": (-12183.3350 / 4 - MU_HF_O_POOR) / 2}), - ], -) -def test_solve_chemical_potentials(reference_energies, expected): - assert solve_chemical_potentials(reference_energies) == pytest.approx(expected) - - -@pytest.mark.parametrize("reference_energies", [{"Hf": HFO2, "O": HFO2}, {"Hf": HFO2, "O": HFO2_3X3X3}, {"Hf": HFO2}]) -def test_solve_chemical_potentials_raises(reference_energies): - with pytest.raises(ValueError): - solve_chemical_potentials(reference_energies) - - -# scopeTrack globals of the m-HfO2 jobs rxCNizLKg7hkrPgAh (Zr_Hf) and mpxeDWYNKPBr6zc8Z (V_O, E_O from O2). +# scopeTrack globals of the m-HfO2 jobs rxCNizLKg7hkrPgAh (Zr_Hf) and mpxeDWYNKPBr6zc8Z (V_O). ZR_HF_SCOPE_TRACK = [ - {"scope": {"global": {"SUM_DELTA_N_TIMES_MU": 0.0, "DELTA_N_BY_SYMBOL": {"Hf": -1, "O": 0, "Zr": 1}}}}, - {"scope": {"global": {"SUM_DELTA_N_TIMES_MU": 811.1902858605226, "DEFECT_FORMATION_ENERGY": 0.3664}}}, + {"scope": {"global": {"DELTA_N_BY_SYMBOL": {"Hf": -1, "O": 0, "Zr": 1}}}}, + {"scope": {"global": {"DEFECT_FORMATION_ENERGY": 0.3664}}}, ] V_O_SCOPE_TRACK = [ - {"scope": {"global": {"DELTA_N_BY_SYMBOL": {"Hf": 0, "O": -1}, "SUM_DELTA_N_TIMES_MU": -MU_O_RICH}}}, + {"scope": {"global": {"DELTA_N_BY_SYMBOL": {"O": -1, "Hf": 0}}}}, {"scope": {"global": {"DEFECT_FORMATION_ENERGY": 6.356}}}, ] +ELEMENTAL = {"O": 0.0, "Hf": 0.0, "Zr": 0.0} +O_RICH = {"O": 0.0, "Hf": -10.6926, "Zr": -10.3396} +O_POOR = {"O": -5.346, "Hf": 0.0} @pytest.mark.parametrize( - "scope_track, reference_energies, expected", - [(ZR_HF_SCOPE_TRACK, {"O": O2, "Hf": HFO2, "Zr": ZRO2}, 0.0134), (V_O_SCOPE_TRACK, {"O": O2, "Hf": HFO2}, 6.356)], + "scope_track, delta_mu, expected", + [(ZR_HF_SCOPE_TRACK, ELEMENTAL, 0.3664), (ZR_HF_SCOPE_TRACK, O_RICH, 0.0134), (V_O_SCOPE_TRACK, O_POOR, 1.010)], ) -def test_get_formation_energy_at_references(scope_track, reference_energies, expected): - chemical_potentials = solve_chemical_potentials(reference_energies) +def test_get_formation_energy_at_chemical_potentials(scope_track, delta_mu, expected): scope = flatten_scope_track(scope_track) - assert get_formation_energy_at_references(scope, chemical_potentials) == pytest.approx(expected, abs=1e-3) + assert get_formation_energy_at_chemical_potentials(scope, delta_mu) == pytest.approx(expected, abs=1e-4) + + +def test_get_formation_energy_at_chemical_potentials_raises_on_a_missing_element(): + with pytest.raises(KeyError): + get_formation_energy_at_chemical_potentials(flatten_scope_track(ZR_HF_SCOPE_TRACK), O_POOR) From a753cfc3f569ddca78f22c0a6bad71075fbaeb64 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Tue, 29 Sep 2026 13:44:44 -0700 Subject: [PATCH 15/18] fix(SOF-7975): label the chemical-potential plot by its combination of delta mu, define the inputs, move the helpers to defect_analysis MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The formation energy at a condition is E_f - sum_i dN_i dmu_i, so the plot's x axis is that combination of dmu, and its label is now built from the job's dN: "Δμ_O (eV)" for an O vacancy, "Δμ_Hf − Δμ_Zr (eV)" for Zr on a Hf site (get_chemical_potential_combination). The conditions are drawn as named points on the line. The shifted value is formation_energy_at_condition in both tables, so formation_energy keeps the meaning it has in results_df (the stored value); the charged table carries both and the plot's x is their difference, replacing formation_energy_shift. The 7.1 markdown defines E_i and dN_i and says where the dmu come from (the host's stability range, formation energies from analyze_convex_hull.ipynb); 7.2 says it uses the stored E_f and that a condition does not move the transition levels. The charged 7.1 cell guards on the records, so a run where every job errored still shows results_df. flatten_scope_track and get_formation_energy_at_chemical_potentials move to defect_analysis.py, the mu_i counterpart of the Fermi-level functions there, and chemical_potentials.py and its test file are removed. The parameter example sits on its own comment line with the values that reproduce 0.013 / 1.010 eV. A new test case pins that an element in the condition but not in the job is ignored. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../workflows/defect_formation_energy.ipynb | 18 +++++--- .../defect_formation_energy_charged.ipynb | 25 ++++++----- .../entity/property/chemical_potentials.py | 21 --------- .../core/entity/property/defect_analysis.py | 35 ++++++++++++++- .../test_property_chemical_potentials.py | 34 -------------- .../entity/test_property_defect_analysis.py | 44 +++++++++++++++++++ 6 files changed, 102 insertions(+), 75 deletions(-) delete mode 100644 src/py/mat3ra/notebooks_utils/core/entity/property/chemical_potentials.py delete mode 100644 tests/py/unit/core/entity/test_property_chemical_potentials.py diff --git a/other/materials_designer/workflows/defect_formation_energy.ipynb b/other/materials_designer/workflows/defect_formation_energy.ipynb index 4033a43b4..9ce9273cb 100644 --- a/other/materials_designer/workflows/defect_formation_energy.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy.ipynb @@ -123,7 +123,8 @@ "outputs": [], "source": [ "SCF_KGRID = [4, 4, 4] # same grid for the pristine and defective cells\n", - "CHEMICAL_POTENTIALS = None # Δμ per element, eV/atom, per condition, e.g. {\"O-rich\": {\"O\": 0.0, \"Hf\": -10.69, \"Zr\": -10.34}, \"O-poor\": {\"O\": -5.35, \"Hf\": 0.0, \"Zr\": 0.0}}\n" + "CHEMICAL_POTENTIALS = None # Δμ per element, eV/atom, per condition, e.g.\n", + "# {\"O-rich\": {\"O\": 0.0, \"Hf\": -10.693, \"Zr\": -10.340}, \"O-poor\": {\"O\": -5.346, \"Hf\": 0.0, \"Zr\": 0.0}}\n" ] }, { @@ -559,7 +560,7 @@ "source": [ "## 7. Retrieve results\n", "### 7.1. Retrieve and visualize defect formation energy\n", - "The job reports $E_f$ at the elemental chemical potentials it used, $\\mu_i = E_i$. `CHEMICAL_POTENTIALS` sets $\\Delta\\mu_i = \\mu_i - E_i \\le 0$ for every element of the job, per named condition; take them from a phase diagram computed separately. The table lists $E_f(\\Delta\\mu) = E_f - \\sum_i \\Delta N_i\\,\\Delta\\mu_i$ for each condition, and the plot draws this line through them." + "The job reports $E_f$ at the elemental chemical potentials it used, $\\mu_i = E_i$: $E_i$ is the energy per atom of the element's reference (e.g. O₂ for O, the bulk metal for Hf), and $\\Delta N_i$ is the number of atoms of $i$ the defect adds (+) or removes (−). `CHEMICAL_POTENTIALS` sets $\\Delta\\mu_i = \\mu_i - E_i \\le 0$ for every element of the job, per named condition. Obtain them separately from the host's stability range, $\\sum_i n_i\\,\\Delta\\mu_i = \\Delta H_f$ over its formula unit $n_i$, with $\\Delta H_f$ its formation energy (e.g. from [Formation Energies and Convex Hull](analyze_convex_hull.ipynb)). The table lists `formation_energy_at_condition`, $E_f(\\Delta\\mu) = E_f - \\sum_i \\Delta N_i\\,\\Delta\\mu_i$, for each condition, and the plot draws it against $-\\sum_i \\Delta N_i\\,\\Delta\\mu_i$ with the conditions marked." ] }, { @@ -571,8 +572,9 @@ "source": [ "import pandas as pd\n", "import plotly.graph_objects as go\n", - "from mat3ra.notebooks_utils.core.entity.property.chemical_potentials import (\n", + "from mat3ra.notebooks_utils.core.entity.property.defect_analysis import (\n", " flatten_scope_track,\n", + " get_chemical_potential_combination,\n", " get_formation_energy_at_chemical_potentials,\n", ")\n", "from mat3ra.notebooks_utils.ipython.entity.property.visualize import visualize_properties\n", @@ -584,13 +586,15 @@ "if CHEMICAL_POTENTIALS:\n", " scope = flatten_scope_track(client.jobs.get(defect_job_id)[\"scopeTrack\"])\n", " chemical_potentials_df = pd.DataFrame.from_dict(CHEMICAL_POTENTIALS, orient=\"index\")\n", - " chemical_potentials_df[\"formation_energy\"] = [\n", + " chemical_potentials_df[\"formation_energy_at_condition\"] = [\n", " get_formation_energy_at_chemical_potentials(scope, delta_mu) for delta_mu in CHEMICAL_POTENTIALS.values()\n", " ]\n", " print(chemical_potentials_df)\n", - " energies = chemical_potentials_df[\"formation_energy\"]\n", - " figure = go.Figure(go.Scatter(x=energies - scope[\"DEFECT_FORMATION_ENERGY\"], y=energies, text=energies.index))\n", - " figure.update_layout(xaxis_title=\"−Σ ΔN_i Δμ_i (eV)\", yaxis_title=\"Defect formation energy (eV)\")\n", + " energies = chemical_potentials_df[\"formation_energy_at_condition\"]\n", + " figure = go.Figure(go.Scatter(x=energies - scope[\"DEFECT_FORMATION_ENERGY\"], y=energies, text=energies.index,\n", + " mode=\"lines+markers+text\", textposition=\"top center\"))\n", + " figure.update_layout(xaxis_title=f\"{get_chemical_potential_combination(scope['DELTA_N_BY_SYMBOL'])} (eV)\",\n", + " yaxis_title=\"Defect formation energy (eV)\")\n", " render_figure(figure)" ] } diff --git a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb index 5eec99343..ac6eac943 100644 --- a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb @@ -133,7 +133,8 @@ "CHARGES = [0] # e.g. [1, 0, -1, -2, -3]\n", "\n", "SCF_KGRID = [4, 4, 4] # same grid for the pristine and defective cells, divided by the supercell scaling\n", - "CHEMICAL_POTENTIALS = None # Δμ per element, eV/atom, per condition, e.g. {\"O-rich\": {\"O\": 0.0, \"Hf\": -10.69, \"Zr\": -10.34}, \"O-poor\": {\"O\": -5.35, \"Hf\": 0.0, \"Zr\": 0.0}}\n", + "CHEMICAL_POTENTIALS = None # Δμ per element, eV/atom, per condition, e.g.\n", + "# {\"O-rich\": {\"O\": 0.0, \"Hf\": -10.693, \"Zr\": -10.340}, \"O-poor\": {\"O\": -5.346, \"Hf\": 0.0, \"Zr\": 0.0}}\n", "\n", "RELAX_PRISTINE_MATERIAL = False\n", "RELAX_DEFECTIVE_MATERIAL = False" @@ -751,7 +752,7 @@ "source": [ "## 7. Retrieve results\n", "### 7.1. Formation energies\n", - "`formation_energy` is $E_f[X^q]$ at $\\mu_e = 0$, as the workflow reports it. `length` is $L = V^{1/3}$ of the supercell, the length scale the image-charge terms in section 7.3 are expressed in. The workflow uses the elemental chemical potentials $\\mu_i = E_i$. `CHEMICAL_POTENTIALS` sets $\\Delta\\mu_i = \\mu_i - E_i \\le 0$ for every element of the job, per named condition; take them from a phase diagram computed separately. A table then lists $E_f(\\Delta\\mu) = E_f - \\sum_i \\Delta N_i\\,\\Delta\\mu_i$ per size, charge and condition, and a plot draws these lines for the largest size, one per charge state." + "`formation_energy` is $E_f[X^q]$ at $\\mu_e = 0$, as the workflow reports it. `length` is $L = V^{1/3}$ of the supercell, the length scale the image-charge terms in section 7.3 are expressed in. The workflow uses the elemental chemical potentials $\\mu_i = E_i$: $E_i$ is the energy per atom of the element's reference (e.g. O₂ for O, the bulk metal for Hf), and $\\Delta N_i$ is the number of atoms of $i$ the defect adds (+) or removes (−). `CHEMICAL_POTENTIALS` sets $\\Delta\\mu_i = \\mu_i - E_i \\le 0$ for every element of the job, per named condition. Obtain them separately from the host's stability range, $\\sum_i n_i\\,\\Delta\\mu_i = \\Delta H_f$ over its formula unit $n_i$, with $\\Delta H_f$ its formation energy (e.g. from [Formation Energies and Convex Hull](analyze_convex_hull.ipynb)). A table then lists `formation_energy_at_condition`, $E_f(\\Delta\\mu) = E_f - \\sum_i \\Delta N_i\\,\\Delta\\mu_i$, per size, charge and condition, and a plot draws it against $-\\sum_i \\Delta N_i\\,\\Delta\\mu_i$ for the largest size, one line per charge state, with the conditions marked." ] }, { @@ -762,8 +763,9 @@ "outputs": [], "source": [ "import plotly.graph_objects as go\n", - "from mat3ra.notebooks_utils.core.entity.property.chemical_potentials import (\n", + "from mat3ra.notebooks_utils.core.entity.property.defect_analysis import (\n", " flatten_scope_track,\n", + " get_chemical_potential_combination,\n", " get_formation_energy_at_chemical_potentials,\n", ")\n", "from mat3ra.notebooks_utils.ipython.plot._plotly import render_figure\n", @@ -780,11 +782,10 @@ " if CHEMICAL_POTENTIALS and properties:\n", " scope = flatten_scope_track(job[\"scopeTrack\"])\n", " for condition, delta_mu in CHEMICAL_POTENTIALS.items():\n", - " formation_energy = get_formation_energy_at_chemical_potentials(scope, delta_mu)\n", " chemical_potential_records.append({\n", " \"scaling\": record[\"scaling\"], \"charge\": record[\"charge\"], \"condition\": condition, **delta_mu,\n", - " \"formation_energy\": formation_energy,\n", - " \"formation_energy_shift\": formation_energy - scope[\"DEFECT_FORMATION_ENERGY\"],\n", + " \"formation_energy\": scope[\"DEFECT_FORMATION_ENERGY\"],\n", + " \"formation_energy_at_condition\": get_formation_energy_at_chemical_potentials(scope, delta_mu),\n", " })\n", "\n", "results_df = pd.DataFrame(results_records)\n", @@ -792,16 +793,18 @@ "if len(successful_df) < len(results_df):\n", " print(f\"⚠️ {len(results_df) - len(successful_df)} of {len(results_df)} job(s) returned no formation energy; \"\n", " \"they are left out of sections 7.2 and 7.3.\")\n", - "if CHEMICAL_POTENTIALS:\n", + "if chemical_potential_records:\n", " chemical_potentials_df = pd.DataFrame(chemical_potential_records)\n", " print(chemical_potentials_df.to_string(index=False))\n", " largest_df = chemical_potentials_df[chemical_potentials_df[\"scaling\"] == chemical_potentials_df[\"scaling\"].max()]\n", " figure = go.Figure()\n", " for charge, group in largest_df.groupby(\"charge\"):\n", - " figure.add_scatter(x=group[\"formation_energy_shift\"], y=group[\"formation_energy\"], text=group[\"condition\"],\n", - " name=f\"q = {charge:+d}\")\n", + " figure.add_scatter(x=group[\"formation_energy_at_condition\"] - group[\"formation_energy\"],\n", + " y=group[\"formation_energy_at_condition\"], text=group[\"condition\"],\n", + " mode=\"lines+markers+text\", textposition=\"top center\", name=f\"q = {charge:+d}\")\n", " figure.update_layout(title=f\"Defect formation energy (n={largest_df['scaling'].iloc[0]})\",\n", - " xaxis_title=\"−Σ ΔN_i Δμ_i (eV)\", yaxis_title=\"Formation energy at the VBM (eV)\")\n", + " xaxis_title=f\"{get_chemical_potential_combination(scope['DELTA_N_BY_SYMBOL'])} (eV)\",\n", + " yaxis_title=\"Formation energy at the VBM (eV)\")\n", " render_figure(figure)\n", "results_df" ] @@ -812,7 +815,7 @@ "metadata": {}, "source": [ "### 7.2. Formation energy versus the electron chemical potential\n", - "Each charge state gives a straight line of slope $q$ over $\\mu_e \\in [0, E_{\\text{gap}}]$; the lower envelope is the charge state the defect actually adopts, and the crossings between the lines are the charge transition levels. The largest supercell with a result for every charge state is used, with $E_{\\text{gap}}$ from its own Band Gap job. The table lists each charge state's formation energy at the VBM and the range of the Fermi level $\\mu_e$ in which it is the stable state." + "Each charge state gives a straight line of slope $q$ over $\\mu_e \\in [0, E_{\\text{gap}}]$; the lower envelope is the charge state the defect actually adopts, and the crossings between the lines are the charge transition levels. The largest supercell with a result for every charge state is used, with $E_{\\text{gap}}$ from its own Band Gap job. The table lists each charge state's formation energy at the VBM and the range of the Fermi level $\\mu_e$ in which it is the stable state. 7.2 and 7.3 use the stored $E_f$ ($\\mu_i = E_i$). A condition shifts every charge state by the same $-\\sum_i \\Delta N_i\\,\\Delta\\mu_i$, so the transition levels do not depend on it." ] }, { diff --git a/src/py/mat3ra/notebooks_utils/core/entity/property/chemical_potentials.py b/src/py/mat3ra/notebooks_utils/core/entity/property/chemical_potentials.py deleted file mode 100644 index 423185576..000000000 --- a/src/py/mat3ra/notebooks_utils/core/entity/property/chemical_potentials.py +++ /dev/null @@ -1,21 +0,0 @@ -from typing import Any, Dict, List - - -def flatten_scope_track(scope_track: List[dict]) -> Dict[str, Any]: - """The global scope of a job's `scopeTrack`, later values overriding earlier ones.""" - return {key: value for item in scope_track for key, value in item["scope"]["global"].items()} - - -def get_formation_energy_at_chemical_potentials(scope: Dict[str, Any], delta_mu: Dict[str, float]) -> float: - """ - Defect formation energy at the chemical potentials mu_i = E_i + delta_mu[i], from its job's flattened scope. - - E_i is the elemental energy per atom the job used, so the stored value is the one at delta_mu = 0: - E_f(delta_mu) = E_f - sum_i dN_i * delta_mu[i]. - - Raises: - KeyError: If `delta_mu` has no value for an element of the job, including one with dN_i = 0. - """ - return scope["DEFECT_FORMATION_ENERGY"] - sum( - count * delta_mu[element] for element, count in scope["DELTA_N_BY_SYMBOL"].items() - ) diff --git a/src/py/mat3ra/notebooks_utils/core/entity/property/defect_analysis.py b/src/py/mat3ra/notebooks_utils/core/entity/property/defect_analysis.py index 9a8766f02..09583f628 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/property/defect_analysis.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/property/defect_analysis.py @@ -1,4 +1,5 @@ -"""Formation energies of charged defects: dependence on the Fermi level and on the supercell size. +"""Formation energies of charged defects: dependence on the Fermi level, on the supercell size and on the chemical +potentials. Pure computation — no API calls, no display. Kept apart from `analysis.py`, which imports pymatgen's phase diagram module at load time; that module needs tqdm, which JupyterLite does not install. @@ -6,7 +7,7 @@ The Fermi level is the electron chemical potential mu_e, measured from the valence band maximum (VBM). """ -from typing import Dict, Sequence +from typing import Any, Dict, List, Sequence import numpy as np import pandas as pd @@ -108,3 +109,33 @@ def evaluate_finite_size_fit(fit: Dict[str, float], inverse_lengths: Sequence[fl """E_inf + a / L + b / L^3 from the coefficients of `fit_finite_size`, at the given 1 / L.""" inverse_lengths_array = np.asarray(inverse_lengths, dtype=float) return fit["E_inf"] + fit["a"] * inverse_lengths_array + fit.get("b", 0.0) * inverse_lengths_array**3 + + +def flatten_scope_track(scope_track: List[dict]) -> Dict[str, Any]: + """The global scope of a job's `scopeTrack`, later values overriding earlier ones.""" + return {key: value for item in scope_track for key, value in item["scope"]["global"].items()} + + +def get_formation_energy_at_chemical_potentials(scope: Dict[str, Any], delta_mu: Dict[str, float]) -> float: + """ + Defect formation energy at the chemical potentials mu_i = E_i + delta_mu[i], from its job's flattened scope. + + E_i is the elemental energy per atom the job used, so the stored value is the one at delta_mu = 0: + E_f(delta_mu) = E_f - sum_i dN_i * delta_mu[i]. + + Raises: + KeyError: If `delta_mu` has no value for an element of the job, including one with dN_i = 0. + """ + return scope["DEFECT_FORMATION_ENERGY"] - sum( + count * delta_mu[element] for element, count in scope["DELTA_N_BY_SYMBOL"].items() + ) + + +def get_chemical_potential_combination(delta_n_by_symbol: Dict[str, int]) -> str: + """-sum_i dN_i * delta_mu_i written out, removed atoms first: "Δμ_Hf − Δμ_Zr" for dN = {Hf: -1, Zr: +1}.""" + terms = [ + f"{'+' if count < 0 else '−'} {abs(count) if abs(count) != 1 else ''}Δμ_{element}" + for element, count in sorted(delta_n_by_symbol.items(), key=lambda item: item[1]) + if count + ] + return " ".join(terms).lstrip("+ ") diff --git a/tests/py/unit/core/entity/test_property_chemical_potentials.py b/tests/py/unit/core/entity/test_property_chemical_potentials.py deleted file mode 100644 index f1c2ec34d..000000000 --- a/tests/py/unit/core/entity/test_property_chemical_potentials.py +++ /dev/null @@ -1,34 +0,0 @@ -"""Unit tests for defect formation energies at chemical potentials relative to the elemental references.""" - -import pytest -from mat3ra.notebooks_utils.core.entity.property.chemical_potentials import ( - flatten_scope_track, - get_formation_energy_at_chemical_potentials, -) - -# scopeTrack globals of the m-HfO2 jobs rxCNizLKg7hkrPgAh (Zr_Hf) and mpxeDWYNKPBr6zc8Z (V_O). -ZR_HF_SCOPE_TRACK = [ - {"scope": {"global": {"DELTA_N_BY_SYMBOL": {"Hf": -1, "O": 0, "Zr": 1}}}}, - {"scope": {"global": {"DEFECT_FORMATION_ENERGY": 0.3664}}}, -] -V_O_SCOPE_TRACK = [ - {"scope": {"global": {"DELTA_N_BY_SYMBOL": {"O": -1, "Hf": 0}}}}, - {"scope": {"global": {"DEFECT_FORMATION_ENERGY": 6.356}}}, -] -ELEMENTAL = {"O": 0.0, "Hf": 0.0, "Zr": 0.0} -O_RICH = {"O": 0.0, "Hf": -10.6926, "Zr": -10.3396} -O_POOR = {"O": -5.346, "Hf": 0.0} - - -@pytest.mark.parametrize( - "scope_track, delta_mu, expected", - [(ZR_HF_SCOPE_TRACK, ELEMENTAL, 0.3664), (ZR_HF_SCOPE_TRACK, O_RICH, 0.0134), (V_O_SCOPE_TRACK, O_POOR, 1.010)], -) -def test_get_formation_energy_at_chemical_potentials(scope_track, delta_mu, expected): - scope = flatten_scope_track(scope_track) - assert get_formation_energy_at_chemical_potentials(scope, delta_mu) == pytest.approx(expected, abs=1e-4) - - -def test_get_formation_energy_at_chemical_potentials_raises_on_a_missing_element(): - with pytest.raises(KeyError): - get_formation_energy_at_chemical_potentials(flatten_scope_track(ZR_HF_SCOPE_TRACK), O_POOR) diff --git a/tests/py/unit/core/entity/test_property_defect_analysis.py b/tests/py/unit/core/entity/test_property_defect_analysis.py index 5a5d981c0..2b63b7cf0 100644 --- a/tests/py/unit/core/entity/test_property_defect_analysis.py +++ b/tests/py/unit/core/entity/test_property_defect_analysis.py @@ -9,8 +9,11 @@ STABLE_TO_COLUMN, evaluate_finite_size_fit, fit_finite_size, + flatten_scope_track, get_charge_state_table, + get_chemical_potential_combination, get_formation_energies_vs_fermi_level, + get_formation_energy_at_chemical_potentials, ) # GaAs As-vacancy formation energies at the VBM (eV), 2x2x2 cell, from the QuantumATK tutorial. @@ -58,3 +61,44 @@ def test_fit_finite_size_fits_the_cubic_term_with_four_sizes(): def test_fit_finite_size_needs_two_sizes(): with pytest.raises(ValueError, match="at least two"): fit_finite_size([10.0], [2.3]) + + +# scopeTrack globals of the m-HfO2 jobs rxCNizLKg7hkrPgAh (Zr_Hf) and mpxeDWYNKPBr6zc8Z (V_O). +ZR_HF_SCOPE_TRACK = [ + {"scope": {"global": {"DELTA_N_BY_SYMBOL": {"Hf": -1, "O": 0, "Zr": 1}}}}, + {"scope": {"global": {"DEFECT_FORMATION_ENERGY": 0.3664}}}, +] +V_O_SCOPE_TRACK = [ + {"scope": {"global": {"DELTA_N_BY_SYMBOL": {"O": -1, "Hf": 0}}}}, + {"scope": {"global": {"DEFECT_FORMATION_ENERGY": 6.356}}}, +] +ELEMENTAL = {"O": 0.0, "Hf": 0.0, "Zr": 0.0} +O_RICH = {"O": 0.0, "Hf": -10.6926, "Zr": -10.3396} +O_POOR = {"O": -5.346, "Hf": 0.0} + + +@pytest.mark.parametrize( + "scope_track, delta_mu, expected", + [ + (ZR_HF_SCOPE_TRACK, ELEMENTAL, 0.3664), + (ZR_HF_SCOPE_TRACK, O_RICH, 0.0134), + (V_O_SCOPE_TRACK, O_POOR, 1.010), + (V_O_SCOPE_TRACK, O_RICH, 6.356), # Zr is not an element of the job, so its delta_mu is ignored + ], +) +def test_get_formation_energy_at_chemical_potentials(scope_track, delta_mu, expected): + scope = flatten_scope_track(scope_track) + assert get_formation_energy_at_chemical_potentials(scope, delta_mu) == pytest.approx(expected, abs=1e-4) + + +def test_get_formation_energy_at_chemical_potentials_raises_on_a_missing_element(): + with pytest.raises(KeyError): + get_formation_energy_at_chemical_potentials(flatten_scope_track(ZR_HF_SCOPE_TRACK), O_POOR) + + +@pytest.mark.parametrize( + "delta_n_by_symbol, expected", + [({"O": -1, "Hf": 0}, "Δμ_O"), ({"Zr": 1, "Hf": -1, "O": 0}, "Δμ_Hf − Δμ_Zr"), ({"O": -2}, "2Δμ_O")], +) +def test_get_chemical_potential_combination(delta_n_by_symbol, expected): + assert get_chemical_potential_combination(delta_n_by_symbol) == expected From 07182b0da9caee82c10a67547d3ef72a6db6a1cf Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Tue, 29 Sep 2026 16:44:25 -0700 Subject: [PATCH 16/18] fix(SOF-7975): flat CHEMICAL_POTENTIALS, print the elemental references the job used, draw one line from delta mu = 0 CHEMICAL_POTENTIALS is now a flat {element: delta mu} dict, empty by default, in both notebooks: empty means no correction, and the job's value from the total energies is reported as is. The results cells print the elemental reference energies per atom the job used, mu_i^0, from TE_CONTRIBUTIONS_BY_SYMBOL[element].total_energy_per_atom in the job's scope, which is what the workflow builds SUM_DELTA_N_TIMES_MU from. With delta mu given, they print E_f - sum_i dN_i delta mu_i next to the job's value and draw the line from delta mu = 0 to that point on the existing combination axis. The charged notebook adds the corrected value to results_df as formation_energy_at_chemical_potentials and draws one line per charge for the largest size. The per-condition table, the condition names, the example line and every HfO2 value are gone from the notebooks. The 7.1 markdown states mu_i = mu_i^0 + delta mu_i and where delta mu comes from, and 7.2 refers to delta mu = 0 instead of the removed E_i. The explicit plot mode is dropped, because plotly draws a two-point trace as lines+markers by default. The test constants lose their condition names. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../workflows/defect_formation_energy.ipynb | 24 ++++++------ .../defect_formation_energy_charged.ipynb | 38 +++++++++---------- .../entity/test_property_defect_analysis.py | 15 ++++---- 3 files changed, 36 insertions(+), 41 deletions(-) diff --git a/other/materials_designer/workflows/defect_formation_energy.ipynb b/other/materials_designer/workflows/defect_formation_energy.ipynb index 9ce9273cb..16af1289d 100644 --- a/other/materials_designer/workflows/defect_formation_energy.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy.ipynb @@ -123,8 +123,7 @@ "outputs": [], "source": [ "SCF_KGRID = [4, 4, 4] # same grid for the pristine and defective cells\n", - "CHEMICAL_POTENTIALS = None # Δμ per element, eV/atom, per condition, e.g.\n", - "# {\"O-rich\": {\"O\": 0.0, \"Hf\": -10.693, \"Zr\": -10.340}, \"O-poor\": {\"O\": -5.346, \"Hf\": 0.0, \"Zr\": 0.0}}\n" + "CHEMICAL_POTENTIALS = {} # Δμ per element, eV/atom; empty = no correction, the job's value from the total energies is reported as is\n" ] }, { @@ -560,7 +559,7 @@ "source": [ "## 7. Retrieve results\n", "### 7.1. Retrieve and visualize defect formation energy\n", - "The job reports $E_f$ at the elemental chemical potentials it used, $\\mu_i = E_i$: $E_i$ is the energy per atom of the element's reference (e.g. O₂ for O, the bulk metal for Hf), and $\\Delta N_i$ is the number of atoms of $i$ the defect adds (+) or removes (−). `CHEMICAL_POTENTIALS` sets $\\Delta\\mu_i = \\mu_i - E_i \\le 0$ for every element of the job, per named condition. Obtain them separately from the host's stability range, $\\sum_i n_i\\,\\Delta\\mu_i = \\Delta H_f$ over its formula unit $n_i$, with $\\Delta H_f$ its formation energy (e.g. from [Formation Energies and Convex Hull](analyze_convex_hull.ipynb)). The table lists `formation_energy_at_condition`, $E_f(\\Delta\\mu) = E_f - \\sum_i \\Delta N_i\\,\\Delta\\mu_i$, for each condition, and the plot draws it against $-\\sum_i \\Delta N_i\\,\\Delta\\mu_i$ with the conditions marked." + "The chemical potentials are $\\mu_i = \\mu_i^0 + \\Delta\\mu_i$: $\\mu_i^0$ are the elemental reference energies per atom the job used, printed for reference, and $\\Delta\\mu_i \\le 0$ come from the host's stability range on the phase diagram (formation energies, e.g. [Formation Energies and Convex Hull](analyze_convex_hull.ipynb)), typed in `CHEMICAL_POTENTIALS`. With $\\Delta\\mu$ given, $E_f$ is corrected by $-\\sum_i \\Delta N_i\\,\\Delta\\mu_i$ and drawn from $\\Delta\\mu = 0$ to the given point." ] }, { @@ -570,7 +569,6 @@ "metadata": {}, "outputs": [], "source": [ - "import pandas as pd\n", "import plotly.graph_objects as go\n", "from mat3ra.notebooks_utils.core.entity.property.defect_analysis import (\n", " flatten_scope_track,\n", @@ -583,16 +581,16 @@ "defect_energy_data = client.properties.get_for_job(defect_job_id)\n", "visualize_properties(defect_energy_data, title=\"Defect Formation Energy\")\n", "\n", + "scope = flatten_scope_track(client.jobs.get(defect_job_id)[\"scopeTrack\"])\n", + "for element, contribution in scope[\"TE_CONTRIBUTIONS_BY_SYMBOL\"].items():\n", + " print(f\"μ⁰_{element} = {contribution['total_energy_per_atom']:.4f} eV/atom\")\n", "if CHEMICAL_POTENTIALS:\n", - " scope = flatten_scope_track(client.jobs.get(defect_job_id)[\"scopeTrack\"])\n", - " chemical_potentials_df = pd.DataFrame.from_dict(CHEMICAL_POTENTIALS, orient=\"index\")\n", - " chemical_potentials_df[\"formation_energy_at_condition\"] = [\n", - " get_formation_energy_at_chemical_potentials(scope, delta_mu) for delta_mu in CHEMICAL_POTENTIALS.values()\n", - " ]\n", - " print(chemical_potentials_df)\n", - " energies = chemical_potentials_df[\"formation_energy_at_condition\"]\n", - " figure = go.Figure(go.Scatter(x=energies - scope[\"DEFECT_FORMATION_ENERGY\"], y=energies, text=energies.index,\n", - " mode=\"lines+markers+text\", textposition=\"top center\"))\n", + " formation_energy = scope[\"DEFECT_FORMATION_ENERGY\"]\n", + " formation_energy_at_chemical_potentials = get_formation_energy_at_chemical_potentials(scope, CHEMICAL_POTENTIALS)\n", + " print(f\"E_f = {formation_energy:.4f} eV from the total energies, \"\n", + " f\"{formation_energy_at_chemical_potentials:.4f} eV at Δμ = {CHEMICAL_POTENTIALS}\")\n", + " figure = go.Figure(go.Scatter(x=[0, formation_energy_at_chemical_potentials - formation_energy],\n", + " y=[formation_energy, formation_energy_at_chemical_potentials]))\n", " figure.update_layout(xaxis_title=f\"{get_chemical_potential_combination(scope['DELTA_N_BY_SYMBOL'])} (eV)\",\n", " yaxis_title=\"Defect formation energy (eV)\")\n", " render_figure(figure)" diff --git a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb index ac6eac943..31a2b71e5 100644 --- a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb @@ -133,8 +133,7 @@ "CHARGES = [0] # e.g. [1, 0, -1, -2, -3]\n", "\n", "SCF_KGRID = [4, 4, 4] # same grid for the pristine and defective cells, divided by the supercell scaling\n", - "CHEMICAL_POTENTIALS = None # Δμ per element, eV/atom, per condition, e.g.\n", - "# {\"O-rich\": {\"O\": 0.0, \"Hf\": -10.693, \"Zr\": -10.340}, \"O-poor\": {\"O\": -5.346, \"Hf\": 0.0, \"Zr\": 0.0}}\n", + "CHEMICAL_POTENTIALS = {} # Δμ per element, eV/atom; empty = no correction, the job's value from the total energies is reported as is\n", "\n", "RELAX_PRISTINE_MATERIAL = False\n", "RELAX_DEFECTIVE_MATERIAL = False" @@ -752,7 +751,7 @@ "source": [ "## 7. Retrieve results\n", "### 7.1. Formation energies\n", - "`formation_energy` is $E_f[X^q]$ at $\\mu_e = 0$, as the workflow reports it. `length` is $L = V^{1/3}$ of the supercell, the length scale the image-charge terms in section 7.3 are expressed in. The workflow uses the elemental chemical potentials $\\mu_i = E_i$: $E_i$ is the energy per atom of the element's reference (e.g. O₂ for O, the bulk metal for Hf), and $\\Delta N_i$ is the number of atoms of $i$ the defect adds (+) or removes (−). `CHEMICAL_POTENTIALS` sets $\\Delta\\mu_i = \\mu_i - E_i \\le 0$ for every element of the job, per named condition. Obtain them separately from the host's stability range, $\\sum_i n_i\\,\\Delta\\mu_i = \\Delta H_f$ over its formula unit $n_i$, with $\\Delta H_f$ its formation energy (e.g. from [Formation Energies and Convex Hull](analyze_convex_hull.ipynb)). A table then lists `formation_energy_at_condition`, $E_f(\\Delta\\mu) = E_f - \\sum_i \\Delta N_i\\,\\Delta\\mu_i$, per size, charge and condition, and a plot draws it against $-\\sum_i \\Delta N_i\\,\\Delta\\mu_i$ for the largest size, one line per charge state, with the conditions marked." + "`formation_energy` is $E_f[X^q]$ at $\\mu_e = 0$, as the workflow reports it. `length` is $L = V^{1/3}$ of the supercell, the length scale the image-charge terms in section 7.3 are expressed in. The chemical potentials are $\\mu_i = \\mu_i^0 + \\Delta\\mu_i$: $\\mu_i^0$ are the elemental reference energies per atom the job used, printed for reference, and $\\Delta\\mu_i \\le 0$ come from the host's stability range on the phase diagram (formation energies, e.g. [Formation Energies and Convex Hull](analyze_convex_hull.ipynb)), typed in `CHEMICAL_POTENTIALS`. With $\\Delta\\mu$ given, $E_f$ is corrected by $-\\sum_i \\Delta N_i\\,\\Delta\\mu_i$, with $\\Delta N_i$ the atoms of $i$ the defect adds (+) or removes (−), listed as `formation_energy_at_chemical_potentials` and drawn from $\\Delta\\mu = 0$ for the largest size, one line per charge state." ] }, { @@ -770,7 +769,7 @@ ")\n", "from mat3ra.notebooks_utils.ipython.plot._plotly import render_figure\n", "\n", - "results_records, chemical_potential_records = [], []\n", + "results_records = []\n", "for record in job_records:\n", " job = client.jobs.get(record[\"job_id\"])\n", " properties = client.properties.get_for_job(record[\"job_id\"], property_name=\"defect_formation_energy\")\n", @@ -779,29 +778,28 @@ " \"final_status\": job.get(\"status\"),\n", " \"formation_energy\": properties[0].get(\"value\") if properties else None,\n", " })\n", - " if CHEMICAL_POTENTIALS and properties:\n", + " if properties:\n", " scope = flatten_scope_track(job[\"scopeTrack\"])\n", - " for condition, delta_mu in CHEMICAL_POTENTIALS.items():\n", - " chemical_potential_records.append({\n", - " \"scaling\": record[\"scaling\"], \"charge\": record[\"charge\"], \"condition\": condition, **delta_mu,\n", - " \"formation_energy\": scope[\"DEFECT_FORMATION_ENERGY\"],\n", - " \"formation_energy_at_condition\": get_formation_energy_at_chemical_potentials(scope, delta_mu),\n", - " })\n", + " if CHEMICAL_POTENTIALS:\n", + " results_records[-1][\"formation_energy_at_chemical_potentials\"] = (\n", + " get_formation_energy_at_chemical_potentials(scope, CHEMICAL_POTENTIALS)\n", + " )\n", "\n", "results_df = pd.DataFrame(results_records)\n", "successful_df = results_df[(results_df[\"final_status\"] == \"finished\") & results_df[\"formation_energy\"].notna()]\n", "if len(successful_df) < len(results_df):\n", " print(f\"⚠️ {len(results_df) - len(successful_df)} of {len(results_df)} job(s) returned no formation energy; \"\n", " \"they are left out of sections 7.2 and 7.3.\")\n", - "if chemical_potential_records:\n", - " chemical_potentials_df = pd.DataFrame(chemical_potential_records)\n", - " print(chemical_potentials_df.to_string(index=False))\n", - " largest_df = chemical_potentials_df[chemical_potentials_df[\"scaling\"] == chemical_potentials_df[\"scaling\"].max()]\n", + "if not successful_df.empty:\n", + " for element, contribution in scope[\"TE_CONTRIBUTIONS_BY_SYMBOL\"].items():\n", + " print(f\"μ⁰_{element} = {contribution['total_energy_per_atom']:.4f} eV/atom\")\n", + "if CHEMICAL_POTENTIALS and not successful_df.empty:\n", + " largest_df = successful_df[successful_df[\"scaling\"] == successful_df[\"scaling\"].max()]\n", " figure = go.Figure()\n", - " for charge, group in largest_df.groupby(\"charge\"):\n", - " figure.add_scatter(x=group[\"formation_energy_at_condition\"] - group[\"formation_energy\"],\n", - " y=group[\"formation_energy_at_condition\"], text=group[\"condition\"],\n", - " mode=\"lines+markers+text\", textposition=\"top center\", name=f\"q = {charge:+d}\")\n", + " for row in largest_df.itertuples():\n", + " figure.add_scatter(x=[0, row.formation_energy_at_chemical_potentials - row.formation_energy],\n", + " y=[row.formation_energy, row.formation_energy_at_chemical_potentials],\n", + " name=f\"q = {row.charge:+d}\")\n", " figure.update_layout(title=f\"Defect formation energy (n={largest_df['scaling'].iloc[0]})\",\n", " xaxis_title=f\"{get_chemical_potential_combination(scope['DELTA_N_BY_SYMBOL'])} (eV)\",\n", " yaxis_title=\"Formation energy at the VBM (eV)\")\n", @@ -815,7 +813,7 @@ "metadata": {}, "source": [ "### 7.2. Formation energy versus the electron chemical potential\n", - "Each charge state gives a straight line of slope $q$ over $\\mu_e \\in [0, E_{\\text{gap}}]$; the lower envelope is the charge state the defect actually adopts, and the crossings between the lines are the charge transition levels. The largest supercell with a result for every charge state is used, with $E_{\\text{gap}}$ from its own Band Gap job. The table lists each charge state's formation energy at the VBM and the range of the Fermi level $\\mu_e$ in which it is the stable state. 7.2 and 7.3 use the stored $E_f$ ($\\mu_i = E_i$). A condition shifts every charge state by the same $-\\sum_i \\Delta N_i\\,\\Delta\\mu_i$, so the transition levels do not depend on it." + "Each charge state gives a straight line of slope $q$ over $\\mu_e \\in [0, E_{\\text{gap}}]$; the lower envelope is the charge state the defect actually adopts, and the crossings between the lines are the charge transition levels. The largest supercell with a result for every charge state is used, with $E_{\\text{gap}}$ from its own Band Gap job. The table lists each charge state's formation energy at the VBM and the range of the Fermi level $\\mu_e$ in which it is the stable state. 7.2 and 7.3 use the job's $E_f$ ($\\Delta\\mu = 0$); $\\Delta\\mu$ shifts every charge state by the same $-\\sum_i \\Delta N_i\\,\\Delta\\mu_i$, so the transition levels do not depend on it." ] }, { diff --git a/tests/py/unit/core/entity/test_property_defect_analysis.py b/tests/py/unit/core/entity/test_property_defect_analysis.py index 2b63b7cf0..93ff292fd 100644 --- a/tests/py/unit/core/entity/test_property_defect_analysis.py +++ b/tests/py/unit/core/entity/test_property_defect_analysis.py @@ -72,18 +72,17 @@ def test_fit_finite_size_needs_two_sizes(): {"scope": {"global": {"DELTA_N_BY_SYMBOL": {"O": -1, "Hf": 0}}}}, {"scope": {"global": {"DEFECT_FORMATION_ENERGY": 6.356}}}, ] -ELEMENTAL = {"O": 0.0, "Hf": 0.0, "Zr": 0.0} -O_RICH = {"O": 0.0, "Hf": -10.6926, "Zr": -10.3396} -O_POOR = {"O": -5.346, "Hf": 0.0} +ZR_HF_DELTA_MU = {"O": 0.0, "Hf": -10.693, "Zr": -10.340} +V_O_DELTA_MU = {"O": -5.346, "Hf": 0.0} @pytest.mark.parametrize( "scope_track, delta_mu, expected", [ - (ZR_HF_SCOPE_TRACK, ELEMENTAL, 0.3664), - (ZR_HF_SCOPE_TRACK, O_RICH, 0.0134), - (V_O_SCOPE_TRACK, O_POOR, 1.010), - (V_O_SCOPE_TRACK, O_RICH, 6.356), # Zr is not an element of the job, so its delta_mu is ignored + (ZR_HF_SCOPE_TRACK, {"O": 0.0, "Hf": 0.0, "Zr": 0.0}, 0.3664), + (ZR_HF_SCOPE_TRACK, ZR_HF_DELTA_MU, 0.0134), + (V_O_SCOPE_TRACK, V_O_DELTA_MU, 1.010), + (V_O_SCOPE_TRACK, ZR_HF_DELTA_MU, 6.356), # Zr is not an element of the job, so its delta_mu is ignored ], ) def test_get_formation_energy_at_chemical_potentials(scope_track, delta_mu, expected): @@ -93,7 +92,7 @@ def test_get_formation_energy_at_chemical_potentials(scope_track, delta_mu, expe def test_get_formation_energy_at_chemical_potentials_raises_on_a_missing_element(): with pytest.raises(KeyError): - get_formation_energy_at_chemical_potentials(flatten_scope_track(ZR_HF_SCOPE_TRACK), O_POOR) + get_formation_energy_at_chemical_potentials(flatten_scope_track(ZR_HF_SCOPE_TRACK), V_O_DELTA_MU) @pytest.mark.parametrize( From 0868e37006c90a4819025062a28d5e7f295031a0 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Tue, 29 Sep 2026 17:06:49 -0700 Subject: [PATCH 17/18] fix(SOF-7975): the charged relaxation as a function next to the defect workflow, the job result read in the library The charged notebook's 6.2 loop no longer carries the fixed-cell relaxation inline: relax_defective_supercell(defective_supercell, scaling, charge), defined in 5.2 next to create_defect_workflow, builds the relaxation workflow for the supercell's k-grid with tot_charge for q != 0, reuses a finished or running relaxation of that name and grid or creates and submits one tagged charge:q, waits, and returns the final structure. The loop body is now: relax if RELAX_DEFECTIVE_MATERIAL, create the defect workflow, create the job, record it. Jobs, prints and their order are unchanged. The 7.1 cells no longer read the job's scopeTrack. get_defect_job_result(api_client, job_id) returns a DefectJobResult (formation energy, dN per element, the elemental reference energies per atom the job used), and get_formation_energy_at_chemical_potentials takes that result. plot_formation_energy_vs_chemical_potentials in defect_plot.py draws each line from (0, E_f) to (x, E_f at delta mu); the neutral notebook passes one line named after DEFECTIVE_NAME, the charged one a line per charge for the largest size. flatten_scope_track is private. The CHEMICAL_POTENTIALS comment gives an example, and the 7.1 markdown says what the cell shows. The y axis of the neutral plot now reads "Formation energy at the VBM (eV)" and the plot has a title, as the charged one does. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../workflows/defect_formation_energy.ipynb | 27 +++---- .../defect_formation_energy_charged.ipynb | 78 ++++++++++--------- .../core/entity/property/defect_analysis.py | 33 ++++++-- .../ipython/entity/property/defect_plot.py | 17 +++- .../entity/test_property_defect_analysis.py | 59 +++++++++----- 5 files changed, 135 insertions(+), 79 deletions(-) diff --git a/other/materials_designer/workflows/defect_formation_energy.ipynb b/other/materials_designer/workflows/defect_formation_energy.ipynb index 16af1289d..b0eef9631 100644 --- a/other/materials_designer/workflows/defect_formation_energy.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy.ipynb @@ -123,7 +123,7 @@ "outputs": [], "source": [ "SCF_KGRID = [4, 4, 4] # same grid for the pristine and defective cells\n", - "CHEMICAL_POTENTIALS = {} # Δμ per element, eV/atom; empty = no correction, the job's value from the total energies is reported as is\n" + "CHEMICAL_POTENTIALS = {} # Δμ per element, eV/atom, e.g. {\"O\": -0.12}\n" ] }, { @@ -559,7 +559,7 @@ "source": [ "## 7. Retrieve results\n", "### 7.1. Retrieve and visualize defect formation energy\n", - "The chemical potentials are $\\mu_i = \\mu_i^0 + \\Delta\\mu_i$: $\\mu_i^0$ are the elemental reference energies per atom the job used, printed for reference, and $\\Delta\\mu_i \\le 0$ come from the host's stability range on the phase diagram (formation energies, e.g. [Formation Energies and Convex Hull](analyze_convex_hull.ipynb)), typed in `CHEMICAL_POTENTIALS`. With $\\Delta\\mu$ given, $E_f$ is corrected by $-\\sum_i \\Delta N_i\\,\\Delta\\mu_i$ and drawn from $\\Delta\\mu = 0$ to the given point." + "The elemental reference energies per atom $\\mu_i^0$ that the job used are printed below. With $\\Delta\\mu_i$ in `CHEMICAL_POTENTIALS` (e.g. from the host's stability range in [Formation Energies and Convex Hull](analyze_convex_hull.ipynb)), $E_f$ is corrected by $-\\sum_i \\Delta N_i\\,\\Delta\\mu_i$ and plotted from $\\Delta\\mu = 0$." ] }, { @@ -569,31 +569,28 @@ "metadata": {}, "outputs": [], "source": [ - "import plotly.graph_objects as go\n", "from mat3ra.notebooks_utils.core.entity.property.defect_analysis import (\n", - " flatten_scope_track,\n", " get_chemical_potential_combination,\n", + " get_defect_job_result,\n", " get_formation_energy_at_chemical_potentials,\n", ")\n", + "from mat3ra.notebooks_utils.ipython.entity.property.defect_plot import plot_formation_energy_vs_chemical_potentials\n", "from mat3ra.notebooks_utils.ipython.entity.property.visualize import visualize_properties\n", "from mat3ra.notebooks_utils.ipython.plot._plotly import render_figure\n", "\n", "defect_energy_data = client.properties.get_for_job(defect_job_id)\n", "visualize_properties(defect_energy_data, title=\"Defect Formation Energy\")\n", "\n", - "scope = flatten_scope_track(client.jobs.get(defect_job_id)[\"scopeTrack\"])\n", - "for element, contribution in scope[\"TE_CONTRIBUTIONS_BY_SYMBOL\"].items():\n", - " print(f\"μ⁰_{element} = {contribution['total_energy_per_atom']:.4f} eV/atom\")\n", + "result = get_defect_job_result(client, defect_job_id)\n", + "for element, energy in result.reference_energies_per_atom.items():\n", + " print(f\"μ⁰_{element} = {energy:.4f} eV/atom\")\n", "if CHEMICAL_POTENTIALS:\n", - " formation_energy = scope[\"DEFECT_FORMATION_ENERGY\"]\n", - " formation_energy_at_chemical_potentials = get_formation_energy_at_chemical_potentials(scope, CHEMICAL_POTENTIALS)\n", - " print(f\"E_f = {formation_energy:.4f} eV from the total energies, \"\n", + " formation_energy_at_chemical_potentials = get_formation_energy_at_chemical_potentials(result, CHEMICAL_POTENTIALS)\n", + " print(f\"E_f = {result.formation_energy:.4f} eV from the total energies, \"\n", " f\"{formation_energy_at_chemical_potentials:.4f} eV at Δμ = {CHEMICAL_POTENTIALS}\")\n", - " figure = go.Figure(go.Scatter(x=[0, formation_energy_at_chemical_potentials - formation_energy],\n", - " y=[formation_energy, formation_energy_at_chemical_potentials]))\n", - " figure.update_layout(xaxis_title=f\"{get_chemical_potential_combination(scope['DELTA_N_BY_SYMBOL'])} (eV)\",\n", - " yaxis_title=\"Defect formation energy (eV)\")\n", - " render_figure(figure)" + " lines = {DEFECTIVE_NAME: (result.formation_energy, formation_energy_at_chemical_potentials)}\n", + " x_label = f\"{get_chemical_potential_combination(result.delta_n_by_symbol)} (eV)\"\n", + " render_figure(plot_formation_energy_vs_chemical_potentials(lines, x_label, \"Defect formation energy\"))" ] } ], diff --git a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb index 31a2b71e5..01e7f2184 100644 --- a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb @@ -133,7 +133,7 @@ "CHARGES = [0] # e.g. [1, 0, -1, -2, -3]\n", "\n", "SCF_KGRID = [4, 4, 4] # same grid for the pristine and defective cells, divided by the supercell scaling\n", - "CHEMICAL_POTENTIALS = {} # Δμ per element, eV/atom; empty = no correction, the job's value from the total energies is reported as is\n", + "CHEMICAL_POTENTIALS = {} # Δμ per element, eV/atom, e.g. {\"O\": -0.12}\n", "\n", "RELAX_PRISTINE_MATERIAL = False\n", "RELAX_DEFECTIVE_MATERIAL = False" @@ -543,12 +543,16 @@ "source": [ "from mat3ra.standata.workflows import WorkflowStandata\n", "from mat3ra.wode.workflows import Workflow\n", + "from mat3ra.notebooks_utils.core.entity.job.api import find_job_for_material\n", "from mat3ra.notebooks_utils.ipython.entity.workflow.visualize import visualize_workflow\n", "from mat3ra.notebooks_utils.workflow import apply_scf_kgrid, patch_workflow_qe_input, set_assignment_value\n", "\n", "defect_workflow_config = WorkflowStandata.filter_by_application(app.name).get_by_name_first_match(\n", " WORKFLOW_SEARCH_TERM\n", ")\n", + "defect_relax_workflow_config = WorkflowStandata.filter_by_application(app.name).get_by_name_first_match(\n", + " \"fixed_cell_relaxation.json\"\n", + ")\n", "\n", "REFERENCE_JOB_FILTER_UNITS = {\n", " \"total_energy\": \"assign-reference-job-filter\",\n", @@ -574,6 +578,28 @@ " return apply_scf_kgrid(workflow, get_scf_kgrid_for_supercell(scaling), material=defective_supercell)\n", "\n", "\n", + "async def relax_defective_supercell(defective_supercell, scaling, charge):\n", + " \"\"\"The defective supercell relaxed at fixed cell in its charge state; a finished or running relaxation is reused.\"\"\"\n", + " kgrid = get_scf_kgrid_for_supercell(scaling)\n", + " relax_workflow = Workflow.create(defect_relax_workflow_config)\n", + " relax_workflow.name = f\"{relax_workflow.name} {defective_supercell.name} q={charge:+d}\"\n", + " apply_scf_kgrid(relax_workflow, kgrid, material=defective_supercell, unit_name=\"pw_relax\")\n", + " if charge:\n", + " patch_workflow_qe_input(relax_workflow, {\"system\": {\"tot_charge\": charge}}, unit_names=[\"pw_relax\"])\n", + " relax_job = find_job_for_material(\n", + " client, defective_supercell.id, relax_workflow.name, ACCOUNT_ID, kgrid=kgrid, unit_name=\"pw_relax\",\n", + " statuses=(\"finished\", \"active\", \"submitted\", \"queued\"),\n", + " )\n", + " if relax_job is None:\n", + " relax_job = create_job_for_materials([defective_supercell], relax_workflow, [f\"charge:{charge}\"])\n", + " submit_jobs(client.jobs, [relax_job[\"_id\"]])\n", + " print(f\"✅ n={scaling} q={charge:+d}: created a relaxation job {relax_job['_id']}\")\n", + " else:\n", + " print(f\"♻️ n={scaling} q={charge:+d}: reusing the relaxation from job {relax_job['_id']}\")\n", + " await wait_for_jobs_to_finish_async(client.jobs, [relax_job[\"_id\"]], poll_interval=POLL_INTERVAL)\n", + " return get_final_structure_for_job(client, relax_job[\"_id\"])\n", + "\n", + "\n", "visualize_workflow(Workflow.create(defect_workflow_config))" ] }, @@ -669,34 +695,14 @@ "outputs": [], "source": [ "import pandas as pd\n", - "from mat3ra.notebooks_utils.core.entity.job.api import find_job_for_material\n", "\n", - "defect_relax_workflow_config = WorkflowStandata.filter_by_application(app.name).get_by_name_first_match(\n", - " \"fixed_cell_relaxation.json\"\n", - ")\n", "job_records = []\n", "for scaling, (pristine_supercell, defective_supercell) in saved_supercell_pairs.items():\n", " kgrid = get_scf_kgrid_for_supercell(scaling)\n", " for charge in CHARGES:\n", " defective_material = defective_supercell\n", " if RELAX_DEFECTIVE_MATERIAL:\n", - " relax_workflow = Workflow.create(defect_relax_workflow_config)\n", - " relax_workflow.name = f\"{relax_workflow.name} {defective_supercell.name} q={charge:+d}\"\n", - " apply_scf_kgrid(relax_workflow, kgrid, material=defective_supercell, unit_name=\"pw_relax\")\n", - " if charge:\n", - " patch_workflow_qe_input(relax_workflow, {\"system\": {\"tot_charge\": charge}}, unit_names=[\"pw_relax\"])\n", - " relax_job = find_job_for_material(\n", - " client, defective_supercell.id, relax_workflow.name, ACCOUNT_ID, kgrid=kgrid, unit_name=\"pw_relax\",\n", - " statuses=(\"finished\", \"active\", \"submitted\", \"queued\"),\n", - " )\n", - " if relax_job is None:\n", - " relax_job = create_job_for_materials([defective_supercell], relax_workflow, [f\"charge:{charge}\"])\n", - " submit_jobs(client.jobs, [relax_job[\"_id\"]])\n", - " print(f\"✅ n={scaling} q={charge:+d}: created a relaxation job {relax_job['_id']}\")\n", - " else:\n", - " print(f\"♻️ n={scaling} q={charge:+d}: reusing the relaxation from job {relax_job['_id']}\")\n", - " await wait_for_jobs_to_finish_async(client.jobs, [relax_job[\"_id\"]], poll_interval=POLL_INTERVAL)\n", - " defective_material = get_final_structure_for_job(client, relax_job[\"_id\"])\n", + " defective_material = await relax_defective_supercell(defective_supercell, scaling, charge)\n", " workflow = create_defect_workflow(scaling, charge, pristine_reference_jobs[scaling])\n", " # Order matters: [0] defective (computed), [1] pristine (reference).\n", " job = create_job_for_materials([defective_material, pristine_supercell], workflow, [f\"charge:{charge}\"])\n", @@ -751,7 +757,7 @@ "source": [ "## 7. Retrieve results\n", "### 7.1. Formation energies\n", - "`formation_energy` is $E_f[X^q]$ at $\\mu_e = 0$, as the workflow reports it. `length` is $L = V^{1/3}$ of the supercell, the length scale the image-charge terms in section 7.3 are expressed in. The chemical potentials are $\\mu_i = \\mu_i^0 + \\Delta\\mu_i$: $\\mu_i^0$ are the elemental reference energies per atom the job used, printed for reference, and $\\Delta\\mu_i \\le 0$ come from the host's stability range on the phase diagram (formation energies, e.g. [Formation Energies and Convex Hull](analyze_convex_hull.ipynb)), typed in `CHEMICAL_POTENTIALS`. With $\\Delta\\mu$ given, $E_f$ is corrected by $-\\sum_i \\Delta N_i\\,\\Delta\\mu_i$, with $\\Delta N_i$ the atoms of $i$ the defect adds (+) or removes (−), listed as `formation_energy_at_chemical_potentials` and drawn from $\\Delta\\mu = 0$ for the largest size, one line per charge state." + "`formation_energy` is $E_f[X^q]$ at $\\mu_e = 0$, as the workflow reports it. `length` is $L = V^{1/3}$ of the supercell, the length scale the image-charge terms in section 7.3 are expressed in. The elemental reference energies per atom $\\mu_i^0$ that the jobs used are printed below. With $\\Delta\\mu_i$ in `CHEMICAL_POTENTIALS` (e.g. from the host's stability range in [Formation Energies and Convex Hull](analyze_convex_hull.ipynb)), $E_f$ is corrected by $-\\sum_i \\Delta N_i\\,\\Delta\\mu_i$, with $\\Delta N_i$ the atoms of $i$ the defect adds (+) or removes (−), in `formation_energy_at_chemical_potentials` and plotted from $\\Delta\\mu = 0$ for the largest size, one line per charge state." ] }, { @@ -761,12 +767,12 @@ "metadata": {}, "outputs": [], "source": [ - "import plotly.graph_objects as go\n", "from mat3ra.notebooks_utils.core.entity.property.defect_analysis import (\n", - " flatten_scope_track,\n", " get_chemical_potential_combination,\n", + " get_defect_job_result,\n", " get_formation_energy_at_chemical_potentials,\n", ")\n", + "from mat3ra.notebooks_utils.ipython.entity.property.defect_plot import plot_formation_energy_vs_chemical_potentials\n", "from mat3ra.notebooks_utils.ipython.plot._plotly import render_figure\n", "\n", "results_records = []\n", @@ -779,10 +785,10 @@ " \"formation_energy\": properties[0].get(\"value\") if properties else None,\n", " })\n", " if properties:\n", - " scope = flatten_scope_track(job[\"scopeTrack\"])\n", + " result = get_defect_job_result(client, record[\"job_id\"])\n", " if CHEMICAL_POTENTIALS:\n", " results_records[-1][\"formation_energy_at_chemical_potentials\"] = (\n", - " get_formation_energy_at_chemical_potentials(scope, CHEMICAL_POTENTIALS)\n", + " get_formation_energy_at_chemical_potentials(result, CHEMICAL_POTENTIALS)\n", " )\n", "\n", "results_df = pd.DataFrame(results_records)\n", @@ -791,19 +797,15 @@ " print(f\"⚠️ {len(results_df) - len(successful_df)} of {len(results_df)} job(s) returned no formation energy; \"\n", " \"they are left out of sections 7.2 and 7.3.\")\n", "if not successful_df.empty:\n", - " for element, contribution in scope[\"TE_CONTRIBUTIONS_BY_SYMBOL\"].items():\n", - " print(f\"μ⁰_{element} = {contribution['total_energy_per_atom']:.4f} eV/atom\")\n", + " for element, energy in result.reference_energies_per_atom.items():\n", + " print(f\"μ⁰_{element} = {energy:.4f} eV/atom\")\n", "if CHEMICAL_POTENTIALS and not successful_df.empty:\n", " largest_df = successful_df[successful_df[\"scaling\"] == successful_df[\"scaling\"].max()]\n", - " figure = go.Figure()\n", - " for row in largest_df.itertuples():\n", - " figure.add_scatter(x=[0, row.formation_energy_at_chemical_potentials - row.formation_energy],\n", - " y=[row.formation_energy, row.formation_energy_at_chemical_potentials],\n", - " name=f\"q = {row.charge:+d}\")\n", - " figure.update_layout(title=f\"Defect formation energy (n={largest_df['scaling'].iloc[0]})\",\n", - " xaxis_title=f\"{get_chemical_potential_combination(scope['DELTA_N_BY_SYMBOL'])} (eV)\",\n", - " yaxis_title=\"Formation energy at the VBM (eV)\")\n", - " render_figure(figure)\n", + " lines = {f\"q = {row.charge:+d}\": (row.formation_energy, row.formation_energy_at_chemical_potentials)\n", + " for row in largest_df.itertuples()}\n", + " x_label = f\"{get_chemical_potential_combination(result.delta_n_by_symbol)} (eV)\"\n", + " title = f\"Defect formation energy (n={largest_df['scaling'].iloc[0]})\"\n", + " render_figure(plot_formation_energy_vs_chemical_potentials(lines, x_label, title))\n", "results_df" ] }, diff --git a/src/py/mat3ra/notebooks_utils/core/entity/property/defect_analysis.py b/src/py/mat3ra/notebooks_utils/core/entity/property/defect_analysis.py index 09583f628..e457e5087 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/property/defect_analysis.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/property/defect_analysis.py @@ -1,13 +1,13 @@ """Formation energies of charged defects: dependence on the Fermi level, on the supercell size and on the chemical potentials. -Pure computation — no API calls, no display. Kept apart from `analysis.py`, which imports pymatgen's +No display. Kept apart from `analysis.py`, which imports pymatgen's phase diagram module at load time; that module needs tqdm, which JupyterLite does not install. The Fermi level is the electron chemical potential mu_e, measured from the valence band maximum (VBM). """ -from typing import Any, Dict, List, Sequence +from typing import Any, Dict, List, NamedTuple, Sequence import numpy as np import pandas as pd @@ -111,14 +111,33 @@ def evaluate_finite_size_fit(fit: Dict[str, float], inverse_lengths: Sequence[fl return fit["E_inf"] + fit["a"] * inverse_lengths_array + fit.get("b", 0.0) * inverse_lengths_array**3 -def flatten_scope_track(scope_track: List[dict]) -> Dict[str, Any]: +class DefectJobResult(NamedTuple): + formation_energy: float + delta_n_by_symbol: Dict[str, int] + reference_energies_per_atom: Dict[str, float] + + +def _flatten_scope_track(scope_track: List[dict]) -> Dict[str, Any]: """The global scope of a job's `scopeTrack`, later values overriding earlier ones.""" return {key: value for item in scope_track for key, value in item["scope"]["global"].items()} -def get_formation_energy_at_chemical_potentials(scope: Dict[str, Any], delta_mu: Dict[str, float]) -> float: +def get_defect_job_result(api_client, job_id: str) -> DefectJobResult: + """The result of a finished Defect Formation Energy job, read from its scope.""" + scope = _flatten_scope_track(api_client.jobs.get(job_id)["scopeTrack"]) + return DefectJobResult( + formation_energy=scope["DEFECT_FORMATION_ENERGY"], + delta_n_by_symbol=scope["DELTA_N_BY_SYMBOL"], + reference_energies_per_atom={ + element: contribution["total_energy_per_atom"] + for element, contribution in scope["TE_CONTRIBUTIONS_BY_SYMBOL"].items() + }, + ) + + +def get_formation_energy_at_chemical_potentials(result: DefectJobResult, delta_mu: Dict[str, float]) -> float: """ - Defect formation energy at the chemical potentials mu_i = E_i + delta_mu[i], from its job's flattened scope. + Defect formation energy at the chemical potentials mu_i = E_i + delta_mu[i]. E_i is the elemental energy per atom the job used, so the stored value is the one at delta_mu = 0: E_f(delta_mu) = E_f - sum_i dN_i * delta_mu[i]. @@ -126,8 +145,8 @@ def get_formation_energy_at_chemical_potentials(scope: Dict[str, Any], delta_mu: Raises: KeyError: If `delta_mu` has no value for an element of the job, including one with dN_i = 0. """ - return scope["DEFECT_FORMATION_ENERGY"] - sum( - count * delta_mu[element] for element, count in scope["DELTA_N_BY_SYMBOL"].items() + return result.formation_energy - sum( + count * delta_mu[element] for element, count in result.delta_n_by_symbol.items() ) diff --git a/src/py/mat3ra/notebooks_utils/ipython/entity/property/defect_plot.py b/src/py/mat3ra/notebooks_utils/ipython/entity/property/defect_plot.py index e57eaf9ae..d3b3f2ef1 100644 --- a/src/py/mat3ra/notebooks_utils/ipython/entity/property/defect_plot.py +++ b/src/py/mat3ra/notebooks_utils/ipython/entity/property/defect_plot.py @@ -4,7 +4,7 @@ needs tqdm, which JupyterLite does not install. """ -from typing import Dict +from typing import Dict, Tuple import numpy as np import pandas as pd @@ -67,3 +67,18 @@ def plot_finite_size_fits(results: pd.DataFrame, fits: Dict[int, Dict[str, float yaxis_title="Formation energy at the VBM (eV)", ) return figure + + +def plot_formation_energy_vs_chemical_potentials( + lines: Dict[str, Tuple[float, float]], x_label: str, title: str +) -> go.Figure: + """One line per name, from (0, E_f at delta_mu = 0) to (x, E_f at delta_mu), with x their difference.""" + figure = go.Figure() + for name, (formation_energy, formation_energy_at_chemical_potentials) in lines.items(): + figure.add_scatter( + x=[0, formation_energy_at_chemical_potentials - formation_energy], + y=[formation_energy, formation_energy_at_chemical_potentials], + name=name, + ) + figure.update_layout(title=title, xaxis_title=x_label, yaxis_title="Formation energy at the VBM (eV)") + return figure diff --git a/tests/py/unit/core/entity/test_property_defect_analysis.py b/tests/py/unit/core/entity/test_property_defect_analysis.py index 93ff292fd..3ae40a13d 100644 --- a/tests/py/unit/core/entity/test_property_defect_analysis.py +++ b/tests/py/unit/core/entity/test_property_defect_analysis.py @@ -1,5 +1,7 @@ """Unit tests for charged-defect formation energy analysis.""" +from unittest.mock import MagicMock + import numpy as np import pytest from mat3ra.notebooks_utils.core.entity.property.defect_analysis import ( @@ -7,11 +9,12 @@ FORMATION_ENERGY_AT_VBM_COLUMN, STABLE_FROM_COLUMN, STABLE_TO_COLUMN, + DefectJobResult, evaluate_finite_size_fit, fit_finite_size, - flatten_scope_track, get_charge_state_table, get_chemical_potential_combination, + get_defect_job_result, get_formation_energies_vs_fermi_level, get_formation_energy_at_chemical_potentials, ) @@ -64,35 +67,55 @@ def test_fit_finite_size_needs_two_sizes(): # scopeTrack globals of the m-HfO2 jobs rxCNizLKg7hkrPgAh (Zr_Hf) and mpxeDWYNKPBr6zc8Z (V_O). -ZR_HF_SCOPE_TRACK = [ - {"scope": {"global": {"DELTA_N_BY_SYMBOL": {"Hf": -1, "O": 0, "Zr": 1}}}}, - {"scope": {"global": {"DEFECT_FORMATION_ENERGY": 0.3664}}}, -] -V_O_SCOPE_TRACK = [ - {"scope": {"global": {"DELTA_N_BY_SYMBOL": {"O": -1, "Hf": 0}}}}, - {"scope": {"global": {"DEFECT_FORMATION_ENERGY": 6.356}}}, -] +HAFNIUM = {"total_energy_per_atom": -2160.2365} +OXYGEN = {"total_energy_per_atom": -437.4523} +ZIRCONIUM = {"total_energy_per_atom": -1349.0462} +ZR_HF_JOB = { + "scopeTrack": [ + {"scope": {"global": {"DELTA_N_BY_SYMBOL": {"Hf": -1, "O": 0, "Zr": 1}}}}, + {"scope": {"global": {"TE_CONTRIBUTIONS_BY_SYMBOL": {"Hf": HAFNIUM, "O": OXYGEN, "Zr": ZIRCONIUM}}}}, + {"scope": {"global": {"DEFECT_FORMATION_ENERGY": 0.3664}}}, + ] +} +V_O_JOB = { + "scopeTrack": [ + {"scope": {"global": {"DELTA_N_BY_SYMBOL": {"O": -1, "Hf": 0}}}}, + {"scope": {"global": {"TE_CONTRIBUTIONS_BY_SYMBOL": {"Hf": HAFNIUM, "O": OXYGEN}}}}, + {"scope": {"global": {"DEFECT_FORMATION_ENERGY": 6.356}}}, + ] +} +ZR_HF_RESULT = DefectJobResult( + 0.3664, {"Hf": -1, "O": 0, "Zr": 1}, {"Hf": -2160.2365, "O": -437.4523, "Zr": -1349.0462} +) +V_O_RESULT = DefectJobResult(6.356, {"O": -1, "Hf": 0}, {"Hf": -2160.2365, "O": -437.4523}) ZR_HF_DELTA_MU = {"O": 0.0, "Hf": -10.693, "Zr": -10.340} V_O_DELTA_MU = {"O": -5.346, "Hf": 0.0} +@pytest.mark.parametrize("job, expected", [(ZR_HF_JOB, ZR_HF_RESULT), (V_O_JOB, V_O_RESULT)]) +def test_get_defect_job_result(job, expected): + api_client = MagicMock() + api_client.jobs.get.return_value = job + assert get_defect_job_result(api_client, "job-1") == expected + api_client.jobs.get.assert_called_once_with("job-1") + + @pytest.mark.parametrize( - "scope_track, delta_mu, expected", + "result, delta_mu, expected", [ - (ZR_HF_SCOPE_TRACK, {"O": 0.0, "Hf": 0.0, "Zr": 0.0}, 0.3664), - (ZR_HF_SCOPE_TRACK, ZR_HF_DELTA_MU, 0.0134), - (V_O_SCOPE_TRACK, V_O_DELTA_MU, 1.010), - (V_O_SCOPE_TRACK, ZR_HF_DELTA_MU, 6.356), # Zr is not an element of the job, so its delta_mu is ignored + (ZR_HF_RESULT, {"O": 0.0, "Hf": 0.0, "Zr": 0.0}, 0.3664), + (ZR_HF_RESULT, ZR_HF_DELTA_MU, 0.0134), + (V_O_RESULT, V_O_DELTA_MU, 1.010), + (V_O_RESULT, ZR_HF_DELTA_MU, 6.356), # Zr is not an element of the job, so its delta_mu is ignored ], ) -def test_get_formation_energy_at_chemical_potentials(scope_track, delta_mu, expected): - scope = flatten_scope_track(scope_track) - assert get_formation_energy_at_chemical_potentials(scope, delta_mu) == pytest.approx(expected, abs=1e-4) +def test_get_formation_energy_at_chemical_potentials(result, delta_mu, expected): + assert get_formation_energy_at_chemical_potentials(result, delta_mu) == pytest.approx(expected, abs=1e-4) def test_get_formation_energy_at_chemical_potentials_raises_on_a_missing_element(): with pytest.raises(KeyError): - get_formation_energy_at_chemical_potentials(flatten_scope_track(ZR_HF_SCOPE_TRACK), V_O_DELTA_MU) + get_formation_energy_at_chemical_potentials(ZR_HF_RESULT, V_O_DELTA_MU) @pytest.mark.parametrize( From 317adae9897de523f72bd0976271971b125a2bce Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Tue, 29 Sep 2026 17:46:50 -0700 Subject: [PATCH 18/18] fix(SOF-7975): generic defaults, PRISTINE_NAME "Si" by Standata first match, SCF_KGRID None reads the pristine job's grid The charged notebook's PRISTINE_NAME is "Si" and is loaded as the sibling notebooks load a material, Materials.get_by_name_first_match (mp-149 for "Si"), falling back to load_material on the uploads folder and the platform when Standata has no match; only the lookup sits in the try, so only its not-found ValueError falls back. The exact-name filter is gone. SCF_KGRID defaults to None in both notebooks. With None, the pristine Total Energy reference is found on any grid, or created without one so the platform chooses it; its grid is read with get_kgrid_of_job, and the Band Gap reference is looked up or created on that grid, so the two share one grid. After the references finish, the notebook reads the Total Energy job's grid back, prints it, and uses it for the defective relaxations, the defect jobs and the relaxation lookups (get_scf_kgrid_for_supercell returns it). An explicit SCF_KGRID behaves as before, divided by the supercell scaling. The neutral notebook applies the grid of the pristine job it found to the defect SCF; with None, its missing-reference error asks to run Total Energy on the pristine first. The reuse line prints the grid of the job it reused, the creation line "k-grid: platform default" when none is set. The markdown states where the grid comes from. get_kgrid_of_job (job/api.py) returns a job unit's grid from its kgrid context, or, for a job created without one, from K_POINTS automatic in workflow.subworkflows[].units[name].input[0].rendered, which production renders at creation. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../workflows/defect_formation_energy.ipynb | 13 ++++--- .../defect_formation_energy_charged.ipynb | 35 +++++++++++-------- .../notebooks_utils/core/entity/job/api.py | 15 ++++++++ tests/py/unit/core/entity/test_job_api.py | 32 +++++++++++++++++ 4 files changed, 76 insertions(+), 19 deletions(-) diff --git a/other/materials_designer/workflows/defect_formation_energy.ipynb b/other/materials_designer/workflows/defect_formation_energy.ipynb index b0eef9631..9db16357f 100644 --- a/other/materials_designer/workflows/defect_formation_energy.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy.ipynb @@ -14,7 +14,7 @@ "- **[0] Defective supercell** — the structure containing the defect(s) (vacancy, substitution, interstitial, or a mix).\n", "- **[1] Pristine supercell** — the defect-free version of the *same* supercell.\n", "\n", - "Only the defective cell is computed here (`pw_scf`). The pristine reference energy is fetched from a previously-finished **Total Energy** job on the same k-grid, so the pristine must already be converged/relaxed and have such a job on the platform. Elemental chemical potentials are taken from Standata elemental reference materials (each needs a total energy too), as in [Formation Energy](formation_energy.ipynb).\n", + "Only the defective cell is computed here (`pw_scf`). The pristine reference energy is fetched from a previously-finished **Total Energy** job, so the pristine must already be converged/relaxed and have such a job on the platform. The k-grid comes from `SCF_KGRID`, or from the pristine reference job when it is None, and is the same for the pristine and defective cells. Elemental chemical potentials are taken from Standata elemental reference materials (each needs a total energy too), as in [Formation Energy](formation_energy.ipynb).\n", "\n", "Formula:\n", "\n", @@ -122,7 +122,7 @@ "metadata": {}, "outputs": [], "source": [ - "SCF_KGRID = [4, 4, 4] # same grid for the pristine and defective cells\n", + "SCF_KGRID = None # e.g. [4, 4, 4]; None takes the pristine reference job's grid\n", "CHEMICAL_POTENTIALS = {} # Δμ per element, eV/atom, e.g. {\"O\": -0.12}\n" ] }, @@ -361,7 +361,7 @@ "metadata": {}, "outputs": [], "source": [ - "from mat3ra.notebooks_utils.core.entity.job.api import find_job_for_material_with_property\n", + "from mat3ra.notebooks_utils.core.entity.job.api import find_job_for_material_with_property, get_kgrid_of_job\n", "from mat3ra.notebooks_utils.core.entity.material.api import get_or_create_material\n", "\n", "saved_defective = Material.create(get_or_create_material(client, defective_material, ACCOUNT_ID))\n", @@ -374,8 +374,11 @@ " raise RuntimeError(\n", " f\"No finished Total Energy job for '{saved_pristine.name}' on k-grid {SCF_KGRID}: \"\n", " \"run Total Energy on the pristine at this k-grid.\"\n", + " if SCF_KGRID\n", + " else f\"No finished Total Energy job for '{saved_pristine.name}': run Total Energy on the pristine first.\"\n", " )\n", - "print(f\"♻️ pristine reference: job {pristine_reference_job['_id']}, k-grid {SCF_KGRID}\")\n", + "kgrid = SCF_KGRID or get_kgrid_of_job(pristine_reference_job)\n", + "print(f\"♻️ pristine reference: job {pristine_reference_job['_id']}, k-grid {kgrid}\")\n", "\n", "# Order matters: [0] defective (computed), [1] pristine (reference).\n", "materials = [saved_defective, saved_pristine]\n" @@ -433,7 +436,7 @@ "print(f\"Loaded workflow: {defect_workflow.name}\")\n", "print(f\"Multi-material: {getattr(defect_workflow, 'isMultiMaterial', False)}\")\n", "\n", - "apply_scf_kgrid(defect_workflow, SCF_KGRID, material=defective_material)\n", + "apply_scf_kgrid(defect_workflow, kgrid, material=defective_material)\n", "set_assignment_value(defect_workflow, \"assign-reference-job-filter\", str({\"$in\": [pristine_reference_job[\"_id\"]]}))\n", "\n", "visualize_workflow(defect_workflow)" diff --git a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb index 01e7f2184..69bcfb402 100644 --- a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb +++ b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb @@ -90,7 +90,7 @@ "\n", "# 3. Material parameters\n", "FOLDER = \"../uploads\"\n", - "PRISTINE_NAME = \"Si, Silicon, FCC (Fd-3m) 3D (Bulk), mp-149\" # Standata name, or one in FOLDER or on the platform\n", + "PRISTINE_NAME = \"Si\" # Name of the material to load from Standata, uploads folder, or platform\n", "VISUALIZATION_REPETITIONS = [1, 1, 1]\n", "\n", "# 4. Workflow parameters\n", @@ -132,7 +132,7 @@ "\n", "CHARGES = [0] # e.g. [1, 0, -1, -2, -3]\n", "\n", - "SCF_KGRID = [4, 4, 4] # same grid for the pristine and defective cells, divided by the supercell scaling\n", + "SCF_KGRID = None # e.g. [4, 4, 4]; None takes the pristine reference job's grid\n", "CHEMICAL_POTENTIALS = {} # Δμ per element, eV/atom, e.g. {\"O\": -0.12}\n", "\n", "RELAX_PRISTINE_MATERIAL = False\n", @@ -288,7 +288,7 @@ "metadata": {}, "source": [ "### 3.3. Define how jobs are created\n", - "Every job is tagged with its charge; the pristine references run neutral, tagged `charge:0`. A reference is reused when a finished `charge:0` job on that pristine cell already reported the property (`total_energy`, `band_gaps`), on the same k-grid, `SCF_KGRID` divided by the supercell scaling. The Defect Formation Energy workflow reads each reference from the one job the notebook names to it, so no other calculation on the same cell can stand in; the relaxation carries no charge tag and is never taken for one." + "Every job is tagged with its charge; the pristine references run neutral, tagged `charge:0`. A reference is reused when a finished `charge:0` job on that pristine cell already reported the property (`total_energy`, `band_gaps`), on `SCF_KGRID` divided by the supercell scaling when `SCF_KGRID` is set. The k-grid comes from `SCF_KGRID`, or from the pristine reference job when it is None, and is the same for the pristine and defective cells. The Defect Formation Energy workflow reads each reference from the one job the notebook names to it, so no other calculation on the same cell can stand in; the relaxation carries no charge tag and is never taken for one." ] }, { @@ -337,12 +337,12 @@ "from mat3ra.notebooks_utils.core.entity.material.api import load_material\n", "from mat3ra.notebooks_utils.ipython.entity.material.visualize import visualize_materials\n", "\n", - "standata_matches = [data for data in Materials.get_by_name(PRISTINE_NAME) if data[\"name\"] == PRISTINE_NAME]\n", - "pristine_material = (\n", - " Material.create(standata_matches[0])\n", - " if standata_matches\n", - " else load_material(client, FOLDER, PRISTINE_NAME, ACCOUNT_ID)\n", - ")\n", + "try:\n", + " standata_material_config = Materials.get_by_name_first_match(PRISTINE_NAME)\n", + "except ValueError:\n", + " pristine_material = load_material(client, FOLDER, PRISTINE_NAME, ACCOUNT_ID)\n", + "else:\n", + " pristine_material = Material.create(standata_material_config)\n", "print(f\"Pristine material: {pristine_material.name} ({len(pristine_material.basis.elements.ids)} atoms)\")\n", "\n", "visualize_materials(pristine_material, repetitions=VISUALIZATION_REPETITIONS, title=\"Pristine material\")" @@ -561,7 +561,9 @@ "\n", "\n", "def get_scf_kgrid_for_supercell(scaling):\n", - " \"\"\"SCF_KGRID divided by the scaling, rounded, never below 1.\"\"\"\n", + " \"\"\"SCF_KGRID divided by the scaling, rounded, never below 1; with no SCF_KGRID, the size's pristine job's grid.\"\"\"\n", + " if SCF_KGRID is None:\n", + " return pristine_kgrids.get(scaling)\n", " return [max(1, round(dimension / scaling)) for dimension in SCF_KGRID]\n", "\n", "\n", @@ -642,10 +644,11 @@ "metadata": {}, "outputs": [], "source": [ - "from mat3ra.notebooks_utils.core.entity.job.api import find_job_for_material_with_property\n", + "from mat3ra.notebooks_utils.core.entity.job.api import find_job_for_material_with_property, get_kgrid_of_job\n", "\n", "pristine_reference_jobs = {scaling: {} for scaling in saved_supercell_pairs}\n", "pristine_job_ids = []\n", + "pristine_kgrids = {}\n", "for scaling, (pristine_supercell, _) in saved_supercell_pairs.items():\n", " kgrid = get_scf_kgrid_for_supercell(scaling)\n", " for property_name, workflow_config in pristine_workflow_configs.items():\n", @@ -653,15 +656,16 @@ " client, pristine_supercell.id, property_name, ACCOUNT_ID, tags=[\"charge:0\"], kgrid=kgrid\n", " )\n", " if job is not None:\n", - " print(f\"♻️ n={scaling}: reusing {property_name} from job {job['name']}, k-grid {kgrid}\")\n", + " print(f\"♻️ n={scaling}: reusing {property_name} from job {job['name']}, k-grid {get_kgrid_of_job(job)}\")\n", " else:\n", " workflow = Workflow.create(workflow_config)\n", " workflow.name = f\"{workflow.name} {pristine_supercell.name}\"\n", " apply_scf_kgrid(workflow, kgrid, material=pristine_supercell)\n", " job = create_job_for_materials([pristine_supercell], workflow, [\"charge:0\"])\n", " pristine_job_ids.append(job[\"_id\"])\n", - " print(f\"✅ n={scaling}: created a {property_name} job {job['_id']}, k-grid {kgrid}\")\n", + " print(f\"✅ n={scaling}: created a {property_name} job {job['_id']}, k-grid: {kgrid or 'platform default'}\")\n", " pristine_reference_jobs[scaling][property_name] = job\n", + " kgrid = kgrid or get_kgrid_of_job(job)\n", "\n", "if pristine_job_ids:\n", " submit_jobs(client.jobs, pristine_job_ids)\n", @@ -676,7 +680,10 @@ "outputs": [], "source": [ "if pristine_job_ids:\n", - " await wait_for_jobs_to_finish_async(client.jobs, pristine_job_ids, poll_interval=POLL_INTERVAL)" + " await wait_for_jobs_to_finish_async(client.jobs, pristine_job_ids, poll_interval=POLL_INTERVAL)\n", + "for scaling, reference_jobs in pristine_reference_jobs.items():\n", + " pristine_kgrids[scaling] = get_kgrid_of_job(client.jobs.get(reference_jobs[\"total_energy\"][\"_id\"]))\n", + " print(f\"n={scaling}: k-grid {get_scf_kgrid_for_supercell(scaling)}\")" ] }, { diff --git a/src/py/mat3ra/notebooks_utils/core/entity/job/api.py b/src/py/mat3ra/notebooks_utils/core/entity/job/api.py index 57c33dcd4..363c5b2f9 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/job/api.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/job/api.py @@ -1,3 +1,4 @@ +import re import urllib.request from typing import Any, Dict, Iterable, List, Optional, Union @@ -171,6 +172,20 @@ def get_kgrid_query(kgrid: Optional[List[int]], unit_name: str = "pw_scf") -> Di return {"workflow.subworkflows.units": {"$elemMatch": {"name": unit_name, "context": kgrid_context}}} +def get_kgrid_of_job(job: Dict[str, Any], unit_name: str = "pw_scf") -> Optional[List[int]]: + """ + K-grid dimensions the job's `unit_name` unit ran on: its `kgrid` context, or, for a job created without one, the + grid the platform rendered into the unit's input, `workflow.subworkflows[].units[name].input[0].rendered`. + """ + units = [unit for subworkflow in job["workflow"]["subworkflows"] for unit in subworkflow["units"]] + unit = next(unit for unit in units if unit["name"] == unit_name) + kgrid_context = next((item for item in unit["context"] if item["name"] == "kgrid"), None) + if kgrid_context: + return kgrid_context["data"]["dimensions"] + match = re.search(r"K_POINTS automatic\s+(\d+)\s+(\d+)\s+(\d+)", unit["input"][0]["rendered"]) + return [int(dimension) for dimension in match.groups()] if match else None + + def find_job_for_material_with_property( api_client: APIClient, material_id: str, diff --git a/tests/py/unit/core/entity/test_job_api.py b/tests/py/unit/core/entity/test_job_api.py index dc3a87934..d676a45b5 100644 --- a/tests/py/unit/core/entity/test_job_api.py +++ b/tests/py/unit/core/entity/test_job_api.py @@ -6,6 +6,7 @@ create_job, find_job_for_material, find_job_for_material_with_property, + get_kgrid_of_job, get_kgrid_query, ) @@ -219,3 +220,34 @@ def test_find_job_for_material_matches_the_kgrid_of_the_unit(): assert job == EXISTING_JOB assert RELAX_KGRID_QUERY.items() <= client.jobs.list.call_args.args[0].items() + + +# The pw_scf unit with a kgrid context, input not rendered yet, and as production job BLmZo5WZfFXKKTb2H stores a job +# created without a grid: no context, the platform's grid rendered into the input. +KGRID_CONTEXT: Dict[str, Any] = { + "name": "kgrid", + "isEdited": True, + "data": {"dimensions": [4, 4, 4], "shifts": [0, 0, 0], "gridMetricType": "KPPRA", "gridMetricValue": 768}, + "extraData": {"materialHash": "041d30e32f91e2eeb14c74298dffd08b"}, +} +RENDERED_INPUT = ( + "CELL_PARAMETERS angstrom\n 0.000000000 0.000000000 5.326038000\nK_POINTS automatic\n1 1 1 0 0 0 \n" +) +UNIT_WITH_KGRID_CONTEXT: Dict[str, Any] = { + "name": "pw_scf", + "context": [KGRID_CONTEXT], + "input": [{"template": {"name": "pw_scf.in"}, "rendered": "", "isManuallyChanged": False}], +} +UNIT_WITH_RENDERED_INPUT: Dict[str, Any] = { + "name": "pw_scf", + "context": [], + "input": [{"template": {"name": "pw_scf.in"}, "rendered": RENDERED_INPUT, "isManuallyChanged": False}], +} + + +@pytest.mark.parametrize( + ("unit", "expected_kgrid"), + [(UNIT_WITH_KGRID_CONTEXT, [4, 4, 4]), (UNIT_WITH_RENDERED_INPUT, [1, 1, 1])], +) +def test_get_kgrid_of_job(unit, expected_kgrid): + assert get_kgrid_of_job({"workflow": {"subworkflows": [{"units": [unit]}]}}) == expected_kgrid