{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Extra exercise - vibrational energy\n",
        "===================================\n",
        "\n",
        "The energy calculated with DFT is the electronic energy at 0K. However, to\n",
        "model catalytic reactions we are usually intereseted in the energy landscape\n",
        "at finite temperature. In this exercise we will calculate the energy\n",
        "contributions from vibrations and see how they affect the splitting of\n",
        "N<sub>2</sub> on Ru.\n",
        "\n",
        "We calculate the vibrational energy in the harmonic approximation using the\n",
        "finite displacement method. For further reading see for example:\n",
        "\n",
        "* [Stoffel et al, Angewandte Chemie Int. Ed., 49, 5242-5266 (2010)][1]\n",
        "* [Togo and Tanaka, Scripta Materialia 108, 1-5 (2015)][2]\n",
        "\n",
        "[1]: https://onlinelibrary.wiley.com/doi/abs/10.1002/anie.200906780\n",
        "[2]: https://www.sciencedirect.com/science/article/pii/S1359646215003127\n",
        "\n",
        "\n",
        "Vibrational energy of the initial and final states\n",
        "--------------------------------------------------\n",
        "\n",
        "a) Calculating vibrations requires tighter convergence than normal energy\n",
        "   calculations. Therefore you should first take your already optimised initial\n",
        "   and final state geometries from the NEB calculations and relax them further\n",
        "   to `fmax=0.01eV/\u00c5` with the QuasiNewton optimiser and an energy cutoff of\n",
        "   500eV. Converge the eigenstates to 1e-8. (Note that for other systems you\n",
        "   might need even tighter convergence!)\n",
        "\n",
        "Submit the structures to the queue. The optimisation should take 10-15 mins\n",
        "for each structure on 8 cores.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "b) Once you have done this you can calculate the vibrations using the\n",
        "   [vibrations module in ASE][3] following\n",
        "   the template below. We only calculate the vibrations of the adsorbate and\n",
        "   assume that the frequencies of the substrate are unchanged - this is a\n",
        "   common assumption. Use 4 displacements to fit the frequencies and the same\n",
        "   calculator parameters as in a).\n",
        "\n",
        "[3]: https://wiki.fysik.dtu.dk/ase/ase/vibrations/vibrations.html\n",
        "\n",
        "Submit the calculations for the initial and the final state to the queue. It\n",
        "will take a while to run, but you can start preparing your analysis (part c\n",
        "and d) or read some of the references in the meantime.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# You can't run this cell - use it as a starting point for your python script!\n",
        "from ase.io import read\n",
        "from ase.vibrations import Vibrations\n",
        "from gpaw import GPAW, PW\n",
        "\n",
        "slab = read('your_structure')\n",
        "calc = GPAW(xc='PBE',\n",
        "            mode=PW(500),\n",
        "            kpts=?,\n",
        "            symmetry={'point_group': False},\n",
        "            txt='vib.txt')\n",
        "slab.calc = calc\n",
        "Uini = slab.get_potential_energy()\n",
        "\n",
        "vib = Vibrations(slab,\n",
        "                 name='vib',\n",
        "                 indices=[?, ?],\n",
        "                 nfree=4)\n",
        "vib.run()\n",
        "vib.summary(log='vib_summary.log')\n",
        "for i in range(6):\n",
        "    vib.write_mode(i)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "from ase.io import read\n",
        "from ase.vibrations import Vibrations\n",
        "from gpaw import GPAW, PW\n",
        "\n",
        "slab = read('tight-2Nads.traj')\n",
        "calc = GPAW(xc='PBE',\n",
        "            mode=PW(500),\n",
        "            kpts=?,\n",
        "            symmetry={'point_group': False},\n",
        "            txt='vib2.txt')\n",
        "slab.calc = calc\n",
        "Ufin = slab.get_potential_energy()\n",
        "\n",
        "vib2 = Vibrations(slab,\n",
        "                  name='vib2',\n",
        "                  indices=[?, ?],\n",
        "                  nfree=4)\n",
        "vib2.run()\n",
        "vib2.summary(log='vib2_summary.log')\n",
        "for i in range(6):\n",
        "    vib2.write_mode(i)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "\n",
        "The final three lines write out the frequencies and `.traj` files with\n",
        "animations of the different modes in order of their energy. Take a look at\n",
        "the vibrational modes in the files. Do you understand why the first mode has\n",
        "a low energy, while the last one has a high energy?\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "c) Analyse the frequencies in the harmonic approximation:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "from ase.thermochemistry import HarmonicThermo\n",
        "\n",
        "T = 300  # Kelvin\n",
        "# An example only - insert your calculated energy levels here - in eV:\n",
        "energies = [0.01, 0.05, 0.10]\n",
        "vibs = HarmonicThermo(energies)\n",
        "Efree = vibs.get_helmholtz_energy(T, verbose=True)\n",
        "print('Free energy at 300 K: ', Efree)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "The `verbose` keyword gives a detailed description of the different\n",
        "contributions to the free energy. For more information on what the different\n",
        "contributions are see the [ASE background webpage][4]\n",
        "(go to the **Harmonic limit** sub-heading).\n",
        "\n",
        "[4]: https://wiki.fysik.dtu.dk/ase/ase/thermochemistry/thermochemistry.html\n",
        "\n",
        "Now try to calculate how the different contributions change with temperature.\n",
        "You can for example make a `for` loop and use the `get_entropy()` and\n",
        "`get_internal_energy()` methods [(see description here)][5].\n",
        "\n",
        "[5]: https://wiki.fysik.dtu.dk/ase/ase/thermochemistry/thermochemistry.html#ase.thermochemistry.IdealGasThermo.get_enthalpy\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "d) Calculate how the vibrational energy affects the overall reaction energy.\n",
        "   How does it change with temperature? Which contribution is important for the\n",
        "   change in reaction energy?\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "e) To make sure that your NEB is converged you should also calculate the\n",
        "   vibrational energy of the transition state. Again, this requires tighter\n",
        "   convergence than we have used in the NEB exercise. This takes a while to run\n",
        "   so to save time, we provide the transition state geometry from a reasonably\n",
        "   converged NEB (i.e. `fmax=0.01`, a cutoff energy of 500eV and eigenstates\n",
        "   converged to 1e-8) in the file `TS.xyz`. Calculate the vibrations with these\n",
        "   parameters. How many imaginary modes do you get and how do they look? What\n",
        "   does this mean?\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [],
      "outputs": [],
      "execution_count": null
    },
    {
      "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
}