Merge pull request #294 from Exabyte-io/feature/SOF-7777

feature/SOF 7777

Author: VsevolodX

other/materials_designer/workflows/valence_band_offset.ipynb

@@ -9,10 +9,17 @@
     "\n",
     "Calculate the valence band offset for a labeled interface using a DFT workflow on the Mat3ra platform.\n",
     "\n",
+    "The notebook supports two modes controlled by `IS_POLAR` in cell 1.2:\n",
+    "\n",
+    "- `IS_POLAR = False`: run the standard VBO workflow and report the workflow VBO together with the average ESP profiles.\n",
+    "- `IS_POLAR = True`: keep the same workflow VBO calculation, then run an additional notebook-side post-process at the end that fits the interface and bulk ESP profiles over the slab regions and plots the polar fit.\n",
+    "\n",
+    "When the polar option is enabled, the workflow itself is unchanged. The extra polar analysis is performed in the final notebook cells using the saved interface, left-slab, and right-slab materials plus the ESP results returned by the job.\n",
+    "\n",
     "<h2 style=\"color:green\">Usage</h2>\n",
     "\n",
     "1. Create and save an interface with labels (for example via `create_interface_with_min_strain_zsl.ipynb`).\n",
-    "1. Set the interface and calculation parameters in cells 1.2 and 1.3 below.\n",
+    "1. Set the interface and calculation parameters in cells 1.2 and 1.3 below, including `IS_POLAR` if the interface is polar.\n",
     "1. Click \"Run\" > \"Run All\" to run all cells.\n",
     "1. Wait for the job to complete.\n",
     "1. Scroll down to view the VBO result.\n",
@@ -26,7 +33,8 @@
     "1. Configure compute: get list of clusters and create compute configuration with selected cluster, queue, and number of processors.\n",
     "1. Create the job with materials and workflow configuration: assemble the job from materials, workflow, project, and compute configuration.\n",
     "1. Submit the job and monitor the status: submit the job and wait for completion.\n",
-    "1. Retrieve results: get and display the valence band offset.\n"
+    "1. Retrieve results: get and display the workflow valence band offset and average ESP profiles.\n",
+    "1. If `IS_POLAR = True`, run the final polar post-process cells to fit the ESP profile and generate the polar-fit plot in the notebook.\n"
    ]
   },
   {
@@ -87,6 +95,8 @@
     "RIGHT_SIDE_PART = InterfacePartsEnum.FILM\n",
     "INTERFACE_SYSTEM_NAME = None  # Used as tag to group the materials. Defaults to shorthand from the loaded interface name\n",
     "\n",
+    "IS_POLAR = False  # Whether the interface is polar, to adjust the VBO calculation method accordingly.\n",
+    "\n",
     "# 4. Workflow parameters\n",
     "APPLICATION_NAME = \"espresso\"\n",
     "WORKFLOW_SEARCH_TERM = \"valence_band_offset.json\"\n",
@@ -370,6 +380,7 @@
     "\n",
     "workflow_config = WorkflowStandata.filter_by_application(app.name).get_by_name_first_match(WORKFLOW_SEARCH_TERM)\n",
     "workflow = Workflow.create(workflow_config)\n",
+    "\n",
     "workflow.name = MY_WORKFLOW_NAME\n",
     "\n",
     "visualize_workflow(workflow)\n"
@@ -630,11 +641,11 @@
     "job_data = client.jobs.get(job_id)\n",
     "workflow = job_data[\"workflow\"]\n",
     "\n",
-    "vbo_value = client.properties.get_for_job(\n",
+    "vbo_value_non_polar = client.properties.get_for_job(\n",
     "    job_id,\n",
     "    property_name=PropertyName.scalar.valence_band_offset.value,\n",
     ")[0]\n",
-    "print(f\"Valence Band Offset (VBO) value: {vbo_value['value']:.3f} eV\")\n",
+    "print(f\"Workflow Valence Band Offset (VBO) value: {float(vbo_value_non_polar['value']):.3f} eV\")\n",
     "\n",
     "avg_esp_unit_ids = {}\n",
     "for subworkflow in workflow[\"subworkflows\"]:\n",
@@ -643,31 +654,181 @@
     "        result_names = [result[\"name\"] for result in unit.get(\"results\", [])]\n",
     "        if \"average_potential_profile\" in result_names:\n",
     "            avg_esp_unit_ids[subworkflow_name] = unit[\"flowchartId\"]\n",
-    "            break\n",
     "\n",
     "ordered_names = [\n",
     "    \"BS + Avg ESP (Interface)\",\n",
     "    \"BS + Avg ESP (interface left)\",\n",
     "    \"BS + Avg ESP (interface right)\",\n",
     "]\n",
     "\n",
+    "avg_esp_results = {}\n",
     "for subworkflow_name in ordered_names:\n",
     "    unit_id = avg_esp_unit_ids[subworkflow_name]\n",
     "    avg_esp_data = client.properties.get_for_job(\n",
     "        job_id,\n",
     "        property_name=\"average_potential_profile\",\n",
     "        unit_id=unit_id,\n",
     "    )[0]\n",
-    "    visualize_properties(avg_esp_data, title=subworkflow_name)\n"
+    "    avg_esp_results[subworkflow_name] = avg_esp_data\n",
+    "    visualize_properties(avg_esp_data, title=subworkflow_name)\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "43",
+   "metadata": {},
+   "source": [
+    "### 8.2. Polar-only notebook post-process\n",
+    "If `IS_POLAR = True`, use the slab boundaries from the notebook materials and fit the ESP profiles directly in the notebook.\n"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "43",
+   "id": "44",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "from scipy.stats import linregress\n",
+    "from mat3ra.prode import PropertyName\n",
+    "from utils.api import get_property_holder_for_job, update_property_holder_value\n",
+    "\n",
+    "\n",
+    "def get_region_indices(x_data, x_min, x_max):\n",
+    "    mask = (x_data >= x_min) & (x_data <= x_max)\n",
+    "    indices = np.where(mask)[0]\n",
+    "    if len(indices) == 0:\n",
+    "        return 0, len(x_data)\n",
+    "    return indices[0], indices[-1] + 1\n",
+    "\n",
+    "\n",
+    "def fit_and_average(x_data, y_data, start_idx, end_idx):\n",
+    "    x_region = x_data[start_idx:end_idx]\n",
+    "    y_region = y_data[start_idx:end_idx]\n",
+    "    if len(x_region) < 2:\n",
+    "        avg = float(np.mean(y_region)) if len(y_region) > 0 else 0.0\n",
+    "        return avg, 0.0, avg\n",
+    "    slope, intercept, _, _, _ = linregress(x_region, y_region)\n",
+    "    x_mid = (x_region[0] + x_region[-1]) / 2.0\n",
+    "    avg_value = slope * x_mid + intercept\n",
+    "    return float(avg_value), float(slope), float(intercept)\n",
+    "\n",
+    "\n",
+    "interface_profile = avg_esp_results[\"BS + Avg ESP (Interface)\"]\n",
+    "left_profile = avg_esp_results[\"BS + Avg ESP (interface left)\"]\n",
+    "right_profile = avg_esp_results[\"BS + Avg ESP (interface right)\"]\n",
+    "\n",
+    "interface_material.to_cartesian()\n",
+    "left_material.to_cartesian()\n",
+    "right_material.to_cartesian()\n",
+    "\n",
+    "interface_z = sorted(coord[2] for coord in interface_material.basis.coordinates.values)\n",
+    "n_left = len(left_material.basis.elements.values)\n",
+    "z_min_1 = interface_z[0]\n",
+    "z_max_1 = interface_z[n_left - 1]\n",
+    "z_min_2 = interface_z[n_left]\n",
+    "z_max_2 = interface_z[-1]\n",
+    "\n",
+    "print(f\"Detected slab 1 boundaries: z = {z_min_1:.3f} to {z_max_1:.3f} A\")\n",
+    "print(f\"Detected slab 2 boundaries: z = {z_min_2:.3f} to {z_max_2:.3f} A\")\n",
+    "\n",
+    "X = np.array(interface_profile[\"xDataArray\"])\n",
+    "Y = np.array(interface_profile[\"yDataSeries\"][1])\n",
+    "X_left = np.array(left_profile[\"xDataArray\"])\n",
+    "Y_left = np.array(left_profile[\"yDataSeries\"][1])\n",
+    "X_right = np.array(right_profile[\"xDataArray\"])\n",
+    "Y_right = np.array(right_profile[\"yDataSeries\"][1])\n",
+    "\n",
+    "slab1_start, slab1_end = get_region_indices(X, z_min_1, z_max_1)\n",
+    "slab2_start, slab2_end = get_region_indices(X, z_min_2, z_max_2)\n",
+    "slab1_start_left, slab1_end_left = get_region_indices(X_left, z_min_1, z_max_1)\n",
+    "slab2_start_right, slab2_end_right = get_region_indices(X_right, z_min_2, z_max_2)\n",
+    "\n",
+    "Va_interface, slope_a, intercept_a = fit_and_average(X, Y, slab1_start, slab1_end)\n",
+    "Vb_interface, slope_b, intercept_b = fit_and_average(X, Y, slab2_start, slab2_end)\n",
+    "Va_bulk_left, _, _ = fit_and_average(X_left, Y_left, slab1_start_left, slab1_end_left)\n",
+    "Vb_bulk_right, _, _ = fit_and_average(X_right, Y_right, slab2_start_right, slab2_end_right)\n",
+    "\n",
+    "vbo_value_polar = (Vb_interface - Va_interface) - (Vb_bulk_right - Va_bulk_left)\n",
+    "\n",
+    "print(f\"Interface ESP slab 1: {Va_interface:.3f} eV\")\n",
+    "print(f\"Interface ESP slab 2: {Vb_interface:.3f} eV\")\n",
+    "print(f\"Bulk ESP left: {Va_bulk_left:.3f} eV\")\n",
+    "print(f\"Bulk ESP right: {Vb_bulk_right:.3f} eV\")\n",
+    "print(f\"Fit slope slab 1: {slope_a:.6f} eV/A\")\n",
+    "print(f\"Fit slope slab 2: {slope_b:.6f} eV/A\")\n",
+    "print(f\"Polar post-process VBO: {vbo_value_polar:.3f} eV\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "45",
+   "metadata": {},
+   "source": [
+    "### 8.3. Plot the ESP profiles with slab regions and fits"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "46",
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "plt.figure(figsize=(10, 6))\n",
+    "plt.plot(X, Y, label=\"Macroscopic Average Potential\", linewidth=2)\n",
+    "plt.axvspan(z_min_1, z_max_1, color=\"red\", alpha=0.2, label=\"Slab 1 Region\")\n",
+    "plt.axvspan(z_min_2, z_max_2, color=\"blue\", alpha=0.2, label=\"Slab 2 Region\")\n",
+    "\n",
+    "if slab1_end > slab1_start:\n",
+    "    x_fit1 = X[slab1_start:slab1_end]\n",
+    "    y_fit1 = slope_a * x_fit1 + intercept_a\n",
+    "    plt.plot(x_fit1, y_fit1, color=\"darkred\", linestyle=\"--\", linewidth=2, label=\"Fit Slab 1\")\n",
+    "\n",
+    "if slab2_end > slab2_start:\n",
+    "    x_fit2 = X[slab2_start:slab2_end]\n",
+    "    y_fit2 = slope_b * x_fit2 + intercept_b\n",
+    "    plt.plot(x_fit2, y_fit2, color=\"darkblue\", linestyle=\"--\", linewidth=2, label=\"Fit Slab 2\")\n",
+    "\n",
+    "plt.axhline(Va_interface, color=\"red\", linestyle=\":\", linewidth=2, label=f\"Avg ESP Slab 1 = {Va_interface:.3f} eV\")\n",
+    "plt.axhline(Vb_interface, color=\"blue\", linestyle=\":\", linewidth=2, label=f\"Avg ESP Slab 2 = {Vb_interface:.3f} eV\")\n",
+    "plt.xlabel(\"z-coordinate (A)\", fontsize=12)\n",
+    "plt.ylabel(\"Macroscopic Average Potential (eV)\", fontsize=12)\n",
+    "plt.title(f\"Polar Interface VBO = {vbo_value_polar:.3f} eV\", fontsize=14, fontweight=\"bold\")\n",
+    "plt.legend(loc=\"best\", fontsize=10)\n",
+    "plt.grid(True, alpha=0.3)\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "47",
+   "metadata": {},
+   "source": [
+    "### 8.4. Update the VBO property with the polar post-process result"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "48",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "if IS_POLAR:\n",
+    "    vbo_property_holder = get_property_holder_for_job(\n",
+    "        client,\n",
+    "        job_id,\n",
+    "        PropertyName.scalar.valence_band_offset.value,\n",
+    "    )\n",
+    "    update_property_holder_value(client, vbo_property_holder[\"_id\"], vbo_value_polar)\n",
+    "    print(f\"Persisted displayed Valence Band Offset (VBO) value: {vbo_value_polar:.3f} eV\")\n",
+    "\n"
+   ]
   }
  ],
  "metadata": {

utils/api.py

@@ -200,6 +200,48 @@ def get_properties_for_job(client: APIClient, job_id: str, property_name: Option
     return [{**prop, "fermiEnergy": fermi_energy} for prop in properties]
 
 
+def get_property_holder_for_job(
+    client: APIClient, job_id: str, property_name: str, unit_id: Optional[str] = None
+) -> dict:
+    """
+    Fetch the first full property holder for a job/property pair.
+
+    Args:
+        client (APIClient): API client instance.
+        job_id (str): Job ID.
+        property_name (str): Property name.
+        unit_id (str, optional): Unit flowchart ID.
+
+    Returns:
+        dict: Full property holder document.
+    """
+    query = {
+        "source.info.jobId": job_id,
+        "data.name": property_name,
+    }
+    if unit_id:
+        query["source.info.unitId"] = unit_id
+    holders = client.properties.list(query=query)
+    if not holders:
+        raise ValueError(f"Property '{property_name}' not found for job '{job_id}'")
+    return holders[0]
+
+
+def update_property_holder_value(client: APIClient, property_holder_id: str, value: float) -> dict:
+    """
+    Update a scalar property's data.value.
+
+    Args:
+        client (APIClient): API client instance.
+        property_holder_id (str): Property holder ID.
+        value (float): New scalar value.
+
+    Returns:
+        dict: Server response payload.
+    """
+    return client.properties.update(property_holder_id, {"$set": {"data.value": value}})
+
+
 def create_job(
     api_client: APIClient,
     materials: List[Union[dict, Material]],