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", 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..fc4e8c6a --- /dev/null +++ b/other/materials_designer/workflows/local/reaction_path_afir_mace.ipynb @@ -0,0 +1,822 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Reaction path and transition state with AFIR + MLFF (MACE)\n", + "\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", + "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: `PRODUCT_NAME` only labels what the search produces and is not checked against it.\n", + "\n", + "

Usage

\n", + "\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", + "## Summary\n", + "\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", + "Runs locally through ASE, without platform compute or 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.mlff import get_mlff_install_profiles\n", + "from mat3ra.notebooks_utils.packages import install_packages\n", + "\n", + "await install_packages(get_mlff_install_profiles(\"mace\"))\n", + "\n", + "from mat3ra.notebooks_utils.pyodide.packages.patches import apply_all_patches\n", + "\n", + "apply_all_patches(\"mace\")" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "### 1.2. Set parameters\n", + "\n", + "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." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "# 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", + "# 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", + "# 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", + "# 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": "5", + "metadata": {}, + "source": [ + "### 1.3. Convergence thresholds and output paths\n", + "\n", + "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": "6", + "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", + "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 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": "7", + "metadata": {}, + "source": [ + "## 2. Load the molecule\n", + "### 2.1. Read from uploads, or fetch the 3D structure from PubChem\n", + "\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 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": "8", + "metadata": {}, + "outputs": [], + "source": [ + "import io\n", + "import os\n", + "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", + "\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", + "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", + "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", + "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", + " 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", + "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": "9", + "metadata": {}, + "source": [ + "### 2.2. View the molecule" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "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", + "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 = MaterialWithBuildMetadata.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", + "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", + "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", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.mlff import create_mlff_calculator\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", + "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", + "# 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=OPTIMIZER_LOG).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 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", + "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=OPTIMIZER_LOG).run(\n", + " fmax=AFIR_FMAX, steps=AFIR_MAX_STEPS_PER_STAGE\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": "17", + "metadata": {}, + "source": [ + "## 5. Recover the physical energy landscape\n", + "### 5.1. Strip the bias\n", + "\n", + "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", + "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", + "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", + "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_MODEL_LABEL}\")\n", + "energy_axes.legend()\n", + "energy_axes.grid(True, linestyle=\":\", alpha=0.6)\n", + "\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", + "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 without the bias, into the product minimum 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=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", + "print(f\"{'bond':<20}{'reactant':>12}{'product':>12}\")\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": "23", + "metadata": {}, + "source": [ + "## 7. Refine the transition state\n", + "\n", + "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", + "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", + "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", + " 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=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=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", + "saddle_force = float(np.abs(transition_state.get_forces()).max())\n", + "\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(\", \".join(f\"{label} = {transition_state.get_distance(*pair):.2f} Å\" for label, pair in TRACKED_PAIRS.items()))" + ] + }, + { + "cell_type": "markdown", + "id": "25", + "metadata": {}, + "source": [ + "## 8. Verify the transition state\n", + "\n", + "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", + "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", + "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()" + ] + }, + { + "cell_type": "markdown", + "id": "27", + "metadata": {}, + "source": [ + "## 9. Confirm which minima the saddle connects\n", + "\n", + "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", + "metadata": {}, + "outputs": [], + "source": [ + "connected_minima = {}\n", + "\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=MODE_FOLLOWING_MAX_STEPS)\n", + " connected_minima[label] = displaced\n", + " energy = (displaced.get_potential_energy() - reactant_energy) * EV_TO_KCAL_PER_MOL\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 TRACKED_PAIRS.values()\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", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "29", + "metadata": {}, + "source": [ + "## 10. Results\n", + "### 10.1. Energy diagram\n", + "\n", + "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", + "metadata": {}, + "outputs": [], + "source": [ + "levels = [\n", + " (MOLECULE_NAME, 0.0),\n", + " (\"transition state\", activation_energy),\n", + " (PRODUCT_LABEL, 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", + "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_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", + "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": "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_LABEL), \"title\": PRODUCT_LABEL},\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", + "Pass the three structures that define the reaction to the environment with `set_materials`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "34", + "metadata": {}, + "outputs": [], + "source": [ + "from mat3ra.notebooks_utils.material import set_materials\n", + "\n", + "structures = (\n", + " (\"reactant\", reactant, f\"{MOLECULE_NAME}, reactant\"),\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)" + ] + }, + { + "cell_type": "markdown", + "id": "35", + "metadata": {}, + "source": [ + "### 11.2. Write the energy profile and the plots\n", + "\n", + "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", + "metadata": {}, + "outputs": [], + "source": [ + "import json\n", + "\n", + "results = {\n", + " \"reaction\": {\"reactant\": MOLECULE_NAME, \"product\": PRODUCT_NAME},\n", + " \"settings\": {\n", + " \"calculator\": MACE_MODEL_LABEL,\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", + " \"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: {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", + " \"distance_angstrom\": {\n", + " label: [round(float(value), 3) for value in distances] for label, distances in tracked_distances.items()\n", + " },\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 +}