diff --git a/other/materials_designer/workflows/Introduction.ipynb b/other/materials_designer/workflows/Introduction.ipynb index 5a9cd46d..332a609f 100644 --- a/other/materials_designer/workflows/Introduction.ipynb +++ b/other/materials_designer/workflows/Introduction.ipynb @@ -70,6 +70,7 @@ "\n", "### 6.5. Defect Energy\n", "#### [6.5.1. Defect formation energy.](defect_formation_energy.ipynb)\n", + "#### [6.5.2. Defect formation energy of charged defects.](defect_formation_energy_charged.ipynb)\n", "\n", "### 6.6. Formation Energy\n", "#### [6.6.1. Compound formation energy.](formation_energy.ipynb)\n", diff --git a/other/materials_designer/workflows/defect_formation_energy_charged.ipynb b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb new file mode 100644 index 00000000..4150a2c7 --- /dev/null +++ b/other/materials_designer/workflows/defect_formation_energy_charged.ipynb @@ -0,0 +1,910 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Defect Formation Energy of Charged Defects\n", + "\n", + "Calculate the formation energy of a defect as a function of its charge state and of the supercell size, with a multi-material DFT workflow on the Mat3ra platform. For the neutral defect in a single supercell see [Defect Formation Energy](defect_formation_energy.ipynb).\n", + "\n", + "One job is created per (supercell size, charge state) pair. Each job takes **two materials, in order**:\n", + "\n", + "- **[0] Defective supercell** — the pristine supercell with `DEFECT_CONFIGS` applied.\n", + "- **[1] Pristine supercell** — the defect-free supercell of the same size.\n", + "\n", + "The workflow reports the raw formation energy of the charged cell,\n", + "\n", + "$$E_f^{\\text{raw}}[X^q] = E_{\\text{tot}}[X^q] - E_{\\text{tot}}[\\text{bulk}] - \\sum_i \\Delta N_i\\, \\mu_i \\quad [\\text{eV}]$$\n", + "\n", + "and this notebook adds the electron-reservoir term, giving the formation energy as a function of the electron chemical potential $\\mu_e$ measured from the valence band maximum:\n", + "\n", + "$$E_f[X^q](\\mu_e) = E_f^{\\text{raw}}[X^q] + q\\,(E_{\\text{VBM}} + \\mu_e), \\qquad 0 \\le \\mu_e \\le E_{\\text{gap}}$$\n", + "\n", + "$E_{\\text{VBM}}$ and $E_{\\text{gap}}$ come from 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", + "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", + "\n", + "

Usage

\n", + "\n", + "1. Set the pristine material, the defects, the supercell sizes and the charge states in cells 1.2 and 1.3 below.\n", + "1. Click \"Run\" > \"Run All\".\n", + "1. Total Energy and Band Gap jobs are created for any pristine supercell that lacks them, then one Defect Formation Energy job per (size, charge).\n", + "1. Scroll down for the results table, the formation energy against the electron chemical potential with the stable charge states, and the finite-size extrapolation.\n", + "\n", + "## Summary\n", + "\n", + "1. Set up the environment and parameters.\n", + "1. Authenticate and initialize API client.\n", + "1. Build a pristine/defective supercell pair per size, resolve elemental references, and save them.\n", + "1. Configure the Defect Formation Energy workflow per charge state and the pristine reference workflows.\n", + "1. Configure compute.\n", + "1. Create, submit, and monitor the pristine reference jobs and one defect job per (size, charge).\n", + "1. Retrieve results: raw formation energies, the dependence on the electron chemical potential, and the finite-size extrapolation." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "## 1. Set up the environment and parameters\n", + "### 1.1. Install packages (JupyterLite)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.packages import install_packages\n", + "\n", + "await install_packages(\"made|api_examples\")" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "### 1.2. Set parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "from datetime import datetime\n", + "from mat3ra.ide.compute import QueueName\n", + "\n", + "# 2. Auth and organization parameters\n", + "ORGANIZATION_NAME = None\n", + "\n", + "# 3. Material parameters\n", + "FOLDER = \"../uploads\"\n", + "PRISTINE_NAME = \"Si\" # defect-free cell, loaded from FOLDER or from Standata by name\n", + "VISUALIZATION_REPETITIONS = [1, 1, 1]\n", + "\n", + "# 4. Workflow parameters\n", + "WORKFLOW_SEARCH_TERM = \"defect_formation_energy.json\"\n", + "APPLICATION_NAME = \"espresso\"\n", + "MY_WORKFLOW_NAME = \"Defect Formation Energy\"\n", + "\n", + "# 5. Compute parameters\n", + "CLUSTER_NAME = None # specify full or partial name i.e. \"cluster-001\" to select\n", + "QUEUE_NAME = QueueName.D\n", + "PPN = 1\n", + "\n", + "# 6. Job parameters\n", + "timestamp = datetime.now().strftime(\"%Y-%m-%d %H:%M\")\n", + "POLL_INTERVAL = 30 # seconds" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "### 1.3. Set specific charge state and supercell size parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "# Defects placed in every supercell, in crystal coordinates of the pristine cell.\n", + "# See create_point_defect.ipynb for the full defect configuration reference.\n", + "DEFECT_CONFIGS = [\n", + " {\"type\": \"vacancy\", \"coordinate\": [0.0, 0.0, 0.0], \"placement_method\": \"closest_site\"},\n", + "]\n", + "\n", + "# Supercell sizes: n -> n x n x n repetitions of the pristine cell. Two or more\n", + "# sizes enable the finite-size extrapolation in section 7.3.\n", + "SUPERCELL_SCALINGS = [1] # e.g. [2, 3]\n", + "\n", + "# Net charge q of the defective supercell in units of e, one job per value.\n", + "CHARGES = [0] # e.g. [1, 0, -1, -2, -3]\n", + "\n", + "# K-grid for the pristine cell. A supercell of size n uses SCF_KGRID / n rounded, which\n", + "# holds the k-point density roughly fixed -- until it bottoms out at 1 and larger cells\n", + "# are all sampled at Gamma. If not set, KPPRA is used by default.\n", + "SCF_KGRID = None # e.g. [8, 8, 8]\n", + "\n", + "# Whose band_gaps property to reuse for the pristine supercells: \"public\" (any owner,\n", + "# highest precision wins), \"curators\" (only curators'), or \"my_account\" (your own, or\n", + "# the organization's when ORGANIZATION_NAME is set). The pristine total energy is not\n", + "# scoped by this -- the workflow finds it by job, so section 6.1 does too.\n", + "PRISTINE_PROPERTY_SOURCE = \"my_account\"" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "## 2. Authenticate and initialize API client\n", + "### 2.1. Authenticate\n", + "Authenticate in the browser and have credentials stored in environment variable `OIDC_ACCESS_TOKEN`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.auth import authenticate\n", + "\n", + "await authenticate()" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "### 2.2. Initialize API client" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.api_client import APIClient\n", + "\n", + "client = APIClient.authenticate()\n", + "client" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "### 2.3. Select account to work under" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "client.list_accounts()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "selected_account = client.my_account\n", + "\n", + "if ORGANIZATION_NAME:\n", + " selected_account = client.get_account(name=ORGANIZATION_NAME)\n", + "\n", + "ACCOUNT_ID = selected_account.id\n", + "print(f\"✅ Selected account ID: {ACCOUNT_ID}, name: {selected_account.name}\")" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "### 2.4. Select project" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "projects = client.projects.list({\"isDefault\": True, \"owner._id\": ACCOUNT_ID})\n", + "project_id = projects[0][\"_id\"]\n", + "print(f\"✅ Using project: {projects[0]['name']} ({project_id})\")" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "## 3. Build the supercell pairs\n", + "### 3.1. Load the pristine cell" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17", + "metadata": {}, + "outputs": [], + "source": [ + "import re\n", + "from mat3ra.made.material import Material\n", + "from mat3ra.standata.materials import Materials\n", + "from mat3ra.notebooks_utils.material import load_material_from_folder\n", + "from mat3ra.notebooks_utils.ipython.entity.material.visualize import visualize_materials\n", + "\n", + "\n", + "def load_pristine(name):\n", + " \"\"\"Look in the uploads folder, then among your platform materials, then in Standata.\"\"\"\n", + " loaded = load_material_from_folder(FOLDER, name)\n", + " if loaded is not None:\n", + " return loaded\n", + " matches = client.materials.list(\n", + " {\"name\": {\"$regex\": re.escape(name), \"$options\": \"i\"}, \"owner._id\": ACCOUNT_ID}\n", + " )\n", + " return Material.create(matches[0] if matches else Materials.get_by_name_first_match(name))\n", + "\n", + "\n", + "pristine = load_pristine(PRISTINE_NAME)\n", + "print(f\"Pristine cell: {pristine.name} ({len(pristine.basis.elements.ids)} atoms)\")\n", + "\n", + "visualize_materials(pristine, repetitions=VISUALIZATION_REPETITIONS, title=\"Pristine cell\")" + ] + }, + { + "cell_type": "markdown", + "id": "18", + "metadata": {}, + "source": [ + "### 3.2. Build a pristine/defective pair per supercell size\n", + "The pristine cell is scaled first and the defects are placed into each supercell, so that every size holds the same isolated defect at the same site. Defect coordinates are given in the pristine cell's crystal basis and are divided by the scaling to reach that site in the supercell.\n", + "\n", + "Atom labels are stripped: they are useful for analysis notebooks but are not compatible with QE input generation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.made.tools.helpers import PointDefectDict, create_multiple_defects, create_supercell\n", + "\n", + "\n", + "def build_pair(scaling):\n", + " supercell = create_supercell(pristine, scaling_factor=[scaling] * 3)\n", + " defect_dicts = [\n", + " PointDefectDict(**{**config, \"coordinate\": [component / scaling for component in config[\"coordinate\"]]})\n", + " for config in DEFECT_CONFIGS\n", + " ]\n", + " defective = create_multiple_defects(material=supercell, defect_dicts=defect_dicts)\n", + " for material, kind in ((supercell, \"pristine\"), (defective, \"defective\")):\n", + " material.basis.set_labels_from_list([])\n", + " material.name = f\"{pristine.name} {scaling}x{scaling}x{scaling} {kind}\"\n", + " return supercell, defective\n", + "\n", + "\n", + "pairs = {scaling: build_pair(scaling) for scaling in SUPERCELL_SCALINGS}\n", + "\n", + "visualize_materials(\n", + " [{\"material\": material, \"title\": material.name} for pair in pairs.values() for material in pair],\n", + " repetitions=VISUALIZATION_REPETITIONS,\n", + " rotation=\"-90x\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "20", + "metadata": {}, + "source": [ + "### 3.3. Resolve Standata elemental reference materials\n", + "Elemental chemical potentials come from Standata materials tagged `elemental` with `metadata.element`. Each must also have a refined `total_energy` -- this is enforced by the workflow at runtime; run [Total Energy](total_energy.ipynb) for any elemental reference that is missing one.\n", + "\n", + "The chemical-potential term covers every species whose count changes, so elements are taken from the union of the pristine and defective structures." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21", + "metadata": {}, + "outputs": [], + "source": [ + "elements = sorted(\n", + " {element for pair in pairs.values() for material in pair for element in material.basis.elements.values}\n", + ")\n", + "elemental_materials_data = client.materials.list(\n", + " {\"tags\": \"elemental\", \"metadata.element\": {\"$in\": elements}},\n", + ")\n", + "available = {data.get(\"metadata\", {}).get(\"element\") for data in elemental_materials_data}\n", + "missing = sorted(set(elements) - available)\n", + "if missing:\n", + " raise RuntimeError(\n", + " f\"Missing elemental reference material(s) for {missing}. \"\n", + " \"Add them from Standata, or tag an existing material 'elemental' with metadata.element set.\"\n", + " )\n", + "print(f\"Resolved elemental reference materials for: {', '.join(elements)}\")" + ] + }, + { + "cell_type": "markdown", + "id": "22", + "metadata": {}, + "source": [ + "### 3.4. Save the supercells to the platform" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.core.entity.material.api import get_or_create_material\n", + "\n", + "saved_pairs = {\n", + " scaling: tuple(Material.create(get_or_create_material(client, material, ACCOUNT_ID)) for material in pair)\n", + " for scaling, pair in pairs.items()\n", + "}\n", + "for scaling, (saved_pristine, saved_defective) in saved_pairs.items():\n", + " print(f\"✅ n={scaling}: pristine {saved_pristine.id} ({len(saved_pristine.basis.elements.ids)} atoms), \"\n", + " f\"defective {saved_defective.id} ({len(saved_defective.basis.elements.ids)} atoms)\")" + ] + }, + { + "cell_type": "markdown", + "id": "24", + "metadata": {}, + "source": [ + "## 4. Configure the workflows\n", + "### 4.1. Select application" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "25", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.ade.application import Application\n", + "from mat3ra.standata.applications import ApplicationStandata\n", + "\n", + "app_config = ApplicationStandata.get_by_name_first_match(APPLICATION_NAME)\n", + "app = Application(**app_config)\n", + "print(f\"Using application: {app.name}\")" + ] + }, + { + "cell_type": "markdown", + "id": "26", + "metadata": {}, + "source": [ + "### 4.2. Load the Defect Formation Energy workflow and apply size and charge\n", + "One workflow per (size, charge): the charge is written into the `&SYSTEM` namelist of the defective-cell SCF, and the k-grid is scaled down with the supercell, which holds the k-point density roughly fixed while the divided grid stays above 1. The first combination is previewed below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "27", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.standata.workflows import WorkflowStandata\n", + "from mat3ra.wode.workflows import Workflow\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\n", + "\n", + "defect_workflow_config = WorkflowStandata.filter_by_application(app.name).get_by_name_first_match(\n", + " WORKFLOW_SEARCH_TERM\n", + ")\n", + "\n", + "\n", + "def scf_kgrid_for(scaling):\n", + " \"\"\"The pristine grid divided by the scaling, rounded, never below 1.\n", + "\n", + " This holds the k-point density roughly fixed across sizes, but only while the result stays\n", + " above 1: a coarse SCF_KGRID at a large scaling bottoms out at a Gamma-only grid, and sizes\n", + " that sample the Brillouin zone differently should not be compared in section 7.3.\n", + " \"\"\"\n", + " return None if SCF_KGRID is None else [max(1, round(dimension / scaling)) for dimension in SCF_KGRID]\n", + "\n", + "\n", + "def defect_workflow_for(scaling, charge):\n", + " workflow = Workflow.create(defect_workflow_config)\n", + " workflow.name = f\"{MY_WORKFLOW_NAME} n={scaling} q={charge:+d}\"\n", + " if charge:\n", + " patch_workflow_qe_input(workflow, {\"system\": {\"tot_charge\": charge}}, unit_names=[\"pw_scf\"])\n", + " return apply_scf_kgrid(workflow, scf_kgrid_for(scaling), material=saved_pairs[scaling][1])\n", + "\n", + "\n", + "visualize_workflow(defect_workflow_for(SUPERCELL_SCALINGS[0], CHARGES[0]))" + ] + }, + { + "cell_type": "markdown", + "id": "28", + "metadata": {}, + "source": [ + "### 4.3. Load the pristine reference workflows\n", + "Each pristine supercell supplies two references: its total energy, which the workflow subtracts, and its band gap, which supplies $E_{\\text{VBM}}$ and the range of $\\mu_e$. Both are reused when the material already carries them and computed otherwise." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "29", + "metadata": {}, + "outputs": [], + "source": [ + "pristine_workflow_configs = {\n", + " property_name: WorkflowStandata.filter_by_application(app.name).get_by_name_first_match(search_term)\n", + " for property_name, search_term in {\"total_energy\": \"total_energy.json\", \"band_gaps\": \"band_gap.json\"}.items()\n", + "}\n", + "print(f\"Loaded pristine reference workflows: {', '.join(pristine_workflow_configs)}\")" + ] + }, + { + "cell_type": "markdown", + "id": "30", + "metadata": {}, + "source": [ + "## 5. Create the compute configuration\n", + "### 5.1. Select cluster" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31", + "metadata": {}, + "outputs": [], + "source": [ + "clusters = client.clusters.list()\n", + "print(f\"Available clusters: {[c['hostname'] for c in clusters]}\")" + ] + }, + { + "cell_type": "markdown", + "id": "32", + "metadata": {}, + "source": [ + "### 5.2. Create compute configuration" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "33", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.ide.compute import Compute\n", + "\n", + "if CLUSTER_NAME:\n", + " cluster = next((c for c in clusters if CLUSTER_NAME in c[\"hostname\"]), None)\n", + "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}\")" + ] + }, + { + "cell_type": "markdown", + "id": "34", + "metadata": {}, + "source": [ + "## 6. Create and run the jobs\n", + "### 6.1. Create and submit the pristine reference jobs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "35", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.api.job import submit_jobs, wait_for_jobs_to_finish_async\n", + "from mat3ra.notebooks_utils.core.entity.property.api import find_property_for_material\n", + "from mat3ra.notebooks_utils.job import create_job\n", + "\n", + "\n", + "def create_job_for(materials, workflow):\n", + " job = create_job(\n", + " api_client=client,\n", + " materials=materials,\n", + " workflow=workflow,\n", + " project_id=project_id,\n", + " owner_id=ACCOUNT_ID,\n", + " prefix=f\"{workflow.name} {timestamp}\",\n", + " compute=compute.to_dict(),\n", + " )\n", + " return job[0] if isinstance(job, list) else job\n", + "\n", + "\n", + "def pristine_reference_exists(material_id, property_name):\n", + " \"\"\"Return the existing reference, or None. Each is checked the way its own consumer resolves it.\n", + "\n", + " The workflow resolves E_PRISTINE through the most recent *finished job* named \"Total Energy\",\n", + " then that job's Quantum ESPRESSO total_energy (shell/fetch_bulk_total_energy) -- a property on\n", + " its own, however it was curated, does not satisfy it. The band gap is read here in section 7.2,\n", + " as a property, so it is looked up as one.\n", + " \"\"\"\n", + " if property_name != \"total_energy\":\n", + " return find_property_for_material(\n", + " client, material_id, property_name, PRISTINE_PROPERTY_SOURCE, owner_id=ACCOUNT_ID\n", + " )\n", + " jobs = client.jobs.list(\n", + " {\"_material._id\": material_id, \"status\": \"finished\", \"workflow.name\": {\"$regex\": \"Total Energy\"}},\n", + " {\"sort\": {\"updatedAt\": -1}, \"limit\": 1},\n", + " )\n", + " if not jobs:\n", + " return None\n", + " # Only the most recent such job counts, and only a Quantum ESPRESSO result satisfies the\n", + " # workflow -- a Total Energy job run under another application carries the same name.\n", + " energies = client.properties.list(\n", + " {\"source.info.jobId\": jobs[0][\"_id\"], \"slug\": \"total_energy\", \"group\": {\"$regex\": \"qe:\"}}\n", + " )\n", + " return jobs[0] if energies else None\n", + "\n", + "\n", + "def reference_matches_kgrid(job, scaling):\n", + " \"\"\"True when the reference job sampled the Brillouin zone the way this size will.\n", + "\n", + " A reference computed at another k-grid is the quiet way these numbers go wrong: E_PRISTINE and\n", + " E_VBM would come from a different sampling than the defective cell they are subtracted from.\n", + " Only an explicitly set SCF_KGRID can be compared -- the KPPRA default adapts per material.\n", + " \"\"\"\n", + " expected = scf_kgrid_for(scaling)\n", + " if expected is None:\n", + " return True\n", + " for subworkflow in job.get(\"workflow\", {}).get(\"subworkflows\", []):\n", + " for unit in subworkflow.get(\"units\", []):\n", + " for context in unit.get(\"context\", {}).get(\"kgrid\", {}), unit.get(\"kgrid\", {}):\n", + " dimensions = (context or {}).get(\"dimensions\")\n", + " if dimensions:\n", + " return list(dimensions) == list(expected)\n", + " return False\n", + "\n", + "\n", + "pristine_job_ids = []\n", + "for scaling, (saved_pristine, _) in saved_pairs.items():\n", + " for property_name, workflow_config in pristine_workflow_configs.items():\n", + " existing = pristine_reference_exists(saved_pristine.id, property_name)\n", + " if existing and not reference_matches_kgrid(existing, scaling):\n", + " print(f\"↻ n={scaling}: existing {property_name} used a different k-grid; recomputing\")\n", + " existing = None\n", + " if existing:\n", + " # Name it: a reference computed at a different k-grid is the one silent way the\n", + " # numbers below go wrong.\n", + " print(f\"♻️ n={scaling}: reusing {property_name} from {existing.get('name', existing['_id'])}\")\n", + " continue\n", + " workflow = Workflow.create(workflow_config)\n", + " # Keep the standard name as the prefix: the defect workflow finds this job by a\n", + " # regex on \"Total Energy\".\n", + " workflow.name = f\"{workflow.name} {saved_pristine.name}\"\n", + " apply_scf_kgrid(workflow, scf_kgrid_for(scaling), material=saved_pristine)\n", + " print(f\" k-grid {scf_kgrid_for(scaling) or 'from KPPRA'}\")\n", + " job = create_job_for([saved_pristine], workflow)\n", + " pristine_job_ids.append(job[\"_id\"])\n", + " print(f\"✅ n={scaling}: created a {property_name} job {job['_id']}\")\n", + "\n", + "if pristine_job_ids:\n", + " submit_jobs(client.jobs, pristine_job_ids)\n", + " print(f\"✅ Submitted {len(pristine_job_ids)} pristine reference job(s).\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "36", + "metadata": {}, + "outputs": [], + "source": [ + "if pristine_job_ids:\n", + " await wait_for_jobs_to_finish_async(client.jobs, pristine_job_ids, poll_interval=POLL_INTERVAL)" + ] + }, + { + "cell_type": "markdown", + "id": "37", + "metadata": {}, + "source": [ + "### 6.2. Create the Defect Formation Energy jobs, one per size and charge" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "38", + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "\n", + "job_records = []\n", + "for scaling in SUPERCELL_SCALINGS:\n", + " saved_pristine, saved_defective = saved_pairs[scaling]\n", + " for charge in CHARGES:\n", + " # Order matters: [0] defective (computed), [1] pristine (reference).\n", + " job = create_job_for([saved_defective, saved_pristine], defect_workflow_for(scaling, charge))\n", + " print(f\"n={scaling} q={charge:+d}: k-grid {scf_kgrid_for(scaling) or 'from KPPRA'}, job {job['_id']}\")\n", + " job_records.append({\n", + " \"scaling\": scaling,\n", + " \"charge\": charge,\n", + " \"atoms\": len(saved_defective.basis.elements.ids),\n", + " \"length\": saved_pristine.lattice.cell_volume ** (1 / 3),\n", + " \"job_id\": job[\"_id\"],\n", + " })\n", + "\n", + "job_ids = [record[\"job_id\"] for record in job_records]\n", + "pd.DataFrame(job_records)" + ] + }, + { + "cell_type": "markdown", + "id": "39", + "metadata": {}, + "source": [ + "### 6.3. Submit the jobs and monitor the statuses" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "40", + "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": "41", + "metadata": {}, + "outputs": [], + "source": [ + "if job_ids:\n", + " await wait_for_jobs_to_finish_async(client.jobs, job_ids, poll_interval=POLL_INTERVAL)" + ] + }, + { + "cell_type": "markdown", + "id": "42", + "metadata": {}, + "source": [ + "## 7. Retrieve results\n", + "### 7.1. Raw formation energies\n", + "`length` is $L = V^{1/3}$ of the supercell, the length scale the image-charge terms in section 7.3 are expressed in." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43", + "metadata": {}, + "outputs": [], + "source": [ + "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", + " results_records.append({\n", + " **record,\n", + " \"final_status\": job.get(\"status\"),\n", + " \"raw_formation_energy\": properties[0].get(\"value\") if properties else None,\n", + " })\n", + "\n", + "results_df = pd.DataFrame(results_records)\n", + "successful_df = results_df[(results_df[\"final_status\"] == \"finished\") & results_df[\"raw_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", + "results_df" + ] + }, + { + "cell_type": "markdown", + "id": "44", + "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 is used, with $E_{\\text{VBM}}$ and $E_{\\text{gap}}$ from its own Band Gap job." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "45", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import plotly.graph_objects as go\n", + "from mat3ra.notebooks_utils.ipython.plot._plotly import render_figure\n", + "\n", + "# Every charge state must be present at the size plotted: the lower envelope is a comparison\n", + "# between them, so a missing one silently reports the wrong stable charge. This section is\n", + "# skipped rather than raised on, so that 7.3 -- which needs neither -- still runs.\n", + "complete_scalings = [\n", + " scaling\n", + " for scaling, group in successful_df.groupby(\"scaling\")\n", + " if set(group[\"charge\"]) == set(CHARGES)\n", + "]\n", + "largest_scaling = max(complete_scalings) if complete_scalings else None\n", + "band_gaps_property = (\n", + " find_property_for_material(\n", + " client, saved_pairs[largest_scaling][0].id, \"band_gaps\", PRISTINE_PROPERTY_SOURCE, owner_id=ACCOUNT_ID\n", + " )\n", + " if largest_scaling is not None\n", + " else None\n", + ")\n", + "fundamental_gap = (\n", + " min(band_gaps_property[\"data\"][\"values\"], key=lambda entry: entry[\"value\"])\n", + " if band_gaps_property\n", + " else None\n", + ")\n", + "valence_band_maximum = fundamental_gap[\"eigenvalueValence\"] if fundamental_gap else None\n", + "\n", + "if largest_scaling is None:\n", + " print(\"⚠️ No supercell size has a result for every charge state in CHARGES, so the lower \"\n", + " \"envelope would compare an incomplete set. Skipping this section -- 7.3 still runs.\")\n", + "elif valence_band_maximum is None:\n", + " print(f\"⚠️ No band gap with a valence-band eigenvalue for the n={largest_scaling} pristine cell: \"\n", + " \"its Band Gap job did not finish, or PRISTINE_PROPERTY_SOURCE excludes it. Skipping this \"\n", + " \"section -- 7.3 still runs, on the raw energies.\")\n", + "else:\n", + " band_gap = fundamental_gap[\"value\"]\n", + " print(f\"Pristine n={largest_scaling}: VBM {valence_band_maximum:.4f} eV, \"\n", + " f\"{fundamental_gap['type']} gap {band_gap:.4f} eV\")\n", + "\n", + " electron_chemical_potential = np.linspace(0, band_gap, 201)\n", + " largest_df = successful_df[successful_df[\"scaling\"] == largest_scaling]\n", + " formation_energies = {\n", + " row.charge: row.raw_formation_energy + row.charge * (valence_band_maximum + electron_chemical_potential)\n", + " for row in largest_df.itertuples()\n", + " }\n", + "\n", + " figure = go.Figure()\n", + " for charge, energies in sorted(formation_energies.items()):\n", + " figure.add_scatter(x=electron_chemical_potential, y=energies, mode=\"lines\", name=f\"q = {charge:+d}\")\n", + " figure.add_scatter(\n", + " x=electron_chemical_potential,\n", + " y=np.min(list(formation_energies.values()), axis=0),\n", + " mode=\"lines\",\n", + " name=\"stable state\",\n", + " line={\"width\": 8, \"color\": \"black\"},\n", + " opacity=0.2,\n", + " )\n", + " figure.update_layout(\n", + " title=f\"Defect formation energy vs electron chemical potential (n={largest_scaling})\",\n", + " xaxis_title=\"Electron chemical potential above VBM (eV)\",\n", + " yaxis_title=\"Formation energy (eV)\",\n", + " )\n", + " render_figure(figure)\n", + "\n", + " stable_charges = np.array(\n", + " [\n", + " min(formation_energies, key=lambda charge: formation_energies[charge][index])\n", + " for index in range(len(electron_chemical_potential))\n", + " ]\n", + " )\n", + " for charge, energies in sorted(formation_energies.items()):\n", + " window = electron_chemical_potential[stable_charges == charge]\n", + " stability = f\"stable for mu_e in [{window.min():.3f}, {window.max():.3f}] eV\" if window.size else \"never stable\"\n", + " print(f\"q = {charge:+d}: E_f = {energies[0]:.4f} eV at mu_e = 0 (VBM), {stability}\")" + ] + }, + { + "cell_type": "markdown", + "id": "46", + "metadata": {}, + "source": [ + "### 7.3. Finite-size extrapolation\n", + "The raw formation energy of a charged cell carries the spurious interaction of the defect with its periodic images, which falls off as $a/L + b/L^3$. Fitting the sizes in `SUPERCELL_SCALINGS` extrapolates to the isolated defect $E_\\infty$, in place of an analytical correction; $b$ is only fitted with four or more sizes, since three would make the system square and the \"fit\" an interpolation. `E_f_at_vbm` is $E_\\infty + q E_{\\text{VBM}}$ -- the extrapolated formation energy at $\\mu_e = 0$, directly comparable to what 7.2 prints." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "47", + "metadata": {}, + "outputs": [], + "source": [ + "def extrapolate(lengths, energies):\n", + " terms = [np.ones_like(lengths), 1 / lengths] + ([1 / lengths**3] if len(lengths) > 3 else [])\n", + " coefficients, *_ = np.linalg.lstsq(np.stack(terms, axis=1), energies, rcond=None)\n", + " return dict(zip((\"E_inf\", \"a\", \"b\"), coefficients))\n", + "\n", + "\n", + "fit_records = []\n", + "figure = go.Figure()\n", + "for charge, group in successful_df.groupby(\"charge\"):\n", + " lengths, energies = group[\"length\"].to_numpy(), group[\"raw_formation_energy\"].to_numpy()\n", + " if len(lengths) < 2:\n", + " print(f\"q = {charge:+d}: E_f = {energies[0]:.4f} eV at n={group['scaling'].iloc[0]}; \"\n", + " \"add sizes to SUPERCELL_SCALINGS to extrapolate to the isolated defect.\")\n", + " continue\n", + " fit = extrapolate(lengths, energies)\n", + " fit_records.append({\"charge\": charge, **fit})\n", + " inverse_lengths = np.linspace(0, 1 / lengths.min(), 50)\n", + " figure.add_scatter(x=1 / lengths, y=energies, mode=\"markers\", name=f\"q = {charge:+d}\")\n", + " figure.add_scatter(\n", + " x=inverse_lengths,\n", + " y=fit[\"E_inf\"] + fit[\"a\"] * inverse_lengths + fit.get(\"b\", 0.0) * inverse_lengths**3,\n", + " mode=\"lines\",\n", + " showlegend=False,\n", + " )\n", + "\n", + "if fit_records:\n", + " figure.update_layout(\n", + " title=\"Finite-size extrapolation of the defect formation energy\",\n", + " xaxis_title=\"1 / L (1/Angstrom), L = V^(1/3)\",\n", + " yaxis_title=\"Raw formation energy (eV)\",\n", + " )\n", + " render_figure(figure)\n", + "\n", + "fit_df = pd.DataFrame(fit_records)\n", + "if fit_records and valence_band_maximum is not None:\n", + " # Same quantity section 7.2 plots at mu_e = 0, now with the image-charge error extrapolated out.\n", + " fit_df[\"E_f_at_vbm\"] = fit_df[\"E_inf\"] + fit_df[\"charge\"] * valence_band_maximum\n", + "fit_df" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "pygments_lexer": "ipython3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/py/mat3ra/notebooks_utils/core/entity/property/api.py b/src/py/mat3ra/notebooks_utils/core/entity/property/api.py index 050c16f4..7a50ba04 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/property/api.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/property/api.py @@ -71,9 +71,15 @@ def update_property_holder_value(client: APIClient, property_holder_id: str, val return client.properties.update(property_holder_id, {"$set": {"data.value": value}}) -def find_total_energy_for_material(client: APIClient, material_id: str, source: str = "my_account") -> Optional[dict]: +def find_property_for_material( + client: APIClient, + material_id: str, + property_name: str, + source: str = "my_account", + owner_id: Optional[str] = None, +) -> Optional[dict]: """ - Find the best-precision total_energy property for a material. Mirrors the + Find the best-precision property of the given name for a material. Mirrors the platform's "Resolve Total Energies for Elemental Materials" subworkflow, which queries properties directly by material and selects by precision -- no job lookup involved. @@ -83,22 +89,26 @@ def find_total_energy_for_material(client: APIClient, material_id: str, source: Args: client (APIClient): API client instance. - material_id (str): Material _id to look up the total_energy property for. - source (str): Source of the total energy property: `my_account` (default), `curators` or + material_id (str): Material _id to look up the property for. + property_name (str): Property name, e.g. `total_energy` or `band_gaps`. + source (str): Source of the property: `my_account` (default), `curators` or `public`. + owner_id (str, optional): Account the `my_account` scope resolves to. Properties inherit their + job's owner, so pass the account the jobs were created under -- an organization's, when + working on its behalf. Defaults to the caller's personal account. Returns: - The best-precision total_energy property, or None if none exists. + The best-precision property holder, or None if none exists. """ material = client.materials.get(material_id) exabyte_id = material.get("exabyteId") if not exabyte_id: return None - query = {"exabyteId": exabyte_id, "slug": "total_energy"} + query = {"exabyteId": exabyte_id, "slug": property_name} if source == "curators": query["owner.slug"] = "curators" elif source == "my_account": - query["owner._id"] = client.my_account.id + query["owner._id"] = owner_id or client.my_account.id elif source != "public": raise ValueError(f"Invalid source: {source!r}. Expected 'public', 'curators', or 'my_account'.") properties = client.properties.list( @@ -108,6 +118,11 @@ def find_total_energy_for_material(client: APIClient, material_id: str, source: return properties[0] if properties else None +def find_total_energy_for_material(client: APIClient, material_id: str, source: str = "my_account") -> Optional[dict]: + """Find the best-precision total_energy property for a material -- `find_property_for_material` for the details.""" + return find_property_for_material(client, material_id, "total_energy", source) + + def get_property_by_subworkflow_and_unit_indicies( endpoint: PropertiesEndpoints, property_name: str, job: dict, subworkflow_index: int, unit_index: int ) -> dict: diff --git a/tests/py/unit/core/entity/test_property_api_find_total_energy_for_material.py b/tests/py/unit/core/entity/test_property_api.py similarity index 75% rename from tests/py/unit/core/entity/test_property_api_find_total_energy_for_material.py rename to tests/py/unit/core/entity/test_property_api.py index a7596ca7..3731e2da 100644 --- a/tests/py/unit/core/entity/test_property_api_find_total_energy_for_material.py +++ b/tests/py/unit/core/entity/test_property_api.py @@ -1,9 +1,9 @@ -"""Unit tests for find_total_energy_for_material.""" +"""Unit tests for core.entity.property.api.""" from unittest.mock import MagicMock import pytest -from mat3ra.notebooks_utils.core.entity.property.api import find_total_energy_for_material +from mat3ra.notebooks_utils.core.entity.property.api import find_property_for_material, find_total_energy_for_material MATERIAL_ID = "material-a" EXABYTE_ID = "exabyte-a" @@ -108,3 +108,34 @@ def test_find_total_energy_for_material_returns_none_when_material_has_no_exabyt client.properties.list.assert_not_called() assert result is None + + +@pytest.mark.parametrize("property_name", ["total_energy", "band_gaps"]) +def test_find_property_for_material_queries_by_property_name(property_name): + client = _client() + + find_property_for_material(client, MATERIAL_ID, property_name) + + client.properties.list.assert_called_once_with( + query={ + "exabyteId": EXABYTE_ID, + "slug": property_name, + "owner._id": OWNER_ACCOUNT_ID, + }, + projection={"sort": {"precision.value": -1}, "limit": 1}, + ) + + +def test_find_property_for_material_scopes_my_account_to_an_explicit_owner(): + client = _client() + + find_property_for_material(client, MATERIAL_ID, "band_gaps", owner_id="org-account") + + client.properties.list.assert_called_once_with( + query={ + "exabyteId": EXABYTE_ID, + "slug": "band_gaps", + "owner._id": "org-account", + }, + projection={"sort": {"precision.value": -1}, "limit": 1}, + )