{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Catalysis: Dissociative adsorbtion of N2 on a metal surface\n", "\n", "This is the rate limiting step for ammonia synthesis.\n", "\n", "**Scientific disclaimer:** These calculations are done on a flat surface. In reality, the process takes place at the foot of an atomic step on the surface. Doing calculations on this more realistic system would be too slow for these exercises. For the same reason, we use a metal slab with only two layers, a realistic calculation would require the double." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## N2 Adsorption on a metal surface\n", "This notebook shows how to calculate the adsorption energy of an N2 molecule on a closepacked Ru surface. The first cell imports some modules from the ASE and GPAW packages " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "from ase import Atoms\n", "from gpaw import GPAW, PW\n", "from ase.constraints import FixAtoms\n", "from ase.optimize import QuasiNewton\n", "from ase.build import fcc111, hcp0001\n", "import numpy as np\n", "from ase.visualize import view\n", "from ase.io import read, write\n", "import time" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Setting up the metal surface\n", "\n", "Ru crystalises in the hcp structure with a lattice constants a = 2.706 \u00c5 and c = 4.282 \u00c5. It is often better to use the lattice constants corresponding to the DFT variant used (here PBE with PAW). We get this from http://oqmd.org. \n", "\n", "We model the surface by a 2 layer slab of metal atoms, and add 5\u00c5 vacuum on each side. \n", "\n", "We visualize the system with ASE GUI, so you can check that everything looks right. This pops up a new window." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "a_Ru = 2.704 # PBE value from OQMD.org; expt value is 2.706\n", "slab = hcp0001('Ru', a=a_Ru, size=(2, 2, 2), vacuum=5.0)\n", "\n", "# Other metals are possible, for example Rhodium\n", "# Rhodium is FCC so you should use fcc111(...) to set up the system (same arguments).\n", "# Remember to set the FCC lattice constant, get it from OQMD.org.\n", "\n", "# a_Rh = 3.793\n", "# slab = fcc111('Rh', a=a_Rh, size=(2, 2, 2), vacuum=5.0)\n", "\n", "view(slab)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To optimise the slab we need a calculator. We use the GPAW calculator in plane wave (PW) mode with the PBE exchange-correlation functional. The convergence with respect to the cutoff energy and k-point sampling should always be checked - see `Convergence.ipynb`for more information on how this can be done. For this exercise an energy cutoff of 350eV and 4x4x1 k-point mesh is chosen to give reasonable results with a limited computation time." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "calc = GPAW(xc='PBE',\n", " mode=PW(350),\n", " kpts={'size': (4, 4, 1), 'gamma': True}, \n", " convergence={'eigenstates': 1e-6})\n", "slab.set_calculator(calc)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The bottom layer of the slab is fixed during optimisation. The structure is optimised until the forces on all atoms are below 0.05eV/\u00c5. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "z = slab.positions[:, 2]\n", "constraint = FixAtoms(mask=(z < z.min() + 1.0))\n", "slab.set_constraint(constraint)\n", "dyn = QuasiNewton(slab, trajectory='Ru.traj')\n", "t = time.time()\n", "dyn.run(fmax=0.05)\n", "print('Wall time: {} min.'.format((time.time() - t) / 60))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The calculation will take ca. 5 minutes. While the calculation is running you can take a look at the output. How many k-points are there in total and how many are there in the irreducible part of the Brillouin zone? What does this mean for the speed of the calculation?\n", "\n", "What are the forces and the energy after each iteration? You can read it directly in the output above, or from the saved .traj file like this:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "iter0 = read('Ru.traj', index=0)\n", "print('Energy: ', iter0.get_potential_energy())\n", "print('Forces: ', iter0.get_forces())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Often you are only interested in the final energy which can be found like this:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "e_slab = slab.get_potential_energy()\n", "print(e_slab)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Making a Nitrogen molecule\n", "We now make an N2 molecule and optimise it in the same unit cell as we used for the slab." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "d = 1.10\n", "molecule = Atoms('2N', positions=[(0., 0., 0.), (0., 0., d)])\n", "molecule.set_cell(slab.get_cell())\n", "molecule.center()\n", "calc_mol = GPAW(xc='PBE', mode=PW(350))\n", "molecule.set_calculator(calc_mol)\n", "dyn2 = QuasiNewton(molecule, trajectory='N2.traj')\n", "dyn2.run(fmax=0.05)\n", "e_N2 = molecule.get_potential_energy()" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "We can calculate the bond length like this:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "d_N2 = molecule.get_distance(0, 1)\n", "print(d_N2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How does this compare with the experimental value?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Adsorbing the molecule\n", "\n", "Now we adsorb the molecule on top of one of the Ru atoms. \n", "\n", "Here, it would be natural to just add the molecule to the slab, and minimize. However, that takes 45 minutes to an hour to converge, **so we cheat to speed up the calculation.**\n", "\n", "The main slowing-down comes from the relaxation of the topmost metal atom where the N2 molecule binds, this atom moves a quarter of an \u00c5ngstr\u00f6m out. Also, the binding length of the molecule changes when it is adsorbed, so we build a new molecule with a better starting guess." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "h = 1.9 # guess at the binding height\n", "d = 1.2 # guess at the binding distance\n", "slab.positions[4, 2] += 0.2 # pre-relax the binding metal atom.\n", "\n", "molecule = Atoms('2N', positions=[(0, 0, 0), (0, 0, d)])\n", "p = slab.get_positions()\n", "molecule.translate(p + (0, 0, h))\n", "slabN2 = slab + molecule \n", "constraint = FixAtoms(mask=(z < z.min() + 1.0))\n", "slabN2.set_constraint(constraint)\n", "view(slabN2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We optimise the structure. Since we have cheated and have a good guess for the initial configuration we prevent that the optimization algorithm takes too large steps." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "slabN2.set_calculator(calc)\n", "dyn = QuasiNewton(slabN2, trajectory='N2Ru-top.traj', maxstep=0.02)\n", "t = time.time()\n", "dyn.run(fmax=0.05)\n", "print('Wall time: {} min.'.format((time.time() - t) / 60))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The calculation will take a while (10-15 minutes). While it is running please follow the guidelines in the **Exercise** section below." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Once the calculation is finished we can calculate the adsorption energy as:\n", "\n", "Eads = Eslab+N2 - (Eslab + EN2)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "print('Adsorption energy:', slabN2.get_potential_energy() - (e_slab + e_N2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Try to calculate the bond length of N2 adsorbed on the surface. Has it changed? What is the distance between the N2 molecule and the surface?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise\n", "\n", "1) Make a new notebook and set up an adsorption configuration where the N2 molecule is lying down with the center of mass above a three-fold hollow site as shown below. Use an adsorption height of 1.7 \u00c5.\n", "\n", " \n", "\n", "Remember that you can read in the `traj` files you have saved, so you don't need to optimise the surface again. \n", "\n", "View the combined system before you optimize the structure to ensure that you created what you intended." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "slab = read('Ru.traj')\n", "view(slab)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that when viewing the structure, you can find the index of the individual atoms in the ``slab`` object by clicking on them.\n", "\n", "You might also find the [`get_center_of_mass()`](https://wiki.fysik.dtu.dk/ase/ase/atoms.html#ase.Atoms.get_center_of_mass) and [`rotate()`](https://wiki.fysik.dtu.dk/ase/ase/atoms.html#ase.Atoms.rotate) methods useful. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now you should optimize the structure as you did before with the N2 molecule standing. The calculation will probably bee too long to run interactively in a Notebook. Prototype it here, then interrupt the calculation and copy-paste the relevant cells into a script.\n", "\n", "Check the number of irreducible k-points and then submit the job as a batch job running on that number of CPU cores." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3) Make a configuration where two N atoms are adsorbed in hollow sites on the surface as shown below\n", "\n", " \n", "\n", "Note that here the two N atoms sit on next-nearest hollow sites. An alternative would be to have them on nearest neighbor sites. If you feel energetic you could investigate that as well. Also, there are two different kinds of hollow sites, they are not completely equivalent!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Optimise the structure and get the final energy. Is it favourable to dissociate N2 on the surface? What is the N-N distance now? What does that mean for catalysis?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [] } ], "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.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }