{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Constrained Calculations - Surface diffusion energy barriers\n\nIn this tutorial, we will calculate the energy barrier that was found\nusing the :mod:`NEB <ase.mep.neb>` method in the `diffusion tutorial`\ntutorial.  Here, we use a simple :class:`~ase.constraints.FixedPlane`\nconstraint that forces the Au atom to relax in the *yz*-plane only:\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from ase.build import add_adsorbate, fcc100\nfrom ase.calculators.emt import EMT\nfrom ase.constraints import FixAtoms, FixedPlane\nfrom ase.optimize import QuasiNewton\n\n# 2x2-Al(001) surface with 3 layers and an\n# Au atom adsorbed in a hollow site:\nslab = fcc100('Al', size=(2, 2, 3))\nadd_adsorbate(slab, 'Au', 1.7, 'hollow')\nslab.center(axis=2, vacuum=4.0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can visualize the structure with ase visualize:\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, (ax1, ax2) = plt.subplots(1, 2)\nplot_atoms(slab, ax1)\nplot_atoms(slab, ax2, rotation='-90x')\nax1.set_title('top view')\nax2.set_title('side view')\nax1.set_axis_off()\nax2.set_axis_off()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Alternatively, you can use also use view directly:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "#   $ from ase.visualize import view\n#   $ view(slab)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now continue fixing the atoms in the slab:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Fix second and third layers:\nmask = [atom.tag > 1 for atom in slab]\n# print(mask)\nfixlayers = FixAtoms(mask=mask)\n\n# Constrain the last atom (Au atom) to move only in the yz-plane:\nplane = FixedPlane(-1, (1, 0, 0))\n\nslab.set_constraint([fixlayers, plane])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we can perform the calculation optimizing the displacement\nof the gold atom along the x-axis.\nWe do structure optimization here using the EMT potential:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Use EMT potential:\nslab.calc = EMT()\n\nfor i in range(5):\n    qn = QuasiNewton(slab, trajectory=f'mep{i}.traj')\n    qn.run(fmax=0.05)\n    # Move gold atom along x-axis:\n    slab[-1].x += slab.get_cell()[0, 0] / 8\n\n# Let's visualize the saved trajectory.\n# Here is code to visualize\n# a side-view of the path (unit cell repeated twice):\n\n\nfrom ase.io import read\n\nconfigs = [read(f'mep{i}.traj', '-1') for i in range(5)]\n\n# for easier visualization, let's repeat the structures\nconfigs_repeated = [config.repeat((2, 1, 1)) for config in configs]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can visualize the structures with `ase.visualize.plot.animate`:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from ase.visualize.plot import animate\n\nanimate(\n    configs_repeated,\n    ax=None,\n    interval=500,  # in ms; same default value as in FuncAnimation\n    rotation=('-90x,0y,0z'),\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's plot the energy and look at the barrier.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# get the potential energies of the structures\nenergies = [config.get_potential_energy() for config in configs]\n# set last energy value to 0 for easier comparison\nenergies = [energy - energies[-1] for energy in energies]\nplt.ylabel(r'$E(i) - E_{\\mathrm{final}}$ (meV)')\nplt.xlabel('Image number')\nplt.plot(range(1, len(energies) + 1), energies)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The barrier is found to\nbe 0.35 eV - exactly as in the `NEB <diffusion tutorial>`\ntutorial.\n\nThe result can also be analysed with the\ncommand :command:`ase gui mep?.traj -n\n-1` (choose :menuselection:`Tools --> NEB`).\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. seealso::\n\n  * :mod:`ase.mep.neb`\n  * :mod:`ase.constraints`\n  * `diffusion tutorial`\n  * :func:`~ase.build.fcc100`\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
}