{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Creating thermodynamic databases\n", "\n", "
Written by Allan Leal (ETH Zurich) on Jan 14th, 2022
\n", "\n", "```{attention}\n", "Always make sure you are using the [latest version of Reaktoro](https://anaconda.org/conda-forge/reaktoro). Otherwise, some new features documented on this website will not work on your machine and you may receive unintuitive errors. Follow these [update instructions](updating_reaktoro_via_conda) to get the latest version of Reaktoro!\n", "```\n", "\n", "In this tutorial, we will demonstrate how you create your own thermodynamic databases in Reaktoro's YAML format. In these databases, you will have blocks of data for each species (and possibly elements as well, if they are not standard chemical elements from the periodic table).\n", "\n", "It is important to understand that Reaktoro expects each chemical species to know how to calculate its own *standard thermodynamic properties* at a given temperature $T$ and pressure $P$. The following standard thermodynamic properties of a species are expected as output in a standard thermodynamic model in Reaktoro:\n", "\n", "* its standard molar Gibbs energy of formation $G^{\\circ}$ (in J/mol)\n", "* its standard molar enthalpy of formation $H^{\\circ}$ (in J/mol)\n", "* its standard molar volume $V^{\\circ}$ (in m³/mol)\n", "* its standard molar isobaric heat capacity $C_{P}^{\\circ}$ (in J/(mol·K))\n", "* the temperature derivative of $V^{\\circ}$ at constant pressure $(\\partial V^{\\circ}/\\partial T)_P$ of the species (in m³/(mol·K))\n", "* the pressure derivative of $V^{\\circ}$ at constant temperature $(\\partial V^{\\circ}/\\partial P)_T$ (in m³/(mol·Pa))\n", "\n", "After computing these properties for a given $(T,P)$, several others can be easily obtained using thermodynamic relationships (e.g. $U^{\\circ} = H^{\\circ} - PV^{\\circ }$).\n", "\n", "Reaktoro already implements several standard thermodynamic models for you:\n", "\n", "* `StandardThermoModelConstant`\n", "* `StandardThermoModelInterpolation`\n", "* `StandardThermoModelHKF`\n", "* `StandardThermoModelMineralHKF`\n", "* `StandardThermoModelWaterHKF`\n", "* `StandardThermoModelMaierKelley`\n", "* `StandardThermoModelHollandPowell`\n", "* `StandardThermoModelNasa`\n", "\n", "We understand that thermodynamic data are not always formulated and organized in terms of standard thermodynamic properties of species. Many thermodynamic databases exist in which reaction properties are given (e.g., equilibrium constant of reaction, enthalpy of reaction, etc.). Reaktoro also implements several *standard thermodynamic models for reactions* to meet this need:\n", "\n", "* `ReactionStandardThermoModelConstLgK`\n", "* `ReactionStandardThermoModelVantHoff`\n", "* `ReactionStandardThermoModelPhreeqcLgK`\n", "* `ReactionStandardThermoModelGemsLgK`\n", "\n", "These models also expect temperature $T$ and pressure $P$ as inputs, and their output should be:\n", "\n", "* the standard molar Gibbs energy of the reaction $\\Delta G^{\\circ}$ (in J/mol)\n", "* the standard molar enthalpy of the reaction $\\Delta H^{\\circ}$ (in J/mol)\n", "\n", "```{note}\n", "In the background, Reaktoro will convert your *thermodynamic reaction model* assigned to a species into a *standard thermodynamic model* for the same species.\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A very simple custom thermodynamic database in YAML format\n", "\n", "Let's consider fictitious gaseous species A+, B-, and AB. The elements constituting these species are A and B. We consider the following reaction among them:\n", "\n", "$$\\mathrm{A}^{+}+\\mathrm{B}^{-}\\rightleftharpoons\\mathrm{AB}\\qquad(\\lg K_{\\mathrm{AB}}=2)$$\n", "\n", "where $\\lg\\equiv\\log_{10}$.\n", "\n", "```{note}\n", "This example can be modified to work with species of other aggregate state as well, not just gaseous. And yes, gases can exist in charged state: plasma!\n", "```\n", "\n", "Let's say we are interested in calculating the equilibrium state of these species at 25 {{deg}}C and 1 atm starting from an initial state with the same temperature and pressure containing 1 mol of AB. The first step is to create a database file in Reaktoro's YAML format containing the details and thermodynamic data of these species:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "```{tab} Old Format\n", "~~~{literalinclude} db-fictitious-simple.yaml\n", ":language: yaml\n", "~~~\n", "```\n", "```{tab} New Format\n", "~~~{literalinclude} db-fictitious-simple-new-format.yaml\n", ":language: yaml\n", "~~~\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "```{note}\n", "Reaktoro now supports a new experimental database format (click the \"New Format\" tab above). In this new format, YAML dictionaries are used instead of YAML lists to store species and element data. You will notice that each new element and species is now preceded by its unique identifier (the element symbol and the species name). This format has some advantages when working with online graphical editors for YAML and JSON files, with elements and species being found more easily in the tree. Furthermore, this format will also have some advantages in future versions of Reaktoro when the user wants to access any parameter associated with a species and change its value to see how the calculated chemical states (via equilibrium or kinetics computations) are affected.\n", "```\n", "\n", "```{attention}\n", "Reaktoro strictly adopts SI units! Thus, molar mass must be provided in kg/mol, energy in J, mass in kg, volume in m3, amount in mol, temperature in K, pressure in Pa, etc.\n", "```\n", "\n", "Remarks on the above database:\n", "\n", "**Remark 1:** Elements A and B in this problem are not on the periodic table (element B here is not boron!). Thus, they need to be entered into the database. They are needed when building the chemical species. Note that their symbols are referred to in the `Elements` field of the species, eg `Elements: 1:A 1:B` for the `AB` species.\n", "\n", "**Remark 2:** Species `A+` and `B-` compute their standard thermodynamic properties with a *constant standard thermodynamic model* implemented in `StandardThermoModelConstant`. For any $(T,P)$, its properties are constants. By setting `G0` to zero it means $G^{\\circ}_i=0$ for any temperature and pressure. We could have assigned values to `H0`, `V0`, `VT0`, `VP0`, and `Cp0` as well. If not provided, they are zero.\n", "\n", "**Remark 3:** The thermodynamic data available for species `AB` is in the form of reaction properties (reaction equilibrium constant). So we used `FormationReaction` above instead of `StandardThermoModel` for this. Reaktoro requires that a reaction that provides thermodynamic properties for a species be written as a *formation reaction*. In this example, the species `AB` should be considered as a product in the reaction and both `A+` and `B-` as reactants. The reaction is therefore expected to be `A+ + B- = AB`, with `AB` on the right side. In the `Reactants` field of `FormationReaction` for species `AB` we provide `1:A+ 1:B-`.\n", "\n", "**Remark 4:** If the reaction has more than one product, the other products are considered reactants with negative coefficients. For example, `H+ + HCO3- = CO2 + H2O`, a forming reaction for `CO2`, also has `H2O` as a product. This reaction can be interpreted as `H+ + HCO3- - H2O = CO2`, which would produce `Reagents: 1:H+ 1:HCO3- -1:H2O` for the `FormationReaction` of `CO2`. You'll see this in the next section where we show another example of constructing a custom database.\n", "\n", "**Remark 5:** The `AB` species calculates its standard thermodynamic properties with a *constant lgK reaction model* implemented in `ReactionStandardThermoModelConstLgK`. For each $(T,P)$ the equilibrium constant is the same. This equilibrium constant $K_{\\mathrm{AB}}$ is thus converted by the model to $\\Delta G_{\\mathrm{AB}}^{\\circ}$ as follows:\n", "\n", "$$\\Delta G_{\\mathrm{AB}}^{\\circ}=-RT\\ln K_{\\mathrm{AB}}$$\n", "\n", "and from the definition of $\\Delta G_{\\mathrm{AB}}^{\\circ}$:\n", "\n", "$$\\Delta G_{\\mathrm{AB}}^{\\circ}\\equiv G_{\\mathrm{AB}}^{\\circ}-G_{\\mathrm{A^{+}}}^{\\circ} -G_{\\mathrm{B^{-}}}^{\\circ}$$\n", "\n", "Reaktoro will calculate $G_{\\mathrm{AB}}^{\\circ}$ using:\n", "\n", "$$G_{\\mathrm{AB}}^{\\circ}=\\Delta G_{\\mathrm{AB}}^{\\circ}+G_{\\mathrm{A^{+}}}^{\\circ}+ G_{\\mathrm{B^{-}}}^{\\circ}.$$\n", "\n", "This is the principle behind Reaktoro's automated transformation of a thermodynamic model for a reaction to a standard thermodynamic model for a species.\n", "\n", "**Remark 6:** The above conversion technique can produce infinite recursion if all species are defined in terms of formation reactions. We leave the mathematical identification of this recursion for you to figure out. Just note that it would make no sense to define the formation reaction for `AB` in terms of `A+` and `B-` and also define formation reactions for `A+` and `B-` in terms of themselves and `AB` . Infinite recursion can be removed by assigning a standard thermodynamic model to some species, the basic/master species. For example, a constant model that sets $G_{\\mathrm{A^{+}}}^{\\circ}$ and $G_{\\mathrm{B^{-}}}^{\\circ}$ to zero for the basic species.\n", "\n", "```{note}\n", "The above convention of formation reactions is only necessary when providing the means to calculate standard thermodynamic properties for a species. When modeling chemical kinetics, Reaktoro has no absolute imposition on how you define your reactions!\n", "```\n", "\n", "Let us now proceed with our chemical equilibrium calculation. We've saved this database file locally, with file name `db-fictious-simple.yaml`, which is used next:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "INITIAL STATE\n", "+-----------------+------------+------+\n", "| Property | Value | Unit |\n", "+-----------------+------------+------+\n", "| Temperature | 298.1500 | K |\n", "| Pressure | 1.0000 | bar |\n", "| Charge: | 0.0000e+00 | mol |\n", "| Element Amount: | | |\n", "| :: A | 1.0000e+00 | mol |\n", "| :: B | 1.0000e+00 | mol |\n", "| Species Amount: | | |\n", "| :: A+ | 1.0000e-16 | mol |\n", "| :: B- | 1.0000e-16 | mol |\n", "| :: AB | 1.0000e+00 | mol |\n", "+-----------------+------------+------+\n", "FINAL STATE\n", "+-----------------+------------+------+\n", "| Property | Value | Unit |\n", "+-----------------+------------+------+\n", "| Temperature | 298.1500 | K |\n", "| Pressure | 1.0000 | bar |\n", "| Charge: | 0.0000e+00 | mol |\n", "| Element Amount: | | |\n", "| :: A | 1.0000e+00 | mol |\n", "| :: B | 1.0000e+00 | mol |\n", "| Species Amount: | | |\n", "| :: A+ | 9.9504e-02 | mol |\n", "| :: B- | 9.9504e-02 | mol |\n", "| :: AB | 9.0050e-01 | mol |\n", "+-----------------+------------+------+\n" ] } ], "source": [ "from reaktoro import *\n", "\n", "db = Database.fromFile(\"db-fictitious-simple.yaml\")\n", "\n", "gases = GaseousPhase(\"A+ B- AB\")\n", "\n", "system = ChemicalSystem(db, gases)\n", "\n", "state = ChemicalState(system)\n", "state.temperature(25.0, \"celsius\")\n", "state.pressure(1.0, \"bar\")\n", "state.set(\"AB\", 1.0, \"mol\")\n", "\n", "print(\"INITIAL STATE\")\n", "print(state)\n", "\n", "equilibrate(state)\n", "\n", "print(\"FINAL STATE\")\n", "print(state)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A more complicated custom thermodynamic database in YAML format\n", "\n", "We now consider a more complicated database in which thermodynamic reaction properties available in the PHREEQC database `phreeqc.dat` were collected and ported to Reaktoro's YAML format. The database is shown below:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "```{tab} Old Format\n", "~~~{literalinclude} db-geo-complicated.yaml\n", ":language: yaml\n", "~~~\n", "```\n", "```{tab} New Format\n", "~~~{literalinclude} db-geo-complicated-new-format.yaml\n", ":language: yaml\n", "~~~\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will now use this database, saved locally under the name `db-geo-complicated.yaml`, to perform a chemical equilibrium calculation in an initial chemical state at 25 {{deg}}C and 1 bar containing 1 kg of H{{_2}}O, 0.1 mol of CO{{_2}}, and 1 mol of mineral calcite (CaCO{{_3}}):" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "FINAL STATE USING CUSTOM DATABASE\n", "+---------------------+-------------+------+\n", "| Property | Value | Unit |\n", "+---------------------+-------------+------+\n", "| Temperature | 298.1500 | K |\n", "| Pressure | 1.0000 | bar |\n", "| Charge: | -1.0061e-16 | mol |\n", "| Element Amount: | | |\n", "| :: H | 1.1102e+02 | mol |\n", "| :: C | 1.1000e+00 | mol |\n", "| :: O | 5.8708e+01 | mol |\n", "| :: Ca | 1.0000e+00 | mol |\n", "| Species Amount: | | |\n", "| :: H2O | 5.5499e+01 | mol |\n", "| :: H+ | 2.2607e-06 | mol |\n", "| :: OH- | 4.4663e-09 | mol |\n", "| :: HCO3- | 1.7870e-02 | mol |\n", "| :: CO2 | 9.1064e-02 | mol |\n", "| :: CO3-2 | 3.7052e-07 | mol |\n", "| :: Ca+2 | 8.9341e-03 | mol |\n", "| :: Calcite :: CaCO3 | 9.9107e-01 | mol |\n", "+---------------------+-------------+------+\n" ] } ], "source": [ "from reaktoro import *\n", "\n", "db = Database.fromFile(\"db-geo-complicated.yaml\")\n", "\n", "solution = AqueousPhase(\"H2O H+ OH- HCO3- CO2 CO3-2 Ca+2\")\n", "mineral = MineralPhase(\"Calcite\")\n", "\n", "system = ChemicalSystem(db, solution, mineral)\n", "\n", "state = ChemicalState(system)\n", "state.temperature(25.0, \"celsius\")\n", "state.pressure(1.0, \"bar\")\n", "state.set(\"H2O\", 1.0, \"kg\")\n", "state.set(\"CO2\", 0.1, \"mol\")\n", "state.set(\"Calcite\", 1.0, \"mol\")\n", "\n", "equilibrate(state)\n", "\n", "print(\"FINAL STATE USING CUSTOM DATABASE\")\n", "print(state)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's now compare the above calculation with one using the PHREEQC database `phreeqc.dat` directly:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "FINAL STATE USING PHREEQC.DAT\n", "+---------------------+-------------+------+\n", "| Property | Value | Unit |\n", "+---------------------+-------------+------+\n", "| Temperature | 298.1500 | K |\n", "| Pressure | 1.0000 | bar |\n", "| Charge: | -9.7145e-17 | mol |\n", "| Element Amount: | | |\n", "| :: H | 1.1101e+02 | mol |\n", "| :: C | 1.1000e+00 | mol |\n", "| :: O | 5.8706e+01 | mol |\n", "| :: Ca | 1.0000e+00 | mol |\n", "| Species Amount: | | |\n", "| :: H2O | 5.5497e+01 | mol |\n", "| :: H+ | 2.2610e-06 | mol |\n", "| :: OH- | 4.4652e-09 | mol |\n", "| :: HCO3- | 1.7871e-02 | mol |\n", "| :: CO2 | 9.1063e-02 | mol |\n", "| :: CO3-2 | 3.7059e-07 | mol |\n", "| :: Ca+2 | 8.9348e-03 | mol |\n", "| :: Calcite :: CaCO3 | 9.9107e-01 | mol |\n", "+---------------------+-------------+------+\n" ] } ], "source": [ "from reaktoro import *\n", "\n", "db = PhreeqcDatabase(\"phreeqc.dat\")\n", "\n", "solution = AqueousPhase(\"H2O H+ OH- HCO3- CO2 CO3-2 Ca+2\")\n", "mineral = MineralPhase(\"Calcite\")\n", "\n", "system = ChemicalSystem(db, solution, mineral)\n", "\n", "state = ChemicalState(system)\n", "state.temperature(25.0, \"celsius\")\n", "state.pressure(1.0, \"bar\")\n", "state.set(\"H2O\", 1.0, \"kg\")\n", "state.set(\"CO2\", 0.1, \"mol\")\n", "state.set(\"Calcite\", 1.0, \"mol\")\n", "\n", "equilibrate(state)\n", "\n", "print(\"FINAL STATE USING PHREEQC.DAT\")\n", "print(state)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The amounts of species above should be comparable. Not identical, because for some species in our custom database, we are not using the same calculation mode for equilibrium constants (e.g. in our custom database we use a certain value for $\\lg K$ for `CO2`, 16.681, while in `phreeqc.dat` an analytic expression is used)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A very simple custom thermodynamic database using Python API\n", "\n", "We will now construct a {{Database}} object using Reaktoro's Python API only (no YAML files!). This is a more complicated way to construct a custom {{Database}} object, but maybe this guide can be helpful to you in ways we cannot foresee now. \n", "\n", "We want to construct a database to model the equilibrium of `H2O`, `H+`, and `OH-`. The elements composing these species are `H` and `O`. \n", "\n", "The code block below implements function `createElementList` that construct an {{ElementList}} object containing the {{Element}} objects that constitute the {{Species}} objects we will construct next:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def createElementList():\n", " elements = ElementList()\n", "\n", " elements.append(Element()\n", " .withName(\"Hydrogen\")\n", " .withSymbol(\"H\")\n", " .withMolarMass(1.0e-3))\n", "\n", " elements.append(Element()\n", " .withName(\"Oxygen\")\n", " .withSymbol(\"O\")\n", " .withMolarMass(16.0e-3))\n", "\n", " return elements" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now define a function `createSpeciesList` that construct the {{Species}} objects that will populate our {{Database}} object:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "\n", "def createSpeciesList(elements: ElementList):\n", " \n", " # Create the primary species (used to define the secondary species)\n", "\n", " primaryspecies = SpeciesList()\n", "\n", " primaryspecies.append(Species()\n", " .withName(\"H+\")\n", " .withFormula(\"H+\")\n", " .withCharge(1.0)\n", " .withElements(ElementalComposition([ (elements.get(\"H\"), 1.0) ]))\n", " .withAggregateState(AggregateState.Aqueous)\n", " .withStandardGibbsEnergy(0.0))\n", "\n", " primaryspecies.append(Species()\n", " .withName(\"OH-\")\n", " .withFormula(\"OH-\")\n", " .withCharge(-1.0)\n", " .withElements(ElementalComposition([ (elements.get(\"O\"), 1.0), (elements.get(\"H\"), 1.0) ]))\n", " .withAggregateState(AggregateState.Aqueous)\n", " .withStandardGibbsEnergy(0.0))\n", "\n", " # Create the secondary species (using formation reactions from primary species)\n", "\n", " secondaryspecies = []\n", "\n", " secondaryspecies.append(Species()\n", " .withName(\"H2O\")\n", " .withFormula(\"H2O\")\n", " .withElements(ElementalComposition([ (elements.get(\"H\"), 2.0), (elements.get(\"O\"), 1.0) ]))\n", " .withAggregateState(AggregateState.Aqueous)\n", " .withFormationReaction(\n", " FormationReaction()\n", " .withReactants([\n", " (primaryspecies.get(\"H+\") , 1.0),\n", " (primaryspecies.get(\"OH-\"), 1.0) ])\n", " .withEquilibriumConstant(14.0)))\n", "\n", " return primaryspecies + secondaryspecies # concatenate both species lists" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now create our auxiliary `createCustomDatabase` function that creates our {{Database}} object:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def createCustomDatabase():\n", " elements = createElementList()\n", " species = createSpeciesList(elements)\n", " return Database(species.data())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are now ready for our chemical equilibrium calculation:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "+-----------------+------------+------+\n", "| Property | Value | Unit |\n", "+-----------------+------------+------+\n", "| Temperature | 298.1500 | K |\n", "| Pressure | 1.0000 | bar |\n", "| Charge: | 0.0000e+00 | mol |\n", "| Element Amount: | | |\n", "| :: H | 1.1111e+02 | mol |\n", "| :: O | 5.5556e+01 | mol |\n", "| Species Amount: | | |\n", "| :: H+ | 1.0008e-07 | mol |\n", "| :: OH- | 1.0008e-07 | mol |\n", "| :: H2O | 5.5556e+01 | mol |\n", "+-----------------+------------+------+\n" ] } ], "source": [ "db = createCustomDatabase()\n", "\n", "solution = AqueousPhase()\n", "\n", "system = ChemicalSystem(db, solution)\n", "\n", "state = ChemicalState(system)\n", "state.set(\"H2O\", 1.0, \"kg\")\n", "\n", "equilibrate(state)\n", "\n", "print(state)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That's it! We managed to create a Reaktoro script entirely from scratch. From database creation to an equilibrium calculation.\n", "\n", "\n" ] } ], "metadata": { "interpreter": { "hash": "e4e8b2f3ae27709963f14fd23a6560d362beea55eaec742263828e04d814e23c" }, "kernelspec": { "display_name": "Python 3.9.7 64-bit ('reaktoro-jupyter-book': conda)", "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.9.9" }, "orig_nbformat": 4 }, "nbformat": 4, "nbformat_minor": 2 }