{
  "cells": [
    {
      "cell_type": "code",
      "metadata": {},
      "source": [],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Introduction to Nudged Elastic Band (NEB) calculations\n",
        "\n",
        "This tutorial describes how to use the NEB method to calculate the diffusion\n",
        "barrier for an Au atom on Al(001). If you are not familiar with the NEB\n",
        "method some relevant references are listed\n",
        "[here.](https://wiki.fysik.dtu.dk/ase/ase/neb.html)\n",
        "\n",
        "The tutorial uses the EMT potential in stead of DFT, as this is a lot faster.\n",
        "It is based on a [tutorial found on the ASE\n",
        "webpage](https://wiki.fysik.dtu.dk/ase/tutorials/neb/diffusion.html#diffusion-tutorial).\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "%matplotlib notebook\n",
        "import matplotlib.pyplot as plt\n",
        "from ase.build import fcc100, add_adsorbate\n",
        "from ase.constraints import FixAtoms\n",
        "from ase.calculators.emt import EMT\n",
        "from ase.optimize import BFGS\n",
        "from ase.visualize import view\n",
        "from ase.mep import NEB\n",
        "from ase.io import read"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "First we set up the initial state and check that it looks ok:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# 2x2-Al(001) surface with 3 layers and an\n",
        "# Au atom adsorbed in a hollow site:\n",
        "slab = fcc100('Al', size=(2, 2, 3))\n",
        "add_adsorbate(slab, 'Au', 1.7, 'hollow')\n",
        "slab.center(axis=2, vacuum=4.0)\n",
        "view(slab)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Then we optimise the structure and save it\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Fix second and third layers:\n",
        "mask = [atom.tag > 1 for atom in slab]\n",
        "slab.set_constraint(FixAtoms(mask=mask))\n",
        "\n",
        "# Use EMT potential:\n",
        "slab.calc = EMT()\n",
        "\n",
        "# Optimise initial state:\n",
        "qn = BFGS(slab, trajectory='initial.traj')\n",
        "qn.run(fmax=0.05)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "We make the final state by moving the Au atom one lattice constant and\n",
        "optimise again.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Optimise final state:\n",
        "slab[-1].x += slab.get_cell()[0, 0] / 2\n",
        "qn = BFGS(slab, trajectory='final.traj')\n",
        "qn.run(fmax=0.05)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Now we make a NEB calculation with 3 images between the inital and final\n",
        "states. The images are initially made as copies of the initial state and the\n",
        "command `interpolate()` makes a linear interpolation between the initial and\n",
        "final state. As always, we check that everything looks ok before we run the\n",
        "calculation.\n",
        "\n",
        "NOTE: The linear interpolation works well in this case but not for e.g.\n",
        "rotations. In this case an improved starting guess can be made with the [IDPP\n",
        "method.](https://wiki.fysik.dtu.dk/ase/tutorials/neb/idpp.html#idpp-tutorial)\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "initial = read('initial.traj')\n",
        "final = read('final.traj')\n",
        "\n",
        "constraint = FixAtoms(mask=[atom.tag > 1 for atom in initial])\n",
        "\n",
        "n_im = 3  # number of images\n",
        "images = [initial]\n",
        "for i in range(n_im):\n",
        "    image = initial.copy()\n",
        "    image.calc = EMT()\n",
        "    image.set_constraint(constraint)\n",
        "    images.append(image)\n",
        "\n",
        "images.append(final)\n",
        "\n",
        "neb = NEB(images)\n",
        "neb.interpolate()\n",
        "view(images)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "qn = BFGS(neb, trajectory='neb.traj')\n",
        "qn.run(fmax=0.05)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "We visualize the final path with:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "view(images)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You can find the barrier by selecting Tools->NEB in the gui (unfortunately,\n",
        "the gui cannot show graphs when started from a notebook), or you can make a\n",
        "script using\n",
        "[NEBTools](https://wiki.fysik.dtu.dk/ase/ase/neb.html#ase.neb.NEBTools),\n",
        "e.g.:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "from ase.mep import NEBTools\n",
        "\n",
        "images = read('neb.traj@-5:')\n",
        "\n",
        "nebtools = NEBTools(images)\n",
        "\n",
        "# Get the calculated barrier and the energy change of the reaction.\n",
        "Ef, dE = nebtools.get_barrier()\n",
        "print('Barrier:', Ef)\n",
        "print('Reaction energy:', dE)\n",
        "\n",
        "# Generate new matplotlib axis - otherwise nebtools plot double.\n",
        "fig, ax = plt.subplots()\n",
        "nebtools.plot_band(ax)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Exercise\n",
        "\n",
        "Now you should make your own NEB using the configuration with N<sub>2</sub>\n",
        "lying down as the initial state and the configuration with two N atoms\n",
        "adsorbed on the surface as the final state. The NEB needs to run in parallel\n",
        "so you should make it as a python script, however you can use the Notebook to\n",
        "test your configurations (but not the parallelisation) if you like and export\n",
        "it as a script in the end.\n",
        "\n",
        "### Parallelisation\n",
        "\n",
        "The NEB should be parallelised over images. An example can be found in [this\n",
        "GPAW tutorial](https://wiki.fysik.dtu.dk/gpaw/tutorialsexercises/moleculardynamics/neb/neb.html). The\n",
        "script enumerates the cpu's and uses this number (``rank``) along with the\n",
        "total number of cpu's (``size``) to distribute the tasks.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# This code is just for illustration\n",
        "from ase.parallel import world\n",
        "n_im = 4              # Number of images\n",
        "n = world.size // n_im      # number of cpu's per image\n",
        "j = 1 + world.rank // n     # image number on this cpu\n",
        "assert n_im * n == world.size"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "For each image we assign a set of cpu's identified by their rank. The rank\n",
        "numbers are given to the calculator associated with this image.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# This code is just for illustration\n",
        "from gpaw import GPAW, PW\n",
        "images = [initial]\n",
        "for i in range(n_im):\n",
        "    ranks = range(i * n, (i + 1) * n)\n",
        "    image = initial.copy()\n",
        "    image.set_constraint(constraint)\n",
        "    if world.rank in ranks:\n",
        "        calc = GPAW(mode=PW(350),\n",
        "                    nbands='130%',\n",
        "                    ...,\n",
        "                    txt=f'{i}.txt',\n",
        "                    communicator=ranks)\n",
        "        image.calc = calc\n",
        "    images.append(image)\n",
        "images.append(final)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "When running the parallel NEB, you should choose the number of CPU cores\n",
        "properly.  Let Ncore = N_im * Nk where N_im is the number of images, and Nk\n",
        "is a divisor of the number of k-points; i.e. if there are 6 irreducible\n",
        "k-point, Nk should be 1, 2, 3 or 6.  Keep the total number of cores to 24 or\n",
        "less, or your job will wait too long in the queue.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "### Input parameters\n",
        "\n",
        "Some suitable parameters for the NEB are given below:\n",
        "\n",
        "* Use the same calculator and constraints as for the initial and final images, but remember to set the `communicator` as described above\n",
        "* Use 6 images. This gives a reasonable description of the energy landscape and can be run e.g. on 12 cores.\n",
        "* Use a spring constant of 1.0 between the images. A lower value will slow the convergence\n",
        "* Relax the initial NEB until `fmax=0.1eV/\u00c5`, then switch on the climbing image and relax until `fmax=0.05eV/\u00c5`.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Once the calculation is done you should check that the final path looks\n",
        "reasonable. What is the N-N distance at the saddle point? Use NEBTools to\n",
        "calculate the barrier. Is N<sub>2</sub> likely to dissociate on the surface\n",
        "at room temperature?\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [],
      "outputs": [],
      "execution_count": null
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "CAMD2022",
      "language": "python",
      "name": "camd2022"
    },
    "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.6.1"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 1
}