{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Constrained minima hopping (global optimization)\n\nThis is an example of a search for a global optimum geometric\nconfiguration using the minima hopping algorithm,\nalong with the Hookean class of constraints.\nThis type of approach is useful in searching for the global\noptimum position of adsorbates on a surface while enforcing\nthat the adsorbates' identity is preserved.\nIn short, this example searches for a global optimum configuration using the\n**minima hopping** algorithm together with **Hookean** constraints to\npreserve molecular identity of an adsorbate.\n\nWe look for the optimum binding configuration of a $\\mathrm{Cu_2}$\ndimer on a fixed Pt(110) surface (Cu/Pt chosen only because EMT supports them).\nReplace the $\\mathrm{Cu_2}$ dimer with, e.g., CO to find its optimal site\nwhile preventing dissociation into separate C and O adsorbates.\n\nTwo **Hookean** constraints are used:\n\n1. **Bond-preserving**: apply a restorative force if the Cu\u2013Cu distance\n   exceeds 2.6 \u00c5 (no force below 2.6 \u00c5), keeping the dimer intact.\n2. **Keep below plane**: apply a downward force if one Cu goes above\n   $z=15$ \u00c5 to avoid the dimer flying off into vacuum.\n\n**Outputs.** The run writes a text log (``hop.log``) and a trajectory of\naccepted minima (``minima.traj``). You can also visualize progress with the\n``mhsummary.py`` utility if available.\n\n## References\n- Minima hopping in ASE usage:\n  :mod:`ase.optimize.minimahopping` (Goedecker, JCP 120, 9911 (2004))\n- Hookean constraints: :mod:`ase.constraints.Hookean`\n- Constrained minima hopping tutorial page in the ASE docs\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Imports and calculator\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from __future__ import annotations\n\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nfrom ase import Atoms\nfrom ase.build import fcc110\nfrom ase.calculators.emt import EMT\nfrom ase.constraints import FixAtoms, Hookean\nfrom ase.io import read\nfrom ase.optimize.minimahopping import MHPlot, MinimaHopping\nfrom ase.visualize.plot import plot_atoms\n\n\n# Make results reproducible across doc builds. We will pass this random\n# number generator to the MinimaHopping algorithm.\nseed = 42\nrng = np.random.default_rng(seed)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build a fixed Pt(110) slab\nA small slab is enough for the demonstration and keeps CI fast.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "slab = fcc110(\n    'Pt', size=(3, 2, 3), vacuum=12.0\n)  # (x, y, z repeats), ~ few hundred atoms at most\nprint(slab.pbc)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This sets the desired periodic boundary conditions to be periodic in\nx- and y-direction.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Fix the whole slab (only the adsorbate will move in this minimal example).\nfix = FixAtoms(mask=slab.symbols == 'Pt')\nslab.set_constraint(fix)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, we can visualize the slab.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots()\nplot_atoms(slab, ax, rotation=('270x,0y,0z'))\nax.set_axis_off()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can also attach a calculator.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Attach an EMT calculator\nslab.calc = EMT()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Add a Cu\\ :sub:`2` dimer above the surface\nPlace the dimer roughly above a surface hollow; exact site is not important\nfor the example since minima hopping will explore.\nStart with a 2.3 \u00c5 Cu\u2013Cu distance (near gas-phase value).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "top_z = slab.positions[:, 2].max()\ncu2 = Atoms('Cu2', positions=[[0.0, 0.0, top_z + 3.0], [2.3, 0.0, top_z + 3.0]])\n\n# Translate laterally to the middle of the cell\ncell = slab.get_cell()\nxy_center = 0.5 * (cell[0] + cell[1])\ncu2.translate([xy_center[0], xy_center[1], 0.0])\n\natoms = slab + cu2\natoms.calc = EMT()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, we can visualize the new structure.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots()\nplot_atoms(atoms, ax, rotation=('270x,0y,0z'))\nax.set_axis_off()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Hookean constraints\nAs mentioned above, we want to add two kinds of Hookean constraints:\n\n1) Preserve the Cu\u2013Cu bond if it stretches beyond 2.6 \u00c5:\n   no force for r <= 2.6 \u00c5, Hookean spring k*(r-rt) beyond.\n2) Keep one Cu below the plane z = 15 \u00c5 using the plane form of Hookean:\n   plane (A,B,C,D) with Ax+By+Cz+D = 0  ->  (0,0,1,-15) gives z=15;\n   a downward force is applied if z > 15.\n\nThe Hookean API (ASE >= 3.22) accepts:\n  Hookean(a1, a2, k, rt=None)\n  - a2 can be an atom index, a fixed point (x,y,z), or a plane (A,B,C,D).\n\nSpring constants here are modest to guide but not dominate the dynamics.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "i_cu0, i_cu1 = atoms.symbols.search('Cu')\n\n# contraint #1\nbond_constraint = Hookean(\n    a1=i_cu0, a2=int(i_cu1), k=5.0, rt=2.6\n)  # eV/\u00c5^2, threshold 2.6 \u00c5\n# contraint #2\nz_plane_constraint = Hookean(\n    a1=i_cu0, a2=(0.0, 0.0, 1.0, -15.0), k=2.0\n)  # plane z=15 \u00c5\n\natoms.set_constraint([fix, bond_constraint, z_plane_constraint])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Run minima hopping (short demo)\nWe keep totalsteps small so the example runs quickly in CI.\nOutputs:\n\n- hop.log: text progress\n- minima.traj: accepted local minima\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mh = MinimaHopping(\n    atoms,\n    T0=1500.0,  # initial MD \"temperature\" (K)\n    Ediff0=1.0,  # initial acceptance window (eV)\n    mdmin=2,  # MD stop criterion\n    logfile='hop.log',\n    minima_traj='minima.traj',\n    rng=rng,\n)\n\n# Run a few steps only for documentation builds.\nmh(totalsteps=10)\nprint(\n    'Minima hopping finished.'\n    \"See 'hop.log' and 'minima.traj' in the working directory.\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Visualization of Results\nNow, we can visualize the new structure.\nFor this we are loading the forth image of the trajectory.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "atoms = read('minima.traj')\n\nfig, ax = plt.subplots()\nplot_atoms(atoms, ax, rotation=('270x,0y,0z'))\nax.set_axis_off()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can also analyze this with\n:class:`~ase.optimize.minimahopping.MHPlot` and save the results.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mhplot = MHPlot()\nmhplot.save_figure('summary.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This will make a summary figure, which should see above.\nAs the search is inherently random, yours will look different\nthan this if you used a different random seed\n(and this will look different each time the documentation is rebuilt).\nIn this figure, you will see on the $E_\\mathrm{pot}$ axes the energy\nlevels of the conformers found. The flat bars represent the energy at the\nend of each local optimization step. The checkmark indicates the local\nminimum was accepted; red arrows indicate it was rejected for the three\npossible reasons. The black path between steps is the potential energy\nduring the molecular dynamics (MD) portion of the step; the dashed line\nis the local optimization on termination of the MD step.\nNote the y axis is broken to allow different energy scales\nbetween the local minima and the space explored in the MD simulations.\nThe $T$ and $E_\\mathrm{diff}$ plots show the values of\nthe self-adjusting parameters as the algorithm progresses.\n\n## Further Examples\nYou can find an example of the implementation of this for real adsorbates\nas well as find suitable parameters for the Hookean constraints:\n\nAndrew Peterson\n:doi:`Global optimization of adsorbate\u2013surface\nstructures while preserving molecular identity <10.1007/s11244-013-0161-8>`\nTop. Catal., Vol. **57**, 40 (2014)\n\n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.13.7"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}