Creating thermodynamic databases

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

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 thermodynamic reaction models to meet this need:

  • ReactionThermoModelConstLgK

  • ReactionThermoModelVantHoff

  • ReactionThermoModelPhreeqcLgK

  • ReactionThermoModelGemsLgK

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 stated: 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:

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-
    ReactionThermoModel:
      ConstLgK:
        lgKr: 2.0

Attention

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.

Notes on the above database:

  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.

  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.

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

  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.

  5. The AB species calculates its standard thermodynamic properties with a constant lgK reaction model implemented in ReactionThermoModelConstLgK. 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.

  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.yml, which is used next:

from reaktoro import *

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

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.15 |    K |
| Pressure        | 100000 |   Pa |
| Charge:         |      0 |  mol |
| Element Amount: |        |      |
| :: A            |      1 |  mol |
| :: B            |      1 |  mol |
| Species Amount: |        |      |
| :: A+           |  1e-16 |  mol |
| :: B-           |  1e-16 |  mol |
| :: AB           |      1 |  mol |
+-----------------+--------+------+
FINAL STATE
+-----------------+-----------+------+
| Property        |     Value | Unit |
+-----------------+-----------+------+
| Temperature     |    298.15 |    K |
| Pressure        |    100000 |   Pa |
| Charge:         |         0 |  mol |
| Element Amount: |           |      |
| :: A            |         1 |  mol |
| :: B            |         1 |  mol |
| Species Amount: |           |      |
| :: A+           | 0.0995037 |  mol |
| :: B-           | 0.0995037 |  mol |
| :: AB           |  0.900496 |  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:

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: H2O
  Formula: H2O
  Elements: 2:H 1:O
  AggregateState: Aqueous
  FormationReaction:
    Reactants: 1:H+ 1:OH-
    ReactionThermoModel:
      ConstLgK:
        lgKr: 14.0
- Name: HCO3-
  Formula: HCO3-
  Elements: 1:H 1:C 3:O
  Charge: -1.0
  AggregateState: Aqueous
  FormationReaction:
    Reactants: 1:CO3-2 1:H+
    ReactionThermoModel:
      VantHoff:
        lgKr: 10.329
        dHr: -14899.224  # in J/mol
- Name: CO2
  Formula: CO2
  Elements: 1:H 1:C 3:O
  AggregateState: Aqueous
  FormationReaction:
    Reactants: 1:CO3-2 2:H+ -1:H2O
    ReactionThermoModel:
      PhreeqcLgK:
        A1: 464.1965
        A2: 0.09344813
        A3: -26986.16
        A4: -165.75951
        A5: 2248628.9
        A6: 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: Calcite
  Formula: CaCO3
  Elements: 1:Ca 1:C 3:O
  AggregateState: Solid
  FormationReaction:
    Reactants: 1:CO3-2 1:Ca+2
    ReactionThermoModel:
      VantHoff:
        lgKr: 8.48
        dHr: 10832.376  # in J/mol

We will now use this database, saved locally under the name db-geo-complicated.yml, to perform a chemical equilibrium calculation in an initial chemical state at 25 °C and 1 bar containing 1 kg of H2O, 0.1 mol of CO2, and 1 mol of mineral calcite (CaCO3):

from reaktoro import *

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

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.15 |    K |
| Pressure        |       100000 |   Pa |
| Charge:         | -9.71445e-17 |  mol |
| Element Amount: |              |      |
| :: H            |      111.017 |  mol |
| :: C            |          1.1 |  mol |
| :: O            |      58.7084 |  mol |
| :: Ca           |            1 |  mol |
| Species Amount: |              |      |
| :: H2O          |      55.4995 |  mol |
| :: H+           |  2.26067e-06 |  mol |
| :: OH-          |  4.46631e-09 |  mol |
| :: HCO3-        |    0.0178697 |  mol |
| :: CO2          |     0.091064 |  mol |
| :: CO3-2        |  3.70518e-07 |  mol |
| :: Ca+2         |   0.00893409 |  mol |
| :: Calcite      |     0.991066 |  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.15 |    K |
| Pressure        |       100000 |   Pa |
| Charge:         | -1.04083e-16 |  mol |
| Element Amount: |              |      |
| :: H            |      111.012 |  mol |
| :: C            |          1.1 |  mol |
| :: O            |      58.7062 |  mol |
| :: Ca           |            1 |  mol |
| Species Amount: |              |      |
| :: H2O          |      55.4973 |  mol |
| :: H+           |  2.26105e-06 |  mol |
| :: OH-          |  4.46515e-09 |  mol |
| :: HCO3-        |    0.0178712 |  mol |
| :: CO2          |    0.0910633 |  mol |
| :: CO3-2        |  3.70591e-07 |  mol |
| :: Ca+2         |   0.00893484 |  mol |
| :: Calcite      |     0.991065 |  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.15 |    K |
| Pressure        |      100000 |   Pa |
| Charge:         |           0 |  mol |
| Element Amount: |             |      |
| :: H            |     111.111 |  mol |
| :: O            |     55.5556 |  mol |
| Species Amount: |             |      |
| :: H+           | 1.00085e-07 |  mol |
| :: OH-          | 1.00085e-07 |  mol |
| :: H2O          |     55.5556 |  mol |
+-----------------+-------------+------+

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