07_Adsorption_energy_calculation_in_QE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adsorption energy calculation using DFT\n",
    "\n",
    "## Adsorption Energy\n",
    "- Adsorption energy refers to the change of energy when an adsorbate is attached to the surface of an adsorbent. \n",
    "- Adsorption energy, ${\\Delta}E_{a}$, can be obtained by calculating the difference between the energy of the adsorbed surface and the sum of the energy for each structure composing the adsorbed surface:\n",
    "<br> ${\\Delta}E_{a} =  E_{ads,slab} -  E_{slab} - E_{ads}$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating Adsorption Energy in QE\n",
    "\n",
    "### Consider O$_2$ adsorption on Platinum (Pt) surface with constrained atom positions\n",
    "\n",
    "For this situation, we will first build the Atoms object for Platinum slab and adsorb an oxygen molecule to the slab.\n",
    "\n",
    "#### Building atoms\n",
    "\n",
    "* In the following script, a Pt slab is created. The bottom two layers of Pt (in z-direction) are fixed. The exact position to divide the atoms is determined by visualizing the ‘input.traj’ file with ase gui command. By clicking on the atoms in the pop-out window, the exact position within a cell will be displayed. After applying the constraint to the system, visualizing ‘input.traj’ file with ase gui command, you should be able to see the bottom two layers of Pt atoms have a cross on them, indicating that they are fixed. Another way to check or to manipulate the constraint is by ‘POSCAR’, where it lists the positions and constraints of all atoms in the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.build import bulk , molecule , surface , add_adsorbate\n",
    "from ase.constraints import FixAtoms\n",
    "from ase.io import read\n",
    "\n",
    "#building Platinum slab\n",
    "\n",
    "bulk = bulk ('Pt')\n",
    "miller_indices = (1, 1, 1)\n",
    "layers = 4\n",
    "slab = surface (bulk , indices = miller_indices , layers = layers , vacuum =10)\n",
    "slab *= (2 ,2 ,1)\n",
    "slab.write('slab.traj')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Build an O2 molecule with ase library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.build import molecule, bulk\n",
    "from ase.data.pubchem import pubchem_atoms_search\n",
    "import numpy as np\n",
    "from shutil import copy\n",
    "import os\n",
    "\n",
    "#building oxygen molecule\n",
    "\n",
    "atoms = molecule('O2')\n",
    "max_pos = np.amax(atoms.get_positions(), axis=0)\n",
    "atoms.set_cell(np.array([10, 10, 10]) + max_pos)\n",
    "atoms.center()\n",
    "\n",
    "atoms.write('O2.traj')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Attach the O$_2$ to the Pt slab using ase.build.add_adsorbate()."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.build import bulk , molecule , surface , add_adsorbate\n",
    "from ase.constraints import FixAtoms\n",
    "from ase.io import read\n",
    "\n",
    "slab = read('slab.traj')\n",
    "adsorbate = molecule('O2')\n",
    "add_adsorbate(slab, adsorbate, height = 2.5, position = (1.3 ,1) )\n",
    "c = FixAtoms ( indices = [ atom . index for atom in slab if atom . position [2] <14 ]) #constrained atoms (fixing the bottom 2 layers)\n",
    "slab.set_constraint(c)\n",
    "slab.write('slab_ads.traj')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### DFT calculation \n",
    "- In the following script, we define a Quantum Espresso (QE) calculator object and attach it to an ase.Atoms object so that ase knows how to get the data. After defining the calculator, a DFT calculation can be run. The potential energy can be obtained by the method: get_potential_energy(). The potential energy result can be stored in a text file e.g. 'energy.txt'.\n",
    "- Details on how to submit these calculations can be seen in <a href=\"https://github.com/medford-group/dft-training/blob/devel_jagriti/05_Running_DFT_calculation_on_PACE.ipynb\">Lecture 5</a>.\n",
    "- Clean slab calculation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.io import read\n",
    "from ase.optimize import BFGSLineSearch \n",
    "from espresso import iEspresso as Espresso\n",
    "import numpy as np\n",
    "\n",
    "image = read('slab.traj')\n",
    "\n",
    "# setup calculator\n",
    "parameters = dict(xc='BEEF-vdW',\n",
    "        kpts=(4, 4, 1), #only need 1 kpt in z-direction\n",
    "        pw=400.,\n",
    "        dw=4000.,\n",
    "        spinpol=False,\n",
    "        convergence={'energy':1e-6,\n",
    "                    'mixing':0.05,\n",
    "                    'maxsteps':1000,\n",
    "                    'diag':'david'},\n",
    "        startingwfc='atomic',\n",
    "        smearing='fd', #fermi-dirac electron smearing\n",
    "        sigma=0.1, #smearing width\n",
    "        dipole={'status':True}, #dipole corrections True turns them on\n",
    "        outdir ='esp.log')\n",
    "\n",
    "calc = Espresso(atoms = image, **parameters)\n",
    "\n",
    "image.set_calculator(calc)\n",
    "\n",
    "relax = BFGSLineSearch(image,\n",
    "                       trajectory='opt.traj',\n",
    "                       logfile='opt.log',\n",
    "                       restart='opt.pckl')\n",
    "\n",
    "relax.run(fmax=0.05)\n",
    "\n",
    "eng = image.get_potential_energy()\n",
    "\n",
    "image.write('converged.traj')\n",
    "\n",
    "with open('slab_energy.txt', 'w') as f:\n",
    "    f.write(str(eng))\n",
    "image.calc.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- O2 on the Pt slab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.io import read\n",
    "from ase.optimize import BFGSLineSearch \n",
    "from espresso import iEspresso as Espresso\n",
    "import numpy as np\n",
    "\n",
    "image = read('slab_ads.traj')\n",
    "\n",
    "# setup calculator\n",
    "parameters = dict(xc='BEEF-vdW',\n",
    "        kpts=(4, 4, 1), #only need 1 kpt in z-direction\n",
    "        pw=400.,\n",
    "        dw=4000.,\n",
    "        spinpol=False,\n",
    "        convergence={'energy':1e-6,\n",
    "                    'mixing':0.05,\n",
    "                    'maxsteps':1000,\n",
    "                    'diag':'david'},\n",
    "        startingwfc='atomic',\n",
    "        smearing='fd', #fermi-dirac electron smearing\n",
    "        sigma=0.1, #smearing width\n",
    "        dipole={'status':True}, #dipole corrections True turns them on\n",
    "        outdir ='esp.log')\n",
    "\n",
    "calc = Espresso(atoms = image, **parameters)\n",
    "\n",
    "image.set_calculator(calc)\n",
    "\n",
    "relax = BFGSLineSearch(image,\n",
    "                       trajectory='opt.traj',\n",
    "                       logfile='opt.log',\n",
    "                       restart='opt.pckl')\n",
    "\n",
    "relax.run(fmax=0.05)\n",
    "\n",
    "eng = image.get_potential_energy()\n",
    "\n",
    "image.write('converged.traj')\n",
    "\n",
    "with open('slab_ads_energy.txt', 'w') as f:\n",
    "    f.write(str(eng))\n",
    "image.calc.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- O2 gas reference energy calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.io import read\n",
    "from ase.optimize import BFGSLineSearch \n",
    "from espresso import iEspresso as Espresso\n",
    "import numpy as np\n",
    "\n",
    "image = read('O2.traj')\n",
    "\n",
    "# setup calculator\n",
    "parameters = dict(xc='BEEF-vdW',\n",
    "        kpts=(4, 4, 1), #only need 1 kpt in z-direction\n",
    "        pw=400.,\n",
    "        dw=4000.,\n",
    "        spinpol=False,\n",
    "        convergence={'energy':1e-6,\n",
    "                    'mixing':0.05,\n",
    "                    'maxsteps':1000,\n",
    "                    'diag':'david'},\n",
    "        startingwfc='atomic',\n",
    "        smearing='fd', #fermi-dirac electron smearing\n",
    "        sigma=0.1, #smearing width\n",
    "        dipole={'status':True}, #dipole corrections True turns them on\n",
    "        outdir ='esp.log')\n",
    "\n",
    "calc = Espresso(atoms = image, **parameters)\n",
    "\n",
    "image.set_calculator(calc)\n",
    "\n",
    "relax = BFGSLineSearch(image,\n",
    "                       trajectory='opt.traj',\n",
    "                       logfile='opt.log',\n",
    "                       restart='opt.pckl')\n",
    "\n",
    "relax.run(fmax=0.05)\n",
    "\n",
    "eng = image.get_potential_energy()\n",
    "\n",
    "image.write('converged.traj')\n",
    "\n",
    "with open('O2_energy.txt', 'w') as f:\n",
    "    f.write(str(eng))\n",
    "image.calc.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Analysis of adsorption energies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from ase.io import read\n",
    "\n",
    "def readFile(path):\n",
    "    f = open(path, 'r')\n",
    "    content = f.readlines()\n",
    "    return content\n",
    "\n",
    "def read_energy(path):\n",
    "    energy = float(readFile(path)[0])\n",
    "    return energy\n",
    "\n",
    "E_slab = read_energy('slab_energy.txt')\n",
    "E_slab_ads = read_energy('slab_ads_energy.txt')\n",
    "E_O2 = read_energy('O2.txt')\n",
    "\n",
    "# E_adsorption is the adsorption energy\n",
    "E_adsorption = E_slab_ads - E_slab - E_O2\n",
    "\n",
    "print('The adsorption energy of O2 on Pt slab is {0:1.3f} eV'.format(E_adsorption))\n"
   ]
  }
 ],
 "metadata": {
  "authors": [
   {
    "name": "Kaylee Tian"
   }
  ],
  "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}