{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "624c1da5",
   "metadata": {},
   "source": [
    "# Chemical equilibrium calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32855fc1",
   "metadata": {},
   "source": [
    "This tutorial demonstrates how to use Reaktoro to perform chemical equilibrium calculations. We start by importing the `reaktoro` package:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36cc6cb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "from reaktoro import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0991e598",
   "metadata": {},
   "source": [
    "## Initializing thermodynamic database"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9215f192",
   "metadata": {},
   "source": [
    "Next, we need a thermodynamic database that enables us to compute the thermodynamic properties of species and\n",
    "reactions. For this, we create an object of class [Database](https://reaktoro.org/cpp/classReaktoro_1_1Database.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "499b3392",
   "metadata": {},
   "outputs": [],
   "source": [
    "db = Database(\"supcrt98.xml\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "381471d4",
   "metadata": {},
   "source": [
    "> For more detailed overview on the functionality of this class, please check the tutorial\n",
    "[**Database**](cl.database.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32bfb26a",
   "metadata": {},
   "source": [
    "## Initializing chemical system"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44912b66",
   "metadata": {},
   "source": [
    "To indicate the phases of interest (as well as their species) that may potentially exist at equilibrium,\n",
    "we create an object of class [ChemicalEditor](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalEditor.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "366db1e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "editor = ChemicalEditor(db)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d09a1fd3",
   "metadata": {},
   "source": [
    "We consider an aqueous phase composed of all aqueous species in the database that can be formed with the given\n",
    "chemical elements below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07c61656",
   "metadata": {},
   "outputs": [],
   "source": [
    "editor.addAqueousPhaseWithElements(\"H O Na Cl C Ca Mg Si\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e5ed40e",
   "metadata": {},
   "source": [
    "> **Note:** This automatic selection of chemical species for a phase can result in a large number of them. This\n",
    "potentially increases the computing cost of the chemical reaction calculations. If you are using Reaktoro in\n",
    "demanding applications, you might want to manually specify the chemical species of each phase in your chemical\n",
    "system. This can be achieved by providing an explicit list of species names, e.g.,\n",
    "`editor.addAqueousPhase(\"H2O(l) H+ OH- CO2(aq)\")`. Note, however, that care is required here to ensure relevant\n",
    "species are not missing. The just\n",
    "given example is a bad one in fact, with important species such as `HCO3-` and `CO3--` missing in the list."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a7362f1",
   "metadata": {},
   "source": [
    "We are interested in a gaseous phase containing exactly the following gases (which may not exist in positive\n",
    "amounts at the end of our equilibrium calculation later):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e96501b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "editor.addGaseousPhase(\"H2O(g) CO2(g)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d850cc8",
   "metadata": {},
   "source": [
    "Finally, we consider some pure minerals that could exist in positive amounts in our equilibrium calculations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48fd8ba4",
   "metadata": {},
   "outputs": [],
   "source": [
    "editor.addMineralPhase(\"Halite\")\n",
    "editor.addMineralPhase(\"Calcite\")\n",
    "editor.addMineralPhase(\"Magnesite\")\n",
    "editor.addMineralPhase(\"Dolomite\")\n",
    "editor.addMineralPhase(\"Quartz\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "801cda63",
   "metadata": {},
   "source": [
    "> See the tutorial [**ChemicalEditor**](cl.chemical-editor.ipynb) for studying further capabilities of\n",
    "[ChemicalEditor](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalEditor.html) class."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d4a5f60",
   "metadata": {},
   "source": [
    "## Initializing chemical system"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8080b83",
   "metadata": {},
   "source": [
    "Next, follows an important step with the creation of the chemical system with the information so far collected in the\n",
    "[ChemicalEditor](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalEditor.html) object `editor`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67b49a8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "system = ChemicalSystem(editor)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ac63714",
   "metadata": {},
   "source": [
    "The [ChemicalSystem](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html) class is one of the most\n",
    "important classes in Reaktoro. It is the class used to computationally represent a chemical system, with all its\n",
    "phases, species, and elements. It is also the class used for computation of thermodynamic properties of phases and\n",
    "species, such as activities, chemical potentials, standard Gibbs energies, enthalpies, phase molar volumes, densities,\n",
    "and many others. Many classes in Reaktoro require an instance of\n",
    "[ChemicalSystem](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html) for their initialization,\n",
    "since any chemical calculation needs to know the definition of the chemical system and the thermodynamic models\n",
    "describing the non-ideal behavior of the phases."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39b2586b",
   "metadata": {},
   "source": [
    "> See [**ChemicalSystem**](cl.chemical-system.ipynb) for the explanation on functionality of the class\n",
    "[ChemicalSystem](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a73ea68c",
   "metadata": {},
   "source": [
    "## Initializing equilibrium problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42c91ad2",
   "metadata": {},
   "source": [
    "We use class [EquilibriumProblem](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumProblem.html) to specify the\n",
    "conditions at which our system should be in equilibrium."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67463b40",
   "metadata": {},
   "outputs": [],
   "source": [
    "problem = EquilibriumProblem(system)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c6aea30",
   "metadata": {},
   "source": [
    "Reaktoro provides the class [EquilibriumProblem](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumProblem.html)\n",
    "for the convenient description of equilibrium conditions. Using this class allows one to set the temperature and\n",
    "pressure at equilibrium, and a recipe that describes a mixture of substances and their amounts, which can be seen\n",
    "as initial conditions for the equilibrium calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a62ce246",
   "metadata": {},
   "outputs": [],
   "source": [
    "problem.setTemperature(70, \"celsius\")\n",
    "problem.setPressure(100, \"bar\")\n",
    "problem.add(\"H2O\", 1.0, \"kg\")\n",
    "problem.add(\"CO2\", 2.0, \"mol\")\n",
    "problem.add(\"NaCl\", 1.0, \"mol\")\n",
    "problem.add(\"CaCO3\", 10.0, \"g\")\n",
    "problem.add(\"MgCO3\", 5.0, \"g\")\n",
    "problem.add(\"Quartz\", 1.0, \"mol\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "856d7018",
   "metadata": {},
   "source": [
    "> **Note:** The substance names above can either be chemical formulas, such as `CaCO3` and `CaCl2`, as well as names of\n",
    "species that can be found in the database, such as `Quartz`. Reaktoro will break down the chemical formulas of the\n",
    "substances and calculate the amount of each chemical element in the system. These element amounts are inputs to the\n",
    "equilibrium calculation. In the future, we will only allow species names to be provided since this is a safer way\n",
    "of preventing unfeasible elemental mass conditions to be imposed (e.g., there are *x* moles of C and *y* moles of\n",
    "O, and distributing these among the species always produce an excess of either C or O)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a616809b",
   "metadata": {},
   "source": [
    "## Equilibration of chemical problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af44533d",
   "metadata": {},
   "source": [
    "We now perform a fast Gibbs energy minimization calculation to compute the chemical equilibrium state of the system\n",
    "at given conditions stored in `problem`. For this, we use the convenient function\n",
    "[equilibrate](https://reaktoro.org/cpp/namespaceReaktoro.html#af2d3b39d3e0b8f9cb5a4d9bbb06b697e):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3e895e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "state = equilibrate(problem)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7e459e9",
   "metadata": {},
   "source": [
    "## Analyzing species amounts"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48d7bc6b",
   "metadata": {},
   "source": [
    "The result of the `equilibrate` call before, `state`, is an object of class [ChemicalState](\n",
    "https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html). This object contains the temperature, pressure,\n",
    "and amounts of the species at the computed chemical equilibrium state. We can access these properties as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3229e647",
   "metadata": {},
   "outputs": [],
   "source": [
    "T = state.temperature()\n",
    "P = state.pressure()\n",
    "n = state.speciesAmounts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b979eec",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"T = {T} K\")\n",
    "print(f\"P = {P} Pa\")\n",
    "print(f\"n (in mol) = \\n{n}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82d9bd74",
   "metadata": {},
   "source": [
    "To print the name of each species and its amount (in mol), we execute the following loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf990a68",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Species names : n (in mol)\")\n",
    "for species in system.species():\n",
    "    name = species.name()\n",
    "    amount = state.speciesAmount(name)\n",
    "    print(f\"{name:>13} = {amount}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77b0d758",
   "metadata": {},
   "source": [
    "To examine the complete information about the chemical state, you can also output it to a file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f95aecad",
   "metadata": {},
   "outputs": [],
   "source": [
    "state.output(\"state.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fceb21f4",
   "metadata": {},
   "source": [
    "## Analyzing chemical properties"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73945441",
   "metadata": {},
   "source": [
    "If you require chemical properties of a system that depend on temperature (*T*), pressure (*P*), and composition (*n*),\n",
    "then [ChemicalProperties](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalProperties.html) class is what you need"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1c68d7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "properties = ChemicalProperties(system)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "609164ee",
   "metadata": {},
   "source": [
    "We can compute the chemical properties of the system at the state of equilibrium we found before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1ee8cae",
   "metadata": {},
   "outputs": [],
   "source": [
    "properties.update(T, P, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d303c357",
   "metadata": {},
   "source": [
    "Alternatively, one can also do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4068667",
   "metadata": {},
   "outputs": [],
   "source": [
    "properties = state.properties()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1e05236",
   "metadata": {},
   "source": [
    "> **Note:** The call above creates a new object of [\n",
    "> ChemicalProperties](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalProperties.html) each time. If you are using\n",
    "> Reaktoro in a simulator that needs the chemical properties of the system at millions/billions of states each time\n",
    "> step, instead of populating many instances of ChemicalProperties, prefer to use\n",
    "> [ChemicalProperties::update](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalProperties.html#af923d85484865039fa56889c1a2f36c9)\n",
    "> method of an existing [ChemicalProperties](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalProperties.html)\n",
    "> object."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3115d720",
   "metadata": {},
   "source": [
    "### Analyzing activities"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83f49009",
   "metadata": {},
   "source": [
    "Once we have computed the chemical properties, we can query for some of them. Below we get the natural log of\n",
    "species activities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a68d1ee2",
   "metadata": {},
   "outputs": [],
   "source": [
    "lna = properties.lnActivities().val"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84bc2425",
   "metadata": {},
   "source": [
    "> **Note:** The use of `.ddT`, `.ddP`, and `.ddn`, instead of `.val`, extracts the derivatives of the activities\n",
    "(or any other chemical property) with respect to *T*, *P*, and *n*, respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07d4021a",
   "metadata": {},
   "source": [
    "To compute the actual activities (not their natural log), and print them one by one, we do"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8460cb1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = numpy.exp(lna)\n",
    "print(\"Species names : activities\")\n",
    "for i, species in enumerate(system.species()):\n",
    "    print(f\"{species.name():>13} : {a[i]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0ae96fa",
   "metadata": {},
   "source": [
    "### Analyzing chemical potentials"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee29d4a4",
   "metadata": {},
   "source": [
    "Similarly, we can inspect chemical potentials:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "249fc66b",
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = properties.chemicalPotentials().val\n",
    "print(\"Species names : potentials (in kJ/mol)\")\n",
    "for i, species in enumerate(system.species()):\n",
    "    print(f\"{species.name():>13} : {mu[i]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "311e4014",
   "metadata": {},
   "source": [
    "### Calculating the pH of the aqueous solution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed5e5d3e",
   "metadata": {},
   "source": [
    "Let's create a pH function that computes the pH of the aqueous solution given the chemical properties of the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb92d5a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "evaluate_pH = ChemicalProperty.pH(system)\n",
    "pH = evaluate_pH(properties)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a9e259e",
   "metadata": {},
   "source": [
    "> **Note:** This will be soon simplified!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80a026fb",
   "metadata": {},
   "source": [
    "Besides its value, which can be obtained by  `pH.val`, `ph` also contains derivatives with respect to the species\n",
    "amounts. It can be accessed by `pH.ddn`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53f7b122",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "print(f\"The pH of the aqueous phase is {pH.val}.\")\n",
    "print(f\"Its sensitivity with respect to speciation is:\")\n",
    "print(\"Species names : ∂(pH)/∂n (in 1/mol)\")\n",
    "for i, species in enumerate(system.species()):\n",
    "    print(f\"{species.name():>13} : {pH.ddn[i]}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}