From 2400cbadf4b6f3fc84e2213a673e5c87a8590add Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Fri, 2 Oct 2026 17:01:38 -0700 Subject: [PATCH 01/17] =?UTF-8?q?SOF-8064:=20name=20the=20seven=20Gr/h-BN?= =?UTF-8?q?=20interfaces=20and=20re-set=20the=20cell=20to=20the=20120?= =?UTF-8?q?=C2=B0=20hexagonal=20setting?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The platform resolves the symbolic K point from a per-lattice-type table in the 120° convention; the ZSL cell comes out at 60°, where [1/3, 1/3, 0] is not K, so the cell is re-set and typed HEX. Co-Authored-By: Claude Opus 5.5 (1M context) --- ...terface_2d_2d_boron_nitride_graphene.ipynb | 55 +++++++++++++++++-- 1 file changed, 49 insertions(+), 6 deletions(-) diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb index 48b18a1a2..cb33e3b1b 100644 --- a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb @@ -8,7 +8,7 @@ "\n", "## Introduction\n", "\n", - "This notebook demonstrates the creation of an interface between two 2D materials, Boron Nitride (BN) and Graphene and shifting the film along the y-axis to create multiple (7) stacking configurations.\n", + "This notebook demonstrates the creation of an interface between two 2D materials, Boron Nitride (BN) and Graphene and shifting the film along the y-axis to create multiple (7) stacking configurations. The three symmetric stackings carry the labels of Jung et al. (2015) and Giovannetti et al. (2007): AA (C over B and N), AB (C over N and a hexagon centre), BA (C over B and a hexagon centre).\n", "\n", "Following the manuscript:\n", "> **Jeil Jung, Ashley M. DaSilva, Allan H. MacDonald & Shaffique Adam**\n", @@ -62,7 +62,21 @@ "MAX_ANGLE_TOLERANCE = 0.02\n", "\n", "# Whether to reduce the resulting interface cell to the primitive cell after the interface creation.\n", - "REDUCE_RESULT_CELL_TO_PRIMITIVE = True" + "REDUCE_RESULT_CELL_TO_PRIMITIVE = True\n", + "\n", + "# Names of the seven shifted interfaces, in the order the shift loop builds them. The registry of\n", + "# each was measured from the built cell: n=2 is BA (one C over B, the other over a hexagon centre),\n", + "# n=4 is AA (C over B and C over N), n=6 is AB (C over N, the other over a hexagon centre); the odd\n", + "# n are bridge positions; n=8 repeats n=2 one period later.\n", + "INTERFACE_NAMES = [\n", + " \"Gr/hBN d3.4 shift 0of6 BA\",\n", + " \"Gr/hBN d3.4 shift 1of6\",\n", + " \"Gr/hBN d3.4 shift 2of6 AA\",\n", + " \"Gr/hBN d3.4 shift 3of6\",\n", + " \"Gr/hBN d3.4 shift 4of6 AB\",\n", + " \"Gr/hBN d3.4 shift 5of6\",\n", + " \"Gr/hBN d3.4 shift 6of6 BA\",\n", + "]" ] }, { @@ -338,6 +352,27 @@ ")" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 3.5. Set the cell to the standard hexagonal setting\n", + "The ZSL interface cell comes out with γ = 60°. The symbolic K point the band-structure notebook uses assumes the standard 120° hexagonal cell, so the cell is re-set here (same four atoms) and each shifted interface is typed HEX." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.made.tools.helpers import create_supercell\n", + "\n", + "interface = create_supercell(interface, supercell_matrix=[[1, 0, 0], [-1, 1, 0], [0, 0, 1]])\n", + "print(f\"{len(interface.basis.elements.ids)} atoms, a = {interface.lattice.a:.4f} Å, \"\n", + " f\"gamma = {interface.lattice.gamma:.1f}°\")" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -372,16 +407,24 @@ "outputs": [], "source": [ "import numpy as np\n", + "from mat3ra.made.tools.analyze.other import get_average_interlayer_distance\n", + "from mat3ra.made.tools.convert.interface_parts_enum import InterfacePartsEnum\n", "from mat3ra.made.tools.modify import interface_displace_part\n", "\n", "a = interface.lattice.a\n", "shifted_interfaces = []\n", - "for n in range(2, 9):\n", + "for index, n in enumerate(range(2, 9)):\n", " shifted_interface = interface_displace_part(\n", " interface=interface,\n", " displacement=[0, n * a / np.sqrt(3) / 2, 0],\n", " use_cartesian_coordinates=True)\n", - " shifted_interfaces.append(shifted_interface)" + " shifted_interface.name = INTERFACE_NAMES[index]\n", + " shifted_interface.lattice.type = \"HEX\"\n", + " shifted_interfaces.append(shifted_interface)\n", + " interlayer_distance = get_average_interlayer_distance(\n", + " shifted_interface, InterfacePartsEnum.SUBSTRATE.value, InterfacePartsEnum.FILM.value)\n", + " print(f\"{shifted_interface.name}: {len(shifted_interface.basis.elements.ids)} atoms, \"\n", + " f\"gamma = {shifted_interface.lattice.gamma:.1f}°, interlayer distance = {interlayer_distance:.3f} Å\")" ] }, { @@ -417,8 +460,8 @@ "from mat3ra.notebooks_utils.material import set_materials\n", "\n", "set_materials(shifted_interfaces)\n", - "for idx, shifted_interface in enumerate(shifted_interfaces):\n", - " download_content_to_file(shifted_interface.to_json(), f\"interface_{idx}.json\")" + "for shifted_interface in shifted_interfaces:\n", + " download_content_to_file(shifted_interface.to_json(), f\"{shifted_interface.name}.json\")" ] } ], From 719f573550f7df50daf78e5c126f01dc27fac07e Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Fri, 2 Oct 2026 17:06:00 -0700 Subject: [PATCH 02/17] SOF-8064: shorter comment on the interface names Co-Authored-By: Claude Opus 5.5 (1M context) --- .../interface_2d_2d_boron_nitride_graphene.ipynb | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb index cb33e3b1b..eb28b2de6 100644 --- a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb @@ -64,10 +64,8 @@ "# Whether to reduce the resulting interface cell to the primitive cell after the interface creation.\n", "REDUCE_RESULT_CELL_TO_PRIMITIVE = True\n", "\n", - "# Names of the seven shifted interfaces, in the order the shift loop builds them. The registry of\n", - "# each was measured from the built cell: n=2 is BA (one C over B, the other over a hexagon centre),\n", - "# n=4 is AA (C over B and C over N), n=6 is AB (C over N, the other over a hexagon centre); the odd\n", - "# n are bridge positions; n=8 repeats n=2 one period later.\n", + "# One name per shift (n = 2..8); the three symmetric registries carry their Jung 2015 / Giovannetti 2007 label,\n", + "# n = 8 repeats n = 2 one period later.\n", "INTERFACE_NAMES = [\n", " \"Gr/hBN d3.4 shift 0of6 BA\",\n", " \"Gr/hBN d3.4 shift 1of6\",\n", From bad76e175b11ccbe6d002e11b69094023e9bc7cc Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Fri, 2 Oct 2026 17:22:10 -0700 Subject: [PATCH 03/17] =?UTF-8?q?SOF-8064:=20simulation=20notebook=20?= =?UTF-8?q?=E2=80=94=20stacking=20energy=20and=20gap=20at=20K=20for=20the?= =?UTF-8?q?=20seven=20Gr/h-BN=20interfaces?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Opus 5.5 (1M context) --- .../specific_examples/Introduction.ipynb | 2 +- ...2d_boron_nitride_graphene_SIMULATION.ipynb | 724 ++++++++++++++++++ 2 files changed, 725 insertions(+), 1 deletion(-) create mode 100644 other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb diff --git a/other/materials_designer/specific_examples/Introduction.ipynb b/other/materials_designer/specific_examples/Introduction.ipynb index 71075fbc4..156267f34 100644 --- a/other/materials_designer/specific_examples/Introduction.ipynb +++ b/other/materials_designer/specific_examples/Introduction.ipynb @@ -24,7 +24,7 @@ "| `P-0D-NRB` | Nanoribbon | *To be added* | — | — |\n", "| `C-2D-HST` | Heterostack | [Si/SiO₂/HfO₂/TiN Heterostructure](heterostructure_silicon_silicon_dioxide_hafnium_dioxide_titanium_nitride.ipynb) | *To be added* | [[3]](#ref3) |\n", "| `C-2D-INT-S` | Interface Simple | *To be added* | — | — |\n", - "| `C-2D-INT-Z` | Interface ZSL | [BN/Graphene 2D–2D Interface](interface_2d_2d_boron_nitride_graphene.ipynb) | *To be added* | [[4]](#ref4) |\n", + "| `C-2D-INT-Z` | Interface ZSL | [BN/Graphene 2D–2D Interface](interface_2d_2d_boron_nitride_graphene.ipynb) | [Gr/h-BN Stacking Energy and Band Gap](interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb) | [[4]](#ref4) |\n", "| `C-2D-INT-Z` | Interface ZSL | [Graphene/SiO₂ 2D–3D Interface](interface_2d_3d_graphene_silicon_dioxide.ipynb) | *To be added* | [[5]](#ref5) |\n", "| `C-2D-INT-Z` | Interface ZSL | [Cu/Cristobalite 3D–3D Interface](interface_3d_3d_copper_cristobalite.ipynb) | *To be added* | [[6]](#ref6) |\n", "| `C-2D-INT-Z` | Interface ZSL | [Graphene/Ni Interface Film XY Position Optimization](optimization_interface_film_xy_position_graphene_nickel.ipynb) | [Gr/Ni(111) Registry and Work of Adhesion](optimization_interface_film_xy_position_graphene_nickel_SIMULATION.ipynb) | [[7]](#ref7) |\n", diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb new file mode 100644 index 000000000..756b8e270 --- /dev/null +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb @@ -0,0 +1,724 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Stacking Energy and Band Gap of Graphene on h-BN\n", + "\n", + "> **Gianluca Giovannetti, Petr A. Khomyakov, Geert Brocks, Paul J. Kelly & Jeroen van den Brink**\n", + "> Substrate-induced band gap in graphene on hexagonal boron nitride: Ab initio density functional calculations. Physical Review B, 76, 073103. 2007.\n", + "> [https://doi.org/10.1103/PhysRevB.76.073103](https://doi.org/10.1103/PhysRevB.76.073103)\n", + "\n", + "Calculate the total energy and the band structure of the seven graphene/h-BN stackings at d = 3.4 Å\n", + "created in the [structure notebook](interface_2d_2d_boron_nitride_graphene.ipynb), using Quantum\n", + "ESPRESSO and the band structure + density of states workflow from Standata, and compare the\n", + "ordering of the sliding energies and the gap at K with Giovannetti et al. (2007).\n", + "\n", + "The seven stackings follow the sliding path of Fig. 7(a) in Jung et al. (2015). That figure is a model\n", + "curve parameterised from RPA calculations; the DFT numbers compared here are Giovannetti's.\n", + "\n", + "

Usage

\n", + "\n", + "1. Create the materials in the [structure notebook](interface_2d_2d_boron_nitride_graphene.ipynb), which saves them to the `uploads` folder under the names used in cell 1.2 below.\n", + "1. Set the materials and the calculation parameters in cells 1.2 and 1.3, or keep the default values: everything runs with them.\n", + "1. Click \"Run\" > \"Run All\" to run all cells.\n", + "1. Wait for the jobs to complete.\n", + "1. Scroll down to view the results; the comparison cell in section 8 prints the verdict.\n", + "\n", + "## Summary\n", + "\n", + "1. Set up the environment and parameters: install packages (JupyterLite only) and configure the materials, workflow, model, compute resources and jobs.\n", + "1. Authenticate and initialize API client: authenticate via browser, initialize the client, then select account and project.\n", + "1. Load the materials by name from the `uploads` folder, print their provenance and save them to the platform.\n", + "1. Configure the workflow: select the application and the model, load the band structure + DOS workflow from Standata, and set the computational parameters for each material.\n", + "1. Configure compute: get the list of clusters and create a compute configuration.\n", + "1. Run one job per material, one after another, re-using a job that already ran under the same workflow name.\n", + "1. Retrieve results: band structures, total energies and the gap at K.\n", + "1. Compare with Giovannetti et al. (2007)." + ] + }, + { + "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 and configurations for the workflow and jobs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "from datetime import datetime\n", + "from mat3ra.ide.compute import QueueName\n", + "\n", + "ORGANIZATION_NAME = None # set to your organization name (full or partial); otherwise, your default one is used\n", + "FOLDER = \"./uploads\"\n", + "\n", + "# Names saved by the structure notebook; the symmetric stackings carry their Jung 2015 / Giovannetti 2007 label\n", + "MATERIALS = {\n", + " \"Gr/hBN d3.4 shift 0of6 BA\": {\"shift\": 0, \"stacking\": \"BA\"},\n", + " \"Gr/hBN d3.4 shift 1of6\": {\"shift\": 1, \"stacking\": None},\n", + " \"Gr/hBN d3.4 shift 2of6 AA\": {\"shift\": 2, \"stacking\": \"AA\"},\n", + " \"Gr/hBN d3.4 shift 3of6\": {\"shift\": 3, \"stacking\": None},\n", + " \"Gr/hBN d3.4 shift 4of6 AB\": {\"shift\": 4, \"stacking\": \"AB\"},\n", + " \"Gr/hBN d3.4 shift 5of6\": {\"shift\": 5, \"stacking\": None},\n", + " \"Gr/hBN d3.4 shift 6of6 BA\": {\"shift\": 6, \"stacking\": \"BA\"},\n", + "}\n", + "\n", + "WORKFLOW_SEARCH_TERM = \"band_structure_dos.json\"\n", + "MY_WORKFLOW_NAME = \"Band Structure + DOS\"\n", + "APPLICATION_NAME = \"espresso\"\n", + "\n", + "CLUSTER_NAME = \"001\" # specify full or partial name i.e. \"cluster-001\" to select\n", + "QUEUE_NAME = QueueName.OR\n", + "PPN = 16 # queue OR on cluster-001 allows at most 16 cores per node\n", + "TIME_LIMIT = \"02:00:00\"\n", + "\n", + "timestamp = datetime.now().strftime(\"%Y-%m-%d %H:%M\")\n", + "POLL_INTERVAL = 60 # seconds" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "### 1.3. Set the DFT parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "MODEL_SUBTYPE = \"lda\"\n", + "FUNCTIONAL = \"pz\" # Giovannetti et al. 2007 use LDA: GGA gives essentially no interlayer binding\n", + "PSEUDOPOTENTIAL_TYPE = \"us\" # GBRV ultrasoft, the only LDA family the platform publishes for B, C and N\n", + "ECUTWFC = 50 # Ry\n", + "ECUTRHO = 400 # Ry, 8x for ultrasoft pseudopotentials\n", + "\n", + "KGRID = [36, 36, 1] # Giovannetti et al. 2007; a multiple of 3 keeps K on the mesh\n", + "SMEARING_SETTINGS = {\"degauss\": 0.001} # Ry; the gaps compared are 30-80 meV\n", + "MODEL_TAG = f\"{FUNCTIONAL}-{PSEUDOPOTENTIAL_TYPE} {ECUTWFC}-{ECUTRHO}Ry k{KGRID[0]} g{SMEARING_SETTINGS['degauss']}\"\n", + "\n", + "SCF_UNIT = \"pw_scf\"\n", + "NSCF_UNIT = \"pw_nscf\"\n", + "BANDS_UNIT = \"pw_bands\"\n", + "N_OCCUPIED_BANDS = 8 # 16 valence electrons: C 4 + 4, B 3, N 5\n", + "\n", + "KPATH_STEPS = 40\n", + "KPATH = [\n", + " {\"point\": \"Γ\", \"steps\": KPATH_STEPS},\n", + " {\"point\": \"K\", \"steps\": KPATH_STEPS},\n", + " {\"point\": \"M\", \"steps\": KPATH_STEPS},\n", + " {\"point\": \"Γ\", \"steps\": 1},\n", + "]" + ] + }, + { + "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" + ] + }, + { + "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. Load the materials\n", + "### 3.1. Load from the uploads folder or the platform, and print provenance\n", + "\n", + "The structures, their geometry and their lattice type all come from the structure notebook; this one\n", + "only loads them by name." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17", + "metadata": {}, + "outputs": [], + "source": [ + "from collections import Counter\n", + "from mat3ra.made.tools.analyze.other import get_average_interlayer_distance\n", + "from mat3ra.made.tools.convert.interface_parts_enum import InterfacePartsEnum\n", + "from mat3ra.notebooks_utils.core.entity.material.api import load_material\n", + "\n", + "materials = {}\n", + "for name, settings in MATERIALS.items():\n", + " material = load_material(client, FOLDER, name, ACCOUNT_ID)\n", + " materials[name] = material\n", + " composition = \"\".join(\n", + " f\"{element}{count}\" for element, count in sorted(Counter(material.basis.elements.values).items()))\n", + " interlayer_distance = get_average_interlayer_distance(\n", + " material, InterfacePartsEnum.SUBSTRATE.value, InterfacePartsEnum.FILM.value)\n", + " print(f\"{name}: {composition}, {material.basis.number_of_atoms} atoms, a = {material.lattice.a:.4f} Å, \"\n", + " f\"gamma = {material.lattice.gamma:.3f}°, interlayer distance = {interlayer_distance:.3f} Å, \"\n", + " f\"{settings['stacking'] or 'bridge'}\")" + ] + }, + { + "cell_type": "markdown", + "id": "18", + "metadata": {}, + "source": [ + "### 3.2. Preview the materials" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.ipython.entity.material.visualize import visualize_materials\n", + "\n", + "visualize_materials([{\"material\": material, \"title\": name} for name, material in materials.items()],\n", + " viewer=\"wave\")" + ] + }, + { + "cell_type": "markdown", + "id": "20", + "metadata": {}, + "source": [ + "### 3.3. Save the materials to the platform" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.made.material import Material\n", + "from mat3ra.notebooks_utils.core.entity.material.api import get_or_create_material\n", + "\n", + "saved_materials = {}\n", + "for name, material in materials.items():\n", + " material.basis.set_labels_from_list([])\n", + " saved_materials[name] = Material.create(get_or_create_material(client, material, ACCOUNT_ID))" + ] + }, + { + "cell_type": "markdown", + "id": "22", + "metadata": {}, + "source": [ + "## 4. Configure the workflow\n", + "### 4.1. Select application and model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.standata.applications import ApplicationStandata\n", + "from mat3ra.standata.model_tree import ModelTreeStandata\n", + "from mat3ra.ade.application import Application\n", + "from mat3ra.mode import ModelFactory\n", + "\n", + "app_config = ApplicationStandata.get_by_name_first_match(APPLICATION_NAME)\n", + "app = Application(**app_config)\n", + "\n", + "model_config = ModelTreeStandata.get_model_by_parameters(type=\"dft\", subtype=MODEL_SUBTYPE, functional=FUNCTIONAL)\n", + "model_config[\"method\"] = {\"type\": \"pseudopotential\", \"subtype\": PSEUDOPOTENTIAL_TYPE}\n", + "model = ModelFactory.create(model_config)\n", + "print(f\"Using application: {app.name}, model: {MODEL_TAG}\")" + ] + }, + { + "cell_type": "markdown", + "id": "24", + "metadata": {}, + "source": [ + "### 4.2. Configure the workflow\n", + "\n", + "One workflow per material, all with the settings of 1.3. Each workflow is named after its material\n", + "and `MODEL_TAG`; section 6 finds a job that already ran by that name." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "25", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.standata.workflows import WorkflowStandata\n", + "from mat3ra.wode.workflows import Workflow\n", + "from mat3ra.wode.context.providers import PlanewaveCutoffsContextProvider, PointsGridDataProvider, \\\n", + " PointsPathDataProvider\n", + "from mat3ra.notebooks_utils.workflow import patch_workflow_qe_input\n", + "\n", + "from copy import deepcopy\n", + "\n", + "workflow_config = WorkflowStandata.filter_by_application(app.name).get_by_name_first_match(WORKFLOW_SEARCH_TERM)\n", + "\n", + "cutoffs_context = PlanewaveCutoffsContextProvider(\n", + " wavefunction=ECUTWFC, density=ECUTRHO, isEdited=True).get_context_item_data()\n", + "grid_context = PointsGridDataProvider(\n", + " material=next(iter(materials.values())), dimensions=KGRID, isEdited=True).get_context_item_data()\n", + "path_context = PointsPathDataProvider(path=KPATH, isEdited=True).get_context_item_data()\n", + "\n", + "workflows = {}\n", + "for name in MATERIALS:\n", + " workflow = Workflow.create(deepcopy(workflow_config))\n", + " workflow.name = f\"{MY_WORKFLOW_NAME} {name} {MODEL_TAG}\"\n", + " subworkflow = workflow.subworkflows[0]\n", + " subworkflow.model = model\n", + "\n", + " for unit_name, contexts in [(SCF_UNIT, [grid_context, cutoffs_context]),\n", + " (NSCF_UNIT, [grid_context, cutoffs_context]),\n", + " (BANDS_UNIT, [path_context, cutoffs_context])]:\n", + " unit = subworkflow.get_unit_by_name(name=unit_name)\n", + " for context in contexts:\n", + " unit.add_context(context)\n", + " subworkflow.set_unit(unit)\n", + " patch_workflow_qe_input(workflow, {\"system\": SMEARING_SETTINGS}, [SCF_UNIT, NSCF_UNIT, BANDS_UNIT])\n", + " workflows[name] = workflow\n", + " print(workflow.name)" + ] + }, + { + "cell_type": "markdown", + "id": "26", + "metadata": {}, + "source": [ + "### 4.3. Preview the workflow" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "27", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.ipython.entity.workflow.visualize import visualize_workflow\n", + "\n", + "visualize_workflow(workflows[next(iter(MATERIALS))])" + ] + }, + { + "cell_type": "markdown", + "id": "28", + "metadata": {}, + "source": [ + "### 4.4. Save the workflows to the collection" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "29", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.core.entity.workflow.api import get_or_create_workflow\n", + "\n", + "for name, workflow in workflows.items():\n", + " saved = Workflow.create(get_or_create_workflow(client, workflow, ACCOUNT_ID))\n", + " print(f\"{name}: workflow {saved.id}\")" + ] + }, + { + "cell_type": "markdown", + "id": "30", + "metadata": {}, + "source": [ + "## 5. Create the compute configuration\n", + "### 5.1. Get list of clusters" + ] + }, + { + "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 the compute configuration for the jobs" + ] + }, + { + "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", + " if cluster is None:\n", + " raise ValueError(f\"Cluster '{CLUSTER_NAME}' not found. Available: {[c['hostname'] for c in clusters]}\")\n", + "else:\n", + " cluster = clusters[0]\n", + "compute = Compute(cluster=cluster, queue=QUEUE_NAME, ppn=PPN, timeLimit=TIME_LIMIT)\n", + "print(f\"Using cluster: {compute.cluster.hostname}, queue: {QUEUE_NAME}, ppn: {PPN}, \"\n", + " f\"time limit: {TIME_LIMIT}\")" + ] + }, + { + "cell_type": "markdown", + "id": "34", + "metadata": {}, + "source": [ + "## 6. Run the jobs one at a time\n", + "\n", + "A job is found by its material and its workflow name, which carries the model tag, and re-used if it\n", + "exists; otherwise it is created and submitted. The jobs run one after another." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "35", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.core.entity.job.api import find_job_for_material\n", + "from mat3ra.notebooks_utils.job import create_job\n", + "from mat3ra.notebooks_utils.api.job import submit_jobs, wait_for_jobs_to_finish_async\n", + "\n", + "job_ids = {}\n", + "for name, workflow in workflows.items():\n", + " job = find_job_for_material(\n", + " client, saved_materials[name].id, workflow.name, ACCOUNT_ID,\n", + " statuses=(\"submitted\", \"queued\", \"active\", \"finished\"),\n", + " )\n", + " if job is None:\n", + " job = create_job(\n", + " api_client=client, materials=[saved_materials[name]], workflow=workflow, project_id=project_id,\n", + " owner_id=ACCOUNT_ID, compute=compute.to_dict(), prefix=f\"{workflow.name} {timestamp}\",\n", + " )\n", + " submit_jobs(client.jobs, [job[\"_id\"]])\n", + " print(f\"✅ {name}: submitted job {job['_id']}\")\n", + " else:\n", + " print(f\"♻️ {name}: reusing job {job['_id']}\")\n", + " await wait_for_jobs_to_finish_async(client.jobs, [job[\"_id\"]], poll_interval=POLL_INTERVAL)\n", + " job_ids[name] = job[\"_id\"]" + ] + }, + { + "cell_type": "markdown", + "id": "36", + "metadata": {}, + "source": [ + "## 7. Retrieve the results\n", + "### 7.1. Band structures" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "37", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.core.entity.property.api import get_properties_for_job\n", + "from mat3ra.notebooks_utils.ipython.entity.property.visualize import visualize_properties\n", + "\n", + "band_structures = {}\n", + "for name, job_id in job_ids.items():\n", + " band_structures[name] = get_properties_for_job(client, job_id, property_name=\"band_structure\")\n", + " visualize_properties(band_structures[name], title=f\"Band Structure: {name}\",\n", + " extra_config={\"material\": materials[name].to_dict()})" + ] + }, + { + "cell_type": "markdown", + "id": "38", + "metadata": {}, + "source": [ + "### 7.2. Total energies" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "39", + "metadata": {}, + "outputs": [], + "source": [ + "total_energies = {}\n", + "for name, job_id in job_ids.items():\n", + " total_energies[name] = get_properties_for_job(client, job_id, \"total_energy\")[0][\"value\"]\n", + " print(f\"{name}: {total_energies[name]:.6f} eV\")" + ] + }, + { + "cell_type": "markdown", + "id": "40", + "metadata": {}, + "source": [ + "### 7.3. Gap at K\n", + "\n", + "The gap is the smallest difference along the path between the lowest unoccupied and the highest\n", + "occupied band; the k-point where it occurs is printed in crystal coordinates, K being (1/3, 1/3, 0)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "41", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "gaps = {}\n", + "for name, band_structure in band_structures.items():\n", + " bands = np.array(band_structure[0][\"yDataSeries\"])\n", + " gaps_along_path = bands[N_OCCUPIED_BANDS] - bands[N_OCCUPIED_BANDS - 1]\n", + " kpoint_index = int(np.argmin(gaps_along_path))\n", + " gaps[name] = 1000 * gaps_along_path[kpoint_index]\n", + " kpoint = band_structure[0][\"xDataArray\"][kpoint_index]\n", + " print(f\"{name}: {gaps[name]:.1f} meV at k = ({kpoint[0]:.4f}, {kpoint[1]:.4f}, {kpoint[2]:.4f})\")" + ] + }, + { + "cell_type": "markdown", + "id": "42", + "metadata": {}, + "source": [ + "## 8. Compare with Giovannetti et al. (2007)\n", + "\n", + "ΔE is the total energy relative to the first entry of `MATERIALS`. The paper's gaps are read off its\n", + "Fig. 4 at d = 3.4 Å; the gaps at its own equilibrium distances are printed for reference." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.plot import plot_series\n", + "\n", + "GIOVANNETTI_GAP_AT_3P4 = {\"AA\": 80, \"AB\": 45, \"BA\": 30} # meV, Fig. 4 read off at d = 3.4 Å, ±5 meV\n", + "GIOVANNETTI_GAP_AT_EQUILIBRIUM = {\"AA\": (3.50, 56), \"AB\": (3.40, 46), \"BA\": (3.22, 53)} # (Å, meV), not compared\n", + "GAP_TOLERANCE = 15 # meV\n", + "\n", + "reference_name = next(iter(MATERIALS))\n", + "series = [\n", + " {\n", + " \"shift\": settings[\"shift\"],\n", + " \"stacking\": settings[\"stacking\"] or \"bridge\",\n", + " \"energy\": 1000 * (total_energies[name] - total_energies[reference_name]),\n", + " \"gap\": gaps[name],\n", + " }\n", + " for name, settings in MATERIALS.items()\n", + "]\n", + "symmetric_stackings = {}\n", + "for item in series:\n", + " if item[\"stacking\"] in GIOVANNETTI_GAP_AT_3P4:\n", + " symmetric_stackings.setdefault(item[\"stacking\"], item)\n", + "\n", + "print(f\"{'shift':>5} {'stacking':>8} {'ΔE (meV)':>9} {'gap (meV)':>10} {'paper gap (meV)':>16}\")\n", + "for item in series:\n", + " print(f\"{item['shift']:>5} {item['stacking']:>8} {item['energy']:>9.2f} {item['gap']:>10.1f} \"\n", + " f\"{GIOVANNETTI_GAP_AT_3P4.get(item['stacking'], ''):>16}\")\n", + "\n", + "plot_series(series=series, x_key=\"shift\", y_key=\"energy\", xlabel=\"Shift along y (sixths of √3 a)\",\n", + " ylabel=\"ΔE (meV per cell)\", title=\"Total energy along the sliding path\")\n", + "plot_series(series=series, x_key=\"shift\", y_key=\"gap\", xlabel=\"Shift along y (sixths of √3 a)\",\n", + " ylabel=\"Gap at K (meV)\", title=\"Gap at K along the sliding path\")\n", + "\n", + "energy_ordering = (symmetric_stackings[\"BA\"][\"energy\"] < symmetric_stackings[\"AB\"][\"energy\"]\n", + " < symmetric_stackings[\"AA\"][\"energy\"])\n", + "energy_extrema = (min(series, key=lambda item: item[\"energy\"])[\"stacking\"] == \"BA\"\n", + " and max(series, key=lambda item: item[\"energy\"])[\"stacking\"] == \"AA\")\n", + "gap_ordering = symmetric_stackings[\"AA\"][\"gap\"] > symmetric_stackings[\"AB\"][\"gap\"] > symmetric_stackings[\"BA\"][\"gap\"]\n", + "gaps_within_tolerance = all(abs(symmetric_stackings[stacking][\"gap\"] - paper_gap) <= GAP_TOLERANCE\n", + " for stacking, paper_gap in GIOVANNETTI_GAP_AT_3P4.items())\n", + "verdict = \"yes\" if energy_ordering and energy_extrema and gap_ordering and gaps_within_tolerance else \"no\"\n", + "\n", + "print(f\"E(BA) < E(AB) < E(AA): {energy_ordering}\")\n", + "print(f\"Minimum at BA, maximum at AA: {energy_extrema}\")\n", + "print(f\"Gap AA > AB > BA: {gap_ordering}\")\n", + "print(f\"Each gap within {GAP_TOLERANCE} meV of Fig. 4: {gaps_within_tolerance}\")\n", + "print(\"Gaps at the paper's equilibrium distances: \" + \", \".join(\n", + " f\"{stacking} {distance:.2f} Å {gap} meV\" for stacking, (distance, gap) in GIOVANNETTI_GAP_AT_EQUILIBRIUM.items()))\n", + "print(\"Offsets from the paper's setup:\")\n", + "print(f\" cell a = {materials[reference_name].lattice.a:.3f} Å vs the paper's 2.445 Å\")\n", + "print(\" one h-BN layer vs four\")\n", + "print(\" GBRV ultrasoft LDA pseudopotentials vs PAW\")\n", + "print(f\" Gaussian smearing of {SMEARING_SETTINGS['degauss']} Ry vs the tetrahedron method\")\n", + "print(\" no dipole correction\")\n", + "print(f\"Reproduces Giovannetti et al. (2007): {verdict}\")" + ] + }, + { + "cell_type": "markdown", + "id": "44", + "metadata": {}, + "source": [ + "## References\n", + "\n", + "[1] G. Giovannetti, P. A. Khomyakov, G. Brocks, P. J. Kelly and J. van den Brink, \"Substrate-induced band gap in graphene on hexagonal boron nitride: Ab initio density functional calculations\", Phys. Rev. B 76, 073103 (2007). https://doi.org/10.1103/PhysRevB.76.073103\n", + "\n", + "[2] J. Jung, A. M. DaSilva, A. H. MacDonald and S. Adam, \"Origin of band gaps in graphene on hexagonal boron nitride\", Nat. Commun. 6, 6308 (2015). https://doi.org/10.1038/ncomms7308" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From aed923576ac07dee1c5927b08cd6bc9170a69a99 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Fri, 2 Oct 2026 17:28:29 -0700 Subject: [PATCH 04/17] =?UTF-8?q?SOF-8064:=20comparison=20cell=20=E2=80=94?= =?UTF-8?q?=20x-axis=20label=20and=20the=20pseudopotential=20offset=20line?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Sonnet 5 --- .../interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb index 756b8e270..5efc59fd1 100644 --- a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb @@ -658,9 +658,9 @@ " print(f\"{item['shift']:>5} {item['stacking']:>8} {item['energy']:>9.2f} {item['gap']:>10.1f} \"\n", " f\"{GIOVANNETTI_GAP_AT_3P4.get(item['stacking'], ''):>16}\")\n", "\n", - "plot_series(series=series, x_key=\"shift\", y_key=\"energy\", xlabel=\"Shift along y (sixths of √3 a)\",\n", + "plot_series(series=series, x_key=\"shift\", y_key=\"energy\", xlabel=\"Shift from BA along y (steps of a/(2√3))\",\n", " ylabel=\"ΔE (meV per cell)\", title=\"Total energy along the sliding path\")\n", - "plot_series(series=series, x_key=\"shift\", y_key=\"gap\", xlabel=\"Shift along y (sixths of √3 a)\",\n", + "plot_series(series=series, x_key=\"shift\", y_key=\"gap\", xlabel=\"Shift from BA along y (steps of a/(2√3))\",\n", " ylabel=\"Gap at K (meV)\", title=\"Gap at K along the sliding path\")\n", "\n", "energy_ordering = (symmetric_stackings[\"BA\"][\"energy\"] < symmetric_stackings[\"AB\"][\"energy\"]\n", @@ -681,7 +681,7 @@ "print(\"Offsets from the paper's setup:\")\n", "print(f\" cell a = {materials[reference_name].lattice.a:.3f} Å vs the paper's 2.445 Å\")\n", "print(\" one h-BN layer vs four\")\n", - "print(\" GBRV ultrasoft LDA pseudopotentials vs PAW\")\n", + "print(\" GBRV ultrasoft LDA pseudopotentials vs the paper's VASP potentials at 600 eV\")\n", "print(f\" Gaussian smearing of {SMEARING_SETTINGS['degauss']} Ry vs the tetrahedron method\")\n", "print(\" no dipole correction\")\n", "print(f\"Reproduces Giovannetti et al. (2007): {verdict}\")" From 4f3f498c09bce51b5cfe63756c9bea8cf22966a5 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Fri, 2 Oct 2026 18:09:38 -0700 Subject: [PATCH 05/17] =?UTF-8?q?SOF-8064:=20review=20=E2=80=94=20GBRV=204?= =?UTF-8?q?0/200=20Ry,=20full=20names,=20path=20steps=20in=20the=20reuse?= =?UTF-8?q?=20tag,=20constants=20where=20used?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Opus 5.5 (1M context) --- ...2d_boron_nitride_graphene_SIMULATION.ipynb | 48 ++++++++++--------- 1 file changed, 26 insertions(+), 22 deletions(-) diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb index 5efc59fd1..a49f504d2 100644 --- a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb @@ -14,7 +14,7 @@ "Calculate the total energy and the band structure of the seven graphene/h-BN stackings at d = 3.4 Å\n", "created in the [structure notebook](interface_2d_2d_boron_nitride_graphene.ipynb), using Quantum\n", "ESPRESSO and the band structure + density of states workflow from Standata, and compare the\n", - "ordering of the sliding energies and the gap at K with Giovannetti et al. (2007).\n", + "ordering of the sliding energies and the direct gap at K with Giovannetti et al. (2007).\n", "\n", "The seven stackings follow the sliding path of Fig. 7(a) in Jung et al. (2015). That figure is a model\n", "curve parameterised from RPA calculations; the DFT numbers compared here are Giovannetti's.\n", @@ -123,19 +123,19 @@ "MODEL_SUBTYPE = \"lda\"\n", "FUNCTIONAL = \"pz\" # Giovannetti et al. 2007 use LDA: GGA gives essentially no interlayer binding\n", "PSEUDOPOTENTIAL_TYPE = \"us\" # GBRV ultrasoft, the only LDA family the platform publishes for B, C and N\n", - "ECUTWFC = 50 # Ry\n", - "ECUTRHO = 400 # Ry, 8x for ultrasoft pseudopotentials\n", + "ECUTWFC = 40 # Ry, GBRV's tested cutoff\n", + "ECUTRHO = 200 # Ry, GBRV's tested charge-density cutoff\n", "\n", "KGRID = [36, 36, 1] # Giovannetti et al. 2007; a multiple of 3 keeps K on the mesh\n", "SMEARING_SETTINGS = {\"degauss\": 0.001} # Ry; the gaps compared are 30-80 meV\n", - "MODEL_TAG = f\"{FUNCTIONAL}-{PSEUDOPOTENTIAL_TYPE} {ECUTWFC}-{ECUTRHO}Ry k{KGRID[0]} g{SMEARING_SETTINGS['degauss']}\"\n", + "KPATH_STEPS = 40\n", + "MODEL_TAG = f\"{FUNCTIONAL}-{PSEUDOPOTENTIAL_TYPE} {ECUTWFC}-{ECUTRHO}Ry k{KGRID[0]} p{KPATH_STEPS} g{SMEARING_SETTINGS['degauss']}\"\n", "\n", "SCF_UNIT = \"pw_scf\"\n", "NSCF_UNIT = \"pw_nscf\"\n", "BANDS_UNIT = \"pw_bands\"\n", - "N_OCCUPIED_BANDS = 8 # 16 valence electrons: C 4 + 4, B 3, N 5\n", + "NUMBER_OF_OCCUPIED_BANDS = 8 # 16 valence electrons: C 4 + 4, B 3, N 5\n", "\n", - "KPATH_STEPS = 40\n", "KPATH = [\n", " {\"point\": \"Γ\", \"steps\": KPATH_STEPS},\n", " {\"point\": \"K\", \"steps\": KPATH_STEPS},\n", @@ -247,7 +247,7 @@ "metadata": {}, "source": [ "## 3. Load the materials\n", - "### 3.1. Load from the uploads folder or the platform, and print provenance\n", + "### 3.1. Load from the uploads folder and print provenance\n", "\n", "The structures, their geometry and their lattice type all come from the structure notebook; this one\n", "only loads them by name." @@ -589,10 +589,11 @@ "id": "40", "metadata": {}, "source": [ - "### 7.3. Gap at K\n", + "### 7.3. Direct gap along the path\n", "\n", - "The gap is the smallest difference along the path between the lowest unoccupied and the highest\n", - "occupied band; the k-point where it occurs is printed in crystal coordinates, K being (1/3, 1/3, 0)." + "The smallest difference along Γ–K–M–Γ between the lowest unoccupied and the highest occupied band;\n", + "the k-point where it occurs is printed in crystal coordinates. For the three symmetric stackings it\n", + "sits at K = (1/3, 1/3, 0)." ] }, { @@ -607,7 +608,7 @@ "gaps = {}\n", "for name, band_structure in band_structures.items():\n", " bands = np.array(band_structure[0][\"yDataSeries\"])\n", - " gaps_along_path = bands[N_OCCUPIED_BANDS] - bands[N_OCCUPIED_BANDS - 1]\n", + " gaps_along_path = bands[NUMBER_OF_OCCUPIED_BANDS] - bands[NUMBER_OF_OCCUPIED_BANDS - 1]\n", " kpoint_index = int(np.argmin(gaps_along_path))\n", " gaps[name] = 1000 * gaps_along_path[kpoint_index]\n", " kpoint = band_structure[0][\"xDataArray\"][kpoint_index]\n", @@ -621,8 +622,9 @@ "source": [ "## 8. Compare with Giovannetti et al. (2007)\n", "\n", - "ΔE is the total energy relative to the first entry of `MATERIALS`. The paper's gaps are read off its\n", - "Fig. 4 at d = 3.4 Å; the gaps at its own equilibrium distances are printed for reference." + "ΔE is the total energy relative to the first entry of `MATERIALS`. The paper's AA and BA gaps are read\n", + "off its Fig. 4 at d = 3.4 Å, and AB is the paper's value at its 3.40 Å equilibrium; the gaps at its\n", + "own equilibrium distances are printed for reference." ] }, { @@ -634,8 +636,10 @@ "source": [ "from mat3ra.notebooks_utils.plot import plot_series\n", "\n", - "GIOVANNETTI_GAP_AT_3P4 = {\"AA\": 80, \"AB\": 45, \"BA\": 30} # meV, Fig. 4 read off at d = 3.4 Å, ±5 meV\n", - "GIOVANNETTI_GAP_AT_EQUILIBRIUM = {\"AA\": (3.50, 56), \"AB\": (3.40, 46), \"BA\": (3.22, 53)} # (Å, meV), not compared\n", + "# meV; AA and BA read off Fig. 4 at d = 3.4 Å (±5 meV), AB is the paper's value at its 3.40 Å equilibrium\n", + "GIOVANNETTI_GAPS_AT_3_4_ANGSTROM = {\"AA\": 80, \"AB\": 46, \"BA\": 30}\n", + "GIOVANNETTI_GAPS_AT_EQUILIBRIUM = {\"AA\": (3.50, 56), \"AB\": (3.40, 46), \"BA\": (3.22, 53)}\n", + "GIOVANNETTI_LATTICE_CONSTANT = 2.445 # Å, graphene LDA\n", "GAP_TOLERANCE = 15 # meV\n", "\n", "reference_name = next(iter(MATERIALS))\n", @@ -650,18 +654,18 @@ "]\n", "symmetric_stackings = {}\n", "for item in series:\n", - " if item[\"stacking\"] in GIOVANNETTI_GAP_AT_3P4:\n", + " if item[\"stacking\"] in GIOVANNETTI_GAPS_AT_3_4_ANGSTROM:\n", " symmetric_stackings.setdefault(item[\"stacking\"], item)\n", "\n", "print(f\"{'shift':>5} {'stacking':>8} {'ΔE (meV)':>9} {'gap (meV)':>10} {'paper gap (meV)':>16}\")\n", "for item in series:\n", " print(f\"{item['shift']:>5} {item['stacking']:>8} {item['energy']:>9.2f} {item['gap']:>10.1f} \"\n", - " f\"{GIOVANNETTI_GAP_AT_3P4.get(item['stacking'], ''):>16}\")\n", + " f\"{GIOVANNETTI_GAPS_AT_3_4_ANGSTROM.get(item['stacking'], ''):>16}\")\n", "\n", "plot_series(series=series, x_key=\"shift\", y_key=\"energy\", xlabel=\"Shift from BA along y (steps of a/(2√3))\",\n", " ylabel=\"ΔE (meV per cell)\", title=\"Total energy along the sliding path\")\n", "plot_series(series=series, x_key=\"shift\", y_key=\"gap\", xlabel=\"Shift from BA along y (steps of a/(2√3))\",\n", - " ylabel=\"Gap at K (meV)\", title=\"Gap at K along the sliding path\")\n", + " ylabel=\"Direct gap along the path (meV)\", title=\"Gap at K along the sliding path\")\n", "\n", "energy_ordering = (symmetric_stackings[\"BA\"][\"energy\"] < symmetric_stackings[\"AB\"][\"energy\"]\n", " < symmetric_stackings[\"AA\"][\"energy\"])\n", @@ -669,17 +673,17 @@ " and max(series, key=lambda item: item[\"energy\"])[\"stacking\"] == \"AA\")\n", "gap_ordering = symmetric_stackings[\"AA\"][\"gap\"] > symmetric_stackings[\"AB\"][\"gap\"] > symmetric_stackings[\"BA\"][\"gap\"]\n", "gaps_within_tolerance = all(abs(symmetric_stackings[stacking][\"gap\"] - paper_gap) <= GAP_TOLERANCE\n", - " for stacking, paper_gap in GIOVANNETTI_GAP_AT_3P4.items())\n", + " for stacking, paper_gap in GIOVANNETTI_GAPS_AT_3_4_ANGSTROM.items())\n", "verdict = \"yes\" if energy_ordering and energy_extrema and gap_ordering and gaps_within_tolerance else \"no\"\n", "\n", "print(f\"E(BA) < E(AB) < E(AA): {energy_ordering}\")\n", "print(f\"Minimum at BA, maximum at AA: {energy_extrema}\")\n", "print(f\"Gap AA > AB > BA: {gap_ordering}\")\n", - "print(f\"Each gap within {GAP_TOLERANCE} meV of Fig. 4: {gaps_within_tolerance}\")\n", + "print(f\"Each gap within {GAP_TOLERANCE} meV of the paper: {gaps_within_tolerance}\")\n", "print(\"Gaps at the paper's equilibrium distances: \" + \", \".join(\n", - " f\"{stacking} {distance:.2f} Å {gap} meV\" for stacking, (distance, gap) in GIOVANNETTI_GAP_AT_EQUILIBRIUM.items()))\n", + " f\"{stacking} {distance:.2f} Å {gap} meV\" for stacking, (distance, gap) in GIOVANNETTI_GAPS_AT_EQUILIBRIUM.items()))\n", "print(\"Offsets from the paper's setup:\")\n", - "print(f\" cell a = {materials[reference_name].lattice.a:.3f} Å vs the paper's 2.445 Å\")\n", + "print(f\" cell a = {materials[reference_name].lattice.a:.3f} Å vs the paper's {GIOVANNETTI_LATTICE_CONSTANT} Å\")\n", "print(\" one h-BN layer vs four\")\n", "print(\" GBRV ultrasoft LDA pseudopotentials vs the paper's VASP potentials at 600 eV\")\n", "print(f\" Gaussian smearing of {SMEARING_SETTINGS['degauss']} Ry vs the tetrahedron method\")\n", From 42016fedefa557ab558e1398d295c9fe924ad9e8 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Fri, 2 Oct 2026 18:10:58 -0700 Subject: [PATCH 06/17] =?UTF-8?q?SOF-8064:=20two=20labels=20=E2=80=94=20di?= =?UTF-8?q?rect=20gap=20along=20the=20path?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Sonnet 5 --- .../interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb index a49f504d2..a780ba8a3 100644 --- a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb @@ -35,7 +35,7 @@ "1. Configure the workflow: select the application and the model, load the band structure + DOS workflow from Standata, and set the computational parameters for each material.\n", "1. Configure compute: get the list of clusters and create a compute configuration.\n", "1. Run one job per material, one after another, re-using a job that already ran under the same workflow name.\n", - "1. Retrieve results: band structures, total energies and the gap at K.\n", + "1. Retrieve results: band structures, total energies and the direct gap along the path.\n", "1. Compare with Giovannetti et al. (2007)." ] }, @@ -665,7 +665,7 @@ "plot_series(series=series, x_key=\"shift\", y_key=\"energy\", xlabel=\"Shift from BA along y (steps of a/(2√3))\",\n", " ylabel=\"ΔE (meV per cell)\", title=\"Total energy along the sliding path\")\n", "plot_series(series=series, x_key=\"shift\", y_key=\"gap\", xlabel=\"Shift from BA along y (steps of a/(2√3))\",\n", - " ylabel=\"Direct gap along the path (meV)\", title=\"Gap at K along the sliding path\")\n", + " ylabel=\"Direct gap along the path (meV)\", title=\"Direct gap along the sliding path\")\n", "\n", "energy_ordering = (symmetric_stackings[\"BA\"][\"energy\"] < symmetric_stackings[\"AB\"][\"energy\"]\n", " < symmetric_stackings[\"AA\"][\"energy\"])\n", From 15e3381b84d2dc67683a3d31539c3e5060b96b35 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 5 Oct 2026 15:55:43 -0700 Subject: [PATCH 07/17] =?UTF-8?q?SOF-8064:=20gap=20tolerance=20relative=20?= =?UTF-8?q?=E2=80=94=2015=20%=20of=20the=20paper's=20value?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Sonnet 5.5 --- ..._2d_2d_boron_nitride_graphene_SIMULATION.ipynb | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb index a780ba8a3..a904c07d3 100644 --- a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb @@ -640,7 +640,7 @@ "GIOVANNETTI_GAPS_AT_3_4_ANGSTROM = {\"AA\": 80, \"AB\": 46, \"BA\": 30}\n", "GIOVANNETTI_GAPS_AT_EQUILIBRIUM = {\"AA\": (3.50, 56), \"AB\": (3.40, 46), \"BA\": (3.22, 53)}\n", "GIOVANNETTI_LATTICE_CONSTANT = 2.445 # Å, graphene LDA\n", - "GAP_TOLERANCE = 15 # meV\n", + "GAP_TOLERANCE_FRACTION = 0.15 # of the paper's value\n", "\n", "reference_name = next(iter(MATERIALS))\n", "series = [\n", @@ -657,10 +657,12 @@ " if item[\"stacking\"] in GIOVANNETTI_GAPS_AT_3_4_ANGSTROM:\n", " symmetric_stackings.setdefault(item[\"stacking\"], item)\n", "\n", - "print(f\"{'shift':>5} {'stacking':>8} {'ΔE (meV)':>9} {'gap (meV)':>10} {'paper gap (meV)':>16}\")\n", + "print(f\"{'shift':>5} {'stacking':>8} {'ΔE (meV)':>9} {'gap (meV)':>10} {'paper gap (meV)':>16} {'deviation':>10}\")\n", "for item in series:\n", + " paper_gap = GIOVANNETTI_GAPS_AT_3_4_ANGSTROM.get(item[\"stacking\"])\n", + " deviation = f\"{(item['gap'] - paper_gap) / paper_gap:+.0%}\" if paper_gap else \"\"\n", " print(f\"{item['shift']:>5} {item['stacking']:>8} {item['energy']:>9.2f} {item['gap']:>10.1f} \"\n", - " f\"{GIOVANNETTI_GAPS_AT_3_4_ANGSTROM.get(item['stacking'], ''):>16}\")\n", + " f\"{paper_gap or '':>16} {deviation:>10}\")\n", "\n", "plot_series(series=series, x_key=\"shift\", y_key=\"energy\", xlabel=\"Shift from BA along y (steps of a/(2√3))\",\n", " ylabel=\"ΔE (meV per cell)\", title=\"Total energy along the sliding path\")\n", @@ -672,14 +674,15 @@ "energy_extrema = (min(series, key=lambda item: item[\"energy\"])[\"stacking\"] == \"BA\"\n", " and max(series, key=lambda item: item[\"energy\"])[\"stacking\"] == \"AA\")\n", "gap_ordering = symmetric_stackings[\"AA\"][\"gap\"] > symmetric_stackings[\"AB\"][\"gap\"] > symmetric_stackings[\"BA\"][\"gap\"]\n", - "gaps_within_tolerance = all(abs(symmetric_stackings[stacking][\"gap\"] - paper_gap) <= GAP_TOLERANCE\n", - " for stacking, paper_gap in GIOVANNETTI_GAPS_AT_3_4_ANGSTROM.items())\n", + "gaps_within_tolerance = all(\n", + " abs(symmetric_stackings[stacking][\"gap\"] - paper_gap) <= GAP_TOLERANCE_FRACTION * paper_gap\n", + " for stacking, paper_gap in GIOVANNETTI_GAPS_AT_3_4_ANGSTROM.items())\n", "verdict = \"yes\" if energy_ordering and energy_extrema and gap_ordering and gaps_within_tolerance else \"no\"\n", "\n", "print(f\"E(BA) < E(AB) < E(AA): {energy_ordering}\")\n", "print(f\"Minimum at BA, maximum at AA: {energy_extrema}\")\n", "print(f\"Gap AA > AB > BA: {gap_ordering}\")\n", - "print(f\"Each gap within {GAP_TOLERANCE} meV of the paper: {gaps_within_tolerance}\")\n", + "print(f\"Each gap within {100 * GAP_TOLERANCE_FRACTION:.0f} % of the paper's: {gaps_within_tolerance}\")\n", "print(\"Gaps at the paper's equilibrium distances: \" + \", \".join(\n", " f\"{stacking} {distance:.2f} Å {gap} meV\" for stacking, (distance, gap) in GIOVANNETTI_GAPS_AT_EQUILIBRIUM.items()))\n", "print(\"Offsets from the paper's setup:\")\n", From a0e1e485c067e486a5311a5e4745bab697894f45 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 5 Oct 2026 17:03:55 -0700 Subject: [PATCH 08/17] =?UTF-8?q?SOF-8064:=20comparison=20cell=20reports,?= =?UTF-8?q?=20it=20does=20not=20grade=20=E2=80=94=20the=20SnO=20shape?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Sonnet 5.5 --- ...2d_boron_nitride_graphene_SIMULATION.ipynb | 36 +++---------------- 1 file changed, 5 insertions(+), 31 deletions(-) diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb index a904c07d3..a625363e7 100644 --- a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb @@ -25,7 +25,7 @@ "1. Set the materials and the calculation parameters in cells 1.2 and 1.3, or keep the default values: everything runs with them.\n", "1. Click \"Run\" > \"Run All\" to run all cells.\n", "1. Wait for the jobs to complete.\n", - "1. Scroll down to view the results; the comparison cell in section 8 prints the verdict.\n", + "1. Scroll down to view the results; section 8 prints the comparison table and plots.\n", "\n", "## Summary\n", "\n", @@ -624,7 +624,9 @@ "\n", "ΔE is the total energy relative to the first entry of `MATERIALS`. The paper's AA and BA gaps are read\n", "off its Fig. 4 at d = 3.4 Å, and AB is the paper's value at its 3.40 Å equilibrium; the gaps at its\n", - "own equilibrium distances are printed for reference." + "own equilibrium distances are printed for reference.\n", + "\n", + "Settings that differ from the paper's: cell a = 2.509 Å (h-BN unstrained, graphene +1.79 %) against 2.445 Å; one h-BN layer against four; GBRV ultrasoft pseudopotentials against VASP at 600 eV; Gaussian smearing against the tetrahedron method; no dipole correction." ] }, { @@ -639,8 +641,6 @@ "# meV; AA and BA read off Fig. 4 at d = 3.4 Å (±5 meV), AB is the paper's value at its 3.40 Å equilibrium\n", "GIOVANNETTI_GAPS_AT_3_4_ANGSTROM = {\"AA\": 80, \"AB\": 46, \"BA\": 30}\n", "GIOVANNETTI_GAPS_AT_EQUILIBRIUM = {\"AA\": (3.50, 56), \"AB\": (3.40, 46), \"BA\": (3.22, 53)}\n", - "GIOVANNETTI_LATTICE_CONSTANT = 2.445 # Å, graphene LDA\n", - "GAP_TOLERANCE_FRACTION = 0.15 # of the paper's value\n", "\n", "reference_name = next(iter(MATERIALS))\n", "series = [\n", @@ -652,11 +652,6 @@ " }\n", " for name, settings in MATERIALS.items()\n", "]\n", - "symmetric_stackings = {}\n", - "for item in series:\n", - " if item[\"stacking\"] in GIOVANNETTI_GAPS_AT_3_4_ANGSTROM:\n", - " symmetric_stackings.setdefault(item[\"stacking\"], item)\n", - "\n", "print(f\"{'shift':>5} {'stacking':>8} {'ΔE (meV)':>9} {'gap (meV)':>10} {'paper gap (meV)':>16} {'deviation':>10}\")\n", "for item in series:\n", " paper_gap = GIOVANNETTI_GAPS_AT_3_4_ANGSTROM.get(item[\"stacking\"])\n", @@ -669,29 +664,8 @@ "plot_series(series=series, x_key=\"shift\", y_key=\"gap\", xlabel=\"Shift from BA along y (steps of a/(2√3))\",\n", " ylabel=\"Direct gap along the path (meV)\", title=\"Direct gap along the sliding path\")\n", "\n", - "energy_ordering = (symmetric_stackings[\"BA\"][\"energy\"] < symmetric_stackings[\"AB\"][\"energy\"]\n", - " < symmetric_stackings[\"AA\"][\"energy\"])\n", - "energy_extrema = (min(series, key=lambda item: item[\"energy\"])[\"stacking\"] == \"BA\"\n", - " and max(series, key=lambda item: item[\"energy\"])[\"stacking\"] == \"AA\")\n", - "gap_ordering = symmetric_stackings[\"AA\"][\"gap\"] > symmetric_stackings[\"AB\"][\"gap\"] > symmetric_stackings[\"BA\"][\"gap\"]\n", - "gaps_within_tolerance = all(\n", - " abs(symmetric_stackings[stacking][\"gap\"] - paper_gap) <= GAP_TOLERANCE_FRACTION * paper_gap\n", - " for stacking, paper_gap in GIOVANNETTI_GAPS_AT_3_4_ANGSTROM.items())\n", - "verdict = \"yes\" if energy_ordering and energy_extrema and gap_ordering and gaps_within_tolerance else \"no\"\n", - "\n", - "print(f\"E(BA) < E(AB) < E(AA): {energy_ordering}\")\n", - "print(f\"Minimum at BA, maximum at AA: {energy_extrema}\")\n", - "print(f\"Gap AA > AB > BA: {gap_ordering}\")\n", - "print(f\"Each gap within {100 * GAP_TOLERANCE_FRACTION:.0f} % of the paper's: {gaps_within_tolerance}\")\n", "print(\"Gaps at the paper's equilibrium distances: \" + \", \".join(\n", - " f\"{stacking} {distance:.2f} Å {gap} meV\" for stacking, (distance, gap) in GIOVANNETTI_GAPS_AT_EQUILIBRIUM.items()))\n", - "print(\"Offsets from the paper's setup:\")\n", - "print(f\" cell a = {materials[reference_name].lattice.a:.3f} Å vs the paper's {GIOVANNETTI_LATTICE_CONSTANT} Å\")\n", - "print(\" one h-BN layer vs four\")\n", - "print(\" GBRV ultrasoft LDA pseudopotentials vs the paper's VASP potentials at 600 eV\")\n", - "print(f\" Gaussian smearing of {SMEARING_SETTINGS['degauss']} Ry vs the tetrahedron method\")\n", - "print(\" no dipole correction\")\n", - "print(f\"Reproduces Giovannetti et al. (2007): {verdict}\")" + " f\"{stacking} {distance:.2f} Å {gap} meV\" for stacking, (distance, gap) in GIOVANNETTI_GAPS_AT_EQUILIBRIUM.items()))" ] }, { From 0558cfc88beec2be071693aec3490543caa21429 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Mon, 5 Oct 2026 19:09:46 -0700 Subject: [PATCH 09/17] =?UTF-8?q?SOF-8064:=20queue=20D,=20one=20core=20?= =?UTF-8?q?=E2=80=94=20a=20four-atom=20cell=20needs=20seconds=20of=20QE?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Sonnet 5.5 --- .../interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb index a625363e7..392a370b3 100644 --- a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb @@ -97,9 +97,9 @@ "APPLICATION_NAME = \"espresso\"\n", "\n", "CLUSTER_NAME = \"001\" # specify full or partial name i.e. \"cluster-001\" to select\n", - "QUEUE_NAME = QueueName.OR\n", - "PPN = 16 # queue OR on cluster-001 allows at most 16 cores per node\n", - "TIME_LIMIT = \"02:00:00\"\n", + "QUEUE_NAME = QueueName.D\n", + "PPN = 1\n", + "TIME_LIMIT = \"01:00:00\"\n", "\n", "timestamp = datetime.now().strftime(\"%Y-%m-%d %H:%M\")\n", "POLL_INTERVAL = 60 # seconds" From 6e22d2af52553445107de3e275d1850221f3be20 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Wed, 7 Oct 2026 18:58:35 -0700 Subject: [PATCH 10/17] =?UTF-8?q?SOF-8064:=20structure=20notebook=20?= =?UTF-8?q?=E2=80=94=20Giovannetti's=20cell:=20graphene=20on=204-layer=20A?= =?UTF-8?q?A'=20h-BN=20at=20a=20=3D=202.445=20=C3=85,=20three=20registries?= =?UTF-8?q?,=20distance=20list?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Opus 5.5 --- ...terface_2d_2d_boron_nitride_graphene.ipynb | 364 +++++------------- 1 file changed, 98 insertions(+), 266 deletions(-) diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb index eb28b2de6..ac49af613 100644 --- a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb @@ -8,18 +8,13 @@ "\n", "## Introduction\n", "\n", - "This notebook demonstrates the creation of an interface between two 2D materials, Boron Nitride (BN) and Graphene and shifting the film along the y-axis to create multiple (7) stacking configurations. The three symmetric stackings carry the labels of Jung et al. (2015) and Giovannetti et al. (2007): AA (C over B and N), AB (C over N and a hexagon centre), BA (C over B and a hexagon centre).\n", + "This notebook creates graphene on four layers of hexagonal boron nitride (h-BN) in a common hexagonal cell, for three registries of graphene on the top h-BN layer and a list of graphene–h-BN distances.\n", "\n", "Following the manuscript:\n", - "> **Jeil Jung, Ashley M. DaSilva, Allan H. MacDonald & Shaffique Adam**\n", - "> **Origin of the band gap in graphene on hexagonal boron nitride**\n", - "> Nature Communications volume 6, Article number: 6308 (2015)\n", - "> [DOI: 10.1038/ncomms7308](https://doi.org/10.1038/ncomms7308)\n", - "\n", - "\n", - "Relicating the materials with profile give in Figure 7. a (top row):\n", - "\n", - "\n" + "> **Gianluca Giovannetti, Petr A. Khomyakov, Geert Brocks, Paul J. Kelly, and Jeroen van den Brink**\n", + "> **Substrate-induced band gap in graphene on hexagonal boron nitride: Ab initio density functional calculations**\n", + "> Physical Review B 76, 073103 (2007)\n", + "> [DOI: 10.1103/PhysRevB.76.073103](https://doi.org/10.1103/PhysRevB.76.073103)" ] }, { @@ -38,42 +33,25 @@ "metadata": {}, "outputs": [], "source": [ - "FILM_MILLER_INDICES = (0, 0, 1)\n", - "FILM_THICKNESS = 1 # in atomic layers\n", - "FILM_TERMINATION_FORMULA = None # if None, the first termination will be used\n", - "FILM_VACUUM = 0.0 # in angstroms\n", - "\n", - "SUBSTRATE_MILLER_INDICES = (0, 0, 1)\n", - "SUBSTRATE_THICKNESS = 1 # in atomic layers\n", - "SUBSTRATE_TERMINATION_FORMULA = None # if None, the first termination will be used\n", - "SUBSTRATE_VACUUM = 0.0 # in angstroms\n", - "\n", - "INTERFACE_DISTANCE = 3.4 # Gap between substrate and film, in Angstrom\n", - "INTERFACE_VACUUM = 20.0 # Vacuum over film, in Angstrom\n", - "\n", - "# Whether to convert materials to conventional cells before creating slabs.\n", - "USE_CONVENTIONAL_CELL = True\n", - "\n", - "# Maximum area for the superlattice search algorithm (the final interface area will be smaller)\n", - "MAX_AREA = 350 # in Angstrom^2\n", - "# Additional fine-tuning parameters (increase values to get more strained matches):\n", - "MAX_AREA_TOLERANCE = 0.09 # in Angstrom^2\n", - "MAX_LENGTH_TOLERANCE = 0.05\n", - "MAX_ANGLE_TOLERANCE = 0.02\n", - "\n", - "# Whether to reduce the resulting interface cell to the primitive cell after the interface creation.\n", - "REDUCE_RESULT_CELL_TO_PRIMITIVE = True\n", - "\n", - "# One name per shift (n = 2..8); the three symmetric registries carry their Jung 2015 / Giovannetti 2007 label,\n", - "# n = 8 repeats n = 2 one period later.\n", - "INTERFACE_NAMES = [\n", - " \"Gr/hBN d3.4 shift 0of6 BA\",\n", - " \"Gr/hBN d3.4 shift 1of6\",\n", - " \"Gr/hBN d3.4 shift 2of6 AA\",\n", - " \"Gr/hBN d3.4 shift 3of6\",\n", - " \"Gr/hBN d3.4 shift 4of6 AB\",\n", - " \"Gr/hBN d3.4 shift 5of6\",\n", - " \"Gr/hBN d3.4 shift 6of6 BA\",\n", + "LATTICE_CONSTANT = 2.445 # Å, graphene LDA (Giovannetti et al. 2007); h-BN is compressed to it\n", + "H_BN_LAYERS = 4\n", + "H_BN_INTERLAYER_DISTANCE = 3.24 # Å, the paper's LDA value\n", + "VACUUM = 15.0 # Å above graphene\n", + "\n", + "# Registry of graphene on the top h-BN layer, Giovannetti et al. 2007 Fig. 1: one C over B and the other over\n", + "# N (a), over N and a hexagon centre (b), over B and a hexagon centre (c). All three run in the paper;\n", + "# uncomment to build the others.\n", + "STACKINGS = [\n", + " \"c\",\n", + " # \"a\",\n", + " # \"b\",\n", + "]\n", + "STACKING_SHIFTS = {\"a\": 0, \"b\": 1, \"c\": -1} # in units of a/√3 along y\n", + "# Graphene–h-BN distance, Å. The paper scans 2.5–3.9; uncomment to build the full set.\n", + "DISTANCES = [\n", + " 3.2,\n", + " 3.4,\n", + " # 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.3, 3.5, 3.6, 3.7, 3.8, 3.9,\n", "]" ] }, @@ -101,7 +79,7 @@ "metadata": {}, "source": [ "### 1.3. Get input materials and assign `substrate` and `film`\n", - "Materials are loaded with `get_data()`. The first material is assigned as substrate and the second as film." + "Bulk h-BN (substrate) and graphene (film) are loaded from Standata." ] }, { @@ -113,8 +91,8 @@ "from mat3ra.standata.materials import Materials\n", "from mat3ra.made.material import Material\n", "\n", - "film = Material.create(Materials.get_by_name_first_match(\"Graphene\"))\n", - "substrate = Material.create(Materials.get_by_name_and_categories(\"BN\", \"2D\"))" + "film = Material.create(Materials.get_by_name_first_match(\"C, Graphene, HEX (P6/mmm) 2D (Monolayer), 2dm-3993\"))\n", + "substrate = Material.create(Materials.get_by_name_first_match(\"BN, Boron Nitride, HEX (P6_3/mmc) 3D (Bulk), mp-7991\"))" ] }, { @@ -139,30 +117,10 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## 2. Configure slabs for interface\n", - "\n", - "### 2.1. Get possible terminations for the slabs" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from mat3ra.made.tools.helpers import get_slab_terminations\n", + "## 2. Prepare the slabs\n", "\n", - "film_slab_terminations = get_slab_terminations(material=film, miller_indices=FILM_MILLER_INDICES)\n", - "substrate_slab_terminations = get_slab_terminations(material=substrate, miller_indices=SUBSTRATE_MILLER_INDICES)\n", - "print(\"Film slab terminations:\", film_slab_terminations)\n", - "print(\"Substrate slab terminations:\", substrate_slab_terminations)\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### 2.2. Visualize slabs for all possible terminations" + "### 2.1. Strain the materials to the common lattice constant\n", + "Both materials are strained in-plane to `LATTICE_CONSTANT`, and h-BN along c to `H_BN_INTERLAYER_DISTANCE` between layers." ] }, { @@ -171,30 +129,23 @@ "metadata": {}, "outputs": [], "source": [ - "from mat3ra.made.tools.helpers import create_slab\n", + "from mat3ra.made.tools.build_components.operations.core.modifications.strain.helpers import create_strain\n", "\n", - "film_slabs = [create_slab(film, miller_indices=FILM_MILLER_INDICES, termination_top=termination, vacuum=0) for termination\n", - " in\n", - " film_slab_terminations]\n", - "\n", - "substrate_slabs = [create_slab(substrate, miller_indices=SUBSTRATE_MILLER_INDICES, termination_top=termination, vacuum=0)\n", - " for termination in\n", - " substrate_slab_terminations]\n", - "\n", - "film_slabs_with_titles = [{\"material\": slab, \"title\": str(termination)} for slab, termination in\n", - " zip(film_slabs, film_slab_terminations)]\n", - "substrate_slabs_with_titles = [{\"material\": slab, \"title\": str(termination)} for slab, termination in\n", - " zip(substrate_slabs, substrate_slab_terminations)]\n", - "\n", - "visualize(film_slabs_with_titles, repetitions=[3, 3, 1], rotation=\"-90x\")\n", - "visualize(substrate_slabs_with_titles, repetitions=[3, 3, 1], rotation=\"-90x\")" + "film_scale = LATTICE_CONSTANT / film.lattice.a\n", + "substrate_scale = LATTICE_CONSTANT / substrate.lattice.a\n", + "substrate_c_scale = 2 * H_BN_INTERLAYER_DISTANCE / substrate.lattice.c\n", + "film = create_strain(film, strain_matrix=[[film_scale, 0, 0], [0, film_scale, 0], [0, 0, 1]])\n", + "substrate = create_strain(\n", + " substrate, strain_matrix=[[substrate_scale, 0, 0], [0, substrate_scale, 0], [0, 0, substrate_c_scale]]\n", + ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### 2.3. Select terminations for the Slabs" + "### 2.2. Set the AA' stacking of h-BN\n", + "The atoms of the upper layer of the bulk h-BN cell are moved by (1/3, 2/3, 0) in crystal coordinates, so that B sits over N in adjacent layers." ] }, { @@ -203,19 +154,21 @@ "metadata": {}, "outputs": [], "source": [ - "from mat3ra.made.tools.helpers import select_slab_termination\n", + "from mat3ra.made.tools.analyze.other import get_atom_indices_with_condition_on_coordinates\n", + "from mat3ra.made.tools.operations.core.unary import translate_atoms\n", "\n", - "film_termination = select_slab_termination(film_slab_terminations, FILM_TERMINATION_FORMULA)\n", - "substrate_termination = select_slab_termination(substrate_slab_terminations, SUBSTRATE_TERMINATION_FORMULA)" + "upper_layer_ids = get_atom_indices_with_condition_on_coordinates(substrate, lambda coordinate: coordinate[2] > 0.5)\n", + "substrate = translate_atoms(\n", + " substrate, atom_ids=upper_layer_ids, vector=[1 / 3, 2 / 3, 0], use_cartesian_coordinates=False\n", + ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### 2.4. Create Substrate and Film Slabs\n", - "Slab Configuration lets define the slab thickness, vacuum, and the Miller indices of the interfacial plane and get the slabs with possible terminations.\n", - "Define the substrate slab cell that will be used as a base for the interface and the film slab cell that will be placed on top of the substrate slab.\n" + "### 2.3. Create Substrate and Film Slabs\n", + "The substrate slab has `H_BN_LAYERS` layers of h-BN, the film slab one layer of graphene." ] }, { @@ -228,20 +181,16 @@ "\n", "substrate_slab_config = SlabConfiguration.from_parameters(\n", " material_or_dict=substrate,\n", - " miller_indices=SUBSTRATE_MILLER_INDICES,\n", - " number_of_layers=SUBSTRATE_THICKNESS,\n", + " miller_indices=(0, 0, 1),\n", + " number_of_layers=H_BN_LAYERS // 2, # two BN layers per bulk cell\n", " vacuum=0.0,\n", - " termination_top_formula=SUBSTRATE_TERMINATION_FORMULA,\n", - " use_conventional_cell=USE_CONVENTIONAL_CELL\n", ")\n", "\n", "film_slab_config = SlabConfiguration.from_parameters(\n", " material_or_dict=film,\n", - " miller_indices=FILM_MILLER_INDICES,\n", - " number_of_layers=FILM_THICKNESS,\n", + " miller_indices=(0, 0, 1),\n", + " number_of_layers=1,\n", " vacuum=0.0,\n", - " termination_top_formula=FILM_TERMINATION_FORMULA,\n", - " use_conventional_cell=USE_CONVENTIONAL_CELL\n", ")\n", "\n", "substrate_slab = SlabBuilder().get_material(substrate_slab_config)\n", @@ -252,150 +201,10 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## 3. Analyze possible interfaces with ZSL Analyzer\n", - "### 3.1. Initialize ZSL Analyzer\n", - "The search algorithm for supercells matching can be tuned by setting its parameters directly, otherwise the default values are used." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from mat3ra.made.tools.analyze.interface import ZSLInterfaceAnalyzer\n", + "## 3. Create the interfaces\n", "\n", - "zsl_analyzer = ZSLInterfaceAnalyzer(\n", - " substrate_slab_configuration=substrate_slab_config,\n", - " film_slab_configuration=film_slab_config,\n", - " max_area=MAX_AREA,\n", - " max_area_ratio_tol=MAX_AREA_TOLERANCE,\n", - " max_length_tol=MAX_LENGTH_TOLERANCE,\n", - " max_angle_tol=MAX_ANGLE_TOLERANCE,\n", - " reduce_result_cell=False # Reduces supercell matrices in analyzer\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### 3.2. Generate matches with strain analyzer\n", - "Matches are sorted by size and strain." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "matches = zsl_analyzer.zsl_match_holders" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### 3.3. Plot matches by area and strain" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from mat3ra.notebooks_utils.ipython.entity.material.plot import plot_strain_vs_area\n", - "\n", - "PLOT_SETTINGS = {\n", - " \"HEIGHT\": 600,\n", - " \"X_SCALE\": \"log\", # or linear\n", - " \"Y_SCALE\": \"log\", # or linear\n", - "}\n", - "\n", - "plot_strain_vs_area(matches, PLOT_SETTINGS)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### 3.4. Select the interface\n", - "\n", - "Select the index for the interface with the lowest strain and the smallest area." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "selected_index = 0\n", - "\n", - "from mat3ra.made.tools.helpers import create_interface_zsl_between_slabs as create_zsl_interface_between_slabs\n", - "\n", - "interface = create_zsl_interface_between_slabs(\n", - " substrate_slab=substrate_slab,\n", - " film_slab=film_slab,\n", - " gap=INTERFACE_DISTANCE,\n", - " vacuum=INTERFACE_VACUUM,\n", - " match_id=selected_index,\n", - " max_area=MAX_AREA,\n", - " max_area_ratio_tol=MAX_AREA_TOLERANCE,\n", - " max_length_tol=MAX_LENGTH_TOLERANCE,\n", - " max_angle_tol=MAX_ANGLE_TOLERANCE,\n", - " reduce_result_cell_to_primitive=REDUCE_RESULT_CELL_TO_PRIMITIVE,\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### 3.5. Set the cell to the standard hexagonal setting\n", - "The ZSL interface cell comes out with γ = 60°. The symbolic K point the band-structure notebook uses assumes the standard 120° hexagonal cell, so the cell is re-set here (same four atoms) and each shifted interface is typed HEX." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from mat3ra.made.tools.helpers import create_supercell\n", - "\n", - "interface = create_supercell(interface, supercell_matrix=[[1, 0, 0], [-1, 1, 0], [0, 0, 1]])\n", - "print(f\"{len(interface.basis.elements.ids)} atoms, a = {interface.lattice.a:.4f} Å, \"\n", - " f\"gamma = {interface.lattice.gamma:.1f}°\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 4. Preview the selected interface and create variants\n", - "\n", - "### 4.1. Preview the selected interface" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "visualize(interface, repetitions=[3, 3, 1])\n", - "visualize(interface, repetitions=[3, 3, 1], rotation=\"-90x\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### 4.2. Shift film along y-axis\n", - "Shifting with the step of a/sqrt(3)/2 angstroms, where a is the lattice constant of the interface." + "### 3.1. Place graphene at each distance and shift it into each registry\n", + "The film is placed on the substrate slab at each of `DISTANCES` and shifted along y by `STACKING_SHIFTS` for each of `STACKINGS`; the registry of each carbon atom is measured on the result." ] }, { @@ -407,29 +216,51 @@ "import numpy as np\n", "from mat3ra.made.tools.analyze.other import get_average_interlayer_distance\n", "from mat3ra.made.tools.convert.interface_parts_enum import InterfacePartsEnum\n", + "from mat3ra.made.tools.helpers import create_interface_zsl_between_slabs as create_zsl_interface_between_slabs\n", "from mat3ra.made.tools.modify import interface_displace_part\n", "\n", - "a = interface.lattice.a\n", - "shifted_interfaces = []\n", - "for index, n in enumerate(range(2, 9)):\n", - " shifted_interface = interface_displace_part(\n", - " interface=interface,\n", - " displacement=[0, n * a / np.sqrt(3) / 2, 0],\n", - " use_cartesian_coordinates=True)\n", - " shifted_interface.name = INTERFACE_NAMES[index]\n", - " shifted_interface.lattice.type = \"HEX\"\n", - " shifted_interfaces.append(shifted_interface)\n", - " interlayer_distance = get_average_interlayer_distance(\n", - " shifted_interface, InterfacePartsEnum.SUBSTRATE.value, InterfacePartsEnum.FILM.value)\n", - " print(f\"{shifted_interface.name}: {len(shifted_interface.basis.elements.ids)} atoms, \"\n", - " f\"gamma = {shifted_interface.lattice.gamma:.1f}°, interlayer distance = {interlayer_distance:.3f} Å\")" + "\n", + "def get_registry(interface):\n", + " elements = np.array(interface.basis.elements.values)\n", + " coordinates = np.array(interface.basis.coordinates.values)\n", + " top_layer = (elements != \"C\") & np.isclose(coordinates[:, 2], coordinates[elements != \"C\", 2].max())\n", + " registry = []\n", + " for carbon in coordinates[elements == \"C\"]:\n", + " in_plane_offsets = (coordinates[top_layer, :2] - carbon[:2] + 0.5) % 1 - 0.5\n", + " atoms_below = elements[top_layer][np.all(np.abs(in_plane_offsets) < 1e-3, axis=1)]\n", + " registry.append(atoms_below[0] if len(atoms_below) else \"hollow\")\n", + " return registry\n", + "\n", + "\n", + "interfaces = []\n", + "for stacking in STACKINGS:\n", + " for distance in DISTANCES:\n", + " # the builder adds the gap to the vacuum above the film\n", + " interface = create_zsl_interface_between_slabs(\n", + " substrate_slab=substrate_slab, film_slab=film_slab, gap=distance, vacuum=VACUUM - distance\n", + " )\n", + " interface = interface_displace_part(\n", + " interface=interface,\n", + " displacement=[0, STACKING_SHIFTS[stacking] * interface.lattice.a / np.sqrt(3), 0],\n", + " use_cartesian_coordinates=True,\n", + " )\n", + " interface.name = f\"Gr/hBN ({stacking}) d{distance:.2f}\"\n", + " interface.lattice.type = \"HEX\"\n", + " interfaces.append(interface)\n", + " interlayer_distance = get_average_interlayer_distance(\n", + " interface, InterfacePartsEnum.SUBSTRATE.value, InterfacePartsEnum.FILM.value\n", + " )\n", + " vacuum = interface.lattice.c * (1 - max(coordinate[2] for coordinate in interface.basis.coordinates.values))\n", + " print(f\"{interface.name}: {len(interface.basis.elements.ids)} atoms, a = {interface.lattice.a:.4f} Å, \"\n", + " f\"gamma = {interface.lattice.gamma:.1f}°, distance = {interlayer_distance:.3f} Å, \"\n", + " f\"vacuum = {vacuum:.2f} Å, C over {' / '.join(get_registry(interface))}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### 4.3. Preview the shifted materials" + "### 3.2. Preview the interfaces" ] }, { @@ -438,14 +269,15 @@ "metadata": {}, "outputs": [], "source": [ - "visualize(shifted_interfaces, repetitions=[3, 3, 1])" + "visualize(interfaces, repetitions=[3, 3, 1])\n", + "visualize(interfaces, repetitions=[3, 3, 1], rotation=\"-90x\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## 5. Pass data to the outside runtime" + "## 4. Pass data to the outside runtime" ] }, { @@ -457,9 +289,9 @@ "from mat3ra.notebooks_utils.io import download_content_to_file\n", "from mat3ra.notebooks_utils.material import set_materials\n", "\n", - "set_materials(shifted_interfaces)\n", - "for shifted_interface in shifted_interfaces:\n", - " download_content_to_file(shifted_interface.to_json(), f\"{shifted_interface.name}.json\")" + "set_materials(interfaces)\n", + "for interface in interfaces:\n", + " download_content_to_file(interface.to_json(), f\"{interface.name}.json\")" ] } ], From 102092b41197ac34338980dfc8449dedb36d7415 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Wed, 7 Oct 2026 19:43:17 -0700 Subject: [PATCH 11/17] SOF-8064: three default distances around the (c) minimum Co-Authored-By: Claude Opus 5.5 --- .../interface_2d_2d_boron_nitride_graphene.ipynb | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb index ac49af613..8348bf116 100644 --- a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb @@ -47,11 +47,11 @@ " # \"b\",\n", "]\n", "STACKING_SHIFTS = {\"a\": 0, \"b\": 1, \"c\": -1} # in units of a/√3 along y\n", - "# Graphene–h-BN distance, Å. The paper scans 2.5–3.9; uncomment to build the full set.\n", + "# Graphene–h-BN distance, Å. The paper scans 2.5–3.9; three points around (c)'s minimum are active, uncomment\n", + "# the rest for the full set.\n", "DISTANCES = [\n", - " 3.2,\n", - " 3.4,\n", - " # 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.3, 3.5, 3.6, 3.7, 3.8, 3.9,\n", + " 3.1, 3.2, 3.3,\n", + " # 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9,\n", "]" ] }, From 3ab4acac8ef206c529958f2e75f08feb5004c24e Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Wed, 7 Oct 2026 19:49:40 -0700 Subject: [PATCH 12/17] =?UTF-8?q?SOF-8064:=20simulation=20notebook=20?= =?UTF-8?q?=E2=80=94=20the=20paper's=20stackings=20=C3=97=20distances,=20t?= =?UTF-8?q?etrahedron=20occupations,=20dipole=20correction?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Opus 5.5 --- ...2d_boron_nitride_graphene_SIMULATION.ipynb | 363 ++++++++++++------ 1 file changed, 256 insertions(+), 107 deletions(-) diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb index 392a370b3..141ed0214 100644 --- a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb @@ -5,37 +5,37 @@ "id": "0", "metadata": {}, "source": [ - "# Stacking Energy and Band Gap of Graphene on h-BN\n", + "# Substrate-Induced Band Gap of Graphene on h-BN\n", "\n", - "> **Gianluca Giovannetti, Petr A. Khomyakov, Geert Brocks, Paul J. Kelly & Jeroen van den Brink**\n", - "> Substrate-induced band gap in graphene on hexagonal boron nitride: Ab initio density functional calculations. Physical Review B, 76, 073103. 2007.\n", - "> [https://doi.org/10.1103/PhysRevB.76.073103](https://doi.org/10.1103/PhysRevB.76.073103)\n", + "> **Gianluca Giovannetti, Petr A. Khomyakov, Geert Brocks, Paul J. Kelly, and Jeroen van den Brink**\n", + "> **Substrate-induced band gap in graphene on hexagonal boron nitride: Ab initio density functional calculations**\n", + "> Physical Review B 76, 073103 (2007)\n", + "> [DOI: 10.1103/PhysRevB.76.073103](https://doi.org/10.1103/PhysRevB.76.073103)\n", "\n", - "Calculate the total energy and the band structure of the seven graphene/h-BN stackings at d = 3.4 Å\n", - "created in the [structure notebook](interface_2d_2d_boron_nitride_graphene.ipynb), using Quantum\n", - "ESPRESSO and the band structure + density of states workflow from Standata, and compare the\n", - "ordering of the sliding energies and the direct gap at K with Giovannetti et al. (2007).\n", - "\n", - "The seven stackings follow the sliding path of Fig. 7(a) in Jung et al. (2015). That figure is a model\n", - "curve parameterised from RPA calculations; the DFT numbers compared here are Giovannetti's.\n", + "Calculate the total energy and the band structure of graphene on four layers of h-BN for the three\n", + "stackings and the list of graphene–h-BN distances created in the\n", + "[structure notebook](interface_2d_2d_boron_nitride_graphene.ipynb), using Quantum ESPRESSO and the\n", + "band structure + density of states workflow from Standata. The results reproduce Giovannetti et al.\n", + "(2007): Fig. 2 (total energy vs distance and the equilibrium distances), Fig. 3 (bands and density of\n", + "states of stacking (c) at its equilibrium distance, gap at K) and Fig. 4 (gap at K vs distance).\n", "\n", "

Usage

\n", "\n", - "1. Create the materials in the [structure notebook](interface_2d_2d_boron_nitride_graphene.ipynb), which saves them to the `uploads` folder under the names used in cell 1.2 below.\n", - "1. Set the materials and the calculation parameters in cells 1.2 and 1.3, or keep the default values: everything runs with them.\n", + "1. Create the materials in the [structure notebook](interface_2d_2d_boron_nitride_graphene.ipynb), which saves them to the `uploads` folder under the names built from `MATERIAL_NAME` in cell 1.2 below.\n", + "1. Set the stackings, distances and calculation parameters in cells 1.2 and 1.3, or keep the default values: the active stackings × distances run with them. Uncomment the rest of `STACKINGS` and `DISTANCES`, here and in the structure notebook, for the paper's full set.\n", "1. Click \"Run\" > \"Run All\" to run all cells.\n", "1. Wait for the jobs to complete.\n", - "1. Scroll down to view the results; section 8 prints the comparison table and plots.\n", + "1. Scroll down to view the results; section 8 prints the comparison with the paper.\n", "\n", "## Summary\n", "\n", - "1. Set up the environment and parameters: install packages (JupyterLite only) and configure the materials, workflow, model, compute resources and jobs.\n", + "1. Set up the environment and parameters: install packages (JupyterLite only) and configure the stackings, distances, workflow, model, compute resources and jobs.\n", "1. Authenticate and initialize API client: authenticate via browser, initialize the client, then select account and project.\n", "1. Load the materials by name from the `uploads` folder, print their provenance and save them to the platform.\n", - "1. Configure the workflow: select the application and the model, load the band structure + DOS workflow from Standata, and set the computational parameters for each material.\n", + "1. Configure the workflow: select the application and the model, load the band structure + DOS workflow from Standata, and set the computational parameters, the tetrahedron occupations and the dipole correction for each material.\n", "1. Configure compute: get the list of clusters and create a compute configuration.\n", "1. Run one job per material, one after another, re-using a job that already ran under the same workflow name.\n", - "1. Retrieve results: band structures, total energies and the direct gap along the path.\n", + "1. Retrieve results: total energies, the direct gap at K, the equilibrium distance of each stacking, the total energy and the gap vs distance, the bands and density of states of (c) at its equilibrium distance, the h-BN gap and the effective mass at K.\n", "1. Compare with Giovannetti et al. (2007)." ] }, @@ -81,16 +81,16 @@ "ORGANIZATION_NAME = None # set to your organization name (full or partial); otherwise, your default one is used\n", "FOLDER = \"./uploads\"\n", "\n", - "# Names saved by the structure notebook; the symmetric stackings carry their Jung 2015 / Giovannetti 2007 label\n", - "MATERIALS = {\n", - " \"Gr/hBN d3.4 shift 0of6 BA\": {\"shift\": 0, \"stacking\": \"BA\"},\n", - " \"Gr/hBN d3.4 shift 1of6\": {\"shift\": 1, \"stacking\": None},\n", - " \"Gr/hBN d3.4 shift 2of6 AA\": {\"shift\": 2, \"stacking\": \"AA\"},\n", - " \"Gr/hBN d3.4 shift 3of6\": {\"shift\": 3, \"stacking\": None},\n", - " \"Gr/hBN d3.4 shift 4of6 AB\": {\"shift\": 4, \"stacking\": \"AB\"},\n", - " \"Gr/hBN d3.4 shift 5of6\": {\"shift\": 5, \"stacking\": None},\n", - " \"Gr/hBN d3.4 shift 6of6 BA\": {\"shift\": 6, \"stacking\": \"BA\"},\n", - "}\n", + "STACKINGS = [\n", + " \"c\",\n", + " # \"a\",\n", + " # \"b\",\n", + "]\n", + "DISTANCES = [\n", + " 3.1, 3.2, 3.3,\n", + " # 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9,\n", + "]\n", + "MATERIAL_NAME = \"Gr/hBN ({stacking}) d{distance:.2f}\" # as saved by the structure notebook\n", "\n", "WORKFLOW_SEARCH_TERM = \"band_structure_dos.json\"\n", "MY_WORKFLOW_NAME = \"Band Structure + DOS\"\n", @@ -99,7 +99,7 @@ "CLUSTER_NAME = \"001\" # specify full or partial name i.e. \"cluster-001\" to select\n", "QUEUE_NAME = QueueName.D\n", "PPN = 1\n", - "TIME_LIMIT = \"01:00:00\"\n", + "TIME_LIMIT = \"04:00:00\"\n", "\n", "timestamp = datetime.now().strftime(\"%Y-%m-%d %H:%M\")\n", "POLL_INTERVAL = 60 # seconds" @@ -123,18 +123,20 @@ "MODEL_SUBTYPE = \"lda\"\n", "FUNCTIONAL = \"pz\" # Giovannetti et al. 2007 use LDA: GGA gives essentially no interlayer binding\n", "PSEUDOPOTENTIAL_TYPE = \"us\" # GBRV ultrasoft, the only LDA family the platform publishes for B, C and N\n", - "ECUTWFC = 40 # Ry, GBRV's tested cutoff\n", + "ECUTWFC = 40 # Ry, GBRV's tested cutoff; the paper's 600 eV is a VASP number\n", "ECUTRHO = 200 # Ry, GBRV's tested charge-density cutoff\n", "\n", "KGRID = [36, 36, 1] # Giovannetti et al. 2007; a multiple of 3 keeps K on the mesh\n", - "SMEARING_SETTINGS = {\"degauss\": 0.001} # Ry; the gaps compared are 30-80 meV\n", + "OCCUPATIONS_SETTINGS = {\"occupations\": \"tetrahedra\"} # the paper's tetrahedron method\n", + "# Dipole correction along z; emaxpos is set per material in 4.2, in the middle of the vacuum above graphene\n", + "DIPOLE_SETTINGS = {\"control\": {\"tefield\": True, \"dipfield\": True}, \"system\": {\"edir\": 3, \"eopreg\": 0.05}}\n", "KPATH_STEPS = 40\n", - "MODEL_TAG = f\"{FUNCTIONAL}-{PSEUDOPOTENTIAL_TYPE} {ECUTWFC}-{ECUTRHO}Ry k{KGRID[0]} p{KPATH_STEPS} g{SMEARING_SETTINGS['degauss']}\"\n", + "MODEL_TAG = f\"{FUNCTIONAL}-{PSEUDOPOTENTIAL_TYPE} {ECUTWFC}-{ECUTRHO}Ry k{KGRID[0]} p{KPATH_STEPS} tetra dip\"\n", "\n", "SCF_UNIT = \"pw_scf\"\n", "NSCF_UNIT = \"pw_nscf\"\n", "BANDS_UNIT = \"pw_bands\"\n", - "NUMBER_OF_OCCUPIED_BANDS = 8 # 16 valence electrons: C 4 + 4, B 3, N 5\n", + "NUMBER_OF_OCCUPIED_BANDS = 20 # 40 valence electrons: C 4×2, B 3×4, N 5×4\n", "\n", "KPATH = [\n", " {\"point\": \"Γ\", \"steps\": KPATH_STEPS},\n", @@ -250,7 +252,7 @@ "### 3.1. Load from the uploads folder and print provenance\n", "\n", "The structures, their geometry and their lattice type all come from the structure notebook; this one\n", - "only loads them by name." + "only loads them by name, one per stacking and distance." ] }, { @@ -260,22 +262,24 @@ "metadata": {}, "outputs": [], "source": [ - "from collections import Counter\n", "from mat3ra.made.tools.analyze.other import get_average_interlayer_distance\n", "from mat3ra.made.tools.convert.interface_parts_enum import InterfacePartsEnum\n", "from mat3ra.notebooks_utils.core.entity.material.api import load_material\n", "\n", + "material_names = {\n", + " stacking: {distance: MATERIAL_NAME.format(stacking=stacking, distance=distance) for distance in sorted(DISTANCES)}\n", + " for stacking in STACKINGS\n", + "}\n", "materials = {}\n", - "for name, settings in MATERIALS.items():\n", - " material = load_material(client, FOLDER, name, ACCOUNT_ID)\n", - " materials[name] = material\n", - " composition = \"\".join(\n", - " f\"{element}{count}\" for element, count in sorted(Counter(material.basis.elements.values).items()))\n", - " interlayer_distance = get_average_interlayer_distance(\n", - " material, InterfacePartsEnum.SUBSTRATE.value, InterfacePartsEnum.FILM.value)\n", - " print(f\"{name}: {composition}, {material.basis.number_of_atoms} atoms, a = {material.lattice.a:.4f} Å, \"\n", - " f\"gamma = {material.lattice.gamma:.3f}°, interlayer distance = {interlayer_distance:.3f} Å, \"\n", - " f\"{settings['stacking'] or 'bridge'}\")" + "for stacking, names in material_names.items():\n", + " for name in names.values():\n", + " material = load_material(client, FOLDER, name, ACCOUNT_ID)\n", + " materials[name] = material\n", + " interlayer_distance = get_average_interlayer_distance(\n", + " material, InterfacePartsEnum.SUBSTRATE.value, InterfacePartsEnum.FILM.value)\n", + " print(f\"{name}: {material.basis.number_of_atoms} atoms, a = {material.lattice.a:.4f} Å, \"\n", + " f\"gamma = {material.lattice.gamma:.3f}°, graphene–h-BN distance = {interlayer_distance:.3f} Å, \"\n", + " f\"stacking ({stacking})\")" ] }, { @@ -387,8 +391,14 @@ " material=next(iter(materials.values())), dimensions=KGRID, isEdited=True).get_context_item_data()\n", "path_context = PointsPathDataProvider(path=KPATH, isEdited=True).get_context_item_data()\n", "\n", + "\n", + "def get_vacuum_center(material):\n", + " heights = [coordinate[2] for coordinate in material.basis.coordinates.values]\n", + " return round((max(heights) + 1 + min(heights)) / 2, 4)\n", + "\n", + "\n", "workflows = {}\n", - "for name in MATERIALS:\n", + "for name, material in materials.items():\n", " workflow = Workflow.create(deepcopy(workflow_config))\n", " workflow.name = f\"{MY_WORKFLOW_NAME} {name} {MODEL_TAG}\"\n", " subworkflow = workflow.subworkflows[0]\n", @@ -401,7 +411,9 @@ " for context in contexts:\n", " unit.add_context(context)\n", " subworkflow.set_unit(unit)\n", - " patch_workflow_qe_input(workflow, {\"system\": SMEARING_SETTINGS}, [SCF_UNIT, NSCF_UNIT, BANDS_UNIT])\n", + " patch_workflow_qe_input(workflow, DIPOLE_SETTINGS, [SCF_UNIT, NSCF_UNIT, BANDS_UNIT])\n", + " patch_workflow_qe_input(workflow, {\"system\": {**OCCUPATIONS_SETTINGS, \"emaxpos\": get_vacuum_center(material)}},\n", + " [SCF_UNIT, NSCF_UNIT, BANDS_UNIT])\n", " workflows[name] = workflow\n", " print(workflow.name)" ] @@ -423,7 +435,7 @@ "source": [ "from mat3ra.notebooks_utils.ipython.entity.workflow.visualize import visualize_workflow\n", "\n", - "visualize_workflow(workflows[next(iter(MATERIALS))])" + "visualize_workflow(next(iter(workflows.values())))" ] }, { @@ -543,7 +555,7 @@ "metadata": {}, "source": [ "## 7. Retrieve the results\n", - "### 7.1. Band structures" + "### 7.1. Total energies" ] }, { @@ -554,13 +566,11 @@ "outputs": [], "source": [ "from mat3ra.notebooks_utils.core.entity.property.api import get_properties_for_job\n", - "from mat3ra.notebooks_utils.ipython.entity.property.visualize import visualize_properties\n", "\n", - "band_structures = {}\n", + "total_energies = {}\n", "for name, job_id in job_ids.items():\n", - " band_structures[name] = get_properties_for_job(client, job_id, property_name=\"band_structure\")\n", - " visualize_properties(band_structures[name], title=f\"Band Structure: {name}\",\n", - " extra_config={\"material\": materials[name].to_dict()})" + " total_energies[name] = get_properties_for_job(client, job_id, \"total_energy\")[0][\"value\"]\n", + " print(f\"{name}: {total_energies[name]:.6f} eV\")" ] }, { @@ -568,7 +578,10 @@ "id": "38", "metadata": {}, "source": [ - "### 7.2. Total energies" + "### 7.2. Direct gap at K\n", + "\n", + "The difference between the lowest unoccupied and the highest occupied band at K, the end of the first\n", + "leg of `KPATH`." ] }, { @@ -578,10 +591,17 @@ "metadata": {}, "outputs": [], "source": [ - "total_energies = {}\n", + "import numpy as np\n", + "\n", + "K_VERTEX_INDEX = KPATH_STEPS\n", + "\n", + "band_structures = {}\n", + "gaps = {}\n", "for name, job_id in job_ids.items():\n", - " total_energies[name] = get_properties_for_job(client, job_id, \"total_energy\")[0][\"value\"]\n", - " print(f\"{name}: {total_energies[name]:.6f} eV\")" + " band_structures[name] = get_properties_for_job(client, job_id, property_name=\"band_structure\")[0]\n", + " bands_at_k = np.array(band_structures[name][\"yDataSeries\"])[:, K_VERTEX_INDEX]\n", + " gaps[name] = 1000 * (bands_at_k[NUMBER_OF_OCCUPIED_BANDS] - bands_at_k[NUMBER_OF_OCCUPIED_BANDS - 1])\n", + " print(f\"{name}: {gaps[name]:.1f} meV\")" ] }, { @@ -589,11 +609,10 @@ "id": "40", "metadata": {}, "source": [ - "### 7.3. Direct gap along the path\n", + "### 7.3. Equilibrium distance of each stacking\n", "\n", - "The smallest difference along Γ–K–M–Γ between the lowest unoccupied and the highest occupied band;\n", - "the k-point where it occurs is printed in crystal coordinates. For the three symmetric stackings it\n", - "sits at K = (1/3, 1/3, 0)." + "The vertex of the parabola through the three distances with the lowest total energy; the gap at K is\n", + "printed for the distance in `DISTANCES` nearest to it." ] }, { @@ -603,81 +622,211 @@ "metadata": {}, "outputs": [], "source": [ - "import numpy as np\n", - "\n", - "gaps = {}\n", - "for name, band_structure in band_structures.items():\n", - " bands = np.array(band_structure[0][\"yDataSeries\"])\n", - " gaps_along_path = bands[NUMBER_OF_OCCUPIED_BANDS] - bands[NUMBER_OF_OCCUPIED_BANDS - 1]\n", - " kpoint_index = int(np.argmin(gaps_along_path))\n", - " gaps[name] = 1000 * gaps_along_path[kpoint_index]\n", - " kpoint = band_structure[0][\"xDataArray\"][kpoint_index]\n", - " print(f\"{name}: {gaps[name]:.1f} meV at k = ({kpoint[0]:.4f}, {kpoint[1]:.4f}, {kpoint[2]:.4f})\")" + "equilibrium_distances = {}\n", + "equilibrium_names = {}\n", + "for stacking, names in material_names.items():\n", + " lowest = sorted(names, key=lambda distance: total_energies[names[distance]])[:3]\n", + " quadratic, linear, _ = np.polyfit(lowest, [total_energies[names[distance]] for distance in lowest], 2)\n", + " equilibrium_distances[stacking] = -linear / (2 * quadratic)\n", + " nearest = min(names, key=lambda distance: abs(distance - equilibrium_distances[stacking]))\n", + " equilibrium_names[stacking] = names[nearest]\n", + " print(f\"({stacking}): equilibrium distance {equilibrium_distances[stacking]:.3f} Å, \"\n", + " f\"gap at K at {nearest:.2f} Å: {gaps[names[nearest]]:.1f} meV\")" ] }, { "cell_type": "markdown", "id": "42", "metadata": {}, + "source": [ + "### 7.4. Total energy vs distance (Fig. 2 of the paper)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.plot import plot_series\n", + "\n", + "for stacking, names in material_names.items():\n", + " reference_energy = total_energies[names[max(names)]]\n", + " series = [\n", + " {\"distance\": distance, \"energy\": total_energies[name] - reference_energy} for distance, name in names.items()\n", + " ]\n", + " plot_series(series=series, x_key=\"distance\", y_key=\"energy\", xlabel=\"Graphene–h-BN distance (Å)\",\n", + " ylabel=f\"E − E({max(names):.2f} Å) (eV per cell)\",\n", + " title=f\"Total energy vs distance, stacking ({stacking})\")" + ] + }, + { + "cell_type": "markdown", + "id": "44", + "metadata": {}, + "source": [ + "### 7.5. Gap at K vs distance (Fig. 4 of the paper)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "45", + "metadata": {}, + "outputs": [], + "source": [ + "for stacking, names in material_names.items():\n", + " series = [{\"distance\": distance, \"gap\": gaps[name]} for distance, name in names.items()]\n", + " plot_series(series=series, x_key=\"distance\", y_key=\"gap\", xlabel=\"Graphene–h-BN distance (Å)\",\n", + " ylabel=\"Direct gap at K (meV)\", title=f\"Gap at K vs distance, stacking ({stacking})\")" + ] + }, + { + "cell_type": "markdown", + "id": "46", + "metadata": {}, + "source": [ + "### 7.6. Bands and density of states of (c) at its equilibrium distance (Fig. 3 of the paper)\n", + "\n", + "The job at the distance nearest to the equilibrium of (c): band structure, density of states, and the\n", + "two bands around the gap near K." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "47", + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib import pyplot as plt\n", + "from mat3ra.notebooks_utils.ipython.entity.property.visualize import visualize_properties\n", + "from mat3ra.notebooks_utils.plot import display_matplotlib_figure\n", + "\n", + "name = equilibrium_names[\"c\"]\n", + "visualize_properties(band_structures[name], title=f\"Band Structure: {name}\",\n", + " extra_config={\"material\": materials[name].to_dict()})\n", + "visualize_properties(get_properties_for_job(client, job_ids[name], property_name=\"density_of_states\"),\n", + " title=f\"Density of States: {name}\", extra_config={\"material\": materials[name].to_dict()})\n", + "\n", + "bands = np.array(band_structures[name][\"yDataSeries\"])\n", + "path_indices = np.arange(K_VERTEX_INDEX - 8, K_VERTEX_INDEX + 9)\n", + "center = bands[NUMBER_OF_OCCUPIED_BANDS - 1 : NUMBER_OF_OCCUPIED_BANDS + 1, K_VERTEX_INDEX].mean()\n", + "figure, axes = plt.subplots(figsize=(6, 5))\n", + "for band in bands[NUMBER_OF_OCCUPIED_BANDS - 1 : NUMBER_OF_OCCUPIED_BANDS + 1]:\n", + " axes.plot(path_indices, band[path_indices], marker=\"o\")\n", + "axes.set_ylim(center - 0.3, center + 0.3)\n", + "axes.set_xlabel(f\"Path index (K at {K_VERTEX_INDEX})\")\n", + "axes.set_ylabel(\"Energy (eV)\")\n", + "axes.set_title(\"Bands around K\")\n", + "display_matplotlib_figure(figure)" + ] + }, + { + "cell_type": "markdown", + "id": "48", + "metadata": {}, + "source": [ + "### 7.7. h-BN gap at K" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "49", + "metadata": {}, + "outputs": [], + "source": [ + "# the h-BN bands adjacent to the graphene pair\n", + "h_bn_gap = bands[NUMBER_OF_OCCUPIED_BANDS + 1, K_VERTEX_INDEX] - bands[NUMBER_OF_OCCUPIED_BANDS - 2, K_VERTEX_INDEX]\n", + "print(f\"{name}: h-BN gap at K {h_bn_gap:.2f} eV\")" + ] + }, + { + "cell_type": "markdown", + "id": "50", + "metadata": {}, + "source": [ + "### 7.8. Effective mass at K\n", + "\n", + "The curvature of the lowest unoccupied band at K, from the parabola through the two path points on\n", + "either side of it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "51", + "metadata": {}, + "outputs": [], + "source": [ + "HBAR_SQUARED_OVER_ELECTRON_MASS = 7.6200 # eV·Å²\n", + "K_STEP = 4 * np.pi / (3 * materials[name].lattice.a) / KPATH_STEPS # Å⁻¹ per path step along Γ–K\n", + "\n", + "# the K–M steps are half as long as the Γ–K steps\n", + "k_offsets = K_STEP * np.array([-2, -1, 0, 0.5, 1])\n", + "quadratic, _, _ = np.polyfit(k_offsets, bands[NUMBER_OF_OCCUPIED_BANDS, K_VERTEX_INDEX - 2 : K_VERTEX_INDEX + 3], 2)\n", + "effective_mass = HBAR_SQUARED_OVER_ELECTRON_MASS / (2 * quadratic)\n", + "print(f\"{name}: effective mass at K {effective_mass:.2e} m_e\")" + ] + }, + { + "cell_type": "markdown", + "id": "52", + "metadata": {}, "source": [ "## 8. Compare with Giovannetti et al. (2007)\n", "\n", - "ΔE is the total energy relative to the first entry of `MATERIALS`. The paper's AA and BA gaps are read\n", - "off its Fig. 4 at d = 3.4 Å, and AB is the paper's value at its 3.40 Å equilibrium; the gaps at its\n", - "own equilibrium distances are printed for reference.\n", + "The paper's equilibrium distances (Fig. 2) and gaps at equilibrium (Fig. 4, p. 3) for each active\n", + "stacking, then the h-BN gap and the effective mass at K for (c), printed beside the values above.\n", "\n", - "Settings that differ from the paper's: cell a = 2.509 Å (h-BN unstrained, graphene +1.79 %) against 2.445 Å; one h-BN layer against four; GBRV ultrasoft pseudopotentials against VASP at 600 eV; Gaussian smearing against the tetrahedron method; no dipole correction." + "Settings that differ from the paper's: GBRV ultrasoft pseudopotentials at 40/200 Ry against VASP at 600 eV." ] }, { "cell_type": "code", "execution_count": null, - "id": "43", + "id": "53", "metadata": {}, "outputs": [], "source": [ - "from mat3ra.notebooks_utils.plot import plot_series\n", + "GIOVANNETTI_EQUILIBRIUM_DISTANCE = {\"a\": 3.50, \"b\": 3.40, \"c\": 3.22} # Å, Fig. 2\n", + "GIOVANNETTI_GAP_AT_EQUILIBRIUM = {\"a\": 56, \"b\": 46, \"c\": 53} # meV, Fig. 4 / p.3\n", + "GIOVANNETTI_H_BN_GAP_AT_K = 4.7 # eV, p.3\n", + "GIOVANNETTI_EFFECTIVE_MASS = 4.7e-3 # m_e, (c), p.4\n", "\n", - "# meV; AA and BA read off Fig. 4 at d = 3.4 Å (±5 meV), AB is the paper's value at its 3.40 Å equilibrium\n", - "GIOVANNETTI_GAPS_AT_3_4_ANGSTROM = {\"AA\": 80, \"AB\": 46, \"BA\": 30}\n", - "GIOVANNETTI_GAPS_AT_EQUILIBRIUM = {\"AA\": (3.50, 56), \"AB\": (3.40, 46), \"BA\": (3.22, 53)}\n", - "\n", - "reference_name = next(iter(MATERIALS))\n", - "series = [\n", - " {\n", - " \"shift\": settings[\"shift\"],\n", - " \"stacking\": settings[\"stacking\"] or \"bridge\",\n", - " \"energy\": 1000 * (total_energies[name] - total_energies[reference_name]),\n", - " \"gap\": gaps[name],\n", - " }\n", - " for name, settings in MATERIALS.items()\n", - "]\n", - "print(f\"{'shift':>5} {'stacking':>8} {'ΔE (meV)':>9} {'gap (meV)':>10} {'paper gap (meV)':>16} {'deviation':>10}\")\n", - "for item in series:\n", - " paper_gap = GIOVANNETTI_GAPS_AT_3_4_ANGSTROM.get(item[\"stacking\"])\n", - " deviation = f\"{(item['gap'] - paper_gap) / paper_gap:+.0%}\" if paper_gap else \"\"\n", - " print(f\"{item['shift']:>5} {item['stacking']:>8} {item['energy']:>9.2f} {item['gap']:>10.1f} \"\n", - " f\"{paper_gap or '':>16} {deviation:>10}\")\n", + "print(f\"{'stacking':>8} {'d_eq (Å)':>9} {'paper':>6} {'deviation':>9} {'gap at d_eq (meV)':>18} {'paper':>6} \"\n", + " f\"{'deviation':>9}\")\n", + "for stacking in STACKINGS:\n", + " distance = equilibrium_distances[stacking]\n", + " paper_distance = GIOVANNETTI_EQUILIBRIUM_DISTANCE[stacking]\n", + " gap = gaps[equilibrium_names[stacking]]\n", + " paper_gap = GIOVANNETTI_GAP_AT_EQUILIBRIUM[stacking]\n", + " print(f\"{stacking:>8} {distance:>9.3f} {paper_distance:>6.2f} \"\n", + " f\"{(distance - paper_distance) / paper_distance:>+9.1%} \"\n", + " f\"{gap:>18.1f} {paper_gap:>6} {(gap - paper_gap) / paper_gap:>+9.1%}\")\n", "\n", - "plot_series(series=series, x_key=\"shift\", y_key=\"energy\", xlabel=\"Shift from BA along y (steps of a/(2√3))\",\n", - " ylabel=\"ΔE (meV per cell)\", title=\"Total energy along the sliding path\")\n", - "plot_series(series=series, x_key=\"shift\", y_key=\"gap\", xlabel=\"Shift from BA along y (steps of a/(2√3))\",\n", - " ylabel=\"Direct gap along the path (meV)\", title=\"Direct gap along the sliding path\")\n", + "print(f\"h-BN gap at K, (c): {h_bn_gap:.2f} eV, paper {GIOVANNETTI_H_BN_GAP_AT_K} eV\")\n", + "print(f\"Effective mass at K, (c): {effective_mass:.2e} m_e, paper {GIOVANNETTI_EFFECTIVE_MASS:.1e} m_e\")\n", "\n", - "print(\"Gaps at the paper's equilibrium distances: \" + \", \".join(\n", - " f\"{stacking} {distance:.2f} Å {gap} meV\" for stacking, (distance, gap) in GIOVANNETTI_GAPS_AT_EQUILIBRIUM.items()))" + "energy_order = [stacking for stacking in [\"c\", \"b\", \"a\"] if stacking in STACKINGS]\n", + "energy_order_holds = all(\n", + " total_energies[material_names[lower][distance]] < total_energies[material_names[higher][distance]]\n", + " for lower, higher in zip(energy_order, energy_order[1:])\n", + " for distance in DISTANCES\n", + ")\n", + "print(f\"E(c) < E(b) < E(a) at every distance: {energy_order_holds}\")" ] }, { "cell_type": "markdown", - "id": "44", + "id": "54", "metadata": {}, "source": [ "## References\n", "\n", - "[1] G. Giovannetti, P. A. Khomyakov, G. Brocks, P. J. Kelly and J. van den Brink, \"Substrate-induced band gap in graphene on hexagonal boron nitride: Ab initio density functional calculations\", Phys. Rev. B 76, 073103 (2007). https://doi.org/10.1103/PhysRevB.76.073103\n", - "\n", - "[2] J. Jung, A. M. DaSilva, A. H. MacDonald and S. Adam, \"Origin of band gaps in graphene on hexagonal boron nitride\", Nat. Commun. 6, 6308 (2015). https://doi.org/10.1038/ncomms7308" + "[1] G. Giovannetti, P. A. Khomyakov, G. Brocks, P. J. Kelly and J. van den Brink, \"Substrate-induced band gap in graphene on hexagonal boron nitride: Ab initio density functional calculations\", Phys. Rev. B 76, 073103 (2007). https://doi.org/10.1103/PhysRevB.76.073103" ] } ], From 7ff17e590f22bb5e19bf2ffe7bf292eec370db8b Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Wed, 7 Oct 2026 19:55:48 -0700 Subject: [PATCH 13/17] =?UTF-8?q?SOF-8064:=20simulation=20notebook=20?= =?UTF-8?q?=E2=80=94=20E(d),=20gap(d),=20bands=20and=20DOS=20at=20the=20eq?= =?UTF-8?q?uilibrium=20distance,=20as=20Giovannetti=202007?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Opus 5.5 --- ...nterface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb index 141ed0214..ecc91d35c 100644 --- a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb @@ -128,7 +128,6 @@ "\n", "KGRID = [36, 36, 1] # Giovannetti et al. 2007; a multiple of 3 keeps K on the mesh\n", "OCCUPATIONS_SETTINGS = {\"occupations\": \"tetrahedra\"} # the paper's tetrahedron method\n", - "# Dipole correction along z; emaxpos is set per material in 4.2, in the middle of the vacuum above graphene\n", "DIPOLE_SETTINGS = {\"control\": {\"tefield\": True, \"dipfield\": True}, \"system\": {\"edir\": 3, \"eopreg\": 0.05}}\n", "KPATH_STEPS = 40\n", "MODEL_TAG = f\"{FUNCTIONAL}-{PSEUDOPOTENTIAL_TYPE} {ECUTWFC}-{ECUTRHO}Ry k{KGRID[0]} p{KPATH_STEPS} tetra dip\"\n", @@ -393,8 +392,7 @@ "\n", "\n", "def get_vacuum_center(material):\n", - " heights = [coordinate[2] for coordinate in material.basis.coordinates.values]\n", - " return round((max(heights) + 1 + min(heights)) / 2, 4)\n", + " return (max(coordinate[2] for coordinate in material.basis.coordinates.values) + 1) / 2\n", "\n", "\n", "workflows = {}\n", @@ -713,9 +711,9 @@ "bands = np.array(band_structures[name][\"yDataSeries\"])\n", "path_indices = np.arange(K_VERTEX_INDEX - 8, K_VERTEX_INDEX + 9)\n", "center = bands[NUMBER_OF_OCCUPIED_BANDS - 1 : NUMBER_OF_OCCUPIED_BANDS + 1, K_VERTEX_INDEX].mean()\n", - "figure, axes = plt.subplots(figsize=(6, 5))\n", + "figure, axes = plt.subplots()\n", "for band in bands[NUMBER_OF_OCCUPIED_BANDS - 1 : NUMBER_OF_OCCUPIED_BANDS + 1]:\n", - " axes.plot(path_indices, band[path_indices], marker=\"o\")\n", + " axes.plot(path_indices, band[path_indices])\n", "axes.set_ylim(center - 0.3, center + 0.3)\n", "axes.set_xlabel(f\"Path index (K at {K_VERTEX_INDEX})\")\n", "axes.set_ylabel(\"Energy (eV)\")\n", From 90a9ed955a758b992580ec1b8f73a5ef4b483235 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Wed, 7 Oct 2026 20:01:19 -0700 Subject: [PATCH 14/17] SOF-8064: effective mass from the gapped-cone relation; 100 steps per path leg Co-Authored-By: Claude Opus 5.5 --- ...2d_boron_nitride_graphene_SIMULATION.ipynb | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb index ecc91d35c..1a8324117 100644 --- a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb @@ -129,7 +129,7 @@ "KGRID = [36, 36, 1] # Giovannetti et al. 2007; a multiple of 3 keeps K on the mesh\n", "OCCUPATIONS_SETTINGS = {\"occupations\": \"tetrahedra\"} # the paper's tetrahedron method\n", "DIPOLE_SETTINGS = {\"control\": {\"tefield\": True, \"dipfield\": True}, \"system\": {\"edir\": 3, \"eopreg\": 0.05}}\n", - "KPATH_STEPS = 40\n", + "KPATH_STEPS = 100\n", "MODEL_TAG = f\"{FUNCTIONAL}-{PSEUDOPOTENTIAL_TYPE} {ECUTWFC}-{ECUTRHO}Ry k{KGRID[0]} p{KPATH_STEPS} tetra dip\"\n", "\n", "SCF_UNIT = \"pw_scf\"\n", @@ -709,12 +709,12 @@ " title=f\"Density of States: {name}\", extra_config={\"material\": materials[name].to_dict()})\n", "\n", "bands = np.array(band_structures[name][\"yDataSeries\"])\n", - "path_indices = np.arange(K_VERTEX_INDEX - 8, K_VERTEX_INDEX + 9)\n", + "path_indices = np.arange(K_VERTEX_INDEX - 12, K_VERTEX_INDEX + 13)\n", "center = bands[NUMBER_OF_OCCUPIED_BANDS - 1 : NUMBER_OF_OCCUPIED_BANDS + 1, K_VERTEX_INDEX].mean()\n", "figure, axes = plt.subplots()\n", "for band in bands[NUMBER_OF_OCCUPIED_BANDS - 1 : NUMBER_OF_OCCUPIED_BANDS + 1]:\n", " axes.plot(path_indices, band[path_indices])\n", - "axes.set_ylim(center - 0.3, center + 0.3)\n", + "axes.set_ylim(center - 0.5, center + 0.5)\n", "axes.set_xlabel(f\"Path index (K at {K_VERTEX_INDEX})\")\n", "axes.set_ylabel(\"Energy (eV)\")\n", "axes.set_title(\"Bands around K\")\n", @@ -748,8 +748,8 @@ "source": [ "### 7.8. Effective mass at K\n", "\n", - "The curvature of the lowest unoccupied band at K, from the parabola through the two path points on\n", - "either side of it." + "For a gapped Dirac cone, E(k) = ±√((Δ/2)² + (ħv k)²), the band-edge mass is m* = ħ²Δ / (2(ħv)²), with Δ\n", + "the gap at K from 7.2 and ħv the slope of the lowest unoccupied band against |k − K| on the Γ–K leg." ] }, { @@ -762,11 +762,10 @@ "HBAR_SQUARED_OVER_ELECTRON_MASS = 7.6200 # eV·Å²\n", "K_STEP = 4 * np.pi / (3 * materials[name].lattice.a) / KPATH_STEPS # Å⁻¹ per path step along Γ–K\n", "\n", - "# the K–M steps are half as long as the Γ–K steps\n", - "k_offsets = K_STEP * np.array([-2, -1, 0, 0.5, 1])\n", - "quadratic, _, _ = np.polyfit(k_offsets, bands[NUMBER_OF_OCCUPIED_BANDS, K_VERTEX_INDEX - 2 : K_VERTEX_INDEX + 3], 2)\n", - "effective_mass = HBAR_SQUARED_OVER_ELECTRON_MASS / (2 * quadratic)\n", - "print(f\"{name}: effective mass at K {effective_mass:.2e} m_e\")" + "k_offsets = K_STEP * np.arange(8, 2, -1)\n", + "hbar_velocity, _ = np.polyfit(k_offsets, bands[NUMBER_OF_OCCUPIED_BANDS, K_VERTEX_INDEX - 8 : K_VERTEX_INDEX - 2], 1)\n", + "effective_mass = HBAR_SQUARED_OVER_ELECTRON_MASS * gaps[name] / 1000 / (2 * hbar_velocity**2)\n", + "print(f\"{name}: ħv = {hbar_velocity:.2f} eV·Å, effective mass at K {effective_mass:.2e} m_e\")" ] }, { From 2f6a178cea644f5172bba64ced135d2c8dd4d762 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Wed, 7 Oct 2026 20:04:06 -0700 Subject: [PATCH 15/17] =?UTF-8?q?SOF-8064:=20tetrahedra=20on=20the=20scf?= =?UTF-8?q?=20and=20nscf=20units=20only=20=E2=80=94=20pw=5Fbands=20keeps?= =?UTF-8?q?=20the=20template's=20smearing?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Sonnet 5.5 --- .../interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb index 1a8324117..63c019b91 100644 --- a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb @@ -410,7 +410,8 @@ " unit.add_context(context)\n", " subworkflow.set_unit(unit)\n", " patch_workflow_qe_input(workflow, DIPOLE_SETTINGS, [SCF_UNIT, NSCF_UNIT, BANDS_UNIT])\n", - " patch_workflow_qe_input(workflow, {\"system\": {**OCCUPATIONS_SETTINGS, \"emaxpos\": get_vacuum_center(material)}},\n", + " patch_workflow_qe_input(workflow, {\"system\": OCCUPATIONS_SETTINGS}, [SCF_UNIT, NSCF_UNIT])\n", + " patch_workflow_qe_input(workflow, {\"system\": {\"emaxpos\": get_vacuum_center(material)}},\n", " [SCF_UNIT, NSCF_UNIT, BANDS_UNIT])\n", " workflows[name] = workflow\n", " print(workflow.name)" From ec8349892c130a3c774a831de37990e1027bb37b Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Wed, 7 Oct 2026 20:48:22 -0700 Subject: [PATCH 16/17] =?UTF-8?q?SOF-8064:=20review=20round=202=20?= =?UTF-8?q?=E2=80=94=20eamp=200=20(no=20external=20field),=20paper-style?= =?UTF-8?q?=20Figs.=202/4,=20=C4=A7v=20from=20both=20sides,=20interpolated?= =?UTF-8?q?=20gap=20at=20d=5Feq,=20no=20grading=20line?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - DIPOLE_SETTINGS sets eamp = 0.0: with tefield and no eamp, pw.x applied its default 0.001 a.u. (0.051 V/Å) on top of the dipole correction. MODEL_TAG is now built from the settings ("… tetrahedra eamp0.0"), so the field-on jobs named "… tetra dip" are never reused. - §8 no longer prints the E(c) < E(b) < E(a) boolean (it printed True for stackings that never ran); the ordering is read off Fig. 2. - Figs. 2 and 4: one axis, one curve per active stacking. Fig. 2 has one common zero, E of (c) at the largest distance, or of the first active stacking when (c) is off. - ħv is the mean of the Γ–K and K–M slope fits over VELOCITY_FIT_RANGE, with step lengths 4π/(3a)/KPATH_STEPS and 2π/(3a)/KPATH_STEPS. - Gap at d_eq interpolated between the grid distances; Fig. 3 stays at the nearest grid point. - Zoom plotted against k − K in Å⁻¹; ZOOM_POINTS, ZOOM_WINDOW and VELOCITY_FIT_RANGE in §1.3; `name` renamed to `equilibrium_name_c`. - Docstrings on get_registry and get_vacuum_center; REGISTRY_TOLERANCE in the structure notebook's parameters. - PPN = 2, queue D's maximum on cluster-001: one core measured 31 min per job. Not applied: scipy.constants for ħ²/mₑ, made's get_atom_indices_by_layer in get_registry, extending the "settings that differ" note. Co-Authored-By: Claude Opus 5.5 --- ...terface_2d_2d_boron_nitride_graphene.ipynb | 4 +- ...2d_boron_nitride_graphene_SIMULATION.ipynb | 114 ++++++++++-------- 2 files changed, 69 insertions(+), 49 deletions(-) diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb index 8348bf116..44fa8b893 100644 --- a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb @@ -47,6 +47,7 @@ " # \"b\",\n", "]\n", "STACKING_SHIFTS = {\"a\": 0, \"b\": 1, \"c\": -1} # in units of a/√3 along y\n", + "REGISTRY_TOLERANCE = 1e-3 # crystal units\n", "# Graphene–h-BN distance, Å. The paper scans 2.5–3.9; three points around (c)'s minimum are active, uncomment\n", "# the rest for the full set.\n", "DISTANCES = [\n", @@ -221,13 +222,14 @@ "\n", "\n", "def get_registry(interface):\n", + " \"\"\"Atom of the top h-BN layer under each carbon atom, or \"hollow\" (Giovannetti et al. 2007 Fig. 1).\"\"\"\n", " elements = np.array(interface.basis.elements.values)\n", " coordinates = np.array(interface.basis.coordinates.values)\n", " top_layer = (elements != \"C\") & np.isclose(coordinates[:, 2], coordinates[elements != \"C\", 2].max())\n", " registry = []\n", " for carbon in coordinates[elements == \"C\"]:\n", " in_plane_offsets = (coordinates[top_layer, :2] - carbon[:2] + 0.5) % 1 - 0.5\n", - " atoms_below = elements[top_layer][np.all(np.abs(in_plane_offsets) < 1e-3, axis=1)]\n", + " atoms_below = elements[top_layer][np.all(np.abs(in_plane_offsets) < REGISTRY_TOLERANCE, axis=1)]\n", " registry.append(atoms_below[0] if len(atoms_below) else \"hollow\")\n", " return registry\n", "\n", diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb index 63c019b91..ad54c1528 100644 --- a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb @@ -98,7 +98,7 @@ "\n", "CLUSTER_NAME = \"001\" # specify full or partial name i.e. \"cluster-001\" to select\n", "QUEUE_NAME = QueueName.D\n", - "PPN = 1\n", + "PPN = 2\n", "TIME_LIMIT = \"04:00:00\"\n", "\n", "timestamp = datetime.now().strftime(\"%Y-%m-%d %H:%M\")\n", @@ -128,9 +128,11 @@ "\n", "KGRID = [36, 36, 1] # Giovannetti et al. 2007; a multiple of 3 keeps K on the mesh\n", "OCCUPATIONS_SETTINGS = {\"occupations\": \"tetrahedra\"} # the paper's tetrahedron method\n", - "DIPOLE_SETTINGS = {\"control\": {\"tefield\": True, \"dipfield\": True}, \"system\": {\"edir\": 3, \"eopreg\": 0.05}}\n", + "# eamp = 0: the sawtooth is the dipole correction only, no external field\n", + "DIPOLE_SETTINGS = {\"control\": {\"tefield\": True, \"dipfield\": True}, \"system\": {\"edir\": 3, \"eamp\": 0.0, \"eopreg\": 0.05}}\n", "KPATH_STEPS = 100\n", - "MODEL_TAG = f\"{FUNCTIONAL}-{PSEUDOPOTENTIAL_TYPE} {ECUTWFC}-{ECUTRHO}Ry k{KGRID[0]} p{KPATH_STEPS} tetra dip\"\n", + "MODEL_TAG = (f\"{FUNCTIONAL}-{PSEUDOPOTENTIAL_TYPE} {ECUTWFC}-{ECUTRHO}Ry k{KGRID[0]} p{KPATH_STEPS} \"\n", + " f\"{OCCUPATIONS_SETTINGS['occupations']} eamp{DIPOLE_SETTINGS['system']['eamp']}\")\n", "\n", "SCF_UNIT = \"pw_scf\"\n", "NSCF_UNIT = \"pw_nscf\"\n", @@ -142,7 +144,11 @@ " {\"point\": \"K\", \"steps\": KPATH_STEPS},\n", " {\"point\": \"M\", \"steps\": KPATH_STEPS},\n", " {\"point\": \"Γ\", \"steps\": 1},\n", - "]" + "]\n", + "\n", + "VELOCITY_FIT_RANGE = (3, 8) # path points from K used for the ħv fit\n", + "ZOOM_POINTS = 12 # path points each side of K in the zoom\n", + "ZOOM_WINDOW = 0.5 # eV each side of the band edges" ] }, { @@ -392,6 +398,7 @@ "\n", "\n", "def get_vacuum_center(material):\n", + " \"\"\"Crystal-z midpoint of the vacuum above the film, QE's `emaxpos` unit.\"\"\"\n", " return (max(coordinate[2] for coordinate in material.basis.coordinates.values) + 1) / 2\n", "\n", "\n", @@ -610,8 +617,9 @@ "source": [ "### 7.3. Equilibrium distance of each stacking\n", "\n", - "The vertex of the parabola through the three distances with the lowest total energy; the gap at K is\n", - "printed for the distance in `DISTANCES` nearest to it." + "The vertex of the parabola through the three distances with the lowest total energy, and the gap at K\n", + "interpolated to it from 7.2. The bands of 7.6 (Fig. 3) are those of the distance in `DISTANCES`\n", + "nearest to it." ] }, { @@ -623,14 +631,17 @@ "source": [ "equilibrium_distances = {}\n", "equilibrium_names = {}\n", + "gap_at_equilibrium = {}\n", "for stacking, names in material_names.items():\n", " lowest = sorted(names, key=lambda distance: total_energies[names[distance]])[:3]\n", " quadratic, linear, _ = np.polyfit(lowest, [total_energies[names[distance]] for distance in lowest], 2)\n", " equilibrium_distances[stacking] = -linear / (2 * quadratic)\n", " nearest = min(names, key=lambda distance: abs(distance - equilibrium_distances[stacking]))\n", " equilibrium_names[stacking] = names[nearest]\n", + " gap_at_equilibrium[stacking] = np.interp(\n", + " equilibrium_distances[stacking], list(names), [gaps[name] for name in names.values()])\n", " print(f\"({stacking}): equilibrium distance {equilibrium_distances[stacking]:.3f} Å, \"\n", - " f\"gap at K at {nearest:.2f} Å: {gaps[names[nearest]]:.1f} meV\")" + " f\"gap at K interpolated to it {gap_at_equilibrium[stacking]:.1f} meV\")" ] }, { @@ -648,16 +659,20 @@ "metadata": {}, "outputs": [], "source": [ - "from mat3ra.notebooks_utils.plot import plot_series\n", + "from matplotlib import pyplot as plt\n", + "from mat3ra.notebooks_utils.plot import display_matplotlib_figure\n", "\n", + "reference_stacking = \"c\" if \"c\" in STACKINGS else STACKINGS[0]\n", + "reference_energy = total_energies[material_names[reference_stacking][max(DISTANCES)]]\n", + "figure, axes = plt.subplots()\n", "for stacking, names in material_names.items():\n", - " reference_energy = total_energies[names[max(names)]]\n", - " series = [\n", - " {\"distance\": distance, \"energy\": total_energies[name] - reference_energy} for distance, name in names.items()\n", - " ]\n", - " plot_series(series=series, x_key=\"distance\", y_key=\"energy\", xlabel=\"Graphene–h-BN distance (Å)\",\n", - " ylabel=f\"E − E({max(names):.2f} Å) (eV per cell)\",\n", - " title=f\"Total energy vs distance, stacking ({stacking})\")" + " axes.plot(list(names), [total_energies[name] - reference_energy for name in names.values()], marker=\"o\",\n", + " label=f\"({stacking})\")\n", + "axes.set_xlabel(\"Graphene–h-BN distance (Å)\")\n", + "axes.set_ylabel(f\"E − E(({reference_stacking}), {max(DISTANCES):.2f} Å) (eV per cell)\")\n", + "axes.set_title(\"Total energy vs distance\")\n", + "axes.legend()\n", + "display_matplotlib_figure(figure)" ] }, { @@ -675,10 +690,14 @@ "metadata": {}, "outputs": [], "source": [ + "figure, axes = plt.subplots()\n", "for stacking, names in material_names.items():\n", - " series = [{\"distance\": distance, \"gap\": gaps[name]} for distance, name in names.items()]\n", - " plot_series(series=series, x_key=\"distance\", y_key=\"gap\", xlabel=\"Graphene–h-BN distance (Å)\",\n", - " ylabel=\"Direct gap at K (meV)\", title=f\"Gap at K vs distance, stacking ({stacking})\")" + " axes.plot(list(names), [gaps[name] for name in names.values()], marker=\"o\", label=f\"({stacking})\")\n", + "axes.set_xlabel(\"Graphene–h-BN distance (Å)\")\n", + "axes.set_ylabel(\"Gap at K (meV)\")\n", + "axes.set_title(\"Gap at K vs distance\")\n", + "axes.legend()\n", + "display_matplotlib_figure(figure)" ] }, { @@ -699,24 +718,27 @@ "metadata": {}, "outputs": [], "source": [ - "from matplotlib import pyplot as plt\n", "from mat3ra.notebooks_utils.ipython.entity.property.visualize import visualize_properties\n", - "from mat3ra.notebooks_utils.plot import display_matplotlib_figure\n", "\n", - "name = equilibrium_names[\"c\"]\n", - "visualize_properties(band_structures[name], title=f\"Band Structure: {name}\",\n", - " extra_config={\"material\": materials[name].to_dict()})\n", - "visualize_properties(get_properties_for_job(client, job_ids[name], property_name=\"density_of_states\"),\n", - " title=f\"Density of States: {name}\", extra_config={\"material\": materials[name].to_dict()})\n", + "equilibrium_name_c = equilibrium_names[\"c\"]\n", + "visualize_properties(band_structures[equilibrium_name_c], title=f\"Band Structure: {equilibrium_name_c}\",\n", + " extra_config={\"material\": materials[equilibrium_name_c].to_dict()})\n", + "visualize_properties(get_properties_for_job(client, job_ids[equilibrium_name_c], property_name=\"density_of_states\"),\n", + " title=f\"Density of States: {equilibrium_name_c}\",\n", + " extra_config={\"material\": materials[equilibrium_name_c].to_dict()})\n", + "\n", + "K_STEP_GAMMA_K = 4 * np.pi / (3 * materials[equilibrium_name_c].lattice.a) / KPATH_STEPS # Å⁻¹ per Γ–K step\n", + "K_STEP_K_M = 2 * np.pi / (3 * materials[equilibrium_name_c].lattice.a) / KPATH_STEPS # Å⁻¹ per K–M step\n", "\n", - "bands = np.array(band_structures[name][\"yDataSeries\"])\n", - "path_indices = np.arange(K_VERTEX_INDEX - 12, K_VERTEX_INDEX + 13)\n", - "center = bands[NUMBER_OF_OCCUPIED_BANDS - 1 : NUMBER_OF_OCCUPIED_BANDS + 1, K_VERTEX_INDEX].mean()\n", + "bands = np.array(band_structures[equilibrium_name_c][\"yDataSeries\"])\n", + "path_steps = np.arange(-ZOOM_POINTS, ZOOM_POINTS + 1)\n", + "k_offsets = path_steps * np.where(path_steps < 0, K_STEP_GAMMA_K, K_STEP_K_M)\n", "figure, axes = plt.subplots()\n", "for band in bands[NUMBER_OF_OCCUPIED_BANDS - 1 : NUMBER_OF_OCCUPIED_BANDS + 1]:\n", - " axes.plot(path_indices, band[path_indices])\n", - "axes.set_ylim(center - 0.5, center + 0.5)\n", - "axes.set_xlabel(f\"Path index (K at {K_VERTEX_INDEX})\")\n", + " axes.plot(k_offsets, band[K_VERTEX_INDEX + path_steps])\n", + "axes.set_ylim(bands[NUMBER_OF_OCCUPIED_BANDS - 1, K_VERTEX_INDEX] - ZOOM_WINDOW,\n", + " bands[NUMBER_OF_OCCUPIED_BANDS, K_VERTEX_INDEX] + ZOOM_WINDOW)\n", + "axes.set_xlabel(\"k − K (Å⁻¹), Γ side < 0 < M side\")\n", "axes.set_ylabel(\"Energy (eV)\")\n", "axes.set_title(\"Bands around K\")\n", "display_matplotlib_figure(figure)" @@ -739,7 +761,7 @@ "source": [ "# the h-BN bands adjacent to the graphene pair\n", "h_bn_gap = bands[NUMBER_OF_OCCUPIED_BANDS + 1, K_VERTEX_INDEX] - bands[NUMBER_OF_OCCUPIED_BANDS - 2, K_VERTEX_INDEX]\n", - "print(f\"{name}: h-BN gap at K {h_bn_gap:.2f} eV\")" + "print(f\"{equilibrium_name_c}: h-BN gap at K {h_bn_gap:.2f} eV\")" ] }, { @@ -750,7 +772,8 @@ "### 7.8. Effective mass at K\n", "\n", "For a gapped Dirac cone, E(k) = ±√((Δ/2)² + (ħv k)²), the band-edge mass is m* = ħ²Δ / (2(ħv)²), with Δ\n", - "the gap at K from 7.2 and ħv the slope of the lowest unoccupied band against |k − K| on the Γ–K leg." + "the gap at K from 7.2 and ħv the slope of the lowest unoccupied band against |k − K|, the mean of the\n", + "Γ–K and the K–M sides of K over `VELOCITY_FIT_RANGE`." ] }, { @@ -761,12 +784,15 @@ "outputs": [], "source": [ "HBAR_SQUARED_OVER_ELECTRON_MASS = 7.6200 # eV·Å²\n", - "K_STEP = 4 * np.pi / (3 * materials[name].lattice.a) / KPATH_STEPS # Å⁻¹ per path step along Γ–K\n", "\n", - "k_offsets = K_STEP * np.arange(8, 2, -1)\n", - "hbar_velocity, _ = np.polyfit(k_offsets, bands[NUMBER_OF_OCCUPIED_BANDS, K_VERTEX_INDEX - 8 : K_VERTEX_INDEX - 2], 1)\n", - "effective_mass = HBAR_SQUARED_OVER_ELECTRON_MASS * gaps[name] / 1000 / (2 * hbar_velocity**2)\n", - "print(f\"{name}: ħv = {hbar_velocity:.2f} eV·Å, effective mass at K {effective_mass:.2e} m_e\")" + "fit_steps = np.arange(VELOCITY_FIT_RANGE[0], VELOCITY_FIT_RANGE[1] + 1)\n", + "conduction_band = bands[NUMBER_OF_OCCUPIED_BANDS]\n", + "hbar_velocity_gamma_k, _ = np.polyfit(K_STEP_GAMMA_K * fit_steps, conduction_band[K_VERTEX_INDEX - fit_steps], 1)\n", + "hbar_velocity_k_m, _ = np.polyfit(K_STEP_K_M * fit_steps, conduction_band[K_VERTEX_INDEX + fit_steps], 1)\n", + "hbar_velocity = (hbar_velocity_gamma_k + hbar_velocity_k_m) / 2\n", + "effective_mass = HBAR_SQUARED_OVER_ELECTRON_MASS * gaps[equilibrium_name_c] / 1000 / (2 * hbar_velocity**2)\n", + "print(f\"{equilibrium_name_c}: ħv = {hbar_velocity_gamma_k:.2f} eV·Å (Γ–K), {hbar_velocity_k_m:.2f} eV·Å (K–M), \"\n", + " f\"mean {hbar_velocity:.2f} eV·Å, effective mass at K {effective_mass:.2e} m_e\")" ] }, { @@ -799,22 +825,14 @@ "for stacking in STACKINGS:\n", " distance = equilibrium_distances[stacking]\n", " paper_distance = GIOVANNETTI_EQUILIBRIUM_DISTANCE[stacking]\n", - " gap = gaps[equilibrium_names[stacking]]\n", + " gap = gap_at_equilibrium[stacking]\n", " paper_gap = GIOVANNETTI_GAP_AT_EQUILIBRIUM[stacking]\n", " print(f\"{stacking:>8} {distance:>9.3f} {paper_distance:>6.2f} \"\n", " f\"{(distance - paper_distance) / paper_distance:>+9.1%} \"\n", " f\"{gap:>18.1f} {paper_gap:>6} {(gap - paper_gap) / paper_gap:>+9.1%}\")\n", "\n", "print(f\"h-BN gap at K, (c): {h_bn_gap:.2f} eV, paper {GIOVANNETTI_H_BN_GAP_AT_K} eV\")\n", - "print(f\"Effective mass at K, (c): {effective_mass:.2e} m_e, paper {GIOVANNETTI_EFFECTIVE_MASS:.1e} m_e\")\n", - "\n", - "energy_order = [stacking for stacking in [\"c\", \"b\", \"a\"] if stacking in STACKINGS]\n", - "energy_order_holds = all(\n", - " total_energies[material_names[lower][distance]] < total_energies[material_names[higher][distance]]\n", - " for lower, higher in zip(energy_order, energy_order[1:])\n", - " for distance in DISTANCES\n", - ")\n", - "print(f\"E(c) < E(b) < E(a) at every distance: {energy_order_holds}\")" + "print(f\"Effective mass at K, (c): {effective_mass:.2e} m_e, paper {GIOVANNETTI_EFFECTIVE_MASS:.1e} m_e\")" ] }, { From bbd3e427f7992cf71d14865bd0633790b3775518 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Wed, 7 Oct 2026 23:56:57 -0700 Subject: [PATCH 17/17] =?UTF-8?q?SOF-8064:=20one=20note=20on=20the=20stack?= =?UTF-8?q?ing=20and=20distance=20lists=20=E2=80=94=20what=20the=20paper?= =?UTF-8?q?=20computes,=20what=20the=20defaults=20run,=20how=20to=20run=20?= =?UTF-8?q?more?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Sonnet 5.5 --- .../interface_2d_2d_boron_nitride_graphene.ipynb | 14 ++++++-------- ...e_2d_2d_boron_nitride_graphene_SIMULATION.ipynb | 3 +++ 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb index 44fa8b893..a20b1e6a7 100644 --- a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene.ipynb @@ -38,22 +38,20 @@ "H_BN_INTERLAYER_DISTANCE = 3.24 # Å, the paper's LDA value\n", "VACUUM = 15.0 # Å above graphene\n", "\n", - "# Registry of graphene on the top h-BN layer, Giovannetti et al. 2007 Fig. 1: one C over B and the other over\n", - "# N (a), over N and a hexagon centre (b), over B and a hexagon centre (c). All three run in the paper;\n", - "# uncomment to build the others.\n", + "# Giovannetti et al. 2007 compute all three stackings at every distance 2.5–3.9 Å; the defaults run\n", + "# stacking (c) at three distances around its minimum. Uncomment entries to compute more\n", + "# (one job per stacking × distance, ~25 min each on queue D with two cores).\n", "STACKINGS = [\n", " \"c\",\n", " # \"a\",\n", " # \"b\",\n", "]\n", - "STACKING_SHIFTS = {\"a\": 0, \"b\": 1, \"c\": -1} # in units of a/√3 along y\n", - "REGISTRY_TOLERANCE = 1e-3 # crystal units\n", - "# Graphene–h-BN distance, Å. The paper scans 2.5–3.9; three points around (c)'s minimum are active, uncomment\n", - "# the rest for the full set.\n", "DISTANCES = [\n", " 3.1, 3.2, 3.3,\n", " # 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9,\n", - "]" + "]\n", + "STACKING_SHIFTS = {\"a\": 0, \"b\": 1, \"c\": -1} # in units of a/√3 along y\n", + "REGISTRY_TOLERANCE = 1e-3 # crystal units" ] }, { diff --git a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb index ad54c1528..adc2cac6a 100644 --- a/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb +++ b/other/materials_designer/specific_examples/interface_2d_2d_boron_nitride_graphene_SIMULATION.ipynb @@ -81,6 +81,9 @@ "ORGANIZATION_NAME = None # set to your organization name (full or partial); otherwise, your default one is used\n", "FOLDER = \"./uploads\"\n", "\n", + "# Giovannetti et al. 2007 compute all three stackings at every distance 2.5–3.9 Å; the defaults run\n", + "# stacking (c) at three distances around its minimum. Uncomment entries to compute more\n", + "# (one job per stacking × distance, ~25 min each on queue D with two cores).\n", "STACKINGS = [\n", " \"c\",\n", " # \"a\",\n",