Creating thermodynamic databases#

Written by Allan Leal (ETH Zurich) on Jan 14th, 2022

Attention

Always make sure you are using the latest version of 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 to get the latest version of Reaktoro!

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).

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:

its standard molar Gibbs energy of formation \(G^{\circ}\) (in J/mol)
its standard molar enthalpy of formation \(H^{\circ}\) (in J/mol)
its standard molar volume \(V^{\circ}\) (in m³/mol)
its standard molar isobaric heat capacity \(C_{P}^{\circ}\) (in J/(mol·K))
the temperature derivative of \(V^{\circ}\) at constant pressure \((\partial V^{\circ}/\partial T)_P\) of the species (in m³/(mol·K))
the pressure derivative of \(V^{\circ}\) at constant temperature \((\partial V^{\circ}/\partial P)_T\) (in m³/(mol·Pa))

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 }\)).

Reaktoro already implements several standard thermodynamic models for you:

StandardThermoModelConstant
StandardThermoModelInterpolation
StandardThermoModelHKF
StandardThermoModelMineralHKF
StandardThermoModelWaterHKF
StandardThermoModelMaierKelley
StandardThermoModelHollandPowell
StandardThermoModelNasa

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:

ReactionStandardThermoModelConstLgK
ReactionStandardThermoModelVantHoff
ReactionStandardThermoModelPhreeqcLgK
ReactionStandardThermoModelGemsLgK

These models also expect temperature \(T\) and pressure \(P\) as inputs, and their output should be:

the standard molar Gibbs energy of the reaction \(\Delta G^{\circ}\) (in J/mol)
the standard molar enthalpy of the reaction \(\Delta H^{\circ}\) (in J/mol)

Note

In the background, Reaktoro will convert your thermodynamic reaction model assigned to a species into a standard thermodynamic model for the same species.

A very simple custom thermodynamic database in YAML format#

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:

\[\mathrm{A}^{+}+\mathrm{B}^{-}\rightleftharpoons\mathrm{AB}\qquad(\lg K_{\mathrm{AB}}=2)\]

where \(\lg\equiv\log_{10}\).

Note

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!

Let’s say we are interested in calculating the equilibrium state of these species at 25 °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:

Old Format

Elements:
- Symbol: A
  MolarMass: 0.001  # in kg/mol
- Symbol: B
  MolarMass: 0.002  # in kg/mol

Species:
- Name: A+
  Formula: A+
  Elements: 1:A
  Charge: 1.0
  AggregateState: Gas
  StandardThermoModel:
    Constant:
      G0: 0.0  # in J/mol
- Name: B-
  Formula: B-
  Elements: 1:B
  Charge: -1.0
  AggregateState: Gas
  StandardThermoModel:
    Constant:
      G0: 0.0  # in J/mol
- Name: AB
  Formula: AB
  Elements: 1:A 1:B
  AggregateState: Gas
  FormationReaction:
    Reactants: 1:A+ 1:B-
    ReactionStandardThermoModel:
      ConstLgK:
        lgKr: 2.0

New Format

Elements:
  A:
    Symbol: A
    MolarMass: 0.001  # in kg/mol
  B:
    Symbol: B
    MolarMass: 0.002  # in kg/mol

Species:
  A+:
    Name: A+
    Formula: A+
    Elements: 1:A
    Charge: 1.0
    AggregateState: Gas
    StandardThermoModel:
      Constant:
        G0: 0.0  # in J/mol
  B-:
    Name: B-
    Formula: B-
    Elements: 1:B
    Charge: -1.0
    AggregateState: Gas
    StandardThermoModel:
      Constant:
        G0: 0.0  # in J/mol
  AB:
    Name: AB
    Formula: AB
    Elements: 1:A 1:B
    AggregateState: Gas
    FormationReaction:
      Reactants: 1:A+ 1:B-
      ReactionStandardThermoModel:
        ConstLgK:
          lgKr: 2.0

Note

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.

Attention

Reaktoro strictly adopts SI units! Thus, molar mass must be provided in kg/mol, energy in J, mass in kg, volume in m³, amount in mol, temperature in K, pressure in Pa, etc.

Remarks on the above database:

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.

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.

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-.

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.

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:

\[\Delta G_{\mathrm{AB}}^{\circ}=-RT\ln K_{\mathrm{AB}}\]

and from the definition of \(\Delta G_{\mathrm{AB}}^{\circ}\):

\[\Delta G_{\mathrm{AB}}^{\circ}\equiv G_{\mathrm{AB}}^{\circ}-G_{\mathrm{A^{+}}}^{\circ} -G_{\mathrm{B^{-}}}^{\circ}\]

Reaktoro will calculate \(G_{\mathrm{AB}}^{\circ}\) using:

\[G_{\mathrm{AB}}^{\circ}=\Delta G_{\mathrm{AB}}^{\circ}+G_{\mathrm{A^{+}}}^{\circ}+ G_{\mathrm{B^{-}}}^{\circ}.\]

This is the principle behind Reaktoro’s automated transformation of a thermodynamic model for a reaction to a standard thermodynamic model for a species.

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.

Note

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!

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:

from reaktoro import *

db = Database.fromFile("db-fictitious-simple.yaml")

gases = GaseousPhase("A+ B- AB")

system = ChemicalSystem(db, gases)

state = ChemicalState(system)
state.temperature(25.0, "celsius")
state.pressure(1.0, "bar")
state.set("AB", 1.0, "mol")

print("INITIAL STATE")
print(state)

equilibrate(state)

print("FINAL STATE")
print(state)

INITIAL STATE
+-----------------+------------+------+
| Property        |      Value | Unit |
+-----------------+------------+------+
| Temperature     |   298.1500 |    K |
| Pressure        |     1.0000 |  bar |
| Charge:         | 0.0000e+00 |  mol |
| Element Amount: |            |      |
| :: A            | 1.0000e+00 |  mol |
| :: B            | 1.0000e+00 |  mol |
| Species Amount: |            |      |
| :: A+           | 1.0000e-16 |  mol |
| :: B-           | 1.0000e-16 |  mol |
| :: AB           | 1.0000e+00 |  mol |
+-----------------+------------+------+
FINAL STATE
+-----------------+------------+------+
| Property        |      Value | Unit |
+-----------------+------------+------+
| Temperature     |   298.1500 |    K |
| Pressure        |     1.0000 |  bar |
| Charge:         | 0.0000e+00 |  mol |
| Element Amount: |            |      |
| :: A            | 1.0000e+00 |  mol |
| :: B            | 1.0000e+00 |  mol |
| Species Amount: |            |      |
| :: A+           | 9.9504e-02 |  mol |
| :: B-           | 9.9504e-02 |  mol |
| :: AB           | 9.0050e-01 |  mol |
+-----------------+------------+------+

A more complicated custom thermodynamic database in YAML format#

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:

Old Format

Species:
- Name: H+
  Formula: H+
  Elements: 1:H
  Charge: 1.0
  AggregateState: Aqueous
  StandardThermoModel:
    Constant:
      G0: 0.0  # in J/mol
- Name: OH-
  Formula: OH-
  Elements: 1:O 1:H
  Charge: -1.0
  AggregateState: Aqueous
  StandardThermoModel:
    Constant:
      G0: 0.0  # in J/mol
- Name: CO3-2
  Formula: CO3-2
  Elements: 1:C 3:O
  Charge: -2.0
  AggregateState: Aqueous
  StandardThermoModel:
    Constant:
      G0: 0.0
- Name: Ca+2
  Formula: Ca+2
  Elements: 1:Ca
  Charge: 2.0
  AggregateState: Aqueous
  StandardThermoModel:
    Interpolation:
      Temperatures: [298.15, 350.0, 400.0]  # in K
      Pressures: [1.0e+5, 10.0e+5]          # in Pa
      G0: [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
      H0: [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
- Name: H2O
  Formula: H2O
  Elements: 2:H 1:O
  AggregateState: Aqueous
  FormationReaction:
    Reactants: 1:H+ 1:OH-
    ReactionStandardThermoModel:
      PhreeqcLgK:
        A1: -293.29227
        A2: -0.1360833
        A3: 10576.913
        A4: 123.73158
        A5: 0.0
        A6: 6.996455e-5
- Name: HCO3-
  Formula: HCO3-
  Elements: 1:H 1:C 3:O
  Charge: -1.0
  AggregateState: Aqueous
  FormationReaction:
    Reactants: 1:CO3-2 1:H+
    ReactionStandardThermoModel:
      VantHoff:
        lgKr: 10.329
        dHr: -14899.224  # in J/mol
- Name: CO2
  Formula: CO2
  Elements: 1:C 2:O
  AggregateState: Aqueous
  FormationReaction:
    Reactants: 1:CO3-2 2:H+ -1:H2O
    ReactionStandardThermoModel:
      ConstLgK:
        lgKr: 16.681
- Name: Calcite
  Formula: CaCO3
  Elements: 1:Ca 1:C 3:O
  AggregateState: Solid
  FormationReaction:
    Reactants: 1:CO3-2 1:Ca+2
    ReactionStandardThermoModel:
      VantHoff:
        lgKr: 8.48
        dHr: 10832.376  # in J/mol

New Format

Species:
  H+:
    Name: H+
    Formula: H+
    Elements: 1:H
    Charge: 1.0
    AggregateState: Aqueous
    StandardThermoModel:
      Constant:
        G0: 0.0  # in J/mol
  OH-:
    Name: OH-
    Formula: OH-
    Elements: 1:O 1:H
    Charge: -1.0
    AggregateState: Aqueous
    StandardThermoModel:
      Constant:
        G0: 0.0  # in J/mol
  CO3-2:
    Name: CO3-2
    Formula: CO3-2
    Elements: 1:C 3:O
    Charge: -2.0
    AggregateState: Aqueous
    StandardThermoModel:
      Constant:
        G0: 0.0
  H2O:
    Name: H2O
    Formula: H2O
    Elements: 2:H 1:O
    AggregateState: Aqueous
    FormationReaction:
      Reactants: 1:H+ 1:OH-
      ReactionStandardThermoModel:
        ConstLgK:
          lgKr: 14.0
  HCO3-:
    Name: HCO3-
    Formula: HCO3-
    Elements: 1:H 1:C 3:O
    Charge: -1.0
    AggregateState: Aqueous
    FormationReaction:
      Reactants: 1:CO3-2 1:H+
      ReactionStandardThermoModel:
        VantHoff:
          lgKr: 10.329
          dHr: -14899.224  # in J/mol
  CO2:
    Name: CO2
    Formula: CO2
    Elements: 1:H 1:C 3:O
    AggregateState: Aqueous
    FormationReaction:
      Reactants: 1:CO3-2 2:H+ -1:H2O
      ReactionStandardThermoModel:
        PhreeqcLgK:
          A1: 464.1965
          A2: 0.09344813
          A3: -26986.16
          A4: -165.75951
          A5: 2248628.9
          A6: 0.0
  Ca+2:
    Name: Ca+2
    Formula: Ca+2
    Elements: 1:Ca
    Charge: 2.0
    AggregateState: Aqueous
    StandardThermoModel:
      Interpolation:
        Temperatures: [298.15, 350.0, 400.0]  # in K
        Pressures: [1.0e+5, 10.0e+5]          # in Pa
        G0: [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
        H0: [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
  Calcite:
    Name: Calcite
    Formula: CaCO3
    Elements: 1:Ca 1:C 3:O
    AggregateState: Solid
    FormationReaction:
      Reactants: 1:CO3-2 1:Ca+2
      ReactionStandardThermoModel:
        VantHoff:
          lgKr: 8.48
          dHr: 10832.376  # in J/mol

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 °C and 1 bar containing 1 kg of H₂O, 0.1 mol of CO₂, and 1 mol of mineral calcite (CaCO₃):

from reaktoro import *

db = Database.fromFile("db-geo-complicated.yaml")

solution = AqueousPhase("H2O H+ OH- HCO3- CO2 CO3-2 Ca+2")
mineral = MineralPhase("Calcite")

system = ChemicalSystem(db, solution, mineral)

state = ChemicalState(system)
state.temperature(25.0, "celsius")
state.pressure(1.0, "bar")
state.set("H2O", 1.0, "kg")
state.set("CO2", 0.1, "mol")
state.set("Calcite", 1.0, "mol")

equilibrate(state)

print("FINAL STATE USING CUSTOM DATABASE")
print(state)

FINAL STATE USING CUSTOM DATABASE
+---------------------+-------------+------+
| Property            |       Value | Unit |
+---------------------+-------------+------+
| Temperature         |    298.1500 |    K |
| Pressure            |      1.0000 |  bar |
| Charge:             | -1.0408e-16 |  mol |
| Element Amount:     |             |      |
| :: H                |  1.1102e+02 |  mol |
| :: C                |  1.1000e+00 |  mol |
| :: O                |  5.8708e+01 |  mol |
| :: Ca               |  1.0000e+00 |  mol |
| Species Amount:     |             |      |
| :: H2O              |  5.5499e+01 |  mol |
| :: H+               |  2.2607e-06 |  mol |
| :: OH-              |  4.4663e-09 |  mol |
| :: HCO3-            |  1.7870e-02 |  mol |
| :: CO2              |  9.1064e-02 |  mol |
| :: CO3-2            |  3.7052e-07 |  mol |
| :: Ca+2             |  8.9341e-03 |  mol |
| :: Calcite :: CaCO3 |  9.9107e-01 |  mol |
+---------------------+-------------+------+

Let’s now compare the above calculation with one using the PHREEQC database phreeqc.dat directly:

from reaktoro import *

db = PhreeqcDatabase("phreeqc.dat")

solution = AqueousPhase("H2O H+ OH- HCO3- CO2 CO3-2 Ca+2")
mineral = MineralPhase("Calcite")

system = ChemicalSystem(db, solution, mineral)

state = ChemicalState(system)
state.temperature(25.0, "celsius")
state.pressure(1.0, "bar")
state.set("H2O", 1.0, "kg")
state.set("CO2", 0.1, "mol")
state.set("Calcite", 1.0, "mol")

equilibrate(state)

print("FINAL STATE USING PHREEQC.DAT")
print(state)

FINAL STATE USING PHREEQC.DAT
+---------------------+-------------+------+
| Property            |       Value | Unit |
+---------------------+-------------+------+
| Temperature         |    298.1500 |    K |
| Pressure            |      1.0000 |  bar |
| Charge:             | -1.0061e-16 |  mol |
| Element Amount:     |             |      |
| :: H                |  1.1101e+02 |  mol |
| :: C                |  1.1000e+00 |  mol |
| :: O                |  5.8706e+01 |  mol |
| :: Ca               |  1.0000e+00 |  mol |
| Species Amount:     |             |      |
| :: H2O              |  5.5497e+01 |  mol |
| :: H+               |  2.2613e-06 |  mol |
| :: OH-              |  4.4650e-09 |  mol |
| :: HCO3-            |  1.7870e-02 |  mol |
| :: CO2              |  9.1064e-02 |  mol |
| :: CO3-2            |  3.7053e-07 |  mol |
| :: Ca+2             |  8.9342e-03 |  mol |
| :: Calcite :: CaCO3 |  9.9107e-01 |  mol |
+---------------------+-------------+------+

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).

A very simple custom thermodynamic database using Python API#

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.

We want to construct a database to model the equilibrium of H2O, H+, and OH-. The elements composing these species are H and O.

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:

def createElementList():
    elements = ElementList()

    elements.append(Element()
        .withName("Hydrogen")
        .withSymbol("H")
        .withMolarMass(1.0e-3))

    elements.append(Element()
        .withName("Oxygen")
        .withSymbol("O")
        .withMolarMass(16.0e-3))

    return elements

We now define a function createSpeciesList that construct the Species objects that will populate our Database object:

def createSpeciesList(elements: ElementList):
    
    # Create the primary species (used to define the secondary species)

    primaryspecies = SpeciesList()

    primaryspecies.append(Species()
        .withName("H+")
        .withFormula("H+")
        .withCharge(1.0)
        .withElements(ElementalComposition([ (elements.get("H"),  1.0) ]))
        .withAggregateState(AggregateState.Aqueous)
        .withStandardGibbsEnergy(0.0))

    primaryspecies.append(Species()
        .withName("OH-")
        .withFormula("OH-")
        .withCharge(-1.0)
        .withElements(ElementalComposition([ (elements.get("O"), 1.0), (elements.get("H"), 1.0) ]))
        .withAggregateState(AggregateState.Aqueous)
        .withStandardGibbsEnergy(0.0))

    # Create the secondary species (using formation reactions from primary species)

    secondaryspecies = []

    secondaryspecies.append(Species()
        .withName("H2O")
        .withFormula("H2O")
        .withElements(ElementalComposition([ (elements.get("H"), 2.0), (elements.get("O"), 1.0) ]))
        .withAggregateState(AggregateState.Aqueous)
        .withFormationReaction(
            FormationReaction()
                .withReactants([
                    (primaryspecies.get("H+") , 1.0),
                    (primaryspecies.get("OH-"), 1.0) ])
                .withEquilibriumConstant(14.0)))

    return primaryspecies + secondaryspecies  # concatenate both species lists

We now create our auxiliary createCustomDatabase function that creates our Database object:

def createCustomDatabase():
    elements = createElementList()
    species = createSpeciesList(elements)
    return Database(species.data())

We are now ready for our chemical equilibrium calculation:

db = createCustomDatabase()

solution = AqueousPhase()

system = ChemicalSystem(db, solution)

state = ChemicalState(system)
state.set("H2O", 1.0, "kg")

equilibrate(state)

print(state)

+-----------------+------------+------+
| Property        |      Value | Unit |
+-----------------+------------+------+
| Temperature     |   298.1500 |    K |
| Pressure        |     1.0000 |  bar |
| Charge:         | 0.0000e+00 |  mol |
| Element Amount: |            |      |
| :: H            | 1.1111e+02 |  mol |
| :: O            | 5.5556e+01 |  mol |
| Species Amount: |            |      |
| :: H+           | 1.0000e-07 |  mol |
| :: OH-          | 1.0000e-07 |  mol |
| :: H2O          | 5.5556e+01 |  mol |
+-----------------+------------+------+

That’s it! We managed to create a Reaktoro script entirely from scratch. From database creation to an equilibrium calculation.

Creating thermodynamic databases

Contents

Creating thermodynamic databases#

A very simple custom thermodynamic database in YAML format#

A more complicated custom thermodynamic database in YAML format#

A very simple custom thermodynamic database using Python API#