diff --git a/other/materials_designer/workflows/defect_formation_energy.ipynb b/other/materials_designer/workflows/defect_formation_energy.ipynb index 10651a132..9db16357f 100644 --- a/other/materials_designer/workflows/defect_formation_energy.ipynb +++ b/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\"))" ] } ], diff --git a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb index ff38759ef..69bcfb402 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, 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.\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", @@ -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", @@ -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", @@ -131,9 +132,11 @@ "\n", "CHARGES = [0] # e.g. [1, 0, -1, -2, -3]\n", "\n", - "SCF_KGRID = None # e.g. [8, 8, 8]\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" + "RELAX_PRISTINE_MATERIAL = False\n", + "RELAX_DEFECTIVE_MATERIAL = False" ] }, { @@ -275,8 +278,8 @@ "else:\n", " cluster = clusters[0]\n", "\n", - "compute = Compute(cluster=cluster, queue=QUEUE_NAME, ppn=PPN)\n", - "print(f\"Using cluster: {compute.cluster.hostname}, queue: {QUEUE_NAME}, ppn: {PPN}\")" + "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}, time limit: {TIME_LIMIT}\")" ] }, { @@ -285,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`), 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 `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." ] }, { @@ -334,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\")" @@ -351,7 +354,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 k-grid `SCF_KGRID`, and the supercells below are built from the relaxed cell." ] }, { @@ -361,29 +364,37 @@ "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", + " 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", + " 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 (job {relax_job['_id']})\")\n", + " f\"{relaxed_pristine_material.lattice.a:.4f} Å relaxed\")\n", " pristine_material = relaxed_pristine_material" ] }, @@ -532,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", @@ -546,8 +561,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; 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", "def create_defect_workflow(scaling, charge, reference_jobs):\n", @@ -563,6 +580,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))" ] }, @@ -605,26 +644,28 @@ "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", " 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 {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 or 'from KPPRA'}\")\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", @@ -639,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)}\")" ] }, { @@ -661,12 +705,15 @@ "\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", + " 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_supercell, 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 = create_job_for_materials([defective_material, pristine_supercell], workflow, [f\"charge:{charge}\"])\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 +764,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. 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." ] }, { @@ -727,6 +774,14 @@ "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.plot._plotly import render_figure\n", + "\n", "results_records = []\n", "for record in job_records:\n", " job = client.jobs.get(record[\"job_id\"])\n", @@ -736,12 +791,28 @@ " \"final_status\": job.get(\"status\"),\n", " \"formation_energy\": properties[0].get(\"value\") if properties else None,\n", " })\n", + " if properties:\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(result, 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 not successful_df.empty:\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", + " 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" ] }, @@ -751,7 +822,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 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/src/py/mat3ra/notebooks_utils/core/entity/job/api.py b/src/py/mat3ra/notebooks_utils/core/entity/job/api.py index f2a6ff443..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 @@ -127,9 +128,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 +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): 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: dict, optional: The matching job, or None if none exists. @@ -147,22 +153,51 @@ 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}, ) 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 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, 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 +205,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): 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. @@ -177,7 +213,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/src/py/mat3ra/notebooks_utils/core/entity/material/api.py b/src/py/mat3ra/notebooks_utils/core/entity/material/api.py index 8cd3b3e03..4848f1a65 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): 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: 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/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..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,12 +1,13 @@ -"""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 +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 Dict, Sequence +from typing import Any, Dict, List, NamedTuple, Sequence import numpy as np import pandas as pd @@ -108,3 +109,52 @@ 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 + + +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_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]. + + 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 result.formation_energy - sum( + count * delta_mu[element] for element, count in result.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/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/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", ] diff --git a/tests/py/unit/core/entity/test_job_api.py b/tests/py/unit/core/entity/test_job_api.py index 8d90484b5..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,8 @@ create_job, find_job_for_material, find_job_for_material_with_property, + get_kgrid_of_job, + get_kgrid_query, ) OWNER_ID = "account-1" @@ -175,3 +177,77 @@ 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] + + +SCF_KGRID_QUERY: Dict[str, Any] = { + "workflow.subworkflows.units": { + "$elemMatch": {"name": "pw_scf", "context": {"$elemMatch": {"name": "kgrid", "data.dimensions": [4, 4, 4]}}} + } +} +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", "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_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=[4, 4, 4]) + + assert job == EXISTING_JOB + assert SCF_KGRID_QUERY.items() <= client.jobs.list.call_args.args[0].items() + + +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=[4, 4, 4], unit_name="pw_relax" + ) + + 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 diff --git a/tests/py/unit/core/entity/test_material_api.py b/tests/py/unit/core/entity/test_material_api.py index 2985ef887..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,6 +360,20 @@ 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" + assert get_kgrid_query([4, 4, 4], "pw_vc-relax").items() <= client.jobs.list.call_args.args[0].items() + + @pytest.mark.parametrize( ("folder_names", "account_names", "name", "expected_name"), [ 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..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,10 +9,14 @@ FORMATION_ENERGY_AT_VBM_COLUMN, STABLE_FROM_COLUMN, STABLE_TO_COLUMN, + DefectJobResult, evaluate_finite_size_fit, fit_finite_size, get_charge_state_table, + get_chemical_potential_combination, + get_defect_job_result, 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 +64,63 @@ 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). +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( + "result, delta_mu, expected", + [ + (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(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(ZR_HF_RESULT, V_O_DELTA_MU) + + +@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