From 65862d33d40032cba5742a08346ae810306a138f Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Tue, 11 Aug 2026 19:53:54 -0700 Subject: [PATCH 1/5] feat: add AFIR NB first implementation (claude) --- .../local/reaction_path_afir_mace.ipynb | 746 ++++++++++++++++++ 1 file changed, 746 insertions(+) create mode 100644 other/materials_designer/workflows/local/reaction_path_afir_mace.ipynb diff --git a/other/materials_designer/workflows/local/reaction_path_afir_mace.ipynb b/other/materials_designer/workflows/local/reaction_path_afir_mace.ipynb new file mode 100644 index 00000000..b2bec2fa --- /dev/null +++ b/other/materials_designer/workflows/local/reaction_path_afir_mace.ipynb @@ -0,0 +1,746 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Discover a reaction path with AFIR + MLFF (MACE)\n", + "\n", + "Find how a molecule rearranges **without knowing the product in advance**: an artificial force (AFIR) pulls two selected atoms together on top of the physical potential energy surface, while a **MACE** foundation model supplies energies and forces, so bonds break and form on the fly.\n", + "\n", + "The example is the **Claisen rearrangement** of allyl vinyl ether (C5H8O) into 4-pentenal — a concerted [3,3]-sigmatropic shift in which a C–O bond breaks and a C–C bond forms in one step. Its ~30 kcal/mol barrier is far too high for a plain relaxation to cross, which is exactly the situation AFIR is built for.\n", + "\n", + "

Usage

\n", + "\n", + "1. Set the molecule, the reacting atom pairs and the MACE model in cells 1.2 and 1.3 (or use default values).\n", + "1. Click \"Run\" > \"Run All\" to run all cells.\n", + "1. Wait for the search to finish (a few minutes on a laptop CPU).\n", + "1. Scroll down to view the discovered path, the transition state and the energy diagram.\n", + "\n", + "## Summary\n", + "\n", + "Load the molecule, relax it with MACE, ramp an artificial force between the two target atoms until the reaction happens, strip the bias to recover the physical energy profile, relax the discovered product, refine the highest point of the path into a true saddle point, verify it with a vibrational analysis and by relaxing along the imaginary mode in both directions, and save the reactant, transition state and product together with the energy profile and the plots.\n", + "\n", + "Everything runs locally through ASE — no platform compute and no authentication." + ] + }, + { + "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|torch|mace\")\n", + "\n", + "from mat3ra.notebooks_utils.pyodide.packages.patches import apply_all_patches\n", + "\n", + "apply_all_patches(\"mace\")" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "### 1.2. Set the reaction\n", + "\n", + "The atom pairs are indices into the loaded structure. Cell 2.1 prints the element of every index so they can be checked before the search starts." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "FOLDER = \"../../uploads\" # molecules are read from here, and the resulting materials are written back here\n", + "RESULTS_FOLDER = \"results\" # energy profile, plots and the raw trajectory\n", + "MOLECULE_NAME = \"allyl vinyl ether\" # looked up in the uploads folder first, then fetched from PubChem by name\n", + "PRODUCT_NAME = \"4-pentenal\"\n", + "\n", + "# The bond that AFIR forces to form: the two terminal CH2 carbons\n", + "BOND_FORMING_PAIR = (4, 5)\n", + "# The bond expected to break in response: the ether oxygen and the allylic CH2\n", + "BOND_BREAKING_PAIR = (0, 1)\n", + "# The bond that becomes the product carbonyl\n", + "CARBONYL_PAIR = (0, 3)" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "### 1.3. AFIR and MACE options\n", + "\n", + "AFIR adds a bias term $E_\\text{bias} = \\alpha \\, r_{ij}$ between the target atoms, i.e. a constant attractive force $\\alpha$ that pulls them together no matter how the physical surface resists. A single value of $\\alpha$ either leaves the molecule stuck in a biased minimum (too weak) or drags it through a badly distorted geometry (too strong), so the force is ramped: each stage starts from the structure the previous one converged to." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "# Artificial force applied between the target atoms, ramped stage by stage (eV/Å)\n", + "AFIR_FORCE_RAMP = [1.0, 2.0, 3.0, 4.0]\n", + "AFIR_MAX_STEPS_PER_STAGE = 150\n", + "AFIR_MAX_DISPLACEMENT = 0.1 # per-step displacement cap, keeps the biased path smooth (Å)\n", + "AFIR_TRAJECTORY_PATH = f\"{RESULTS_FOLDER}/afir_path.traj\"\n", + "\n", + "# Maximum force on any atom at convergence (eV/Å)\n", + "RELAXATION_FMAX = 0.03\n", + "AFIR_FMAX = 0.05\n", + "SADDLE_FMAX = 0.02\n", + "\n", + "# Modes below this magnitude are the translations and rotations of a free molecule (cm⁻¹)\n", + "IMAGINARY_MODE_THRESHOLD = 50\n", + "# Displacement along the imaginary mode used to leave the saddle point (Å)\n", + "REACTION_MODE_DISPLACEMENT = 0.3\n", + "\n", + "MACE_MODEL_FAMILY = \"off\" # \"off\": MACE-OFF23, trained on organic molecules; \"mp\": MACE-MP-0, trained on inorganic crystals\n", + "MACE_MODEL_SIZE = \"medium\" # choose between \"small\", \"medium\" and \"large\"\n", + "MACE_DEFAULT_DTYPE = \"float64\" # float64 is recommended for geometry optimization and vibrations\n", + "MACE_DEVICE = \"cpu\" # hardware target: \"cpu\" or \"cuda\" (GPU)" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "## 2. Load the molecule\n", + "### 2.1. Read from uploads, or fetch the 3D structure from PubChem\n", + "\n", + "PubChem serves an optimized 3D conformer for most small molecules, which makes any named molecule a valid starting point for the search." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "import io\n", + "import os\n", + "from urllib.parse import quote\n", + "\n", + "from ase.io import read, write\n", + "\n", + "PUBCHEM_STRUCTURE_URL = \"https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/name/{name}/SDF?record_type=3d\"\n", + "\n", + "\n", + "def fetch_pubchem_structure(name):\n", + " url = PUBCHEM_STRUCTURE_URL.format(name=quote(name))\n", + " try:\n", + " from pyodide.http import open_url\n", + "\n", + " return open_url(url).read()\n", + " except ImportError:\n", + " from urllib.request import urlopen\n", + "\n", + " return urlopen(url).read().decode()\n", + "\n", + "\n", + "molecule_path = os.path.join(FOLDER, MOLECULE_NAME.replace(\" \", \"_\") + \".xyz\")\n", + "\n", + "if not os.path.exists(molecule_path):\n", + " write(molecule_path, read(io.StringIO(fetch_pubchem_structure(MOLECULE_NAME)), format=\"sdf\"))\n", + " print(f\"Fetched {MOLECULE_NAME} from PubChem, saved to {molecule_path}\")\n", + "\n", + "molecule = read(molecule_path)\n", + "\n", + "print(f\"{MOLECULE_NAME}: {molecule.get_chemical_formula()}\")\n", + "print(\"atoms: \" + \", \".join(f\"{index}:{symbol}\" for index, symbol in enumerate(molecule.get_chemical_symbols())))\n", + "print(f\"bond to form {BOND_FORMING_PAIR}: {molecule.get_distance(*BOND_FORMING_PAIR):.2f} Å\")\n", + "print(f\"bond to break {BOND_BREAKING_PAIR}: {molecule.get_distance(*BOND_BREAKING_PAIR):.2f} Å\")" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "### 2.2. View the molecule" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.made.material import Material\n", + "from mat3ra.made.tools.convert import from_ase\n", + "from mat3ra.notebooks_utils.ipython.entity.material.visualize import ViewersEnum, visualize_materials as visualize\n", + "\n", + "VACUUM = 5.0 # padding around the molecule, in Å, so it can be handled as a material\n", + "\n", + "\n", + "def to_material(atoms, name):\n", + " boxed_atoms = atoms.copy()\n", + " boxed_atoms.center(vacuum=VACUUM)\n", + " material = Material.create(from_ase(boxed_atoms))\n", + " material.name = name\n", + " return material\n", + "\n", + "\n", + "visualize([{\"material\": to_material(molecule, MOLECULE_NAME), \"title\": MOLECULE_NAME}], viewer=ViewersEnum.wave)" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "## 3. Relax the reactant with MACE\n", + "### 3.1. Create the ASE calculator\n", + "\n", + "MACE-OFF23 is trained on organic molecules and is the appropriate foundation model here; it is distributed under the Academic Software License, which does not permit commercial use. MACE-MP-0 is trained on inorganic crystals from Materials Project trajectories — it is selectable through `MACE_MODEL_FAMILY` for comparison, but it does not break the C–O bond for this reaction." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "from mace.calculators import mace_mp, mace_off\n", + "\n", + "mace_foundation_model = {\"off\": mace_off, \"mp\": mace_mp}[MACE_MODEL_FAMILY]\n", + "calculator = mace_foundation_model(\n", + " model=MACE_MODEL_SIZE,\n", + " default_dtype=MACE_DEFAULT_DTYPE,\n", + " device=MACE_DEVICE,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "13", + "metadata": {}, + "source": [ + "### 3.2. Relax the reactant" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ + "from ase.optimize import BFGS\n", + "\n", + "reactant = molecule.copy()\n", + "reactant.calc = calculator\n", + "\n", + "BFGS(reactant, logfile=None).run(fmax=RELAXATION_FMAX)\n", + "reactant_energy = reactant.get_potential_energy()\n", + "\n", + "print(f\"Relaxed {MOLECULE_NAME}: {reactant_energy:.3f} eV\")" + ] + }, + { + "cell_type": "markdown", + "id": "15", + "metadata": {}, + "source": [ + "## 4. Push the reaction with an artificial force\n", + "\n", + "The bias is applied with ASE's `ExternalForce` constraint, which adds exactly the AFIR term: a constant force of $\\alpha$ along the vector connecting the two target atoms, with an energy contribution $\\alpha \\, r_{ij}$. The physical forces still come from MACE, so the molecule is free to respond by breaking whatever bond is in the way." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16", + "metadata": {}, + "outputs": [], + "source": [ + "from ase.constraints import ExternalForce\n", + "from ase.io.trajectory import Trajectory\n", + "\n", + "structure = reactant.copy()\n", + "structure.calc = calculator\n", + "\n", + "os.makedirs(RESULTS_FOLDER, exist_ok=True)\n", + "trajectory = Trajectory(AFIR_TRAJECTORY_PATH, \"w\", structure)\n", + "trajectory.write()\n", + "\n", + "for artificial_force in AFIR_FORCE_RAMP:\n", + " structure.set_constraint(ExternalForce(*BOND_FORMING_PAIR, -artificial_force))\n", + " BFGS(structure, trajectory=trajectory, maxstep=AFIR_MAX_DISPLACEMENT, logfile=None).run(\n", + " fmax=AFIR_FMAX, steps=AFIR_MAX_STEPS_PER_STAGE\n", + " )\n", + " print(\n", + " f\"α = {artificial_force:.1f} eV/Å → \"\n", + " f\"d{BOND_FORMING_PAIR} = {structure.get_distance(*BOND_FORMING_PAIR):.2f} Å, \"\n", + " f\"d{BOND_BREAKING_PAIR} = {structure.get_distance(*BOND_BREAKING_PAIR):.2f} Å\"\n", + " )\n", + "\n", + "structure.set_constraint()" + ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, + "source": [ + "## 5. Recover the physical energy landscape\n", + "### 5.1. Strip the bias\n", + "\n", + "Every structure along the biased path is re-evaluated with the bias removed, which turns the AFIR trajectory into a physical energy profile. Its maximum is the first estimate of the transition state." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "EV_TO_KCAL_PER_MOL = 23.060548\n", + "\n", + "images = read(AFIR_TRAJECTORY_PATH, index=\":\")\n", + "unbiased_energies = []\n", + "for image in images:\n", + " image.set_constraint()\n", + " image.calc = calculator\n", + " unbiased_energies.append(image.get_potential_energy())\n", + "\n", + "path_energies = (np.array(unbiased_energies) - reactant_energy) * EV_TO_KCAL_PER_MOL\n", + "transition_state_guess_index = int(np.argmax(path_energies))\n", + "\n", + "print(f\"AFIR path: {len(images)} structures\")\n", + "print(\n", + " f\"Highest point at step {transition_state_guess_index}: \"\n", + " f\"{path_energies[transition_state_guess_index]:.1f} kcal/mol above the reactant\"\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "19", + "metadata": {}, + "source": [ + "### 5.2. Plot the discovered path" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20", + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib import pyplot as plt\n", + "from mat3ra.notebooks_utils.plot import display_matplotlib_figure\n", + "\n", + "forming_distances = [image.get_distance(*BOND_FORMING_PAIR) for image in images]\n", + "breaking_distances = [image.get_distance(*BOND_BREAKING_PAIR) for image in images]\n", + "\n", + "path_figure, (energy_axes, distance_axes) = plt.subplots(2, 1, figsize=(8, 7), sharex=True)\n", + "\n", + "energy_axes.plot(path_energies, color=\"#2b5c8f\", linewidth=2)\n", + "energy_axes.axvline(transition_state_guess_index, color=\"#c93b3b\", linestyle=\"--\", label=\"Transition state guess\")\n", + "energy_axes.set_ylabel(\"Energy relative to reactant (kcal/mol)\")\n", + "energy_axes.set_title(f\"AFIR path: {MOLECULE_NAME} (MACE-{MACE_MODEL_FAMILY.upper()}, {MACE_MODEL_SIZE})\")\n", + "energy_axes.legend()\n", + "energy_axes.grid(True, linestyle=\":\", alpha=0.6)\n", + "\n", + "distance_axes.plot(forming_distances, color=\"#2e8b57\", label=f\"forming {BOND_FORMING_PAIR}\")\n", + "distance_axes.plot(breaking_distances, color=\"#c93b3b\", label=f\"breaking {BOND_BREAKING_PAIR}\")\n", + "distance_axes.axvline(transition_state_guess_index, color=\"#c93b3b\", linestyle=\"--\")\n", + "distance_axes.set_xlabel(\"AFIR step\")\n", + "distance_axes.set_ylabel(\"Distance (Å)\")\n", + "distance_axes.legend()\n", + "distance_axes.grid(True, linestyle=\":\", alpha=0.6)\n", + "\n", + "path_figure.tight_layout()\n", + "display_matplotlib_figure(path_figure)" + ] + }, + { + "cell_type": "markdown", + "id": "21", + "metadata": {}, + "source": [ + "## 6. Relax the discovered product\n", + "\n", + "The last structure of the biased path is relaxed with the bias removed, which lets it settle into the product basin the search reached." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22", + "metadata": {}, + "outputs": [], + "source": [ + "product = images[-1].copy()\n", + "product.calc = calculator\n", + "\n", + "BFGS(product, logfile=None).run(fmax=RELAXATION_FMAX)\n", + "reaction_energy = (product.get_potential_energy() - reactant_energy) * EV_TO_KCAL_PER_MOL\n", + "\n", + "print(f\"Product energy: {reaction_energy:.1f} kcal/mol relative to the reactant\\n\")\n", + "print(f\"{'bond':<20}{'reactant':>12}{'product':>12}\")\n", + "for label, pair in (\n", + " (\"forming\", BOND_FORMING_PAIR),\n", + " (\"breaking\", BOND_BREAKING_PAIR),\n", + " (\"carbonyl\", CARBONYL_PAIR),\n", + "):\n", + " print(f\"{label + ' ' + str(pair):<20}{reactant.get_distance(*pair):>10.2f} Å{product.get_distance(*pair):>10.2f} Å\")" + ] + }, + { + "cell_type": "markdown", + "id": "23", + "metadata": {}, + "source": [ + "## 7. Refine the transition state\n", + "\n", + "The maximum of the AFIR path is a point on a biased trajectory, not a stationary point of the physical surface: the forces there are still large and its energy overestimates the barrier. The dimer method follows the lowest-curvature mode uphill and all remaining modes downhill, converging on the first-order saddle point nearest that guess — that saddle is what sets the activation energy. It is started along the reaction direction: the forming pair closing, the breaking pair opening." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24", + "metadata": {}, + "outputs": [], + "source": [ + "from ase.mep import DimerControl, MinModeAtoms, MinModeTranslate\n", + "\n", + "transition_state = images[transition_state_guess_index].copy()\n", + "transition_state.calc = calculator\n", + "print(f\"Maximum force at the AFIR guess: {np.abs(transition_state.get_forces()).max():.2f} eV/Å\")\n", + "\n", + "reaction_direction = np.zeros_like(transition_state.positions)\n", + "for pair, sign in ((BOND_FORMING_PAIR, 1.0), (BOND_BREAKING_PAIR, -1.0)):\n", + " unit_vector = transition_state.positions[pair[1]] - transition_state.positions[pair[0]]\n", + " unit_vector /= np.linalg.norm(unit_vector)\n", + " reaction_direction[pair[0]] += sign * unit_vector\n", + " reaction_direction[pair[1]] -= sign * unit_vector\n", + "reaction_direction /= np.linalg.norm(reaction_direction)\n", + "\n", + "dimer_control = DimerControl(\n", + " initial_eigenmode_method=\"displacement\",\n", + " displacement_method=\"vector\",\n", + " logfile=None,\n", + ")\n", + "dimer = MinModeAtoms(transition_state, dimer_control)\n", + "dimer.displace(displacement_vector=0.05 * reaction_direction, mask=[True] * len(transition_state))\n", + "\n", + "MinModeTranslate(dimer, logfile=None).run(fmax=SADDLE_FMAX, steps=200)\n", + "\n", + "transition_state_energy = transition_state.get_potential_energy()\n", + "activation_energy = (transition_state_energy - reactant_energy) * EV_TO_KCAL_PER_MOL\n", + "\n", + "print(f\"Maximum force at the saddle point: {np.abs(transition_state.get_forces()).max():.3f} eV/Å\")\n", + "print(f\"Activation energy: {activation_energy:.1f} kcal/mol\")\n", + "print(\n", + " f\"d{BOND_FORMING_PAIR} = {transition_state.get_distance(*BOND_FORMING_PAIR):.2f} Å, \"\n", + " f\"d{BOND_BREAKING_PAIR} = {transition_state.get_distance(*BOND_BREAKING_PAIR):.2f} Å\"\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "25", + "metadata": {}, + "source": [ + "## 8. Verify the transition state\n", + "\n", + "A first-order saddle point has exactly one imaginary vibrational frequency, and its mode is the reaction coordinate. A free molecule also has six translational and rotational modes at (numerically) near-zero frequency, which appear as small imaginary values — `IMAGINARY_MODE_THRESHOLD` separates those from a genuine reaction mode." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "26", + "metadata": {}, + "outputs": [], + "source": [ + "from ase.vibrations import Vibrations\n", + "\n", + "vibrations = Vibrations(transition_state, name=\"transition_state_vibrations\")\n", + "vibrations.run()\n", + "vibrations.summary()\n", + "\n", + "frequencies = vibrations.get_frequencies()\n", + "imaginary_mode_indices = [\n", + " index\n", + " for index, frequency in enumerate(frequencies)\n", + " if np.iscomplex(frequency) and abs(frequency.imag) > IMAGINARY_MODE_THRESHOLD\n", + "]\n", + "\n", + "print(f\"\\nImaginary modes above {IMAGINARY_MODE_THRESHOLD} cm⁻¹: {len(imaginary_mode_indices)}\")\n", + "for index in imaginary_mode_indices:\n", + " print(f\" mode {index}: {abs(frequencies[index].imag):.0f}i cm⁻¹\")\n", + "\n", + "reaction_mode = vibrations.get_mode(imaginary_mode_indices[0])\n", + "vibrations.clean()" + ] + }, + { + "cell_type": "markdown", + "id": "27", + "metadata": {}, + "source": [ + "## 9. Confirm which minima the saddle connects\n", + "\n", + "Displacing along the imaginary mode in both directions and relaxing shows what the transition state actually joins: one side must fall back to the reactant, the other into the product." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "28", + "metadata": {}, + "outputs": [], + "source": [ + "connected_minima = {}\n", + "\n", + "print(f\"{'direction':<12}{'energy, kcal/mol':>18}{'d' + str(BOND_FORMING_PAIR):>14}{'d' + str(BOND_BREAKING_PAIR):>14}\")\n", + "for sign, label in ((1.0, \"forward\"), (-1.0, \"reverse\")):\n", + " displaced = transition_state.copy()\n", + " displaced.positions += sign * REACTION_MODE_DISPLACEMENT * reaction_mode / np.linalg.norm(reaction_mode)\n", + " displaced.calc = calculator\n", + " BFGS(displaced, logfile=None).run(fmax=RELAXATION_FMAX, steps=400)\n", + " connected_minima[label] = displaced\n", + " energy = (displaced.get_potential_energy() - reactant_energy) * EV_TO_KCAL_PER_MOL\n", + " print(\n", + " f\"{label:<12}{energy:>18.1f}\"\n", + " f\"{displaced.get_distance(*BOND_FORMING_PAIR):>12.2f} Å\"\n", + " f\"{displaced.get_distance(*BOND_BREAKING_PAIR):>12.2f} Å\"\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "29", + "metadata": {}, + "source": [ + "## 10. Results\n", + "### 10.1. Energy diagram\n", + "\n", + "The experimental activation energy for this reaction is 30.6 kcal/mol in the gas phase, and the rearrangement is strongly exothermic. MACE reproduces the reaction energy well; the barrier is overestimated, as foundation models trained on near-equilibrium structures generally are in the bond-breaking region." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "30", + "metadata": {}, + "outputs": [], + "source": [ + "EXPERIMENTAL_ACTIVATION_ENERGY = 30.6 # kcal/mol, gas phase\n", + "\n", + "levels = [\n", + " (MOLECULE_NAME, 0.0),\n", + " (\"transition state\", activation_energy),\n", + " (PRODUCT_NAME, reaction_energy),\n", + "]\n", + "\n", + "diagram_figure, axes = plt.subplots(figsize=(7, 4.5))\n", + "axes.plot(range(len(levels)), [energy for _, energy in levels], linestyle=\"--\", color=\"#999999\")\n", + "for position, (label, energy) in enumerate(levels):\n", + " axes.hlines(energy, position - 0.25, position + 0.25, color=\"#2b5c8f\", linewidth=4)\n", + " axes.annotate(f\"{energy:.1f}\", (position, energy), textcoords=\"offset points\", xytext=(0, 10), ha=\"center\")\n", + "axes.axhline(EXPERIMENTAL_ACTIVATION_ENERGY, color=\"#c93b3b\", linestyle=\":\", label=\"experimental barrier\")\n", + "axes.set_xticks(range(len(levels)))\n", + "axes.set_xticklabels([label for label, _ in levels])\n", + "axes.set_ylabel(\"Energy relative to reactant (kcal/mol)\")\n", + "axes.set_title(f\"{MOLECULE_NAME} → {PRODUCT_NAME} (MACE-{MACE_MODEL_FAMILY.upper()}, {MACE_MODEL_SIZE})\")\n", + "axes.legend()\n", + "axes.grid(True, axis=\"y\", linestyle=\":\", alpha=0.6)\n", + "diagram_figure.tight_layout()\n", + "display_matplotlib_figure(diagram_figure)\n", + "\n", + "print(f\"Activation energy: {activation_energy:.1f} kcal/mol (experiment: {EXPERIMENTAL_ACTIVATION_ENERGY})\")\n", + "print(f\"Reaction energy: {reaction_energy:.1f} kcal/mol\")" + ] + }, + { + "cell_type": "markdown", + "id": "31", + "metadata": {}, + "source": [ + "### 10.2. View the reactant, transition state and product" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "32", + "metadata": {}, + "outputs": [], + "source": [ + "visualize(\n", + " [\n", + " {\"material\": to_material(reactant, MOLECULE_NAME), \"title\": MOLECULE_NAME},\n", + " {\"material\": to_material(transition_state, \"Transition state\"), \"title\": \"Transition state\"},\n", + " {\"material\": to_material(product, PRODUCT_NAME), \"title\": PRODUCT_NAME},\n", + " ],\n", + " viewer=ViewersEnum.wave,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "33", + "metadata": {}, + "source": [ + "## 11. Save the results\n", + "### 11.1. Hand the structures back as materials\n", + "\n", + "The three structures that define the reaction are passed to the environment with `set_materials`, each carrying the search settings and the resulting energies as metadata." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "34", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.material import set_materials\n", + "\n", + "reaction_metadata = {\n", + " \"type\": \"reaction_path\",\n", + " \"method\": {\n", + " \"name\": \"AFIR\",\n", + " \"artificial_force_ramp\": {\"values\": AFIR_FORCE_RAMP, \"units\": \"eV/angstrom\"},\n", + " \"bond_forming_pair\": list(BOND_FORMING_PAIR),\n", + " \"bond_breaking_pair\": list(BOND_BREAKING_PAIR),\n", + " },\n", + " \"engine\": {\"name\": \"ASE\", \"optimizer\": \"BFGS\", \"saddle_search\": \"Dimer\"},\n", + " \"calculator\": {\n", + " \"name\": f\"MACE-{MACE_MODEL_FAMILY.upper()}\",\n", + " \"parameters\": {\"model\": MACE_MODEL_SIZE, \"default_dtype\": MACE_DEFAULT_DTYPE},\n", + " },\n", + " \"reactant\": MOLECULE_NAME,\n", + " \"product\": PRODUCT_NAME,\n", + " \"activation_energy\": {\"value\": round(float(activation_energy), 2), \"units\": \"kcal/mol\"},\n", + " \"reaction_energy\": {\"value\": round(float(reaction_energy), 2), \"units\": \"kcal/mol\"},\n", + " \"imaginary_frequency\": {\n", + " \"value\": round(float(abs(frequencies[imaginary_mode_indices[0]].imag)), 1),\n", + " \"units\": \"cm-1\",\n", + " },\n", + "}\n", + "\n", + "structures = (\n", + " (\"reactant\", reactant, f\"{MOLECULE_NAME}, reactant\"),\n", + " (\"transition_state\", transition_state, f\"{MOLECULE_NAME} to {PRODUCT_NAME}, transition state\"),\n", + " (\"product\", product, f\"{PRODUCT_NAME}, product\"),\n", + ")\n", + "\n", + "materials = []\n", + "for role, atoms, name in structures:\n", + " material = to_material(atoms, name)\n", + " material.metadata.reaction_path = {**reaction_metadata, \"role\": role}\n", + " materials.append(material)\n", + "\n", + "set_materials(materials, FOLDER)" + ] + }, + { + "cell_type": "markdown", + "id": "35", + "metadata": {}, + "source": [ + "### 11.2. Write the energy profile and the plots\n", + "\n", + "The trajectory of the search is already in `RESULTS_FOLDER`; this adds the profile behind the plots, the key energies and the two figures, so the run can be reported without re-running it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "36", + "metadata": {}, + "outputs": [], + "source": [ + "import json\n", + "\n", + "results = {\n", + " \"reaction_path\": reaction_metadata,\n", + " \"distances_angstrom\": {\n", + " role: {\n", + " \"forming\": round(float(atoms.get_distance(*BOND_FORMING_PAIR)), 3),\n", + " \"breaking\": round(float(atoms.get_distance(*BOND_BREAKING_PAIR)), 3),\n", + " \"carbonyl\": round(float(atoms.get_distance(*CARBONYL_PAIR)), 3),\n", + " }\n", + " for role, atoms, _ in structures\n", + " },\n", + " \"afir_path\": {\n", + " \"transition_state_guess_index\": transition_state_guess_index,\n", + " \"energy_kcal_per_mol\": [round(float(value), 4) for value in path_energies],\n", + " \"forming_distance_angstrom\": [round(float(value), 3) for value in forming_distances],\n", + " \"breaking_distance_angstrom\": [round(float(value), 3) for value in breaking_distances],\n", + " },\n", + "}\n", + "\n", + "with open(os.path.join(RESULTS_FOLDER, \"reaction_path.json\"), \"w\") as file:\n", + " json.dump(results, file, indent=2)\n", + "\n", + "path_figure.savefig(os.path.join(RESULTS_FOLDER, \"afir_path.png\"), dpi=140)\n", + "diagram_figure.savefig(os.path.join(RESULTS_FOLDER, \"energy_diagram.png\"), dpi=140)\n", + "\n", + "print(f\"Saved to {RESULTS_FOLDER}/: \" + \", \".join(sorted(os.listdir(RESULTS_FOLDER))))" + ] + }, + { + "cell_type": "markdown", + "id": "37", + "metadata": {}, + "source": [ + "## References\n", + "\n", + "[1] AFIR method: S. Maeda, K. Morokuma, \"Communications: A systematic method for locating transition structures of A+B → X type reactions\", J. Chem. Phys. 132, 241102 (2010). https://doi.org/10.1063/1.3457903 \n", + "[2] AFIR review: S. Maeda, K. Ohno, K. Morokuma, \"Systematic exploration of the mechanism of chemical reactions: the global reaction route mapping (GRRM) strategy\", Phys. Chem. Chem. Phys. 15, 3683 (2013). https://doi.org/10.1039/C3CP44063J \n", + "[3] MACE-OFF23 organic foundation models: D. P. Kovács et al., arXiv:2312.15211. https://arxiv.org/abs/2312.15211 \n", + "[4] MACE-MP-0 materials foundation model: I. Batatia et al., arXiv:2401.00096. https://arxiv.org/abs/2401.00096 \n", + "[5] Dimer method: G. Henkelman, H. Jónsson, \"A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives\", J. Chem. Phys. 111, 7010 (1999). https://doi.org/10.1063/1.480097 \n", + "[6] Experimental Claisen barrier: F. W. Schuler, G. W. Murphy, \"The Kinetics of the Rearrangement of Vinyl Allyl Ether\", J. Am. Chem. Soc. 72, 3155 (1950). https://doi.org/10.1021/ja01163a096 \n", + "[7] ASE optimizers, constraints and vibrations: https://wiki.fysik.dtu.dk/ase/ase/optimize.html \n", + "[8] PubChem compound \"allyl vinyl ether\" (CID 221523): https://pubchem.ncbi.nlm.nih.gov/compound/221523 " + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 3b088bc8ad6515ae9284b6a5834b36d97f6395ba Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Tue, 11 Aug 2026 19:54:05 -0700 Subject: [PATCH 2/5] update: intro nb --- other/materials_designer/workflows/Introduction.ipynb | 3 +++ 1 file changed, 3 insertions(+) diff --git a/other/materials_designer/workflows/Introduction.ipynb b/other/materials_designer/workflows/Introduction.ipynb index a97425ba..51bef958 100644 --- a/other/materials_designer/workflows/Introduction.ipynb +++ b/other/materials_designer/workflows/Introduction.ipynb @@ -89,6 +89,9 @@ "### 7.3. Vibrational Frequency (NWChem)\n", "#### [7.3.1. Vibrational frequency calculation.](homo_lumo_frequency.ipynb)\n", "\n", + "### 7.4. Reaction Path Discovery with MLFF\n", + "#### [7.4.1. AFIR reaction path discovery (MACE).](local/reaction_path_afir_mace.ipynb)\n", + "\n", "\n", "## 8. Electronics\n", "\n", From bc9ea7748b7b18660ff2d33641c71b967df76d05 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Wed, 12 Aug 2026 10:46:23 -0700 Subject: [PATCH 3/5] update: save data to files --- .../local/reaction_path_afir_mace.ipynb | 141 ++++++++++-------- 1 file changed, 81 insertions(+), 60 deletions(-) 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 b2bec2fa..b8914ee7 100644 --- a/other/materials_designer/workflows/local/reaction_path_afir_mace.ipynb +++ b/other/materials_designer/workflows/local/reaction_path_afir_mace.ipynb @@ -41,9 +41,10 @@ "metadata": {}, "outputs": [], "source": [ + "from mat3ra.notebooks_utils.mlff import get_mlff_install_profiles\n", "from mat3ra.notebooks_utils.packages import install_packages\n", "\n", - "await install_packages(\"made|api_examples|torch|mace\")\n", + "await install_packages(get_mlff_install_profiles(\"mace\"))\n", "\n", "from mat3ra.notebooks_utils.pyodide.packages.patches import apply_all_patches\n", "\n", @@ -112,9 +113,14 @@ "IMAGINARY_MODE_THRESHOLD = 50\n", "# Displacement along the imaginary mode used to leave the saddle point (Å)\n", "REACTION_MODE_DISPLACEMENT = 0.3\n", - "\n", - "MACE_MODEL_FAMILY = \"off\" # \"off\": MACE-OFF23, trained on organic molecules; \"mp\": MACE-MP-0, trained on inorganic crystals\n", - "MACE_MODEL_SIZE = \"medium\" # choose between \"small\", \"medium\" and \"large\"\n", + "# Two relaxed structures count as different minima if a target distance differs by more than this (Å)\n", + "MINIMUM_SEPARATION = 0.5\n", + "\n", + "# \"mp\": MACE-MP models shipped with the platform, the only option that works in JupyterLite;\n", + "# \"off\": MACE-OFF23 for organic molecules, downloaded on first use, so local runs only\n", + "MACE_MODEL_FAMILY = \"mp\"\n", + "MACE_MODEL = \"medium\" # choose between \"small\", \"medium\" and \"large\"\n", + "MACE_DISPERSION = False # D3 dispersion correction, not needed for this intramolecular rearrangement\n", "MACE_DEFAULT_DTYPE = \"float64\" # float64 is recommended for geometry optimization and vibrations\n", "MACE_DEVICE = \"cpu\" # hardware target: \"cpu\" or \"cuda\" (GPU)" ] @@ -187,7 +193,7 @@ "metadata": {}, "outputs": [], "source": [ - "from mat3ra.made.material import Material\n", + "from mat3ra.made.tools.build_components import MaterialWithBuildMetadata\n", "from mat3ra.made.tools.convert import from_ase\n", "from mat3ra.notebooks_utils.ipython.entity.material.visualize import ViewersEnum, visualize_materials as visualize\n", "\n", @@ -197,7 +203,7 @@ "def to_material(atoms, name):\n", " boxed_atoms = atoms.copy()\n", " boxed_atoms.center(vacuum=VACUUM)\n", - " material = Material.create(from_ase(boxed_atoms))\n", + " material = MaterialWithBuildMetadata.create(from_ase(boxed_atoms))\n", " material.name = name\n", " return material\n", "\n", @@ -213,7 +219,9 @@ "## 3. Relax the reactant with MACE\n", "### 3.1. Create the ASE calculator\n", "\n", - "MACE-OFF23 is trained on organic molecules and is the appropriate foundation model here; it is distributed under the Academic Software License, which does not permit commercial use. MACE-MP-0 is trained on inorganic crystals from Materials Project trajectories — it is selectable through `MACE_MODEL_FAMILY` for comparison, but it does not break the C–O bond for this reaction." + "With `MACE_MODEL_FAMILY = \"mp\"` the calculator is built from the MACE-MP models shipped with the platform (`packages/models`), so nothing is downloaded at run time and the notebook runs in JupyterLite.\n", + "\n", + "Those models are trained on inorganic crystal trajectories. For this organic rearrangement they close the C–C bond but do not break the C–O bond, so the search ends at a cyclic structure rather than at 4-pentenal. `MACE_MODEL_FAMILY = \"off\"` selects MACE-OFF23, which is trained on organic molecules and reproduces the published mechanism — it is downloaded on first use, so it works locally but not in the browser, and its Academic Software License does not permit commercial use." ] }, { @@ -223,14 +231,22 @@ "metadata": {}, "outputs": [], "source": [ - "from mace.calculators import mace_mp, mace_off\n", - "\n", - "mace_foundation_model = {\"off\": mace_off, \"mp\": mace_mp}[MACE_MODEL_FAMILY]\n", - "calculator = mace_foundation_model(\n", - " model=MACE_MODEL_SIZE,\n", - " default_dtype=MACE_DEFAULT_DTYPE,\n", - " device=MACE_DEVICE,\n", - ")" + "from mat3ra.notebooks_utils.mlff import create_mlff_calculator\n", + "\n", + "if MACE_MODEL_FAMILY == \"off\":\n", + " from mace.calculators import mace_off\n", + "\n", + " calculator = mace_off(model=MACE_MODEL, default_dtype=MACE_DEFAULT_DTYPE, device=MACE_DEVICE)\n", + "else:\n", + " calculator = create_mlff_calculator(\n", + " \"mace\",\n", + " {\n", + " \"model\": MACE_MODEL,\n", + " \"dispersion\": MACE_DISPERSION,\n", + " \"default_dtype\": MACE_DEFAULT_DTYPE,\n", + " \"device\": MACE_DEVICE,\n", + " },\n", + " )" ] }, { @@ -250,10 +266,13 @@ "source": [ "from ase.optimize import BFGS\n", "\n", + "# ASE sends optimizer logs to /dev/null when no logfile is given, and Pyodide cannot flush it; a buffer works in both\n", + "OPTIMIZER_LOG = io.StringIO()\n", + "\n", "reactant = molecule.copy()\n", "reactant.calc = calculator\n", "\n", - "BFGS(reactant, logfile=None).run(fmax=RELAXATION_FMAX)\n", + "BFGS(reactant, logfile=OPTIMIZER_LOG).run(fmax=RELAXATION_FMAX)\n", "reactant_energy = reactant.get_potential_energy()\n", "\n", "print(f\"Relaxed {MOLECULE_NAME}: {reactant_energy:.3f} eV\")" @@ -288,7 +307,7 @@ "\n", "for artificial_force in AFIR_FORCE_RAMP:\n", " structure.set_constraint(ExternalForce(*BOND_FORMING_PAIR, -artificial_force))\n", - " BFGS(structure, trajectory=trajectory, maxstep=AFIR_MAX_DISPLACEMENT, logfile=None).run(\n", + " BFGS(structure, trajectory=trajectory, maxstep=AFIR_MAX_DISPLACEMENT, logfile=OPTIMIZER_LOG).run(\n", " fmax=AFIR_FMAX, steps=AFIR_MAX_STEPS_PER_STAGE\n", " )\n", " print(\n", @@ -365,7 +384,7 @@ "energy_axes.plot(path_energies, color=\"#2b5c8f\", linewidth=2)\n", "energy_axes.axvline(transition_state_guess_index, color=\"#c93b3b\", linestyle=\"--\", label=\"Transition state guess\")\n", "energy_axes.set_ylabel(\"Energy relative to reactant (kcal/mol)\")\n", - "energy_axes.set_title(f\"AFIR path: {MOLECULE_NAME} (MACE-{MACE_MODEL_FAMILY.upper()}, {MACE_MODEL_SIZE})\")\n", + "energy_axes.set_title(f\"AFIR path: {MOLECULE_NAME} (MACE-MP, {MACE_MODEL})\")\n", "energy_axes.legend()\n", "energy_axes.grid(True, linestyle=\":\", alpha=0.6)\n", "\n", @@ -401,7 +420,7 @@ "product = images[-1].copy()\n", "product.calc = calculator\n", "\n", - "BFGS(product, logfile=None).run(fmax=RELAXATION_FMAX)\n", + "BFGS(product, logfile=OPTIMIZER_LOG).run(fmax=RELAXATION_FMAX)\n", "reaction_energy = (product.get_potential_energy() - reactant_energy) * EV_TO_KCAL_PER_MOL\n", "\n", "print(f\"Product energy: {reaction_energy:.1f} kcal/mol relative to the reactant\\n\")\n", @@ -448,17 +467,21 @@ "dimer_control = DimerControl(\n", " initial_eigenmode_method=\"displacement\",\n", " displacement_method=\"vector\",\n", - " logfile=None,\n", + " logfile=OPTIMIZER_LOG,\n", + " eigenmode_logfile=OPTIMIZER_LOG,\n", ")\n", "dimer = MinModeAtoms(transition_state, dimer_control)\n", "dimer.displace(displacement_vector=0.05 * reaction_direction, mask=[True] * len(transition_state))\n", "\n", - "MinModeTranslate(dimer, logfile=None).run(fmax=SADDLE_FMAX, steps=200)\n", + "MinModeTranslate(dimer, logfile=OPTIMIZER_LOG).run(fmax=SADDLE_FMAX, steps=200)\n", "\n", "transition_state_energy = transition_state.get_potential_energy()\n", "activation_energy = (transition_state_energy - reactant_energy) * EV_TO_KCAL_PER_MOL\n", + "saddle_force = float(np.abs(transition_state.get_forces()).max())\n", "\n", - "print(f\"Maximum force at the saddle point: {np.abs(transition_state.get_forces()).max():.3f} eV/Å\")\n", + "print(f\"Maximum force at the saddle point: {saddle_force:.3f} eV/Å\")\n", + "if saddle_force > SADDLE_FMAX:\n", + " print(f\"⚠️ Not converged to {SADDLE_FMAX} eV/Å — this structure is not a transition state and the numbers below say nothing about the reaction.\")\n", "print(f\"Activation energy: {activation_energy:.1f} kcal/mol\")\n", "print(\n", " f\"d{BOND_FORMING_PAIR} = {transition_state.get_distance(*BOND_FORMING_PAIR):.2f} Å, \"\n", @@ -500,7 +523,10 @@ "for index in imaginary_mode_indices:\n", " print(f\" mode {index}: {abs(frequencies[index].imag):.0f}i cm⁻¹\")\n", "\n", - "reaction_mode = vibrations.get_mode(imaginary_mode_indices[0])\n", + "if not imaginary_mode_indices:\n", + " print(\"⚠️ No imaginary mode above the threshold — this structure is not a transition state.\")\n", + "\n", + "reaction_mode = vibrations.get_mode(imaginary_mode_indices[0] if imaginary_mode_indices else 0)\n", "vibrations.clean()" ] }, @@ -528,14 +554,25 @@ " displaced = transition_state.copy()\n", " displaced.positions += sign * REACTION_MODE_DISPLACEMENT * reaction_mode / np.linalg.norm(reaction_mode)\n", " displaced.calc = calculator\n", - " BFGS(displaced, logfile=None).run(fmax=RELAXATION_FMAX, steps=400)\n", + " BFGS(displaced, logfile=OPTIMIZER_LOG).run(fmax=RELAXATION_FMAX, steps=400)\n", " connected_minima[label] = displaced\n", " energy = (displaced.get_potential_energy() - reactant_energy) * EV_TO_KCAL_PER_MOL\n", " print(\n", " f\"{label:<12}{energy:>18.1f}\"\n", " f\"{displaced.get_distance(*BOND_FORMING_PAIR):>12.2f} Å\"\n", " f\"{displaced.get_distance(*BOND_BREAKING_PAIR):>12.2f} Å\"\n", - " )" + " )\n", + "\n", + "connects_two_minima = any(\n", + " abs(connected_minima[\"forward\"].get_distance(*pair) - connected_minima[\"reverse\"].get_distance(*pair))\n", + " > MINIMUM_SEPARATION\n", + " for pair in (BOND_FORMING_PAIR, BOND_BREAKING_PAIR)\n", + ")\n", + "print(\n", + " \"\\n✅ The imaginary mode connects two distinct minima.\"\n", + " if connects_two_minima\n", + " else \"\\n⚠️ Both directions relax to the same structure — the saddle does not connect a reactant and a product.\"\n", + ")" ] }, { @@ -573,13 +610,14 @@ "axes.set_xticks(range(len(levels)))\n", "axes.set_xticklabels([label for label, _ in levels])\n", "axes.set_ylabel(\"Energy relative to reactant (kcal/mol)\")\n", - "axes.set_title(f\"{MOLECULE_NAME} → {PRODUCT_NAME} (MACE-{MACE_MODEL_FAMILY.upper()}, {MACE_MODEL_SIZE})\")\n", + "axes.set_title(f\"{MOLECULE_NAME} → {PRODUCT_NAME} (MACE-MP, {MACE_MODEL})\")\n", "axes.legend()\n", "axes.grid(True, axis=\"y\", linestyle=\":\", alpha=0.6)\n", "diagram_figure.tight_layout()\n", "display_matplotlib_figure(diagram_figure)\n", "\n", - "print(f\"Activation energy: {activation_energy:.1f} kcal/mol (experiment: {EXPERIMENTAL_ACTIVATION_ENERGY})\")\n", + "saddle_note = \"\" if saddle_force <= SADDLE_FMAX else \" ⚠️ no transition state was found, see 7\"\n", + "print(f\"Activation energy: {activation_energy:.1f} kcal/mol (experiment: {EXPERIMENTAL_ACTIVATION_ENERGY}){saddle_note}\")\n", "print(f\"Reaction energy: {reaction_energy:.1f} kcal/mol\")" ] }, @@ -616,7 +654,7 @@ "## 11. Save the results\n", "### 11.1. Hand the structures back as materials\n", "\n", - "The three structures that define the reaction are passed to the environment with `set_materials`, each carrying the search settings and the resulting energies as metadata." + "The three structures that define the reaction are passed to the environment with `set_materials`." ] }, { @@ -628,42 +666,13 @@ "source": [ "from mat3ra.notebooks_utils.material import set_materials\n", "\n", - "reaction_metadata = {\n", - " \"type\": \"reaction_path\",\n", - " \"method\": {\n", - " \"name\": \"AFIR\",\n", - " \"artificial_force_ramp\": {\"values\": AFIR_FORCE_RAMP, \"units\": \"eV/angstrom\"},\n", - " \"bond_forming_pair\": list(BOND_FORMING_PAIR),\n", - " \"bond_breaking_pair\": list(BOND_BREAKING_PAIR),\n", - " },\n", - " \"engine\": {\"name\": \"ASE\", \"optimizer\": \"BFGS\", \"saddle_search\": \"Dimer\"},\n", - " \"calculator\": {\n", - " \"name\": f\"MACE-{MACE_MODEL_FAMILY.upper()}\",\n", - " \"parameters\": {\"model\": MACE_MODEL_SIZE, \"default_dtype\": MACE_DEFAULT_DTYPE},\n", - " },\n", - " \"reactant\": MOLECULE_NAME,\n", - " \"product\": PRODUCT_NAME,\n", - " \"activation_energy\": {\"value\": round(float(activation_energy), 2), \"units\": \"kcal/mol\"},\n", - " \"reaction_energy\": {\"value\": round(float(reaction_energy), 2), \"units\": \"kcal/mol\"},\n", - " \"imaginary_frequency\": {\n", - " \"value\": round(float(abs(frequencies[imaginary_mode_indices[0]].imag)), 1),\n", - " \"units\": \"cm-1\",\n", - " },\n", - "}\n", - "\n", "structures = (\n", " (\"reactant\", reactant, f\"{MOLECULE_NAME}, reactant\"),\n", " (\"transition_state\", transition_state, f\"{MOLECULE_NAME} to {PRODUCT_NAME}, transition state\"),\n", " (\"product\", product, f\"{PRODUCT_NAME}, product\"),\n", ")\n", "\n", - "materials = []\n", - "for role, atoms, name in structures:\n", - " material = to_material(atoms, name)\n", - " material.metadata.reaction_path = {**reaction_metadata, \"role\": role}\n", - " materials.append(material)\n", - "\n", - "set_materials(materials, FOLDER)" + "set_materials([to_material(atoms, name) for _, atoms, name in structures], FOLDER)" ] }, { @@ -673,7 +682,7 @@ "source": [ "### 11.2. Write the energy profile and the plots\n", "\n", - "The trajectory of the search is already in `RESULTS_FOLDER`; this adds the profile behind the plots, the key energies and the two figures, so the run can be reported without re-running it." + "The trajectory of the search is already in `RESULTS_FOLDER`; this adds the profile behind the plots, the settings and the resulting energies, so the run can be reported without re-running it." ] }, { @@ -686,7 +695,19 @@ "import json\n", "\n", "results = {\n", - " \"reaction_path\": reaction_metadata,\n", + " \"reaction\": {\"reactant\": MOLECULE_NAME, \"product\": PRODUCT_NAME},\n", + " \"settings\": {\n", + " \"calculator\": f\"MACE-MP {MACE_MODEL}\",\n", + " \"bond_forming_pair\": list(BOND_FORMING_PAIR),\n", + " \"bond_breaking_pair\": list(BOND_BREAKING_PAIR),\n", + " \"artificial_force_ramp_ev_per_angstrom\": AFIR_FORCE_RAMP,\n", + " },\n", + " \"activation_energy_kcal_per_mol\": round(float(activation_energy), 2),\n", + " \"reaction_energy_kcal_per_mol\": round(float(reaction_energy), 2),\n", + " \"imaginary_frequency_cm\": (\n", + " round(float(abs(frequencies[imaginary_mode_indices[0]].imag)), 1) if imaginary_mode_indices else None\n", + " ),\n", + " \"transition_state_found\": bool(saddle_force <= SADDLE_FMAX and connects_two_minima),\n", " \"distances_angstrom\": {\n", " role: {\n", " \"forming\": round(float(atoms.get_distance(*BOND_FORMING_PAIR)), 3),\n", From c91f992690d04521d356f2670954ec3ab9067a84 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Wed, 12 Aug 2026 13:55:26 -0700 Subject: [PATCH 4/5] update: better description --- .../local/reaction_path_afir_mace.ipynb | 270 +++++++++++------- 1 file changed, 173 insertions(+), 97 deletions(-) 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 b8914ee7..6b852d8c 100644 --- a/other/materials_designer/workflows/local/reaction_path_afir_mace.ipynb +++ b/other/materials_designer/workflows/local/reaction_path_afir_mace.ipynb @@ -5,24 +5,35 @@ "id": "0", "metadata": {}, "source": [ - "# Discover a reaction path with AFIR + MLFF (MACE)\n", + "# Reaction path and transition state with AFIR + MLFF (MACE)\n", "\n", - "Find how a molecule rearranges **without knowing the product in advance**: an artificial force (AFIR) pulls two selected atoms together on top of the physical potential energy surface, while a **MACE** foundation model supplies energies and forces, so bonds break and form on the fly.\n", + "Discover a reaction path and its transition state locally with ASE. An artificial force (AFIR) pulls two selected atoms together while a MACE foundation model supplies energies and forces, so bonds break and form during the optimization.\n", "\n", - "The example is the **Claisen rearrangement** of allyl vinyl ether (C5H8O) into 4-pentenal — a concerted [3,3]-sigmatropic shift in which a C–O bond breaks and a C–C bond forms in one step. Its ~30 kcal/mol barrier is far too high for a plain relaxation to cross, which is exactly the situation AFIR is built for.\n", + "Example: the Claisen rearrangement of allyl vinyl ether (C5H8O) to 4-pentenal, a concerted [3,3]-sigmatropic shift in which a C–O bond breaks and a C–C bond forms. The barrier of ~30 kcal/mol is too high to cross by relaxation.\n", + "\n", + "The inputs are the molecule, the atom pair to push together and the force ramp. The product, the barrier and the transition state are results, not inputs.\n", "\n", "

Usage

\n", "\n", - "1. Set the molecule, the reacting atom pairs and the MACE model in cells 1.2 and 1.3 (or use default values).\n", + "1. Set the molecule, the atoms to push together and the MACE model in cells 1.2 to 1.5 (or use default values).\n", "1. Click \"Run\" > \"Run All\" to run all cells.\n", - "1. Wait for the search to finish (a few minutes on a laptop CPU).\n", + "1. Wait for the search to finish.\n", "1. Scroll down to view the discovered path, the transition state and the energy diagram.\n", "\n", "## Summary\n", "\n", - "Load the molecule, relax it with MACE, ramp an artificial force between the two target atoms until the reaction happens, strip the bias to recover the physical energy profile, relax the discovered product, refine the highest point of the path into a true saddle point, verify it with a vibrational analysis and by relaxing along the imaginary mode in both directions, and save the reactant, transition state and product together with the energy profile and the plots.\n", + "2. Load the molecule from the uploads folder or from PubChem.\n", + "3. Relax it with MACE to get the reactant.\n", + "4. Ramp an artificial force between the target atoms until the molecule reacts.\n", + "5. Remove the bias to recover the physical energy profile. Its maximum is the first estimate of the transition state.\n", + "6. Relax the end of the path to get the product.\n", + "7. Refine that estimate to a saddle point with the dimer method.\n", + "8. Verify the saddle point has one imaginary frequency.\n", + "9. Follow that mode in both directions to identify the minima it connects.\n", + "10. Plot the energy diagram and view the structures.\n", + "11. Save the structures, the energy profile and the plots.\n", "\n", - "Everything runs locally through ASE — no platform compute and no authentication." + "Runs locally through ASE, without platform compute or authentication." ] }, { @@ -56,9 +67,11 @@ "id": "3", "metadata": {}, "source": [ - "### 1.2. Set the reaction\n", + "### 1.2. Inputs that steer the search\n", + "\n", + "AFIR adds a bias term $E_\\text{bias} = \\alpha \\, r_{ij}$ between the target atoms: a constant attractive force $\\alpha$. Too small a force stops in a biased minimum, too large a force distorts the geometry, so it is ramped, each stage starting from the structure the previous one converged to.\n", "\n", - "The atom pairs are indices into the loaded structure. Cell 2.1 prints the element of every index so they can be checked before the search starts." + "Atom pairs are indices into the loaded structure. Cell 2.1 prints the element of each index." ] }, { @@ -68,17 +81,15 @@ "metadata": {}, "outputs": [], "source": [ - "FOLDER = \"../../uploads\" # molecules are read from here, and the resulting materials are written back here\n", - "RESULTS_FOLDER = \"results\" # energy profile, plots and the raw trajectory\n", - "MOLECULE_NAME = \"allyl vinyl ether\" # looked up in the uploads folder first, then fetched from PubChem by name\n", - "PRODUCT_NAME = \"4-pentenal\"\n", + "# Starting geometry: read from the uploads folder, or fetched from PubChem under this name\n", + "MOLECULE_NAME = \"allyl vinyl ether\"\n", "\n", - "# The bond that AFIR forces to form: the two terminal CH2 carbons\n", + "# The only pair AFIR pulls together; this choice selects which reaction is explored.\n", + "# Two terminal CH2 carbons, whose bond closes the ring of the [3,3] shift.\n", "BOND_FORMING_PAIR = (4, 5)\n", - "# The bond expected to break in response: the ether oxygen and the allylic CH2\n", - "BOND_BREAKING_PAIR = (0, 1)\n", - "# The bond that becomes the product carbonyl\n", - "CARBONYL_PAIR = (0, 3)" + "\n", + "# Strength of that pull, ramped stage by stage until the molecule reacts (eV/Å)\n", + "AFIR_FORCE_RAMP = [1.0, 2.0, 3.0, 4.0]" ] }, { @@ -86,9 +97,9 @@ "id": "5", "metadata": {}, "source": [ - "### 1.3. AFIR and MACE options\n", + "### 1.3. The model\n", "\n", - "AFIR adds a bias term $E_\\text{bias} = \\alpha \\, r_{ij}$ between the target atoms, i.e. a constant attractive force $\\alpha$ that pulls them together no matter how the physical surface resists. A single value of $\\alpha$ either leaves the molecule stuck in a biased minimum (too weak) or drags it through a badly distorted geometry (too strong), so the force is ramped: each stage starts from the structure the previous one converged to." + "The family determines whether the model covers this chemistry, see 3.1." ] }, { @@ -98,48 +109,93 @@ "metadata": {}, "outputs": [], "source": [ - "# Artificial force applied between the target atoms, ramped stage by stage (eV/Å)\n", - "AFIR_FORCE_RAMP = [1.0, 2.0, 3.0, 4.0]\n", - "AFIR_MAX_STEPS_PER_STAGE = 150\n", - "AFIR_MAX_DISPLACEMENT = 0.1 # per-step displacement cap, keeps the biased path smooth (Å)\n", - "AFIR_TRAJECTORY_PATH = f\"{RESULTS_FOLDER}/afir_path.traj\"\n", + "# What the MACE model was trained on: \"organic\" for molecules, \"inorganic\" for crystals and surfaces\n", + "MACE_MODEL_FAMILY = \"organic\"\n", + "MACE_MODEL = \"medium\" # choose between \"small\", \"medium\" and \"large\"\n", + "MACE_DISPERSION = False # D3 dispersion correction, not needed for this intramolecular rearrangement\n", + "MACE_DEFAULT_DTYPE = \"float64\" # float64 is recommended for geometry optimization and vibrations\n", + "MACE_DEVICE = \"cpu\" # hardware target: \"cpu\" or \"cuda\" (GPU)" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "### 1.4. Convergence and verification thresholds\n", "\n", + "These set precision and run time, not which reaction is found." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ "# Maximum force on any atom at convergence (eV/Å)\n", "RELAXATION_FMAX = 0.03\n", "AFIR_FMAX = 0.05\n", "SADDLE_FMAX = 0.02\n", "\n", + "AFIR_MAX_STEPS_PER_STAGE = 150\n", + "AFIR_MAX_DISPLACEMENT = 0.1 # per-step displacement cap, keeps the biased path smooth (Å)\n", + "\n", "# Modes below this magnitude are the translations and rotations of a free molecule (cm⁻¹)\n", "IMAGINARY_MODE_THRESHOLD = 50\n", "# Displacement along the imaginary mode used to leave the saddle point (Å)\n", "REACTION_MODE_DISPLACEMENT = 0.3\n", "# Two relaxed structures count as different minima if a target distance differs by more than this (Å)\n", - "MINIMUM_SEPARATION = 0.5\n", + "MINIMUM_SEPARATION = 0.5" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "### 1.5. Labels, reported distances and output locations\n", "\n", - "# \"mp\": MACE-MP models shipped with the platform, the only option that works in JupyterLite;\n", - "# \"off\": MACE-OFF23 for organic molecules, downloaded on first use, so local runs only\n", - "MACE_MODEL_FAMILY = \"mp\"\n", - "MACE_MODEL = \"medium\" # choose between \"small\", \"medium\" and \"large\"\n", - "MACE_DISPERSION = False # D3 dispersion correction, not needed for this intramolecular rearrangement\n", - "MACE_DEFAULT_DTYPE = \"float64\" # float64 is recommended for geometry optimization and vibrations\n", - "MACE_DEVICE = \"cpu\" # hardware target: \"cpu\" or \"cuda\" (GPU)" + "None of these affect the result. `PRODUCT_NAME` labels whatever the search produces; it is not a target and is not checked against the product, the bond distances in section 6 are. `BOND_BREAKING_PAIR` is reported along the path and sets the initial direction of the saddle search in section 7; AFIR does not act on it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "PRODUCT_NAME = \"4-pentenal\" # label for the discovered product, used in plots and saved names\n", + "\n", + "# Reported along the path: the ether oxygen and the allylic CH2 expected to come apart\n", + "BOND_BREAKING_PAIR = (0, 1)\n", + "# Reported only: the bond that becomes the product carbonyl\n", + "CARBONYL_PAIR = (0, 3)\n", + "\n", + "FOLDER = \"../../uploads\" # molecules are read from here, and the resulting materials are written back here\n", + "RESULTS_FOLDER = \"results\" # energy profile, plots and the raw trajectory\n", + "AFIR_TRAJECTORY_PATH = f\"{RESULTS_FOLDER}/afir_path.traj\"" ] }, { "cell_type": "markdown", - "id": "7", + "id": "11", "metadata": {}, "source": [ "## 2. Load the molecule\n", "### 2.1. Read from uploads, or fetch the 3D structure from PubChem\n", "\n", - "PubChem serves an optimized 3D conformer for most small molecules, which makes any named molecule a valid starting point for the search." + "PubChem provides an optimized 3D conformer for most small molecules, so any molecule name is a valid starting point.\n", + "\n", + "The connectivity listing below is what the atom indices in 1.2 and 1.5 refer to: hydrogens are summarized as a count, so each heavy atom is identified by its neighbours. In this molecule atoms 4 and 5 are the only carbons with a single carbon neighbour and two hydrogens, i.e. the two terminal CH2 groups, and atom 1 is the CH2 attached to the ether oxygen." ] }, { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -148,6 +204,7 @@ "from urllib.parse import quote\n", "\n", "from ase.io import read, write\n", + "from ase.neighborlist import NeighborList, natural_cutoffs\n", "\n", "PUBCHEM_STRUCTURE_URL = \"https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/name/{name}/SDF?record_type=3d\"\n", "\n", @@ -172,15 +229,33 @@ "\n", "molecule = read(molecule_path)\n", "\n", - "print(f\"{MOLECULE_NAME}: {molecule.get_chemical_formula()}\")\n", - "print(\"atoms: \" + \", \".join(f\"{index}:{symbol}\" for index, symbol in enumerate(molecule.get_chemical_symbols())))\n", - "print(f\"bond to form {BOND_FORMING_PAIR}: {molecule.get_distance(*BOND_FORMING_PAIR):.2f} Å\")\n", + "BOND_TOLERANCE = 1.1 # scales covalent radii when deciding which atoms count as bonded, for the listing below\n", + "\n", + "neighbor_list = NeighborList(natural_cutoffs(molecule, mult=BOND_TOLERANCE), self_interaction=False, bothways=True)\n", + "neighbor_list.update(molecule)\n", + "symbols = molecule.get_chemical_symbols()\n", + "selected_pairs = {BOND_FORMING_PAIR: \"forming\", BOND_BREAKING_PAIR: \"breaking\"}\n", + "\n", + "print(f\"{MOLECULE_NAME}: {molecule.get_chemical_formula()}\n", + "\")\n", + "print(f\"{'index':>5} {'atom':<5}bonded to\")\n", + "for index, symbol in enumerate(symbols):\n", + " if symbol == \"H\":\n", + " continue\n", + " neighbors = sorted(int(neighbor) for neighbor in neighbor_list.get_neighbors(index)[0])\n", + " heavy_neighbors = \", \".join(f\"{neighbor}:{symbols[neighbor]}\" for neighbor in neighbors if symbols[neighbor] != \"H\")\n", + " hydrogen_count = sum(1 for neighbor in neighbors if symbols[neighbor] == \"H\")\n", + " role = next((f\"{role} pair {pair}\" for pair, role in selected_pairs.items() if index in pair), \"\")\n", + " print(f\"{index:>5} {symbol:<5}{heavy_neighbors}{f' + {hydrogen_count}H' if hydrogen_count else '':<10}{role:>22}\")\n", + "\n", + "print(f\"\n", + "bond to form {BOND_FORMING_PAIR}: {molecule.get_distance(*BOND_FORMING_PAIR):.2f} Å\")\n", "print(f\"bond to break {BOND_BREAKING_PAIR}: {molecule.get_distance(*BOND_BREAKING_PAIR):.2f} Å\")" ] }, { "cell_type": "markdown", - "id": "9", + "id": "13", "metadata": {}, "source": [ "### 2.2. View the molecule" @@ -189,7 +264,7 @@ { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -213,45 +288,46 @@ }, { "cell_type": "markdown", - "id": "11", + "id": "15", "metadata": {}, "source": [ "## 3. Relax the reactant with MACE\n", "### 3.1. Create the ASE calculator\n", "\n", - "With `MACE_MODEL_FAMILY = \"mp\"` the calculator is built from the MACE-MP models shipped with the platform (`packages/models`), so nothing is downloaded at run time and the notebook runs in JupyterLite.\n", + "Models are loaded from the ones shipped with the platform (`packages/models`), so nothing is downloaded at run time and local and JupyterLite runs use the same weights.\n", + "\n", + "`MACE_MODEL_FAMILY = \"organic\"` selects MACE-OFF23, trained on organic molecules, distributed under the Academic Software License, which does not permit commercial use. `\"inorganic\"` selects MACE-MP-0, trained on Materials Project crystals, MIT licensed.\n", "\n", - "Those models are trained on inorganic crystal trajectories. For this organic rearrangement they close the C–C bond but do not break the C–O bond, so the search ends at a cyclic structure rather than at 4-pentenal. `MACE_MODEL_FAMILY = \"off\"` selects MACE-OFF23, which is trained on organic molecules and reproduces the published mechanism — it is downloaded on first use, so it works locally but not in the browser, and its Academic Software License does not permit commercial use." + "With `\"inorganic\"` this reaction closes the C–C bond without breaking the C–O bond and ends at a cyclic structure, and section 7 reports that no transition state was found. Energies of the two families use different references and are not comparable." ] }, { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "16", "metadata": {}, "outputs": [], "source": [ "from mat3ra.notebooks_utils.mlff import create_mlff_calculator\n", - "\n", - "if MACE_MODEL_FAMILY == \"off\":\n", - " from mace.calculators import mace_off\n", - "\n", - " calculator = mace_off(model=MACE_MODEL, default_dtype=MACE_DEFAULT_DTYPE, device=MACE_DEVICE)\n", - "else:\n", - " calculator = create_mlff_calculator(\n", - " \"mace\",\n", - " {\n", - " \"model\": MACE_MODEL,\n", - " \"dispersion\": MACE_DISPERSION,\n", - " \"default_dtype\": MACE_DEFAULT_DTYPE,\n", - " \"device\": MACE_DEVICE,\n", - " },\n", - " )" + "from mat3ra.notebooks_utils.pyodide.packages.mace import MODEL_FAMILY_LABELS\n", + "\n", + "MACE_MODEL_LABEL = f\"{MODEL_FAMILY_LABELS[MACE_MODEL_FAMILY]} ({MACE_MODEL})\"\n", + "\n", + "calculator = create_mlff_calculator(\n", + " \"mace\",\n", + " {\n", + " \"family\": MACE_MODEL_FAMILY,\n", + " \"model\": MACE_MODEL,\n", + " \"dispersion\": MACE_DISPERSION,\n", + " \"default_dtype\": MACE_DEFAULT_DTYPE,\n", + " \"device\": MACE_DEVICE,\n", + " },\n", + ")" ] }, { "cell_type": "markdown", - "id": "13", + "id": "17", "metadata": {}, "source": [ "### 3.2. Relax the reactant" @@ -260,7 +336,7 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "18", "metadata": {}, "outputs": [], "source": [ @@ -280,18 +356,18 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "19", "metadata": {}, "source": [ "## 4. Push the reaction with an artificial force\n", "\n", - "The bias is applied with ASE's `ExternalForce` constraint, which adds exactly the AFIR term: a constant force of $\\alpha$ along the vector connecting the two target atoms, with an energy contribution $\\alpha \\, r_{ij}$. The physical forces still come from MACE, so the molecule is free to respond by breaking whatever bond is in the way." + "The bias is applied as an ASE `ExternalForce` constraint, which is the AFIR term: a constant force $\\alpha$ along the vector between the target atoms, with energy $\\alpha \\, r_{ij}$. Physical forces come from MACE, so which bond breaks in response is determined by the model." ] }, { "cell_type": "code", "execution_count": null, - "id": "16", + "id": "20", "metadata": {}, "outputs": [], "source": [ @@ -321,19 +397,19 @@ }, { "cell_type": "markdown", - "id": "17", + "id": "21", "metadata": {}, "source": [ "## 5. Recover the physical energy landscape\n", "### 5.1. Strip the bias\n", "\n", - "Every structure along the biased path is re-evaluated with the bias removed, which turns the AFIR trajectory into a physical energy profile. Its maximum is the first estimate of the transition state." + "Each structure along the biased path is re-evaluated without the bias, turning the AFIR trajectory into a physical energy profile. Its maximum is the first estimate of the transition state." ] }, { "cell_type": "code", "execution_count": null, - "id": "18", + "id": "22", "metadata": {}, "outputs": [], "source": [ @@ -360,7 +436,7 @@ }, { "cell_type": "markdown", - "id": "19", + "id": "23", "metadata": {}, "source": [ "### 5.2. Plot the discovered path" @@ -369,7 +445,7 @@ { "cell_type": "code", "execution_count": null, - "id": "20", + "id": "24", "metadata": {}, "outputs": [], "source": [ @@ -384,7 +460,7 @@ "energy_axes.plot(path_energies, color=\"#2b5c8f\", linewidth=2)\n", "energy_axes.axvline(transition_state_guess_index, color=\"#c93b3b\", linestyle=\"--\", label=\"Transition state guess\")\n", "energy_axes.set_ylabel(\"Energy relative to reactant (kcal/mol)\")\n", - "energy_axes.set_title(f\"AFIR path: {MOLECULE_NAME} (MACE-MP, {MACE_MODEL})\")\n", + "energy_axes.set_title(f\"AFIR path: {MOLECULE_NAME}, {MACE_MODEL_LABEL}\")\n", "energy_axes.legend()\n", "energy_axes.grid(True, linestyle=\":\", alpha=0.6)\n", "\n", @@ -402,18 +478,18 @@ }, { "cell_type": "markdown", - "id": "21", + "id": "25", "metadata": {}, "source": [ "## 6. Relax the discovered product\n", "\n", - "The last structure of the biased path is relaxed with the bias removed, which lets it settle into the product basin the search reached." + "The last structure of the biased path is relaxed without the bias, into the product minimum the search reached." ] }, { "cell_type": "code", "execution_count": null, - "id": "22", + "id": "26", "metadata": {}, "outputs": [], "source": [ @@ -435,18 +511,18 @@ }, { "cell_type": "markdown", - "id": "23", + "id": "27", "metadata": {}, "source": [ "## 7. Refine the transition state\n", "\n", - "The maximum of the AFIR path is a point on a biased trajectory, not a stationary point of the physical surface: the forces there are still large and its energy overestimates the barrier. The dimer method follows the lowest-curvature mode uphill and all remaining modes downhill, converging on the first-order saddle point nearest that guess — that saddle is what sets the activation energy. It is started along the reaction direction: the forming pair closing, the breaking pair opening." + "The maximum of the AFIR path lies on a biased trajectory, not on a stationary point of the physical surface: forces there are large and its energy overestimates the barrier. The dimer method follows the lowest-curvature mode uphill and the remaining modes downhill to the nearest first-order saddle point, which defines the activation energy. The search is started along the reaction direction, with the forming pair closing and the breaking pair opening." ] }, { "cell_type": "code", "execution_count": null, - "id": "24", + "id": "28", "metadata": {}, "outputs": [], "source": [ @@ -491,18 +567,18 @@ }, { "cell_type": "markdown", - "id": "25", + "id": "29", "metadata": {}, "source": [ "## 8. Verify the transition state\n", "\n", - "A first-order saddle point has exactly one imaginary vibrational frequency, and its mode is the reaction coordinate. A free molecule also has six translational and rotational modes at (numerically) near-zero frequency, which appear as small imaginary values — `IMAGINARY_MODE_THRESHOLD` separates those from a genuine reaction mode." + "A first-order saddle point has exactly one imaginary frequency, and its mode is the reaction coordinate. A free molecule also has six translational and rotational modes near zero frequency, which appear as small imaginary values; `IMAGINARY_MODE_THRESHOLD` separates them from a reaction mode." ] }, { "cell_type": "code", "execution_count": null, - "id": "26", + "id": "30", "metadata": {}, "outputs": [], "source": [ @@ -532,18 +608,18 @@ }, { "cell_type": "markdown", - "id": "27", + "id": "31", "metadata": {}, "source": [ "## 9. Confirm which minima the saddle connects\n", "\n", - "Displacing along the imaginary mode in both directions and relaxing shows what the transition state actually joins: one side must fall back to the reactant, the other into the product." + "Displacement along the imaginary mode in both directions, followed by relaxation, identifies the two minima the saddle point joins: one returns to the reactant, the other reaches the product." ] }, { "cell_type": "code", "execution_count": null, - "id": "28", + "id": "32", "metadata": {}, "outputs": [], "source": [ @@ -577,19 +653,19 @@ }, { "cell_type": "markdown", - "id": "29", + "id": "33", "metadata": {}, "source": [ "## 10. Results\n", "### 10.1. Energy diagram\n", "\n", - "The experimental activation energy for this reaction is 30.6 kcal/mol in the gas phase, and the rearrangement is strongly exothermic. MACE reproduces the reaction energy well; the barrier is overestimated, as foundation models trained on near-equilibrium structures generally are in the bond-breaking region." + "The experimental activation energy in the gas phase is 30.6 kcal/mol. Foundation models trained on near-equilibrium structures overestimate barriers in the bond-breaking region; reaction energies are reproduced more closely." ] }, { "cell_type": "code", "execution_count": null, - "id": "30", + "id": "34", "metadata": {}, "outputs": [], "source": [ @@ -610,7 +686,7 @@ "axes.set_xticks(range(len(levels)))\n", "axes.set_xticklabels([label for label, _ in levels])\n", "axes.set_ylabel(\"Energy relative to reactant (kcal/mol)\")\n", - "axes.set_title(f\"{MOLECULE_NAME} → {PRODUCT_NAME} (MACE-MP, {MACE_MODEL})\")\n", + "axes.set_title(f\"{MOLECULE_NAME} → {PRODUCT_NAME}, {MACE_MODEL_LABEL}\")\n", "axes.legend()\n", "axes.grid(True, axis=\"y\", linestyle=\":\", alpha=0.6)\n", "diagram_figure.tight_layout()\n", @@ -623,7 +699,7 @@ }, { "cell_type": "markdown", - "id": "31", + "id": "35", "metadata": {}, "source": [ "### 10.2. View the reactant, transition state and product" @@ -632,7 +708,7 @@ { "cell_type": "code", "execution_count": null, - "id": "32", + "id": "36", "metadata": {}, "outputs": [], "source": [ @@ -648,19 +724,19 @@ }, { "cell_type": "markdown", - "id": "33", + "id": "37", "metadata": {}, "source": [ "## 11. Save the results\n", "### 11.1. Hand the structures back as materials\n", "\n", - "The three structures that define the reaction are passed to the environment with `set_materials`." + "Pass the three structures that define the reaction to the environment with `set_materials`." ] }, { "cell_type": "code", "execution_count": null, - "id": "34", + "id": "38", "metadata": {}, "outputs": [], "source": [ @@ -677,18 +753,18 @@ }, { "cell_type": "markdown", - "id": "35", + "id": "39", "metadata": {}, "source": [ "### 11.2. Write the energy profile and the plots\n", "\n", - "The trajectory of the search is already in `RESULTS_FOLDER`; this adds the profile behind the plots, the settings and the resulting energies, so the run can be reported without re-running it." + "Write the settings, the resulting energies and the profile behind the plots next to the trajectory already in `RESULTS_FOLDER`." ] }, { "cell_type": "code", "execution_count": null, - "id": "36", + "id": "40", "metadata": {}, "outputs": [], "source": [ @@ -697,7 +773,7 @@ "results = {\n", " \"reaction\": {\"reactant\": MOLECULE_NAME, \"product\": PRODUCT_NAME},\n", " \"settings\": {\n", - " \"calculator\": f\"MACE-MP {MACE_MODEL}\",\n", + " \"calculator\": MACE_MODEL_LABEL,\n", " \"bond_forming_pair\": list(BOND_FORMING_PAIR),\n", " \"bond_breaking_pair\": list(BOND_BREAKING_PAIR),\n", " \"artificial_force_ramp_ev_per_angstrom\": AFIR_FORCE_RAMP,\n", @@ -735,7 +811,7 @@ }, { "cell_type": "markdown", - "id": "37", + "id": "41", "metadata": {}, "source": [ "## References\n", From 143fd323b51fb1caebcf9c2f9d0f64242733cd19 Mon Sep 17 00:00:00 2001 From: VsevolodX Date: Wed, 12 Aug 2026 17:48:25 -0700 Subject: [PATCH 5/5] update: generalize --- .../local/reaction_path_afir_mace.ipynb | 293 ++++++++---------- 1 file changed, 136 insertions(+), 157 deletions(-) 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 6b852d8c..fc4e8c6a 100644 --- a/other/materials_designer/workflows/local/reaction_path_afir_mace.ipynb +++ b/other/materials_designer/workflows/local/reaction_path_afir_mace.ipynb @@ -11,12 +11,14 @@ "\n", "Example: the Claisen rearrangement of allyl vinyl ether (C5H8O) to 4-pentenal, a concerted [3,3]-sigmatropic shift in which a C–O bond breaks and a C–C bond forms. The barrier of ~30 kcal/mol is too high to cross by relaxation.\n", "\n", - "The inputs are the molecule, the atom pair to push together and the force ramp. The product, the barrier and the transition state are results, not inputs.\n", + "The inputs are the molecule, the atom pair to push together and the force ramp. The product, the barrier and the transition state are results, not inputs: `PRODUCT_NAME` only labels what the search produces and is not checked against it.\n", "\n", "

Usage

\n", "\n", - "1. Set the molecule, the atoms to push together and the MACE model in cells 1.2 to 1.5 (or use default values).\n", - "1. Click \"Run\" > \"Run All\" to run all cells.\n", + "1. Set the molecule and the MACE model in cell 1.2 (or use the default values).\n", + "1. Run the cells up to 2.1: it lists the atoms with their neighbours and reports the pairs currently selected. Nothing before section 3 uses the model, so this takes seconds.\n", + "1. Set the atom pairs in 1.2 from that listing. AFIR pushes exactly these atoms together, so indices left over from another molecule search a different reaction.\n", + "1. Click \"Run\" > \"Run All\" to run the remaining cells.\n", "1. Wait for the search to finish.\n", "1. Scroll down to view the discovered path, the transition state and the energy diagram.\n", "\n", @@ -67,11 +69,9 @@ "id": "3", "metadata": {}, "source": [ - "### 1.2. Inputs that steer the search\n", + "### 1.2. Set parameters\n", "\n", - "AFIR adds a bias term $E_\\text{bias} = \\alpha \\, r_{ij}$ between the target atoms: a constant attractive force $\\alpha$. Too small a force stops in a biased minimum, too large a force distorts the geometry, so it is ramped, each stage starting from the structure the previous one converged to.\n", - "\n", - "Atom pairs are indices into the loaded structure. Cell 2.1 prints the element of each index." + "Atom indices refer to the loaded structure and are listed by cell 2.1. AFIR adds a bias term $E_\\text{bias} = \\alpha \\, r_{ij}$ between the atoms of `BOND_FORMING_PAIR`: a constant attractive force $\\alpha$. Too small a force stops in a biased minimum, too large a force distorts the geometry, so it is ramped, each stage starting from the structure the previous one converged to." ] }, { @@ -81,34 +81,23 @@ "metadata": {}, "outputs": [], "source": [ - "# Starting geometry: read from the uploads folder, or fetched from PubChem under this name\n", + "# Input structure: a file named after the molecule in FOLDER, otherwise fetched from PubChem by that name\n", + "FOLDER = \"../../uploads\"\n", "MOLECULE_NAME = \"allyl vinyl ether\"\n", "\n", - "# The only pair AFIR pulls together; this choice selects which reaction is explored.\n", - "# Two terminal CH2 carbons, whose bond closes the ring of the [3,3] shift.\n", - "BOND_FORMING_PAIR = (4, 5)\n", + "# Atom indices, as listed by cell 2.1\n", + "BOND_FORMING_PAIR = (4, 5) # the only pair AFIR pulls together; this choice selects the reaction\n", + "BOND_BREAKING_PAIR = (0, 1) # expected to come apart; also directs the saddle search. None if not known\n", + "REPORTED_PAIRS = {\"carbonyl\": (0, 3)} # further distances to follow, {} for none\n", "\n", - "# Strength of that pull, ramped stage by stage until the molecule reacts (eV/Å)\n", - "AFIR_FORCE_RAMP = [1.0, 2.0, 3.0, 4.0]" - ] - }, - { - "cell_type": "markdown", - "id": "5", - "metadata": {}, - "source": [ - "### 1.3. The model\n", + "# Artificial force between the target atoms, ramped stage by stage until the molecule reacts (eV/Å)\n", + "AFIR_FORCE_RAMP = [1.0, 2.0, 3.0, 4.0]\n", + "\n", + "# Name of the expected product, used to label results. None if unknown\n", + "PRODUCT_NAME = \"4-pentenal\"\n", + "# Measured activation energy to draw on the diagram (kcal/mol). None to leave it out\n", + "EXPERIMENTAL_ACTIVATION_ENERGY = 30.6\n", "\n", - "The family determines whether the model covers this chemistry, see 3.1." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6", - "metadata": {}, - "outputs": [], - "source": [ "# What the MACE model was trained on: \"organic\" for molecules, \"inorganic\" for crystals and surfaces\n", "MACE_MODEL_FAMILY = \"organic\"\n", "MACE_MODEL = \"medium\" # choose between \"small\", \"medium\" and \"large\"\n", @@ -119,18 +108,18 @@ }, { "cell_type": "markdown", - "id": "7", + "id": "5", "metadata": {}, "source": [ - "### 1.4. Convergence and verification thresholds\n", + "### 1.3. Convergence thresholds and output paths\n", "\n", - "These set precision and run time, not which reaction is found." + "These set precision and run time, not which reaction is found. Each search writes to its own folder, named after the molecule and the pair pushed together, so several molecules or several pairs can be explored side by side." ] }, { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "6", "metadata": {}, "outputs": [], "source": [ @@ -141,47 +130,33 @@ "\n", "AFIR_MAX_STEPS_PER_STAGE = 150\n", "AFIR_MAX_DISPLACEMENT = 0.1 # per-step displacement cap, keeps the biased path smooth (Å)\n", + "SADDLE_MAX_STEPS = 200\n", + "MODE_FOLLOWING_MAX_STEPS = 400\n", "\n", "# Modes below this magnitude are the translations and rotations of a free molecule (cm⁻¹)\n", "IMAGINARY_MODE_THRESHOLD = 50\n", "# Displacement along the imaginary mode used to leave the saddle point (Å)\n", "REACTION_MODE_DISPLACEMENT = 0.3\n", - "# Two relaxed structures count as different minima if a target distance differs by more than this (Å)\n", - "MINIMUM_SEPARATION = 0.5" - ] - }, - { - "cell_type": "markdown", - "id": "9", - "metadata": {}, - "source": [ - "### 1.5. Labels, reported distances and output locations\n", - "\n", - "None of these affect the result. `PRODUCT_NAME` labels whatever the search produces; it is not a target and is not checked against the product, the bond distances in section 6 are. `BOND_BREAKING_PAIR` is reported along the path and sets the initial direction of the saddle search in section 7; AFIR does not act on it." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "10", - "metadata": {}, - "outputs": [], - "source": [ - "PRODUCT_NAME = \"4-pentenal\" # label for the discovered product, used in plots and saved names\n", - "\n", - "# Reported along the path: the ether oxygen and the allylic CH2 expected to come apart\n", - "BOND_BREAKING_PAIR = (0, 1)\n", - "# Reported only: the bond that becomes the product carbonyl\n", - "CARBONYL_PAIR = (0, 3)\n", - "\n", - "FOLDER = \"../../uploads\" # molecules are read from here, and the resulting materials are written back here\n", - "RESULTS_FOLDER = \"results\" # energy profile, plots and the raw trajectory\n", + "# Two relaxed structures count as different minima if a tracked distance differs by more than this (Å)\n", + "MINIMUM_SEPARATION = 0.5\n", + "# Scales covalent radii when deciding which atoms count as bonded, for the listing in 2.1\n", + "BOND_TOLERANCE = 1.1\n", + "\n", + "# Distances followed along the path and reported for every structure\n", + "TRACKED_PAIRS = {\"forming\": BOND_FORMING_PAIR}\n", + "if BOND_BREAKING_PAIR:\n", + " TRACKED_PAIRS[\"breaking\"] = BOND_BREAKING_PAIR\n", + "TRACKED_PAIRS.update(REPORTED_PAIRS)\n", + "\n", + "PRODUCT_LABEL = PRODUCT_NAME or \"product\"\n", + "SEARCH_NAME = f\"{MOLECULE_NAME.replace(' ', '_')}_{BOND_FORMING_PAIR[0]}-{BOND_FORMING_PAIR[1]}\"\n", + "RESULTS_FOLDER = f\"results/{SEARCH_NAME}\" # one folder per search\n", "AFIR_TRAJECTORY_PATH = f\"{RESULTS_FOLDER}/afir_path.traj\"" ] }, { "cell_type": "markdown", - "id": "11", + "id": "7", "metadata": {}, "source": [ "## 2. Load the molecule\n", @@ -189,13 +164,15 @@ "\n", "PubChem provides an optimized 3D conformer for most small molecules, so any molecule name is a valid starting point.\n", "\n", - "The connectivity listing below is what the atom indices in 1.2 and 1.5 refer to: hydrogens are summarized as a count, so each heavy atom is identified by its neighbours. In this molecule atoms 4 and 5 are the only carbons with a single carbon neighbour and two hydrogens, i.e. the two terminal CH2 groups, and atom 1 is the CH2 attached to the ether oxygen." + "The connectivity listing below is what the atom indices in 1.2 refer to: hydrogens are summarized as a count, so each heavy atom is identified by its neighbours. In this molecule atoms 4 and 5 are the only carbons with a single carbon neighbour and two hydrogens, i.e. the two terminal CH2 groups, and atom 1 is the CH2 attached to the ether oxygen.\n", + "\n", + "The selected pairs are then resolved against the structure and reported, so indices left over from another molecule show up here rather than after the search." ] }, { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -228,34 +205,51 @@ " print(f\"Fetched {MOLECULE_NAME} from PubChem, saved to {molecule_path}\")\n", "\n", "molecule = read(molecule_path)\n", - "\n", - "BOND_TOLERANCE = 1.1 # scales covalent radii when deciding which atoms count as bonded, for the listing below\n", + "symbols = molecule.get_chemical_symbols()\n", "\n", "neighbor_list = NeighborList(natural_cutoffs(molecule, mult=BOND_TOLERANCE), self_interaction=False, bothways=True)\n", "neighbor_list.update(molecule)\n", - "symbols = molecule.get_chemical_symbols()\n", - "selected_pairs = {BOND_FORMING_PAIR: \"forming\", BOND_BREAKING_PAIR: \"breaking\"}\n", + "neighbors_of = {index: sorted(int(n) for n in neighbor_list.get_neighbors(index)[0]) for index in range(len(molecule))}\n", + "selected_pairs = {pair: label for label, pair in TRACKED_PAIRS.items()}\n", + "\n", "\n", - "print(f\"{MOLECULE_NAME}: {molecule.get_chemical_formula()}\n", - "\")\n", + "def describe(index):\n", + " hydrogen_count = sum(1 for neighbor in neighbors_of[index] if symbols[neighbor] == \"H\")\n", + " return f\"{index}:{symbols[index]}\" + (f\" +{hydrogen_count}H\" if hydrogen_count else \"\")\n", + "\n", + "\n", + "print(f\"{MOLECULE_NAME}: {molecule.get_chemical_formula()}\\n\")\n", "print(f\"{'index':>5} {'atom':<5}bonded to\")\n", "for index, symbol in enumerate(symbols):\n", " if symbol == \"H\":\n", " continue\n", - " neighbors = sorted(int(neighbor) for neighbor in neighbor_list.get_neighbors(index)[0])\n", - " heavy_neighbors = \", \".join(f\"{neighbor}:{symbols[neighbor]}\" for neighbor in neighbors if symbols[neighbor] != \"H\")\n", - " hydrogen_count = sum(1 for neighbor in neighbors if symbols[neighbor] == \"H\")\n", - " role = next((f\"{role} pair {pair}\" for pair, role in selected_pairs.items() if index in pair), \"\")\n", - " print(f\"{index:>5} {symbol:<5}{heavy_neighbors}{f' + {hydrogen_count}H' if hydrogen_count else '':<10}{role:>22}\")\n", + " heavy_neighbors = \", \".join(f\"{n}:{symbols[n]}\" for n in neighbors_of[index] if symbols[n] != \"H\")\n", + " hydrogen_count = sum(1 for n in neighbors_of[index] if symbols[n] == \"H\")\n", + " bonded = heavy_neighbors + (f\" + {hydrogen_count}H\" if hydrogen_count else \"\")\n", + " role = next((f\"{label} pair {pair}\" for pair, label in selected_pairs.items() if index in pair), \"\")\n", + " print(f\"{index:>5} {symbol:<5}{bonded:<22}{role}\")\n", + "\n", + "out_of_range = {label: pair for label, pair in TRACKED_PAIRS.items() if max(pair) >= len(molecule)}\n", + "if out_of_range:\n", + " raise ValueError(\n", + " f\"{out_of_range} out of range for {MOLECULE_NAME}, which has {len(molecule)} atoms. \"\n", + " \"Set the pairs in 1.2 from the listing above.\"\n", + " )\n", + "\n", + "print()\n", + "for label, pair in TRACKED_PAIRS.items():\n", + " connection = f\"{describe(pair[0])} — {describe(pair[1])}\"\n", + " print(f\"{label:<10}{connection:<24}{molecule.get_distance(*pair):>8.2f} Å\")\n", "\n", - "print(f\"\n", - "bond to form {BOND_FORMING_PAIR}: {molecule.get_distance(*BOND_FORMING_PAIR):.2f} Å\")\n", - "print(f\"bond to break {BOND_BREAKING_PAIR}: {molecule.get_distance(*BOND_BREAKING_PAIR):.2f} Å\")" + "if any(symbols[index] == \"H\" for index in BOND_FORMING_PAIR):\n", + " print(\"\\n⚠️ The forming pair includes a hydrogen: intended for a hydrogen transfer, otherwise pick two heavy atoms.\")\n", + "if BOND_FORMING_PAIR[1] in neighbors_of[BOND_FORMING_PAIR[0]]:\n", + " print(\"\\n⚠️ The forming pair is already bonded. AFIR pushes together atoms that are not bonded yet.\")" ] }, { "cell_type": "markdown", - "id": "13", + "id": "9", "metadata": {}, "source": [ "### 2.2. View the molecule" @@ -264,7 +258,7 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -288,7 +282,7 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "11", "metadata": {}, "source": [ "## 3. Relax the reactant with MACE\n", @@ -304,7 +298,7 @@ { "cell_type": "code", "execution_count": null, - "id": "16", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -327,7 +321,7 @@ }, { "cell_type": "markdown", - "id": "17", + "id": "13", "metadata": {}, "source": [ "### 3.2. Relax the reactant" @@ -336,7 +330,7 @@ { "cell_type": "code", "execution_count": null, - "id": "18", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -356,7 +350,7 @@ }, { "cell_type": "markdown", - "id": "19", + "id": "15", "metadata": {}, "source": [ "## 4. Push the reaction with an artificial force\n", @@ -367,7 +361,7 @@ { "cell_type": "code", "execution_count": null, - "id": "20", + "id": "16", "metadata": {}, "outputs": [], "source": [ @@ -386,18 +380,15 @@ " BFGS(structure, trajectory=trajectory, maxstep=AFIR_MAX_DISPLACEMENT, logfile=OPTIMIZER_LOG).run(\n", " fmax=AFIR_FMAX, steps=AFIR_MAX_STEPS_PER_STAGE\n", " )\n", - " print(\n", - " f\"α = {artificial_force:.1f} eV/Å → \"\n", - " f\"d{BOND_FORMING_PAIR} = {structure.get_distance(*BOND_FORMING_PAIR):.2f} Å, \"\n", - " f\"d{BOND_BREAKING_PAIR} = {structure.get_distance(*BOND_BREAKING_PAIR):.2f} Å\"\n", - " )\n", + " distances = \", \".join(f\"{label} = {structure.get_distance(*pair):.2f} Å\" for label, pair in TRACKED_PAIRS.items())\n", + " print(f\"α = {artificial_force:.1f} eV/Å → {distances}\")\n", "\n", "structure.set_constraint()" ] }, { "cell_type": "markdown", - "id": "21", + "id": "17", "metadata": {}, "source": [ "## 5. Recover the physical energy landscape\n", @@ -409,7 +400,7 @@ { "cell_type": "code", "execution_count": null, - "id": "22", + "id": "18", "metadata": {}, "outputs": [], "source": [ @@ -436,7 +427,7 @@ }, { "cell_type": "markdown", - "id": "23", + "id": "19", "metadata": {}, "source": [ "### 5.2. Plot the discovered path" @@ -445,15 +436,15 @@ { "cell_type": "code", "execution_count": null, - "id": "24", + "id": "20", "metadata": {}, "outputs": [], "source": [ "from matplotlib import pyplot as plt\n", "from mat3ra.notebooks_utils.plot import display_matplotlib_figure\n", "\n", - "forming_distances = [image.get_distance(*BOND_FORMING_PAIR) for image in images]\n", - "breaking_distances = [image.get_distance(*BOND_BREAKING_PAIR) for image in images]\n", + "TRACKED_COLORS = {\"forming\": \"#2e8b57\", \"breaking\": \"#c93b3b\"}\n", + "tracked_distances = {label: [image.get_distance(*pair) for image in images] for label, pair in TRACKED_PAIRS.items()}\n", "\n", "path_figure, (energy_axes, distance_axes) = plt.subplots(2, 1, figsize=(8, 7), sharex=True)\n", "\n", @@ -464,8 +455,8 @@ "energy_axes.legend()\n", "energy_axes.grid(True, linestyle=\":\", alpha=0.6)\n", "\n", - "distance_axes.plot(forming_distances, color=\"#2e8b57\", label=f\"forming {BOND_FORMING_PAIR}\")\n", - "distance_axes.plot(breaking_distances, color=\"#c93b3b\", label=f\"breaking {BOND_BREAKING_PAIR}\")\n", + "for label, distances in tracked_distances.items():\n", + " distance_axes.plot(distances, color=TRACKED_COLORS.get(label), label=f\"{label} {TRACKED_PAIRS[label]}\")\n", "distance_axes.axvline(transition_state_guess_index, color=\"#c93b3b\", linestyle=\"--\")\n", "distance_axes.set_xlabel(\"AFIR step\")\n", "distance_axes.set_ylabel(\"Distance (Å)\")\n", @@ -478,7 +469,7 @@ }, { "cell_type": "markdown", - "id": "25", + "id": "21", "metadata": {}, "source": [ "## 6. Relax the discovered product\n", @@ -489,7 +480,7 @@ { "cell_type": "code", "execution_count": null, - "id": "26", + "id": "22", "metadata": {}, "outputs": [], "source": [ @@ -501,17 +492,13 @@ "\n", "print(f\"Product energy: {reaction_energy:.1f} kcal/mol relative to the reactant\\n\")\n", "print(f\"{'bond':<20}{'reactant':>12}{'product':>12}\")\n", - "for label, pair in (\n", - " (\"forming\", BOND_FORMING_PAIR),\n", - " (\"breaking\", BOND_BREAKING_PAIR),\n", - " (\"carbonyl\", CARBONYL_PAIR),\n", - "):\n", + "for label, pair in TRACKED_PAIRS.items():\n", " print(f\"{label + ' ' + str(pair):<20}{reactant.get_distance(*pair):>10.2f} Å{product.get_distance(*pair):>10.2f} Å\")" ] }, { "cell_type": "markdown", - "id": "27", + "id": "23", "metadata": {}, "source": [ "## 7. Refine the transition state\n", @@ -522,7 +509,7 @@ { "cell_type": "code", "execution_count": null, - "id": "28", + "id": "24", "metadata": {}, "outputs": [], "source": [ @@ -533,7 +520,8 @@ "print(f\"Maximum force at the AFIR guess: {np.abs(transition_state.get_forces()).max():.2f} eV/Å\")\n", "\n", "reaction_direction = np.zeros_like(transition_state.positions)\n", - "for pair, sign in ((BOND_FORMING_PAIR, 1.0), (BOND_BREAKING_PAIR, -1.0)):\n", + "direction_pairs = [(BOND_FORMING_PAIR, 1.0)] + ([(BOND_BREAKING_PAIR, -1.0)] if BOND_BREAKING_PAIR else [])\n", + "for pair, sign in direction_pairs:\n", " unit_vector = transition_state.positions[pair[1]] - transition_state.positions[pair[0]]\n", " unit_vector /= np.linalg.norm(unit_vector)\n", " reaction_direction[pair[0]] += sign * unit_vector\n", @@ -549,7 +537,7 @@ "dimer = MinModeAtoms(transition_state, dimer_control)\n", "dimer.displace(displacement_vector=0.05 * reaction_direction, mask=[True] * len(transition_state))\n", "\n", - "MinModeTranslate(dimer, logfile=OPTIMIZER_LOG).run(fmax=SADDLE_FMAX, steps=200)\n", + "MinModeTranslate(dimer, logfile=OPTIMIZER_LOG).run(fmax=SADDLE_FMAX, steps=SADDLE_MAX_STEPS)\n", "\n", "transition_state_energy = transition_state.get_potential_energy()\n", "activation_energy = (transition_state_energy - reactant_energy) * EV_TO_KCAL_PER_MOL\n", @@ -559,15 +547,12 @@ "if saddle_force > SADDLE_FMAX:\n", " print(f\"⚠️ Not converged to {SADDLE_FMAX} eV/Å — this structure is not a transition state and the numbers below say nothing about the reaction.\")\n", "print(f\"Activation energy: {activation_energy:.1f} kcal/mol\")\n", - "print(\n", - " f\"d{BOND_FORMING_PAIR} = {transition_state.get_distance(*BOND_FORMING_PAIR):.2f} Å, \"\n", - " f\"d{BOND_BREAKING_PAIR} = {transition_state.get_distance(*BOND_BREAKING_PAIR):.2f} Å\"\n", - ")" + "print(\", \".join(f\"{label} = {transition_state.get_distance(*pair):.2f} Å\" for label, pair in TRACKED_PAIRS.items()))" ] }, { "cell_type": "markdown", - "id": "29", + "id": "25", "metadata": {}, "source": [ "## 8. Verify the transition state\n", @@ -578,7 +563,7 @@ { "cell_type": "code", "execution_count": null, - "id": "30", + "id": "26", "metadata": {}, "outputs": [], "source": [ @@ -608,7 +593,7 @@ }, { "cell_type": "markdown", - "id": "31", + "id": "27", "metadata": {}, "source": [ "## 9. Confirm which minima the saddle connects\n", @@ -619,30 +604,28 @@ { "cell_type": "code", "execution_count": null, - "id": "32", + "id": "28", "metadata": {}, "outputs": [], "source": [ "connected_minima = {}\n", "\n", - "print(f\"{'direction':<12}{'energy, kcal/mol':>18}{'d' + str(BOND_FORMING_PAIR):>14}{'d' + str(BOND_BREAKING_PAIR):>14}\")\n", + "columns = \"\".join(f\"{label:>14}\" for label in TRACKED_PAIRS)\n", + "print(f\"{'direction':<12}{'energy, kcal/mol':>18}{columns}\")\n", "for sign, label in ((1.0, \"forward\"), (-1.0, \"reverse\")):\n", " displaced = transition_state.copy()\n", " displaced.positions += sign * REACTION_MODE_DISPLACEMENT * reaction_mode / np.linalg.norm(reaction_mode)\n", " displaced.calc = calculator\n", - " BFGS(displaced, logfile=OPTIMIZER_LOG).run(fmax=RELAXATION_FMAX, steps=400)\n", + " BFGS(displaced, logfile=OPTIMIZER_LOG).run(fmax=RELAXATION_FMAX, steps=MODE_FOLLOWING_MAX_STEPS)\n", " connected_minima[label] = displaced\n", " energy = (displaced.get_potential_energy() - reactant_energy) * EV_TO_KCAL_PER_MOL\n", - " print(\n", - " f\"{label:<12}{energy:>18.1f}\"\n", - " f\"{displaced.get_distance(*BOND_FORMING_PAIR):>12.2f} Å\"\n", - " f\"{displaced.get_distance(*BOND_BREAKING_PAIR):>12.2f} Å\"\n", - " )\n", + " distances = \"\".join(f\"{displaced.get_distance(*pair):>12.2f} Å\" for pair in TRACKED_PAIRS.values())\n", + " print(f\"{label:<12}{energy:>18.1f}{distances}\")\n", "\n", "connects_two_minima = any(\n", " abs(connected_minima[\"forward\"].get_distance(*pair) - connected_minima[\"reverse\"].get_distance(*pair))\n", " > MINIMUM_SEPARATION\n", - " for pair in (BOND_FORMING_PAIR, BOND_BREAKING_PAIR)\n", + " for pair in TRACKED_PAIRS.values()\n", ")\n", "print(\n", " \"\\n✅ The imaginary mode connects two distinct minima.\"\n", @@ -653,7 +636,7 @@ }, { "cell_type": "markdown", - "id": "33", + "id": "29", "metadata": {}, "source": [ "## 10. Results\n", @@ -665,16 +648,14 @@ { "cell_type": "code", "execution_count": null, - "id": "34", + "id": "30", "metadata": {}, "outputs": [], "source": [ - "EXPERIMENTAL_ACTIVATION_ENERGY = 30.6 # kcal/mol, gas phase\n", - "\n", "levels = [\n", " (MOLECULE_NAME, 0.0),\n", " (\"transition state\", activation_energy),\n", - " (PRODUCT_NAME, reaction_energy),\n", + " (PRODUCT_LABEL, reaction_energy),\n", "]\n", "\n", "diagram_figure, axes = plt.subplots(figsize=(7, 4.5))\n", @@ -682,24 +663,26 @@ "for position, (label, energy) in enumerate(levels):\n", " axes.hlines(energy, position - 0.25, position + 0.25, color=\"#2b5c8f\", linewidth=4)\n", " axes.annotate(f\"{energy:.1f}\", (position, energy), textcoords=\"offset points\", xytext=(0, 10), ha=\"center\")\n", - "axes.axhline(EXPERIMENTAL_ACTIVATION_ENERGY, color=\"#c93b3b\", linestyle=\":\", label=\"experimental barrier\")\n", + "if EXPERIMENTAL_ACTIVATION_ENERGY:\n", + " axes.axhline(EXPERIMENTAL_ACTIVATION_ENERGY, color=\"#c93b3b\", linestyle=\":\", label=\"measured barrier\")\n", "axes.set_xticks(range(len(levels)))\n", "axes.set_xticklabels([label for label, _ in levels])\n", "axes.set_ylabel(\"Energy relative to reactant (kcal/mol)\")\n", - "axes.set_title(f\"{MOLECULE_NAME} → {PRODUCT_NAME}, {MACE_MODEL_LABEL}\")\n", + "axes.set_title(f\"{MOLECULE_NAME} → {PRODUCT_LABEL}, {MACE_MODEL_LABEL}\")\n", "axes.legend()\n", "axes.grid(True, axis=\"y\", linestyle=\":\", alpha=0.6)\n", "diagram_figure.tight_layout()\n", "display_matplotlib_figure(diagram_figure)\n", "\n", "saddle_note = \"\" if saddle_force <= SADDLE_FMAX else \" ⚠️ no transition state was found, see 7\"\n", - "print(f\"Activation energy: {activation_energy:.1f} kcal/mol (experiment: {EXPERIMENTAL_ACTIVATION_ENERGY}){saddle_note}\")\n", + "measured = f\" (measured: {EXPERIMENTAL_ACTIVATION_ENERGY})\" if EXPERIMENTAL_ACTIVATION_ENERGY else \"\"\n", + "print(f\"Activation energy: {activation_energy:.1f} kcal/mol{measured}{saddle_note}\")\n", "print(f\"Reaction energy: {reaction_energy:.1f} kcal/mol\")" ] }, { "cell_type": "markdown", - "id": "35", + "id": "31", "metadata": {}, "source": [ "### 10.2. View the reactant, transition state and product" @@ -708,7 +691,7 @@ { "cell_type": "code", "execution_count": null, - "id": "36", + "id": "32", "metadata": {}, "outputs": [], "source": [ @@ -716,7 +699,7 @@ " [\n", " {\"material\": to_material(reactant, MOLECULE_NAME), \"title\": MOLECULE_NAME},\n", " {\"material\": to_material(transition_state, \"Transition state\"), \"title\": \"Transition state\"},\n", - " {\"material\": to_material(product, PRODUCT_NAME), \"title\": PRODUCT_NAME},\n", + " {\"material\": to_material(product, PRODUCT_LABEL), \"title\": PRODUCT_LABEL},\n", " ],\n", " viewer=ViewersEnum.wave,\n", ")" @@ -724,7 +707,7 @@ }, { "cell_type": "markdown", - "id": "37", + "id": "33", "metadata": {}, "source": [ "## 11. Save the results\n", @@ -736,7 +719,7 @@ { "cell_type": "code", "execution_count": null, - "id": "38", + "id": "34", "metadata": {}, "outputs": [], "source": [ @@ -744,8 +727,8 @@ "\n", "structures = (\n", " (\"reactant\", reactant, f\"{MOLECULE_NAME}, reactant\"),\n", - " (\"transition_state\", transition_state, f\"{MOLECULE_NAME} to {PRODUCT_NAME}, transition state\"),\n", - " (\"product\", product, f\"{PRODUCT_NAME}, product\"),\n", + " (\"transition_state\", transition_state, f\"{MOLECULE_NAME}, transition state\"),\n", + " (\"product\", product, f\"{PRODUCT_LABEL} from {MOLECULE_NAME}\"),\n", ")\n", "\n", "set_materials([to_material(atoms, name) for _, atoms, name in structures], FOLDER)" @@ -753,7 +736,7 @@ }, { "cell_type": "markdown", - "id": "39", + "id": "35", "metadata": {}, "source": [ "### 11.2. Write the energy profile and the plots\n", @@ -764,7 +747,7 @@ { "cell_type": "code", "execution_count": null, - "id": "40", + "id": "36", "metadata": {}, "outputs": [], "source": [ @@ -774,8 +757,7 @@ " \"reaction\": {\"reactant\": MOLECULE_NAME, \"product\": PRODUCT_NAME},\n", " \"settings\": {\n", " \"calculator\": MACE_MODEL_LABEL,\n", - " \"bond_forming_pair\": list(BOND_FORMING_PAIR),\n", - " \"bond_breaking_pair\": list(BOND_BREAKING_PAIR),\n", + " \"tracked_pairs\": {label: list(pair) for label, pair in TRACKED_PAIRS.items()},\n", " \"artificial_force_ramp_ev_per_angstrom\": AFIR_FORCE_RAMP,\n", " },\n", " \"activation_energy_kcal_per_mol\": round(float(activation_energy), 2),\n", @@ -785,18 +767,15 @@ " ),\n", " \"transition_state_found\": bool(saddle_force <= SADDLE_FMAX and connects_two_minima),\n", " \"distances_angstrom\": {\n", - " role: {\n", - " \"forming\": round(float(atoms.get_distance(*BOND_FORMING_PAIR)), 3),\n", - " \"breaking\": round(float(atoms.get_distance(*BOND_BREAKING_PAIR)), 3),\n", - " \"carbonyl\": round(float(atoms.get_distance(*CARBONYL_PAIR)), 3),\n", - " }\n", + " role: {label: round(float(atoms.get_distance(*pair)), 3) for label, pair in TRACKED_PAIRS.items()}\n", " for role, atoms, _ in structures\n", " },\n", " \"afir_path\": {\n", " \"transition_state_guess_index\": transition_state_guess_index,\n", " \"energy_kcal_per_mol\": [round(float(value), 4) for value in path_energies],\n", - " \"forming_distance_angstrom\": [round(float(value), 3) for value in forming_distances],\n", - " \"breaking_distance_angstrom\": [round(float(value), 3) for value in breaking_distances],\n", + " \"distance_angstrom\": {\n", + " label: [round(float(value), 3) for value in distances] for label, distances in tracked_distances.items()\n", + " },\n", " },\n", "}\n", "\n", @@ -811,7 +790,7 @@ }, { "cell_type": "markdown", - "id": "41", + "id": "37", "metadata": {}, "source": [ "## References\n",