Skip to content

ChemistryLab.jlChemistry you can script

Formulas, species, reactions, thermodynamic databases and equilibrium โ€” for aqueous geochemistry and cement chemistry.

ChemistryLab

What it does โ€‹

ChemistryLab handles chemical formulas, species and reactions as first-class objects, reads thermodynamic data from ThermoFun and Cemdata, and solves equilibrium by Gibbs-energy minimization. Kinetics and chemo-mechanical coupling build on the same objects.

It is written for work that has to be reproducible and scripted: aqueous geochemistry, cement chemistry, and any problem where speciation, a database and a solver have to be driven from code rather than from a dialog box.

A first calculation โ€‹

Calcite in water. The species come from a thermodynamic database, four of them are declared as the independent basis, and the equilibrium state follows from a Gibbs-energy minimization:

julia
using ChemistryLab, DynamicQuantities
using OptimaSolver          # the default equilibrium back-end

species = speciation(build_species(datapath("cemdata18-thermofun.json"); verbose = false),
                     split("Cal H2O@ CO2");
                     aggregate_state = [AS_AQUEOUS],
                     exclude_species = split("H2@ O2@ CH4@"))
byname = Dict(symbol(s) => s for s in species)

system = ChemicalSystem(collect(values(byname)),
                        [byname[s] for s in split("H2O@ H+ CO3-2 Ca+2")])

state = ChemicalState(system)
set_quantity!(state, "Cal", 1e-3u"mol")     # 1 mmol of calcite
set_quantity!(state, "H2O@", 1.0u"kg")      # in 1 kg of water
V = volume(state)
set_quantity!(state, "H+",  1e-4u"mol/L" * V.liquid)
set_quantity!(state, "OH-", 1e-10u"mol/L" * V.liquid)

equilibrated = equilibrate(state)
โ”Œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”
โ”‚  Loading database: data/cemdata18-thermofun.json  โ”‚
โ””โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”˜
โ”Œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”
โ”‚  Building species  โ”‚
โ””โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”˜

Progress:  44%|โ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–                                                                                         |  ETA: 0:00:00
Progress:  84%|โ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–‹                         |  ETA: 0:00:00
Progress: 100%|โ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆโ–ˆ| Time: 0:00:00
โ”Œ Warning: Verbosity toggle: missing_second_order_ad
โ”‚  The selected optimization algorithm requires second order derivatives, but `SecondOrder` ADtype was not provided. So a `SecondOrder` with AutoForwardDiff() for both inner and outer will be created, this can be suboptimal and not work in some cases so an explicit `SecondOrder` ADtype is recommended.
โ”” @ OptimizationBase ~/.julia/packages/OptimizationBase/BEgWO/src/cache.jl:116

The state carries everything derived from it โ€” pH, phase volumes, individual amounts:

julia
using Printf
@printf("pH                = %.2f\n", pH(equilibrated))
@printf("dissolved Ca(2+)  = %.3e mol\n", ustrip(moles(equilibrated, "Ca+2")))
@printf("remaining calcite = %.3e mol\n", ustrip(moles(equilibrated, "Cal")))
pH                = 9.53
dissolved Ca(2+)  = 1.552e-04 mol
remaining calcite = 8.389e-04 mol

Getting started takes the same problem more slowly, and shows how the solubility constant is obtained analytically before any solver is involved.

Reading paths โ€‹

The documentation is organized by the question each chapter answers. Theory explains why a calculation is the right one, the Manual describes how each kind of object is written, the Tutorials drive a calculation from start to finish, and the Applications show complete cases and what the modeling choices cost in numbers. Three entry points follow from it, depending on what the reader already knows.

New to chemical thermodynamics. The chapter Theory is written for this reader and is best read in the order it states, beginning with Thermochemistry and Standard states, which fix the notation and the conventions every other page relies on. The Manual then introduces the objects one at a time, from Species onwards, and the tutorial Chemical Equilibrium puts them together.

Knowing what is to be computed. Getting started and the tutorial Chemical Equilibrium are enough to write a first calculation. The aqueous applications, beginning with COโ‚‚ dissolution and the carbonate system, are small enough to check by hand and are the place to test one's understanding before a larger system.

Holding the analysis of a cement. The route from an oxide analysis to a chemical system is laid out in Bogue calculation and Choosing the species list; the binders of EN 197-1 are then worked in increasing order of difficulty, from A CEM I from its clinker phases to the composite cements, and The binders, and what distinguishes them says which of the package's models each family requires.

Where this comes from โ€‹

This package would not exist without two bodies of work that came before it, and it is worth saying so on the first page rather than in a footnote.

GEM-Selektor and GEMS3K, from the Paul Scherrer Institute and Empa (Kulik et al., 2013), established the Gibbs energy minimization approach used here, and much of the vocabulary with it โ€” the phase stability index this package computes as is the same quantity as their , reached from the same KKT conditions. CEMDATA18 (Lothenbach et al., 2019), the thermodynamic database behind every cement calculation in this manual, is their laboratory's work and ships here unchanged; the zeolite extension is transcribed from two further papers by the same group. Nothing here would produce a number without it.

Reaktoro (Leal et al., 2017), by Allan Leal and contributors, is both an ancestor and a reference: parts of the thermodynamics and kinetics here are Julia ports of its C++ implementation, and it is the code this package checks itself against throughout โ€” see the validation page. Where a result differs, the burden of proof has been on us.

Both are mature, carefully built and widely used, and both address a wider range of problems than this package attempts. What ChemistryLab tries to add is narrower: a Julia-native formulation in which an equilibrium comes with a proof of its optimality rather than a converged iterate, differentiable end to end so that a calibration can be posed as an optimization, and with the cementitious special cases โ€” cement chemist notation, Bogue, the oxide-budget entry route for a glass, the binder families of EN 197-1 โ€” treated as first-class rather than as an application layer.

That is an addition to their work, not a substitute for it.