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. 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, as configured here, 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

No single back end certifies every case measured below, 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 met to 3 % only, 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 where it can be (the default)
state_eq = equilibrate(state; certify = false)  # one back end, polished
state_eq = equilibrate(state; certify = true)   # certified, or an error
eq, cert = equilibrate_certified(state)         # when the proof itself is wanted

Pass a solver explicitly — equilibrate(state, OptimaOptimizer()) — to use that one back end without the search. Its answer is still polished, as described next, and certificate = Ref{Any}() receives the certificate of what is returned, nothing when none could be computed.

Trace elements and the certificate

Each conservation row is judged on its own scale: the balance of row is met when   , what the row holds, so that an element present as a trace is balanced to the same relative accuracy as a major one. A phase made of such a trace, which cannot hold the amounts in moles the search compares it with, is judged on the scale of what it can hold once the search has sent it away and taken it back (OptimaSolver ≥ 0.8.7): a portlandite suspension whose aluminum is a mol seed of ettringite certifies, with the gypsum where a seed of mol puts it. Before that, no route certified it, the aluminum row ended several times out and the gypsum 14 % off at 1.74 mmol.

When no route certifies, equilibrate and equilibrate_certified return the answer with the smallest KKT error and warn, once per session, naming the row furthest out, judged as the certificate judges it. That answer is not an equilibrium. Ask for one with certify = true, which raises UncertifiedEquilibrium instead, the answer and its certificate in the exception:

eq = try
    equilibrate(state; certify = true)
catch e
    e isa UncertifiedEquilibrium || rethrow()
    @show e.certificate.balance_relative    # the worst row, relative to what it holds
    e.state                                  # the answer with the smallest KKT error
end

optimality_certificate audits any state against its budget: balance is the worst row in moles, balance_relative the worst relative to what it holds.

One equilibrium, whatever the back end

The back ends do not all solve the same equations. Ipopt minimizes the scalar n⋅μ(n), the interior point of OptimaSolver steers on the gradient μ, and the dual Newton solves the conditions of equilibrium themselves, μ(n) = −Aᵀy on the species present. The gradient of n⋅μ(n) is μ + Jᵀn, where J = ∂μ/∂n, and Jᵀn = 0 is the Gibbs–Duhem relation: where the activity model satisfies it, n⋅μ(n) is the Gibbs energy and the three agree; where it does not, the minimum of n⋅μ(n) is another composition. The ideal model, Pitzer (Pitzer, 1975), Davies on ions and the Debye–Hückel form with one ion size satisfy it; the B-dot model with its ion sizes and linear terms, Davies with a neutral solute and SIT do not. On calcite and carbon dioxide in a sodium chloride solution, Ipopt's answer meets the element balance to 1e-15 mol under every model; its dissolved calcium is the equilibrium's to 2e-6 under the ideal model, and 7e-4 off under Davies, whose dissolved carbon dioxide carries a salting-out term with no partner.

So every back end's answer is polished: the dual Newton is started from it, and the composition it certifies is the one returned. The route then decides how fast the answer is reached, not which answer it is, and a derivative taken through it is that of the composition returned. Without OptimaSolver nothing can polish, and a back end that minimizes n⋅μ(n) refuses a model that breaks the Gibbs–Duhem relation, with an error that names it, rather than return a composition that is not the equilibrium.

What the certified route does first

With every activity at one, the equilibrium is a linear program over the pure phases, minimize gᵀn subject to A n = b, n ≥ 0, and the certified route solves it before anything else.

• A budget that no non-negative amounts of the declared species can meet is refused at once, with the reason in words: the certificate has budget_feasible = false, route = :infeasible, and unplaceable, for instance "the budget asks for −0.001 mol of Ca+2, and every declared species holds it with a non-negative coefficient". There is no equilibrium to search for, and nothing is.

• Otherwise the vertex of the program, with every other species raised to the amount its multipliers give it, is the first start of the search. On cold cement pastes it certifies twenty to fifty times faster than the recipe as given; a state already near its answer loses nothing.

cert.route says which start the answer came from (:lp_start, :state, :ideal_mixing, :ideal, :continuation or :restart) and cert.n_dual_solves how many solves the search ran. lp_start = false turns the program off.

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.

Shown, not executed

The blocks in this section are shown rather than executed. 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 (Leal et al., 2017) 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 the minerals in Cemdata18 do not depend on pressure — V⁰(1 bar) = V⁰(100 bar) to the last bit for Cal — water's falls by 0.44 % over 100 bar, and the partial volumes of a few aqueous ions vary, OH⁻'s by 8 % over 100 bar. The relative lever (∂V/∂P)·P/V of water alone is 4.5e-5 at 1 bar, some 225 bar to change its volume by one percent, and a paste, whose solids do not move, has less. That is the physics, not a solver limitation: the volume of a nearly 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, a C library, does not operate on them. The equilibrium is solved once at the primal values, by the back end chosen and polished, and the sensitivities come from the conditions of equilibrium at the answer, 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

db = Dict(symbol(s) => s for s in
          build_species(datapath("slop98-inorganic-thermofun.json"); verbose = false))
cs = ChemicalSystem([db[s] for s in split("H2O@ H+ OH- CO2@ HCO3- CO3-2 Ca+2 Cal")])
i_h2o, i_cal, i_co2, i_ca =
    (findfirst(s -> symbol(s) == x, cs.species) for x in ("H2O@", "Cal", "CO2@", "Ca+2"))

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

Each mole of carbon dioxide added dissolves 0.16 mol of calcite at this point, where 3.7 mmol are dissolved. The certified route, equilibrate(state) or equilibrate_certified, returns the same derivative to fifteen digits: it lifts its answer the same way, with the active set frozen and the constraint's unknowns included, so the titrant a prescribed pH needs or the temperature an adiabatic solve reaches comes back with its derivative too. Both routes lift one level of duals at a time, so that a derivative nested in another (a gradient differentiated again for a Hessian, as an inversion needs) is exact at every level, and both lift whatever carries the duals: the amounts, the temperature, the budget, the standard potentials of the species or the parameters of the activity model. The conditions lifted are those the polished answer satisfies, so the derivative is that of the map the solve returns, which is what a derivative by finite differences of that map would approach. Without OptimaSolver, a single back end falls back to a sensitivity solved in Float64, for the amounts of the state only and one level deep, and refuses the rest rather than return derivatives that are not there.

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; the finite-difference reference gives 2×10⁻⁹, the size of its truncation error.

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	any start, default
Val(:log)	log nᵢ	refining a composition already close to the answer

In logarithmic variables, the optimality condition of species i is its amount nᵢ times the gap between its potential and the combination of component potentials its formula implies, so a species held at the regularization floor contributes almost nothing to the optimality test. A starting state in which most species sit at that floor, such as one built from a recipe, then already passes the test of the interior point: from it, an explicit solver in the logarithmic space returns the start or stops without converging. Called without a solver, equilibrate goes through the certified route, whose answer does not depend on this choice; with an explicit solver, the logarithmic space is to be started from a solved state:

state_eq = equilibrate(state, IpoptOptimizer())                        # linear, from the recipe
state_eq_log = equilibrate(state_eq, IpoptOptimizer(); 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

Solver options are passed as keyword arguments and forwarded to the optimizer's solve. Ipopt reads reltol as its convergence tolerance tol; it has no counterpart for abstol, which it ignores:

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

OptimaSolver, the default back end once it is loaded, ignores these keywords: its tolerances are the options of the optimizer, OptimaOptimizer(; tol = 1e-12). The answer of equilibrate_certified is judged by its certificate, whatever the back end's tolerance.

Activity models

All activity models inherit from AbstractActivityModel. The three in the table below cover ideal behavior through to the extended Debye-Hückel level used by standard geochemical codes; PitzerActivityModel and SITActivityModel extend the range to concentrated solutions, each with its own parameter set.

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)

activity_model_range returns the bound of this table for a model, and nothing where the table gives no number. None of these formulas announces that it has left its range, so the comparison has to be made: check ionic_strength(eq) <= activity_model_range(model) on any answer that matters. When equilibrate_certified gives up on an answer whose ionic strength is past the bound, its message says so. It is not a formality: measured on blended cement pastes taken to full reaction, the limiting law with a B-dot term (å = 0) did not certify, and with an ion size per ion, HKFActivityModel(), the same budgets did.

---

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ᵢ + ln(P/P°)

state_eq = equilibrate(state)   # DiluteSolutionModel is the default

---

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

Implements the extended Debye-Hückel model of Helgeson (1969) and 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 — Table 3 of Helgeson et al. (1981) (27 common ions).

4. REJ_CHARGE_DEFAULT — fallback by formal charge, (Xu et al., 2012) Table H.1-1.

5. model.å_default (default: the ion size of NaCl, 3.72 Å, (Helgeson et al., 1981) Table 2).

å_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 (Kulik et al., 2013), 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, that of NaCl, overriding the per-species tables
model = HKFActivityModel(å = literature_value("Helgeson1981", "nacl_distance_of_closest_approach"))

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

The model Cemdata18 prescribes

Cemdata18 was derived with, and prescribes, an extended Debye–Hückel equation with one ion-size parameter and one B-dot common to every ion, set by the dominant electrolyte of the solution (Lothenbach et al., 2019; Appendix C, Eq. C.1): 3.67 Å and 0.123 kg/mol for KOH, 3.31 Å and 0.098 for NaOH, with the same B-dot on the neutral species. It is cemdata18_activity_model:

model = cemdata18_activity_model(:KOH)    # or :NaOH

Lothenbach et al. (2019) state it applicable up to about 1 mol/kg of ionic strength, which activity_model_range returns, and which the certificate of equilibrate_certified compares with the ionic strength of each answer (cert.within_activity_range).

Reproducing one particular GEM-Selektor run

A GEM-Selektor run of a Portland cement with CEMDATA18, whose printed activity coefficients are the fixture of A CEM I at equilibrium, with every solid solution declared, was configured differently: fitting its coefficients by charge class identifies the Debye–Hückel limiting law, å = 0, with a B-dot of about 0.0976 and none on the neutral species. That is a configuration of that run, not what Cemdata18 prescribes, and the two models give different coefficients at the ionic strength of that run, 0.2097 mol/kg:

using ChemistryLab
I = 0.2097                                   # mol/kg, the ionic strength of that run
log10γ(m, z) = -m.A * z^2 * sqrt(I) / (1 + m.B * something(m.å, 0.0) * sqrt(I)) + m.Ḃ * I
prescribed = cemdata18_activity_model(:KOH; temperature_dependent = false)
that_run = HKFActivityModel(å = 0.0, Ḃ = 0.0976, Kₙ = 0.0)
[(model = name, γ1 = round(10^log10γ(m, 1); digits = 3), γ2 = round(10^log10γ(m, 2); digits = 3))
 for (name, m) in (("Cemdata18, Eq. C.1", prescribed), ("that GEMS run", that_run))]

2-element Vector{@NamedTuple{model::String, γ1::Float64, γ2::Float64}}:
 (model = "Cemdata18, Eq. C.1", γ1 = 0.75, γ2 = 0.265)
 (model = "that GEMS run", γ1 = 0.611, γ2 = 0.121)

The run printed γ = 0.611 for the monovalent ions and 0.121 for the divalent ones, which the second line reproduces and the first does not. To compare with that run, use its model; to follow Cemdata18, use the first:

model = HKFActivityModel(å = 0.0, Ḃ = Ḃ_gems, Kₙ = 0.0)   # that run
model = cemdata18_activity_model(:KOH)                   # the paper

The package default, HKFActivityModel(), is a third model, with an ion size per ion: do not mix it with either 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 PitzerActivityModel.

---

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. activity_coefficients reads the coefficients of a model without a formula per ion from its closure, as    with the molalities built on the closure's own floor, and it does so only on the molality scale:

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_hg = SolidSolutionPhase("C3(AF)S0.84H",
    [dict["C3AFS0.84H4.32"], dict["C3FS0.84H4.32"]])

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 extensions, filters and solid solutions 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 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_SO4_OH and AFt_SO4_CO3 carry CEMDATA18's non-ideal parameters, read from data/literature/Lothenbach2019.json, with two instances each because both unmix. Hydrogarnet and Hydrotalcite are deliberate alternatives to the CEMDATA18 phase model, not reproductions of it: GEM-Selektor treats C3AH6, C3FH6 and hydrotalcite as pure phases, 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["C3AFS0.84H4.32"], dict["C3FS0.84H4.32"], ...],
    ["H2O@", "Al+3", ...];           # primaries
    solid_solutions = [ss_hg],       # 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
RegularSolutionModel	ln γᵢ from one Wᵢⱼ per pair	Any number of end-members; convexity from mixing_convexity
SubregularSolutionModel	ln γᵢ from Wᵢⱼ ≠ Wⱼᵢ per pair and Wᵢⱼₖ per triple	Any number of end-members, asymmetric, after Helffrich and Wood (1989)
MulticomponentRedlichKisterModel	a series in xᵢ − xⱼ per pair, a term per triple	Any number of end-members, any order, after Redlich and Kister (1948)
VanLaarModel	one Wᵢⱼ per pair, one size αᵢ per end-member	Any number of end-members, asymmetric through the sizes, after Holland and Powell (2003)
SublatticeModel	ln aᵢ = Σₛ mₛ ln y_{s,σₛ(i)}	Ideal mixing on sites; CNASH_ss of Myers et al. (2014), CSH3T of Kulik (2011)
CompoundEnergyModel	ideal site mixing, plus the reference surface of the end-member energies and regular interactions on a site	The CASH+ C-S-H of Kulik et al. (2022); needs the end-members' standard Gibbs energies

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("C3(AF)S0.84H", [em_al, em_fe])   # IdealSolidSolutionModel() by default
cs = ChemicalSystem([...]; solid_solutions=[ss])
state_eq = equilibrate(state)

Non-ideal binary: Redlich-Kister

# Cemdata18's AFm SO4/OH binary, from its dimensionless Guggenheim parameters
# (data/literature/Lothenbach2019.json): a = α R T, C4AH13 first.
p  = literature_row("Lothenbach2019", "guggenheim_parameters", "AFm SO4/OH")
RT = R_GAS * 298.15
rk = RedlichKisterModel(a0 = p.alpha0 * RT, a1 = p.alpha1 * RT)   # a2 defaults to 0.0
# It unmixes, so it needs two coexisting compositions:
ss = SolidSolutionPhase("AFm_SO4_OH", [dict["C4AH13"], dict["monosulphate12"]];
                        model = rk, instances = 2)

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), the symmetric multi-component Margules model (RegularSolutionModel), or an asymmetric one (SubregularSolutionModel, MulticomponentRedlichKisterModel, VanLaarModel).

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.