Workshop_3-Assessment-DFT-2022.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-4ba158673a862ad8",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "# CH413 Computational Workshop 3: \n",
    "# Predicting the stability of an oxygen-covered Palladium(100) surface"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-7161d605498db50d",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "**ATTENTION** This notebook contains the homework assignment that will be assessed (Part 3).\n",
    "\n",
    "This workshop has 3 parts.\n",
    "* Part 1: Calculating the adsorption energy of oxygen on diamond(110): \n",
    "* Part 2: Calculating stability phase diagrams for different oxygen-covered diamond phases \n",
    "* Part 3: Assessed Assignment (this note beook, which is on the oxygen-covered Pd(100) surfaces).\n",
    "\n",
    "The assessment contains a number of tasks, which are very similar to the workshop. You need to answer them in order to pass this part of the course. You have until **January 12th 2023 12pm (noon)** to complete this assessment, but **do not leave it to the last minute** Some of the DFT calculations that are required are very slow, so this will take some time. \n",
    "\n",
    "To save existing calculations, you can assign variables to the numbers you've already calculated. This way you don't have to recalculate everything if you take breaks.\n",
    "\n",
    "\n",
    "### IMPORTANT - how to use this notebook\n",
    "* This notebook is read-only, so **you do need to make a copy first (make sure to include your last name)**!\n",
    "* Write your answers in the designated markup boxes after the ``# YOUR ANSWER`` line\n",
    "* You need to **fill in the missing code parts as instructed below.**\n",
    "* You can create new code boxes, but discussion comments need to be written in the designated comment cells. This is required for full grading.\n",
    "* The notebook should not raise any errors, otherwise you won't be able to submit it via the system. So, double-check for the source of errors and ask for help if needed.\n",
    "\n",
    "\n",
    "In this assessed work, you will do the phase diagram for investigating the oxygenated Pd(100) surfaces based on the tutorial work covered in the Parts 1-2 of the workshop. Each task is detailed below in individual parts.\n",
    "\n",
    "<b>If you have any issues, contact Bora Karasulu (bora.karasulu@warwick.ac.uk)</b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-2dec7871f68ec1cf",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "## Part 3: Assessed part of the Workshop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-a43b312d71a7d776",
     "locked": true,
     "points": 30,
     "schema_version": 3,
     "solution": false,
     "task": true
    }
   },
   "source": [
    "## TASK 1: Convergence testing\n",
    "\n",
    "\n",
    "As discussed in the lecture, we should check the convergence of this total energy by varying the basis set cutoff and by varying the k-points grid. In this assessed work, you should use projector-augmented plane-waves (viz. PW module from the GPAW package) instead of a finite-difference/numerical basis set (as used in the workshop).\n",
    "\n",
    "1. Recalculate the energy of the clean Pd(100) surface using different settings to study the convergence of the total energy, first with basis set cutoff energy, and then with k-point.\n",
    "    * Vary the basis set cutoff energy between 300-1000 eV, Note: stick with 1x1x1 k-grid for computational efficiency.\n",
    "    * Vary the k-grid between [(1,1,1), (2,2,1), (3,3,1), (4,4,1)], while sticking with a single cutoff energy value (probably one that gives a good precision:time ratio) from the previous step.\n",
    "    \n",
    "    \n",
    "2. Generate two plots of \n",
    "* $E_{tot}$ vs. basis cutoff energy (300, 400, 500, 600, 700, 800) eV\n",
    "and \n",
    "* $E_{tot}$ vs. k point (1,2,3,4)\n",
    "* along with the time required for each setting\n",
    "\n",
    "\n",
    "\n",
    "Below you'll find some prepared functions to help you plot the results. **Please correctly label the axes.**\n",
    "\n",
    "\n",
    "\n",
    "**Attention:** Some of the calculations with higher-accuracy settings may run a few minutes and require a bit of patience\n",
    "\n",
    "You can add as many cells as you think you'll need, but please fill in your discussion in the labeled comments cell\n",
    "\n",
    "**Note:** Use only single point energies (i.e. do not optimise the geometries with each setting you pick) in order to limit the computational efforts and the time needed to finish the tasks. However, please keep in mind that you need to utilise optimised geometries when producing final plots in Task 2 (e.g. phase diagrams)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Importing all the necessary Python packages / ASE modules here\n",
    "from ase import Atoms\n",
    "from ase.visualize import view\n",
    "from ase.io import read\n",
    "from ase.build import surface, add_adsorbate, molecule, fcc100\n",
    "from ase.optimize import BFGS\n",
    "from gpaw import GPAW,PW\n",
    "\n",
    "from imolecule.atoms_to_pybel import atoms_to_pybel as a2p\n",
    "from openbabel import pybel\n",
    "pybel.ipython_3d = True\n",
    "\n",
    "#Needed for time keeping\n",
    "from time import time\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task1.1: Test the convergence of total energy with respect to the plane-wave cutof energy and k-point grid for the clean Pd(100)-fcc surface.** [10pts] \n",
    "\n",
    "You will first you need to create the model for that. You can use the *ase.build* tools (check <a href=https://wiki.fysik.dtu.dk/ase/ase/build/surface.html#ase.build.surface>here</a>). For calculation time considerations, you should use a 3x3x3 surface slab, with a vacuum thickness of 4 Angstroem in +z/-z directions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now create and visualise the Pd(100) surface slab model here. See ase.build webpage for examples\n",
    "from ase.build import fcc100, surface\n",
    "\n",
    "a2p(clean_Pd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculations for the clean Pd(100) surface.\n",
    "\n",
    "# Plane-wave basis-set cutoff energy to investigate...\n",
    "basis = [300,400,500,600,700,800] #If you think necessary, feel free to explore the points in between.\n",
    "#NOTE: use 1x1x1 k-point grid to speed up the calculations, in real calculations you will need to use the one that you find from the k-pt convergence test, which you will do next.\n",
    "\n",
    "# Use a for loop to repeat the following process for each basis cutoff...\n",
    "energies=[] #array for total energy values\n",
    "times=[] #array for timings\n",
    "\n",
    "#To track the elapsed time \n",
    "#init_time=time() #in the beginning of a step\n",
    "#final_time=time() #at the end of a step, i.e. [current time - init time]\n",
    "\n",
    "for base in basis:\n",
    "    #FILL IN\n",
    "    #print and store energy and time at each loop step\n",
    "   \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "####Basis-set Convergence Plotting#####\n",
    "import numpy as np\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib notebook \n",
    "\n",
    "# Plot the total energies against the Basis-set Cutoff Energy.\n",
    "#Relevant keywords/functions:\n",
    "#plt.ylim(low, high)\n",
    "### xlabel and ylabel\n",
    "\n",
    "##TASK: Try to plot the two curves on the same plot (with proper labelling of the axes)\n",
    "#For this you can use the following functions (Google them for more info!!):\n",
    "#fig, ax1 = plt.subplots()\n",
    "#ax2 = ax1.twinx()  \n",
    "\n",
    "#Please don't forget to label the axes properly and add a legend.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task1.2 Now let's do the same for k-points convergence** [10pts]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now test the convergence of total E against the K-point grid. \n",
    "#Use kxkx1, as we are dealing with a surface slab, so we need only one k-point along the z direction.\n",
    "#Note: to speed up the calculations you should use a small cutoff energy in the calculations (e.g. 400eV)\n",
    "\n",
    "ks = [1,2,3,4] # feel free to also explore further points , if you deem necessary\n",
    "\n",
    "# Use a for loop to repeat the following process for each setting...\n",
    "energies=[] #array for total energy values\n",
    "times=[] #array for timings\n",
    "\n",
    "####FILL IN\n",
    "for k in ks:\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "####K point Convergence Plotting#####\n",
    "\n",
    "####FILL IN\n",
    "import numpy as np\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib notebook \n",
    "\n",
    "\n",
    "#plt.ylim(low, high)\n",
    "### xlabel and ylabel\n",
    "\n",
    "##TASK: Try to plot the two curves on the same plot (with proper labelling of the axes)\n",
    "#For this you can use the following functions (Google them for more info!!):\n",
    "#fig, ax1 = plt.subplots()\n",
    "#ax2 = ax1.twinx()  \n",
    "#Please don't forget to label the axes properly and add a legend.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-58c40bec15334afd",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "**Task 1.3: DISCUSSION OF RESULTS for TASK 1** [10pts]\n",
    "\n",
    "Discuss the convergence behavior and discuss in detail if any of the chosen settings would be numerically precise enough to predict the adsorption energy of different oxygen-Pd phases up to a precision of $\\pm$0.10 eV. (Note: Discuss the precision, not the method accuracy.) Also comment on the time required to get the energy result vs. precision achieved, and whether the going for more prcise setting is justified in terms of time required.\n",
    "\n",
    "Enter discussion text below where it says \"YOUR ANSWER\" by double clicking. \n",
    "Use Markdown and Latex syntax, where necessary. Google 'Markdown cheatsheet'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-a5ee57da5e97df1e",
     "locked": false,
     "points": 10,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "source": [
    "#### YOUR ANSWER: PLEASE DOUBLE CLICK HERE AND ENTER TEXT USING LATEX and Markdown\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-108559ed3c2f29bd",
     "locked": true,
     "points": 40,
     "schema_version": 3,
     "solution": false,
     "task": true
    }
   },
   "source": [
    "## Task 2: Phase Diagram of various oxygen-covered Pd(110) surfaces\n",
    "\n",
    "In this part we are interested in the different Pd-O phases containing various oxygen-containing Pd(100) surfaces. To obtain the oxygenated surfaces, you can use the *add_adsorbate* function in the *ase.build* package. You can place the oxygens  1.6A from the surface (which does not much matter in the end as the structure will be optimised, but giving a reasonable value will speed-up the geometry opt convergence).\n",
    "\n",
    "\n",
    "Part 2 tasks:\n",
    "\n",
    "1. Build and visualise the clean and oxygen terminated Pd(100) phases and then optimise each structure (and also the O2 molecule) using the settings that you have identified in Task 1 :\n",
    "   * phase1: oxygen on 'ontop' position (with 1 oxygen atom per unit cell)\n",
    "   * phase2: oxygen on 'hollow' position (with 1 oxygen atom per unit cell)\n",
    "   * phase3: oxygen on 'bridge' position (with 1 oxygen atom per unit cell)\n",
    "   * phase4: two oxygens on 'hollow' and 'ontop' positions (with 2 oxygen atoms per unit cell)\n",
    "   * phase5: two oxygens on 'bridge' and 'hollow' position (with 2 oxygen atoms per unit cell)\n",
    "   \n",
    "   \n",
    "2. Calculate the adsorption energies, surface area, gradient (i.e., number of oxygen atoms per surface area) for each for all clean and oxygenated phases.\n",
    "\n",
    "\n",
    "3. Plot the phase diagram of adsorption free energy ($\\Delta G^{ad}$) vs. change in chemical potential ($\\Delta\\mu_O$). \n",
    "   \n",
    "   *Note*: In all cases, the number of Pd atoms is the same, so the contribution from the change in bulk atoms can be ignored.\n",
    "\n",
    "\n",
    "4. Discussion of the results from Task 2.3: Identify the most stable phases for each region of chemical potential, in a comment cell.\n",
    "\n",
    "\n",
    "5. Translate the phase diagram for a change in chemical potential to a Temperature vs. Pressure phase diagram in the temperature range of 100 K to 1200K and pressure range of 10E-9 atm to 1atm.\n",
    "\n",
    "   Do this using following expression: \n",
    "   $$ \\Delta\\mu(T,p)= a\\cdot T^4+b\\cdot T^3+c\\cdot T^2+d\\cdot T+ e + \\frac{1}{2}k_BT\\ln{\\frac{p}{p^0}}$$\n",
    "   where $T^0=298.15$ K and $p^0=1$ atm and $a=-4.231\\cdot 10^{-12}$ eV/K, $b=6.505\\cdot 10^{-9}$ eV/K, $c=-3.401\\cdot 10^{-6}$ eV/K, $d=-1.259\\cdot 10^{-3}$ eV, $e=-8.883\\cdot 10^{-2}$ eV\n",
    "   \n",
    "   Below is given a function that outputs the chemical potential corresponding to a given pair of T [in K] and p [in atm]\n",
    "\n",
    "\n",
    "6. Discuss your results from Task2.5 results in detail and answer the questions (listed at the end of the notebook).\n",
    "   \n",
    "\n",
    "\n",
    "**NOTE:** If the geometry optimisations are taking too much time (as we are using only one core), you can use the preoptimised structures (optimsied with 600 eV cutoff and 2x2x1 k-point grid) provided for all phases in the 'Part3_DFT/structures/' folder (viz *phaseX-opt.xyz*), also see *ase.io.read* tool. Even if you decide to  use these preoptimised structures, you need to show how to build each model in your notebook. \n",
    "\n",
    "**NOTE2:** If you are still having issues with calculation times with the converged settings you determined in *Task 1*, feel free to use 600eV cutoff and 2x2x1 K-point grid to accelerate the geometry optimisations of the preoptimised geometries, but make sure that you clearly indicate this is your notebook and discuss what adverse effects this might have on the phase diagrams you calculate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib notebook \n",
    "\n",
    "from ase.build import fcc100, add_adsorbate\n",
    "from ase.optimize import BFGS\n",
    "from gpaw import PW,GPAW\n",
    "import os\n",
    "import ase.io\n",
    "\n",
    "#Visualisation related\n",
    "from imolecule.atoms_to_pybel import atoms_to_pybel as a2p\n",
    "from openbabel import pybel\n",
    "pybel.ipython_3d = True\n",
    "from ase.visualize import view\n",
    "\n",
    "#a2p(INSERT STRUCTURE)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-9bc38a7c3749e02d",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def deltamu_Tp(T,p):\n",
    "    \"\"\"\n",
    "    takes T in Kelvin and p in atm and converts to Delta mu in eV\n",
    "    The function is a polynomial fit based on the data provided here:\n",
    "    https://janaf.nist.gov/tables/O-029.html\n",
    "    \n",
    "    Use this function to translate T and P values into delta mu values.\n",
    "    \"\"\"\n",
    "    from math import log\n",
    "    a = -4.231E-12\n",
    "    b =  6.505E-09\n",
    "    c = -3.401E-06\n",
    "    d = -1.259E-03\n",
    "    e = -8.883E-02\n",
    "    kB = 8.6173E-05\n",
    "    \n",
    "    mu_Tp0 = a*(T**4)+b*(T**3)+c*(T**2)+d*T+e\n",
    "    mu_Tp = mu_Tp0 + 0.5*T*kB*log(p) \n",
    "    return mu_Tp\n",
    "\n",
    "#### Example: The value of chemical potential at 600K and standard pressure is \n",
    "print(deltamu_Tp(600,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So let's get started:\n",
    "\n",
    "**Task 2.1 Building models for molecular O2 and all Pd phases and geometry optimisation calculations** [20pts]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#DEFINE a function to do a geometry optimisation using BFGS algorithm\n",
    "from ase.optimize import BFGS #we use a geometry optimisation algorithm from ASE package (among many others available)\n",
    "\n",
    "\n",
    "def geom_opt(atoms):\n",
    "    #FILL IN\n",
    "    return atoms\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now write a function that returns an ASE calculator object, as you will call this function several times while optimising each structure\n",
    "def mycalc():\n",
    "    calc = GPAW(\n",
    "        #FILLIN\n",
    "    )\n",
    "    return calc\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first do the molecular oxygen system, you need to first set up a peridic O2 model (with a proper vacuum padding, e.g. 3A or larger in each direction), and then setup a calculator and do the actual geometry optimisation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up the molecular gas-phase oxygen system.\n",
    "\n",
    "# Re-create the unit cell from the workshop (there is no reason to change this).\n",
    "# Therefore, unit cell will be 6 Angstrom cube containing a single O2 molecule, with bond length 1 Angstrom.\n",
    "# You could also use the ase.build.molecule('O2') to create an O2 model.\n",
    "\n",
    "a = 6.0\n",
    "cell = [[a,0,0],[0,a,0],[0,0,a]]\n",
    "o2 = Atoms('O2', positions = [[0,0,0],[0,0,1.0]])\n",
    "\n",
    "# Set the cell to the O2 molecule.\n",
    "o2.set_cell(cell)\n",
    "# Ensure that there is periodicity in all 3 directions (this is a gas).\n",
    "o2.set_pbc([True,True,True])\n",
    "# Centre the molecule in the cell.\n",
    "o2.center()\n",
    "print(o2)\n",
    "\n",
    "# Visualize O2.\n",
    "#a2p(o2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now optimise the O2 molecule (using the two functions you defined above), and store the final energy\n",
    "\n",
    "o2.set_calculator(mycalc())\n",
    "etot_o2 = o2.get_potential_energy()\n",
    "print('etot_o2 before opt: ', etot_o2) #before the geometry optimisation\n",
    "\n",
    "o2=geom_opt(o2)\n",
    "etot_o2 = o2.get_potential_energy()\n",
    "print('etot_o2 after opt: ', etot_o2) #after the geometry optimisation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now build and optimise the clean_Pd surface using the two funcions you defined above, and store the final energy. \n",
    "from ase.build import fcc100  #you can use this function as you did in Task 1\n",
    "\n",
    "#NOTE: these DFT calculations are slow! So depending on whether you you pipe the output to a text file in your calculation, e.g. GPAW(txt=FILE) option), \n",
    "#and on your DFT settings (basis-set and k-points), it may take a few minutes to get something printed out on the screen as output, as \n",
    "\n",
    "#print energies before and after geometry optimisation, store the final energy.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's visualise the trajectory\n",
    "import ase.io\n",
    "from imolecule.atoms_to_pybel import atoms_to_pybel as a2p\n",
    "from openbabel import pybel\n",
    "pybel.ipython_3d = True\n",
    "from ase.visualize import view\n",
    "\n",
    "\n",
    "steps=ase.io.read('opt.traj',index=':') #read the trajectory from the .traj file, index option is needed for reading all time steps, by default last step is read.\n",
    "view(steps,viewer='ngl')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now optimise the first oxygenated-Pd surface (phase1, ontop), and store the final energy. \n",
    "#As discussed above you can read in the preoptimised structures to speed up optimisations, \n",
    "#but you still need to show how to build the model from scratch using add_adsorbate function from ase.build package (see webpage for details).\n",
    "#print energies before and after geometry optimisation, store the final energy.\n",
    "\n",
    "#Better to copy the clean surface model to another object first, and add the adsorbate(s) to this new phase object, e.g.\n",
    "#phase1=clean_Pd.copy()\n",
    "\n",
    "from ase.build import add_adsorbate\n",
    "#phase1=ase.io.read('Part3_DFT/structures/phase1-opt.xyz') #or one can use the pre-optimised structure provided\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's visualise the trajectory\n",
    "import ase.io\n",
    "steps=ase.io.read('opt.traj',index=':') #read the trajectory from the .traj file, index option is needed for reading all time steps, by default last step is read.\n",
    "view(steps,viewer='ngl')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now optimise the oxygenated-Pd surface (phase2, hollow), and store the final energy. \n",
    "#As discussed above you can read in the preoptimised structures , but you still need to show how to build the model from scratch using add_adsorbate function from ase.build package.\n",
    "#Similar settings as phase 1, now we consider a different adsorption 'position' in add_adsorbate\n",
    "\n",
    "#phase2=clean_Pd.copy()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's visualise the trajectory\n",
    "import ase.io\n",
    "steps=ase.io.read('opt.traj',index=':') #read the trajectory from the .traj file, index option is needed for reading all time steps, by default last step is read.\n",
    "view(steps,viewer='ngl')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now optimise the first oxygenated-Pd surface (phase3, bridge), and store the final energy. \n",
    "#As discussed above you can read in the preoptimised structures , but you still need to show how to build the model from scratch using add_adsorbate function from ase.build package.\n",
    "\n",
    "#phase3=clean_Pd.copy()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's visualise the trajectory\n",
    "import ase.io\n",
    "steps=ase.io.read('opt.traj',index=':') #read the trajectory from the .traj file, index option is needed for reading all time steps, by default last step is read.\n",
    "view(steps,viewer='ngl')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now optimise the oxygenated-Pd surface (phase4, holloow+ontop), and store the final energy. \n",
    "#As discussed above you can read in the preoptimised structures , but you still need to show how to build the model from scratch using add_adsorbate function from ase.build package.\n",
    "#HINT: for adding two oxygen in the unit cell, call add_adsorbate twice with different 'position' settings.\n",
    "\n",
    "#phase4=clean_Pd.copy()\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's visualise the trajectory\n",
    "import ase.io\n",
    "steps=ase.io.read('opt.traj',index=':') #read the trajectory from the .traj file, index option is needed for reading all time steps, by default last step is read.\n",
    "view(steps,viewer='ngl')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now optimise the oxygenated-Pd surface (phase5, holloow+bridge), and store the final energy. \n",
    "#As discussed above you can read in the preoptimised structures , but you still need to show how to build the model from scratch using add_adsorbate function from ase.build package.\n",
    "#HINT: for adding two oxygen in the unit cell, call add_adsorbate twice with different 'position' settings.\n",
    "\n",
    "#phase5=clean_Pd.copy()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's visualise the trajectory\n",
    "import ase.io\n",
    "steps=ase.io.read('opt.traj',index=':') #read the trajectory from the .traj file, index option is needed for reading all time steps, by default last step is read.\n",
    "view(steps,viewer='ngl')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#OPTIONAL: \n",
    "#If you would like to store your optimised structures to use again at a later point (i.e. after reconnecting to the Notebook Server), you can save them using:\n",
    "#clean_Pd.write('clean_opt.xyz',format='extxyz')\n",
    "#phase1.write('phase1-opt.xyz',format='extxyz')\n",
    "#...\n",
    "\n",
    "#and then read them again using e.g.\n",
    "#phase1=ase.io.read('phase1-opt.xyz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 2.2: Adsorption energy, surface area and gradient (Number of O atoms per surface area) calculations for all phases** [10pts]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Task 2.2 #Adsorption energy calculations for all phases\n",
    "\n",
    "# Define a function to calculate the adsorption energies associated with each phase (1,2,4,5).\n",
    "def get_Eads(etot_o2, etot_clean, etot, phase):\n",
    "    #FILL IN\n",
    "    return e_ads\n",
    "\n",
    "#use get_Eads to calculate and print e_ads for each phase.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define a function that returns the gradient (no of oxygen atoms per surface area) for a given structure\n",
    "#and then collect the gradients for all phases\n",
    "import numpy as np\n",
    "def getGrad(atoms):\n",
    "    #FILL IN\n",
    "    return (N_O/area)\n",
    "\n",
    "\n",
    "# Collect the gradients for all phases using the function.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task2.3  Plot the phase diagram (G vs. $\\Delta\\mu_O$)** [10pts]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#Task2.3  ###Phase diagram plot G vs. mu\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib notebook \n",
    "#this gives us an interactive graph with a cross hair\n",
    "\n",
    "#colour the graph with respect to the stable phase in a given chemical potential range\n",
    "#don't forget to add proper labels to axes and the inset, and a legend.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-6d8e733edc1e8caf",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "**Task 2.4 Identify the most stable phases for each region of chemical potential (mu_O2). Note and discuss this in the next comment cell** [10pts]\n",
    "\n",
    "Enter discussion text below where it says \"YOUR ANSWER\" by double clicking. \n",
    "\n",
    "mu range -X.0 - -Y.0 : Phase X is most stable\n",
    "etc...\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-3dbf517970042989",
     "locked": false,
     "points": 10,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "source": [
    "#### YOUR ANSWER: PLEASE DOUBLE CLICK HERE AND ENTER TEXT USING LATEX and Markdown\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 2.5: Now translate the phase diagram into *T* vs. *p* diagram** [10 pts]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Task 2.5 ####Translate phase diagram into T vs. p diagram\n",
    "#HINT: First, generate values of mu for a range of temperatures (100, 200, etc. steps of hundred) and \n",
    "#pressures (10E-9,10E-8, etc. logarithmic)\n",
    "\n",
    "#Generate a two-dimensional numpy array of mu values mu = np.zeros[(N_temps,N_pressures)] and fill it with mu values.\n",
    "#For example:\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib notebook \n",
    "\n",
    "#x, y = np.meshgrid([T1,T2,T3,...], [....])\n",
    "#z = deltamu_Tp(x,y)\n",
    "##Now assign each value of mu with a number identifier for the phase (1,2,3,...)\n",
    "\n",
    "\n",
    "#Check out following examples to learn how to plot a contour plot\n",
    "#   https://matplotlib.org/examples/pylab_examples/contour_demo.html\n",
    "#   https://matplotlib.org/gallery/images_contours_and_fields/contour_corner_mask.html\n",
    "    \n",
    "\n",
    "#Colour and label the different regions of the phase diagram according to the most stable phase.\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-7b363988c25c6c2c",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "**Task 2.6 DISCUSSION OF RESULTS** [10pts]\n",
    "\n",
    "   * What is the most stable phase at ambient pressure and temperature conditions?\n",
    "   * What happens when we heat the surface (assuming no kinetic barriers)?\n",
    "   * What happens when we leave the temperature constant and increase the pressure?\n",
    "   * In Lecture 5, we have seen that also the bulk oxide phase is a part of the Pd-O phase diagram. Considering the heat of formation ($ΔH_f$) given in the Lecture 5 (slide 19), under which conditions would this bulk oxide become relevant?\n",
    "   * Discuss possible reasons why your calculated phase digrams differ from the published ones (shown on the Lecture slide 20, taken from Rogal et al., Phys Rev Lett 98, 046101 (2007)) \n",
    "   * (Optional) if you have decided to use the preoptimised structures and the preset settings (600eV and 2x2x1, rather than the converged values), discuss the possible effects of your choice on the resulting phase diagrams\n",
    "\n",
    "\n",
    "Enter discussion text below where it says \"YOUR ANSWER\" by double clicking. \n",
    "Use Markdown and Latex syntax, where necessary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-4cb78d8045852001",
     "locked": false,
     "points": 10,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "source": [
    "#### YOUR ANSWER: PLEASE DOUBLE CLICK HERE AND ENTER TEXT USING LATEX and Markdown\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Bonus task (not graded):** Normally, you should compute the enthalpy of formation of the bulk oxide phase with the same settings you used for calculating the other phases for consistency. If you are up for an additional task, you can try to obtain the experimentally-observed tetragonal PdO phase from Materials Project Database (Entry ID: 1336, https://materialsproject.org/materials/mp-1336?chemsys=Pd-O), and optimise it with the same settings and add the calculated $ΔH_f$ value into the phase diagrams."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Create Assignment",
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}