{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Phase Diagram\n",
    "In this example, we will generate a 1D and 2D phase diagram for a CO/Pt(111) system.\n",
    "\n",
    "## Topics Covered\n",
    "- Create ``StatMech`` objects \n",
    "- Initialize ``Reaction`` objects to describe the formation reaction of CO/Pt(111) species\n",
    "- Generate a 1D phase diagram by varying T\n",
    "- Generate a 2D phase diagram by varying T and P\n",
    "- Save the ``PhaseDiagram`` object as a ``JSON`` file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create Species for Phase Diagram\n",
    "We will be considering six CO/Pt(111) configurations. The configurations have CO adsorbed in different sites and different coverages.\n",
    "\n",
    "<img src=\"images/configurations.png\" width=600>\n",
    "\n",
    "First, we initialize the species as a dictionary to enable easy `Reaction` initialization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from ase.build import molecule\n",
    "\n",
    "from pmutt.statmech import StatMech, presets\n",
    "\n",
    "species = {\n",
    "    'CO': StatMech(name='CO', atoms=molecule('CO'),\n",
    "                   potentialenergy=-14.8021,\n",
    "                   vib_wavenumbers=[2121.2], symmetrynumber=1, \n",
    "                   **presets['idealgas']),\n",
    "    'Pt': StatMech(name='Pt', potentialenergy=-383.161235,\n",
    "                   **presets['electronic']),\n",
    "    'CO(S) 1/16ML fcc': StatMech(\n",
    "            name='CO(S) 1/16ML fcc', potentialenergy=-399.48282843,\n",
    "            vib_wavenumbers=[1731.942697, 349.970617, 322.15111,\n",
    "                             319.114152, 161.45669],\n",
    "            **presets['harmonic']),\n",
    "    'CO(S) 1/16ML br': StatMech(\n",
    "            name='CO(S) 1/16ML br', potentialenergy=-399.464095,\n",
    "            vib_wavenumbers=[1831.626557, 394.436054, 388.098645,\n",
    "                             373.063005, 203.887416, 52.987012],\n",
    "            **presets['harmonic']),\n",
    "    'CO(S) 1/16ML top': StatMech(\n",
    "            name='CO(S) 1/16ML top', potentialenergy=-399.39545350,\n",
    "            vib_wavenumbers=[2045.797559, 489.514815, 396.498284,\n",
    "                             393.395406, 56.058884, 52.157548],\n",
    "            **presets['harmonic']),\n",
    "    'CO(S) 1/8ML': StatMech(\n",
    "            name='CO(S) 1/8ML', potentialenergy=-415.67626828,\n",
    "            vib_wavenumbers=[2047.452988, 1730.209946, 482.24755,\n",
    "                             394.675312, 392.79586, 354.078848,\n",
    "                             323.143303, 320.375056, 162.356233,\n",
    "                             158.239412, 60.269377, 51.362263],\n",
    "            **presets['harmonic']),\n",
    "    'CO(S) 3/16ML': StatMech(\n",
    "            name='CO(S) 3/16ML', potentialenergy=-431.867618,\n",
    "            vib_wavenumbers=[2049.767728, 1746.427506, 1733.474666,\n",
    "                             478.755939, 391.899407, 389.661616,\n",
    "                             354.568306, 352.532192, 325.154407,\n",
    "                             322.578758, 319.593333, 315.883097,\n",
    "                             163.2316, 162.672434, 158.815096,\n",
    "                             157.87804, 59.576319, 50.284495],\n",
    "            **presets['harmonic']),\n",
    "    'CO(S) 1/2ML': StatMech(\n",
    "            name='CO(S) 1/2ML', potentialenergy=-512.817507,\n",
    "            vib_wavenumbers=[2072.099888, 2053.332551, 2052.632444,\n",
    "                             2052.501762, 1835.620624, 1824.088854,\n",
    "                             1823.712945, 1823.531493, 481.148383,\n",
    "                             480.426246, 480.187182, 479.70589, \n",
    "                             414.42128, 411.357815, 411.091615,\n",
    "                             406.851876, 404.128284, 403.391877, \n",
    "                             402.879585, 401.452804, 401.134231,\n",
    "                             397.539281, 394.569066, 394.101234,\n",
    "                             393.933956, 390.740547, 390.173637,\n",
    "                             389.805187, 388.420025, 387.427067,\n",
    "                             383.620218, 383.348263, 201.654999,\n",
    "                             200.123762, 196.698042, 195.736534, \n",
    "                             75.269065, 72.94012, 70.402739,\n",
    "                             68.651958, 65.289743, 64.556735,\n",
    "                             63.904694, 60.051442, 58.698334,\n",
    "                             55.589005, 52.608038, 43.883525],\n",
    "            **presets['harmonic']),\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create Reactions for Phase Diagram\n",
    "The reactions will be initialized and put in a list. Notice that the stoichiometric coefficient of CO changes for higher coverages. If you are unfamiliar with initializing reactions, see the Reactions example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pmutt.reaction import Reaction\n",
    "\n",
    "reactions=[\n",
    "    Reaction.from_string('Pt = Pt', species), # Clean surface\n",
    "    Reaction.from_string('Pt + CO = CO(S) 1/16ML fcc', species),\n",
    "    Reaction.from_string('Pt + CO = CO(S) 1/16ML br', species),\n",
    "    Reaction.from_string('Pt + CO = CO(S) 1/16ML top', species),\n",
    "    Reaction.from_string('Pt + 2CO = CO(S) 1/8ML', species),\n",
    "    Reaction.from_string('Pt + 3CO = CO(S) 3/16ML', species),\n",
    "    Reaction.from_string('Pt + 8CO = CO(S) 1/2ML', species)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create PhaseDiagram Object\n",
    "Now we have everything we need to create the ``PhaseDiagram`` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pmutt.reaction.phasediagram import PhaseDiagram\n",
    "\n",
    "phase_diagram = PhaseDiagram(reactions=reactions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a 1D Phase Diagram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 75348 characters]\n",
      "text/plain": [
       "<Figure size 900x600 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "T = np.linspace(300, 1000, 200) # K\n",
    "fig1, ax1 = phase_diagram.plot_1D(x_name='T', x_values=T, P=1., G_units='kJ/mol')\n",
    "\n",
    "'''Plotting adjustments'''\n",
    "# Set colors to lines\n",
    "colors = ('#000080', '#0029FF', '#00D5FF', '#7AFF7D', \n",
    "          '#FFE600', '#FF4A00', '#800000')\n",
    "for color, line in zip(colors, ax1.get_lines()):\n",
    "    line.set_color(color)\n",
    "\n",
    "# Set labels to lines\n",
    "labels = ('0 ML', '1/16 ML (fcc)', '1/16 ML (bridge)', \n",
    "          '1/16 ML (top)', '1/8 ML', '3/16 ML', '1/2 ML')\n",
    "handles, _ = ax1.get_legend_handles_labels()\n",
    "ax1.get_legend().remove()\n",
    "ax1.legend(handles[::-1], labels[::-1], loc=0,\n",
    "           title='CO/Pt Configuration')\n",
    "ax1.set_xlabel('Temperature (K)')\n",
    "\n",
    "fig1.set_dpi(150.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a 2D Phase Diagram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 35480 characters]\n",
      "text/plain": [
       "<Figure size 900x600 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Generate Pressure range\n",
    "T = np.linspace(300, 1000) # K\n",
    "P = np.logspace(-3, 0) # bar\n",
    "\n",
    "fig2, ax2, c2, cbar2 = phase_diagram.plot_2D(x1_name='T', x1_values=T, \n",
    "                                             x2_name='P', x2_values=P)\n",
    "\n",
    "'''Plotting adjustments'''\n",
    "# Change y axis to use log scale\n",
    "ax2.set_yscale('log')\n",
    "# Add axis labels\n",
    "ax2.set_xlabel('Temperature (K)')\n",
    "ax2.set_ylabel('CO Pressure (bar)')\n",
    "# Change color scheme\n",
    "plt.set_cmap('jet')\n",
    "# Add labels\n",
    "cbar2.ax.set_yticklabels(labels)\n",
    "\n",
    "fig2.set_dpi(150.)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
