{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Balance simple reactions\n",
    "Here we present eleven different sets of reactants and products, not necessarily based on real reactions.\n",
    "The molecules can be passed as formulas, or as PLAMS Molecule objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scm.plams.tools.reaction import ReactionEquation, balance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Aspirin in water"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 C9H8O4 + 1 H2O => 1 C2H4O2 + 1 C7H6O3\n"
     ]
    }
   ],
   "source": [
    "reactants = [\"C9H8O4\", \"H2O\"]\n",
    "products = [\"C2H4O2\", \"C7H6O3\"]\n",
    "\n",
    "reaction = balance(reactants, products)\n",
    "print(reaction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "``reaction`` is of type ``ReactionEquation`` and has the following attributes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reaction.coeffs=array([1, 1, 1, 1])\n",
      "reaction.message='Success'\n",
      "reaction.reactants=['C9H8O4', 'H2O']\n",
      "reaction.products=['C2H4O2', 'C7H6O3']\n"
     ]
    }
   ],
   "source": [
    "print(f\"{reaction.coeffs=}\")\n",
    "print(f\"{reaction.message=}\")\n",
    "print(f\"{reaction.reactants=}\")\n",
    "print(f\"{reaction.products=}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Specify reactants/products as Species\n",
    "\n",
    "An alternative to the string representation of reactants and products is to use Species, which also allows you to set the molecular charges and minimum coefficients (see later examples):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 C9H8O4 + 1 H2O => 1 C2H4O2 + 1 C7H6O3\n"
     ]
    }
   ],
   "source": [
    "from scm.plams.tools.reaction import Species as S\n",
    "\n",
    "reactants = [S(\"C9H8O4\", charge=0), S(\"H2O\", charge=0)]\n",
    "products = [S(\"C2H4O2\", charge=0), S(\"C7H6O3\", charge=0)]\n",
    "reaction = balance(reactants, products)\n",
    "print(reaction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initialize species from SMILES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 CH4 + 2 O2 => 1 CO2 + 2 H2O\n"
     ]
    }
   ],
   "source": [
    "from scm.plams.tools.reaction import Species as S\n",
    "\n",
    "reactants = [S.from_smiles(\"C\"), S.from_smiles(\"O=O\")]\n",
    "products = [S.from_smiles(\"O=C=O\"), S.from_smiles(\"O\")]\n",
    "reaction = balance(reactants, products)\n",
    "print(reaction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also format the reaction with the SMILES strings (requires that all reactants and products have the ``smiles`` attribute set):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 C + 2 O=O => 1 O=C=O + 2 O\n"
     ]
    }
   ],
   "source": [
    "print(f\"{reaction:smiles}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Superfluous reactants and products (multiple possible solutions)\n",
    "\n",
    "When there are multiple possible solutions, a single one is chosen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 C9H8O4 + 1 H2O => 1 C2H4O2 + 1 C7H6O3\n",
      "Some coefficients are zero:\n",
      "reaction.coeffs=array([1, 1, 0, 0, 1, 1, 0])\n"
     ]
    }
   ],
   "source": [
    "reactants = [\"C9H8O4\", \"H2O\", \"HO\", \"HO2\"]\n",
    "products = [\"C2H4O2\", \"C7H6O3\", \"H2O2\"]\n",
    "reaction = balance(reactants, products)\n",
    "print(reaction)\n",
    "print(\"Some coefficients are zero:\")\n",
    "print(f\"{reaction.coeffs=}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exampe that is not easily balanced"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9 C9O4H6 + 25 OH => 23 C2O2H3 + 5 C7O3H2\n"
     ]
    }
   ],
   "source": [
    "reactants = [\"C9O4H6\", \"OH\"]\n",
    "products = [\"C2O2H3\", \"C7O3H2\"]\n",
    "reaction = balance(reactants, products)\n",
    "print(reaction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reaction that cannot be balanced"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inconsistent system: No solution\n",
      "reaction.coeffs=None\n",
      "reaction.message='Empty nullspace'\n"
     ]
    }
   ],
   "source": [
    "reactants = [\"FeS2\", \"HNO3\"]\n",
    "products = [\"Fe2S3O12\", \"NO\", \"H2SO4\"]\n",
    "reaction = balance(reactants, products)\n",
    "print(reaction)\n",
    "print(f\"{reaction.coeffs=}\")\n",
    "print(f\"{reaction.message=}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sympy method\n",
    "\n",
    "By default, a native method is used to compute the nullspace. Optionally, this can be done with sympy. The result should be the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 C9H8O4 + 1 H2O => 1 C2H4O2 + 1 C7H6O3\n"
     ]
    }
   ],
   "source": [
    "reactants = [\"C9H8O4\", \"H2O\"]\n",
    "products = [\"C2H4O2\", \"C7H6O3\"]\n",
    "\n",
    "reaction = balance(reactants, products, method=\"sympy\")\n",
    "print(reaction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Balance reactions with charged species\n",
    "\n",
    "Charges are printed within square brackets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 OH[-1] + 1 H3O[+1] => 2 H2O\n"
     ]
    }
   ],
   "source": [
    "from scm.plams.tools.reaction import Species as S\n",
    "\n",
    "reactants = [S(\"OH\", charge=-1), S(\"H3O\", charge=+1)]\n",
    "products = [S(\"H2O\", charge=0)]\n",
    "reaction = balance(reactants, products)\n",
    "print(reaction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the charges cannot be balanced, no solution will be found:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inconsistent system: No solution\n"
     ]
    }
   ],
   "source": [
    "from scm.plams.tools.reaction import Species as S\n",
    "\n",
    "reactants = [S(\"OH\", charge=-1), S(\"H3O\", charge=0)]  # \"neutral\" H3O!\n",
    "products = [S(\"H2O\", charge=0)]\n",
    "reaction = balance(reactants, products)\n",
    "print(reaction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The charges are also inferred from the SMILES strings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 H[+1] + 1 HO[-1] => 1 H2O\n"
     ]
    }
   ],
   "source": [
    "from scm.plams.tools.reaction import Species as S\n",
    "\n",
    "reactants = [S.from_smiles(\"[H+]\"), S.from_smiles(\"[OH-]\")]\n",
    "products = [S.from_smiles(\"O\")]\n",
    "reaction = balance(reactants, products)\n",
    "print(reaction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 [H+] + 1 [OH-] => 1 O\n"
     ]
    }
   ],
   "source": [
    "print(f\"{reaction:smiles}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reaction with many possible products: Set minimum coefficients\n",
    "\n",
    "You can set minimum values for the coefficients with the ``min_coeff`` attribute for a ``Species``. This is useful if you have a list of many possible products, and want to find a balanced reaction that includes one or more specific products:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6 H2O + 1 C9H8O4 => 1 C7H6O3 + 3 H2O2 + 1 CH4O + 1 CH4\n"
     ]
    }
   ],
   "source": [
    "from scm.plams.tools.reaction import Species as S\n",
    "\n",
    "reactant_smiles = [\"O\", \"CC(=O)Oc1ccccc1C(=O)O\"]\n",
    "product_smiles = [\n",
    "    \"CC(=O)O\",\n",
    "    \"O=C(O)c1ccccc1O\",\n",
    "    \"OO\",\n",
    "    \"O=C(O)O\",\n",
    "    \"O=C=O\",\n",
    "    \"CO\",\n",
    "    \"C\",\n",
    "    \"O=C(O)C1=CC(O)C=CC1=O\",\n",
    "]\n",
    "\n",
    "reactants = [S.from_smiles(x, min_coeff=1) for x in reactant_smiles]  # both reactants must appear\n",
    "products = [S.from_smiles(x) for x in product_smiles]\n",
    "products[1].min_coeff = 1  # the second product must appear\n",
    "products[2].min_coeff = 3  # the third product must have a coefficient of at least 3\n",
    "reaction = balance(reactants, products)\n",
    "print(reaction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6 O + 1 CC(=O)Oc1ccccc1C(=O)O => 1 O=C(O)c1ccccc1O + 3 OO + 1 CO + 1 C\n"
     ]
    }
   ],
   "source": [
    "print(f\"{reaction:smiles}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Direct usage of the ReactionEquation class\n",
    "\n",
    "You can also directly use the ReactionEquation class. This can be useful and save some time if you want to loop over many possible values of ``min_coeffs``, for example.\n",
    "\n",
    "In this case, do not use ``Species`` but provide formulas or PLAMS Molecules directly, use the ``set_charges`` method to set charges, and set the ``min_coeffs`` array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "\n",
    "reactant_smiles = [\"O\", \"CC(=O)Oc1ccccc1C(=O)O\"]\n",
    "product_smiles = [\n",
    "    \"CC(=O)O\",\n",
    "    \"O=C(O)c1ccccc1O\",\n",
    "    \"OO\",\n",
    "    \"O=C(O)O\",\n",
    "    \"O=C=O\",\n",
    "    \"CO\",\n",
    "    \"C\",\n",
    "    \"O=C(O)C1=CC(O)C=CC1=O\",\n",
    "    \"O=C(O)C1=CC(O)C(O)C=C1O\",\n",
    "    \"O=C(O)C(C=CO)=C(O)C=CO\",\n",
    "    \"O=C(CO)Oc1ccccc1\",\n",
    "    \"Oc1ccccc1\",\n",
    "    \"O=C1C=CC=CC1\",\n",
    "    \"OC1=CC(O)C=CC1\",\n",
    "    \"C=CC=CC=C=O\",\n",
    "    \"C1=C=CC=CC=1\",\n",
    "    \"O=C=CO\",\n",
    "    \"O=CO\",\n",
    "    \"O=C=C(O)O\",\n",
    "    \"O=C=C=O\",\n",
    "    \"O=C1OC1=O\",\n",
    "    \"O=COc1ccccc1\",\n",
    "    \"O=C1CC=CC(O)C1\",\n",
    "    \"C=CC(O)C(O)C=C=O\",\n",
    "    \"OC1=CC=CC(O)C1=C(O)O\",\n",
    "    \"O=C(O)C12C(=O)C1C=CC2O\",\n",
    "    \"O=C(O)CO\",\n",
    "    \"O=CC(=O)O\",\n",
    "    \"O=C1C=CC(O)=C(C1)C(=O)O\",\n",
    "    \"C=CC=CC(=O)O\",\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting loop over products..\n",
      "                  CC(=O)O: 1 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 CC(=O)O + 1 O=C(O)C1=CC=CC=C1O\n",
      "          O=C(O)c1ccccc1O: 1 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 CC(=O)O + 1 O=C(O)C1=CC=CC=C1O\n",
      "                       OO: 2 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 CC(=O)O + 1 OO + 1 O=COC1C=CC=CC=1\n",
      "                  O=C(O)O: 1 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 CC(=O)O + 1 O=C(O)O + 1 C1=C=CC=CC=1\n",
      "                    O=C=O: 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 CC(=O)O + 1 O=C=O + 1 C1=C=CC=CC=1\n",
      "                       CO: 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 CO + 1 C1=C=CC=CC=1 + 1 O=C1OC1=O\n",
      "                        C: 1 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 O=C(O)C1=CC=CC=C1O + 1 O=C=O + 1 C\n",
      "    O=C(O)C1=CC(O)C=CC1=O: 2 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 C + 1 O=C(O)C1=CC(O)C=CC1=O + 1 O=CO\n",
      "  O=C(O)C1=CC(O)C(O)C=C1O: 3 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 C + 1 O=C(O)C1=CC(O)C(O)C=C1O + 1 O=CO\n",
      "   O=C(O)C(C=CO)=C(O)C=CO: 3 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 C + 1 O=C(O)C(C=CO)=C(O)C=CO + 1 O=CO\n",
      "         O=C(CO)Oc1ccccc1: 1 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 O=C(CO)OC1C=CC=CC=1 + 1 O=CO\n",
      "                Oc1ccccc1: 1 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 CC(=O)O + 1 O=C=O + 1 OC1C=CC=CC=1\n",
      "             O=C1C=CC=CC1: 1 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 CC(=O)O + 1 O=C=O + 1 O=C1C=CC=CC1\n",
      "           OC1=CC(O)C=CC1: 2 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 CC(=O)O + 1 O=C=O + 1 OC1=CC(O)C=CC1\n",
      "              C=CC=CC=C=O: 1 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 CC(=O)O + 1 O=C=O + 1 C=CC=CC=C=O\n",
      "             C1=C=CC=CC=1: 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 CC(=O)O + 1 O=C=O + 1 C1=C=CC=CC=1\n",
      "                   O=C=CO: 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 O=C=CO + 1 O=COC1C=CC=CC=1\n",
      "                     O=CO: 1 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 O=C(CO)OC1C=CC=CC=1 + 1 O=CO\n",
      "                O=C=C(O)O: 1 O + 2 CC(=O)OC1=CC=CC=C1C(=O)O => 2 O=C(CO)OC1C=CC=CC=1 + 1 O=C=C(O)O\n",
      "                  O=C=C=O: 2 CC(=O)OC1=CC=CC=C1C(=O)O => 2 O=C(CO)OC1C=CC=CC=1 + 1 O=C=C=O\n",
      "                O=C1OC1=O: 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 CO + 1 C1=C=CC=CC=1 + 1 O=C1OC1=O\n",
      "             O=COc1ccccc1: 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 O=C=CO + 1 O=COC1C=CC=CC=1\n",
      "           O=C1CC=CC(O)C1: 2 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 CC(=O)O + 1 O=C=O + 1 O=C1CC=CC(O)C1\n",
      "         C=CC(O)C(O)C=C=O: 2 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 CO + 1 O=C=C=O + 1 C=CC(O)C(O)C=C=O\n",
      "     OC1=CC=CC(O)C1=C(O)O: 2 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 CC(=O)O + 1 OC1=CC=CC(O)C1=C(O)O\n",
      "   O=C(O)C12C(=O)C1C=CC2O: 2 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 C + 1 O=CO + 1 O=C(O)C12C(=O)C1C=CC2O\n",
      "                 O=C(O)CO: 1 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 O=COC1C=CC=CC=1 + 1 O=C(O)CO\n",
      "                O=CC(=O)O: 1 O + 2 CC(=O)OC1=CC=CC=C1C(=O)O => 2 O=C(CO)OC1C=CC=CC=1 + 1 O=CC(=O)O\n",
      "  O=C1C=CC(O)=C(C1)C(=O)O: 2 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 C + 1 O=CO + 1 O=C1C=CC(O)=C(C1)C(=O)O\n",
      "             C=CC=CC(=O)O: 2 O + 1 CC(=O)OC1=CC=CC=C1C(=O)O => 1 CC(=O)O + 1 O=C=CO + 1 C=CC=CC(=O)O\n"
     ]
    }
   ],
   "source": [
    "from scm.plams import from_smiles\n",
    "\n",
    "reactants = [from_smiles(smiles) for smiles in reactant_smiles]  # PLAMS Molecules\n",
    "products = [from_smiles(smiles) for smiles in product_smiles]  # PLAMS Molecules\n",
    "\n",
    "# Create the Reaction object with all the molecules\n",
    "reaction = ReactionEquation(reactants, products)\n",
    "\n",
    "print(\"Starting loop over products..\")\n",
    "nmols = len(reactants) + len(products)\n",
    "nreactants = len(reactants)\n",
    "for iprod, _ in enumerate(products):\n",
    "    print(f\"{product_smiles[iprod]:>25s}: \", end=\"\")\n",
    "    min_coeffs = numpy.zeros(nmols)\n",
    "    min_coeffs[nreactants + iprod] = 1\n",
    "    reaction.balance(min_coeffs=min_coeffs)\n",
    "    print(f\"{reaction:smiles}\")"
   ]
  }
 ],
 "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}