Simplified Clinker Dissolution

This example illustrates ChemistryLab's equilibrium workflow applied to cement chemistry: computing the hydration state of a simplified Portland clinker at a given water-to-cement ratio. ChemistryLab computes thermodynamic equilibrium by minimising the Gibbs free energy of the system subject to element-conservation constraints. The workflow always follows the same four steps:

Build a ChemicalSystem (species + stoichiometric matrix).
Create an initial ChemicalState (temperature, pressure, initial amounts).
Call equilibrate (or use EquilibriumSolver explicitly).
Inspect the resulting ChemicalState.

Minimal workflow

The convenience function equilibrate handles everything with sensible defaults. The example below reproduces a simplified clinker dissolution calculation.

using ChemistryLab
using DynamicQuantities

substances = build_species("../../../data/cemdata18-thermofun.json")
input_species = split("C3S C2S C3A C4AF Gp Anh Portlandite Jennite H2O@ ettringite monosulphate12 C3AH6 C3FH6 C4FH13")
species = speciation(substances, input_species; aggregate_state=[AS_AQUEOUS])

cs = ChemicalSystem(species, CEMDATA_PRIMARIES)

state = ChemicalState(cs)

# Clinker + gypsum composition (mass fractions, total = 1)
compo = ["C3S" => 0.678, "C2S" => 0.166, "C3A" => 0.04, "C4AF" => 0.072, "Gp" => 0.028]
c     = sum(last.(compo))
wc    = 0.4          # water-to-cement ratio
w     = wc * c
mtot  = c + w

for x in compo
    set_quantity!(state, x.first, x.second / mtot * u"kg")
end
set_quantity!(state, "H2O@", w / mtot * u"kg")

# pH-neutral seed for H⁺ and OH⁻
V = volume(state)
set_quantity!(state, "H+",  1e-7u"mol/L" * V.liquid)
set_quantity!(state, "OH-", 1e-7u"mol/L" * V.liquid)

state_eq = equilibrate(state)

ChemicalState{Species, AbstractReaction, DynamicQuantities.Quantity{Float64, DynamicQuantities.Dimensions{DynamicQuantities.FRInt32}}}
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│             T : 298.15 K                                                                         │
│             P : 1.0 bar                                                                          │
╞══════════════════════════════════════════════════════════════════════════════════════════════════╡
│    # liquid #│             n [mol]│               m [g]│             V [cm³]│           c [mol/L]│
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│   tot. liquid│             4.56714│             82.3078│             82.4392│                    │
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│          H2O@│              4.5639│             82.2186│             82.4623│             55.3608│
│           OH-│          0.00203298│            0.034575│         -0.00957097│           0.0246604│
│          Ca+2│         0.000846461│           0.0339245│          -0.0156077│           0.0102677│
│         CaOH+│         0.000347624│           0.0198441│          0.00200317│          0.00421672│
│         AlO2-│          7.28378e-6│         0.000429594│          6.89581e-5│          8.83533e-5│
│       CaSiO3@│          3.11949e-6│          0.00036236│                 0.0│          3.78399e-5│
│        HSiO3-│          1.08066e-7│          8.33079e-6│          4.89091e-7│          1.31085e-6│
│      Ca(SO4)@│          6.85274e-8│           9.3289e-6│          3.22122e-7│          8.31247e-7│
│         SO4-2│          4.76342e-8│          4.57555e-6│          6.15321e-7│          5.77809e-7│
│    Ca(HSiO3)+│          4.10167e-8│          4.80584e-6│         -2.76311e-7│          4.97538e-7│
│         FeO2-│          3.85934e-8│          3.39016e-6│          1.74269e-8│          4.68143e-7│
│        SiO3-2│          3.56706e-8│          2.71389e-6│                 0.0│          4.32689e-7│
│         SiO2@│          1.28592e-8│          7.72619e-7│          2.06561e-7│          1.55984e-7│
│        AlO2H@│          8.08311e-9│          4.84886e-7│          1.05157e-7│          9.80493e-8│
│        FeO2H@│          7.41608e-9│          6.58926e-7│          5.34649e-8│          8.99582e-8│
│      AlSiO5-3│          6.93203e-9│          9.36251e-7│         -2.36602e-7│          8.40866e-8│
│            H+│          4.70026e-9│          4.73786e-9│                 0.0│          5.70149e-8│
│           O2@│          3.24895e-9│           1.0396e-7│          9.90957e-8│          3.94102e-8│
│          FeO+│          2.89581e-9│          2.08047e-7│         -1.21684e-7│          3.51266e-8│
│          AlO+│          2.66921e-9│          1.14724e-7│         8.20344e-10│          3.23779e-8│
│         HSO4-│          2.51208e-9│          2.43833e-7│          8.75239e-8│           3.0472e-8│
│      Si4O10-4│           1.7528e-9│          4.77341e-7│                 0.0│          2.12618e-8│
│        AlOH+2│          1.55229e-9│          6.82831e-8│         -4.23389e-9│          1.88295e-8│
│         FeOH+│           1.4532e-9│          1.05868e-7│         -2.42871e-8│          1.76275e-8│
│        FeOH+2│          1.42639e-9│          1.03916e-7│         -3.61455e-8│          1.73024e-8│
│          Fe+2│          1.24226e-9│          6.93742e-8│         -2.81254e-8│          1.50688e-8│
│          Al+3│          1.05834e-9│          2.85556e-8│         -4.78823e-8│          1.28378e-8│
│      Fe(SO4)@│           1.0129e-9│          1.53861e-7│          1.69344e-9│          1.22866e-8│
│      Al(SO4)+│         9.42937e-10│          1.16017e-7│         -5.67557e-9│           1.1438e-8│
│          Fe+3│         8.96998e-10│          5.00929e-8│         -3.38978e-8│          1.08807e-8│
│      Fe(SO4)+│         8.16264e-10│          1.23991e-7│         -2.15121e-9│           9.9014e-9│
│     Al(SO4)2-│         7.97514e-10│           1.7473e-7│          2.48121e-8│          9.67396e-9│
│           H2@│         7.89022e-10│          1.59067e-9│          1.99341e-8│          9.57096e-9│
│     Fe(SO4)2-│         6.91611e-10│           1.7149e-7│           2.1084e-8│          8.38935e-9│
│     Fe(HSO4)+│         6.85214e-10│          1.04775e-7│          1.28864e-8│          8.31175e-9│
│         SO3-2│         6.43064e-10│          5.14818e-8│         -2.64684e-9│          7.80046e-9│
│    Fe(HSO4)+2│         5.83228e-10│          8.91808e-8│          1.35301e-9│          7.07464e-9│
│         HSO3-│         5.65897e-10│          4.58745e-8│          1.86502e-8│          6.86442e-9│
│           HS-│         2.02334e-10│          6.69078e-9│          4.08907e-9│          2.45434e-9│
│          H2S@│         1.93539e-10│          6.59503e-9│          6.76442e-9│          2.34766e-9│
│        S2O3-2│         1.79909e-10│          2.01709e-8│           4.9641e-9│          2.18232e-9│
╞══════════════════════════════════════════════════════════════════════════════════════════════════╡
│     # solid #│             n [mol]│               m [g]│             V [cm³]│                    │
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│    tot. solid│              6.0689│             917.692│             387.255│                    │
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│   Portlandite│             2.89109│             214.207│             95.5796│                    │
│       Jennite│             2.85524│             546.425│             223.851│                    │
│monosulphate12│            0.118057│             73.4911│             36.6095│                    │
│         C3FH6│            0.107551│             46.8929│             16.7012│                    │
│         C3AH6│           0.0969547│             36.6761│             14.5143│                    │
│        C4FH13│          3.22011e-8│          1.99064e-5│          9.20758e-6│                    │
│           C2S│          1.46253e-8│          2.51902e-6│          7.57446e-7│                    │
│    ettringite│          1.34178e-8│          1.68403e-5│          9.48675e-6│                    │
│            Gp│          5.83184e-9│          1.00403e-6│           4.3558e-7│                    │
│           Anh│          5.53484e-9│           7.5348e-7│          2.54271e-7│                    │
│           C3S│          1.76753e-9│          4.03551e-7│          1.29348e-7│                    │
│          C4AF│         7.72976e-10│          3.75632e-7│          1.00643e-7│                    │
│           C3A│         6.80971e-10│          1.83992e-7│          6.07542e-8│                    │
╞══════════════════════════════════════════════════════════════════════════════════════════════════╡
│     # TOTAL #│             n [mol]│               m [g]│             V [cm³]│                    │
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│              │              10.636│              1000.0│             469.695│                    │
╞══════════════════════════════════════════════════════════════════════════════════════════════════╡
│            pH : 12.3921                                                                          │
│           pOH : 1.608                                                                            │
│      porosity : 0.175517                                                                         │
│    saturation : 1.0                                                                              │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

Quick shortcut

If your ChemicalState was constructed from a ChemicalSystem built with CEMDATA_PRIMARIES, calling equilibrate(state) with no other arguments is usually sufficient for cement-chemistry problems.

Inspecting the equilibrium state

The returned ChemicalState carries all derived thermodynamic quantities:

println("pH      = ", pH(state_eq))
println("pOH     = ", pOH(state_eq))
println("porosity   = ", porosity(state_eq))
println("saturation = ", saturation(state_eq))

Phase volumes and mole amounts are accessible via named tuples:

v = volume(state_eq)
println("V liquid = ", v.liquid)
println("V solid  = ", v.solid)
println("V total  = ", v.total)

m = moles(state_eq)
println("n liquid = ", m.liquid)
println("n solid  = ", m.solid)

Individual species amounts (in mol):

cs_eq = state_eq.system
for (i, sp) in enumerate(cs_eq.species)
    n_i = state_eq.n[i]
    println(rpad(symbol(sp), 20), ustrip(n_i), " mol")
end

Controlling the solver

Variable space: `:linear` vs `:log`

equilibrate accepts a variable_space keyword that selects the optimisation variable space:

`variable_space`	Variables	Recommended when
`Val(:linear)`	mole amounts `nᵢ ≥ 0`	most systems, default
`Val(:log)`	`log nᵢ`	systems spanning many orders of magnitude

state_eq_log = equilibrate(state; variable_space=Val(:log))

Convergence

Solving a system of equations in chemistry can be a difficult undertaking. The orders of magnitude can vary greatly, and convergence is not guaranteed.

Tolerances

Tighter tolerances are passed directly as keyword arguments and forwarded to the underlying Ipopt solver:

state_eq_tight = equilibrate(state; abstol=1e-12, reltol=1e-12)

Using `EquilibriumSolver` explicitly

For batch calculations where many different initial states share the same system and activity model, construct an EquilibriumSolver once and reuse it:

using OptimizationIpopt

opt = IpoptOptimizer(
    acceptable_tol        = 1e-12,
    dual_inf_tol          = 1e-12,
    acceptable_iter       = 1000,
    constr_viol_tol       = 1e-12,
    warm_start_init_point = "no",
)

solver = EquilibriumSolver(
    cs,
    DiluteSolutionModel(),
    opt;
    variable_space = Val(:linear),
    abstol  = 1e-10,
    reltol  = 1e-10,
)

Once built, solver is called with any compatible ChemicalState:

using OptimizationIpopt #hide

opt = IpoptOptimizer( #hide
    acceptable_tol        = 1e-12, #hide
    dual_inf_tol          = 1e-12, #hide
    acceptable_iter       = 1000, #hide
    constr_viol_tol       = 1e-12, #hide
    warm_start_init_point = "no", #hide
) #hide

solver = EquilibriumSolver( #hide
    cs, #hide
    DiluteSolutionModel(), #hide
    opt; #hide
    variable_space = Val(:linear), #hide
    abstol  = 1e-10, #hide
    reltol  = 1e-10, #hide
) #hide
state_eq2 = solve(solver, state)

Performance

The potential function μ(n, p) is compiled once during EquilibriumSolver construction. Repeated calls to solve(solver, ...) with different states reuse it, avoiding redundant compilation overhead.

Activity models

All activity models inherit from AbstractActivityModel. The only built-in model is DiluteSolutionModel, which implements:

Phase	Law	Expression
Solvent (H₂O)	Raoult	`ln a = ln xₛ`
Aqueous solutes	Henry	`ln a = ln(cᵢ / c°)`, `c° = 1 mol/L`
Crystals	Pure solid	`ln a = 0`
Gas	Ideal mixture	`ln a = ln xᵢ`

To implement a custom activity model, define a new subtype and extend activity_model:

struct MyModel <: AbstractActivityModel
    # model parameters
end

function ChemistryLab.activity_model(::ChemicalSystem, ::MyModel)
    # return a closure lna(n, p) -> Vector{Float64}
    return (_, _) -> begin
        # compute and return log-activities
    end
end

Pass your model to equilibrate or EquilibriumSolver:

state_eq = equilibrate(state; model=MyModel(...))