notebooks/2017-12-15-Plotting a Pourbaix Diagram.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting a pourbaix diagram using the legacy Materials Project API"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**NOTE** this notebook uses the LEGACY Materials Project API, which is frozen at database release `v2021.03.13`. For an updated experience and access to the latest data, we recommend the new API. See [Materials Project API reference](https://next-gen.materialsproject.org/api) for more information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*If you use this infrastructure, please consider citing the following work:*\n",
    "\n",
    "[K. A. Persson, B. Waldwick, P. Lazic and G. Ceder, Phys. Rev. B, 2012, 85, 235438.](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.85.235438)\n",
    "\n",
    "[A. K. Singh, L. Zhou, A. Shinde, S. K. Suram, J. H. Montoya, D. Winston, J. M. Gregoire and K. A. Persson, Chem. Mater., 2017, DOI: acs.chemmater.7b03980.](http://pubs.acs.org/doi/10.1021/acs.chemmater.7b03980)\n",
    "\n",
    "**Notebook Author:** Joseph Montoya, Ryan Kingsbury\n",
    "\n",
    "**Last Update** 2023-03-10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----------\n",
    "\n",
    "The Materials Project REST interface includes functionality to construct pourbaix diagrams from computed entries.  Note that the Pourbaix infrastructure is still undergoing revision, but now includes a simplified interface that enables MP and pymatgen users to fetch entries that have been processed according to the [Materials Project Aqueous Compatibility scheme](https://docs.materialsproject.org/methodology/materials-methodology/aqueous-stability-pourbaix).  Thus, users can reproduce web Pourbaix diagrams in two or three steps in pymatgen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncomment the subsequent lines in this cell to install dependencies for Google Colab.\n",
    "# !pip install pymatgen==2022.7.19\n",
    "from __future__ import annotations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**NOTE** - make sure to paste your **legacy** API key below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "MP_API_KEY = \"<your-legacy-API-key-here>\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import necessary tools from pymatgen\n",
    "from pymatgen.analysis.pourbaix_diagram import PourbaixDiagram, PourbaixPlotter\n",
    "from pymatgen.ext.matproj import MPRester\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "# Initialize the MP Rester\n",
    "mpr = MPRester()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To retrieve entries necessary to construct a Pourbaix Diagram use `MPRester.get_pourbaix_entries(LIST_OF_ELEMENTS)` with a list of elements comprising your chemical system.  It is not necessary to include 'O' and 'H' in your list, as they are added automatically.  This function also contains all of necessary preprocessing to ensure the PourbaixEntries are compatible with the pymatgen `PourbaixDiagram` constructor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/shyue/miniconda3/envs/mavrl/lib/python3.11/site-packages/mp_api/client/mprester.py:368: UserWarning: mpcontribs-client not installed. Install the package to query MPContribs data, or construct pourbaix diagrams: 'pip install mpcontribs-client'\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "ename": "AttributeError",
     "evalue": "'NoneType' object has no attribute 'query_contributions'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[4], line 2\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;66;03m# Get all pourbaix entries corresponding to the Cu-O-H chemical system.\u001b[39;00m\n\u001b[0;32m----> 2\u001b[0m entries \u001b[38;5;241m=\u001b[39m \u001b[43mmpr\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget_pourbaix_entries\u001b[49m\u001b[43m(\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mCu\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/miniconda3/envs/mavrl/lib/python3.11/site-packages/mp_api/client/mprester.py:851\u001b[0m, in \u001b[0;36mMPRester.get_pourbaix_entries\u001b[0;34m(self, chemsys, solid_compat, use_gibbs)\u001b[0m\n\u001b[1;32m    845\u001b[0m chemsys \u001b[38;5;241m=\u001b[39m \u001b[38;5;28msorted\u001b[39m(e\u001b[38;5;241m.\u001b[39mcapitalize() \u001b[38;5;28;01mfor\u001b[39;00m e \u001b[38;5;129;01min\u001b[39;00m chemsys)\n\u001b[1;32m    847\u001b[0m \u001b[38;5;66;03m# Get ion entries first, because certain ions have reference\u001b[39;00m\n\u001b[1;32m    848\u001b[0m \u001b[38;5;66;03m# solids that aren't necessarily in the chemsys (Na2SO4)\u001b[39;00m\n\u001b[1;32m    849\u001b[0m \n\u001b[1;32m    850\u001b[0m \u001b[38;5;66;03m# download the ion reference data from MPContribs\u001b[39;00m\n\u001b[0;32m--> 851\u001b[0m ion_data \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget_ion_reference_data_for_chemsys\u001b[49m\u001b[43m(\u001b[49m\u001b[43mchemsys\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    853\u001b[0m \u001b[38;5;66;03m# build the PhaseDiagram for get_ion_entries\u001b[39;00m\n\u001b[1;32m    854\u001b[0m ion_ref_comps \u001b[38;5;241m=\u001b[39m [\n\u001b[1;32m    855\u001b[0m     Ion\u001b[38;5;241m.\u001b[39mfrom_formula(d[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdata\u001b[39m\u001b[38;5;124m\"\u001b[39m][\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mRefSolid\u001b[39m\u001b[38;5;124m\"\u001b[39m])\u001b[38;5;241m.\u001b[39mcomposition \u001b[38;5;28;01mfor\u001b[39;00m d \u001b[38;5;129;01min\u001b[39;00m ion_data\n\u001b[1;32m    856\u001b[0m ]\n",
      "File \u001b[0;32m~/miniconda3/envs/mavrl/lib/python3.11/site-packages/mp_api/client/mprester.py:989\u001b[0m, in \u001b[0;36mMPRester.get_ion_reference_data_for_chemsys\u001b[0;34m(self, chemsys)\u001b[0m\n\u001b[1;32m    956\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mget_ion_reference_data_for_chemsys\u001b[39m(\u001b[38;5;28mself\u001b[39m, chemsys: \u001b[38;5;28mstr\u001b[39m \u001b[38;5;241m|\u001b[39m \u001b[38;5;28mlist\u001b[39m) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m \u001b[38;5;28mlist\u001b[39m[\u001b[38;5;28mdict\u001b[39m]:\n\u001b[1;32m    957\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"Download aqueous ion reference data used in the construction of Pourbaix diagrams.\u001b[39;00m\n\u001b[1;32m    958\u001b[0m \n\u001b[1;32m    959\u001b[0m \u001b[38;5;124;03m    Use this method to examine the ion reference data and to add additional\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    987\u001b[0m \u001b[38;5;124;03m            compounds and aqueous species, Wiley, New York (1978)'}}\u001b[39;00m\n\u001b[1;32m    988\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 989\u001b[0m     ion_data \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget_ion_reference_data\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    991\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(chemsys, \u001b[38;5;28mstr\u001b[39m):\n\u001b[1;32m    992\u001b[0m         chemsys \u001b[38;5;241m=\u001b[39m chemsys\u001b[38;5;241m.\u001b[39msplit(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m-\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
      "File \u001b[0;32m~/miniconda3/envs/mavrl/lib/python3.11/site-packages/mp_api/client/mprester.py:948\u001b[0m, in \u001b[0;36mMPRester.get_ion_reference_data\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    919\u001b[0m \u001b[38;5;129m@lru_cache\u001b[39m\n\u001b[1;32m    920\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mget_ion_reference_data\u001b[39m(\u001b[38;5;28mself\u001b[39m) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m \u001b[38;5;28mlist\u001b[39m[\u001b[38;5;28mdict\u001b[39m]:\n\u001b[1;32m    921\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"Download aqueous ion reference data used in the construction of Pourbaix diagrams.\u001b[39;00m\n\u001b[1;32m    922\u001b[0m \n\u001b[1;32m    923\u001b[0m \u001b[38;5;124;03m    Use this method to examine the ion reference data and to add additional\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    946\u001b[0m \u001b[38;5;124;03m            compounds and aqueous species, Wiley, New York (1978)'}}\u001b[39;00m\n\u001b[1;32m    947\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 948\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcontribs\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mquery_contributions\u001b[49m(  \u001b[38;5;66;03m# type: ignore\u001b[39;00m\n\u001b[1;32m    949\u001b[0m         query\u001b[38;5;241m=\u001b[39m{\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mproject\u001b[39m\u001b[38;5;124m\"\u001b[39m: \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mion_ref_data\u001b[39m\u001b[38;5;124m\"\u001b[39m},\n\u001b[1;32m    950\u001b[0m         fields\u001b[38;5;241m=\u001b[39m[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124midentifier\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mformula\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdata\u001b[39m\u001b[38;5;124m\"\u001b[39m],\n\u001b[1;32m    951\u001b[0m         paginate\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m,\n\u001b[1;32m    952\u001b[0m     )\u001b[38;5;241m.\u001b[39mget(\n\u001b[1;32m    953\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdata\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    954\u001b[0m     )\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'NoneType' object has no attribute 'query_contributions'"
     ]
    }
   ],
   "source": [
    "# Get all pourbaix entries corresponding to the Cu-O-H chemical system.\n",
    "entries = mpr.get_pourbaix_entries([\"Cu\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pourbaix diagrams can be constructed using `PourbaixDiagram(RETRIEVED_ENTRIES)` as below.  Note that a `comp_dict` keyword argument may also be supplied to the `PourbaixDiagram` constructor if a fixed composition for a multi-element pourbaix diagram is desired."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct the PourbaixDiagram object\n",
    "pbx = PourbaixDiagram(entries)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `PourbaixDiagram` includes a number of useful functions for determining stable species and stability of entries relative to a given pourbaix facet (i.e. as a function of pH and V)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get an entry stability as a function of pH and V\n",
    "entry = [e for e in entries if e.entry_id == \"mp-1692\"][0]\n",
    "print(\n",
    "    \"CuO's potential energy per atom relative to the most\",\n",
    "    f\"stable decomposition product is {pbx.get_decomposition_energy(entry, pH=7, V=-0.2):0.2f} eV/atom\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This suggests that CuO, for example, has a large driving force for decomposition at neutral pH and mildly reducing conditions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see this in more detail, we can plot the pourbaix diagram.  The `PourbaixPlotter` object is also initialized using an instance of the `PourbaixDiagram` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotter = PourbaixPlotter(pbx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotter.get_pourbaix_plot().show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can customize the appearance of the plot using keyword arguments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotter.get_pourbaix_plot(limits=[[0,14], [-2,2]], # pH and E limits\n",
    "                          title=\"Customized Pourbaix Plot\",\n",
    "                            label_domains=True, # set True to label phases\n",
    "                            label_fontsize=10, # font size for the phase labels\n",
    "                            show_neutral_axes=False, # show or hide axes at pH=0, E=0\n",
    "                            show_water_lines=False, # show or hide lines for water stability\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The PourbaixPlotter object can also plot the relative stability of an entry across the pourbaix diagram.  To do this, use the `PourbaixPlotter.plot_entry_stability` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt = plotter.plot_entry_stability(entry)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting k-nary systems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pymatgen also supports binary/ternary pourbaix diagrams with fixed compositions of the non-H or O elements.  This is achieved by finding all possible combinations of entries that fulfill the composition constraint and treating them as individual entries in pourbaix space.  Note that you can supply a composition dictionary to further tune the pourbaix diagram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get all pourbaix entries corresponding to the Bi-V-O-H chemical system.\n",
    "entries = mpr.get_pourbaix_entries([\"Bi\", \"V\"])\n",
    "# Construct the PourbaixDiagram object\n",
    "pbx = PourbaixDiagram(\n",
    "    entries,\n",
    "    comp_dict={\"Bi\": 0.5, \"V\": 0.5},\n",
    "    conc_dict={\"Bi\": 1e-8, \"V\": 1e-8},\n",
    "    filter_solids=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the `filter_solids` keyword argument in the `PourbaixDiagram` instantiation above tells the constructor whether to filter solids by phase stability on the compositional phase diagram for solids. The filtering process significantly reduces the time it takes to generate all of the combined entries for a binary or ternary pourbaix diagram though, so it may be prudent to use in those cases. Pourbaix Diagrams generated with and without this argument may look significantly different in the OER and HER regions, since highly oxidized materials (e. g. Bi$_2$O$_5$) or highly reduced materials (e. g. most hydrides) may not be stable on the compositional phase diagram. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct the pourbaix analyzer\n",
    "plotter = PourbaixPlotter(pbx)\n",
    "plt = plotter.get_pourbaix_plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Getting the heatmaps for a solid entry in these cases is a bit more involved, because many of the regions of the pourbaix diagram include multiphase entries.  Here's an example for this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bivo4_entry = [entry for entry in entries if entry.entry_id == \"mp-545850\"][0]\n",
    "plt = plotter.plot_entry_stability(bivo4_entry)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If a compound is provided that doesn't meet the constraints of the composition (in this case, equal parts Bi and V), you will get a `ValueError`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bio2_entry = [entry for entry in entries if entry.entry_id == \"mp-557993\"][0]\n",
    "plt = plotter.plot_entry_stability(bio2_entry)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}