{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# NEB and Dimer method for Self-diffusion on the Al(110) surface\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In this exercise, we will find minimum-energy paths and transition states\nusing the :mod:`Nudged Elastic Band <ase.mep.neb>` method. We will illustrate\nhow NEB can be used in ASE to compute and compare three different\ndiffusion pathways for an Al atom\non a Al(110) surface. Finally, another method\nfor finding the transition state (i.e. the highest-energy state), the Dimer\nmethod, will also be explored.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Initialize the system\n\nAl(110) surface can be generated with ASE code\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from math import sqrt\n\nimport numpy as np\n\nfrom ase import Atom, Atoms\n\na = 4.0614\nb = a / sqrt(2)\nh = b / 2"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Create  :class:`~ase.Atoms` object and\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "initial = Atoms(\n    'Al2',\n    positions=[(0, 0, 0), (a / 2, b / 2, -h)],\n    cell=(a, b, 2 * h),\n    pbc=(1, 1, 0),\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Multiply the unit cell to make it larger in x,y,z\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "initial *= (4, 4, 2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "You can visualize the surface using\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=('-60x, 10y,0z'))\nax.set_axis_off()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then, lets add the Al addatom which will be the one moving on the surface.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "initial.append(Atom('Al', (a / 2, b / 2, 3 * h)))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Center the cell in vacuum along the z axis\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "initial.center(vacuum=4.0, axis=2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Visualize the new atom in the cell\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots()\nplot_atoms(initial, ax, rotation=('-60x, 10y,0z'))\nax.set_axis_off()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Perform a NEB calculation\n\nThe adatom can jump along the rows (into the picture) or across\nthe rows (to the right inthe picture).\nWe are going to compute this motion to find out which of the two jump\nwill have the largest energy barrier.\nTo do this, you need to create an image with the atoms at\ntheir final position.\nFirst copy the initial :class:`~ase.Atoms` object\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "final = initial.copy()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then move the last atom of the :class:`~ase.Atoms` object \"final\"\n(the one atom we just added before) of +b along the second positional array\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "final.positions[-1, 1] += b"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Visualize the new atom in the cell\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots()\nplot_atoms(final, ax, rotation=('-60x, 10y,0z'))\nax.set_axis_off()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let us fix the atoms that are not moving by creating a constraint\nand setting this constraint to the images.\nTo do this, we create a mask of boolean array that select\nfixed atoms (the two bottom layers):\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from ase.calculators.emt import EMT\nfrom ase.constraints import FixAtoms\n\nmask = initial.positions[:, 2] - min(initial.positions[:, 2]) < 1.5 * h\nconstraint = FixAtoms(mask=mask)\nprint(mask)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Set the :class:`~ase.constraints.FixAtoms` to the :class:`~ase.Atoms` objects,\nand in the same loop, set the calculator (in this example we use EMT,\nbut you can use any calculator supported by ASE)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "initial.calc = EMT()\ninitial.set_constraint(constraint)\nfinal.calc = EMT()\nfinal.set_constraint(constraint)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Use :class:`~ase.calculators.emt` calculator and\n:class:`~ase.optimize.QuasiNewton` Algorithm to optimize\nthe geometry of the initial and final states\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from ase.optimize import QuasiNewton\n\nQuasiNewton(initial).run(fmax=0.05)\nQuasiNewton(final).run(fmax=0.05)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then, construct a list of images by copying the first image several time in\nan array and append to this list the final image\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "images = [initial]\nfor i in range(5):\n    images.append(initial.copy())\nimages.append(final)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Because the .copy() method does not copy the calculator, you need to set a\nnew one for the created images\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for image in images:\n    image.calc = EMT()\n    image.set_constraint(constraint)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Create a Nudged Elastic Band (:class:`~ase.mep import NEB`) object\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from ase.mep import NEB\n\nneb = NEB(images)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Make a starting guess for the minimum energy path by performing a linear\ninterpolation from the initial to the final image\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "neb.interpolate()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Perform the NEB calculation minimizing the force below 0.05 eV/A\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from ase.optimize import MDMin\n\nminimizer = MDMin(neb)\nminimizer.run(fmax=0.05)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Visualize the minimum energy path (MEP) in side view to see the motion.\nHere, we look at the surface slab in yz direction.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for image in images:\n    fig, ax = plt.subplots()\n    plot_atoms(image, ax, rotation=('-90x, 90y,0z'))\n    ax.set_axis_off()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the variation of potential energy\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "potential_energies = [image.get_potential_energy() for image in images]\n\nfig, ax = plt.subplots()\nplt.plot(\n    range(len(potential_energies)),\n    potential_energies - potential_energies[0],\n    marker='+',\n)\nplt.xlabel('Image number')\nplt.ylabel('Potential energy (eV)')\n\ndiff = np.max(potential_energies) - potential_energies[0]\nprint(f'The energy barrier is {diff:.4f} eV.')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "You can visualize the NEB path using ASE GUI after saving the\ntajectory in a file\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from ase.io import write\n\nwrite('neb_path.traj', images, format='traj')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Otherwise, you can use ``ase gui neb_path.traj`` command in your terminal and\nvisualize the energy curve by plotting ``i, E[i] - E[1]``.\nYou now can answer those questions :\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "* How is the shape of the potential (symmetric/asymmetric) and does this make\n  sense for this process (when looking at the moving adatom in the\n  simulation)?\n* What is the energy barrier?\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Beyond your first NEB calculation\n\nYou now can redo the same process to find the energy barrier\nto cross one row.\nThe following code will produce the result (by making use of the previously\ninitialized code), though we encourage you to try by yourself.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "final = initial.copy()\nfinal.positions[-1, 0] += a"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the images\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for image in [initial, final]:\n    fig, ax = plt.subplots()\n    plot_atoms(image, ax, rotation=('-90x, 0y, 0z'))\n    ax.set_axis_off()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Construct a list of images:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "images = [initial]\nfor i in range(5):\n    images.append(initial.copy())\nimages.append(final)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Make a mask of zeros and ones that select fixed atoms (the\ntwo bottom layers):\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mask = initial.positions[:, 2] - min(initial.positions[:, 2]) < 1.5 * h\nconstraint = FixAtoms(mask=mask)\nprint(mask)\n\nfor image in images:\n    # Let all images use an EMT calculator:\n    image.calc = EMT()\n    image.set_constraint(constraint)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Relax the initial and final states:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "QuasiNewton(initial).run(fmax=0.05)\nQuasiNewton(final).run(fmax=0.05)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Create a Nudged Elastic Band:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "neb = NEB(images)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Make a starting guess for the minimum energy path (a straight line\nfrom the initial to the final state):\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "neb.interpolate()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Relax the NEB path:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "minimizer = MDMin(neb)\nminimizer.run(fmax=0.05)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Visualize the MEP in side view to see the motion.\nWe are now viewing the xz direction of the cell.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for image in images:\n    fig, ax = plt.subplots()\n    plot_atoms(image, ax, rotation=('-90x, 0y, 0z'))\n    ax.set_axis_off()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the variation of potential energy\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "potential_energies = [image.get_potential_energy() for image in images]\n\nfig, ax = plt.subplots()\nplt.plot(\n    range(len(potential_energies)),\n    potential_energies - potential_energies[0],\n    marker='+',\n)\nplt.xlabel('Image number')\nplt.ylabel('Potential energy (eV)')\n\ndiff = np.max(potential_energies) - potential_energies[0]\nprint(f'The energy barrier is {diff:.4f} eV.')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Finding the third mechanism\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "A third diffusion process can be found: Diffusion by an exchange process.\nYou can read more about it in the paper listed :mod:`here <ase.mep.dimer>`.\nFind the barrier for this process, and compare the energy barrier\nwith the two other ones.\nThe following code will produce the result (by making use of the previously\ninitialized code), though we encourage you to try by yourself.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "a = 4.0614\nb = a / sqrt(2)\nh = b / 2\ninitial = Atoms(\n    'Al2',\n    positions=[(0, 0, 0), (a / 2, b / 2, -h)],\n    cell=(a, b, 2 * h),\n    pbc=(1, 1, 0),\n)\ninitial *= (2, 2, 2)\ninitial.append(Atom('Al', (a / 2, b / 2, 3 * h)))\ninitial.center(vacuum=4.0, axis=2)\n\n\nfinal = initial.copy()\n# move adatom to row atom 14\nfinal.positions[-1, :] = initial.positions[14]\n# Move row atom 14 to the next row\nfinal.positions[14, :] = initial.positions[-1] + [a, b, 0]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Visualize the initial and final images\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for image in [initial, final]:\n    fig, ax = plt.subplots()\n    plot_atoms(image, ax, rotation=('-60x, 10y,0z'))\n    ax.set_axis_off()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Construct a list of images:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "images = [initial]\nfor i in range(5):\n    images.append(initial.copy())\nimages.append(final)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Make a mask of zeros and ones that select fixed atoms (the\ntwo bottom layers):\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mask = initial.positions[:, 2] - min(initial.positions[:, 2]) < 1.5 * h\nconstraint = FixAtoms(mask=mask)\nprint(mask)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let all images use an EMT calculator:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for image in images:\n    image.calc = EMT()\n    image.set_constraint(constraint)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Relax the initial and final states:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "QuasiNewton(initial).run(fmax=0.05)\nQuasiNewton(final).run(fmax=0.05)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Create a Nudged Elastic Band:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "neb = NEB(images)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Make a starting guess for the minimum energy path (a straight line\nfrom the initial to the final state):\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "neb.interpolate()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Relax the NEB path:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "minimizer = MDMin(neb)\nminimizer.run(fmax=0.05)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Visualize the MEP in side view to see the motion\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for image in images:\n    fig, ax = plt.subplots()\n    plot_atoms(image, ax, rotation=('-60x, 10y,0z'))\n    ax.set_axis_off()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the variation of potential energy\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "potential_energies = [image.get_potential_energy() for image in images]\n\nfig, ax = plt.subplots()\nplt.plot(\n    range(len(potential_energies)),\n    potential_energies - potential_energies[0],\n    marker='+',\n)\nplt.xlabel('Image number')\nplt.ylabel('Potential energy (eV)')\n\ndiff = np.max(potential_energies) - potential_energies[0]\nprint(f'The energy barrier is {diff:.4f} eV.')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. hint::\n  When opening a trajectory with :program:`ase gui` with calculated energies,\n  the default plot window shows the energy versus frame number.\n  To get a better feel of the energy barrier in an NEB calculation;\n  choose  :menuselection:`Tools --> NEB`.\n  This will give a smooth curve   of the energy as a function of the NEB path\n  length, with the slope at each point estimated from the force.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Performing Dimer-method calculation\n\nIn the NEB calculations above we knew the final states, so all we had\nto do was to calculate the path between the initial state\nand the final state.\nBut in some cases we do not know the final state.\nThen the :mod:`Dimer method <ase.mep.dimer>` can be used to\nfind the transition state.\nThe result of a Dimer calculation will hence not be the complete particle\ntrajectory as in the NEB output, but rather the configuration of\nthe transition-state image.\n\nThe following code will find the transition-state image of the\njump along the row.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from ase.io import Trajectory\nfrom ase.mep import DimerControl, MinModeAtoms, MinModeTranslate\n\na = 4.0614\nb = a / sqrt(2)\nh = b / 2\ninitial = Atoms(\n    'Al2',\n    positions=[(0, 0, 0), (a / 2, b / 2, -h)],\n    cell=(a, b, 2 * h),\n    pbc=(1, 1, 0),\n)\ninitial *= (2, 2, 2)\ninitial.append(Atom('Al', (a / 2, b / 2, 3 * h)))\ninitial.center(vacuum=4.0, axis=2)\n\n\ninitial_copy = initial.copy()\n\nN = len(initial)  # number of atoms\n\n\n# Make a mask of zeros and ones that select fixed atoms - the two\n# bottom layers:\nmask = initial.positions[:, 2] - min(initial.positions[:, 2]) < 1.5 * h\nconstraint = FixAtoms(mask=mask)\ninitial.set_constraint(constraint)\n\n# Calculate using EMT:\ninitial.calc = EMT()\n\n# Relax the initial state:\nQuasiNewton(initial).run(fmax=0.05)\ne0 = initial.get_potential_energy()\n\n# To save the trajectory file\ntraj = Trajectory('dimer_along.traj', 'w', initial)\ntraj.write()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Making dimer mask list:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "d_mask = [False] * (N - 1) + [True]\n\n# Set up the dimer:\nd_control = DimerControl(\n    initial_eigenmode_method='displacement',\n    displacement_method='vector',\n    logfile=None,\n    mask=d_mask,\n)\nd_atoms = MinModeAtoms(initial, d_control)\n\n# Displacement settings:\ndisplacement_vector = np.zeros((N, 3))\n# Strength of displacement along y axis = along row:\ndisplacement_vector[-1, 1] = 0.001\n# The direction of the displacement is set by the a in\n# displacement_vector[-1, a], where a can be 0 for x, 1 for y and 2 for z.\nd_atoms.displace(displacement_vector=displacement_vector)\n\n# Converge to a saddle point:\ndim_rlx = MinModeTranslate(d_atoms, trajectory=traj, logfile=None)\ndim_rlx.run(fmax=0.001)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Visualize the Initial state and the saddle point in side view to see\nthe change\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for image in [initial, initial_copy]:\n    fig, ax = plt.subplots()\n    plot_atoms(image, ax, rotation=('-90x, 90y,0z'))\n    ax.set_axis_off()\n\n\ndiff = initial.get_potential_energy() - e0\nprint(f'The energy barrier is {diff:.4f} eV.')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "* Compare the transition-state images of the NEB and Dimer as viewed\n  in the GUI. Are they identical?\n* What is the energy barrier? How does it compare to the one found in the\n  NEB calculation?\n* Do the same as above for the jump across the row and the exchange\n  process by copying and modifying the Dimer script, while remembering\n  that you have to  give the relevant atoms a kick in a meaningful direction.\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
}