SurfacePourbaixTheory.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Theory for Constructing Mg Surface Pourbaix Diagram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* This Notebook is based on this [paper](https://journals.aps.org/prapplied/abstract/10.1103/PhysRevApplied.1.014001) which depicts the approach for generating the Pourbaix diagram.\n",
    "\n",
    "* We construct the Pourbaix diagram of multiple Mg surfaces, not only focusing on the (0001) basal plane, but also on the (10m10) and (11m20) surfaces.\n",
    "\n",
    "* Not only different on-surface adsorption surfaces are considered, but also the subsurface adsorption surfaces are also included here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. The most essential lessons we could learn from the aforementioned paper are the following two equations:\n",
    "\n",
    "    * <span style=\"color:orange\"> $\\Delta G_{\\alpha} = \\Delta G^{\\circ}_{\\alpha} - \\Sigma n_i \\mu_i$</span>, which depicts the surface adsorption (excess) free energy with respect to the corresponding free energy at the standard condition. \n",
    "      \n",
    "      Here $\\mu_i = \\tilde{\\mu}_i - \\tilde{\\mu}_i^{\\circ}$ refers to the calibrated chemical potential of speices $i$.\n",
    "\n",
    "    * In terms of water system (where $\\mu_H$ and $\\mu_O$ matter), the calibatered chemical potential of H and O are important:\n",
    "\n",
    "      <span style=\"color:cyan\"> $\\mu_H = \\tilde{\\mu}_H - \\tilde{\\mu}_H^{\\circ} = (\\mu_e - \\mu_e^{SHE}) - k_BT~ln10~pH$ </span> and \n",
    "\n",
    "      <span style=\"color:pink\"> $\\mu_O = \\Delta G_f^{\\circ}(H_2O) - 2\\mu_H$ </span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Consider an arbitrary Mg surface in any adsorption state, which has the chemical symbol as $Mg^{\\beta}*O_xH_y$ ($\\beta$ refers to different Mg surfaces)\n",
    "\n",
    "   * The adsorption free energy could be expressed as: <span style=\"color:orange\">$\\Delta G[Mg^{\\beta}*O_xH_y] = \\Delta G^{\\circ}[Mg^{\\beta}*O_xH_y] - \\Sigma_i n_i \\mu_i$</span>\n",
    "  \n",
    "   * The first term <span style=\"color:red\">*$\\Delta G^{\\circ}[Mg^{\\beta}*O_xH_y]$*</span> refers to the standard adsorption free energy, which equals to:\n",
    "\n",
    "      $\\Delta G^{\\circ}[Mg^{\\beta}*O_xH_y] = E[Mg^{\\beta}*O_xH_y] - N_{Mg}\\tilde{\\mu}^{\\circ}_{Mg} - x\\tilde{\\mu}^{\\circ}_O - y\\tilde{\\mu}^{\\circ}_H$\n",
    "\n",
    "      $= E[Mg^{\\beta}*O_xH_y] - N_{Mg} * E_{Mg}^{Bulk} - \\frac{x}{2}E_{O_2} - \\frac{y}{2}E_{H_2}$ \n",
    "\n",
    "   * Then consider the chemical potential biased term:\n",
    "\n",
    "      <span style=\"color:cyan\">$\\Sigma_i n_i \\mu_i = N_{Mg}\\mu_{Mg} + x\\mu_O + y\\mu_H = N_{Mg}(E_{Mg}^{\\beta} - E_{Mg}^{Bulk}) + x\\Delta G_f^{\\circ}(H_2O) + (y-2x)\\mu_H$</span>\n",
    "\n",
    "   * The final expression of the arbitrary Mg surface is \n",
    "      \n",
    "      $\\Delta G[Mg^{\\beta}*O_xH_y] = E[Mg^{\\beta}*O_xH_y] - E[Mg^{\\beta}*] - \\frac{x}{2}E_{O_2} - \\frac{y}{2}E_{H_2} - x\\Delta G_f^{\\circ}(H_2O)$\n",
    "\n",
    "      $+ (2x -y) [(\\mu_e - \\mu_e^{SHE}) - k_BT~ln10~pH]$\n",
    "\n",
    "   * Consider the Fermi energy basis (potential shift), higher chemcial potential of electrons refers to a negative potential (and vice versa), \n",
    "  \n",
    "     thus <span style=\"color:pink\">$U_{SHE} = (\\mu_e^{SHE} - \\mu_e)/e$</span> and <span style=\"color:yellow\">$\\mu_e - \\mu_e^{SHE} = -eU_{SHE}$</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. In terms of the symmetric adsorption structure, the previous equation should be slightly modified with considering the prefactor 2:\n",
    "\n",
    "   * We hope to study the adsorption state of a Mg surface as $Mg^{\\beta}*O_xH_y$, but construct a symmetric system with a chemical symbol as $Mg^{\\beta}*O_{2x}H_{2y}$.\n",
    "\n",
    "   * For this symmetric system, nothing sepcial will change for the final free energy expression, except the prefactors:\n",
    "\n",
    "      $\\Delta G[Mg^{\\beta}*O_{2x}H_{2y}] = E[Mg^{\\beta}*O_{2x}H_{2y}] - E[Mg^{\\beta}*] - xE_{O_2} - yE_{H_2} - 2x\\Delta G_f^{\\circ}(H_2O)$\n",
    "\n",
    "      $+ 2(2x -y) [(\\mu_e - \\mu_e^{SHE}) - k_BT~ln10~pH]$\n",
    "\n",
    "   * Then we just need half of this surface excess free energy from:\n",
    "\n",
    "      $\\Delta G[Mg^{\\beta}*O_{x}H_{y}] = \\frac{1}{2}\\Delta G[Mg^{\\beta}*O_{2x}H_{2y}]= $ <span style=\"color:pink\"> $\\frac{E[Mg^{\\beta}*O_{2x}H_{2y}] - E[Mg^{\\beta}*]}{2}$</span> $ - \\frac{x}{2}E_{O_2} - \\frac{y}{2}E_{H_2} - x\\Delta G_f^{\\circ}(H_2O)$\n",
    "\n",
    "      $+ (2x -y) [(\\mu_e - \\mu_e^{SHE}) - k_BT~ln10~pH]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. <span style=\"color:cyan\">How about a slab within a continuum media (implicit solvent model)?</span>\n",
    "\n",
    "   * Initial feel of the contribution of this solvation system to our surface Pourbaix model is to replace the energy of the pristine slab (in vacuum) and the energy of the complex slab (in vacuum) with the energies in solvent.\n",
    "\n",
    "   * Recall the way we determine $\\mu_H$, we consider the ionization process of the hydrogen atom happens in the vacuum (reservoir), then the gas phase hydronium is solvated in the liquid phase (where we include the dissolving energy). If this is the case, we don not need to consider the solvation effect of the reference gas molecule in the implicit solvent case.\n",
    "\n",
    "   * To sum up, in this implicit slab model, we consider the reservoir is still in vacuum, but our slab is fully covered by the implicit solvent:\n",
    "\n",
    "      $\\Delta G^{sol}[Mg^{\\beta}*O_{x}H_{y}] = $ <span style=\"color:yellow\"> $\\frac{E^{sol}[Mg^{\\beta}*O_{2x}H_{2y}] - E^{sol}[Mg^{\\beta}*]}{2}$</span> $ - \\frac{x}{2}E_{O_2} - \\frac{y}{2}E_{H_2} - x\\Delta G_f^{\\circ}(H_2O)$\n",
    "\n",
    "      $+ (2x -y) [(\\mu_e - \\mu_e^{SHE}) - k_BT~ln10~pH]$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Revise the Implicit Solvation Model implemented in the VASPsol\n",
    "   \n",
    "   * This part is basically based on this [paper](https://aip.scitation.org/doi/10.1063/1.4865107), there are also some other useful papers, such as this [paper](https://aip.scitation.org/doi/10.1063/1.5132354).\n",
    "  \n",
    "   * The free energy *A* of the combined solute/solvent system could be expressed as a sum of two terms, a universal functional *F* of the total electron density and the thermodynamically avearged atomic densities of the <span style=\"color:red\">solvent species ?</span> $\\left\\{N_i(\\vec{r})\\right\\}$, and a term describing the electrostatic energy contribution:\n",
    "  \n",
    "       <span style=\"color:orange\">$A=F\\left[n_{\\mathrm{tot}},\\left\\{N_i(\\vec{r})\\right\\}\\right]+\\int d^3 r V_{\\mathrm{ext}}(\\vec{r})\\left(\\sum_i Z_i N_i(\\vec{r})-n_{\\mathrm{tot}}(\\vec{r})\\right)$</span>\n",
    "\n",
    "       with $n_{\\mathrm{tot}} = n_{\\mathrm{solute}} + n_{\\mathrm{solv}}$, and $V_{\\mathrm{ext}}(\\vec{r})$ is the external potential due to the solvent nuclei.\n",
    "\n",
    "    * Using the variational principle for the free energy of the electron-nuclear system at a fixed external potential to find the ground state free energy of the system:\n",
    "\n",
    "       <span style=\"color:orange\">$A_0=\\min _{n_{\\mathrm{tot}},\\left\\{N_i(\\vec{r})\\right\\}}\\left\\{F\\left[n_{\\mathrm{tot}},\\left\\{N_i(\\vec{r})\\right\\}\\right]\\right.\\left.+\\int d^3 r V_{\\mathrm{ext}}(\\vec{r})\\left(\\sum_i Z_i N_i(\\vec{r})-n_{\\mathrm{tot}}(\\vec{r})\\right)\\right\\}$</span>\n",
    "\n",
    "   * To make it amenable to a computational treatment, the free energy is first minimized over the solvent electron density and then over the solute electron density.\n",
    "\n",
    "   * Firstly minimizing the free energy with respect to the solvent electron density (erase the variation degree of the solvent electron density), we obtain\n",
    "      \n",
    "       $\\tilde{A}=  G\\left[n_{\\text {solute }}(\\vec{r}),\\left\\{N_i(\\vec{r})\\right\\}, V_{\\text {ext }}(\\vec{r})\\right]  -\\int d^3 r V_{\\text {ext }}(\\vec{r}) n_{\\text {solute }}(\\vec{r})$\n",
    "       where\n",
    "       \n",
    "       <span style=\"color:cyan\"> $G {\\left[n_{\\text {solute }}(\\vec{r}),\\left\\{N_i(\\vec{r})\\right\\}, V_{\\mathrm{ext}}(\\vec{r})\\right] } =\\min _{n_{\\text {solv }}}\\left\\{F\\left[n_{\\mathrm{tot}},\\left\\{N_i(\\vec{r})\\right\\}\\right]\\right. \\left.+\\int d^3 r V_{\\mathrm{ext}}(\\vec{r})\\left(\\sum_i Z_i N_i(\\vec{r})-n_{\\text {solv }}(\\vec{r})\\right)\\right\\}$ </span>\n",
    "\n",
    "    * The functional G could be decomposed into two parts: part $i)$ $A_{KS}$ is the usual KS density functional for the solute and part $ii)$ $A_{diel}$ is the term that encapsulates all the interactions of the solute with the solvent and the internal energy of the solvent.\n",
    "\n",
    "       <span style=\"color:cyan\">$\\begin{aligned}\n",
    "        & G\\left[n_{\\text {solute }}(\\vec{r}),\\left\\{N_i(\\vec{r})\\right\\}, V_{\\text {ext }}(\\vec{r})\\right] \\\\\n",
    "        &= A_{\\mathrm{KS}}\\left[n_{\\text {solute }}(\\vec{r}), V_{\\text {ext }}(\\vec{r})\\right] \\\\\n",
    "        &+A_{\\text {diel }}\\left[n_{\\text {solute }}(\\vec{r}),\\left\\{N_i(\\vec{r})\\right\\}, V_{\\text {ext }}(\\vec{r})\\right]\n",
    "        \\end{aligned}$</span>\n",
    "\n",
    "   * To further simplify the free energy expression, the $A_{diel}$ could be minized with respect to the average atomic densities of the solvent $\\left\\{N_i(\\vec{r})\\right\\}$. Finally, the ground state free energy of the solute/solvent combined system could be expressed as:\n",
    "    \n",
    "       <span style=\"color:orange\">$ A_0= \\min _{n_{\\text {solue }}(\\vec{r})}\\left\\{A_{\\mathrm{KS}}\\left[n_{\\text {solute }}(\\vec{r}), V_{\\text {ext }}(\\vec{r})\\right]\\right.\\left.-\\int d^3 r V_{\\mathrm{ext}}(\\vec{r}) n_{\\text {solute }}(\\vec{r})+\\tilde{A}_{\\text {diel }}\\left[n_{\\text {solute }}(\\vec{r}), V_{\\mathrm{ext}}(\\vec{r})\\right]\\right\\}$</span>\n",
    "\n",
    "   * All the solvent effect are contained in the functional $\\tilde{A}_{\\text {diel}}$, which is obtained by the minimazation over the solvent electron density and the thermodynamically average atomic densities of the solvent. Therefore, the functional $\\tilde{A}_{\\text {diel}}$ describe a continuum model for the solvent, whose equilibrium structure is fully determined by the properties of the solute.\n",
    "\n",
    "   * Approximation on the <span style=\"color:red\">$\\tilde{A}_{\\text {diel}}$</span>: $i)$ electrostaic interaction between the solute and the solvent affects the equilibrium polarization of the solvent dipoles. And the solvent polarization can be described by the local relative permittivity of the solvent $\\epsilon(\\vec{r})$. $ii)$ cavtion and dispersion may also play important roles, $A_{\\mathrm{cav}}=\\tau \\int d^3 r|\\nabla S|$, with $\\tau$ as the effective surface tension parameter and $S(\\vec{r})$ is the cavity shape function as $S\\left(n_{\\text {solute }}(\\vec{r})\\right)=\\frac{1}{2} \\operatorname{erfc}\\left\\{\\frac{\\log \\left(n_{\\text {solute }} / n_{\\mathrm{c}}\\right)}{\\sigma \\sqrt{2}}\\right\\}$.\n",
    "\n",
    "   * Decoupling the electrostatic term from the KS functional and combining it within the interaction term and the cavitation term ($\\phi(\\vec{r})$ is the combined electrostatic potential due to the electronic and nuclear charges of the soluye system in a polarizable medium):\n",
    "  \n",
    "       <span style=\"color:orange\">$A\\left[n_{\\text {solute }}(\\vec{r}), \\phi(\\vec{r})\\right]= A_{\\mathrm{TXC}}\\left[n_{\\text {solute }}(\\vec{r})\\right]  +\\int d^3 r \\phi(\\vec{r})\\left(N_{\\text {solute }}(\\vec{r})-n_{\\text {solute }}(\\vec{r})\\right) -\\int d^3 r \\epsilon(\\vec{r}) \\frac{|\\nabla \\phi|^2}{8 \\pi} +A_{\\text{cav }}$</span>\n",
    "\n",
    "   * Mapping the relative permittivity from the real space to the electronic density of the solute, a diffuse cavity functional is determined as:\n",
    "\n",
    "       $\\epsilon\\left(n_{\\text {solute }}(\\vec{r})\\right)=1+\\left(\\epsilon_{\\mathrm{b}}-1\\right) S\\left(n_{\\text {solute }}(\\vec{r})\\right)$\n",
    "\n",
    "   * First we can minimise the $A\\left[n_{\\text {solute }}(\\vec{r}),\\phi(\\vec{r})\\right]$ with respect to $\\phi(\\vec{r})$, which will lead to a generalized Poisson eqaution:\n",
    "\n",
    "       <span style=\"color:cyan\">$\\nabla \\cdot {\\left[\\epsilon\\left(n_{\\text {solute }}(\\vec{r})\\right) \\nabla \\phi\\right] } =-4 \\pi\\left\\{N_{\\text {solute }}(\\vec{r})-n_{\\text {solute }}(\\vec{r})\\right\\}$</span>\n",
    "\n",
    "    * If minimise the $A\\left[n_{\\text {solute }}(\\vec{r}),\\phi(\\vec{r})\\right]$ with respect to $n_{\\text {solute}}$, on top of the typcial KS local part potential, we have the following terms:\n",
    "\n",
    "       $V_{\\text {solv }}=\\frac{d \\epsilon\\left(n_{\\text {solute }}(\\vec{r})\\right)}{d n_{\\text {solute }}(\\vec{r})} \\frac{|\\nabla \\phi|^2}{8 \\pi}+\\tau \\frac{d|\\nabla S|}{d n_{\\text {solute }}(\\vec{r})}$\n",
    "\n",
    "    * The expression of the $A\\left[n_{\\text {solute }}(\\vec{r}),\\phi(\\vec{r})\\right]$ is then expanded by the VASPsol author, with the inclusion of the interaction between the electrolyte(solvent) ions with the solute electrostatic potential $\\phi(\\vec{r})$, and a nonelectrostatic contribution to the free energy from the mobile ions in the electrolyte ($A_\\text{ion} = k_B T S_{\\text {ion }}=-\\int \\sum_i c_i \\frac{z_i q \\phi}{k_B T} d^3 r$):\n",
    "\n",
    "        <span style=\"color:cyan\">$\\begin{aligned}\n",
    "        A[n(\\vec{r}), \\phi(\\vec{r})]= & A_{\\mathrm{TXC}}[n(\\vec{r})]+\\int \\phi(\\vec{r}) \\rho_s(\\vec{r}) d^3 r-\\int \\epsilon(\\vec{r}) \\frac{|\\nabla \\phi|^2}{8 \\pi} d^3 r \\\\\n",
    "        & +\\int \\frac{1}{2} \\phi(\\vec{r}) \\rho_{\\mathrm{ion}}(\\vec{r}) d^3 r+A_{\\mathrm{cav}}+A_{\\mathrm{ion}}\n",
    "        \\end{aligned}$</span>\n",
    "\n",
    "        with $\\rho_s(\\vec{r})$ as the sum of the solute electronic and nuclear charge densities $\\rho_{\\mathrm{s}}(\\vec{r})=n(\\vec{r})+N(\\vec{r})$ and the ion charge density of the electrolyte as $\\rho_{\\text {ion }}(\\vec{r})=\\sum_i q z_i c_i(\\vec{r})$.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.7 ('work')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.9.7"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "2f96c6429f727abf40bda0f738b748f4f96354aa996dc3cc220cc83619beaf19"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}