notebooks/2019-03-11-Interface Reactions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interface Reactions\n",
    "\n",
    "This notebook shows how to\n",
    "* Obtain information about interface reactions between two solid substances in contact, and\n",
    "* Plot reaction energy as a function of mixing ratio.\n",
    "\n",
    "We use the [Materials Project API](https://github.com/materialsproject/mapidoc) to obtain energies of compounds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymatgen.analysis.interface_reactions import InterfacialReactivity\n",
    "from pymatgen.analysis.phase_diagram import GrandPotentialPhaseDiagram, PhaseDiagram\n",
    "from pymatgen.core import Composition, Element\n",
    "from pymatgen.ext.matproj import MPRester\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "# Initialize the REST API interface. You may need to put your own API key in as an arg.\n",
    "mpr = MPRester()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, set the values of the two reactants. Optionally, simulate the case where the reaction system is in contact with an elemental reservoir.\n",
    "\n",
    "Because the methodology here is to generate a pseudo-binary phase stability diagram of two reactants as a function of mixing ratio, the addition of an elemental reservoir implies construction of a so-called _grand potential phase diagram_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Chemical formulae for two solid reactants.\n",
    "reactant1 = \"LiCoO2\"\n",
    "reactant2 = \"Li3PS4\"\n",
    "\n",
    "# Is the system open to an elemental reservoir?\n",
    "grand = True\n",
    "\n",
    "if grand:\n",
    "    # Element in the elemental reservoir.\n",
    "    open_el = \"Co\"\n",
    "    # Relative chemical potential vs. pure substance. Must be non-positive.\n",
    "    relative_mu = -1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, compile the critical reaction information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "efa040645ef248228b61d35fd1d6cbec",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Retrieving ThermoDoc documents:   0%|          | 0/950 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Get the compositions of the reactants\n",
    "comp1 = Composition(reactant1)\n",
    "comp2 = Composition(reactant2)\n",
    "\n",
    "# Gather all elements involved in the chemical system.\n",
    "elements = [e.symbol for e in comp1.elements + comp2.elements]\n",
    "if grand:\n",
    "    elements.append(open_el)\n",
    "elements = list(set(elements))  # Remove duplicates\n",
    "\n",
    "# Get all entries in the chemical system\n",
    "entries = mpr.get_entries_in_chemsys(elements)\n",
    "\n",
    "# Build a phase diagram using these entries.\n",
    "pd = PhaseDiagram(entries)\n",
    "\n",
    "# For an open system, include the grand potential phase diagram.\n",
    "if grand:\n",
    "    # Get the chemical potential of the pure subtance.\n",
    "    mu = pd.get_transition_chempots(Element(open_el))[0]\n",
    "    # Set the chemical potential in the elemental reservoir.\n",
    "    chempots = {open_el: relative_mu + mu}\n",
    "    # Build the grand potential phase diagram\n",
    "    gpd = GrandPotentialPhaseDiagram(entries, chempots)\n",
    "    # Create InterfacialReactivity object.\n",
    "    interface = InterfacialReactivity(\n",
    "        comp1,\n",
    "        comp2,\n",
    "        gpd,\n",
    "        norm=True,\n",
    "        include_no_mixing_energy=True,\n",
    "        pd_non_grand=pd,\n",
    "        use_hull_energy=False,\n",
    "    )\n",
    "else:\n",
    "    interface = InterfacialReactivity(\n",
    "        comp1,\n",
    "        comp2,\n",
    "        pd,\n",
    "        norm=True,\n",
    "        include_no_mixing_energy=False,\n",
    "        pd_non_grand=None,\n",
    "        use_hull_energy=False,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From here, you can plot reaction energy versus mixing ratio as below. Note that the mixing ratio is between the normalized compositions of the reactants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "Li0.25 Co0.25 O0.5 has elements not in the phase diagram S, O, Li, P",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[5], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m plt \u001b[38;5;241m=\u001b[39m \u001b[43minterface\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mplot\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/miniconda3/envs/mavrl/lib/python3.11/site-packages/pymatgen/analysis/interface_reactions.py:203\u001b[0m, in \u001b[0;36mInterfacialReactivity.plot\u001b[0;34m(self, backend)\u001b[0m\n\u001b[1;32m    191\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m    192\u001b[0m \u001b[38;5;124;03mPlots reaction energy as a function of mixing ratio x in self.c1 - self.c2\u001b[39;00m\n\u001b[1;32m    193\u001b[0m \u001b[38;5;124;03mtie line.\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    200\u001b[0m \u001b[38;5;124;03m    Plot of reaction energies as a function of mixing ratio\u001b[39;00m\n\u001b[1;32m    201\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m    202\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m backend\u001b[38;5;241m.\u001b[39mlower() \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mplotly\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[0;32m--> 203\u001b[0m     fig \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_get_plotly_figure\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    204\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m backend\u001b[38;5;241m.\u001b[39mlower() \u001b[38;5;129;01min\u001b[39;00m [\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mmatplotlib\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mmpl\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mplt\u001b[39m\u001b[38;5;124m\"\u001b[39m]:\n\u001b[1;32m    205\u001b[0m     fig \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_get_matplotlib_figure()\n",
      "File \u001b[0;32m~/miniconda3/envs/mavrl/lib/python3.11/site-packages/pymatgen/analysis/interface_reactions.py:335\u001b[0m, in \u001b[0;36mInterfacialReactivity._get_plotly_figure\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    333\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_get_plotly_figure\u001b[39m(\u001b[38;5;28mself\u001b[39m) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m Figure:\n\u001b[1;32m    334\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"Returns a Plotly figure of reaction kinks diagram.\"\"\"\u001b[39;00m\n\u001b[0;32m--> 335\u001b[0m     kinks \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mmap\u001b[39m(\u001b[38;5;28mlist\u001b[39m, \u001b[38;5;28mzip\u001b[39m(\u001b[38;5;241m*\u001b[39m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget_kinks\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m))\n\u001b[1;32m    336\u001b[0m     _, x, energy, reactions, _ \u001b[38;5;241m=\u001b[39m kinks\n\u001b[1;32m    338\u001b[0m     lines \u001b[38;5;241m=\u001b[39m Scatter(\n\u001b[1;32m    339\u001b[0m         x\u001b[38;5;241m=\u001b[39mx,\n\u001b[1;32m    340\u001b[0m         y\u001b[38;5;241m=\u001b[39menergy,\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    344\u001b[0m         hoverinfo\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mnone\u001b[39m\u001b[38;5;124m\"\u001b[39m,\n\u001b[1;32m    345\u001b[0m     )\n",
      "File \u001b[0;32m~/miniconda3/envs/mavrl/lib/python3.11/site-packages/pymatgen/analysis/interface_reactions.py:143\u001b[0m, in \u001b[0;36mInterfacialReactivity.get_kinks\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    133\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mget_kinks\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;28mtuple\u001b[39m[\u001b[38;5;28mint\u001b[39m, \u001b[38;5;28mfloat\u001b[39m, \u001b[38;5;28mfloat\u001b[39m, Reaction, \u001b[38;5;28mfloat\u001b[39m]]:\n\u001b[1;32m    134\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m    135\u001b[0m \u001b[38;5;124;03m    Finds all the kinks in mixing ratio where reaction products changes\u001b[39;00m\n\u001b[1;32m    136\u001b[0m \u001b[38;5;124;03m    along the tie-line of composition self.c1 and composition self.c2.\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    141\u001b[0m \u001b[38;5;124;03m        energy per mol of formula in kJ/mol).\u001b[39;00m\n\u001b[1;32m    142\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 143\u001b[0m     c1_coord \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpd\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpd_coords\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcomp1\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    144\u001b[0m     c2_coord \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mpd\u001b[38;5;241m.\u001b[39mpd_coords(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcomp2)\n\u001b[1;32m    146\u001b[0m     n1 \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcomp1\u001b[38;5;241m.\u001b[39mnum_atoms\n",
      "File \u001b[0;32m~/miniconda3/envs/mavrl/lib/python3.11/site-packages/pymatgen/analysis/phase_diagram.py:512\u001b[0m, in \u001b[0;36mPhaseDiagram.pd_coords\u001b[0;34m(self, comp)\u001b[0m\n\u001b[1;32m    500\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m    501\u001b[0m \u001b[38;5;124;03mThe phase diagram is generated in a reduced dimensional space\u001b[39;00m\n\u001b[1;32m    502\u001b[0m \u001b[38;5;124;03m(n_elements - 1). This function returns the coordinates in that space.\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    509\u001b[0m \u001b[38;5;124;03m    The coordinates for a given composition in the PhaseDiagram's basis\u001b[39;00m\n\u001b[1;32m    510\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m    511\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mset\u001b[39m(comp\u001b[38;5;241m.\u001b[39melements) \u001b[38;5;241m-\u001b[39m \u001b[38;5;28mset\u001b[39m(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39melements):\n\u001b[0;32m--> 512\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mcomp\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m has elements not in the phase diagram \u001b[39m\u001b[38;5;132;01m{\u001b[39;00m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m, \u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;241m.\u001b[39mjoin(\u001b[38;5;28mmap\u001b[39m(\u001b[38;5;28mstr\u001b[39m,\u001b[38;5;250m \u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39melements))\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m    513\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m np\u001b[38;5;241m.\u001b[39marray([comp\u001b[38;5;241m.\u001b[39mget_atomic_fraction(el) \u001b[38;5;28;01mfor\u001b[39;00m el \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39melements[\u001b[38;5;241m1\u001b[39m:]])\n",
      "\u001b[0;31mValueError\u001b[0m: Li0.25 Co0.25 O0.5 has elements not in the phase diagram S, O, Li, P"
     ]
    }
   ],
   "source": [
    "plt = interface.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import OrderedDict\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "critical_rxns = [\n",
    "    OrderedDict(\n",
    "        [\n",
    "            (\"Atomic fraction\", round(ratio, 3)),\n",
    "            (\"Reaction equation\", rxn),\n",
    "            (\"E$_{rxt}$ per mol equation (kJ/mol)\", round(rxn_energy, 1)),\n",
    "            (\"E$_{rxt}$ per reactant atom (eV/atom)\", round(reactivity, 3)),\n",
    "        ]\n",
    "    )\n",
    "    for _, ratio, reactivity, rxn, rxn_energy in interface.get_kinks()\n",
    "]\n",
    "interface_reaction_table = pd.DataFrame(critical_rxns)\n",
    "interface_reaction_table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also obtain the mixing ratio between the original compositions, i.e. mol fraction of the first reactant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "interface_reaction_table[\"Molar fraction\"] = pd.Series(\n",
    "    np.round(interface.get_critical_original_kink_ratio(), 3)\n",
    ")\n",
    "interface_reaction_table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the reaction equations are `Reaction` objects suitable for structured analysis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rxn = critical_rxns[2][\"Reaction equation\"]\n",
    "print(rxn)\n",
    "print(type(rxn))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get interface reaction information for reactants LiCoO2 and Li3PS4 in open system to Co.\n",
    "kinks_from_API = mpr.get_interface_reactions(\n",
    "    \"LiCoO2\", \"Li3PS4\", open_el=\"Co\", relative_mu=-1, use_hull_energy=False\n",
    ")\n",
    "\n",
    "# Get inforamtion for the second critical reaction.\n",
    "print(kinks_from_API[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The critical reaction information from REST API should be the same as in the previous table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "critical_rxns_from_API = [\n",
    "    OrderedDict(\n",
    "        [\n",
    "            (\"Atomic fraction\", round(reaction[\"ratio_atomic\"], 3)),\n",
    "            (\"Molar fraction\", round(reaction[\"ratio_molar\"], 3)),\n",
    "            (\"Reaction equation\", reaction[\"rxn_str\"]),\n",
    "            (\n",
    "                \"E$_{rxt}$ per mol equation (kJ/mol)\",\n",
    "                round(float(reaction[\"rxn_energy_sigdig_kJmol\"]), 1),\n",
    "            ),\n",
    "            (\"E$_{rxt}$ per reactant atom (eV/atom)\", round(reaction[\"energy\"], 3)),\n",
    "        ]\n",
    "    )\n",
    "    for reaction in kinks_from_API\n",
    "]\n",
    "pd.DataFrame(critical_rxns_from_API)"
   ]
  }
 ],
 "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
}