SOF-8067: charged defect notebook — relax the defective cell per charge, reuse references only on the same k-grid

other/materials_designer/workflows/defect_formation_energy.ipynb

@@ -14,13 +14,13 @@
     "- **[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, 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",
     "$$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",
+    "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",
@@ -122,22 +122,8 @@
    "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",
-    "\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",
-    "PRISTINE_TOTAL_ENERGY_SOURCE = \"my_account\"\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"
    ]
   },
   {
@@ -375,17 +361,24 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "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",
-    "from mat3ra.notebooks_utils.core.entity.property.api import find_total_energy_for_material\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",
+    "pristine_reference_job = find_job_for_material_with_property(\n",
+    "    client, saved_pristine.id, \"total_energy\", ACCOUNT_ID, kgrid=SCF_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 {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",
+    "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"
@@ -432,9 +425,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 patch_workflow_qe_input\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",
@@ -444,18 +436,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",
-    "\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",
+    "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)"
    ]
@@ -579,7 +561,8 @@
    "metadata": {},
    "source": [
     "## 7. Retrieve results\n",
-    "### 7.1. Retrieve and visualize defect formation energy"
+    "### 7.1. Retrieve and visualize defect formation energy\n",
+    "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$."
    ]
   },
   {
@@ -589,10 +572,28 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "from mat3ra.notebooks_utils.core.entity.property.defect_analysis import (\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\")"
+    "visualize_properties(defect_energy_data, title=\"Defect Formation Energy\")\n",
+    "\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_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",
+    "    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\"))"
    ]
   }
  ],