Solving an equilibrium

Before this page

The tutorial Chemical Equilibrium. The snippets below reuse its calcite system cs and its initial state state, and they are shown rather than executed.

equilibrate works with its defaults on most systems, and the tutorial Chemical Equilibrium relies on them throughout. This page gathers what those defaults are and how each is changed: the back ends and the way the answer is certified, the constraints that replace a fixed temperature and pressure, the derivatives of an equilibrium with respect to its inputs, the settings of the solver, the activity models of the aqueous phase and the declaration of solid solutions.

Choosing a solver

ChemistryLab provides two solver extensions. Load whichever fits your workflow:

Extension packages	Default solver	When to use
Optimization, OptimizationIpopt	IpoptOptimizer	general-purpose, robust
OptimaSolver	OptimaOptimizer	preferred when available

When both are loaded, OptimaSolver provides the default single back end, and on a cement equilibrium it is the faster of the two: 0.10 s against 0.28 s for Ipopt on a CEM I paste at w/c = 0.45. An earlier version of this page quoted "3 to 26 times faster"; that range is not reproduced — the measured factor there is 2.8, and on a small calcite system Ipopt is both faster (0.015 s against 0.038 s) and far more accurate on the element balance (3e-12 relative against 3e-2). Which back end is quicker depends on the system, so neither ordering is worth stating as a rule.

What is worth stating is that neither of them returns a KKT point on a cement. Both satisfy the element balance to 10⁻¹⁴–10⁻⁹ and agree on the pH to three decimals, while their stationarity residual sits at 133 and 136 in RT units. That does not matter for a pH, and it matters a great deal for a trace species or a saturation index — which is why equilibrate certifies by default, at a cost of about 1.4 s on the same paste.

equilibrate(state) does not pick one of them — it proves the answer

Neither back end is reliable on its own, so the one-argument equilibrate solves by every route available and keeps the answer optimality_certificate proves optimal. Measured, with the element balance judged row by row against each row's own budget:

case	interior point	dual Newton	certified route
calcite in water, 10–40 °C	3.0e-2	3.0e-12	certified
calcite + 1 mmol CO₂	1.0e-3	6.3e-13	certified
calcite + 50 mmol CO₂	1.3e-16	8.5e-14	certified
pure water	1.3e-16	1.3e-16	certified
CEM I paste, w/c 0.45 and 0.60	1.5e-14	2.0e-14	certified
CEM I paste, w/c 0.30	7.2e-16	3.3e-10, not certified	certified

Ten of ten certified through both routes; seven through the interior point alone, nine through the dual Newton alone. The interior point's 3 % on the first row is not a tolerance: it is a charge balance wrong in the second digit, and the reason it cannot be improved is that the fraction-to-boundary rule caps its step at 15 % of a correction the next iteration re-poses — traced over twenty-three iterations with the residual frozen at 3.0e-6 and ‖dn‖ decaying at 1 − α.

state_eq = equilibrate(state)                   # certified (the default)
state_eq = equilibrate(state; certify = false)  # one back end, as before
eq, cert = equilibrate_certified(state)         # when the proof itself is wanted

Pass a solver explicitly — equilibrate(state, OptimaOptimizer()) — to use that one back end and nothing else.

Explicit solver (always works)

Pass the solver as the second positional argument:

using Optimization, OptimizationIpopt
state_eq = equilibrate(state, IpoptOptimizer())

using OptimaSolver
state_eq = equilibrate(state, OptimaOptimizer())

Default shortcut

When only one extension is loaded:

using Optimization, OptimizationIpopt
state_eq = equilibrate(state)   # → IpoptOptimizer

using OptimaSolver
state_eq = equilibrate(state)   # → OptimaOptimizer

With both loaded, OptimaSolver always wins:

using Optimization, OptimizationIpopt
using OptimaSolver
state_eq = equilibrate(state)   # → OptimaOptimizer (priority)

Constraints other than fixed T and P

A closed system at given temperature and pressure is one case among several. A vessel that exchanges no heat has its temperature determined by the reaction; a titration holds the pH and lets the amount of acid follow. Both are equilibrium problems, and both are stated by saying what is held and what is unknown.

Needs OptimaSolver 0.5

The blocks in this section are shown rather than executed: the documentation environment resolves OptimaSolver from the registry, and the parameter block these constraints ride on arrived in 0.5.0. The figures quoted below are the ones test/equilibrium_constraints.jl asserts.

using ChemistryLab, DynamicQuantities, OptimaSolver

# The same calcite system, but adiabatic: no heat leaves, so T is an unknown
st = ChemicalState(cs)
set_quantity!(st, "Cal",  1e-2u"mol")
set_quantity!(st, "H2O@", 1.0u"kg")
set_quantity!(st, "H+",   1e-2u"mol")
set_quantity!(st, "OH-",  1e-10u"mol")

des   = DualEquilibriumSolver(cs, DiluteSolutionModel())
eq_ad = solve(des, st; constraint = Adiabatic())

temperature(eq_ad)      # 298.1996 K — the reaction warms it by 0.0496 K
pH(eq_ad)               # 6.4323, against 6.4329 at fixed temperature

constraint	held	unknown	vehicle
FixedTP (default)	T, P	—	—
Adiabatic	the enthalpy of the initial state	T	parameter
FixedEnthalpy	a prescribed H	T	parameter
FixedVolume	a prescribed V	P	parameter
SealedVolume	the volume of the initial state	P	parameter
FixedpH	−log₁₀ a(H⁺)	the titrant amount	column
FixedActivity	a of any species	the titrant amount	column

The two vehicles

They are not interchangeable, and the distinction is the same one Reaktoro draws.

A prescribed property — an enthalpy, a volume — adds one parameter and one equation to the solver's own square system. Nothing loops around the equilibrium solve: the temperature is an unknown beside the amounts and the element potentials, and the enthalpy balance is one more row of the same Newton system. So an adiabatic solve costs one equation, not a solve per trial temperature.

A prescribed chemical potential — a pH, an activity — adds one column to the conservation matrix instead: an unknown amount of a titrant the system may draw on. The linear rows become A n − A[:, titrant] q = b, so the system is open to that one substance and closed to everything else, and q comes back as part of the answer — it is the reagent consumed, which is what a titration measures.

q = Ref(Float64[])
eq_ph = solve(des, st; constraint = FixedpH(7.0), parameters = q)
q[][1]      # 3.4673e-3 mol of H+ had to be added to hold pH 7

On 10 mmol of calcite in a kilogram of water, whose free pH is 9.90, holding the pH takes more acid the lower the target, and at pH 6 the calcite is gone entirely:

pH held	H⁺ added (mol)	calcite left (mol)
6	1.689e-2	0
7	3.467e-3	7.063e-3
8	8.648e-4	9.149e-3
9	2.504e-4	9.728e-3

What is checked, and what is refused

The adiabatic answer is validated two ways. Against physics: H⁺ + OH⁻ → H₂O comes out at −55.85 kJ/mol at every amount tested, against the accepted −55.8, with nothing fitted to it — the enthalpies come from the database and the temperature is an unknown of the system. And against itself: solving at a fixed temperature equal to the one the adiabatic solve found returns the same composition to 1e-12.

A volume constraint on a condensed system is refused, with the lever it measured named in the error. The molar volumes of water and of the minerals in the shipped databases are exactly pressure-independent — V⁰(1 bar) = V⁰(100 bar) to the last bit for H2O@ and Cal — and only a few aqueous ions vary, OH⁻ by 8 % over 100 bar. The relative lever (∂V/∂P)·P/V is then about 1e-6, meaning some 9 600 bar to change the volume by one percent. That is the physics, not a solver limitation: the volume of an incompressible condensed system is fixed by its composition. Declare a gas phase to give the pressure something to do. To report the volume change of a sealed specimen at fixed pressure — what a hydrating binder actually needs — use porosity with a reference state instead.

Constraints need the certifying solver

The prescribed property or potential is an unknown of the dual Newton's own system, so a constraint other than FixedTP requires OptimaSolver and a system with an aqueous phase and H2O@. Asking for one otherwise raises, rather than being silently ignored. equilibrate(state; constraint = ...) routes through the certified path, so the answer still comes with its proof.

Differentiating an equilibrium

A ChemicalState carries whatever number type its amounts have, so a composition built from ForwardDiff.Dual values propagates through the speciation, and pH, pOH, porosity and saturation come back as duals too.

Crossing the solve works as well, and without asking any solver to iterate on dual numbers — Ipopt is a C library and never could. The equilibrium is solved once at the primal values, and the sensitivities come from the optimality conditions, the implicit-function-theorem route:

restricted to the species actually present. One factorization serves every partial derivative, and the result is exact — no step size to choose.

using ChemistryLab, DynamicQuantities, ForwardDiff, OptimaSolver

f(x) = begin
    n = Any[fill(0.0u"mol", length(cs.species))...]
    n[i_h2o] = 55.5u"mol";  n[i_cal] = 0.05u"mol";  n[i_co2] = x * u"mol"
    eq = equilibrate(ChemicalState(cs, n), OptimaOptimizer())
    ustrip(us"mol", eq.n[i_ca])
end

ForwardDiff.derivative(f, 0.01)     # → 0.15193

Why the complementarity conditions cannot be skipped

The stationarity conditions are ∇G − Aᵀy − z = 0, A n = b, nᵢzᵢ = 0, with z ≥ 0 the stability multipliers. The last block partitions the species, and dropping it does not degrade the answer gently — it destroys it. On calcite + CO₂ in water with a gas phase declared, the unreduced system puts the whole response into the absent gas species (n = 5×10⁻¹¹ mol), returning a sensitivity that satisfies the element balance to 4×10⁻¹⁶ and means nothing.

No back-end returns z, so the active set is recovered internally: a species negligible on the scale of the system that nonetheless takes a leading share of the response is pinned, and the system re-solved.

Verified against Reaktoro

On calcite + CO₂ + water, ∂n/∂(CO₂) from this route agrees with the package's own finite differences to 9×10⁻⁵ — the finite-difference truncation error — and with Reaktoro 2.13 reading the same Cemdata18 file, over the same eleven species, under the same (ideal) activity model:

species	ChemistryLab (AD)	Reaktoro (FD)	rel. diff
H₂O	−0.181336	−0.181397	3.4×10⁻⁴
Ca²⁺	+0.151920	+0.151987	4.4×10⁻⁴
HCO₃⁻	+0.333308	+0.333427	3.6×10⁻⁴
CO₂(aq)	+0.818655	+0.818600	6.7×10⁻⁵
Ca(HCO₃)⁺	+0.029333	+0.029338	1.6×10⁻⁴
calcite	−0.181251	−0.181324	4.0×10⁻⁴

The equilibrium amounts agree to the same order (Ca²⁺ 3.5281×10⁻³ against 3.52902×10⁻³). Reaktoro's own spread across h ∈ {10⁻³, 10⁻⁴, 10⁻⁵} is 7.2×10⁻⁴, so the residual difference sits below the oracle's truncation error on every species.

The absent gas species gets exactly zero from the active-set treatment, against 2×10⁻⁹ by finite differences.

A cross-code comparison has three knobs, not one

Database, species list and activity model all have to match, and each is worth tens of percent here. Reading Cemdata18 in both codes but leaving Reaktoro on its HKF activity model against this package's default DiluteSolutionModel() moves ∂Ca²⁺/∂(CO₂) from +0.1520 to +0.2179 — a 35 % gap that says nothing about either code.

The species list matters just as much. Dropping the aqueous calcium complexes leaves free Ca²⁺ as the only aqueous home for calcium, so ∂Ca²⁺/∂(CO₂) and ∂calcite/∂(CO₂) mirror each other exactly. Restoring CaOH⁺, Ca(CO₃)@ and Ca(HCO₃)⁺ breaks that mirror — the difference is what Ca(HCO₃)⁺ takes up — and both codes break it the same way. The element balance closes to 2×10⁻¹⁶ either way.

Controlling the solver

Variable space: :linear vs :log

equilibrate accepts a variable_space keyword that selects the optimization 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)

Activity models

All activity models inherit from AbstractActivityModel. Three built-in models are provided, covering ideal behavior through to the extended Debye-Hückel level used by standard geochemical codes.

Choosing a model

Model	Formula	Valid range	Parameters needed
DiluteSolutionModel	Raoult / Henry	I ≪ 1 mol/kg	none
HKFActivityModel	B-dot extended Debye-Hückel	I ≲ 1 mol/kg	A, B, Ḃ (defaults at 25 °C)
DaviesActivityModel	Davies equation	I ≲ 0.5 mol/kg	A, b (defaults at 25 °C)

---

DiluteSolutionModel (ideal dilute solution)

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ᵢ

state_eq = equilibrate(state)   # DiluteSolutionModel is the default

---

HKFActivityModel (extended Debye-Hückel B-dot)

Implements the extended Debye-Hückel model of Helgeson (Helgeson, 1969) and Helgeson, Kirkham & Flowers (Helgeson et al., 1981), identical to the model used by PHREEQC (Parkhurst and Appelo, 2013) and EQ3/6.

Ion activity coefficient:

log₁₀ γᵢ = −A zᵢ² √I / (1 + B åᵢ √I)  +  Ḃ I

Neutral aqueous species (salting-out):

log₁₀ γᵢ = Kₙ I

Water activity is computed from the osmotic coefficient via Gibbs-Duhem (not Raoult), which is accurate up to I ≈ 1 mol/kg.

Ionic radius lookup (priority order):

1. model.å — one common radius for every ion, when given. Short-circuits the rest of the chain.

2. sp[:å] — explicit value set in species properties.

3. REJ_HKF — Helgeson et al. (1981) Table 3 (27 common ions) (Helgeson et al., 1981).

4. REJ_CHARGE_DEFAULT — fallback by formal charge (Xu et al., 2011).

5. model.å_default (default: 3.72 Å).

å_default does not impose a common ionic radius

It is the last resort of the chain above, reached only for a charge that neither table covers — in practice |z| ≥ 5. Setting it changes essentially nothing for a real solution. Pass å to impose one common radius, which is what GEM-Selektor, PHREEQC's -gamma and most published cement models actually use.

Usage:

# Fixed A, B at 25 °C / 1 bar (fast — suitable for isothermal calculations)
state_eq = equilibrate(state; model=HKFActivityModel())

# Temperature-dependent A and B (recomputed from T, P at each equilibrium solve)
state_eq = equilibrate(state; model=HKFActivityModel(temperature_dependent=true))

# Custom parameters
model = HKFActivityModel(A=0.52, B=0.33, Ḃ=0.04)

# One common ion size of 3.72 Å, overriding the per-species tables
model = HKFActivityModel(å = 3.72)

# å = 0 collapses the denominator to 1: the Debye-Hückel limiting law plus Ḃ I
model = HKFActivityModel(å = 0.0)

Reproducing a GEM-Selektor CEMDATA18 run

CEMDATA18 (Lothenbach et al., 2019) carries no ion-size parameter, so a GEM-Selektor run of a Portland cement starts from å = 0 and carries the whole non-ideality in the B-dot term, with no salting-out on the neutral species. For a KOH-dominated pore solution that is

model = HKFActivityModel(å = 0.0, Ḃ = 0.097637, Kₙ = 0.0)

which reproduces the activity coefficients such a run reports to 0.25 % on the monovalent ions and 1.2 % on the divalent ones. The package defaults are a different and more defensible model — the limiting law has no validity at I ≈ 0.2 mol/kg — and give divalent coefficients about twice as large, so the two must not be mixed in one comparison.

The A and B parameters depend on the water dielectric constant and density and can be computed explicitly via hkf_debye_huckel_params:

ab = hkf_debye_huckel_params(298.15, 1e5)   # → (A=0.5114, B=0.3288)

Valid range

The B-dot model is reliable for I ≲ 1 mol/kg. For higher ionic strengths (brines, evaporites), use the Pitzer model (planned future extension).

---

DaviesActivityModel (Davies equation)

Simpler alternative with no species-specific ionic radii (Davies, 1962). Suitable when ionic radii data are unavailable or for rapid screening calculations.

Ion activity coefficient:

log₁₀ γᵢ = −A zᵢ² (√I / (1 + √I)  −  b I)

Water activity uses the Raoult (mole fraction) approximation.

state_eq = equilibrate(state; model=DaviesActivityModel())

# Temperature-dependent A
state_eq = equilibrate(state; model=DaviesActivityModel(temperature_dependent=true))

---

Custom activity models

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

struct MyModel <: AbstractActivityModel
    # model parameters
end

function ChemistryLab.activity_model(cs::ChemicalSystem, ::MyModel)
    # Precompute species indices and constants here (called once)
    idx_solvent = only(cs.idx_solvent)
    # ...

    # Return a closure lna(n, p) -> Vector compatible with ForwardDiff
    function lna(n::AbstractVector, p)
        # p contains at minimum: p.ΔₐG⁰overRT, p.T, p.P, p.ϵ
        # n is dimensionless mole vector, same indexing as cs.species
        out = zeros(eltype(n), length(n))
        # ... fill log-activities ...
        return out
    end
    return lna
end

Pass your model to equilibrate or EquilibriumSolver:

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

A custom model should also declare its solute concentration scale, so that activity_coefficients divides by the right thing:

ChemistryLab.concentration_scale(::MyModel) = :molality

Solid solutions

Pure crystalline species have activity ln a = 0. Solid solutions are mineral phases with variable composition (e.g. C-S-H, AFm, hydrogarnet), where the activity of each end-member depends on its mole fraction within the phase.

Defining end-members and phases

End-member species must carry aggregate_state = AS_CRYSTAL. SolidSolutionPhase automatically requalifies any end-member whose class is not already SC_SSENDMEMBER, so database species with SC_COMPONENT can be passed directly.

Workflow A — pass database species directly:

using ChemistryLab

substances = build_species(datapath("cemdata18-thermofun.json"); verbose = false)
dict = Dict(symbol(s) => s for s in substances)

# SolidSolutionPhase requalifies SC_COMPONENT → SC_SSENDMEMBER automatically
ss_afm = SolidSolutionPhase("AFm",
    [dict["monosulphate12"], dict["monocarbonate"]])

Workflow B — automated via build_solid_solutions and a TOML file:

# Load all phases defined in the TOML
ss_phases = build_solid_solutions(datapath("solid_solutions.toml"), dict)

See the manual page Database Interoperability for the TOML format and the pre-built data/solid_solutions.toml file shipped with ChemistryLab.

What the shipped file is, and is not

Two of its entries reproduce CEMDATA18 (Lothenbach et al., 2019) phases of the same name. CSHQ is the six-end-member C-S-H of Kulik's downscaled solid solution model (Kulik, 2011); its KSiOH and NaSiOH members carry the uptake of potassium and sodium, whose records name Robie & Hemingway (Robie and Hemingway, 1995) among their sources, and they are what fixes the pore-solution pH of a Portland cement — leaving them out strands the alkalis in solution. C3(AF)S0.84H is the Fe-siliceous hydrogarnet.

AFm, Hydrogarnet and Hydrotalcite are deliberate alternatives to the CEMDATA18 phase model, not reproductions of it: GEM-Selektor treats monocarbonate, C3AH6, C3FH6 and hydrotalcite as pure phases, its AFm solid solution is C4AH13 + monosulphate12, and its hydrotalcite solid solution is Mg3AlC0.5OH + Mg3FeC0.5OH at Mg:Al = 3. Reproducing a published GEM-Selektor result means declaring the phases in the script, as above, rather than taking this file wholesale.

Then pass solid_solutions as a keyword to ChemicalSystem:

cs = ChemicalSystem(
    [H2O_sp, dict["monosulphate12"], dict["monocarbonate"], ...],
    ["H2O@", "Al+3", ...];           # primaries
    solid_solutions = [ss_afm],      # or solid_solutions = ss_phases
)

Activity models for solid solutions

Model	Formula	Notes
IdealSolidSolutionModel	ln aᵢ = ln xᵢ	Default, any number of end-members
RedlichKisterModel	ln aᵢ = ln xᵢ + ln γᵢ (Margules)	Binary only (2 end-members), parameters in J/mol

The solid-solution activity is computed inside the aqueous activity closure — no separate activity model is needed. The existing equilibrate(state) call handles solid solutions automatically.

Ideal solid solution

ss = SolidSolutionPhase("AFm", [em_ms, em_mc])   # IdealSolidSolutionModel() by default
cs = ChemicalSystem([...]; solid_solutions=[ss])
state_eq = equilibrate(state)

Non-ideal binary: Redlich-Kister

# Interaction parameters for monosulfoaluminate-monocarboaluminate (example values)
rk = RedlichKisterModel(a0 = 3000.0, a1 = 500.0)          # a2 defaults to 0.0
# or 3-parameter:  RedlichKisterModel(a0 = 3000.0, a1 = 500.0, a2 = 50.0)
ss = SolidSolutionPhase("AFm", [em_ms, em_mc]; model=rk)

Activity coefficients (Guggenheim / ThermoCalc convention):

Valid range

RedlichKisterModel requires exactly 2 end-members. For ternary or higher-order solid solutions, use the ideal model (IdealSolidSolutionModel).

Integration with aqueous models

Solid-solution activities are computed independently of the aqueous activity model. You can combine HKFActivityModel() for the aqueous phase with any solid-solution model — the same equilibrate call handles both.

Where to go next

The theory behind each option is in the chapter Theory: the certificate in Proving that an answer is the answer, the aqueous models in Activity models and the mixing models in Solid solutions. What the choice of model costs in numbers is measured in What the choice of activity model costs and Solid solution models, in numbers.