{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Dissociation of a molecule using the NEB method\n\nIn this tutorial we provide an illustrative\nexample of a nudged-elastic band (NEB) calculation.\nFor more information on the NEB technique, see :mod:`ase.mep.neb`.\n\nWe consider the dissociation of a nitrogen molecule on the Cu (111) surface.\n\nThe first step is to find the relaxed structures of\nthe initial and final states.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.animation as animation\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nfrom ase import Atoms\nfrom ase.build import add_adsorbate, fcc111\nfrom ase.calculators.emt import EMT\nfrom ase.constraints import FixAtoms\nfrom ase.io import read, write\nfrom ase.mep import NEB\nfrom ase.optimize import QuasiNewton\nfrom ase.optimize.fire import FIRE\nfrom ase.visualize.plot import plot_atoms"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First, we create the initial state: An N\\ :sub:`2`\\  molecule on a\nCu(111) slab\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Set up a (4 x 4) two layer Cu (111) slab\nslab = fcc111('Cu', size=(4, 4, 2))\nslab.set_pbc((1, 1, 0))\n\n# Add the N2 Molecule,  oriented at 60 degrees:\nd = 1.10  # N2 bond length\nn2_mol = Atoms(\n    'N2', positions=[[0.0, 0.0, 0.0], [0.5 * 3**0.5 * d, 0.5 * d, 0.0]]\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then we put the adsorbate onto the surface\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "add_adsorbate(slab, n2_mol, height=1.0, position='fcc')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And add a calculator for the forces and energies:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "slab.calc = EMT()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We don't want to worry about the Cu degrees of freedom,\nso fix these atoms:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mask = [atom.symbol == 'Cu' for atom in slab]\nslab.set_constraint(FixAtoms(mask=mask))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then we relax the structure and write the trajectory into a file\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "relax = QuasiNewton(slab)\nrelax.run(fmax=0.05)\nprint('initial state:', slab.get_potential_energy())\nwrite('N2.traj', slab)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now the final state.\nMove the second N atom to a neighboring hollow site and relax.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "slab[-1].position[0] = slab[-2].position[0] + 0.25 * slab.cell[0, 0]\nslab[-1].position[1] = slab[-2].position[1]\n\nrelax.run()\nprint('final state:  ', slab.get_potential_energy())\nwrite('2N.traj', slab)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Having obtained these structures we set up an NEB\ncalculation with 9 images.  Using :func:`~neb.interpolate()`\nprovides a guess for the path between the initial\nand final states.  We perform the relaxation of the images\nand obtain the intermediate steps.\n\nFirst, we read the previous configurations\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "initial = read('N2.traj')\nfinal = read('2N.traj')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then, we make 9 images (note the use of copy) and, as before,\nwe fix the Cu atoms\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "configs = [initial.copy() for i in range(8)] + [final]\n\nconstraint = FixAtoms(mask=[atom.symbol != 'N' for atom in initial])\nfor config in configs:\n    config.calc = EMT()\n    config.set_constraint(constraint)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we make the NEB object, and call its interpolate() method to guess\nthe intermediate steps\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "band = NEB(configs)\nband.interpolate()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then we relax the configurations in the NEB object. We can use same FIRE\noptimization class we might apply to optimize an Atoms object; NEB presents\nappropriate degrees of freedom and gradients for this to optimise multiple\nconfigurations simultaneously.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "relax = FIRE(band)\nrelax.run()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, we compare intermediate steps to the initial energy\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "energy_initial = initial.get_potential_energy()\nn2_distances = []\nenergy_differences = []\nfor i, config in enumerate(configs):\n    n2_distance = np.linalg.norm(config[-2].position - config[-1].position)\n    delta_energy = config.get_potential_energy() - energy_initial\n    n2_distances.append(n2_distance)\n    energy_differences.append(delta_energy)\n    print(\n        f'{i:>3}\\td = {n2_distance:>4.3f}AA '\n        f'\\tdelta energy = {delta_energy:>5.3f}eV'\n    )"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "After the calculation is complete, the energy difference\nwith respect to the initial state is given for each image,\nas well as the distance between the N atoms.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "energy_differences = [e * 1.0e03 for e in energy_differences]  # in meV\n\nfig, axs = plt.subplots(\n    nrows=2,\n    gridspec_kw={'height_ratios': [0.5, 1.0]},\n)\n\naxs[0].plot(n2_distances, energy_differences)\nscat = axs[0].scatter(n2_distances[0], energy_differences[0], s=10.0**2.0)\naxs[0].set_ylabel(r'$E(i) - E_{\\mathrm{initial}}$ (meV)')\naxs[0].set_xlabel(r'$d_{\\mathrm{N}-\\mathrm{N}} (\\AA)$')\n\n# Plot the atomic structure, focussing on the two Nitrogen atoms\nplot_atoms(config, axs[1], rotation='0x', show_unit_cell=0)\naxs[1].set_axis_off()\naxs[1].set_ylim(0.0, 5.642)\naxs[1].set_xlim(0.0, 14.833)\n\n\ndef animate(i):\n    scat.set_offsets((n2_distances[i], energy_differences[i]))\n\n    # Remove the previous atomic plot\n    [p.remove() for p in axs[1].patches]\n    plot_atoms(configs[i], axs[1], rotation='0x', show_unit_cell=0)\n    axs[1].set_xlim(0.0, 14.833)\n    axs[1].set_ylim(0.0, 5.642)\n    axs[1].set_axis_off()\n    return (scat,)\n\n\nani = animation.FuncAnimation(\n    fig, animate, repeat=True, frames=len(configs) - 1, interval=200\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
}