Thermodynamic Backends

There are excellent and mature open-source software for modeling chemical systems for which development and maintenance has been ongoing for several years to a few decades. Wouldn’t it be great if Reaktoro could use the best these software can offer in terms of thermodynamic databases and activity models as thermodynamic backends, while relying on Reaktoro’s numerical algorithms for the intense chemical reaction calculations? Well, it turns out this is already supported for two widely used geochemical modeling codes: PHREEQC and GEMS.


PHREEQC is a computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations developed by USGS, United States.


GEMS is a Gibbs energy minimization software for geochemical modeling developed at Paul Scherrer Institute, Switzerland.


Neither PHREEQC nor GEMS are used in Reaktoro for solving the underlying mathematical problems for chemical equilibrium and kinetics. These backends act as providers of thermodynamic properties of species, phases, and reactions (e.g., activity coefficients, activities, standard chemical potentials, phase molar volume and enthalpy, equilibrium constant of reaction). Whenever Reaktoro’s numerical algorithms need these properties, the thermodynamic backend is invoked to retrieve them in a very efficient way by directly using its API. It seems complicated, but rest assure, the usage is rather simple as you’ll see below!


The code below demonstrates the combined use of Reaktoro and PHREEQC to perform a chemical equilibrium calculation in which PHREEQC thermodynamic data and activity models are used together with Reaktoro’s Gibbs energy minimization algorithm.

from reaktoro import *

# This string defines a PHREEQC script problem.
# This problem was taken from the official PHREEQC example named ex1.
ex1 = r'''(
TITLE Example 1.--Add uranium and speciate seawater.
        units   ppm
        pH      8.22
        pe      8.451
        density 1.023
        temp    25.0
        redox   O(0)/O(-2)
        Ca              412.3
        Mg              1291.8
        Na              10768.0
        K               399.1
        Fe              0.002
        Mn              0.0002  pe
        Si              4.28
        Cl              19353.0
        Alkalinity      141.682 as HCO3
        S(6)            2712.0
        N(5)            0.29    gfw   62.0
        N(-3)           0.03    as    NH4
        U               3.3     ppb   N(5)/N(-3)
        O(0)            1.0     O2(g) -0.7
        U       U+4     0.0     238.0290     238.0290
        U(4)    U+4     0.0     238.0290
        U(5)    UO2+    0.0     238.0290
        U(6)    UO2+2   0.0     238.0290
        #primary master species for U
        #is also secondary master species for U(4)
        U+4 = U+4
                log_k          0.0
        U+4 + 4 H2O = U(OH)4 + 4 H+
                log_k          -8.538
                delta_h        24.760 kcal
        U+4 + 5 H2O = U(OH)5- + 5 H+
                log_k          -13.147
                delta_h        27.580 kcal
        #secondary master species for U(5)
        U+4 + 2 H2O = UO2+ + 4 H+ + e-
                log_k          -6.432
                delta_h        31.130 kcal
        #secondary master species for U(6)
        U+4 + 2 H2O = UO2+2 + 4 H+ + 2 e-
                log_k          -9.217
                delta_h        34.430 kcal
        UO2+2 + H2O = UO2OH+ + H+
                log_k          -5.782
                delta_h        11.015 kcal
        2UO2+2 + 2H2O = (UO2)2(OH)2+2 + 2H+
                log_k          -5.626
                delta_h        -36.04 kcal
        3UO2+2 + 5H2O = (UO2)3(OH)5+ + 5H+
                log_k          -15.641
                delta_h        -44.27 kcal
        UO2+2 + CO3-2 = UO2CO3
                log_k          10.064
                delta_h        0.84 kcal
        UO2+2 + 2CO3-2 = UO2(CO3)2-2
                log_k          16.977
                delta_h        3.48 kcal
        UO2+2 + 3CO3-2 = UO2(CO3)3-4
                log_k          21.397
                delta_h        -8.78 kcal
        UO2 + 4 H+ = U+4 + 2 H2O
        log_k          -3.490
        delta_h        -18.630 kcal

# Initialize a Phreeqc instance with the official phreeqc.dat database file
phreeqc = Phreeqc('../../databases/phreeqc/phreeqc.dat')

# Execute a PHREEQC script defining a geochemical problem.
# Here this script is actually embedded into a string named `ex1`.
# However, `ex1` could also be a string containing the path to a script file.
# Method execute will automatically identify when the contents are embedded in
# the string and when the string is actually a path to a script file.

# Initialize a ChemicalSystem instance using the current state of the Phreeqc
# instance. This will allow the use of both PHREEQC thermodynamic data and
# PHREEQC activity models in the subsequent equilibrium calculations using
# Reaktoro's algorithms.
system = ChemicalSystem(phreeqc)

# Initialize an ChemicalState instance using the current state of the
# Phreeqc instance.
state = phreeqc.state(system)

# Output the equilibrium state calculated by PHREEQC to a file.

# Define an equilibrium problem in which the current state is mixed with 1
# mmol of HCl.
problem = EquilibriumProblem(system)
problem.add('HCl', 1.0, 'mmol')

# Set Hessian of Gibbs energy to an approximation, since PHREEQC does not 
# compute molar derivatives of activities 
options = EquilibriumOptions()
options.hessian = GibbsHessian.Approximation

# Calculate the new equilibrium state of the system.
# This will use both PHREEQC thermodynamic data and PHREEQC activity models.
state = ChemicalState(system)
equilibrate(state, problem, options)

# Print the new equilibrium state and check with pH is more acidic now.
Python and C++ files for this demo:

GEMS Backend

Similarly, the code below briefly demonstrates how Reaktoro and GEMS can be used together. You’ll need first to prepare your chemical system definition using GEM-Selektor, the graphical user interface of GEMS. In this step, you’ll be able to select which GEMS’ supported thermodynamic database you want to use as well as the activity models for each phase (aqueous, gaseous, solid solutions). Next, export the GEMS project files to disk, and use it in Reaktoro as shown below.

from reaktoro import *

# **Note:**
# This demo should be executed in the same directory where the script is located.
# Example:
# ~~~
# cd Reaktoro/demos/python
# python
# ~~~

# Use an exported project file from GEMS to initialize a Gems object,
gems = Gems("../resources/gems/CalciteBC-dat.lst")

# and then use it to construct the ChemicalSystem object.
system = ChemicalSystem(gems)

# Create a ChemicalState object that contains the temperature, pressure,
# and amounts of species stored in the exported GEMS file.
state = gems.state(system)

# Output the equilibrium state calculated by GEMS to a file.

# Perturb the equilibrium state calculated by GEMS
state.setSpeciesAmount("CO2@", 0.1)

# and then equilibrate the modified chemical state using Reaktoro's methods.

# Output the updated equilibrium state to a file.
Python and C++ files for this demo:

What about more thermodynamic backends?

Are there other chemical reaction modeling software that you think could be integrated with Reaktoro as thermodynamic backends? Let us know by creating a new issue at Reaktoro’s GitHub Issues.


It would be great if you could contribute to expanding the list of supported Reaktoro’s thermodynamic backends. Contributions can be made in several forms, ranging from direct code contribution to financing a project in which one or more experts will implement this.