{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# NEB with IDPP: |subst|\n\n.. |subst| replace:: Image Dependent Pair Potential\n                     for improved interpolation of NEB initial guess\n\nReference: S. Smidstrup, A. Pedersen, K. Stokbro and H. Jonsson,\n:doi:`Improved initial guess for minimum energy path calculations\n<10.1063/1.4878664>`,\nJ. Chem. Phys. 140, 214106 (2014).\n\nUse of the Nudged Elastic Band\n(NEB) method for transition state search\nis dependent upon generating an initial guess\nfor the images lying between the initial and final states. The most\nsimple approach is to use linear interpolation of the\natomic coordinates. However, this can be problematic as the quality\nof the interpolated path can ofter be far from the real one.\nThe implication being\nthat a lot of time is spent in the NEB routine optimising the shape of\nthe path, before the transition state is homed-in upon.\n\nThe image dependent pair potential (IDPP) is a method that has been\ndeveloped to provide an improvement to the initial guess for the NEB path.\nThe IDPP method uses the bond distance between the atoms involved in\nthe transition state to create target structures for the images, rather\nthan interpolating the atomic positions.\nBy defining an objective function in terms\nof the distances between atoms, the NEB algorithm is used with this\nimage dependent pair potential to create the initial guess for the\nfull NEB calculation.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>The examples below utilise the EMT calculator for illustrative purposes, the\n   results should not be over interpreted.</p></div>\n\nThis tutorial includes example NEB calculations for two different systems.\nFirst, it starts with a simple NEB of Ethane comparing IDPP\nto the standard linear approach.\nThe second example is for a N atom on a Pt step edge.\n\n## Example 1: Ethane\n\nThis example illustrates the use of the IDPP interpolation scheme to\ngenerate an initial guess for rotation of a methyl group around the CC bond.\n\n\n### 1.1 Generate Initial and Final State\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from ase.build import molecule\nfrom ase.calculators.emt import EMT\nfrom ase.mep import NEB\nfrom ase.optimize.fire import FIRE\n\n# Create the molecule.\ninitial = molecule('C2H6')\n# Attach calculators\ninitial.calc = EMT()\n# Relax the molecule\nrelax = FIRE(initial)\nrelax.run(fmax=0.05)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's look at the relaxed molecule.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n\nfrom ase.visualize.plot import plot_atoms\n\nfig, ax = plt.subplots()\nplot_atoms(initial, ax, rotation=('90x,0y,90z'))\nax.set_axis_off()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we can create the final state.\nSince we want to look at the rotation\nof the methyl group, we switch the position of the\nHydrogen atoms on one methyl group.\nThen, we setup and run the NEB calculation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "final = initial.copy()\nfinal.positions[2:5] = initial.positions[[3, 4, 2]]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### 1.2 Linear Interpolation Approach\n\nGenerate blank images.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "images = [initial]\n\nfor i in range(9):\n    images.append(initial.copy())\n\nfor image in images:\n    image.calc = EMT()\n\nimages.append(final)\n\n# Run linear interpolation.\nneb = NEB(images)\nneb.interpolate()\n\n# Run NEB calculation.\nqn = FIRE(neb, trajectory='ethane_linear.traj')\nqn.run(fmax=0.05)\n# You can add a logfile to the FIRE optimizer by adding\n# e.g. `logfile='ethane_linear.log` to save the printed output."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Using the standard linear interpolation approach,\nas in the following example, we can see\nthat 47 iterations are required to find the transition state.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### 1.3 Image Dependent Pair Potential\n\nHowever if we modify our script slightly and use the IDPP method to\nfind the initial guess, we can see that the number of iterations\nrequired to find the transition state is reduced to 7.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Optimise molecule.\ninitial = molecule('C2H6')\ninitial.calc = EMT()\nrelax = FIRE(initial, logfile='opt.log')\nrelax.run(fmax=0.05)\n\n# Create final state.\nfinal = initial.copy()\nfinal.positions[2:5] = initial.positions[[3, 4, 2]]\n\n# Generate blank images.\nimages = [initial]\n\nfor i in range(9):\n    images.append(initial.copy())\n\nfor image in images:\n    image.calc = EMT()\n\nimages.append(final)\n\n# Run IDPP interpolation.\nneb = NEB(images)\nneb.interpolate('idpp')\n\n# Run NEB calculation.\nqn = FIRE(neb, trajectory='ethane_idpp.traj')\nqn.run(fmax=0.05)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Clearly, if one was using a full DFT calculator one can\npotentially gain a significant time improvement.\n\n## Example 2: N Diffusion over a Step Edge\n\nOften we are interested in generating an initial guess for a surface reaction.\n\n### 2.1 Generate Initial and Final State\nThe first part of this example\nillustrates how we can optimise our initial and final state structures\nbefore using the IDPP interpolation to generate our initial guess\nfor the NEB calculation:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n\nfrom ase import Atoms\nfrom ase.calculators.emt import EMT\nfrom ase.constraints import FixAtoms\nfrom ase.lattice.cubic import FaceCenteredCubic\nfrom ase.mep import NEB\nfrom ase.optimize.fire import FIRE as FIRE\n\n# Set the number of images you want.\nnimages = 5\n\n# Some algebra to determine surface normal and the plane of the surface.\nd3 = [2, 1, 1]\na1 = np.array([0, 1, 1])\nd1 = np.cross(a1, d3)\na2 = np.array([0, -1, 1])\nd2 = np.cross(a2, d3)\n\n# Create the slab.\nslab = FaceCenteredCubic(\n    directions=[d1, d2, d3], size=(2, 1, 2), symbol=('Pt'), latticeconstant=3.9\n)\n\n# Add some vacuum to the slab.\nuc = slab.get_cell()\nuc[2] += [0.0, 0.0, 10.0]  # There are ten layers of vacuum.\nuc = slab.set_cell(uc, scale_atoms=False)\n\n# Some positions needed to place the atom in the correct place.\nx1 = 1.379\nx2 = 4.137\nx3 = 2.759\ny1 = 0.0\ny2 = 2.238\nz1 = 7.165\nz2 = 6.439\n\n\n# Add the adatom to the list of atoms and set constraints of surface atoms.\nslab += Atoms('N', [((x2 + x1) / 2, y1, z1 + 1.5)])\nFixAtoms(mask=slab.symbols == 'Pt')\n\n# Optimise the initial state: atom below step.\ninitial = slab.copy()\ninitial.calc = EMT()\nrelax = FIRE(initial, logfile='opt.log')\nrelax.run(fmax=0.05)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now let's visualize this.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots()\nplot_atoms(initial, ax, rotation=('0x,0y,0z'))\nax.set_axis_off()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now create and optimise the final state by\nmoving the atom above the step.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "slab[-1].position = (x3, y2 + 1.0, z2 + 3.5)\nfinal = slab.copy()\nfinal.calc = EMT()\nrelax = FIRE(final, logfile='opt.log')\nrelax.run(fmax=0.05)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now let's visualize this.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots()\nplot_atoms(final, ax, rotation=('0x,0y,0z'))\nax.set_axis_off()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### 2.2 Image Dependent Pair Potential\n\nNow we are ready to setup the NEB with the IDPP interpolation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Create a list of images for interpolation.\nimages = [initial]\nfor i in range(nimages):\n    images.append(initial.copy())\n\nfor image in images:\n    image.calc = EMT()\n\nimages.append(final)\n\n# Carry out idpp interpolation.\nneb = NEB(images)\nneb.interpolate('idpp')\n\n# Run NEB calculation.\nqn = FIRE(neb, trajectory='N_diffusion.traj')\nqn.run(fmax=0.05)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### 2.3 Linear Interpolation Approach\n\nTo again illustrate the potential speedup, the following script\nuses the linear interpolation.\nThis takes more iterations to find a transition\nstate, compared to using the IDPP interpolation.\nWe start from the initial and final state we generated above.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Create a list of images for interpolation.\nimages = [initial]\nfor i in range(nimages):\n    images.append(initial.copy())\n\nfor image in images:\n    image.calc = EMT()\n\nimages.append(final)\n\n# Carry out linear interpolation.\nneb = NEB(images)\nneb.interpolate()\n\n# Run NEB calculation.\nqn = FIRE(neb, trajectory='N_diffusion_lin.traj')\nqn.run(fmax=0.05)"
      ]
    }
  ],
  "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
}