diff --git a/other/experiments/jupyterlite/relax_structure_with_mlff.ipynb b/other/experiments/jupyterlite/relax_structure_with_mlff.ipynb index 470169162..3f77085b2 100644 --- a/other/experiments/jupyterlite/relax_structure_with_mlff.ipynb +++ b/other/experiments/jupyterlite/relax_structure_with_mlff.ipynb @@ -82,7 +82,7 @@ "source": [ "from mat3ra.notebooks_utils.packages import install_packages\n", "from mat3ra.notebooks_utils.primitive.environment import is_pyodide_environment\n", - "from mat3ra.notebooks_utils.mlff import get_mlff_install_profiles\n", + "from mat3ra.notebooks_utils.calculators.mlff import get_mlff_install_profiles\n", "\n", "profiles = get_mlff_install_profiles(MLFF_NAME)\n", "await install_packages(profiles)\n", @@ -157,7 +157,7 @@ "from mat3ra.made.tools.convert import to_ase\n", "from ase.optimize import BFGS\n", "\n", - "from mat3ra.notebooks_utils.mlff import create_mlff_calculator\n", + "from mat3ra.notebooks_utils.calculators.mlff import create_mlff_calculator\n", "from mat3ra.notebooks_utils.ipython.plot._plotly import progress_callback\n", "\n", "calculator = create_mlff_calculator(MLFF_NAME, MLFF_SETTINGS[MLFF_NAME])\n", diff --git a/other/materials_designer/specific_examples/Introduction.ipynb b/other/materials_designer/specific_examples/Introduction.ipynb index 4795f3c84..bed9c1ffc 100644 --- a/other/materials_designer/specific_examples/Introduction.ipynb +++ b/other/materials_designer/specific_examples/Introduction.ipynb @@ -27,7 +27,7 @@ "| `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 | [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) | *To be added* | [[7]](#ref7) |\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", "| `C-2D-INT-T` | Interface Twisted | *To be added* | — | — |\n", "| `C-2D-INT-C` | Interface Commensurate Lattice | [Twisted Commensurate MoS₂ Bilayer](interface_bilayer_twisted_commensurate_lattices_molybdenum_disulfide.ipynb) | [Twisted MoS₂ Bilayer Band Structure](interface_bilayer_twisted_commensurate_lattices_molybdenum_disulfide_SIMULATION.ipynb) | [[8]](#ref8) |\n", "| `C-2D-MLT` | Multi-Layer | *To be added* | — | — |\n", diff --git a/other/materials_designer/specific_examples/optimization_interface_film_xy_position_graphene_nickel.ipynb b/other/materials_designer/specific_examples/optimization_interface_film_xy_position_graphene_nickel.ipynb index f8ad43b81..60aaf1527 100644 --- a/other/materials_designer/specific_examples/optimization_interface_film_xy_position_graphene_nickel.ipynb +++ b/other/materials_designer/specific_examples/optimization_interface_film_xy_position_graphene_nickel.ipynb @@ -56,6 +56,8 @@ "# Material selection\n", "SUBSTRATE_NAME = \"Nickel\"\n", "FILM_NAME = \"Graphene\"\n", + "# canonical name the companion SIMULATION notebook loads by\n", + "BASE_MATERIAL_NAME = \"Graphene_Nickel_interface\"\n", "\n", "# Slab parameters\n", "FILM_MILLER_INDICES = (0, 0, 1)\n", @@ -196,7 +198,7 @@ " reduce_result_cell_to_primitive=REDUCE_RESULT_CELL_TO_PRIMITIVE,\n", ")\n", "\n", - "interface_material.name = f\"{FILM_NAME}_{SUBSTRATE_NAME}_interface\"\n", + "interface_material.name = BASE_MATERIAL_NAME\n", "\n", "# Visualize interface\n", "visualize_materials([interface_material], repetitions=STRUCTURE_REPETITIONS)\n", @@ -264,7 +266,7 @@ "metadata": {}, "outputs": [], "source": [ - "from mat3ra.utils.jupyterlite.plot import plot_2d_heatmap, plot_3d_surface\n", + "from mat3ra.notebooks_utils.ipython.plot._plotly import plot_2d_heatmap, plot_3d_surface\n", "\n", "x_values, y_values = xy_matrix\n", "# Plot energy landscape\n", @@ -317,7 +319,7 @@ "id": "16", "metadata": {}, "source": [ - "# 4. Save optimized material" + "# 4. Save the base and optimized materials" ] }, { @@ -330,6 +332,9 @@ "from mat3ra.notebooks_utils.io import download_content_to_file\n", "from mat3ra.notebooks_utils.material import set_materials\n", "\n", + "set_materials(interface_material)\n", + "download_content_to_file(interface_material.to_json(), f\"{interface_material.name}.json\")\n", + "\n", "optimized_material.name = f\"{interface_material.name}_optimized_xy\"\n", "set_materials(optimized_material)\n", "download_content_to_file(optimized_material.to_json(), f\"{interface_material.name}_optimized_xy.json\")" diff --git a/other/materials_designer/specific_examples/optimization_interface_film_xy_position_graphene_nickel_SIMULATION.ipynb b/other/materials_designer/specific_examples/optimization_interface_film_xy_position_graphene_nickel_SIMULATION.ipynb new file mode 100644 index 000000000..3fb5e1926 --- /dev/null +++ b/other/materials_designer/specific_examples/optimization_interface_film_xy_position_graphene_nickel_SIMULATION.ipynb @@ -0,0 +1,1109 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Graphene/Ni(111) Interface: Registry, Separation and Work of Adhesion\n", + "\n", + "## 0. Introduction\n", + "\n", + "This notebook reproduces the structure and energetics of graphene on Ni(111) following the review:\n", + "\n", + "> **Arjun Dahal, Matthias Batzill**\n", + "> \"Graphene–nickel interfaces: a review\"\n", + "> Nanoscale, 6(5), 2548. (2014)\n", + "> [DOI: 10.1039/c3nr05279f](https://doi.org/10.1039/c3nr05279f)\n", + "\n", + "The review's structural facts (its section 2.1): graphene locks into a 1×1 registry on Ni(111);\n", + "LEED I–V and ion scattering identify the adsorbed structure as one carbon **atop** a first-layer Ni\n", + "and the other in the **fcc hollow**, 0.211 nm above the surface with a 0.005 nm buckling in which\n", + "the atop carbon sits further out. Its computed numbers come from\n", + "[Lahiri et al., New J. Phys. 13, 025001 (2011)](https://doi.org/10.1088/1367-2630/13/2/025001)\n", + "(open access), whose Table 1 is the quantitative target here:\n", + "\n", + "| interface | work of adhesion (J/m²) | separation (Å) |\n", + "|---|---|---|\n", + "| fcc (atop + fcc hollow) | 0.81 | 2.16 |\n", + "| hcp (atop + hcp hollow) | 0.77 | 2.17 |\n", + "| hollow (fcc + hcp hollows) | 0.31 | 3.26 |\n", + "\n", + "The four candidate registries, in the review's own Fig. 1:\n", + "\n", + "\"The\n", + "\n", + "The bridge registry (d) is not quantified in either paper — it is included here as an extra point\n", + "beyond the published set.\n", + "\n", + "The published calculation (Lahiri et al., section 2.2) used **LDA, spin-polarized, with geometry\n", + "relaxation** — five Ni layers with the bottom two fixed — because \"GGA does not provide an adequate\n", + "description of Ni–graphene bonding\". This notebook follows that recipe in two tiers:\n", + "\n", + "- **Fast (here, in minutes):** each registry relaxed with the\n", + " [MACE-MP](https://github.com/ACEsuit/mace) machine-learned force field (+D3), with the bottom\n", + " substrate layers fixed as in the paper, heights only; same-cell references give the work of\n", + " adhesion. MACE is PBE-trained and misses the paper's numbers on this interface: chemisorption\n", + " several times too weak, the separation short, the atop carbon buckled the wrong way. What it\n", + " delivers in minutes is the registry set, the two-branch (chemisorbed / dispersion-bound) energy\n", + " landscape and the starting geometries for the precise tier; its table prints beside the paper's\n", + " so the gap shows.\n", + "- **Precise (platform jobs):** the paper's functional — **LDA** (pz, ultrasoft), spin-polarized,\n", + " **fixed-cell relaxation**, no dispersion correction — for each registry plus the two same-cell\n", + " references the work of adhesion needs; the relaxed geometry is read back and compared too.\n", + "\n", + "**Prerequisite:** run\n", + "[optimization_interface_film_xy_position_graphene_nickel.ipynb](optimization_interface_film_xy_position_graphene_nickel.ipynb)\n", + "first — it creates and saves the base interface material this notebook loads.\n", + "\n", + "## 1. Prepare the Environment\n", + "### 1.1. Install Packages\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.calculators.mlff import get_mlff_install_profiles\n", + "from mat3ra.notebooks_utils.packages import install_packages\n", + "\n", + "await install_packages(get_mlff_install_profiles(\"mace\"))\n", + "\n", + "from mat3ra.notebooks_utils.pyodide.packages.patches import apply_all_patches\n", + "\n", + "apply_all_patches(\"mace\")" + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "### 1.2. Set Parameters\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "from datetime import datetime\n", + "from mat3ra.ide.compute import QueueName\n", + "\n", + "# 2. Auth and organization parameters\n", + "ORGANIZATION_NAME = None\n", + "\n", + "# 3. Material parameters\n", + "FOLDER = \"./uploads\"\n", + "BASE_MATERIAL_NAME = \"Graphene_Nickel_interface\" # created by the companion structure notebook\n", + "\n", + "# 4. MLFF parameters\n", + "MACE_MODEL_FAMILY = \"MACE-MP-0\"\n", + "MACE_MODEL = \"large\"\n", + "MACE_DISPERSION = True\n", + "MACE_DEFAULT_DTYPE = \"float64\"\n", + "MACE_DEVICE = \"cpu\"\n", + "\n", + "# 5. Separation scan, in Angstrom — brackets both published minima (2.16 and 3.26 A)\n", + "Z_SCAN_START = 1.8\n", + "Z_SCAN_STOP = 4.3\n", + "Z_SCAN_STEP = 0.25\n", + "CHEMISORBED_BELOW = 2.6 # boundary between the chemisorbed and dispersion-bound branches\n", + "\n", + "# 6. Relaxation — the paper's scheme; the buckling is one of the published numbers\n", + "FMAX = 0.02 # eV/A\n", + "FROZEN_SUBSTRATE_LAYERS = 2\n", + "\n", + "# 7. Workflow parameters\n", + "WORKFLOW_SEARCH_TERM = \"fixed_cell_relaxation.json\"\n", + "APPLICATION_NAME = \"espresso\"\n", + "MY_WORKFLOW_NAME = \"Fixed-cell Relaxation (Gr/Ni registry)\"\n", + "\n", + "# Method parameters — the published setup (Lahiri et al., section 2.2) where the platform can\n", + "# express it: LDA, spin-polarized, relaxed, no dispersion correction.\n", + "PSEUDOPOTENTIAL_TYPE = \"us\"\n", + "FUNCTIONAL = \"pz\"\n", + "MODEL_SUBTYPE = \"lda\"\n", + "ECUTWFC = 40 # GBRV's published pair\n", + "ECUTRHO = 200\n", + "SCF_KGRID = [12, 12, 1] # multiple of 3 keeps K on the mesh; dense for a metal\n", + "STARTING_MAGNETIZATION = {\"Ni\": 0.7} # near the bulk moment\n", + "\n", + "# SCF settings for a spin-polarized metal slab\n", + "SMEARING = \"mv\"\n", + "DEGAUSS = 0.01 # Ry\n", + "ADDITIONAL_PARAMETERS = {\n", + " \"electrons\": {\n", + " \"mixing_mode\": \"local-TF\",\n", + " \"mixing_beta\": 0.2,\n", + " \"electron_maxstep\": 200,\n", + " },\n", + "}\n", + "\n", + "# 8. Compute parameters — a spin-polarized relaxation on a 12x12x1 grid needs a full node\n", + "CLUSTER_NAME = None\n", + "QUEUE_NAME = QueueName.OF\n", + "PPN = 40\n", + "TIME_LIMIT = \"04:00:00\"\n", + "\n", + "# 9. Job parameters\n", + "timestamp = datetime.now().strftime(\"%Y-%m-%d %H:%M\")\n", + "POLL_INTERVAL = 30\n" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "## 2. Load the Base Interface\n", + "\n", + "The base interface is created by the companion structure notebook and saved into `uploads/`.\n", + "It is required — this notebook does not substitute another material.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.made.material import Material\n", + "from mat3ra.made.tools.analyze.other import get_average_interlayer_distance\n", + "from mat3ra.made.tools.convert import to_ase\n", + "from mat3ra.made.tools.convert.interface_parts_enum import InterfacePartsEnum\n", + "from mat3ra.made.tools.modify import interface_get_part\n", + "from mat3ra.notebooks_utils.material import load_material_from_folder\n", + "from mat3ra.notebooks_utils.ipython.entity.material.visualize import visualize_materials as visualize\n", + "\n", + "base_interface = load_material_from_folder(FOLDER, BASE_MATERIAL_NAME)\n", + "if base_interface is None:\n", + " raise RuntimeError(\n", + " f\"'{BASE_MATERIAL_NAME}' not found in {FOLDER} — run \"\n", + " \"optimization_interface_film_xy_position_graphene_nickel.ipynb first.\"\n", + " )\n", + "\n", + "film_part = interface_get_part(base_interface, part=InterfacePartsEnum.FILM)\n", + "substrate_part = interface_get_part(base_interface, part=InterfacePartsEnum.SUBSTRATE)\n", + "film_elements = set(film_part.basis.elements.values)\n", + "substrate_elements = set(substrate_part.basis.elements.values)\n", + "substrate_indices = [i for i, label in enumerate(base_interface.basis.labels.values)\n", + " if label == InterfacePartsEnum.SUBSTRATE.value]\n", + "n_carbon = len(film_part.basis.elements.values)\n", + "measured_gap = get_average_interlayer_distance(\n", + " to_ase(base_interface), InterfacePartsEnum.SUBSTRATE.value, InterfacePartsEnum.FILM.value)\n", + "\n", + "print(f\"Material: {base_interface.name}\")\n", + "print(f\"Atoms: {len(base_interface.basis.elements.values)} \"\n", + " f\"({n_carbon} film C, {len(substrate_indices)} substrate Ni)\")\n", + "print(f\"Film-substrate separation as built: {measured_gap:.3f} A\")\n", + "\n", + "visualize([{\"material\": base_interface, \"title\": base_interface.name}], repetitions=[3, 3, 1], rotation=\"-90x\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "## 3. Place the Film at the High-Symmetry Registries\n", + "\n", + "The registries are defined by where carbon atoms sit relative to the Ni(111) surface sites:\n", + "**top** (above a first-layer Ni), **hcp hollow** (above a second-layer Ni), **fcc hollow**\n", + "(above a third-layer Ni), and **bridge** — the C–C bond midpoint sits over a first-layer Ni\n", + "(Fig. 1d), so neither carbon lands on a named site.\n", + "The sites are measured from the structure itself — the top three Ni layers — and the film is\n", + "translated so one carbon sublattice lands on each site in turn.\n", + "\n", + "### 3.1. Measure the Registry Shifts\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from mat3ra.made.tools.analyze.other import get_closest_site_id_from_coordinate_and_element\n", + "from mat3ra.made.tools.helpers import SurfaceSiteAnalyzer, get_film_buckling, get_film_site_occupation\n", + "from mat3ra.made.tools.modify import interface_displace_part\n", + "\n", + "surface = SurfaceSiteAnalyzer(material=substrate_part)\n", + "film_z = float(np.mean([c[2] for c in film_part.coordinates_array])) # crystal (fractional) coordinate\n", + "anchor = get_closest_site_id_from_coordinate_and_element(film_part, [1 / 3, 2 / 3, film_z], \"C\")\n", + "film_cartesian = film_part.clone()\n", + "film_cartesian.to_cartesian()\n", + "carbon_xy = [np.array(c[:2]) for c in film_cartesian.coordinates_array]\n", + "\n", + "ANCHOR_SITE = {\"atop_fcc\": \"fcc\", \"atop_hcp\": \"atop\", \"hollow\": \"hcp\"}\n", + "displacements = {label: surface.get_displacement_to_site(carbon_xy[anchor], site, use_cartesian_coordinates=True) for label, site in ANCHOR_SITE.items()}\n", + "displacements[\"bridge\"] = surface.get_displacement_to_site(np.mean(carbon_xy, axis=0), \"atop\", use_cartesian_coordinates=True)\n", + "\n", + "REGISTRY_SITES = {\"atop_fcc\": {\"atop\", \"fcc\"}, \"atop_hcp\": {\"atop\", \"hcp\"}, \"hollow\": {\"fcc\", \"hcp\"}}\n", + "for label, shift in displacements.items():\n", + " occupied = get_film_site_occupation(interface_displace_part(base_interface, displacement=list(shift)), surface)\n", + " print(f\"{label:<10} carbons on {sorted(str(site) for site in occupied.values())} shift (A): {np.round(shift[:2], 3) + 0.0}\")\n", + " if label in REGISTRY_SITES:\n", + " assert set(occupied.values()) == REGISTRY_SITES[label], f\"{label}: carbons on {sorted(map(str, occupied.values()))}\"\n", + " else:\n", + " assert set(occupied.values()) == {None}, f\"{label}: carbons on {sorted(map(str, occupied.values()))}\"\n" + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "### 3.2. Preview Each Registry\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "preview = []\n", + "for label, shift in displacements.items():\n", + " m = interface_displace_part(base_interface, displacement=list(shift))\n", + " m.name = f\"{BASE_MATERIAL_NAME} {label}\"\n", + " preview.append({\"material\": m, \"title\": label})\n", + "\n", + "visualize(preview, repetitions=[2, 2, 1])" + ] + }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "## 4. Fast Tier: Relax Each Registry with MACE\n", + "\n", + "Each registry is bracketed by a rigid scan, then **relaxed** — positions move along z only, with\n", + "the deepest substrate layers fixed, for the interface and both same-cell references (bare Ni slab,\n", + "free-standing graphene) alike. For the paper's symmetric registries this equals full relaxation,\n", + "since in-plane forces vanish by symmetry; for the bridge, an in-plane saddle, it is what keeps the\n", + "point defined. Relaxing all three under the same constraint turns total energies into a work of\n", + "adhesion: W = (E_slab + E_graphene − E_interface) / A. Every atom's xy is held fixed by the\n", + "z-only constraint, so no film can slide into a neighbouring registry in this tier; the occupation\n", + "is still re-measured because the platform tier below relaxes every coordinate freely, where a\n", + "slide is possible, and there a structure that lands in a different registry is dropped rather\n", + "than reported under the wrong name. Distances follow the paper's convention: the averaged carbon\n", + "height above the averaged top-Ni height; buckling is the height difference between the two\n", + "carbons, positive when the atop carbon sits further out.\n", + "\n", + "### 4.1. Build the MACE Calculator\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "import importlib.util\n", + "\n", + "from mat3ra.notebooks_utils.calculators.mlff import create_mlff_calculator\n", + "\n", + "dispersion_available = importlib.util.find_spec(\"torch_dftd\") is not None\n", + "dispersion_active = MACE_DISPERSION and dispersion_available\n", + "if MACE_DISPERSION and not dispersion_available:\n", + " print(\"torch-dftd is not available here: the fast tier runs WITHOUT dispersion — the\")\n", + " print(\"GGA-level picture the manuscript describes as inadequate for this interface.\")\n", + "\n", + "calculator = create_mlff_calculator(\n", + " \"mace\",\n", + " {\n", + " \"family\": MACE_MODEL_FAMILY,\n", + " \"model\": MACE_MODEL,\n", + " \"dispersion\": dispersion_active,\n", + " \"default_dtype\": MACE_DEFAULT_DTYPE,\n", + " \"device\": MACE_DEVICE,\n", + " },\n", + ")\n" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "### 4.2. Relax the Same-Cell References\n", + "\n", + "The bare Ni slab and free-standing graphene, relaxed under the same z-only constraint as the\n", + "interface, are the two references the work of adhesion needs.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.made.tools.calculate import calculate_adhesion_energy, calculate_total_energy\n", + "from mat3ra.made.tools.helpers import get_atom_indices_by_layer\n", + "from mat3ra.notebooks_utils.workflows.relaxation import relax_material\n", + "\n", + "EV_PER_A2_TO_J_PER_M2 = 16.0217663\n", + "layers = get_atom_indices_by_layer(base_interface)\n", + "frozen = [i for layer in layers[:FROZEN_SUBSTRATE_LAYERS] for i in layer if i in substrate_indices]\n", + "\n", + "substrate_layers = get_atom_indices_by_layer(substrate_part)\n", + "slab_relaxed = relax_material(substrate_part, calculator, fmax=FMAX,\n", + " fixed_atom_indices=[i for layer in substrate_layers[:FROZEN_SUBSTRATE_LAYERS] for i in layer],\n", + " along_z_only=True)\n", + "film_relaxed = relax_material(film_part, calculator, fmax=FMAX, along_z_only=True)" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "### 4.3. Scan, Relax and Measure Each Registry\n", + "\n", + "For each registry: a rigid z-scan brackets the minimum on each branch (chemisorbed / dispersion-\n", + "bound), the bracketed geometry is relaxed, and the relaxed structure is checked against its\n", + "registry's expected sites before its work of adhesion, separation and buckling are recorded.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "distances = np.arange(Z_SCAN_START, Z_SCAN_STOP + 1e-9, Z_SCAN_STEP)\n", + "\n", + "scan_results = {}\n", + "for label, shift in displacements.items():\n", + " energies = np.array([\n", + " calculate_total_energy(\n", + " interface_displace_part(base_interface, displacement=list(shift + np.array([0.0, 0.0, d - measured_gap]))),\n", + " calculator,\n", + " )\n", + " for d in distances\n", + " ])\n", + "\n", + " # a bracketed minimum: the lowest scanned point of a branch that is not a window edge\n", + " starts = {}\n", + " for branch, in_branch in ((\"chem\", distances < CHEMISORBED_BELOW), (\"phys\", distances >= CHEMISORBED_BELOW)):\n", + " branch_idx = np.where(in_branch)[0]\n", + " j = int(np.argmin(energies[branch_idx]))\n", + " if 0 < j < len(branch_idx) - 1: # interior of the branch, so j is a real bracketed minimum\n", + " starts[branch] = float(distances[branch_idx[j]])\n", + " if not starts:\n", + " scan_results[label] = {\"energies\": energies, \"chem\": None, \"relaxed\": None}\n", + " print(f\"{label:<10} unbound in this window\" + (\"\" if dispersion_active else \" (dispersion inactive)\"))\n", + " continue\n", + "\n", + " start = starts.get(\"chem\", starts.get(\"phys\"))\n", + " displaced = interface_displace_part(base_interface, displacement=list(shift + np.array([0.0, 0.0, start - measured_gap])))\n", + " relaxed = relax_material(displaced, calculator, fmax=FMAX, fixed_atom_indices=frozen, along_z_only=True)\n", + "\n", + " occupied = set(get_film_site_occupation(relaxed).values())\n", + " buckling = get_film_buckling(relaxed)\n", + " if label in REGISTRY_SITES and occupied != REGISTRY_SITES[label]:\n", + " scan_results[label] = {\"energies\": energies, \"chem\": starts.get(\"chem\"), \"relaxed\": None}\n", + " print(f\"{label:<10} relaxed onto {occupied}: not a {label} result, dropped\")\n", + " continue\n", + " scan_results[label] = {\"energies\": energies, \"chem\": starts.get(\"chem\"), \"relaxed\": {\n", + " \"w_adh\": calculate_adhesion_energy(relaxed, slab_relaxed, film_relaxed, calculator) * EV_PER_A2_TO_J_PER_M2,\n", + " \"separation\": get_average_interlayer_distance(\n", + " to_ase(relaxed), InterfacePartsEnum.SUBSTRATE.value, InterfacePartsEnum.FILM.value),\n", + " \"buckling\": buckling,\n", + " \"material\": relaxed,\n", + " }}\n", + " r = scan_results[label][\"relaxed\"]\n", + " buckling_cell = \" — \" if buckling is None else f\"{buckling:+.3f}\"\n", + " print(f\"{label:<10} relaxed: d = {r['separation']:5.2f} A buckling = {buckling_cell} A \"\n", + " f\"W_adh = {r['w_adh']:.2f} J/m^2\")" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "### 4.4. Plot the Rigid Scans\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17", + "metadata": {}, + "outputs": [], + "source": [ + "import plotly.graph_objects as go\n", + "\n", + "reference = min(float(r[\"energies\"].min()) for r in scan_results.values())\n", + "fig = go.Figure()\n", + "for label, r in scan_results.items():\n", + " fig.add_trace(go.Scatter(x=distances, y=(r[\"energies\"] - reference) * 1000 / n_carbon,\n", + " mode=\"lines+markers\", name=label))\n", + "fig.update_layout(\n", + " title=\"Rigid-scan energy vs. separation (bracketing only; the table below is relaxed)\",\n", + " xaxis_title=\"plane distance (A)\",\n", + " yaxis_title=\"energy above the deepest scanned point (meV / C atom)\",\n", + ")\n", + "fig.show()\n" + ] + }, + { + "cell_type": "markdown", + "id": "18", + "metadata": {}, + "source": [ + "### 4.5. Compare the Fast Tier with the Paper\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19", + "metadata": {}, + "outputs": [], + "source": [ + "# Lahiri et al. (2011), Table 1\n", + "PAPER = {\"atop_fcc\": (0.81, 2.16), \"atop_hcp\": (0.77, 2.17), \"hollow\": (0.31, 3.26)}\n", + "PAPER_BUCKLING_FCC = 0.03 # A, atop_fcc, computed (Lahiri et al. ref 35)\n", + "\n", + "rows = {label: r[\"relaxed\"] for label, r in scan_results.items() if r[\"relaxed\"]}\n", + "print(f\"{'registry':<10}{'W_adh':>7}{'paper':>7} {'d':>5}{'paper':>7} buckling\")\n", + "for label, r in sorted(rows.items(), key=lambda kv: -kv[1][\"w_adh\"]):\n", + " w, d = PAPER.get(label, (\"—\", \"—\"))\n", + " buckling_cell = \" — \" if r[\"buckling\"] is None else f\"{r['buckling']:+.3f}\"\n", + " print(f\"{label:<10}{r['w_adh']:>7.2f}{w:>7} {r['separation']:>5.2f}{d:>7} {buckling_cell}\")\n", + "for label in set(scan_results) - set(rows):\n", + " w, d = PAPER.get(label, (\"—\", \"—\"))\n", + " print(f\"{label:<10}no result here — paper: {w} J/m^2 at {d} A\")" + ] + }, + { + "cell_type": "markdown", + "id": "20", + "metadata": {}, + "source": [ + "## 5. Precise Tier: the Paper's LDA, Relaxed, on the Platform\n", + "\n", + "One **fixed-cell relaxation** per selected registry, starting from the MACE-relaxed geometry, at the\n", + "paper's functional — LDA, spin-polarized, no dispersion correction — plus the two same-cell\n", + "references the work of adhesion needs. Each job's final structure is read back, so separation and\n", + "buckling are compared as well as the energy. The graphene reference runs without spin polarization:\n", + "it is non-magnetic, and a symmetric spin-polarized solution has the same energy.\n", + "\n", + "Divergences from Lahiri et al.: 4 Ni layers, not 5; 20 Å of vacuum, not 90; the platform relaxes\n", + "every atom, where the paper held the bottom two Ni layers; plane-wave ultrasoft pseudopotentials,\n", + "not all-electron LCAO.\n", + "\n", + "A default run selects one registry (three jobs). An **empty** list skips the platform tier entirely,\n", + "which is what the automated test does: relaxations take longer than a browser test may wait.\n", + "\n", + "### 5.1. Select the Registries for the Precise Tier\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21", + "metadata": {}, + "outputs": [], + "source": [ + "DFT_REGISTRY_NAMES = [\n", + " \"atop_fcc\",\n", + " # \"atop_hcp\",\n", + " # \"hollow\",\n", + " # \"bridge\",\n", + "]\n" + ] + }, + { + "cell_type": "markdown", + "id": "22", + "metadata": {}, + "source": [ + "### 5.2. Authenticate\n", + "\n", + "Authenticate in the browser and have credentials stored in environment variable\n", + "`OIDC_ACCESS_TOKEN`.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.auth import authenticate\n", + "\n", + "await authenticate()" + ] + }, + { + "cell_type": "markdown", + "id": "24", + "metadata": {}, + "source": [ + "### 5.3. Initialize the API Client\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "25", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.api_client import APIClient\n", + "\n", + "client = APIClient.authenticate()\n", + "client" + ] + }, + { + "cell_type": "markdown", + "id": "26", + "metadata": {}, + "source": [ + "### 5.4. Select the Account\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "27", + "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": "28", + "metadata": {}, + "source": [ + "### 5.5. Select the Project\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "29", + "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": "30", + "metadata": {}, + "source": [ + "### 5.6. Save the Registry Materials to the Platform\n", + "\n", + "Each selected registry's relaxed geometry, plus the same-cell substrate and film references, is\n", + "saved under the platform's own species labelling.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.core.entity.material.api import get_or_create_material\n", + "\n", + "dft_materials, reference_materials = {}, {}\n", + "if DFT_REGISTRY_NAMES:\n", + " for label in DFT_REGISTRY_NAMES:\n", + " relaxed = scan_results[label][\"relaxed\"]\n", + " if relaxed is None:\n", + " print(f\"{label:<16} skipped: no relaxed structure from the fast tier\")\n", + " continue\n", + " # QE species names must match between input blocks: drop the film/substrate labels.\n", + " m = relaxed[\"material\"].clone()\n", + " m.basis.set_labels_from_list(None)\n", + " m.name = f\"{BASE_MATERIAL_NAME} {label} d{relaxed['separation']:.2f} relaxed\"\n", + " saved = Material.create(get_or_create_material(client, m, ACCOUNT_ID))\n", + " dft_materials[label] = saved\n", + " print(f\"{label:<16} -> '{saved.name}' ({len(saved.basis.elements.values)} atoms)\")\n", + " for name, part in ((\"substrate\", substrate_part), (\"film\", film_part)) if dft_materials else ():\n", + " m = part.clone()\n", + " m.basis.set_labels_from_list(None)\n", + " m.name = f\"{BASE_MATERIAL_NAME} {name} reference\"\n", + " saved = Material.create(get_or_create_material(client, m, ACCOUNT_ID))\n", + " reference_materials[name] = saved\n", + " print(f\"{name + ' ref':<16} -> '{saved.name}' ({len(saved.basis.elements.values)} atoms)\")\n", + "else:\n", + " print(\"DFT tier skipped: no registries selected.\")" + ] + }, + { + "cell_type": "markdown", + "id": "32", + "metadata": {}, + "source": [ + "### 5.7. Select the Application\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "33", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.standata.applications import ApplicationStandata\n", + "from mat3ra.ade.application import Application\n", + "\n", + "app_config = ApplicationStandata.get_by_name_first_match(APPLICATION_NAME)\n", + "app = Application(**app_config)\n", + "print(f\"Using application: {app.name}\")" + ] + }, + { + "cell_type": "markdown", + "id": "34", + "metadata": {}, + "source": [ + "### 5.8. Create the Workflow and Preview It\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "35", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.standata.workflows import WorkflowStandata\n", + "from mat3ra.wode.workflows import Workflow\n", + "from mat3ra.notebooks_utils.ipython.entity.workflow.visualize import visualize_workflow\n", + "\n", + "workflow_config = WorkflowStandata.filter_by_application(app.name).get_by_name_first_match(WORKFLOW_SEARCH_TERM)\n", + "workflow = Workflow.create(workflow_config)\n", + "workflow.name = MY_WORKFLOW_NAME\n", + "\n", + "visualize_workflow(workflow)" + ] + }, + { + "cell_type": "markdown", + "id": "36", + "metadata": {}, + "source": [ + "### 5.9. Set the DFT Model\n", + "\n", + "The paper's functional: LDA, no dispersion correction.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "37", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.mode import ModelFactory\n", + "from mat3ra.standata.model_tree import ModelTreeStandata\n", + "\n", + "# The paper's functional: LDA, no dispersion correction.\n", + "model_config = ModelTreeStandata.get_model_by_parameters(\n", + " type=\"dft\",\n", + " subtype=MODEL_SUBTYPE,\n", + " functional=FUNCTIONAL,\n", + ")\n", + "model_config[\"method\"] = {\"type\": \"pseudopotential\", \"subtype\": PSEUDOPOTENTIAL_TYPE}\n", + "model = ModelFactory.create(model_config)\n" + ] + }, + { + "cell_type": "markdown", + "id": "38", + "metadata": {}, + "source": [ + "### 5.10. Configure the Workflow's Cutoffs, K-grid and Magnetization\n", + "\n", + "The published settings applied to the relaxation unit; a Ni moment only where there is Ni.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "39", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.workflow import apply_planewave_cutoffs, apply_scf_kgrid, patch_workflow_qe_input\n", + "\n", + "RELAX_UNIT = \"pw_relax\"\n", + "\n", + "if dft_materials:\n", + " reference_material = next(iter(dft_materials.values()))\n", + " if reference_material.basis.elements.values[0] != \"Ni\":\n", + " raise RuntimeError(\"Expected Ni as the first species — the magnetization index assumes it\")\n", + " for subworkflow in workflow.subworkflows:\n", + " subworkflow.model = model\n", + " apply_planewave_cutoffs(workflow, ECUTWFC, ECUTRHO, unit_name=RELAX_UNIT)\n", + " apply_scf_kgrid(workflow, SCF_KGRID, material=reference_material, unit_name=RELAX_UNIT)\n", + " system = {\n", + " \"degauss\": DEGAUSS,\n", + " \"smearing\": SMEARING,\n", + " \"nspin\": 2,\n", + " \"starting_magnetization(1)\": STARTING_MAGNETIZATION[\"Ni\"],\n", + " }\n", + " patch_workflow_qe_input(workflow, {\"system\": system, **ADDITIONAL_PARAMETERS}, unit_names=[RELAX_UNIT])" + ] + }, + { + "cell_type": "markdown", + "id": "40", + "metadata": {}, + "source": [ + "### 5.11. Configure and Save the Workflows\n", + "\n", + "The film reference needs its own workflow object — same settings, without the substrate's\n", + "magnetization — so each of the three (interface, substrate, film) is saved separately.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "41", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.core.entity.workflow.api import get_or_create_workflow\n", + "\n", + "saved_workflows = {}\n", + "if dft_materials:\n", + " film_workflow = Workflow.create(WorkflowStandata.filter_by_application(app.name)\n", + " .get_by_name_first_match(WORKFLOW_SEARCH_TERM))\n", + " film_workflow.name = f\"{MY_WORKFLOW_NAME} film\"\n", + " for subworkflow in film_workflow.subworkflows:\n", + " subworkflow.model = model\n", + " apply_planewave_cutoffs(film_workflow, ECUTWFC, ECUTRHO, unit_name=RELAX_UNIT)\n", + " apply_scf_kgrid(film_workflow, SCF_KGRID, material=reference_material, unit_name=RELAX_UNIT)\n", + " system = {\"degauss\": DEGAUSS, \"smearing\": SMEARING, \"nspin\": 1}\n", + " patch_workflow_qe_input(film_workflow, {\"system\": system, **ADDITIONAL_PARAMETERS}, unit_names=[RELAX_UNIT])\n", + "\n", + " workflows = {\"interface\": workflow, \"substrate\": workflow, \"film\": film_workflow}\n", + " for key, built in workflows.items():\n", + " saved_workflows[key] = Workflow.create(get_or_create_workflow(client, built, ACCOUNT_ID))\n", + " print(f\"{key:<12} -> workflow {saved_workflows[key].id}\")" + ] + }, + { + "cell_type": "markdown", + "id": "42", + "metadata": {}, + "source": [ + "### 5.12. List the Available Clusters\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43", + "metadata": {}, + "outputs": [], + "source": [ + "if dft_materials:\n", + " print(f\"Available clusters: {[c['hostname'] for c in client.clusters.list()]}\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "44", + "metadata": {}, + "source": [ + "### 5.13. Create the Compute Configuration\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "45", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.ide.compute import Compute\n", + "\n", + "compute = None\n", + "if dft_materials:\n", + " clusters = client.clusters.list()\n", + " matching = [c for c in clusters if not CLUSTER_NAME or CLUSTER_NAME in c[\"hostname\"]]\n", + " if not matching:\n", + " raise RuntimeError(f\"No cluster matching {CLUSTER_NAME!r} is available; registered: \"\n", + " f\"{[c['hostname'] for c in clusters]}\")\n", + " compute = Compute(cluster=matching[0], queue=QUEUE_NAME, ppn=PPN, timeLimit=TIME_LIMIT)\n", + " print(f\"Using cluster: {compute.cluster.hostname}, queue: {QUEUE_NAME}, ppn: {PPN}, time limit: {TIME_LIMIT}\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "46", + "metadata": {}, + "source": [ + "### 5.14. Create One Job per Registry\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "47", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.utils.namespace import dict_to_namespace_recursive\n", + "from mat3ra.notebooks_utils.job import create_job\n", + "\n", + "jobs, reference_jobs = {}, {}\n", + "if dft_materials:\n", + " for label, m in dft_materials.items():\n", + " job_response = create_job(\n", + " api_client=client,\n", + " materials=[m],\n", + " workflow=workflows[\"interface\"],\n", + " project_id=project_id,\n", + " owner_id=ACCOUNT_ID,\n", + " prefix=f\"{MY_WORKFLOW_NAME} {label} {timestamp}\",\n", + " compute=compute.to_dict(),\n", + " )\n", + " jobs[label] = dict_to_namespace_recursive(job_response)._id\n", + " print(f\"{label:<16} -> job {jobs[label]}\")\n", + " for name, m in reference_materials.items():\n", + " job_response = create_job(\n", + " api_client=client,\n", + " materials=[m],\n", + " workflow=workflows[name],\n", + " project_id=project_id,\n", + " owner_id=ACCOUNT_ID,\n", + " prefix=f\"{MY_WORKFLOW_NAME} {name} reference {timestamp}\",\n", + " compute=compute.to_dict(),\n", + " )\n", + " reference_jobs[name] = dict_to_namespace_recursive(job_response)._id\n", + " print(f\"{name + ' reference':<16} -> job {reference_jobs[name]}\")" + ] + }, + { + "cell_type": "markdown", + "id": "48", + "metadata": {}, + "source": [ + "### 5.15. Submit the Jobs\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "49", + "metadata": {}, + "outputs": [], + "source": [ + "for label, job_id in {**jobs, **reference_jobs}.items():\n", + " client.jobs.submit(job_id)\n", + " print(f\"Submitted {label}: {job_id}\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "50", + "metadata": {}, + "source": [ + "### 5.16. Wait for the Jobs to Finish\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "51", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.api.job import wait_for_jobs_to_finish_async\n", + "\n", + "all_job_ids = list(jobs.values()) + list(reference_jobs.values())\n", + "if all_job_ids:\n", + " await wait_for_jobs_to_finish_async(client.jobs, all_job_ids, poll_interval=POLL_INTERVAL)\n", + "else:\n", + " print(\"Nothing to wait for — the DFT tier was skipped.\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "52", + "metadata": {}, + "source": [ + "### 5.17. Retrieve and Compare the DFT Results\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "53", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.made.tools.analyze.other import get_surface_area\n", + "from mat3ra.made.tools.modify import interface_label_parts_by_elements\n", + "from mat3ra.notebooks_utils.core.entity.material.api import get_final_structure_for_job\n", + "from mat3ra.notebooks_utils.core.entity.property.api import get_properties_for_job\n", + "\n", + "area = get_surface_area(to_ase(base_interface))\n", + "\n", + "dft_results = {}\n", + "if jobs:\n", + " reference_energies = {name: get_properties_for_job(client, job_id, \"total_energy\")[-1][\"value\"]\n", + " for name, job_id in reference_jobs.items()}\n", + " separated = reference_energies[\"substrate\"] + reference_energies[\"film\"]\n", + " for label, job_id in jobs.items():\n", + " energy = get_properties_for_job(client, job_id, \"total_energy\")[-1][\"value\"]\n", + " material = interface_label_parts_by_elements(get_final_structure_for_job(client, job_id), substrate_elements)\n", + " separation = get_average_interlayer_distance(\n", + " to_ase(material), InterfacePartsEnum.SUBSTRATE.value, InterfacePartsEnum.FILM.value)\n", + "\n", + " occupied = set(get_film_site_occupation(material).values())\n", + " buckling = get_film_buckling(material)\n", + "\n", + " drifted = label in REGISTRY_SITES and occupied != REGISTRY_SITES[label]\n", + " if drifted:\n", + " print(f\"! {label}: relaxed onto {sorted(str(s) for s in occupied)}, not {sorted(REGISTRY_SITES[label])} — kept below, marked\")\n", + " dft_results[label] = {\"energy\": energy, \"w_adh\": (separated - energy) / area * EV_PER_A2_TO_J_PER_M2,\n", + " \"separation\": separation, \"buckling\": buckling, \"drifted\": drifted, \"sites\": occupied}\n", + "\n", + " dft_cells = {}\n", + " for label, r in dft_results.items():\n", + " if r[\"drifted\"]:\n", + " sites = \"/\".join(sorted(str(s) for s in r[\"sites\"]))\n", + " dft_cells[label] = f\"{label}→{sites}\"\n", + " else:\n", + " dft_cells[label] = label\n", + " dft_width = max(len(\"registry\"), *(len(c) for c in dft_cells.values()))\n", + " print(f\"{'registry':<{dft_width}}{'E (eV)':>12}{'W_adh':>8}{'paper':>7} {'d':>5}{'paper':>7} buckling\")\n", + " for label, r in sorted(dft_results.items(), key=lambda kv: -kv[1][\"w_adh\"]):\n", + " w, d = PAPER.get(label, (\"—\", \"—\"))\n", + " buckling_cell = \" — \" if r[\"buckling\"] is None else f\"{r['buckling']:+.3f}\"\n", + " print(f\"{dft_cells[label]:<{dft_width}}{r['energy']:>12.4f}{r['w_adh']:>8.2f}{w:>7} {r['separation']:>5.2f}{d:>7} {buckling_cell}\")" + ] + }, + { + "cell_type": "markdown", + "id": "54", + "metadata": {}, + "source": [ + "## 6. Compare with the Article\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "55", + "metadata": {}, + "outputs": [], + "source": [ + "# Lahiri et al. (2011), Table 1, beside what each tier here computes.\n", + "LABELS = (\"atop_fcc\", \"atop_hcp\", \"hollow\", \"bridge\")\n", + "\n", + "final_cells = {}\n", + "for label in LABELS:\n", + " dft = dft_results.get(label)\n", + " if dft and dft[\"drifted\"]:\n", + " sites = \"/\".join(sorted(str(s) for s in dft[\"sites\"]))\n", + " final_cells[label] = f\"{label}→{sites}\"\n", + " else:\n", + " final_cells[label] = label\n", + "\n", + "final_width = max(len(\"registry\"), *(len(c) for c in final_cells.values()))\n", + "print(f\"{'registry':<{final_width}}{'W_adh':>7}{'MACE':>7}{'DFT':>7} {'d':>5}{'MACE':>7}{'DFT':>7} \"\n", + " f\"{'buckling':>8}{'MACE':>9}{'DFT':>9}\")\n", + "for label in LABELS:\n", + " w, d = PAPER.get(label, (None, None))\n", + " b_paper = PAPER_BUCKLING_FCC if label == \"atop_fcc\" else None\n", + " mace, dft = rows.get(label), dft_results.get(label)\n", + "\n", + " w_cell = f\"{w:.2f}\" if w is not None else \"—\"\n", + " d_cell = f\"{d:.2f}\" if d is not None else \"—\"\n", + " b_paper_cell = f\"{b_paper:.2f}\" if b_paper is not None else \"—\"\n", + " mace_w_adh_cell = f\"{mace['w_adh']:.2f}\" if mace else \"—\"\n", + " dft_w_adh_cell = f\"{dft['w_adh']:.2f}\" if dft else \"—\"\n", + " mace_separation_cell = f\"{mace['separation']:.2f}\" if mace else \"—\"\n", + " dft_separation_cell = f\"{dft['separation']:.2f}\" if dft else \"—\"\n", + " mace_buckling = mace[\"buckling\"] if mace else None\n", + " dft_buckling = dft[\"buckling\"] if dft else None\n", + " mace_buckling_cell = \" — \" if mace_buckling is None else f\"{mace_buckling:+.3f}\"\n", + " dft_buckling_cell = \" — \" if dft_buckling is None else f\"{dft_buckling:+.3f}\"\n", + "\n", + " print(f\"{final_cells[label]:<{final_width}}\"\n", + " f\"{w_cell:>7}{mace_w_adh_cell:>7}{dft_w_adh_cell:>7} \"\n", + " f\"{d_cell:>5}{mace_separation_cell:>7}{dft_separation_cell:>7} \"\n", + " f\"{b_paper_cell:>8}{mace_buckling_cell:>9}{dft_buckling_cell:>9}\")" + ] + }, + { + "cell_type": "markdown", + "id": "56", + "metadata": {}, + "source": [ + "## References\n", + "\n", + "[1] Arjun Dahal, Matthias Batzill, \"Graphene-nickel interfaces: a review\",\n", + "Nanoscale 6(5), 2548 (2014). [DOI: 10.1039/c3nr05279f](https://doi.org/10.1039/c3nr05279f)\n", + "\n", + "[2] Jayeeta Lahiri, Travis S. Miller, Andrew J. Ross, Lyudmyla Adamska, Ivan I. Oleynik,\n", + "Matthias Batzill, \"Graphene growth and stability at nickel surfaces\", New J. Phys. 13, 025001\n", + "(2011). [DOI: 10.1088/1367-2630/13/2/025001](https://doi.org/10.1088/1367-2630/13/2/025001)\n", + "\n", + "[3] mat3ra-made: https://github.com/mat3ra/made\n", + "\n", + "[4] MACE-MP-0 foundation models: https://github.com/ACEsuit/mace\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/other/materials_designer/workflows/local/reaction_path_afir_mace.ipynb b/other/materials_designer/workflows/local/reaction_path_afir_mace.ipynb index fc5adbf66..3b517c0ba 100644 --- a/other/materials_designer/workflows/local/reaction_path_afir_mace.ipynb +++ b/other/materials_designer/workflows/local/reaction_path_afir_mace.ipynb @@ -54,7 +54,7 @@ "metadata": {}, "outputs": [], "source": [ - "from mat3ra.notebooks_utils.mlff import get_mlff_install_profiles\n", + "from mat3ra.notebooks_utils.calculators.mlff import get_mlff_install_profiles\n", "from mat3ra.notebooks_utils.packages import install_packages\n", "\n", "await install_packages(get_mlff_install_profiles(\"mace\"))\n", @@ -302,7 +302,7 @@ "metadata": {}, "outputs": [], "source": [ - "from mat3ra.notebooks_utils.mlff import create_mlff_calculator\n", + "from mat3ra.notebooks_utils.calculators.mlff import create_mlff_calculator\n", "\n", "MACE_MODEL_LABEL = f\"{MACE_MODEL_FAMILY} ({MACE_MODEL})\"\n", "\n", diff --git a/other/materials_designer/workflows/local/relaxation_mlff_mace.ipynb b/other/materials_designer/workflows/local/relaxation_mlff_mace.ipynb index 0e3865ade..7de30d2b3 100644 --- a/other/materials_designer/workflows/local/relaxation_mlff_mace.ipynb +++ b/other/materials_designer/workflows/local/relaxation_mlff_mace.ipynb @@ -40,7 +40,7 @@ "metadata": {}, "outputs": [], "source": [ - "from mat3ra.notebooks_utils.mlff import get_mlff_install_profiles\n", + "from mat3ra.notebooks_utils.calculators.mlff import get_mlff_install_profiles\n", "from mat3ra.notebooks_utils.packages import install_packages\n", "\n", "await install_packages(get_mlff_install_profiles(\"mace\"))\n", @@ -244,7 +244,7 @@ "outputs": [], "source": [ "from mat3ra.made.tools.convert import to_ase\n", - "from mat3ra.notebooks_utils.mlff import create_mlff_calculator\n", + "from mat3ra.notebooks_utils.calculators.mlff import create_mlff_calculator\n", "\n", "calculator = create_mlff_calculator(\n", " \"mace\",\n", diff --git a/src/py/mat3ra/notebooks_utils/calculators/__init__.py b/src/py/mat3ra/notebooks_utils/calculators/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/src/py/mat3ra/notebooks_utils/mlff.py b/src/py/mat3ra/notebooks_utils/calculators/mlff.py similarity index 100% rename from src/py/mat3ra/notebooks_utils/mlff.py rename to src/py/mat3ra/notebooks_utils/calculators/mlff.py diff --git a/src/py/mat3ra/notebooks_utils/core/entity/material/api.py b/src/py/mat3ra/notebooks_utils/core/entity/material/api.py index 287beb0ce..9891585e9 100644 --- a/src/py/mat3ra/notebooks_utils/core/entity/material/api.py +++ b/src/py/mat3ra/notebooks_utils/core/entity/material/api.py @@ -3,7 +3,9 @@ from mat3ra.api_client import APIClient from mat3ra.made.material import Material +from mat3ra.prode import PropertyName +from ..property.api import get_properties_for_job from .analysis import get_slab_bulk_crystal, resolve_bulk_query_from_crystal ORDERED_ENTITY_SET_TYPE = "ordered" @@ -32,6 +34,14 @@ def get_or_create_material(api_client: APIClient, material, owner_id: str) -> di return created +def get_final_structure_for_job(api_client: APIClient, job_id: str) -> Material: + """Fetch the relaxed structure a job reported as its `final_structure` property.""" + properties = get_properties_for_job(api_client, job_id, PropertyName.non_scalar.final_structure.value) + if not properties: + raise RuntimeError(f"Job {job_id} reported no 'final_structure'") + return Material.create(api_client.materials.get(properties[-1]["materialId"])) + + def get_bulk_material(api_client: APIClient, slab_material: Material, owner_id: str) -> Material: """ Resolves the platform bulk material a slab was built from, owned by the given account. diff --git a/src/py/mat3ra/notebooks_utils/workflow.py b/src/py/mat3ra/notebooks_utils/workflow.py index e5de259fe..ef43ca99f 100644 --- a/src/py/mat3ra/notebooks_utils/workflow.py +++ b/src/py/mat3ra/notebooks_utils/workflow.py @@ -2,7 +2,7 @@ from typing import Dict, List, Optional from mat3ra.wode import Workflow -from mat3ra.wode.context.providers import PointsGridDataProvider +from mat3ra.wode.context.providers import PlanewaveCutoffsContextProvider, PointsGridDataProvider FORTRAN_NUMBER_PATTERN = re.compile(r"^[+-]?(?:\d+(?:\.\d*)?|\.\d+)(?:[de][+-]?\d+)?$", re.IGNORECASE) @@ -88,3 +88,17 @@ def apply_scf_kgrid( if first_only: break return workflow + + +def apply_planewave_cutoffs(workflow: Workflow, wavefunction, density, *, unit_name: str = "pw_relax") -> Workflow: + """Attaches an edited planewave cutoffs context to units named `unit_name`.""" + context = PlanewaveCutoffsContextProvider( + wavefunction=wavefunction, density=density, isEdited=True + ).get_context_item_data() + for subworkflow in workflow.subworkflows: + if unit_name not in [unit.name for unit in subworkflow.units]: + continue + unit = subworkflow.get_unit_by_name(name=unit_name) + unit.add_context(context) + subworkflow.set_unit(unit) + return workflow diff --git a/src/py/mat3ra/notebooks_utils/workflows/__init__.py b/src/py/mat3ra/notebooks_utils/workflows/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/src/py/mat3ra/notebooks_utils/workflows/relaxation.py b/src/py/mat3ra/notebooks_utils/workflows/relaxation.py new file mode 100644 index 000000000..e8bdd4f0c --- /dev/null +++ b/src/py/mat3ra/notebooks_utils/workflows/relaxation.py @@ -0,0 +1,56 @@ +from typing import Optional, Sequence + +from ase.constraints import FixAtoms, FixedLine +from ase.optimize import BFGS +from mat3ra.made.material import Material +from mat3ra.made.tools.convert import to_ase +from mat3ra.made.tools.third_party import ASECalculator + +Z_DIRECTION = [0, 0, 1] + + +def relax_material( + material: Material, + calculator: ASECalculator, + fmax: float = 0.05, + max_steps: int = 300, + fixed_atom_indices: Optional[Sequence[int]] = None, + along_z_only: bool = False, + logfile: Optional[str] = "-", +) -> Material: + """ + Relax atomic positions with an ASE calculator (e.g. from `create_mlff_calculator`) at fixed + cell, optionally holding atoms fixed or allowing motion along z only. + + Args: + material: The structure to relax; labels and build metadata are preserved in the result. + calculator: Any ASE calculator. + fmax: Force convergence criterion, eV/Angstrom. + max_steps: Optimizer step limit. + fixed_atom_indices: Atoms held fixed. + along_z_only: Restrict every atom's motion to the z direction. + logfile: ASE optimizer log target; "-" is stdout, None silences it. + + Raises: + RuntimeError: when the optimizer stops before the forces fall below `fmax`. + """ + atoms = to_ase(material) + constraints = [] + if fixed_atom_indices: + constraints.append(FixAtoms(indices=list(fixed_atom_indices))) + if along_z_only: + constraints.append(FixedLine(list(range(len(atoms))), direction=Z_DIRECTION)) + if constraints: + atoms.set_constraint(constraints) + atoms.calc = calculator + converged = BFGS(atoms, logfile=logfile).run(fmax=fmax, steps=max_steps) + if not converged: + raise RuntimeError(f"Relaxation of '{material.name}' did not reach fmax={fmax} eV/A within {max_steps} steps.") + + relaxed = material.clone() + was_in_crystal_units = relaxed.basis.is_in_crystal_units + relaxed.to_cartesian() + relaxed.set_coordinates(atoms.positions.tolist()) + if was_in_crystal_units: + relaxed.to_crystal() + return relaxed diff --git a/tests/py/unit/core/entity/test_material_api.py b/tests/py/unit/core/entity/test_material_api.py index b1cfbebb6..1802094d2 100644 --- a/tests/py/unit/core/entity/test_material_api.py +++ b/tests/py/unit/core/entity/test_material_api.py @@ -5,10 +5,12 @@ import pytest from mat3ra.notebooks_utils.core.entity.material.api import ( find_material_set, + get_final_structure_for_job, get_or_create_materials_set, list_materials_by_set, list_materials_in_set, ) +from mat3ra.standata.materials import Materials OWNER_ID = "account-1" MATERIAL_SET_NAME = "H2+H" @@ -224,3 +226,28 @@ def test_get_or_create_materials_set_requires_one_material(): is_ordered=False, ) client.materials.list.assert_not_called() + + +JOB_ID = "job-1" +FINAL_STRUCTURE_MATERIAL_ID = "m-final-structure" + + +@pytest.mark.parametrize( + ("properties", "error"), + [ + ([{"materialId": FINAL_STRUCTURE_MATERIAL_ID}], None), + ([], "reported no 'final_structure'"), + ], +) +def test_get_final_structure_for_job(properties, error): + client = MagicMock() + client.properties.get_for_job.return_value = properties + client.materials.get.return_value = Materials.get_by_name_first_match("Silicon") + if error: + with pytest.raises(RuntimeError, match=error): + get_final_structure_for_job(client, JOB_ID) + return + material = get_final_structure_for_job(client, JOB_ID) + assert material.basis.elements.values == ["Si", "Si"] + client.properties.get_for_job.assert_called_once_with(JOB_ID, "final_structure") + client.materials.get.assert_called_once_with(FINAL_STRUCTURE_MATERIAL_ID) diff --git a/tests/py/unit/test_workflow_utils.py b/tests/py/unit/test_workflow_utils.py index 00a488c57..2a472fd40 100644 --- a/tests/py/unit/test_workflow_utils.py +++ b/tests/py/unit/test_workflow_utils.py @@ -1,6 +1,6 @@ import pytest from mat3ra.made.material import Material -from mat3ra.notebooks_utils.workflow import apply_scf_kgrid, patch_workflow_qe_input +from mat3ra.notebooks_utils.workflow import apply_planewave_cutoffs, apply_scf_kgrid, patch_workflow_qe_input from mat3ra.standata.workflows import WorkflowStandata from mat3ra.wode.workflows import Workflow @@ -88,3 +88,13 @@ def test_apply_scf_kgrid_updates_pw_scf_context(): # KPPRA is per reciprocal atom, and the ratios come from the lattice -- both via `material`. assert kgrid_item["data"]["gridMetricValue"] == 4 * 4 * 1 * 2 assert kgrid_item["data"]["reciprocalVectorRatios"] == [1.0, 1.0, 0.5] + + +@pytest.mark.parametrize("wavefunction,density", [(40, 200), (60, 480)]) +def test_apply_planewave_cutoffs_updates_pw_relax_context(wavefunction, density): + workflow = _relax_workflow() + apply_planewave_cutoffs(workflow, wavefunction, density, unit_name="pw_relax") + unit = workflow.subworkflows[0].get_unit_by_name(name="pw_relax") + cutoffs_item = next(item for item in unit.context if item.get("name") == "cutoffs") + assert cutoffs_item["data"]["wavefunction"] == float(wavefunction) + assert cutoffs_item["data"]["density"] == float(density) diff --git a/tests/py/unit/workflows/__init__.py b/tests/py/unit/workflows/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/py/unit/workflows/test_relaxation.py b/tests/py/unit/workflows/test_relaxation.py new file mode 100644 index 000000000..c643efbd5 --- /dev/null +++ b/tests/py/unit/workflows/test_relaxation.py @@ -0,0 +1,80 @@ +import numpy as np +import pytest +from ase.calculators.emt import EMT +from mat3ra.made.material import Material +from mat3ra.made.tools.calculate import calculate_total_energy +from mat3ra.made.tools.helpers import create_slab +from mat3ra.notebooks_utils.workflows.relaxation import relax_material +from mat3ra.standata.materials import Materials + +# A plain slab, not an interface: relax_material's contract is about constraints (fixed atoms, +# along_z_only, non-convergence), not about Gr/Ni physics, and the interface path is already +# covered end to end by other/materials_designer/specific_examples/ +# optimization_interface_film_xy_position_graphene_nickel_SIMULATION.ipynb and by made's own tests. +# The top atom is displaced deliberately (not left at its built position) so every case tests the +# contract against a real, comfortable force margin rather than however close standata's Ni +# lattice constant happens to sit to EMT's own equilibrium. +LAYER_TOLERANCE = 0.5 # Angstrom: heights closer than this belong to the same layer + +MATERIAL = create_slab( + crystal=Material.create(Materials.get_by_name_first_match("Nickel")), + miller_indices=(1, 0, 0), + number_of_layers=4, + vacuum=10.0, +) + +_cartesian = MATERIAL.clone() +_cartesian.to_cartesian() +_z = [c[2] for c in _cartesian.coordinates_array] +BOTTOM_LAYER = [i for i, z in enumerate(_z) if z - min(_z) < LAYER_TOLERANCE] +TOP_ATOM = max(range(len(_z)), key=lambda i: _z[i]) + + +def _displaced(axis: int, amount: float) -> Material: + displaced = _cartesian.clone() + coordinates = displaced.coordinates_array + coordinates[TOP_ATOM][axis] += amount + displaced.set_coordinates(coordinates) + displaced.to_crystal() + return displaced + + +Z_DISPLACED = _displaced(2, 0.3) # out-of-plane: a real force for the along_z_only case +XY_DISPLACED = _displaced(0, 0.3) # in-plane: a real force to hold still or let drift back + +CALCULATOR = EMT() +RELAX = {"fmax": 0.1, "max_steps": 50, "logfile": None} + +CASES = [ + # (material, fixed_atom_indices, along_z_only, xy_unchanged) + (Z_DISPLACED, BOTTOM_LAYER, True, True), + (XY_DISPLACED, BOTTOM_LAYER, False, False), # in-plane force free to act: the atom drifts back + (XY_DISPLACED, BOTTOM_LAYER, True, True), # same force, held to z: the atom cannot drift +] + + +def _cartesian_positions(material: Material) -> np.ndarray: + cartesian = material.clone() + cartesian.to_cartesian() + return np.array(cartesian.coordinates_array) + + +@pytest.mark.parametrize("material, fixed_atom_indices, along_z_only, xy_unchanged", CASES) +# measured (pytest --durations=0): module 4.34s total, slowest case 0.05s +def test_relax_material(material, fixed_atom_indices, along_z_only, xy_unchanged): + relaxed = relax_material( + material, CALCULATOR, fixed_atom_indices=fixed_atom_indices, along_z_only=along_z_only, **RELAX + ) + assert calculate_total_energy(relaxed, CALCULATOR) < calculate_total_energy(material, CALCULATOR) + assert relaxed.name == material.name + assert relaxed.basis.labels.values == material.basis.labels.values + assert relaxed.basis.is_in_crystal_units == material.basis.is_in_crystal_units + + before, after = _cartesian_positions(material), _cartesian_positions(relaxed) + assert np.allclose(after[fixed_atom_indices], before[fixed_atom_indices]) + assert np.allclose(after[:, :2], before[:, :2], atol=1e-6) == xy_unchanged + + +def test_relax_material_invalid(): + with pytest.raises(RuntimeError): + relax_material(MATERIAL, CALCULATOR, fmax=1e-6, max_steps=1, logfile=None)