{
  "cells": [
    {
      "cell_type": "code",
      "metadata": {},
      "source": [],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Convergence checks\n",
        "==================\n",
        "\n",
        "In this notebook we look at the adsorption energy and height of a nitrogen\n",
        "atom on a Ru(0001) surface in the hcp site.  We check for convergence with\n",
        "respect to:\n",
        "\n",
        "* number of layers\n",
        "* number of k-points in the BZ\n",
        "* plane-wave cutoff energy\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Nitrogen atom\n",
        "-------------\n",
        "\n",
        "First step is an isolated nitrogen atom which has a magnetic moment of 3.\n",
        "More information:\n",
        "[Atoms][1] and [GPAW parameters][2].\n",
        "\n",
        "[1]: https://wiki.fysik.dtu.dk/ase/ase/atoms.html#ase.Atoms\n",
        "[2]: https://wiki.fysik.dtu.dk/gpaw/documentation/basic.html#parameters\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "ecut = 400.0\n",
        "vacuum = 4.0\n",
        "from ase import Atoms\n",
        "from gpaw import GPAW, PW, Davidson\n",
        "nitrogen = Atoms('N', magmoms=[3])\n",
        "nitrogen.center(vacuum=4.0)\n",
        "nitrogen.calc = GPAW(txt='N.txt',\n",
        "                     mode=PW(ecut),\n",
        "                     eigensolver=Davidson(niter=2))\n",
        "en = nitrogen.get_potential_energy()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "print(en, 'eV')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Clean slab\n",
        "----------\n",
        "\n",
        "We use the [ase.build.hcp0001()][3] function to build the Ru(0001) surface.\n",
        "\n",
        "[3]: https://wiki.fysik.dtu.dk/ase/ase/build/surface.html#ase.build.hcp0001\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "nlayers = 2\n",
        "a = 2.72\n",
        "c = 1.58 * a\n",
        "from ase.build import hcp0001\n",
        "slab = hcp0001('Ru', a=a, c=c, size=(1, 1, nlayers), vacuum=vacuum)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "from ase.visualize import view\n",
        "view(slab, repeat=(3, 3, 2))"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "nkpts = 7\n",
        "slab.calc = GPAW(txt='Ru.txt',\n",
        "                 mode=PW(ecut),\n",
        "                 eigensolver=Davidson(niter=2),\n",
        "                 kpts={'size': (nkpts, nkpts, 1), 'gamma': True},\n",
        "                 xc='PBE')\n",
        "eru = slab.get_potential_energy()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "# N/Ru(0001):\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "import numpy as np\n",
        "height = 1.1\n",
        "nslab = hcp0001('Ru', a=a, c=c, size=(1, 1, nlayers))\n",
        "# Calculate the coordianates of the N-atoms:\n",
        "z = slab.positions[:, 2].max() + height\n",
        "x, y = np.dot([2 / 3, 2 / 3], slab.cell[:2, :2])\n",
        "nslab.append('N')\n",
        "nslab.positions[-1] = [x, y, z]\n",
        "nslab.center(vacuum=vacuum, axis=2)  # 2: z-axis"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Alternatively, you can just use the [add_adsorbate()][4] function:\n",
        "\n",
        "[4]: https://wiki.fysik.dtu.dk/ase/ase/build/surface.html#ase.build.add_adsorbate\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "height = 1.1\n",
        "nslab = hcp0001('Ru', a=a, c=c, size=(1, 1, nlayers))\n",
        "from ase.build import add_adsorbate\n",
        "add_adsorbate(nslab, 'N', position='hcp', height=height)\n",
        "nslab.center(vacuum=vacuum, axis=2)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "view(nslab, repeat=(3, 3, 2))"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "from ase.io import write\n",
        "write('nru2.png', nslab.repeat((3, 3, 1)))"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "![rnu](nru2.png)\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "nslab.calc = GPAW(txt='NRu.txt',\n",
        "                  mode=PW(ecut),\n",
        "                  eigensolver=Davidson(niter=2),\n",
        "                  poissonsolver={'dipolelayer': 'xy'},\n",
        "                  kpts={'size': (nkpts, nkpts, 1), 'gamma': True},\n",
        "                  xc='PBE')\n",
        "enru0 = nslab.get_potential_energy()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "print('Unrelaxed adsoption energy:', enru0 - eru - en, 'eV')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "nslab.get_forces()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "The force on the N-atom is quite big.  Let's freeze the surface and relax the\n",
        "adsorbate.  We use\n",
        "[ase.optimize.BFGSLineSearch][5] and\n",
        "[ase.constraints.FixAtoms][6]\n",
        "for this task.\n",
        "\n",
        "[5]: https://wiki.fysik.dtu.dk/ase/ase/optimize.html#module-ase.optimize\n",
        "[6]: https://wiki.fysik.dtu.dk/ase/ase/constraints.html#ase.constraints.FixAtoms\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# This cell will take a few minutes to finish ...\n",
        "from ase.constraints import FixAtoms\n",
        "from ase.optimize import BFGSLineSearch\n",
        "nslab.constraints = FixAtoms(indices=list(range(nlayers)))\n",
        "optimizer = BFGSLineSearch(nslab, trajectory='NRu.traj')\n",
        "optimizer.run(fmax=0.01)\n",
        "height = nslab.positions[-1, 2] - nslab.positions[:-1, 2].max()\n",
        "print('Height:', height, 'Ang')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "enru = nslab.get_potential_energy()\n",
        "print('Relaxed adsorption energy:', enru - eru - en, 'eV')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "In order to make it easy to check for convergence of the adsorption energy\n",
        "and height we write a little function that does all of the stuff above taking\n",
        "`nlayers`, `nkpts` and `ecut` as input parameters.\n",
        "\n",
        "The `adsorb()` function is shown below for completenes, but you should not\n",
        "use it inside this notebook.  Instead, please take a look at the\n",
        "[check_convergence.py][7]\n",
        "script that also contains the definition of the `adsorb()` function.  The\n",
        "script will do a bunch of calculations with different parameters and store\n",
        "the results in a database file (`convergence.db`) that we analyse below ...\n",
        "\n",
        "[7]: https://gitlab.com/gpaw/gpaw/blob/master/doc/summerschools/summerschool18/catalysis/check_convergence.py\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "def adsorb(db, height=1.2, nlayers=3, nkpts=7, ecut=400):\n",
        "    \"\"\"Adsorb nitrogen in hcp-site on Ru(0001) surface.\n",
        "\n",
        "    Do calculations for N/Ru(0001), Ru(0001) and a nitrogen atom\n",
        "    if they have not already been done.\n",
        "\n",
        "    db: Database\n",
        "        Database for collecting results.\n",
        "    height: float\n",
        "        Height of N-atom above top Ru-layer.\n",
        "    nlayers: int\n",
        "        Number of Ru-layers.\n",
        "    nkpts: int\n",
        "        Use a (nkpts * nkpts) Monkhorst-Pack grid that includes the\n",
        "        Gamma point.\n",
        "    ecut: float\n",
        "        Cutoff energy for plane waves.\n",
        "\n",
        "    Returns height.\n",
        "    \"\"\"\n",
        "\n",
        "    name = f'Ru{nlayers}-{nkpts}x{nkpts}-{ecut:.0f}'\n",
        "\n",
        "    parameters = dict(mode=PW(ecut),\n",
        "                      eigensolver=Davidson(niter=2),\n",
        "                      poissonsolver={'dipolelayer': 'xy'},\n",
        "                      kpts={'size': (nkpts, nkpts, 1), 'gamma': True},\n",
        "                      xc='PBE')\n",
        "\n",
        "    # N/Ru(0001):\n",
        "    slab = hcp0001('Ru', a=a, c=c, size=(1, 1, nlayers))\n",
        "    z = slab.positions[:, 2].max() + height\n",
        "    x, y = np.dot([2 / 3, 2 / 3], slab.cell[:2, :2])\n",
        "    slab.append('N')\n",
        "    slab.positions[-1] = [x, y, z]\n",
        "    slab.center(vacuum=vacuum, axis=2)  # 2: z-axis\n",
        "\n",
        "    # Fix first nlayer atoms:\n",
        "    slab.constraints = FixAtoms(indices=list(range(nlayers)))\n",
        "\n",
        "    id = db.reserve(name=f'N/{nlayers}Ru(0001)', nkpts=nkpts, ecut=ecut)\n",
        "    if id is not None:  # skip calculation if already done\n",
        "        slab.calc = GPAW(txt='N' + name + '.txt',\n",
        "                         **parameters)\n",
        "        optimizer = BFGSLineSearch(slab, logfile='N' + name + '.opt')\n",
        "        optimizer.run(fmax=0.01)\n",
        "        height = slab.positions[-1, 2] - slab.positions[:-1, 2].max()\n",
        "        db.write(slab, id=id,\n",
        "                 name=f'N/{nlayers}Ru(0001)', nkpts=nkpts, ecut=ecut,\n",
        "                 height=height)\n",
        "\n",
        "    # Clean surface (single point calculation):\n",
        "    id = db.reserve(name=f'{nlayers}Ru(0001)', nkpts=nkpts, ecut=ecut)\n",
        "    if id is not None:\n",
        "        del slab[-1]  # remove nitrogen atom\n",
        "        slab.calc = GPAW(txt=name + '.txt',\n",
        "                         **parameters)\n",
        "        slab.get_forces()\n",
        "        db.write(slab, id=id,\n",
        "                 name=f'{nlayers}Ru(0001)', nkpts=nkpts, ecut=ecut)\n",
        "\n",
        "    # Nitrogen atom:\n",
        "    id = db.reserve(name='N-atom', ecut=ecut)\n",
        "    if id is not None:\n",
        "        # Create spin-polarized nitrogen atom:\n",
        "        molecule = Atoms('N', magmoms=[3])\n",
        "        molecule.center(vacuum=4.0)\n",
        "        # Remove parameters that make no sense for an isolated atom:\n",
        "        del parameters['kpts']\n",
        "        del parameters['poissonsolver']\n",
        "        # Calculate energy:\n",
        "        molecule.calc = GPAW(txt=name + '.txt', **parameters)\n",
        "        molecule.get_potential_energy()\n",
        "        db.write(molecule, id=id, name='N-atom', ecut=ecut)\n",
        "\n",
        "    return height\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Read more about ASE databases\n",
        "[here](https://wiki.fysik.dtu.dk/ase/ase/db/db.html#module-ase.db).\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "from ase.db import connect\n",
        "db = connect('convergence.db')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "assert len(db) == 72\n",
        "# We assume that the calculations have already been done for you"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# This cell should not take any time because the adsorb() function\n",
        "# is clever enough to skip calculations already in the database.\n",
        "h = 1.2\n",
        "for n in range(1, 10):  # layer\n",
        "    h = adsorb(db, h, n, 7, 400)\n",
        "for k in range(4, 18):  # k-points\n",
        "    h = adsorb(db, h, 2, k, 400)\n",
        "for ecut in range(350, 801, 50):  # plane-wave cutoff\n",
        "    h = adsorb(db, h, 2, 7, ecut)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You can inspect database file with the\n",
        "[command line tool](https://wiki.fysik.dtu.dk/ase/ase/db/db.html#ase-db)\n",
        "`ase db` like this:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "!ase db convergence.db -c ++ -L 0  # show all columns (-c ++); show all rows (-L 0)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "!ase db convergence.db formula=Ru2N,nkpts=7 -c ecut,height -s ecut"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Now we can analyse the results of the convergence tests.  We extract the\n",
        "result from the database with a little helper function `select()`:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "def select(nlayers, nkpts, ecut):\n",
        "    \"\"\"Extract adsorption energy and height from database.\"\"\"\n",
        "    en = db.get(N=1, Ru=0, ecut=ecut).energy\n",
        "    eru = db.get(N=0, Ru=nlayers, nkpts=nkpts, ecut=ecut).energy\n",
        "    row = db.get(N=1, Ru=nlayers, nkpts=nkpts, ecut=ecut)\n",
        "    enru = row.energy\n",
        "    return row.height, enru - eru - en\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "%matplotlib notebook\n",
        "from matplotlib import pyplot as plt"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "n = np.arange(1, 10)\n",
        "h, e = np.array([select(nlayers, 7, 400) for nlayers in n]).T\n",
        "fig, axs = plt.subplots(2, 1, sharex=True)\n",
        "fig.subplots_adjust(hspace=0)\n",
        "axs[0].plot(n, h)\n",
        "axs[1].plot(n, e)\n",
        "axs[0].set_ylabel('height [\u00c5]')\n",
        "axs[1].set_ylabel('ads. energy [eV]')\n",
        "axs[1].set_xlabel('number of layers')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "k = np.arange(4, 18)\n",
        "h, e = np.array([select(2, nkpts, 400) for nkpts in k]).T\n",
        "fig, axs = plt.subplots(2, 1, sharex=True)\n",
        "fig.subplots_adjust(hspace=0)\n",
        "axs[0].plot(k, h)\n",
        "axs[1].plot(k, e)\n",
        "axs[0].set_ylabel('height [\u00c5]')\n",
        "axs[1].set_ylabel('ads. energy [eV]')\n",
        "axs[1].set_xlabel('number of k-points')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "x = np.arange(350, 801, 50)\n",
        "h, e = np.array([select(2, 7, ecut) for ecut in x]).T\n",
        "fig, axs = plt.subplots(2, 1, sharex=True)\n",
        "fig.subplots_adjust(hspace=0)\n",
        "axs[0].plot(x, h)\n",
        "axs[1].plot(x, e)\n",
        "axs[0].set_ylabel('height [\u00c5]')\n",
        "axs[1].set_ylabel('ads. energy [eV]')\n",
        "axs[1].set_xlabel('plane-wave cutoff energy [eV]')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Conclusion\n",
        "----------\n",
        "\n",
        "For accurate calculations you would need:\n",
        "\n",
        "* a plane-wave cutoff of 600 eV\n",
        "* 5 layers of Ru\n",
        "* 9x9 Monkhorst-Pack grid for BZ sampling (for a 1x1 unit cell)\n",
        "\n",
        "For our quick'n'dirty calculations we will use 350 eV, 2 layers and a\n",
        "4x4 $\\\\Gamma$-centered Monkhorst-Pack grid (for a 2x2 unit cell).\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
}