{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bfab7943",
   "metadata": {},
   "source": [
    "# Performing a chemical equilibrium calculation with customized activity models\n",
    "\n",
    "This tutorial demonstrates how to use Reaktoro to perform a chemical equilibrium calculation with customized activity\n",
    "models. We use an example of equilibration of  H<sub>2</sub>O-NaCl-CO<sub>2</sub> system when 1 kg of H<sub>2</sub>O,\n",
    "100 g of CO<sub>2</sub>, and 1 mol of NaCl are mixed at temperature 60 &deg;C and pressure 300 bar.\n",
    "First, we import everything from the `reaktoro` package by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99fc7009",
   "metadata": {},
   "outputs": [],
   "source": [
    "from reaktoro import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1a020b2",
   "metadata": {},
   "source": [
    "To indicate phases and corresponding to them species, as well as models used to evaluate activities of the aqueous\n",
    "species, we use the class [ChemicalEditor](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalEditor.html). Here,\n",
    "the default database SUPCRT98 is used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3713c2ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "editor = ChemicalEditor()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82dbf412",
   "metadata": {},
   "source": [
    "## Specifying different phases and corresponding models\n",
    "\n",
    "Specifying the phases and their species is not enough to fully describe a chemical system in the *computational\n",
    "sense*. Every phase in Reaktoro has two associated models: a *thermodynamic model* and a *chemical model*. These\n",
    "denominations are not standard in the literature, but they are useful in the differentiation of two needed types\n",
    "of models for a phase.\n",
    "\n",
    "* A *thermodynamic model* is a model for the calculation of *standard thermodynamic properties* of the species.\n",
    "Examples include standard Gibbs energies, or standard chemical potentials, standard molar volumes, standard heat\n",
    "capacities, standard enthalpies, and so forth. These models are functions of *temperature* and *pressure* only.\n",
    "\n",
    "* A *chemical model* is a model that describes the *non-ideal behavior* of phases. These models not only depend on\n",
    "temperature and pressure, like the thermodynamic models, but also on the amounts of the species in the phase. To be\n",
    "more precise, on the concentrations of these species, which can be calculated from the amounts of the species.\n",
    "\n",
    "One can define a chemical system with many phases, each phase containing one or more species. This does not mean\n",
    "that all phases and their species exist at positive amounts! To be precise, this means that the chemical\n",
    "calculations, equilibrium or kinetics, are capable of deciding if a phase in the chemical system should exist at\n",
    "positive amounts for some given conditions (e.g., temperature, pressure, overall composition).\n",
    "\n",
    "By selecting as many phases as possible, with the possibilities constrained by the *thermodynamic database* being\n",
    "used, one can increase the confidence level of the estimated chemical states. Note, however, that accurate and\n",
    "realistic estimates depend on many more factors than just the selection of potential phases, such as the *choice of\n",
    "thermodynamic models for non-ideal phases*. Furthermore, note that adding too many phases and species to the\n",
    "definition of the chemical system can result in *more computationally expensive* chemical calculations. In critical\n",
    "performance applications, for instance, when combining chemical reactions and fluid flow and species transport\n",
    "modeling, restricting the number of phases and species might be necessary for achieving feasible simulation times.\n",
    "The modeler is responsible to decide to which extent the number of phases and species can be compromised for\n",
    "efficiency reasons at the expense of chemical realism."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cbada43",
   "metadata": {},
   "source": [
    "### Specifying aqueous phase\n",
    "\n",
    "To initialize [AqueousPhase](https://reaktoro.org/cpp/classReaktoro_1_1AqueousPhase.html#details),\n",
    "we can add the list of exact names of aqueous species we wish to be simulated in computations.\n",
    "\n",
    "> **Note**: Alternative methods to specify species in aqueous phase are described in documentation of the\n",
    "[ChemicalEditor](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalEditor.html) class or tutorial\n",
    "[**ChemicalEditor**](cl.chemical-editor.ipynb).\n",
    "\n",
    "The default model to calculate the activities of solvent water and ionic species is the HKF model.\n",
    "Note that activity models are also needed for *neutral species* then an ideal model is used, in which their activity\n",
    "coefficients are one. For some neutral aqueous species, such as CO<sub>2</sub>(aq), we specify the Drummond model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46e304b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "editor.addAqueousPhase([\"H2O(l)\", \"H+\", \"OH-\", \"Na+\", \"Cl-\", \"HCO3-\", \"CO2(aq)\", \"CO3--\"]) \\\n",
    "    .setActivityModelDrummondCO2()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fa70629",
   "metadata": {},
   "source": [
    "Alternatively, one can select ideal activity model by\n",
    "[AqueousPhase::setActivityModelIdeal](https://reaktoro.org/cpp/classReaktoro_1_1AqueousPhase.html#aef6fbb39539c771554b628de53edeab7)\n",
    "method for neutral aqueous species by providing the corresponding name, e.g.,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "119b2b9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "editor.addAqueousPhase([\"H2O(l)\", \"H+\", \"OH-\", \"Na+\", \"Cl-\", \"HCO3-\", \"CO2(aq)\", \"CO3--\"]) \\\n",
    "    .setActivityModelIdeal(\"CO2(aq)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d38888cf",
   "metadata": {},
   "source": [
    "Finally, we can set the activity model of neutral aqueous species to be the Setschenow one, where not only name of\n",
    "the species but also the Setschenow constant must be provided in method\n",
    "[AqueousPhase::setActivityModelSetschenow](https://reaktoro.org/cpp/classReaktoro_1_1AqueousPhase.html#ab1d6cc44d10a01c7bc1380d6cedfff79)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02625d15",
   "metadata": {},
   "outputs": [],
   "source": [
    "editor.addAqueousPhase([\"H2O(l)\", \"H+\", \"OH-\", \"Na+\", \"Cl-\", \"HCO3-\", \"CO2(aq)\", \"CO3--\", \"NaCl(aq)\"]) \\\n",
    "    .setChemicalModelDebyeHuckel() \\\n",
    "    .setActivityModelSetschenow(\"NaCl(aq)\", 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63269cbc",
   "metadata": {},
   "source": [
    "To choose, for instance, the Pitzer equation of state for aqueous species, we can use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17a8efc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "editor.addAqueousPhase([\"H2O(l)\", \"H+\", \"OH-\", \"Na+\", \"Cl-\", \"HCO3-\", \"CO2(aq)\", \"CO3--\"]) \\\n",
    "    .setChemicalModelPitzerHMW() \\\n",
    "    .setActivityModelDrummondCO2()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cc4a419",
   "metadata": {},
   "source": [
    "The instance of the [AqueousPhase](https://reaktoro.org/cpp/classReaktoro_1_1AqueousPhase.html#details) class can\n",
    "also be constructed explicitly and have its default *chemical model* changed to the Debye-Huckel model, for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de563d71",
   "metadata": {},
   "outputs": [],
   "source": [
    "aqueous_phase = editor.addAqueousPhaseWithElementsOf(\"H2O NaCl CO2\")\n",
    "aqueous_phase.setChemicalModelDebyeHuckel()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b3bbcac",
   "metadata": {},
   "source": [
    "Method\n",
    "[ChemicalEditor::addAqueousPhaseWithElementsOf](\n",
    "https://reaktoro.org/cpp/classReaktoro_1_1ChemicalEditor.html#a23e44f994b87c650a949226ddc195710)\n",
    "used above is practical in those occasions, where it would be preferable to\n",
    "just specify a few compound or substance names, *not necessarily named as in the database*, and then let\n",
    "[ChemicalEditor](https://reaktoro.org/cpp/classReaktorobl_1_1ChemicalEditor.html) to automatically select all chemical\n",
    "species that could be formed out of the combination of those compounds."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00dd447a",
   "metadata": {},
   "source": [
    "### Specifying gaseous phase"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "faa2bd15",
   "metadata": {},
   "source": [
    "Similarly, [GaseousPhase](https://reaktoro.org/cpp/classReaktoro_1_1GaseousPhase.html) can be defined either from\n",
    "list of provided elements or from substance names that are parsed by the\n",
    "[ChemicalEditor](https://reaktoro.org/cpp/classReaktorobl_1_1ChemicalEditor.html) to generate all possible\n",
    "species that can be combined from elements used in those lists. However, for this particular example,\n",
    "the water vapor, H<sub>2</sub>O(g), and gaseous/supercritical carbon dioxide, CO<sub>2</sub>(g), suffices to\n",
    "represent the gaseous phase:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c4afc07",
   "metadata": {},
   "outputs": [],
   "source": [
    "editor.addGaseousPhase([\"H2O(g)\", \"CO2(g)\"]) \\\n",
    "    .setChemicalModelSpycherPruessEnnis()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac443a5c",
   "metadata": {},
   "source": [
    "Here, method [GaseousPhase::setChemicalModelSpycherPruessEnnis](https://reaktoro.org/cpp/classReaktoro_1_1FluidPhase.html#a7fc8a1bc87f24b31939ce3295cd92e8f)\n",
    "sets Spycher et al. (2003) equation of state. This model only\n",
    "supports the gaseous species H<sub>2</sub>O(g) and CO<sub>2</sub>(g). Any other species will result in a runtime\n",
    "error. Alternately, the Spycher and Reed (1988) equation of state can be set (only for H<sub>2</sub>O(g),\n",
    "CO<sub>2</sub>(g), and CH<sub>4</sub>(aq)).\n",
    "If no model is explicitly specified, the Peng-Robinson equation of state is chosen by default to calculate the\n",
    "thermodynamic and chemical properties of this GaseousPhase object."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "822daf9b",
   "metadata": {},
   "source": [
    "### Specifying mineral phase\n",
    "\n",
    "In most geochemical modeling applications, one or more mineral phases are needed. In most cases, these mineral\n",
    "phases are *pure mineral phases*, i.e., they contain only one mineral species. If more than one minerals are\n",
    "present, they are often called *solid solutions*. Defining a pure mineral phase or a solid solution phase is\n",
    "similar to defining any other phase type. The code below demonstrates addition of pure mineral phases halite\n",
    "NaCl:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44ba6354",
   "metadata": {},
   "outputs": [],
   "source": [
    "editor.addMineralPhase(\"Halite\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4b90bed",
   "metadata": {},
   "source": [
    "Definition of solid solution can performed as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3bb2d90",
   "metadata": {},
   "outputs": [],
   "source": [
    "editor.addMineralPhase([\"Calcite\", \"Magnesite\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe19e025",
   "metadata": {},
   "source": [
    "### Creating chemical system\n",
    "\n",
    "Next, we create chemical system using the class [ChemicalSystem](\n",
    "https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html), which helps to specify system's attributes and\n",
    "properties:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a896b105",
   "metadata": {},
   "outputs": [],
   "source": [
    "system = ChemicalSystem(editor)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d94621de",
   "metadata": {},
   "source": [
    "Once system is set, the equilibrium problem using [EquilibriumProblem](\n",
    "https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumProblem.html#details) class must be initialized by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24d51ffb",
   "metadata": {},
   "outputs": [],
   "source": [
    "problem = EquilibriumProblem(system)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11baafd8",
   "metadata": {},
   "source": [
    "We use class [EquilibriumProblem](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumProblem.html) to specify\n",
    "the conditions at which our system should be in equilibrium.\n",
    "For the equilibrium calculation, we have set temperature and pressure with optional units.\n",
    "\n",
    "> **Note**: The default values are 25 &deg;C for the temperature and 1 bar for pressure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4def28c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "problem.setTemperature(60, \"celsius\")\n",
    "problem.setPressure(300, \"bar\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65bf5b32",
   "metadata": {},
   "source": [
    "Additionally, we have to add amount of the solutions used in the equilibrium calculations. For NaCl-brine,\n",
    "we mix 1 kg of water with 1 mol of salt. Plus, we take 100 g of CO<sub>2</sub>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "087a6154",
   "metadata": {},
   "outputs": [],
   "source": [
    "problem.add(\"H2O\", 1, \"kg\")\n",
    "problem.add(\"NaCl\", 1, \"mol\")\n",
    "problem.add(\"CO2\", 100, \"g\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00b93b6a",
   "metadata": {},
   "source": [
    "The units above can be changed, or even suppressed. If not provided, default units are used, such as K for\n",
    "temperatures, Pa for pressures, and mol for amounts. The\n",
    "[EquilibriumProblem::add](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumProblem.html#a18d8ff0e8b8fc66f1eec0127057c7d54)\n",
    "method supports both amount and mass units, such as `mmol`,  `umol`, `g`, `mg`, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79717958",
   "metadata": {},
   "source": [
    "To provide computational representation of the state of a multiphase chemical system resulting from  equilibration\n",
    "process, class [ChemicalState](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html) must be used.\n",
    "Function [equilibrate](https://reaktoro.org/cpp/namespaceReaktoro.html#af2d3b39d3e0b8f9cb5a4d9bbb06b697e)\n",
    "equilibrates a chemical state instance with an equilibrium problem."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e74b102e",
   "metadata": {},
   "source": [
    "The code below uses the definition of the equilibrium problem stored in the object `problem` to perform the\n",
    "equilibrium calculation with utility method [equilibrate](https://reaktoro.org/cpp/namespaceReaktoro.html#af2d3b39d3e0b8f9cb5a4d9bbb06b697e). The result of the calculation is the object `state`,\n",
    "an instance of [ChemicalState](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html) class, which is used\n",
    "to store the chemical state (i.e., the temperature, pressure, and molar amounts of all species) of the system at\n",
    "prescribed equilibrium conditions. The [ChemicalState](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html)\n",
    "class also provides methods for querying thermodynamic properties of the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f9c6abe",
   "metadata": {},
   "outputs": [],
   "source": [
    "state = equilibrate(problem)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "234db9b0",
   "metadata": {},
   "source": [
    "**Note:** Method [equilibrate](https://reaktoro.org/cpp/namespaceReaktoro.html#af2d3b39d3e0b8f9cb5a4d9bbb06b697e)\n",
    "is not the optimal method for performing a sequence of equilibrium calculations\n",
    "(e.g., when coupling Reaktoro with other codes for simulating fluid flow, species transport, etc.). In situations,\n",
    "where many equilibrium calculations need to be performed and sufficient initial guesses are available each time\n",
    "(e.g., the equilibrium state at a previous time step serving as an initial guess for the new equilibrium\n",
    "calculation), use class [EquilibriumSolver](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumSolver.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2667403b",
   "metadata": {},
   "source": [
    "The chemical state can be printed to console by `print(state)` command or saved to a file as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac1ba42c",
   "metadata": {},
   "outputs": [],
   "source": [
    "state.output(\"state.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40f74902",
   "metadata": {},
   "source": [
    "The output information contains details about the equilibrium state of the defined chemical system for the\n",
    "given equilibrium conditions, including, for example, the *amounts*, *masses*, *mole fractions*, *activities*,\n",
    "*activity coefficients*, and *chemical potentials* of the chemical species in each phase. It will also list\n",
    "properties of each phase, such as *density*, *molar volume*, *volume fraction*, as well specific properties of some\n",
    "phases (e.g., *ionic strength*, *pH*, *pe* for the aqueous phase)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2788771a",
   "metadata": {},
   "source": [
    "### Options for the equilibrium calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1294c8e",
   "metadata": {},
   "source": [
    "To customize options for the equilibrium calculations, class\n",
    "[EquilibriumOptions](https://reaktoro.org/cpp/structReaktoro_1_1EquilibriumOptions.html) can be used.\n",
    "For instance, [NonlinearOptions](https://reaktoro.org/cpp/structReaktoro_1_1NonlinearOptions.html) contains\n",
    "information about the nonlinear solver used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a22a394",
   "metadata": {},
   "outputs": [],
   "source": [
    "options = EquilibriumOptions()\n",
    "options.optimum.output.active = True\n",
    "options.epsilon = 1e-50"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7674f07f",
   "metadata": {},
   "source": [
    "Here, we set the field of [OutputterOptions](https://reaktoro.org/cpp/structReaktoro_1_1OutputterOptions.html)\n",
    "class to be `True` to determine whether the intermediate values of equilibrium simulations must be output\n",
    "to the console.\n",
    "Then, we set the parameter *epsilon* (used for the numerical representation of a zero molar amount) to be equal to\n",
    "1e-50. The molar amount of the *i*th species is considered zero if *n[i] < epsilon min b*, where *b* is\n",
    "the vector of element molar amounts."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dbb8a5c",
   "metadata": {},
   "source": [
    "To encounter provided above options, we use function [equilibrate](\n",
    "https://reaktoro.org/cpp/namespaceReaktoro.html#a908245bfa7d236d8c556241dc87f489e)\n",
    "that accepts not only the instance of equilibrium problem but also the specified `options`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8b97023",
   "metadata": {},
   "outputs": [],
   "source": [
    "state = equilibrate(problem, options)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}