{ "cells": [ { "cell_type": "markdown", "id": "da94eeff", "metadata": {}, "source": [ "# Functionality of ChemicalSystem class\n", "\n", "In this tutorial, we provide an explanation on the functionality of the class\n", "[ChemicalSystem](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html)\n", "that represents a chemical system and its attributes and properties. Below,\n", "we provide a tutorial of the methods that can be used to study the characteristics of the considered chemical system.\n", "\n", "A chemical system is a description of the phases of interest in the modeling problem and the chemical species that\n", "compose those phases. For example, when modeling the chemistry of an aqueous solution, only one phase should be\n", "enough, an *aqueous phase*. If one is interested in modeling the solubility of gases in an aqueous solution,\n", "then it makes sense to also define a *gaseous phase* with one or more gases. When modeling aqueous solutions and\n", "minerals, under a variety of temperature, pressure, and chemical conditions, it might make sense to define the\n", "chemical system with many *mineral phases*.\n", "Assume that we have defined `system`, an instance of\n", "[ChemicalSystem](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html) object, by the code below:" ] }, { "cell_type": "code", "execution_count": null, "id": "371269be", "metadata": {}, "outputs": [], "source": [ "# Import the reaktoro Python package\n", "from reaktoro import *\n", "\n", "# Initialize a thermodynamic database\n", "db = Database(\"supcrt98.xml\")\n", "\n", "# Define the chemical system\n", "editor = ChemicalEditor(db)\n", "editor.addAqueousPhaseWithElements(\n", " \"H O Na Cl C\"\n", ") # add aqueous phase by the all possible combination of provided elements\n", "editor.addGaseousPhase([\"CO2(g)\"]) # add one gaseous species\n", "\n", "# Construct the chemical system\n", "system = ChemicalSystem(editor)" ] }, { "cell_type": "markdown", "id": "17c55300", "metadata": {}, "source": [ "> Check previous tutorials to learn the steps above! For example,\n", "> [Basics of equilibrium calculation](eq.equilibrium-basics.ipynb),\n", "> [Equilibrium calculation of carbonate species](eq.equilibrium-carbonates.ipynb),\n", "> [Equilibrium calculations using equilibrium solver](eq.co2-brine-using-equilibrium-solver.ipynb), or\n", "> [Custom activity model for equilibrium calculations](eq.custom-activity-models.ipynb).\n", "\n", "The most general print-out of the chemical system can be done by" ] }, { "cell_type": "code", "execution_count": null, "id": "7a26f7c0", "metadata": {}, "outputs": [], "source": [ "print(system)" ] }, { "cell_type": "markdown", "id": "12a05f72", "metadata": {}, "source": [ "However, to access some limited information, such as the list of species, phases, or elements separately, the following\n", "code can be used" ] }, { "cell_type": "code", "execution_count": null, "id": "8f06a336", "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "code", "execution_count": null, "id": "113aea10", "metadata": {}, "outputs": [], "source": [ "print(f\"List with {system.numSpecies()} species:\")\n", "for species in system.species():\n", " print(species.name())" ] }, { "cell_type": "code", "execution_count": null, "id": "ff1db8d8", "metadata": {}, "outputs": [], "source": [ "print(f\"List with {system.numPhases()} phases:\")\n", "for phase in system.phases():\n", " print(phase.name())" ] }, { "cell_type": "code", "execution_count": null, "id": "5ffbfd9b", "metadata": {}, "outputs": [], "source": [ "print(f\"List with {(system.numElements())} elements:\")\n", "for element in system.elements():\n", " print(element.name())" ] }, { "cell_type": "markdown", "id": "bfa65d4e", "metadata": {}, "source": [ "To output additional information about phases, for instance, one can use" ] }, { "cell_type": "code", "execution_count": null, "id": "f7a06e97", "metadata": {}, "outputs": [], "source": [ "print(\"List of phases with number of species in each phase:\")\n", "for phase in system.phases():\n", " print(f\" * Phase {phase.name()} contains {phase.numSpecies()} species\")" ] }, { "cell_type": "markdown", "id": "86e8793b", "metadata": {}, "source": [ "[ChemicalSystem](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html)\n", "provides the formula matrix (whose entry *(j, i)* is given by the number of atoms of its *j*th element in\n", "its *i*th species). To access it, one need to use" ] }, { "cell_type": "code", "execution_count": null, "id": "69c5be43", "metadata": {}, "outputs": [], "source": [ "matrix = system.formulaMatrix()\n", "print(f\"Formula matrix of the size {matrix.shape[0]} x {matrix.shape[1]}:\\n\", matrix)" ] }, { "cell_type": "markdown", "id": "a6869cf0", "metadata": {}, "source": [ "For instance, the matrix printed above is the matrix of the size *6 x 23*, i.e.,\n", "\n", "\\begin{bmatrix}\n", "1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\\\\n", "0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 2 & 2 & 2 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 \\\\\n", "1 & 2 & 3 & 0 & 1 & 2 & 3 & 4 & 0 & 0 & 1 & 2 & 3 & 0 & 1 & 2 & 2 & 0 & 0 & 1 & 2 & 1 & 2 \\\\\n", "0 & 0 & -2 & -1 & -1 & -1 & -1 & -1 & 1 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & -1 & 1 & 0 & 0 & 0 & -1 & 0 \\\\\n", "\\end{bmatrix}.\n", "\n", "Assume the enumeration starts from 0 (which is the case in Python language). Here, the row with index 1 corresponds to\n", "the element Cl, which is only present in species\n", "Cl-,\n", "ClO-,\n", "ClO2-,\n", "ClO3-,\n", "ClO4-\n", "(with indices 3, 4, 5, 6, and 7),\n", "HCl(aq),\n", "HClO(aq),\n", "HClO_2(aq)\n", "(with indices 13, 14, and 15),\n", "and\n", "NaCl(aq) (with index 18).\n", "Since only one atom of Cl contributes to each species, the second row contains only 1 in non-zero values.\n", "\n", "To get the index of certain element, phase, or species, functions `index__()` or `index__WithError()` (the latter\n", "results in system throwing an exception if the element does not exist), where instead of `__` one can used `Element`,\n", "`Phase`, or `Spieces`:" ] }, { "cell_type": "code", "execution_count": null, "id": "f82b3bd4", "metadata": {}, "outputs": [], "source": [ "print(\"Index of the element H: \", system.indexElement(\"H\"))\n", "print(\"Index of the phase Aqueous: \", system.indexPhase(\"Aqueous\"))\n", "print(\"Index of the species Cl-: \", system.indexSpecies(\"Cl-\"))" ] }, { "cell_type": "markdown", "id": "175c10e1", "metadata": {}, "source": [ "When working with a set of species, one can request a set of corresponding indices. Let us collect all the species\n", "with chlorine and retrieve indices of corresponding species:" ] }, { "cell_type": "code", "execution_count": null, "id": "7517baa3", "metadata": {}, "outputs": [], "source": [ "species = [\n", " \"Cl-\",\n", " \"ClO-\",\n", " \"ClO2-\",\n", " \"ClO3-\",\n", " \"ClO4-\",\n", " \"HCl(aq)\",\n", " \"HClO(aq)\",\n", " \"HClO2(aq)\",\n", " \"NaCl(aq)\",\n", "]\n", "print(\"Indices of species with Cl: \", system.indicesSpecies(species))" ] }, { "cell_type": "markdown", "id": "5f9fad16", "metadata": {}, "source": [ "They must correspond with positions of the non-zero elements in the row 1 of formula matrix discussed-above.\n", "Having the instance of chemical system, we can calculate the molar amounts of the elements (in units of mol), when the\n", "molar amount of species is provided:" ] }, { "cell_type": "code", "execution_count": null, "id": "eb875529", "metadata": {}, "outputs": [], "source": [ "n = np.ones(system.numSpecies())\n", "elements_amount = system.elementAmounts(n)\n", "hydrogen_amount = system.elementAmount(2, n)\n", "print(\n", " \"Element amounts (in mol) provided 1 molal for all species: \",\n", " elements_amount,\n", ")\n", "print(\n", " \"Hydrogen amounts (in mol) provided 1 molal for all species: \",\n", " hydrogen_amount,\n", ")" ] }, { "cell_type": "markdown", "id": "a221c74c", "metadata": {}, "source": [ "To study the chemical system even further, one can access the class\n", "[ThermoProperties](https://reaktoro.org/cpp/classReaktoro_1_1ThermoProperties.html) by providing temperature and\n", "pressure, i.e.," ] }, { "cell_type": "code", "execution_count": null, "id": "97ef7cba", "metadata": {}, "outputs": [], "source": [ "T = 60\n", "P = 100\n", "thermo_properties = system.properties(T, P)" ] }, { "cell_type": "markdown", "id": "ffd5e915", "metadata": {}, "source": [ "Object `thermo_properties` contains information about\n", "* standard partial molar Gibbs energies of the species (in units of J/mol),\n", "* standard partial molar enthalpies of the species (in units of J/mol),\n", "* standard partial molar volumes of the species (in units of m3/mol),\n", "* standard partial molar entropies of the species (in units of J/(mol*K)),\n", "* standard partial molar internal energies of the species (in units of J/mol),\n", "* standard partial molar Helmholtz energies of the species (in units of J/mol),\n", "* standard partial molar isobaric heat capacities of the species (in units of J/(mol*K)),\n", "* standard partial molar isochoric heat capacities of the species (in units of J/(mol*K)).\n", "\n", "For instance, partial molar Gibbs energies or enthalpies can be accessed as follows:" ] }, { "cell_type": "code", "execution_count": null, "id": "38370e8e", "metadata": {}, "outputs": [], "source": [ "print(\"List of standard partial molar Gibbs energies of the species:\")\n", "for energies, species in zip(\n", " thermo_properties.standardPartialMolarGibbsEnergies().val,\n", " system.species()\n", "):\n", " print(f\"\\u03B4G ({species.name():>10}) = {energies}\")" ] }, { "cell_type": "code", "execution_count": null, "id": "9cb4e113", "metadata": {}, "outputs": [], "source": [ "print(\"List of standard partial molar enthalpies of the species:\")\n", "for enthalpies, species in zip(\n", " thermo_properties.standardPartialMolarEnthalpies().val,\n", " system.species()\n", "):\n", " print(f\"\\u03B4H ({species.name():>10}) = {enthalpies}\")" ] }, { "cell_type": "markdown", "id": "d0fe6523", "metadata": {}, "source": [ "Alternatively, by providing the vector with molar amounts of the species (in units of mol) class\n", "[ChemicalProperties](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalProperties.html) can be accessed:" ] }, { "cell_type": "code", "execution_count": null, "id": "22c06661", "metadata": {}, "outputs": [], "source": [ "chemical_properties = system.properties(T, P, n)" ] }, { "cell_type": "markdown", "id": "92001e17", "metadata": {}, "source": [ "This object contains various chemical properties, such as mole fractions, the logarithm of activities, chemical\n", "potentials of the species, in addition to thermodynamic properties listed above:" ] }, { "cell_type": "code", "execution_count": null, "id": "1f27c980", "metadata": {}, "outputs": [], "source": [ "print(\"Chemical potentials of the species:\")\n", "for potential, species, index in zip(\n", " chemical_properties.chemicalPotentials().val,\n", " system.species(),\n", " list(range(1, system.numSpecies()+1))\n", "):\n", " print(f\"\\u03BC_{index} ({species.name():>10}) = {potential}\")" ] }, { "cell_type": "code", "execution_count": null, "id": "98b68474", "metadata": { "lines_to_next_cell": 2 }, "outputs": [], "source": [ "print(\"Logarithms of activities of the species:\")\n", "for activity, species, index in zip(\n", " chemical_properties.lnActivities().val,\n", " system.species(),\n", " list(range(1, system.numSpecies()+1))\n", "):\n", " print(f\"ln (a_{index}) ({species.name():>10}) = {activity}\")" ] }, { "cell_type": "markdown", "id": "0c672ac6", "metadata": {}, "source": [ "### Definition of chemical system using GEMs\n", "\n", "Chemical systems can be initialize using alternative backends, such as GEMs and PHREEQC, using corresponding classes\n", "[Gems](https://reaktoro.org/cpp/classReaktoro_1_1Gems.html) and\n", "[Phreeqc](https://reaktoro.org/cpp/classReaktoro_1_1Phreeqc.html)\n", "`gems = Gems(\"some-gems-project-file.lst\")` or `phreeqc = Phreeqc(\"phreeqc.dat\")`,\n", "which allows to define instance of [ChemicalSystem](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html)\n", "class by respective `system = ChemicalSystem(gems)` or `system = ChemicalSystem(phreeqc)`." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" } }, "nbformat": 4, "nbformat_minor": 5 }