ChemistryLab.jl

Appearance

Equilibrium ​

Activity models ​

ChemistryLab.REJ_CHARGE_DEFAULT Constant
julia
REJ_CHARGE_DEFAULT::Dict{Int,Float64}

Fallback effective electrostatic radii åᵢ [Å] indexed by formal charge, from ToughReact V2 (Xu et al. 2011, Table A2; after Helgeson et al. 1981).

Used by HKFActivityModel with priority 3 in the radius lookup chain: sp[:å] > REJ_HKF > REJ_CHARGE_DEFAULT > model.å_default.

See also: REJ_HKF, HKFActivityModel.

ChemistryLab.REJ_HKF Constant
julia
REJ_HKF::Dict{String,Float64}

Effective electrostatic radii åᵢ [Å] for aqueous ions from Helgeson, Kirkham & Flowers (1981), Am. J. Sci. 281, Table 3.

Keys are PHREEQC-format formula strings (e.g. "Na+", "Ca+2", "SO4-2"). Used by HKFActivityModel with priority 2 in the radius lookup chain: sp[:å] > REJ_HKF > REJ_CHARGE_DEFAULT > model.å_default.

See also: REJ_CHARGE_DEFAULT, HKFActivityModel.

ChemistryLab.AbstractActivityModel Type
julia
abstract type AbstractActivityModel end

Base type for all activity models: the correction that turns a concentration into a chemical potential.

What an activity model is for

The chemical potential of a species is

and the whole of the equilibrium calculation — build_potentials, the Gibbs minimization, the certificate — depends on nothing else about the solution. So an activity model is the single place where a real solution stops being ideal, and it has to answer two questions, not one:

  1. the solute activity, a_i = γ_i m_i (or γ_i c_i), where γ_i corrects for the fact that an ion in a solution of other ions is not in the same state as one alone at the same concentration;

  2. the solvent activity a_w, which is not a small correction on the way to the same answer. It is what a hydrate reaction consumes, and in a cement paste with little mixing water it is the quantity that decides how far the reaction goes.

The second is the one to look at when choosing a model here, because the three built-in models treat it in three different ways and only one of them is defensible past a dilute solution. See Activity models for the comparison.

The interface

A concrete subtype must implement

julia
activity_model(cs::ChemicalSystem, model::YourModel) -> lna    # lna(n, p)
concentration_scale(model::YourModel)                          # :molality | :molarity

lna(n, p) receives the full mole vector n, indexed like cs.species, and the parameter tuple p carrying at least ϵ and usually T, P and ΔₐG⁰overRT; it returns ln aᵢ for every species — solutes, solvent, pure crystals (0), gases, and solid-solution end-members. Three properties are required rather than nice to have:

  • it is differentiated by ForwardDiff at every Newton step, so the output element type must follow n and any regularization must be smooth;

  • it must call _solid_solution_lna!, or the end-members of a solid solution silently get ln a = 0;

  • concentration_scale has no fallback: without it the aqueous accessors raise a MethodError rather than guessing a convention.

One discipline is worth knowing before adding a model: the scalar formula lives in _log10γ_ion and _log10γ_neutral, which the closure calls and activity_coefficients calls, so a model states its algebra once and the accessor cannot drift from what the solver used. A model whose coefficients are not a per-species function of (z, å, I) — an ion-interaction model, where the sum runs over pairs — does not fit that pair of helpers and supplies its own vector-valued path instead.

See also: DiluteSolutionModel, HKFActivityModel, DaviesActivityModel, concentration_scale.

ChemistryLab.DaviesActivityModel Type
julia
struct DaviesActivityModel{T<:Real} <: AbstractActivityModel

The Davies equation ((Davies, 1962)): Debye-Hückel with the ion size taken out, so that nothing per-species has to be known.

Formulas

\ln a_i = \ln 10 \cdot \log_{10}\gamma_i + \ln m_i on the molality scale, and for the solvent Raoult's law, \ln a_w = \ln x_w.

The physics, and the one thing to know before using it

Davies replaces 1 + B å √I with 1 + √I, which amounts to fixing one ion size for everything (about 3 Å at 25 °C), and adds b I to bend the curve back up. The gain is that no per-species datum is needed, which is why it is the model to reach for when a species list contains ions no radius table covers. The loss is that every ion of the same charge is now identical, so the model cannot distinguish Na⁺ from K⁺ at all.

Its water activity does not come from its own activity coefficients

The solutes are non-ideal and the solvent is treated as ideal: a_w = x_w counts molecules. So this model does not satisfy the Gibbs-Duhem relation between its own γᵢ and its a_w — the two halves are not derived from one excess Gibbs energy. In a dilute solution the error is negligible, because both are near their ideal values anyway. In a hydrating cement paste, where the water activity is what decides how far the reaction goes, the two halves disagree and the answer inherits the disagreement.

Use HKFActivityModel, whose a_w comes from an osmotic coefficient, whenever the water activity is part of the question.

Inputs, their defaults, and where each default comes from

fielddefaultunitprovenance
A0.5114(kg/mol)^½(Helgeson et al., 1981) Table 1 at 25 °C / 1 bar, and reproduced by hkf_debye_huckel_params from this package's water model
b0.3kg/molpart of the published equation — Davies fixed it, it is not a free parameter of this implementation
bₙ0.1kg/mola generic salting-out coefficient for neutral species. No source recorded in this package
temperature_dependentfalse—recompute A from p.T, p.P at every call

Valid range

I ≲ 0.5 mol/kg for the activity coefficients — Davies is usually quoted as useful to 0.1 and tolerable to 0.5 — and I ≲ 0.1 mol/kg for anything that depends on the water activity.

Examples

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

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

See also: HKFActivityModel, DiluteSolutionModel, Activity models.

ChemistryLab.DaviesActivityModel Method
julia
DaviesActivityModel(; A=0.5114, b=0.3, bₙ=0.1, temperature_dependent=false)

Construct a DaviesActivityModel.

ChemistryLab.DiluteSolutionModel Type
julia
struct DiluteSolutionModel <: AbstractActivityModel

The ideal dilute solution: every activity coefficient is exactly 1, and an activity is a concentration.

phaselawexpression
solventRaoultln a = ln x_w, the mole fraction
aqueous solutesHenryln a = ln(cᵢ/c°), c° = 1 mol/L
pure crystals—ln a = 0
gasideal mixtureln a = ln xᵢ
solid-solution end-membersideal mixingln a = ln xᵢ within the phase

The physics, and where it runs out

Ideality means an ion does not notice the others. That is true in the limit of infinite dilution and stops being true as soon as the ionic cloud around an ion has a measurable effect on its energy — that is, from a few millimolal upwards for a charged species. Take I ≲ 0.01 mol/kg as the range in which the answer is the answer, and treat anything above as a screening calculation.

Its water activity is a mole fraction, and that is the real limit

a_w = x_w is Raoult's law, which counts molecules and knows nothing about what they are. A cement pore solution at x_w = 0.99 gets a_w = 0.99 whatever it holds in solution, and a paste short of mixing water gets an a_w that follows the amount of water and not its state. If the water activity matters for what is being computed — and in a hydrating paste it decides how far the reaction goes — this model is the wrong one, and HKFActivityModel is the least that will do.

It is nonetheless the default, for two reasons that are about the solve and not about the chemistry: it is exact in the dilute limit, and it makes the log-activity linear in ln n, which is the best-conditioned objective the minimizer will ever see. Start here, then change the model and see whether the answer moves.

See also: HKFActivityModel, DaviesActivityModel, Activity models.

ChemistryLab.HKFActivityModel Type
julia
struct HKFActivityModel{T<:Real} <: AbstractActivityModel

The extended Debye-Hückel model with a B-dot term — the workhorse of aqueous geochemistry, and the model PHREEQC and EQ3/6 use ((Helgeson, 1969), (Helgeson et al., 1981), (Parkhurst and Appelo, 2013)).

Formulas

Ionic species (z ≠ 0), on the molality scale with m° = 1 mol/kg:

Neutral aqueous species (z = 0), the Setschenow form:

and in both cases \ln a_i = \ln 10 \cdot \log_{10}\gamma_i + \ln m_i.

The water activity comes from the osmotic coefficient φ, so that it is consistent with the same parameters rather than being a separate assumption:

with σ(x) = (3/x³)(x − 2\ln(1+x) − 1/(1+x) + 1) (_hkf_sigma, Helgeson et al. 1981 Eqs. 132–137).

The physics each term carries

  • −A z² √I — an ion polarizes the solution around it, and the resulting ionic cloud screens its charge and lowers its energy. Hence γ < 1: an ion in an electrolyte is more stable than one alone at the same molality. The √I dependence is not fitted, it is the Debye-Hückel limiting law, exact as I → 0, and A is fixed by the dielectric constant and density of water.

  • 1 + B å √I — the cloud cannot approach closer than the ion's own size, which cuts the screening off. This is the one term that carries something about the identity of the ion (through å), and it is what extends a law valid at millimolal to roughly 1 mol/kg.

  • + Ḃ I — empirical, linear, positive, and the reason γ turns back upwards at high ionic strength. It is not a physical term that was derived and then measured: (Anderson and Crerar, 1993) (§17.7.1, pp. 445-446) describe it as a deviation function, defined by Helgeson as the difference between the observed activity coefficient of an electrolyte — NaCl — and what the Debye-Hückel expression predicts for it. So it carries short-range ion-solvent and ion-ion interaction and whatever the first two terms failed to capture, together, in one number fitted to one salt. That is why this model has a ceiling rather than an asymptote, and why the ceiling is somewhere around a molal rather than at a sharp value.

  • K_n I for a neutral species — salting out. Water engaged around ions is water unavailable to solvate a neutral molecule, so its activity rises with I and its solubility falls. CO₂(aq) is the case that matters here.

Inputs, their defaults, and where each default comes from

fielddefaultunitprovenance
A0.5114(kg/mol)^½(Helgeson et al., 1981) Table 1 at 25 °C / 1 bar — and reproduced to 1e-3 by hkf_debye_huckel_params from this package's own water model, which is the check in test/activities.jl
B0.3288Å⁻¹(kg/mol)^½same
Ḃ0.041kg/molthe value conventionally carried for a NaCl-dominated solution at 25 °C. What the term is has a source ((Anderson and Crerar, 1993) §17.7.1); this particular number has none recorded in this package, so treat it as a convention rather than a measurement
Kₙ0.1kg/mola generic salting-out coefficient, no source recorded; overridden per species by sp[:Kₙ], which is how CO₂(aq) gets its own
å_default3.72Ålast resort, reached only for a charge no table covers (`
ånothingÅone common radius for every ion, overriding the tables. å = 0 collapses the denominator and gives the limiting law plus Ḃ I
temperature_dependentfalse—recompute A and B from p.T, p.P at every call (needs T and P in p)

Per-ion radii come from REJ_HKF ((Helgeson et al., 1981) Table 3) and, failing that, from REJ_CHARGE_DEFAULT ((Xu et al., 2011) Table A2).

Three of these defaults are conventions, not data

Ḃ, Kₙ and å_default are numbers this package carries without a source to point at. They are in the range everyone uses and they are almost certainly right, but the honest statement is that they have not been traced, and each is a keyword away from being replaced. A, B and the radius tables are traceable, and A and B are additionally derived here rather than tabulated.

What is approximated in the water activity

The γᵢ use a per-ion radius; the osmotic coefficient uses one charge-weighted mean radius å_eff = Σ mᵢzᵢ²åᵢ / Σ mᵢzᵢ². So a_w is not exactly the Gibbs-Duhem integral of the γᵢ this model returns when the ions differ in size, and the inconsistency is measurable: the Gibbs-Duhem residual on 0.3 mol/kg NaCl is a few parts in a thousand (test/activities.jl asserts < 5e-3), where an exactly consistent model would sit at solver tolerance. Å-level differences between Na⁺ (1.91) and Cl⁻ (1.81) are enough to produce it.

That is acceptable for a pore solution at 0.1–0.5 mol/kg and it is the first thing that breaks as the solution concentrates. It is also structural: no choice of Ḃ repairs it, because the defect is in using a mean radius at all.

Ionic radius lookup

The effective radius åᵢ is resolved in order:

  1. model.å — a common radius for every ion, when given. Short-circuits the rest of the chain, so a per-species table entry cannot silently override it.

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

  3. REJ_HKF — Helgeson et al. (1981) Table 3, keyed by PHREEQC formula.

  4. REJ_CHARGE_DEFAULT — fallback by formal charge.

  5. model.å_default — reached only for a charge no table covers, i.e. |z| ≥ 5. It is not a way to impose a common radius; pass å for that.

Valid range

I ≲ 1 mol/kg for the activity coefficients, and less than that for the water activity, for the reason above. Beyond a few mol/kg the Ḃ I term is doing work it was never fitted to do, and nothing in the formula announces it: the model goes on returning finite, plausible numbers. A cement pore solution normally sits at 0.1–0.5 mol/kg and is comfortably inside; a paste short of mixing water is not — see solvent_fraction and Activity models.

Examples

julia
# Default model at 25 °C / 1 bar (fixed A, B)
model = HKFActivityModel()

# Temperature-dependent A and B (recomputed at each solve)
model_tdep = HKFActivityModel(temperature_dependent=true)

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

See also: DaviesActivityModel, hkf_debye_huckel_params, activity_coefficients, ionic_strength.

ChemistryLab.HKFActivityModel Method
julia
HKFActivityModel(; A=0.5114, B=0.3288, Ḃ=0.041, Kₙ=0.1, å_default=3.72,
                   å=nothing, temperature_dependent=false) -> HKFActivityModel

Construct an HKFActivityModel with the given parameters.

Default values are from Helgeson et al. (1981), Table 1, at 25 °C / 1 bar.

å imposes one common effective radius on every charged aqueous species, overriding the per-species tables. Use it to reproduce a published model that was run with a single ion-size parameter — which is what GEM-Selektor, PHREEQC's -gamma and most cement models do. Note that å_default does not do this: it is only the last resort of the lookup chain and is reached only for charges no table covers. å = 0 gives the Debye-Hückel limiting law plus the B-dot term.

Examples

julia
# The package default: per-species radii from REJ_HKF, EQ3/6 NaCl B-dot.
HKFActivityModel()

# One common radius of 3.72 Å for every ion.
HKFActivityModel(å = 3.72)

# The Debye-Hückel limiting law with a KOH-background B-dot, which is what a
# GEM-Selektor CEMDATA18 run of a Portland cement uses: CEMDATA18 carries no
# ion-size parameter, so GEMS starts from å = 0, and the B-dot term is not
# applied to neutral species.
HKFActivityModel(å = 0.0, Ḃ = 0.097637, Kₙ = 0.0)
ChemistryLab._excess_ln_gamma Method
julia
_excess_ln_gamma(model, k, x, T) -> Real

Return the excess log-activity coefficient ln γₖ for end-member k (1-based index) of a solid solution with mole-fraction vector x at temperature T (K).

AD-compatible: all branches preserve ForwardDiff.Dual through computations on x.

Methods:

ChemistryLab._hkf_sigma Method
julia
_hkf_sigma(x) -> Real

Compute the σ function used in the osmotic coefficient formula (Helgeson et al. 1981, Eq. 132–137):

julia
σ(x) = (3/x³)(x − 2 ln(1+x) − 1/(1+x) + 1)

For |x| < 1e-3, a Taylor series 1 − (3/2)x + (9/5)x² is used to avoid catastrophic cancellation. Both branches agree to O(x³) at the threshold, so the gradient is continuous.

AD-compatible: branching is on ForwardDiff.value(x), not on the Dual itself.

ChemistryLab._solid_solution_lna! Method
julia
_solid_solution_lna!(out, _n, ss_groups, ss_models, T, ϵ)

Fill out[i] with ln aᵢ = ln xᵢ + ln γᵢ for all solid-solution end-members.

ss_groups[k] and ss_models[k] describe the k-th solid-solution phase. T is the temperature in K (only relevant for non-ideal models). ϵ is a regularization floor to avoid log(0).

ForwardDiff-compatible.

ChemistryLab.activity_model Method
julia
activity_model(cs::ChemicalSystem, model::DaviesActivityModel) -> Function

Return a closure lna(n, p) -> Vector computing log-activities for the Davies (1962) model. No species-specific ionic radii are required.

AD-compatible: all closure computations accept ForwardDiff.Dual inputs.

ChemistryLab.activity_model Method
julia
activity_model(cs::ChemicalSystem, ::DiluteSolutionModel) -> Function

Return a closure lna(n, p) -> Vector{Float64} computing the vector of log-activities for the dilute ideal solution model.

The returned function has signature lna(n, p) where:

  • n: dimensionless mole vector (same indexing as cs.species)

  • p: NamedTuple containing at least ϵ (floor value to avoid log(0))

Solid-solution end-members (class SC_SSENDMEMBER) receive ln aᵢ = ln xᵢ where xᵢ = nᵢ / Σnⱼ within the same solid-solution phase.

All quantities are dimensionless — units are stripped at construction time.

ChemistryLab.activity_model Method
julia
activity_model(cs::ChemicalSystem, model::HKFActivityModel) -> Function

Return a closure lna(n, p) -> Vector computing log-activities for the extended Debye-Hückel (B-dot) model of Helgeson (1969).

The closure captures all species indices and ionic radii at construction time. Inside lna:

  • Solutes: molality convention, B-dot formula for ions, salting-out for neutrals.

  • Solvent: osmotic coefficient from Gibbs-Duhem (σ-function).

  • Crystals: ln a = 0 (pure solid).

  • Gas: ideal mixture ln a = ln(xᵢ).

If model.temperature_dependent=true, p must contain T (K) and P (Pa) — both are provided automatically by _build_params.

AD-compatible: all closure computations accept ForwardDiff.Dual inputs.

ChemistryLab.build_potentials Method
julia
build_potentials(cs::ChemicalSystem, model::AbstractActivityModel) -> Function

Return a closure μ(n, p) -> Vector{Float64} computing dimensionless chemical potentials μ_i / RT for all species.

  

The returned function is compatible with SciML solvers:

  • n: dimensionless mole vector

  • p: NamedTuple containing:

    • ΔₐG⁰overRT: vector of standard Gibbs energies of formation divided by RT

    • ϵ: regularization floor (e.g. 1e-30)

All quantities are dimensionless — caller is responsible for stripping units from ΔₐG⁰overRT before passing them in p.

Examples

julia
julia> cs = ChemicalSystem([
           Species("H2O"; aggregate_state=AS_AQUEOUS, class=SC_AQSOLVENT),
           Species("Na+"; aggregate_state=AS_AQUEOUS, class=SC_AQSOLUTE),
       ]);

julia> μ = build_potentials(cs, DiluteSolutionModel());

julia> n = [55.5, 0.1];

julia> p = (ΔₐG⁰overRT = [-95.6, -105.6], ϵ = 1e-30);

julia> length(μ(n, p)) == 2
true
ChemistryLab.concentration_scale Method
julia
concentration_scale(model::AbstractActivityModel) -> Symbol

The concentration scale of the model's solute standard state, :molality or :molarity.

An activity coefficient is only defined relative to a scale: a = γ m on the molality scale, a = γ c / c° on the molarity scale. Nothing in the numbers says which one a given model used — DiluteSolutionModel is on the molarity scale but takes ρ = 1 kg/L, so its activities coincide numerically with molalities — so the scale has to be asked rather than inferred. It stops being a formality as soon as the density departs from 1 kg/L, and it is the origin of the 0.0013 pH offset that model carries.

See also: activity_coefficients, molalities.

ChemistryLab.hkf_debye_huckel_params Method
julia
hkf_debye_huckel_params(T_K, P_Pa) -> NamedTuple{(:A, :B)}

Compute the Debye-Hückel A and B parameters from the water density ρ [g/cm³] and dielectric constant εᵣ at temperature T_K (K) and pressure P_Pa (Pa).

Formulas (Helgeson et al. 1981):

julia
A(T,P) = 1.824829238×10⁶ × √ρ / (εᵣ T)^(3/2)    [(kg/mol)^(1/2)]
B(T,P) = 50.29158649      × √ρ / √(εᵣ T)          [Å⁻¹ (kg/mol)^(1/2)]

where ρ is in g/cm³.

Uses water_thermo_props (HGK equation of state) and water_electro_props_jn (Johnson-Norton dielectric constant).

AD-compatible (ForwardDiff-safe). Returns a NamedTuple (A=..., B=...).

Examples

julia
julia> p = hkf_debye_huckel_params(298.15, 1e5);

julia> isapprox(p.A, 0.5114; rtol=1e-3)
true

julia> isapprox(p.B, 0.3288; rtol=1e-3)
true

Ion-interaction (Pitzer) model ​

The virial expansion of the excess Gibbs energy: a coefficient per ion pair and per triplet, so γ and the osmotic coefficient come from one function and the Gibbs-Duhem relation holds by construction. See Activity models §6 for the derivation and The Pitzer model for what the shipped parameter set can be used with. Nothing here carries a default: the parameters are caller input.

ChemistryLab.PitzerActivityModel Type
julia
PitzerActivityModel(; parameters, temperature_dependent = false)

The Pitzer ion-interaction activity model.

parameters is a PitzerParameters and has no default: this model cannot be reached without one, which is the point. Whether the set is complete is not a property of the set alone but of the set relative to a species list, so the check happens when the model meets a ChemicalSystem: every cation-anion pair the system contains must have a beta0 entry, and the error names the pairs that do not.

Why this model rather than an extended Debye-Hückel one

γ and the osmotic coefficient are derived from one excess Gibbs energy — a virial expansion in the molalities, with one coefficient per ion pair and one per triplet — so the Gibbs-Duhem relation between the solutes and the solvent holds by construction rather than approximately. That is what DaviesActivityModel does not do and HKFActivityModel does only up to a mean-radius approximation. See Activity models.

What is not implemented

  • the higher-order electrostatic terms and , which correct the interaction of two ions of unsymmetrical charge (a 1+ with a 2+, say). They are exactly zero for a symmetrical pair, so a single 1-1 or 2-2 electrolyte is unaffected; in a mixture of Na⁺ and Ca²⁺ they are a real omission, and the docstring says so rather than the code pretending otherwise;

  • any temperature dependence of the interaction parameters themselves. A published set is fitted at one temperature — 25 °C for the set shipped here — and this model does not extrapolate it. temperature_dependent governs only , which comes from the water model.

Example

julia
p = build_pitzer_parameters(datapath("pitzer-reardon1990.toml"))
model = PitzerActivityModel(; parameters = p)

See also: PitzerParameters, build_pitzer_parameters, pitzer_origin.

ChemistryLab.PitzerParameters Type
julia
PitzerParameters(; beta0, beta1, beta2, Cphi, theta, psi, lambda,
                   alpha1_22 = 1.4, alpha1 = 2.0, alpha2 = 12.0, b = 1.2,
                   origin = Dict{Tuple{String, String}, String}())

The interaction parameters of a Pitzer model, keyed by species symbol.

Every keyword is mandatory, origin and the four shape constants excepted: a Pitzer model is only as good as the set it is given, and a default would be a number invented on the caller's behalf. Omitting one raises UndefKeywordError before anything is computed.

The tables

keywordkeysmeaning
beta0, beta1, beta2(cation, anion)the second virial coefficient of that pair, as three terms of one ionic-strength dependence
Cphi(cation, anion)the third virial coefficient of that pair
theta(ion, ion) of like chargeinteraction between two cations, or two anions
psi(ion, ion, ion)a triplet: two anions with a cation, or two cations with an anion
lambda(neutral, ion)a neutral solute with an ion — salting out, in Pitzer's form

theta, psi and lambda entries are looked up in any order of their keys. A missing theta, psi or lambda is taken as zero, which is the convention of the literature these tables come from; a missing beta0 for a pair the system actually contains is an error, raised when the model is attached to a ChemicalSystem — see PitzerActivityModel.

The shape constants

alpha1 = 2.0, alpha1_22 = 1.4 (used when both ions carry a charge of magnitude 2 or more), alpha2 = 12.0 and b = 1.2 kg^½·mol^−½ are not fitted parameters: they fix the functional form, and the published β tables were fitted assuming them. Changing one invalidates the table it is used with. They are keywords so that a set fitted with a different convention can be used, not so that they can be tuned.

Provenance

origin maps a pair to a free-text note, and the loader build_pitzer_parameters fills it from the data file. It exists because a published set can contain values that were estimated by analogy rather than measured, and a caller reading beta0 has otherwise no way of telling the two apart. pitzer_origin reads it back.

See also: PitzerActivityModel, build_pitzer_parameters.

ChemistryLab.activity_model Method
julia
activity_model(cs::ChemicalSystem, model::PitzerActivityModel) -> Function

Return the closure lna(n, p) of the Pitzer model for cs.

The completeness of model.parameters is checked here, against the species cs actually contains, and a missing cation-anion pair raises rather than defaulting to ideal behavior.

ChemistryLab.pitzer_origin Method
julia
pitzer_origin(p::PitzerParameters, cation, anion) -> String

What the parameters of that pair are: "fitted", "estimated:<analog>", or "unrecorded" when the set carries no note.

A published Pitzer set can contain values obtained by analogy with a chemically similar ion rather than by fitting a measurement — Reardon's cement set does so for every silicate, aluminate and ferrate pair, which are exactly the ions a cement assemblage needs. This accessor is how that shows up in a calculation instead of staying in a paper's footnote.

ChemistryLab.build_pitzer_parameters Method
julia
build_pitzer_parameters(toml_file) -> PitzerParameters

Read a Pitzer interaction-parameter set from a TOML file.

Nothing about this is automatic: the file has to be named, which is the whole point. data/pitzer-reardon1990.toml ships with the package and is reached as

julia
p = build_pitzer_parameters(datapath("pitzer-reardon1990.toml"))

File format

toml
[meta]
name = "..."
temperature_K = 298.15
speciation = "dissociated"        # informational; read by the caller, not here

[[binary]]                        # one per cation-anion pair
cation = "Na+"
anion  = "Cl-"
beta0  = 0.0765
beta1  = 0.2664
beta2  = 0.0
Cphi   = 0.0013
origin = "fitted"                 # or "estimated:HSO4-"

[[theta]]                         # like-charge pair, either order
i = "Na+"
j = "K+"
value = -0.012

[[psi]]                           # triplet, any order
i = "Na+"
j = "Cl-"
k = "SO4-2"
value = 0.0014

[[lambda]]                        # neutral with an ion
neutral = "H2CO3@"
ion = "Na+"
value = 0.1
origin = "fitted"

beta1, beta2 and Cphi default to zero within a [[binary]] entry; beta0 does not, because an entry without it describes nothing. The shape constants alpha1, alpha1_22, alpha2 and b may be given under [meta] and otherwise take the conventional values documented in PitzerParameters.

What it refuses

A [[binary]] entry without cation, anion or beta0, and a duplicate pair — silently keeping the last of two conflicting entries would be worse than stopping. Species names are not checked against any database here: whether the set covers the system is decided when the model is attached to it, which is where the error can name the pairs that are missing.

See also: PitzerParameters, PitzerActivityModel, pitzer_origin.

Aqueous properties ​

Molalities, ionic strength, activity coefficients and the activity-convention pH, read back off a solved state. See Reading the aqueous properties back for the two traps these exist to avoid.

ChemistryLab.ELECTRON Constant
julia
ELECTRON :: Species

The electron, e⁻, carrying the conventional standard state      .

This is a convention, exactly as it is for H⁺, and not a measurement: no thermodynamic database tabulates a free electron in solution. Fixing its standard state at zero is what makes a half-reaction's well defined, and every redox potential computed from it inherits that convention. The choice is the usual one in geochemistry, so the numbers here are comparable with published half-reaction constants.

Used to balance half-reactions:

julia
julia> r = Reaction([byname["SO4-2"], byname["H+"], ELECTRON,
                     byname["HS-"], byname["H2O@"]]);

julia> r.equation
"SO₄²⁻ + 9H⁺ + 8e⁻ = 4H₂O@ + HS⁻"

See also: pe, Eh, FixedpE.

ChemistryLab.SOLVENT_FRACTION_FLOOR Constant
julia
SOLVENT_FRACTION_FLOOR

The mole fraction of solvent below which an aqueous phase is no longer a solution, and every quantity derived from it is meaningless. See solvent_fraction.

0.5 is not a modeling choice but a floor no real solution comes near: it is about 28 mol of solute per kilogram of water, past saturation for anything. A value below it means the solve has driven the water into the solids, not that the solution is concentrated.

ChemistryLab.Eh Method
julia
Eh(state, model; couple = "SO4-2" => "HS-") -> Float64

The redox potential in volts, from the same half-reaction as pe through the Nernst relation

with the Faraday constant. At 25 °C the factor is 0.05916 V per pe unit.

The caveat of pe applies unchanged: this is the potential of the couple named, and a paste whose couples are not at mutual equilibrium has no single Eh.

ChemistryLab._homotopy_rung Method
julia
_homotopy_rung(cs, A, i_w, n0, model, λ, start, ϵ, verbose, atol, rtol)
    -> Union{ChemicalState, Nothing}

Solve one rung of the continuation, and accept it only if it conserves mass.

nothing means "refuse this rung", which the walk answers by taking a smaller step. Two things can go wrong at a rung, and only one of them raises: a back end can throw, or it can return an answer that is not on the constraint surface. The second is the dangerous one, because the walk would then carry that composition forward as the start of every later rung. Measured on a CEM I paste: an accepted rung off the balance by moles takes the whole walk with it, and equilibrate_certified ends on an answer with an element balance of 6.7 mol — every hydrate at zero, and a table of amounts that reads like a result.

The test is the mixed one, |rᵢ| ≤ atol + rtol · scaleᵢ, and it has to be mixed. A purely relative test is meaningless on this problem because two conservation rows carry a legitimately negligible budget: electroneutrality is exactly zero, and a cement recipe is routinely given a carbon trace of 1e-9 mol. Measured at λ = 0.01 on that paste, a residual of 1.7e-10 mol on the charge row scores 168 against its own budget and a residual of 1.1e-10 mol on the carbon row scores 10.6 — both physically nothing, both rejected, and the walk then never reached its first three targets while still appearing to work. A row holding 1e-11 mol cannot be balanced better than the solver's absolute floor, and asking it to be is a category error.

scaleᵢ = max(|bλᵢ|, Σⱼ |Aᵢⱼ| nⱼ): the row's own budget, or the amount of that component actually being moved around when the budget is near zero — the natural yardstick for a conservation row that nets to nothing. atol sits above the accuracy the interior point itself reaches (about 3e-6 mol on this class of problem), because a rung is a guess and not an answer. Neither number certifies anything: they are there to reject a rung that has wandered, and the certificate judges the result afterwards.

ChemistryLab.activities Method
julia
activities(state::ChemicalState, model::AbstractActivityModel)
    -> OrderedDict{String,Float64}

Activity of every species — exp of log_activities.

Examples

julia
a = activities(eq, HKFActivityModel())
a["H2O@"]                    # water activity, e.g. 0.9937
ChemistryLab.activity_coefficients Method
julia
activity_coefficients(state::ChemicalState, model::AbstractActivityModel)
    -> OrderedDict{String,Float64}

Activity coefficient γᵢ of every aqueous species, from the model's formula.

Solutes are evaluated as γᵢ = 10^(log₁₀ γᵢ) with the model's own expression — −A zᵢ² √I/(1 + B åᵢ √I) + Ḃ I for the ions of HKFActivityModel, Kₙ I for its neutrals, and identically 1 for DiluteSolutionModel, which is ideal on its own (molarity) scale. The solvent is reported as γ_w = a_w / x_w.

Not computed as a ratio of activity to concentration. That ratio agrees for an abundant solute — and the tests check that it does — but it diverges for a species parked at the solver's lower bound, whose log-activity is dominated by the closures' + ϵ term: it returns values of order 1e300 for a charge class whose only members are trace. The formula depends on the ionic strength and the charge alone, so it is exact for a trace species and for a major one alike.

Throws if the system has no aqueous phase.

Examples

julia
γ = activity_coefficients(eq, HKFActivityModel(å = 0.0, Ḃ = 0.097637, Kₙ = 0.0))
γ["K+"], γ["Ca+2"]           # 0.6098, 0.1199 on a CEM I pore solution
γ["H2O@"]                    # the solvent, as a_w / x_w

See also: ionic_strength, concentration_scale, activities.

ChemistryLab.half_reaction Method
julia
half_reaction(state, oxidized, reduced) -> Reaction

The balanced half-reaction taking oxidized to reduced, over the species the system already contains plus H⁺, H₂O and ELECTRON.

Nothing is transcribed: the coefficients come from the element and charge balance, so the reaction is the one this system's species actually support.

julia
julia> half_reaction(eq, "SO4-2", "HS-").equation
"SO₄²⁻ + 9H⁺ + 8e⁻ = 4H₂O@ + HS⁻"
ChemistryLab.homotopy_initial_state Method
julia
homotopy_initial_state(state::ChemicalState; model = DiluteSolutionModel(),
                       steps = ..., ϵ = 1e-16, max_bisections = 6,
                       balance_atol = 1e-5, balance_rtol = 1e-3,
                       verbose = false)
    -> Union{ChemicalState, Nothing}

An automatic initial approximation obtained by continuation in the amount of solute.

Scales every species except the aqueous solvent by a factor λ, and walks λ through steps up to 1, solving at each value from the answer to the previous one. At λ = 1 the composition is exactly state, so the returned composition respects the same element balance; it is a starting point, not a certified equilibrium.

steps is a suggestion, not a schedule. A rung is accepted only if it conserves mass, and a refused rung is retaken by halving the distance back to the last λ that worked, up to max_bisections times per target. Both halves of that matter: a rung that has wandered off the balance would otherwise be carried forward as the start of every later rung, and a rung that simply cannot be taken in one jump can be taken in two.

The acceptance test is |rᵢ| ≤ balance_atol + balance_rtol · scaleᵢ on the element-balance residual, row by row, with scaleᵢ = max(|bᵢ|, Σⱼ |Aᵢⱼ| nⱼ). It has to be mixed rather than relative: two conservation rows carry a legitimately negligible budget — electroneutrality is exactly zero, and a cement recipe is routinely given a carbon trace of 1e-9 mol — so a relative test rejects residuals of 1e-10 mol as if they were failures. Neither tolerance certifies anything; balance_atol sits above the accuracy the interior point itself reaches (~3e-6 mol here) because a rung is a guess, and the certificate judges the result afterwards.

This is what makes a realistic cement solvable from a cold start: with all the mass in the reactants and every product at the ϵ floor, no back end reaches the optimum, because an interior-point method started against the boundary of a system spanning sixteen orders of magnitude in amount cannot take a useful step. At small λ the same system is dilute and well away from that boundary.

It requires nothing from the caller beyond the initial state: no guess at which phases will form, no knowledge of the answer. equilibrate_certified calls it automatically when its ordinary starting points fail to certify, so ordinary use never needs it.

The walk is done under model, which defaults to DiluteSolutionModel and should normally be left there even when the target is a non-ideal model. Measured on the CEM I paste: walking under the extended Debye-Hückel model with a common ion size of zero runs away to an ionic strength of 18 mol/kg, because that model's coefficients fall steeply with I, which raises solubility, which raises I. The ideal model has no such feedback, walks cleanly, and its answer is a good starting point for the non-ideal one — which is how equilibrate_certified uses it.

Returns nothing if the walk produced nothing usable.

Differentiability is unaffected. Under ForwardDiff, equilibrate_certified strips to the primal state, solves in Float64 and attaches the sensitivity through the implicit function theorem, so this never sees a Dual and the derivative does not depend on how the starting point was found.

Examples

julia
# What `equilibrate_certified` does for you when the cold start fails:
guess = homotopy_initial_state(state)
eq, cert = equilibrate_certified(guess; model = HKFActivityModel())

See also: equilibrate_certified.

ChemistryLab.ionic_strength Method
julia
ionic_strength(state::ChemicalState; kind = :effective, ϵ = 1e-16) -> Float64

Molality-basis ionic strength I = ½ Σⱼ mⱼ zⱼ², in mol/kg.

kind selects which of the two ionic strengths is meant, and they are different quantities:

  • :effective (the default) sums over the speciated free ions — the composition as solved, with every complex and ion pair counted at its own charge. A neutral pair such as Ca(SO4)@ contributes nothing. This is the one every activity model in this package is a function of, and the one to pass to any of them.

  • :stoichiometric sums as if every complex were fully dissociated into the components of state.system, so Ca(SO4)@ contributes as Ca²⁺ + SO₄²⁻. This is the analytical ionic strength of the recipe rather than of the solution, and it is what some correlations for salting-out and for diffusivity are fitted against.

The gap between them measures how much of the salt is associated. On a Portland cement pore solution it is small — 0.2121 against 0.2136 mol/kg, 0.7 % — because little is paired at that ionic strength; on a sulfate brine it is not.

Neither is a property of the activity model: both are properties of the composition. Compare the effective one before comparing anything else — an ionic strength that disagrees between two codes means they are not describing the same solution, whatever their volumes happen to agree on.

The stoichiometric sum uses the decomposition of each species over the system's primaries (state.system.SM.A), counting |νₚ| zₚ² for each charged primary p. That is the same decomposition the mass balance uses, so it needs no separate table of dissociation reactions; note that it attributes OH⁻, which CEMDATA18 writes as H₂O − H⁺, one unit of charge through the H⁺ component, which is the right count.

Throws if the system has no aqueous phase, or on an unknown kind.

Examples

julia
ionic_strength(eq)                          # 0.2121 mol/kg — free ions
ionic_strength(eq; kind = :stoichiometric)  # 0.2136 — fully dissociated

See also: molalities, activity_coefficients.

ChemistryLab.log_activities Method
julia
log_activities(state::ChemicalState, model::AbstractActivityModel;
               ϵ = 1e-16, kelvin_shift = 0.0)
    -> OrderedDict{String,Float64}

Natural log of the activity of every species, in the model's own convention.

This is the vector the Gibbs energy is built from: μᵢ/RT = ΔₐG⁰ᵢ/RT + ln aᵢ. Crystals of a pure phase get ln a = 0, solid-solution end-members ln a = ln xᵢ plus their excess term, the solvent its osmotic contribution, and solutes the log of their concentration in the model's scale plus ln γᵢ.

A species at the solver's lower bound has its value dominated by the closures' + ϵ regularization; read activity_coefficients rather than dividing these by a concentration.

Examples

julia
lna = log_activities(eq, model)
lna["H2O@"]                                  # ln a_w
exp(lna["Portlandite"])                      # 1.0 for a pure phase that is present

The capillary shift

kelvin_shift is added to the solvent's log-activity, and it is 0.0 by default so that nothing changes for a caller who does not ask.

It exists because an activity model computes the activity of water from the composition of the solution and knows nothing about the pore that holds it. Under CapillaryWater the solve is posed with a shifted solvent potential — the Kelvin term — and that shift is returned to the caller in the parameters reference, not stored in the state. Reading the activities back without it therefore reports the chemical value and not the pore value: after a solve posed at a_w = 0.90 this function returns 0.999995 unless the shift is passed in.

The two lowerings are independent and their chemical potentials add, so the activities multiply:

with kelvin_shift the logarithm of the second factor — a negative number, since a meniscus lowers the activity. water_activity is the accessor for the composed value.

Examples

julia
q = Ref{Vector{Float64}}()
eq, cert = equilibrate_certified(state; constraint = CapillaryWater(law; reference = fresh),
                                 parameters = q)
lna = log_activities(eq, model; kelvin_shift = q[][1])   # the pore water

See also: activities, activity_coefficients, water_activity, CapillaryWater.

ChemistryLab.molalities Method
julia
molalities(state::ChemicalState; ϵ = 1e-16) -> OrderedDict{String,Float64}

Molality mᵢ = nᵢ / (n_w Mw) of every aqueous solute, in mol per kg of solvent.

The solvent itself is not a solute and is omitted. ϵ floors the amounts the same way the activity closures do, so the values match what the solver saw; species at the floor come back at a molality of order ϵ rather than zero.

Throws if the system has no aqueous phase.

Examples

julia
m = molalities(eq)
m["K+"]                      # mol/kg of water
sum(values(m))               # total solute molality

See also: ionic_strength, activity_coefficients, concentration_scale.

ChemistryLab.pH Method
julia
pH(state::ChemicalState, model::AbstractActivityModel) -> Float64

−log₁₀ a(H⁺) — the pH in the activity convention of model.

This is what GEM-Selektor and Reaktoro report, and it is not what the one-argument pH returns: that one is −log₁₀ c(H⁺) with the concentration taken over the computed liquid volume, and in an alkaline solution it is reconstructed from OH⁻ through pKw. The two differ by the activity coefficient and by the scale conversion. On a Portland cement pore solution at I ≈ 0.2 mol/kg, with γ(H⁺) ≈ 0.61, the gap is about 0.21 units — large enough to be mistaken for a modeling error when comparing against another code.

Returns NaN if the system carries no H+.

Examples

julia
pH(eq)                       # 13.310 — concentration convention
pH(eq, model)                # 13.099 — activity convention, comparable to GEMS

See also: pOH, activity_coefficients, FixedpH.

ChemistryLab.pOH Method
julia
pOH(state::ChemicalState, model::AbstractActivityModel) -> Float64

−log₁₀ a(OH⁻) — the pOH in the activity convention of model.

See pH for why this differs from the one-argument pOH.

Returns NaN if the system carries no OH-.

ChemistryLab.pe Method
julia
pe(state, model; couple = "SO4-2" => "HS-") -> Float64

The electron activity of the pore solution as   , read off one redox couple.

There is no electron species to read an activity from, so pe is inferred from a half-reaction: the couple's two members are balanced over H⁺, H₂O and the electron, the reaction's is computed from the standard Gibbs energies the database carries, and the electron activity is what remains. For a half-reaction with n electrons on the oxidized side,

couple names the oxidized and the reduced member, in that order. The default is sulfate/sulfide, which is the couple a slag-blended cement buffers.

Different couples need not agree, and the disagreement is a result

A single pe exists only if every couple is at mutual equilibrium. In a real paste they are not: sulfate reduction is slow enough to be frozen on the time scale of hydration, so the sulfur couple and the iron couple can report potentials hundreds of millivolts apart. Computing pe from two couples and comparing them measures how far the assumption of a single redox state is from holding — which is worth doing before trusting either.

See also: Eh, half_reaction, FixedpE.

ChemistryLab.saturation_indices Function
julia
saturation_indices(state, model; ϵ = 1e-16) -> OrderedDict{String, <:Real}

LogSI for every species at state: log₁₀(IAP/K) of the reaction that forms it from the system's primaries.

Zero for a phase at equilibrium with the solution, negative for one that is undersaturated, positive for one that should have precipitated. It is the quantity GEM-Selektor prints as LogSI, and the one an optimality_certificate summarizes into a single worst violation without saying which phase that is.

No fitting is involved. The row labels of the conservation matrix are the primary species, so a component's element potential is that primary's chemical potential, y_c = μ_c/RT, and

julia
LogSI_s = [Σ_c A_cs y_c − μ_s/RT] / ln 10

which is saturation_ratio written for the formation reaction s = Σ_c A_cs (primary c). A conservation row that labels no species — the charge row — contributes nothing, since the coefficient of any neutral phase there is zero.

Two things to know before reading the numbers:

  • The check is built in. Every phase actually present at an equilibrium must come out at LogSI = 0; measured on a CEM I paste, the twelve present solids land within 1.2e-12. If they do not, the state is not an equilibrium and no other index in the result means anything.

  • A solid-solution end-member's index is relative to its current mole fraction, since its activity is ln x. For an end-member at the solver's lower bound that is a statement about a vanishing phase, not about whether the solid solution would form.

Examples

julia
si = saturation_indices(eq, model)
si["hydrotalcite"]                     # +5.58: absent, and it should not be
[k for (k, v) in si if v > 1e-4]       # everything supersaturated

See also: optimality_certificate, saturation_ratio, log_activities.

ChemistryLab.solvent_fraction Method
julia
solvent_fraction(state) -> Float64

Mole fraction of the aqueous solvent within the aqueous phase, n_w / Σ_aqueous n. One for pure water, and the number that says whether there is still a solution to speak of.

Every quantity built on the aqueous phase — molality, ionic strength, activity, pH — is defined per kilogram of solvent, and the dual solver parameterizes its interior variables by the solvent's chemical potential. All of that presumes the solvent is the phase, not one species in it. A real electrolyte, even a concentrated one, keeps x_w above about 0.9: seawater is 0.99, a saturated NaCl brine 0.90.

Below SOLVENT_FRACTION_FLOOR the formulation has no ground left. Measured on a sealed cement paste driven under its stoichiometric water demand — w/c = 0.28 on the mix of the w/c example — the Gibbs minimum consumes the free water down to 6e-9 mol, x_w falls to 0.21, and the ionic strength comes out at 409 mol/kg against a Debye-Huckel model valid to about one. The amounts such a solve reports are not equilibrium values; the answer is that the question was posed outside the model's domain.

See also: molalities, ionic_strength, equilibrate_certified, which checks this on the answer it returns.

ChemistryLab.water_activity Method
julia
water_activity(state::ChemicalState, model::AbstractActivityModel;
               ϵ = 1e-16, kelvin_shift = 0.0) -> Float64

The activity of water in state, composed of both lowerings.

The first factor is what the activity model computes from the composition. The second is the capillary term, which is not a property of the composition: it depends on the pore the water sits in, so it has to be supplied. Pass the shift that CapillaryWater returned through its parameters reference, or compute one from a retention law with water_activity(r, S; V_m, T) and take its logarithm.

Left at its default the function reports the chemical water activity alone, which is the right answer for a solution in a container and the wrong one for a solution in a gel pore.

Examples

julia
water_activity(eq, model)                          # chemical only
water_activity(eq, model; kelvin_shift = log(0.9)) # held at RH 90 % as well

See also: log_activities, CapillaryWater, kelvin_activity, PoreHumidity.

Solid solutions ​

ChemistryLab.AbstractSolidSolutionModel Type
julia
abstract type AbstractSolidSolutionModel end

Base type for activity models within a solid solution phase. Each concrete subtype implements _excess_ln_gamma(model, k, x, T).

ChemistryLab.AbstractSolidSolutionPhase Type
julia
abstract type AbstractSolidSolutionPhase end

Base type for solid solution phases. Concrete subtypes group a set of end-member AbstractSpecies and associate them with an AbstractSolidSolutionModel.

ChemistryLab.IdealSolidSolutionModel Type
julia
struct IdealSolidSolutionModel <: AbstractSolidSolutionModel

Ideal (Temkin) solid solution mixing: ln γᵢ = 0, hence ln aᵢ = ln xᵢ.

Valid for any number of end-members.

Example

julia
julia> m = IdealSolidSolutionModel()
IdealSolidSolutionModel()
ChemistryLab.RedlichKisterModel Type
julia
struct RedlichKisterModel{T<:Real} <: AbstractSolidSolutionModel

Binary Redlich-Kister (asymmetric Margules) model for non-ideal solid solutions. Requires exactly 2 end-members per solid solution phase.

Parameters a0, a1, a2 are in J/mol and divided by RT inside the activity closure.

Activity coefficients (Guggenheim / ThermoCalc convention):

julia
ln γ₁ = (x₂²/RT)[a₀ + a₁(3x₁ − x₂) + a₂(x₁ − x₂)(5x₁ − x₂)]
ln γ₂ = (x₁²/RT)[a₀ − a₁(3x₂ − x₁) + a₂(x₂ − x₁)(5x₂ − x₁)]

AD-compatible: all computations propagate ForwardDiff.Dual numbers.

Examples

julia
julia> m = RedlichKisterModel(a0 = 4000.0, a1 = 500.0)
RedlichKisterModel{Float64}(4000.0, 500.0, 0.0)

julia> m.a0
4000.0
ChemistryLab.RedlichKisterModel Method
julia
RedlichKisterModel(; a0=0.0, a1=0.0, a2=0.0)

Keyword constructor for RedlichKisterModel. Parameters are promoted to a common type.

ChemistryLab.RegularSolutionModel Type
julia
struct RegularSolutionModel{T<:Real} <: AbstractSolidSolutionModel

Symmetric regular (multi-component Margules) solid solution.

The gap this fills is arity. RedlichKisterModel is the general binary form and is restricted to two end-members, while IdealSolidSolutionModel takes any number but no interaction at all. A C-S-H with six end-members, or the CNASH and ECSH families of CEMDATA18, had no non-ideal option here.

Excess Gibbs energy, with one symmetric interaction parameter per pair:

julia
G^ex = Σ_{i<j} W_ij x_i x_j
ln γ_k = (1/RT) [ Σ_{j≠k} W_kj x_j  −  Σ_{i<j} W_ij x_i x_j ]

W is in J/mol, symmetric, with a zero diagonal; only the off-diagonal entries are read and W[i,j] is used for the pair (i,j). W_ij > 0 is repulsive — it favors unmixing — and W_ij = 0 recovers ideal mixing.

For two end-members this is exactly RedlichKisterModel(a0 = W₁₂), which the test suite checks; the point of it is n > 2.

AD-compatible: all computations propagate ForwardDiff.Dual numbers.

Examples

julia
julia> m = RegularSolutionModel([0.0 4000.0; 4000.0 0.0]);

julia> m.W[1, 2]
4000.0

References

  • Guggenheim, E.A. (1937). Trans. Faraday Soc. 33, 151–159.
ChemistryLab.RegularSolutionModel Method
julia
RegularSolutionModel(W::AbstractMatrix) -> RegularSolutionModel

Construct a RegularSolutionModel from a symmetric matrix of interaction parameters in J/mol. Raises if W is not square or not symmetric; the diagonal is ignored. Integer entries are promoted to a floating-point type.

ChemistryLab.SolidSolutionPhase Type
julia
struct SolidSolutionPhase{T<:AbstractSpecies, M<:AbstractSolidSolutionModel}
        <: AbstractSolidSolutionPhase

A solid-solution phase consisting of end_members (species with AS_CRYSTAL aggregate state) mixing according to model.

End-members are automatically requalified to SC_SSENDMEMBER at construction time, so database species with SC_COMPONENT can be passed directly.

Construction

Use the keyword constructor:

julia
SolidSolutionPhase(name, end_members; model = IdealSolidSolutionModel())

Validation at construction time:

  • All end-members must have aggregate_state == AS_CRYSTAL.

  • RedlichKisterModel requires exactly 2 end-members.

  • The mixing energy must be convex, unless instances > 1 or check_convexity = false.

A miscibility gap: instances

instances is how many coexisting compositions the declaration may hold. It is 1 for every phase the shipped data describes, and it is 1 because those models are convex: a convex mixing energy has one minimum, so one composition describes the phase.

Inside a spinodal it does not. Where d²g/dx² < 0 the Gibbs minimum is the common-tangent pair — two compositions of the same substance, coexisting — and a formulation carrying one amount per species cannot write that down. So instances = 2 asks ChemicalSystem for a second copy of each end-member, under a derived symbol (monosulphate12#2) sharing the same thermodynamic record, and the minimization is free to put material in either lobe or in both.

This is how GEM-Selektor represents the same thing: CEMDATA18 ships the AFm and AFt binaries under two names each, so that the user can declare them twice. The difference here is only that the duplication is asked for by a keyword rather than carried in the database.

instances > 1 is refused for a convex model, and that is not a formality: two instances of a convex phase are degenerate, every split of the amount between them having the same energy, so the minimum becomes a flat manifold and the optimizer is asked to choose a point on it for no reason. Inside a spinodal the common-tangent pair is unique and the degeneracy does not arise.

julia
# The published AFm sulfate/hydroxide parameters, whose spinodal is
# x in [0.631, 0.914] at 25 C. With one instance this is refused; with two it is
# the case the model was written for.
SolidSolutionPhase("AFm_SO4_OH", [c4ah13, monosulphate];
                   model = RedlichKisterModel(a0 = 20_000.0), instances = 2)

Example

julia
julia> em1 = Species("Ca2SiO4"; aggregate_state=AS_CRYSTAL, class=SC_COMPONENT);

julia> em2 = Species("Ca3Si2O7"; aggregate_state=AS_CRYSTAL, class=SC_COMPONENT);

julia> ss = SolidSolutionPhase("CSH", [em1, em2])
SolidSolutionPhase{Species{Int64}, IdealSolidSolutionModel}
  name: CSH
  end-members (2): Ca2SiO4, Ca3Si2O7
  model: IdealSolidSolutionModel

julia> class(end_members(ss)[1])
SC_SSENDMEMBER::Class = 5
ChemistryLab.SolidSolutionPhase Method
julia
SolidSolutionPhase(name, end_members; model=IdealSolidSolutionModel())

Construct and validate a SolidSolutionPhase.

End-members whose class is not already SC_SSENDMEMBER are automatically requalified via with_class. Passing database species with SC_COMPONENT therefore works directly, without a prior call to with_class.

ChemistryLab._mixing_energy Method
julia
_mixing_energy(A0, A1, A2) -> Function

Molar Gibbs energy of mixing of a binary, in units of RT:

The ideal part is convex everywhere — its second derivative is 1/x + 1/(1-x) — so every miscibility gap is the excess term's doing.

ChemistryLab._rk_coefficients Method
julia
_rk_coefficients(model, T) -> Union{Nothing, NTuple{3,Float64}}

The three Redlich-Kister coefficients of a binary mixing model, in units of RT, or nothing when the model has no excess term this form can express.

RegularSolutionModel is the one-parameter case, a₀ = W₁₂ with a₁ = a₂ = 0, which is why the two share this. IdealSolidSolutionModel — and any model a caller adds — returns nothing: an ideal mixture is convex everywhere, so every construction below is vacuous for it.

Factored out because three functions need exactly this and had three copies of it, which is two opportunities for them to disagree.

ChemistryLab.common_tangent Function
julia
common_tangent(model; T = 298.15, tol = 1e-12, maxit = 100)
    -> Union{Nothing, Tuple{Float64,Float64}}

The two compositions a binary solid solution separates into inside a miscibility gap, or nothing when its mixing energy is convex.

What this computes, and why it is not a minimization

Where the molar Gibbs energy of mixing is concave the equilibrium is not one composition but two, and they are the pair at which a single straight line is tangent to twice — equivalently, at which both end-members have equal chemical potentials in the two phases:

In terms of alone that is

two equations in two unknowns, solved here by Newton with ForwardDiff for the derivatives. It does not involve the rest of the chemical system at all: the pair depends only on the mixing model and the temperature, which is what makes it computable in microseconds and usable as the starting point of a full equilibrium.

This is the approach PHREEQC takes for binary solid solutions, after (Glynn and Reardon, 1990), and it is a different thing from asking a global minimization to discover the split. A minimization started from two identical compositions sits on a stationary point: both instances satisfy every first-order condition jointly, so there is no downhill direction to follow, and it stays there however unstable the state is. Handing it the pair removes the question.

The binodal contains the spinodal

spinodal_interval reports where  , which is where the phase is unstable. The pair returned here is wider: between the two the phase is metastable rather than unstable, and a minimization sees only the tangent. So the spinodal edges bracket the search from inside, and that is where Newton starts.

Verification

For a symmetric model   , hence    , and the common-tangent condition collapses to  :

The test suite checks the returned pair against that equation rather than against a stored number.

Returns nothing for an ideal model, for a phase with more than two end-members (where a one-dimensional construction is not the right object), and when Newton does not converge — never a guess.

See also: spinodal_interval, SolidSolutionPhase.

ChemistryLab.common_tangent Method
julia
common_tangent(phase::SolidSolutionPhase; T = 298.15) -> Union{Nothing, Tuple}

As above for a declared phase.

ChemistryLab.end_members Method
julia
end_members(ss::SolidSolutionPhase) -> Vector{<:AbstractSpecies}

Return the end-member species of the solid solution phase.

ChemistryLab.miscibility_split Function
julia
miscibility_split(model, x̄; T = 298.15) -> NamedTuple

How a binary of overall composition x̄ separates inside its miscibility gap.

Returns (; x_alpha, x_beta, f_alpha, f_beta, Δg) — the two coexisting compositions, the mole fraction of the binary in each, and the molar Gibbs energy the separation releases, in J/mol. Returns nothing when the model is convex, and a single phase (f_alpha = 1, Δg = 0) when x̄ lies outside the pair.

The construction

The compositions come from common_tangent and do not depend on x̄: inside a gap the two phases in equilibrium always have the same pair of compositions, only their proportions change. Those proportions are then the lever rule,

which is mass balance and nothing more:   .

Δg is the distance from the curve down to the common tangent at x̄,

so it measures, in J/mol of binary, how much a single-composition answer overstates the Gibbs energy — that is, how wrong it is.

This is the construction of (Glynn and Reardon, 1990), the one PHREEQC uses for a binary solid solution.

What this does and does not settle

Given x̄, everything above is exact and costs microseconds. Obtaining x̄ from a full aqueous equilibrium inside a gap is the part a minimization over two declared instances does not currently deliver: the symmetric state is a stationary point, and the two-instance problem carries a near-null direction that more iterations make worse rather than better (measured: the element balance degrades from 1.5e-01 to 4.5e+00 between 200 and 5000 iterations).

So read x̄ as the overall composition you have — from a single-phase solve, from an analysis, or as a scan — and this as the exact answer for it.

See also: common_tangent, spinodal_interval.

ChemistryLab.miscibility_split Method
julia
miscibility_split(phase::SolidSolutionPhase, x̄; T = 298.15)

As above for a declared phase.

ChemistryLab.model Method
julia
model(ss::SolidSolutionPhase) -> AbstractSolidSolutionModel

Return the activity model of the solid solution phase.

ChemistryLab.name Method
julia
name(ss::SolidSolutionPhase) -> String

Return the name of the solid solution phase.

ChemistryLab.spinodal_interval Method
julia
spinodal_interval(model, n_members; T = 298.15) -> Union{Nothing, Tuple{Float64, Float64}}

The interval of composition over which a binary mixing energy is concave, or nothing when it is convex throughout.

A solid solution exists as one homogeneous phase only where its molar Gibbs energy of mixing is convex. Where d²g/dx² < 0 — inside the spinodal — the minimum of G is not a single composition but two coexisting ones: the phase unmixes, and the equilibrium is a miscibility gap.

For a binary with the package's Redlich-Kister convention,

and the second derivative is evaluated on a grid rather than in closed form, so that the same routine covers a₂ and the symmetric RegularSolutionModel without a separate derivation. The classical symmetric result is recovered as a check: d²g/dx² at x = 1/2 is 4 − 2A₀, so a regular solution unmixes above W = 2RT.

T is the temperature the parameters are read at; they are stored in J/mol and the criterion is a/RT, so a model that is convex at 25 °C may not be at 5 °C.

Returns nothing for an ideal model, whose second derivative is 1/x + 1/(1-x) and therefore positive everywhere, and for a phase with more than two end-members, where a one-dimensional scan is not the right test — see the note in SolidSolutionPhase.

See also: RedlichKisterModel, RegularSolutionModel.

Certified equilibrium ​

Newton on the KKT system in element-potential space, and the certificate that proves a composition optimal. See Proving that an answer is the answer for the derivation, and equilibrate_certified for the route equilibrate takes by default.

ChemistryLab.DualEquilibriumSolver Type
julia
DualEquilibriumSolver(system, model; tol, maxit, max_active_updates, si_tol, verbose)

Equilibrium by Newton on the KKT system, in element-potential space.

Where EquilibriumSolver minimizes G by an interior-point method over all species and stops on MaxIters for a cement, this solves the stationarity conditions directly and returns a result that optimality_certificate can prove optimal — the Gibbs problem being convex, the KKT conditions are sufficient.

Writing u = −Aᵀy for the element potentials, an aqueous species obeys the mass-action law aᵢ = exp(uᵢ − gᵢ) and a pure phase is present exactly when uᵢ = gᵢ, absent when undersaturated: the classical phase-stability criterion.

Starting from the answer of an EquilibriumSolver is the intended use — the interior-point method reaches a neighborhood, this reaches the conditions.

See also: optimality_certificate, speciated_states.

ChemistryLab._dual_phases Method
julia
_dual_problem(des, p, n0) -> DualNewtonProblem

Package the chemistry as the convex program OptimaSolver solves. Built per solve because the reference potentials Δ_a G⁰/RT depend on temperature and pressure.

ChemistryLab._kkt_error Method
julia
_kkt_error(cert) -> Float64

How far a composition is from satisfying the KKT conditions: the worst of the three residuals the certificate reports, in one number.

All of them, and not the stationarity alone. A composition can be stationary to 1e-3 while violating mass conservation by moles — measured, an answer with stationarity 2.4e-3, an element balance off by 6.7 mol and a phase supersaturated by 45 — and that is not a near-answer, it is not an answer to this problem at all. Ranking on stationarity alone lets such a point beat a candidate that conserves mass, which is how a multi-start search can discard the good answer it just computed.

The constraint's own residual counts too, when there is one. Leaving it out would rank a candidate that minimizes the Gibbs energy while violating the very equation that makes it a constrained answer above one that satisfies both — the same hole OptimaSolver records for the kinetic step, where "a march that should have stopped at saturation dissolved everything and was proved optimal". Read with hasproperty because a certificate also arrives from the kinetic route, which builds its own.

worst_supersaturation is clamped at zero because a negative value is not an error: it means every absent phase is undersaturated, as optimality requires.

ChemistryLab._split_starts Method
julia
_split_starts(model, nmembers) -> Vector{Vector{Float64}}

Where to look for the other lobe of a phase that may unmix, as mole fractions over its members.

The tangent-plane search in OptimaSolver probes the corners of the composition simplex and refines by successive substitution, which converges to the stationary point nearest its start. A corner is usually on the right side of the barrier and sometimes is not; where it is not, the iteration walks back to the phase's own composition and reports nothing, though splitting would lower the energy.

Measured on the AFm sulfate/hydroxide binary with the published Redlich-Kister parameters (spinodal [0.631, 0.914], binodal [0.4999, 0.9700]):

phase sits atcorners alonewith the binodal handed over
x = 0.5268, where a CEM I settles+3.99e-02, foundthe same verdict
x = 0.95+2.4e-16, missed+1.23e-01, trial x = 0.444
x = 0.98, outside the binodalstablestable

Both flagged compositions are metastable, so being metastable is not by itself what defeats the corners — sitting in the lobe they lead back into is, and that is not knowable in advance. Hence a start supplied unconditionally rather than a rule for when to supply one.

common_tangent computes the binodal from the mixing model alone, in microseconds and with no reference to the rest of the system, and the search then refines it with the chemical potentials the system actually has — which is the pair that matters. That division is the whole point: the model's binodal is a good place to look, not the answer. Extra starts can only raise the maximum the search returns, so they never take a verdict away and, as the last row shows, do not invent one.

Empty for anything but a binary with a gap, which is the only case common_tangent is defined for.

ChemistryLab._ss_models Method
julia
_ss_models(des) -> Dict{Int, Any}

The mixing model of each entry of des.ss_groups, by position.

ss_groups is built by walking system.solid_solutions and skipping any declaration whose end-members are not all in the species list, so the two lists are the same length only when nothing was skipped. Rebuilding the correspondence the same way is the only way to be sure a model is matched to its own group; a positional zip would silently pair a phase with someone else's model the first time a declaration is dropped.

ChemistryLab.optimality_certificate Method
julia
optimality_certificate(des, state; b = nothing, ϵ = 1e-16, floor = 1e-25)
    -> (; stationarity, balance, worst_supersaturation, n_interior,
         n_absent_component, param_residual, worst_violation_split,
         split_phases, split_trials, optimal)

Check the KKT conditions at a composition, independently of how it was obtained.

For a convex problem these conditions are sufficient, so optimal = true is a proof of global optimality. Use it to audit any solver — including EquilibriumSolver, whose interior-point iteration reports MaxIters on a cement equilibrium and cannot say whether the point it returns is the answer.

The three quantities are the stationarity of the interior species, the component balance, and the worst saturation index among absent phases (negative when every one of them is undersaturated, as optimality requires).

worst_violation_split extends that last test to the phases that are present: Michelsen's tangent-plane distance, which asks whether a mixing phase would lower the Gibbs energy by separating into two compositions. It is -Inf when no phase could be tested, negative when every one of them is stable, and positive when one wants to unmix — split_phases then names them and split_trials carries, per phase, the composition it wants to split into. That composition is what equilibrate_split seeds a second instance with; it exists nowhere else, being a property of the full system and not of the mixing model alone.

ChemistryLab.solve_certified Method
julia
solve_certified(des, starts; b = nothing, ϵ = 1e-16, floor = 1e-25)
    -> (state, certificate)

Solve from each starting composition in starts and return the first answer optimality_certificate proves optimal. If none is proved, return the one with the smallest KKT error, together with its certificate, so the caller sees what it is getting.

Why trying more than one start is the rigorous thing to do, not a fudge

The certificate is an oracle: for a convex problem it does not rank answers, it decides them. Given an oracle, running several routes and keeping a proved answer is not guesswork — the proof is the same proof whichever route produced the point, and nothing about it depends on having predicted the winner. What would be a fudge is picking a route by taste and reporting its output unproved.

It is also necessary, because no single route dominates. Measured on an LC³ equilibrium at three degrees of reaction, with the dual Newton started from each interior-point back end in turn:

degree of reactionfrom OptimaOptimizerfrom IpoptOptimizer
0.05certified, 1.8e-12certified, 9.1e-13
0.25not certified, 7.2certified, 4.2e-12
1.00certified, 1.2e-11not certified, 3.5e-3

Each back end solves what the other misses. With both offered, all three are proved.

Example

julia
using Optimization, OptimizationIpopt      # for IpoptOptimizer

starts = [
    SciMLBase.solve(EquilibriumSolver(cs, model, OptimaOptimizer()), st; b = b),
    SciMLBase.solve(EquilibriumSolver(cs, model, IpoptOptimizer()), st; b = b),
]
eq, cert = solve_certified(des, starts; b = b)
cert.optimal || @warn "no start produced a certifiable answer" cert

The starts are supplied by the caller rather than built here, so this adds no dependency: whichever back ends are loaded are the ones available.

CommonSolve.solve Method
julia
SciMLBase.solve(des::DualEquilibriumSolver, state; b = nothing, ϵ = 1e-16)
    -> ChemicalState

Equilibrium composition. state supplies the temperature, the pressure and the starting guess; b the component totals, defaulting to those of state.

The amounts are returned as computed, without a floor: an amount of e⁻³⁰⁰ IS the mass-action answer for a species that is not there, and raising it to ϵ falsifies its activity by hundreds of RT units — which is enough to make optimality_certificate report a residual of 74 on a composition solved to 5e-12.

ChemistryLab._MAX_RESTARTS Constant
julia
_MAX_RESTARTS

How many times equilibrate_certified may restart from its own answer before giving up. One round is what the measured cases need; the bound exists so a case that improves by a hair every round cannot loop.

ChemistryLab._REPAIR_FRACTION Constant
julia
_REPAIR_FRACTION

What fraction of the amount the recipe could make of a missing phase _repair_start puts in. Far enough off the boundary for the active-set loop to work with, small enough not to pretend it knows the answer.

ChemistryLab._LazyStarts Type
julia
_LazyStarts

The starting points offered to solve_certified, each back end solved only when the search actually asks for it, and cached once it has been.

solve_certified returns at the first start that certifies, so solving every registered back end before the search begins pays for answers the search may never look at.

Measured honestly, the gain depends entirely on whether the first start certifies. On the CEM I of docs/src/examples/cem1_solid_solutions.md it does not, the later starts are genuinely consumed, and this buys nothing: 107.9 s against 105.1 s eager, within noise. The saving appears only where the first back end settles the problem, which is the common case for the smaller systems. So this is a correctness-of-effort change — never compute a start no one reads — and not a cement optimization; the expensive cement case needs a different answer.

Cached, and that is not optional: the route search offers the same starts again after the ideal pre-solve and again after the continuation, so recomputing them each time would more than undo the gain.

Iterating yields each back end's answer in registration order, skipping any that threw, and finally tail — the caller's own state, which is the only start left if every back end failed.

ChemistryLab._amount_number_type Method
julia
_certified_dual_route(state, model, b, ϵ, verbose, constraint, parameters, kwargs)

nothing for a real-valued composition; the certified answer with its derivative attached for one carrying ForwardDiff.Dual amounts.

Dispatched on the element type, positionally, so the choice is the type system's and neither path pays for the other. Two paths are needed for a mathematical reason, not for want of a generic element type: making the component totals generic would let duals flow into the solver, and what came back would be the derivative of the algorithm — an active set decided by sign tests, a line search with branches, an iteration count that varies with the data — rather than the derivative of the solution. The map b ↦ n*(b) is smooth only piecewise, and on each piece the implicit function theorem gives its derivative at the solution with the active set frozen. That is how OptimaSolver computes its own Sensitivity, and how Optima does upstream.

ChemistryLab._check_solvent Method
julia
_check_solvent(eq)

Say so when the answer has no solution left to be an answer about.

A certificate proves that a composition minimizes the Gibbs energy of the problem as posed. It says nothing about whether the problem was posed inside the model's domain, and there is one way to leave it that produces an ordinary-looking ChemicalState: let the solids take all the water. Every aqueous quantity is then computed per kilogram of a solvent that is not there.

Measured on a sealed cement paste below its stoichiometric water demand, at w/c = 0.28: the free water goes to 6e-9 mol, the solvent falls to a fifth of its own aqueous phase, and the ionic strength is reported as 409 mol/kg by a Debye-Huckel model valid to about one. Nothing in the certificate objects, because nothing is wrong with the minimization — the mix simply does not contain enough water to be a solution chemistry problem.

Warned, not raised, under the default flag: the composition of the solids in such a solve still carries the mass-balance information a caller may legitimately want, and it is the caller who knows whether that is what they asked for. Under STRICT_CONVERGENCE[] it raises, like any other answer that is not one.

ChemistryLab._dual_applicable Method
julia
_dual_applicable(system) -> Bool

Whether DualEquilibriumSolver can be built for system: it needs an aqueous phase, and H2O@ among the species, because it parameterizes the interior variables by the solvent's chemical potential.

ChemistryLab._ideal_start Method
julia
_ideal_start(state, model, bfix, ϵ, constraint, verbose; kwargs...)
    -> Union{ChemicalState, Nothing}

A certified answer to the same problem under ideal activities, to be used as a starting point. nothing when that solve does not certify either, or when model is already the ideal one.

The easier question is the useful one here. Without activity coefficients the residual does not depend on the composition through a second, non-linear path, so the solve is far better conditioned and certifies where the non-ideal model does not; and the phases it finds are the same ones — they differ in amount, not in identity — so the non-ideal solve that starts from it begins with the correct active set instead of discovering it.

That discovery is what was fragile. On a CEM I at w/c = 0.5 with the eight distinct CEMDATA18 solid solutions, eighty phases sit at the bound in the cold state and the active-set search decides its route on comparisons of nearly equal quantities: 100/sum(oxides) summed over a Dict and over an OrderedDict differ by one ulp, and that was enough to choose between an equilibrium certified to 1.1e-14 and a failure with an element balance of 71 mol. Started from the ideal answer, both reach the same certified composition.

Any failure of the inner solve is swallowed: this builds a starting point, and a caller who asked for a result is entitled to the outer verdict rather than to an error raised inside a heuristic.

ChemistryLab._instance_pairs Method
julia
_instance_pairs(cs) -> (twin, untwin)

The species index of each #2 twin, by the index of the species it copies, and the inverse map.

Both directions, because Michelsen's trial can be reported on either instance of a pair and the pair has to be recovered from whichever it names. Empty when the system was not declared with instances = 2, which is how equilibrate_split knows there is nowhere to split into.

ChemistryLab._keep_better Method
julia
_keep_better(eq, cert, eq2, cert2) -> (eq, cert)

Keep the better of two answers: a certificate of optimality beats none, and otherwise the smaller KKT error wins. A round that buys nothing changes nothing, which is what lets the restart loop run without ever making the answer worse.

The optimality flag is compared first, in both directions. Ranking on the residuals alone would let an uncertified point displace a certified one, and no residual is worth trading a proof for.

Among uncertified answers the comparison is _kkt_error — the worst of the three residuals — and not the stationarity alone, which is the same ranking solve_certified uses internally. The distinction is not academic: an answer stationary to 2.4e-3 whose element balance was off by 6.7 mol beat every candidate the continuation produced, because those were stationary to only 1e-2 while conserving mass. The search computed a usable answer and discarded it.

ChemistryLab._repair_round Method
julia
_repair_round(eq, cert, model, bfix, ϵ, solve_from, verbose)
    -> (eq, cert, improved)

One round of "the certificate names a missing phase, so put it in and solve again". Returns the better of the two answers and whether the round bought anything, so the caller can stop as soon as it does not.

solve_from is the search, injected: given a starting composition it returns (state, certificate). Passing it in rather than closing over the caller's solver is what makes this round testable on its own — the situation it exists for, a back end that converges onto the wrong active set, needs a system of some 135 species to arise, while the round's logic needs eight.

ChemistryLab._repair_start Method
julia
_repair_start(eq, model, bfix, ϵ) -> Union{ChemicalState, Nothing}

Build a starting point that makes the phases the certificate says are missing actually present. nothing when there are none.

This is what turns a diagnosis into a repair. optimality_certificate reports a positive worst supersaturation when a phase sits at the lower bound while the solution is supersaturated with respect to it — a genuine KKT failure on a convex problem, so the active set is wrong and the answer is not the answer. saturation_indices says which phase, and the obvious move is then to put it in and solve again.

The failure it exists for is a phase swap, which an active-set loop that admits one phase at a time cannot perform. Measured on a CEM I paste, reported from one machine while another certified the same source: all 0.02515 mol of magnesium sat in brucite with hydrotalcite absent and supersaturated by 5.58 log units, and admitting the hydrotalcite requires dissolving the brucite entirely and taking aluminum back from the hydrogarnet in the same step. The solve was otherwise impeccable — stationarity 1.5e-16, element balance 1.8e-14 — which is exactly what a converged-onto-the-wrong-active-set answer looks like.

The amount each missing phase is given is what the recipe could make of it, min_c b_c / A_cs over the components it consumes, scaled by _REPAIR_FRACTION. That is a chemical bound, not a guess at the answer: a carbonate in a system with 1e-9 mol of carbon is offered 1e-9 mol and no more. What matters is only that the phase starts well away from the boundary, since being at the boundary is what the active-set loop cannot recover from. The element balance of the result is not this function's business — b is fixed once by the caller and every start is projected onto it, so a start that over-spends the budget costs nothing.

ChemistryLab._seed_split! Method
julia
_seed_split!(n, trials, twin, untwin, share) -> Bool

Move material from the fuller instance of each flagged phase into the emptier one, at the composition Michelsen's analysis asks for. Mutates n and returns whether anything moved.

The element budget is not touched, and that is the whole design

The transfer takes the incipient composition t.x out of one instance and puts the same t.x into the other, so the total of every end-member across the pair is unchanged and A n is exactly what it was. The obvious alternative — remove at the DONOR's composition and add at the trial's — moves the same number of moles but a different mixture, so it silently rewrites the element budget the solver is about to be measured against. Two end-members of one binary are different substances; C4AH13 and monosulphate12 do not have the same sulfur.

move is then bounded so that no end-member of the donor goes negative, which is what makes share a request rather than a command.

From the fuller into the emptier

Not the other way round, and not from whichever instance the trial happened to name. The material is typically all in one of the two, and moving a share of an instance holding 1.5e-4 mol while its twin holds 5.1e-2 is a perturbation of three parts in a thousand — a seed that cannot move the answer is indistinguishable from no seed at all.

ChemistryLab._within_domain Method
julia
_within_domain(eq) -> Bool

Whether the answer is inside the domain the model is written for: an aqueous system whose solvent has been eaten by the solids is not.

This is a ranking question, not a diagnosis. _check_solvent already reports such a state on the answer that is returned; what it could not do is stop one from being chosen. Measured on a CEM I with eight solid solutions: a wandering iterate came back with the solvent at x_w = 0.033 and an element balance of 71 mol, the continuation's answer stood at 160, and the smaller number won — so the search settled on a composition in which the water had gone into the solids, and every route that conserved mass was discarded behind it.

Ranking on residuals alone cannot separate the two: both are large, and one is merely larger. The distinction that matters is not how big the residual is but whether the point is a composition this model can describe at all.

ChemistryLab.equilibrate_certified Method
julia
equilibrate_certified(state; model, ϵ, b, verbose, autostart) -> (state, certificate)

Equilibrium composition together with a proof of its global optimality, obtained by solving from every registered back end and keeping the answer optimality_certificate proves optimal.

The starting point is found, not asked for

When no back end certifies from the state as given, an initial approximation is computed by continuation — homotopy_initial_state — and every back end is run again from it. This is what makes a realistic cement solvable without the caller knowing anything about the answer: from the cold state of a CEM I paste (all the mass in the reactants, every product at the ϵ floor) no route reaches the optimum, and with the continuation the same call certifies.

It costs nothing in the ordinary case, because it only runs when nothing else certified. autostart = false declines it, which is what the coupled kinetic step does: there the caller already supplies the previous instant as a warm start, and a handful of extra solves inside an implicit ODE step would be paid at every step.

certificate.optimal == true is a proof, valid because the Gibbs minimization is convex when the mixing terms are — ideal mixing and any activity model whose excess Gibbs energy is convex in the amounts. It is not a proof for a model that is not, and none of the activity models that ship here have been shown to violate it; HKFActivityModel, DaviesActivityModel and the Redlich–Kister solid solutions are used within their stated ranges.

When no route yields a proof, the answer with the smallest KKT error is returned, its certificate says so, and a warning names the residual. That is the honest outcome, and it is not the same thing as a failure: on a low-water cement the returned composition satisfies the element balance to 1e-15 and has a supersaturated phase left out, which the certificate reports as worst_supersaturation > 0.

Requires OptimaSolver (the dual Newton lives there). Systems without an aqueous phase, or without H2O@, cannot use the dual route; for those, this falls back to the plain equilibrate and returns nothing as the certificate.

Example

julia
using ChemistryLab, OptimaSolver
eq, cert = equilibrate_certified(state)
cert.optimal          # true — proved globally optimal
cert.balance          # element balance residual
cert.worst_supersaturation   # negative: every absent phase undersaturated
ChemistryLab.equilibrate_path Method
julia
equilibrate_path(state, budgets; model, kwargs...) -> (states, certificates)

A sequence of certified equilibria, each one started from the last that certified.

budgets is any iterable of element budgets — the vectors equilibrate_certified takes as b. The first is solved from state; every later one is solved from the previous answer, and a previous answer is reused only once the certificate has accepted it. Where none has yet, state is used again.

Why this exists

A cement equilibrium is hard to start cold and easy to start warm, and the gap is not marginal. Measured on the 135-species paste of scripts/ionic_hydration.jl:

cold start, the full multi-start cascade15.2 s
warm start from a neighboring answer0.19 s

Eighty to one. So a sweep that rebuilds its state at every point pays the cold price at every point, and — worse — can fail at one while both of its neighbors certify, which is a starting point and not an infeasibility. Both blended-binder sweeps in this package's documentation had such a point before they were written this way.

What it does and does not change

For a convex problem the minimum is unique, so walking to it cannot change what is found — only whether the search finds it. That premise is checked rather than assumed: SolidSolutionPhase refuses a mixing model whose energy has a spinodal, so a system that was constructed at all is convex unless the refusal was explicitly waived. Waive it and this becomes a genuine choice of branch, because inside a gap the starting point decides which lobe the answer lands in — see common_tangent.

The certificate still decides every point. A refused point is returned like any other, with its certificate, and does not become the next start.

Examples

julia
budgets = [budget_at(f) for f in 0.0:0.05:0.30]
states, certs = equilibrate_path(fresh, budgets; model = HKFActivityModel())
all(c.optimal for c in certs)      # every point proved, not merely converged

See also: equilibrate_certified, common_tangent.

ChemistryLab.equilibrate_split Method
julia
equilibrate_split(state; model, b, maxpasses = 3, share = 0.5, kwargs...)
    -> (state, certificate)

The certified equilibrium of a system whose mixing phases may unmix, found by giving a phase that wants to split a second composition to split into.

The problem this solves

Inside a miscibility gap the Gibbs minimum of a mixing phase is two coexisting compositions, not one. A formulation carrying one amount per species describes that by declaring the phase twice — SolidSolutionPhase(...; instances = 2) — and when the element balance pins the phase's overall composition inside the gap, declaring it is enough: the minimization separates the two instances onto the common-tangent pair by itself, and the certificate proves it. Measured on a calcite/magnesite binary, where 0.025 mol of each fixes x̄ = 1/2 whatever the energetics say, the instances land on the binodal to within 1e-3 and in the proportions the lever rule asks for.

This function is for the case where nothing pins it. In a cement paste the AFm composition is free — the sulfate has ettringite to go to and the hydroxide is abundant — so two instances started at the same composition stay there. The symmetric state satisfies every first-order condition jointly, so it is a stationary point of the minimization, and no descent direction leads away from it however unstable it is. The composition is typically metastable rather than unstable — outside the spinodal, inside the binodal — where reaching the pair needs a finite jump and not a gradient step.

What it does

Michelsen's stability analysis, which is what the certificate already runs on every present mixing phase, does not only answer whether a phase splits: the trial composition that minimizes the tangent-plane distance is an estimate of the incipient phase, and that is what seeds the second instance here. The pass is then repeated until the certificate accepts or stops improving.

The seed comes from the analysis of the full system — the trial composition is computed with the chemical potentials the pore solution actually has — and not from the mixing model alone. common_tangent gives the binodal of the isolated binary, which is a different pair and is the wrong place to start from: seeding there was measured to collapse straight back to one composition.

What it returns, and what it promises

The best certified answer found, or — if none certifies — the last one, exactly as equilibrate_certified would. A pass is kept only when the certificate's KKT error improves — the worst of stationarity, element balance, supersaturation and the constraint residual, not one of them — so the result is never worse than the answer without splitting.

share is how much of the phase's amount is moved into the incipient instance on each pass; maxpasses bounds the work.

This is where convexity has already been given up

A phase that unmixes has a concave mixing energy, so G is not convex and cert.optimal no longer proves a global minimum — it proves a KKT point whose present phases are additionally stable against splitting, which is strictly more than stationarity gives. SolidSolutionPhase refuses such a model unless instances > 1 is asked for, so a system reaching this function was built deliberately.

See also: equilibrate_certified, common_tangent, miscibility_split.

Constraints ​

What is held fixed while the Gibbs energy is minimized. See Constraints other than fixed T and P.

ChemistryLab.Adiabatic Type
julia
Adiabatic()

Enthalpy conserved at the value the initial state carries, temperature unknown. This is a reaction in a vessel that exchanges no heat: the enthalpy released by the reaction goes into raising the temperature.

Equivalent to FixedEnthalpy(enthalpy(state)), and preferable because it cannot be given an enthalpy from a different element basis by mistake.

ChemistryLab.CapillaryWater Type
julia
CapillaryWater(retention, V_ref)
CapillaryWater(retention; reference)

Hydration arrested by self-desiccation: the water left in the pore space is held at reduced activity, and the reaction stops when that activity no longer supports the hydrates.

This is the only constraint here that adds physics rather than a boundary condition. The others say what is held fixed; this one says that water in a fine pore is not the same water as bulk water of the same composition. Its chemical potential is lower by the Kelvin term, RT ln(a_w/a_w^chem) = -2 γ V_m / r, so hydrates that consume water become less stable as the pore space empties. A sealed paste therefore stops hydrating with water still in it — which is what Powers' α_max = w/c / 0.42 describes empirically, and what powers_alpha_max carries into the kinetic rate laws.

retention is a WaterRetention: the material's own relation between the degree of saturation of its pore space and the activity of the water in it. It is measured, not assumed, and it is the caller's to supply.

V_ref is the volume the pore space is referred to — the fresh paste, exactly as for the two-argument porosity. Pass the fresh state and the volume is taken from it. The degree of saturation the constraint uses is then S = V_liquid / (V_ref - V_solid), both volumes recomputed from the current composition at every evaluation.

What it does to the solve

One unknown, q[1] = ln(a_w/a_w^chem) ≤ 0, the Kelvin shift; one equation, the retention law evaluated at the current saturation. Mass is conserved — unlike FixedActivity, no titrant column is added, because a sealed specimen exchanges water with nothing. The shift reaches the caller through the parameters keyword, as the titrant amount does.

What the certificate then proves

Less than usual, and the difference matters. For the unconstrained problem G is convex, so cert.optimal is a proof of a global minimum. A composition- dependent shift of ln a_w is not derived from a convex G for an arbitrary retention law, so under CapillaryWater the certificate proves a KKT point of the constrained problem — stationarity, mass balance, no absent phase supersaturated, and the capillary closure satisfied — but not global optimality. The multi-start route still runs, so agreement across starts is evidence; it is not the proof the fixed-(T, P) route gives.

Reading the shift back

The shift has to be passed to the accessors

The shift lives in the solver's parameter block and not in the activity model, because it is a property of the pore rather than of the composition. So log_activities and everything built on it — activities, pH(state, model), saturation_indices — report the chemical water activity by default, which this constraint barely moves. Both log_activities and water_activity take a kelvin_shift keyword for the composed value:

julia
q = Ref(Float64[])
eq, cert = equilibrate_certified(state; constraint = c, parameters = q)
water_activity(eq, model; kelvin_shift = only(q[]))        # what the solver saw
log_activities(eq, model; kelvin_shift = only(q[]))["H2O@"]

Measured on a calcite system with the shift set to ln 0.90: the chemical log-activity comes back at -5.3e-6, i.e. a_w = 0.999995, while the activity the solve was posed with is 0.90. A caller who reads the default and concludes that nothing happened has read the wrong number — which is why the keyword exists and why it is documented here rather than discovered.

Examples

julia
fresh = fresh_paste(0.30)                     # the volume reference
r     = VanGenuchten(; a = 37.5479e6, m = 1 / 2.1684)   # a measured isotherm
q = Ref(Float64[])
eq, cert = equilibrate_certified(
    state; constraint = CapillaryWater(r; reference = fresh), parameters = q,
)

See also: WaterRetention, kelvin_activity, powers_alpha_max, porosity.

ChemistryLab.EquilibriumConstraint Type
julia
EquilibriumConstraint

What is held fixed while the Gibbs energy is minimized. One of FixedTP, FixedEnthalpy, Adiabatic, FixedVolume or SealedVolume.

ChemistryLab.FixedActivity Type
julia
FixedActivity(species, a; titrant = species)

Activity of species prescribed, the system open to titrant.

titrant is the substance whose amount adjusts to hold the activity — by default the species itself. Its amount is an unknown of the same system, so the answer is both the equilibrium composition and how much titrant it took to get there, which is what a titration measures.

The element balance is then A x − A[:, titrant] q = b: the system is open to that one substance and closed to every other.

ChemistryLab.FixedEh Type
julia
FixedEh(Eh; couple = "SO4-2" => "HS-", titrant = "O2@")

Equilibrium at a prescribed redox potential in volts. Converted to a pe through the Nernst relation at the state's own temperature, and handled by FixedpE from there — so everything said there applies, including the warning about prescribing a potential on a closed system.

Eh is a quantity with the dimensions of an electric potential, so that the unit is stated rather than assumed.

ChemistryLab.FixedEnthalpy Type
julia
FixedEnthalpy(H)

Enthalpy prescribed, temperature unknown. H is a quantity with the dimensions of energy, referred to the same element basis as enthalpy — so in practice it comes from enthalpy(state) of another state on the same budget, never from a table.

ChemistryLab.FixedTP Type
julia
FixedTP()

Temperature and pressure both held at the values the state carries. The default, and the only constraint the interior-point back ends can honor.

ChemistryLab.FixedVolume Type
julia
FixedVolume(V)

Volume prescribed, pressure unknown. V is a quantity with the dimensions of volume, and only the species carrying a molar volume contribute — see volume.

Needs a system whose volume depends on pressure

A condensed system's does not. The molar volumes of water and of the minerals in the shipped databases are exactly pressure-independent, so prescribing the volume of a paste or of an aqueous solution is an equation the pressure cannot satisfy, and the constructor refuses it rather than diverging — the error names the lever it measured. Declare a gas phase to make the constraint meaningful. To report the volume change of a sealed specimen at fixed pressure, which is what a hydrating binder needs, use porosity(state, reference) instead.

ChemistryLab.FixedpE Type
julia
FixedpE(pe; couple = "SO4-2" => "HS-", titrant = "O2@")

Equilibrium at a prescribed electron activity, pe = −log₁₀ a(e⁻), with the oxidation state of the system free to follow.

The vehicle is the implicit titrant of FixedpH, with one difference that is forced on it: there is no electron species, so nothing can have its activity prescribed directly. What is prescribed instead is the half-reaction of couple, whose log K fixes the electron activity once the other members' activities are known:

summed over the couple with the products positive and the reactants negative, n being the number of electrons. That is one linear equation in the log-activities, which is exactly what the titrant mechanism solves.

titrant is the substance the system may draw on to reach the prescribed potential — O2@ by default, so that oxidizing means adding oxygen. It must be a species of the system.

When to prescribe a potential, and when not to

The default in this package is to conserve the oxidation state and let the potential come out as a result, which is what a sealed paste does and what pe then reports. Prescribe one only when the system is open to a redox buffer that is genuinely imposed from outside — a measured Eh, a controlled atmosphere, an electrode. A prescribed potential on a closed paste is a statement about a system that is not the one being modeled.

See also: FixedEh, pe, half_reaction.

ChemistryLab.FixedpH Type
julia
FixedpH(pH; titrant = "H+")

pH prescribed, the system open to titrant. Shorthand for FixedActivity("H+", 10^-pH; titrant = titrant).

titrant = "H+" means acid is added or removed. To titrate with a base instead, name it: FixedpH(12.5; titrant = "OH-"). The answer carries the amount, so the titrant is not a numerical device — it is the reagent consumed.

Which pH is prescribed

This constrains −log₁₀ a(H⁺) in the activity model's own convention, and it holds it exactly: the residual is that quantity, and it comes out at the target to five decimals. pH reports something slightly different — −log₁₀ c(H⁺) with the concentration taken over the computed liquid volume — and with DiluteSolutionModel the two differ by about 0.0013, which is log₁₀ of the ρ ≈ 1 kg/L approximation that model makes when it converts molality to concentration. The gap grows with ionic strength; use HKFActivityModel or DaviesActivityModel where it matters.

ChemistryLab.SaturatedCuring Type
julia
SaturatedCuring(V_ref; titrant = "H2O@")
SaturatedCuring(; reference, titrant = "H2O@")

A specimen cured under water: free to draw in whatever the reaction's own volume loss empties, so the pore space never desiccates.

This is the mirror image of CapillaryWater, and the pair is the two boundary conditions a paste can be cured under. Sealed, the volume that chemical shrinkage empties becomes gas-filled porosity, the saturation falls, the water activity falls with it and the reaction slows — which is what CapillaryWater expresses. Immersed, that volume is refilled from the bath, the specimen stays saturated, and the water activity stays the composition's own.

V_ref is the volume the specimen occupies, normally the fresh paste's: pass the fresh state as reference and it is taken from it. The constraint holds

the total volume of the system, solids and solution together — external dimensions unchanged, the deficit made up by water from outside. Every species carrying no standard molar volume would contribute zero to that sum in silence, so the constructor refuses such a system and names them.

What it does to the solve

One unknown, q[1], the amount of water imbibed; one column −A[:, titrant] in the conservation rows, so the system is open to water and closed to everything else; one equation, the volume closure above — which is linear in the composition, unlike the retention law of CapillaryWater.

The answer includes the chemical shrinkage

q[1] is not a numerical device. It is the water the specimen took up, which is exactly what a chemical-shrinkage measurement reports, and it comes back through the parameters keyword the way a titrant amount does:

julia
q = Ref(Float64[])
fresh = fresh_paste(0.40)
eq, cert = equilibrate_certified(
    state; constraint = SaturatedCuring(; reference = fresh), parameters = q,
)
only(q[])                       # moles of water drawn in

This is a volume condition, not an activity condition

The tempting way to write "cured under water" is FixedActivity("H2O@", 1.0) — and it is wrong. A cement pore solution has a water activity near 0.98 from its dissolved salts alone, so prescribing 1 would draw water in until the solution was dilute enough to reach it, which never happens: the constraint would imbibe without bound. What a bath fixes is not the activity inside the specimen, it is the availability: the pore space stays full. That is a volume statement, and this is it.

What the certificate proves under it

The volume closure is linear and the conservation rows stay affine, so the problem remains a convex minimization on an affine set and cert.optimal keeps its usual meaning — a global minimum of G subject to the budget and the closure. This is a stronger guarantee than CapillaryWater gives, whose composition-dependent activity shift is not derived from a convex G.

See also: CapillaryWater, powers_alpha_max, porosity.

ChemistryLab.SealedVolume Type
julia
SealedVolume()

Volume held at the value the initial state occupies, pressure unknown: a rigid sealed vessel. Equivalent to FixedVolume(volume(state).total).

ChemistryLab._constraint_blocks Method
julia
_constraint_blocks(constraint, des, state, p, n0) -> NamedTuple

The parameter block a constraint contributes: nq, the callbacks gq, hq and cq, the starting guess q0, the difference-step scale qscale, and apply, which writes the parameter that was found back onto the resulting state.

nq = 0 for FixedTP, in which case the solver takes its plain path.

ChemistryLab._molar_volumes Method
julia
_molar_volumes(system, T, P) -> Vector{Float64}

Standard molar volume of every species at T, P, in m³/mol, and zero for a species that has none.

The zero is why CapillaryWater refuses a system with a missing molar volume before it starts: such a species contributes nothing to the volume balance, silently, and the saturation the whole coupling rests on would be wrong with nothing to show for it.

ChemistryLab._pressure_lever Method
julia
_pressure_lever(system, n, T, P) -> Real

Relative sensitivity of the system's volume to pressure, (∂V/∂P)·P/V, by a central difference over one percent of P.

A volume constraint prescribes V and solves for P, so it needs this to be non-negligible. It usually is not for a condensed system: in the databases shipped here the molar volumes of water and of the minerals do not depend on pressure at all — V⁰(1 bar) = V⁰(100 bar) exactly for H2O@ and Cal — and only the partial molar volumes of a few aqueous ions vary, OH- by 8 % over 100 bar. A kilogram of water with 0.25 mol of ions then has a lever of about 1e-6, meaning some 10 000 bar to change the volume by one percent. Newton on that residual takes an enormous step and the pressure leaves the domain of the equation of state.

Which is the physics, not a solver defect: the volume of an incompressible condensed system is fixed by its composition, and pressure has no purchase on it. A gas phase gives it one.

ChemistryLab._redox_terms Method
julia
_redox_terms(des, state, couple) -> (terms, logK, n_e)

The half-reaction of couple as a list of index => coefficient over the system's species, with products positive and reactants negative, together with its log₁₀ K at the state's temperature and its electron count.

The electron is dropped from terms: it has no index, which is the whole reason a redox constraint cannot be written as a prescribed activity.

ChemistryLab._species_index Method
julia
_species_index(des, s) -> Int

Position of a species in the solver's system, by symbol or by object.

ChemistryLab._titrant_blocks Method
julia
_titrant_blocks(des, state, p, terms, i_titrant, ln_target)

The implicit-titrant blocks for a prescribed linear combination of log-activities, Σᵢ νᵢ ln aᵢ = ln_target, where terms is a vector of index => coefficient.

One species with coefficient 1 is a prescribed activity, which is what FixedActivity and FixedpH need. A whole half-reaction is what a prescribed redox potential needs, since no electron activity can be read off a species that does not exist — see FixedpE.

The vehicle is the same either way: one unknown (the titrant amount), one column −A[:, titrant] added to the conservation rows, and one equation.

ChemistryLab._titrant_blocks Method
julia
_titrant_blocks(des, state, p, i_species, i_titrant, ln_a_target)

The parameter block of a prescribed chemical potential: one unknown, the titrant amount; one linear column, its formula; one equation, ln aᵢ = ln a_target.

ChemistryLab._total_enthalpy Method
julia
_total_enthalpy(system, n, T, P)

Σᵢ nᵢ ΔₐH⁰ᵢ(T, P) in joules, as a bare number, over the species that carry an enthalpy of formation. The element type follows n, T and P, so this differentiates.

ChemistryLab._total_volume Method
julia
_total_volume(system, n, T, P)

Σᵢ nᵢ V⁰ᵢ(T, P) in cubic meters, as a bare number, over the species that carry a molar volume.

Water retention ​

The relation between how much water a pore space still holds and how tightly it holds it — the constitutive input CapillaryWater needs, and the source of the humidity PoreHumidity hands to a rate law. It is measured, not assumed, so nothing here carries a default value.

ChemistryLab.FunctionRetention Type
julia
FunctionRetention(f)

A bare function as a retention law: f(S) is the water activity at degree of saturation S.

This is the escape hatch for a curve that is neither a measured table nor van Genuchten — a closed form from a pore-size distribution, an isotherm fitted with someone else's expression, or a constant, which is how the tutorial's negative control imposes one humidity and watches nothing happen.

Because the value returned is an activity, water_activity ignores its V_m and T: there is no Kelvin conversion to make. A law stated as a capillary pressure belongs in VanGenuchten instead, and passing a pressure here would be read as an activity of several million.

Nothing validates f

TabulatedRetention checks its table and VanGenuchten checks its parameters, both in inner constructors. f is opaque, so it is trusted: it should return a value in (0, 1], and it should not increase as the pore space dries. CapillaryWater does test it at S = 1 and refuses a law that is already out of range there, which catches a sign error or a percentage but not a curve that misbehaves in the middle.

CapillaryWater and PoreHumidity wrap a bare function in this type themselves, so CapillaryWater(S -> ...; reference = fresh) needs no explicit construction.

Examples

julia
held = FunctionRetention(_ -> 0.90)               # a paste held at RH 90 %
water_activity(held, 0.5; V_m = 1.807e-5, T = 298.15)   # 0.9, the arguments unused

See also: WaterRetention, TabulatedRetention, VanGenuchten.

ChemistryLab.TabulatedRetention Type
julia
TabulatedRetention(; S, a_w) -> TabulatedRetention

A measured desorption isotherm, as two vectors: degrees of saturation and the water activities (relative humidities) at which they were observed.

Interpolation is linear in ln a_w against S, which is the variable the capillary coupling actually uses and the one in which a Kelvin curve is closest to straight. Both keywords are mandatory.

S must be strictly increasing and a_w non-decreasing with it — a retention curve that falls as the pore space fills is not a retention curve, and an inner constructor refuses it rather than interpolating nonsense. Outside the tabulated range the value is clamped to the nearest end and a warning is issued once: an isotherm measured between S = 0.3 and S = 1 says nothing about S = 0.05, and extrapolating a logarithm there produces confident absurdity.

Examples

julia
# RH measured at four saturations on a hardened paste
r = TabulatedRetention(; S = [0.30, 0.50, 0.75, 1.00], a_w = [0.44, 0.66, 0.85, 1.00])
r(0.60)      # interpolated

See also: WaterRetention, VanGenuchten.

ChemistryLab.VanGenuchten Type
julia
VanGenuchten(; a, m) -> VanGenuchten

The van Genuchten water-retention form, as a capillary pressure

julia
p_c(S) = a * (S^(-1/m) - 1)^(1 - m)

converted to a water activity by the Kelvin relation a_w = exp(-p_c V_m / (R T)), with V_m and T supplied by the constraint that uses it.

Both parameters are mandatory keywords: a is a capillary pressure scale (Pa, or a pressure quantity) and m is the dimensionless shape exponent, 0 < m < 1. There are no defaults — the pair belongs to a fitted material, not to the package.

Published parameters, and the two conventions

(Baroghel-Bouny et al., 1999) fit exactly this expression to measured water-vapor desorption isotherms, and write it with b = 1/m:

julia
p_c(S) = a (S^(-b) - 1)^(1 - 1/b)

so a b from that literature becomes m = 1/b here. Their Table 5, for materials whose mixes and porosities are in their Tables 1 and 4:

mixmaterialW/Ca (MPa)bm = 1/b
COcement paste0.3437.54792.16840.46117
CHpaste, 10 % silica fume0.1996.28371.95400.51177
BOconcrete0.4818.62372.27480.43960
BHconcrete, 10 % silica fume0.2646.93642.06010.48541

Examples

julia
# The ordinary cement paste of Baroghel-Bouny et al. (1999), their mix CO
r = VanGenuchten(; a = 37.5479e6, m = 1 / 2.1684)

A pair belongs to a material, not to this package

The table above is quoted with its source so it can be checked, not so it can be copied blindly. Those four fits are one cement, one curing history and one age; m and a both move with w/c, with silica fume and with the aggregate. Fit your own isotherm where you have one, and cite what you used.

See also: WaterRetention, TabulatedRetention.

ChemistryLab.WaterRetention Type
julia
abstract type WaterRetention

A relation between the degree of saturation of the pore space and the activity of the water left in it.

The interface is one method: a subtype is callable, r(S) -> a_w, mapping a degree of saturation S ∈ [0, 1] to a water activity in (0, 1]. S = 1 is a saturated pore space, where the water is held by nothing and a_w = 1.

Two ways to obtain one:

And one way to build your own: a pore size distribution plus kelvin_activity.

No defaults, deliberately

Every named law takes its parameters as keyword arguments without default values, so omitting one raises UndefKeywordError on the spot. A retention curve is a measurement on a specific material at a specific age; a default would be a number invented on the user's behalf and then silently believed.

See also: CapillaryWater, kelvin_activity.

ChemistryLab.capillary_pressure Method
julia
capillary_pressure(r::VanGenuchten, S) -> Real

The capillary pressure in pascals at saturation S, before the Kelvin relation turns it into an activity. Exposed because it is the form the retention literature tabulates and fits.

ChemistryLab.kelvin_activity Method
julia
kelvin_activity(r; γ, V_m, T) -> Real

Activity of water held in a pore of radius r by its own meniscus, a_w = exp(-2 γ V_m / (r R T)).

This is the whole of the physics the capillary coupling adds. Water in a fine pore is at a lower chemical potential than bulk water of the same composition, so hydrates that consume water become less stable, and below some pore size the reaction stops — which is what self-desiccation is.

Every parameter is a keyword without a default: the surface tension γ belongs to the liquid and the temperature, and the molar volume V_m to the liquid, so neither is the package's to assume.

Examples

julia
julia> using DynamicQuantities

julia> a = kelvin_activity(4.2e-9u"m"; γ = 0.072u"N/m", V_m = 1.8e-5u"m^3/mol", T = 298.15u"K");

julia> round(a; digits = 2)      # the self-desiccation plateau of a sealed paste
0.78

See also: kelvin_radius, WaterRetention.

ChemistryLab.kelvin_radius Method
julia
kelvin_radius(a_w; γ, V_m, T) -> Real

The pore radius, in meters, whose meniscus holds water at activity a_w. The inverse of kelvin_activity.

Useful for reading a retention curve back as a pore size, which is the form in which the literature usually discusses it: a_w = 0.78 — the relative humidity a sealed high-performance paste settles at — is 4.2 nm, the gel-pore scale.

Examples

julia
julia> using DynamicQuantities

julia> r = kelvin_radius(0.78; γ = 0.072u"N/m", V_m = 1.8e-5u"m^3/mol", T = 298.15u"K");

julia> round(r * 1e9; digits = 1)      # nanometers
4.2

See also: kelvin_activity.

ChemistryLab.water_activity Method
julia
water_activity(r::WaterRetention, S; V_m, T) -> Real

Activity of the water a retention law leaves at degree of saturation S.

For a law that is already stated as an activity — TabulatedRetention, or a measured isotherm wrapped in FunctionRetention — this is the law itself and V_m and T are ignored. For a law stated as a capillary pressure, VanGenuchten, it is the Kelvin relation a_w = exp(-p_c V_m / RT), which is why the molar volume and the temperature have to be supplied: a pressure curve carries no temperature, and turning it into an activity does.

Both are keywords without defaults, for the reason given in WaterRetention.

Examples

julia
co = VanGenuchten(; a = 37.5479e6, m = 1 / 2.1684)
water_activity(co, 0.786; V_m = 1.807e-5, T = 298.15)     # ≈ 0.80

See also: capillary_pressure, kelvin_activity, CapillaryWater.

Problem and solver ​

ChemistryLab.EquilibriumProblem Type
julia
EquilibriumProblem

Definition of a chemical equilibrium problem.

Fields

  • b: conservation vector (elemental abundances).

  • A: stoichiometric matrix (conservation matrix).

  • μ: chemical potential function μ(n, p).

  • u0: initial guess for species amounts.

  • p: coefficients for the potential function (default: nothing).

  • lb: lower bounds for species amounts.

  • ub: upper bounds for species amounts.

The problem solves for the species distribution that minimizes the Gibbs energy subject to mass conservation constraints A * n = b.

ChemistryLab.EquilibriumProblem Method
julia
EquilibriumProblem(A, μ, u0; b=A*u0, p=nothing, lb=fill(Tu(1e-16), length(u0)), ub=maximum(abs.(A))/minimum(abs.(A[.!iszero.(A)]))*sum(u0)*one.(u0))

Construct an EquilibriumProblem with the given stoichiometric matrix A, chemical potential function μ, and initial guess u0.

Arguments

  • A: stoichiometric matrix (conservation matrix).

  • μ: chemical potential function μ(n, p).

  • u0: initial guess for species amounts.

  • b: conservation vector (elemental abundances). Defaults to A * u0.

  • p: coefficients for the potential function. Defaults to nothing.

  • lb: lower bounds for species amounts. Defaults to fill(Tu(1e-16), length(u0)).

  • ub: upper bounds for species amounts. Defaults to maximum(abs.(A))/minimum(abs.(A[.!iszero.(A)]))*sum(u0)*one.(u0).

Returns

An EquilibriumProblem instance.

ChemistryLab._concrete_float Method
julia
_concrete_float(x) -> AbstractArray

Narrow a numeric array to a concrete floating-point element type, returning an already-concrete one untouched.

ChemicalSystem stores its stoichiometry as Matrix{Real} whenever integer and rational coefficients coexist, which a cement's does — C3AFS0.84H4.32 and its kind. That is right for the chemistry and wrong for the solver: an abstract element type boxes every entry and turns mul!(res, A, x) into the generic fallback with a dynamic dispatch per element, on a product evaluated at every objective and constraint call. SciMLBase says so out loud, warning that "arrays or dicts to store parameters of different types can hurt performance" as soon as such an array reaches the problem's parameters, and the warning is correct.

Applies to the default b = A * u0 as well, which inherits the same abstract element type from A.

Floating point rather than the promoted exact type (Rational{Int} here) because the numeric pipeline downstream is Float64 throughout — u0, the bounds and DualEquilibriumSolver.A, which has always converted — so keeping rationals would only pay for a rational-to-float conversion at every entry of every product. The conversion is lossy for a non-dyadic rational: 2//5 becomes 0.4, which is not equal to it. That is the same rounding the rest of the solver already applies, and the exact matrix is untouched in system.SM.A, where the stoichiometry belongs.

A caller who passes a concrete array keeps exactly what they passed, identically: an exact Matrix{Rational{Int}}, a Matrix{Int}, or a Matrix{<:Dual} for someone differentiating through it.

ChemistryLab.EXACT_HESSIAN Constant
julia
EXACT_HESSIAN

Whether to hand the back-end the exact Gibbs Hessian diagonal ∂μ/∂n computed by ForwardDiff, instead of letting it approximate one. Default false.

It is off by default only because it currently trades one defect for another, and both are measured:

pure watercalcite + CO₂
false (back-end approximation)[H⁺]/[OH⁻] = 3.78 ✗worst ×19.8
true (exact ∂μ/∂n)[H⁺]/[OH⁻] = 1.000003 ✓worst ×3751 ✗

Turn it on for an aqueous-only system, where it makes the water autoprotolysis come out right. Leave it off when a pure phase is present: the exact curvature of such a phase is zero, the interior-point iteration then stalls essentially at its starting point, and the whole speciation is wrong.

That stall is the open problem. It is not a matter of iterating longer (the answer is identical at 300 and at 200 000 iterations) nor of the near-singular Hessian entry (capping the inverse curvature over five orders of magnitude changes nothing, because the solve stops before the barrier has decayed). The C++ Optima this back-end is ported from offers a Nullspace linear solver in addition to the Rangespace one implemented here — the latter carries the warning that it suits diagonal Hessians only, and it is the one that inverts H. Porting the nullspace path, which never inverts H, is the identified next step.

ChemistryLab.NONCONVERGED Constant
julia
NONCONVERGED :: Ref{Int}

Running count of equilibrium solves that returned a non-success retcode.

With STRICT_CONVERGENCE[] = false — the default — a non-converged solve is a @warn at maxlog = 1 and its result is used anyway. Over the thousands of steps of a coupled kinetics run that is one warning for an arbitrary number of bad speciations, and it never reached the failure count reported by integrate, which only saw solves that actually threw.

Reset it with ChemistryLab.NONCONVERGED[] = 0 before a run and read it after; integrate does exactly that and reports the total.

ChemistryLab.NULLSPACE_STEP Constant
julia
NULLSPACE_STEP

Whether the back-end computes its Newton step by the nullspace method rather than the Schur complement. Default true.

The Schur complement forms S = A H⁻¹ Aᵀ, so it needs H invertible — and a pure phase has unit activity, hence exactly zero curvature. The nullspace method writes dn = dnₚ + Z dz with Z a basis of null(A) and solves (Zᵀ H Z) dz = −Zᵀ(ex + H dnₚ), in which H appears only as a product. This is the route the C++ Optima takes by default, and its Rangespace counterpart — the Schur complement — is documented there as suitable for invertible diagonal Hessians only.

It is what makes the water autoprotolysis come out right: pure water gives [H⁺]/[OH⁻] = 1.0 and pKw = 13.9994, against 3.78 and 13.9897 through the Schur complement, with no change to mixed solid/aqueous systems.

ChemistryLab.STRICT_CONVERGENCE Constant
julia
STRICT_CONVERGENCE

Whether a non-converged equilibrium solve raises (true) or warns (false, the default).

The default is not strict, and deliberately so: the back-end's convergence flag is unreliable in both directions on these problems. It reports MaxIters on points that are numerically excellent — pure water comes back flagged while giving [H⁺]/[OH⁻] = 1.000003 — and reports success on points that are not the minimum. Raising on the flag alone would reject good answers and would still miss the bad ones, so it is offered as an opt-in for callers who want the strictest possible reading.

ChemistryLab._DEFAULT_SOLVER_FACTORY Constant
julia
_DEFAULT_SOLVER_FACTORY

Internal Ref{Union{Nothing, Function}} — populated by extension __init__ functions to register a default solver factory.

  • OptimizationIpoptExt.__init__: registers only if nothing is set (low priority).

  • OptimaSolverExt.__init__: always overrides (high priority).

Result: OptimaSolver wins whenever loaded, regardless of load order.

ChemistryLab._DUAL_AVAILABLE Constant
julia
_DUAL_AVAILABLE

Whether the KKT solver and its certificate are loaded. Set by OptimaSolverExt.__init__, and false otherwise: the dual Newton lives in OptimaSolver, so with only the Ipopt extension loaded there is no certifying route and equilibrate must take its plain path rather than raise.

Registering a back end is not the same question. OptimizationIpopt registers a factory, which makes it a usable STARTING POINT for the certified route, but it cannot certify anything by itself.

ChemistryLab._EXPLORING_STARTS Constant
julia
_EXPLORING_STARTS

Set while a multi-start route is computing or trying starting points, so the diagnostics of a candidate are not reported as diagnostics of the answer.

A start that does not converge is ordinary and expected: equilibrate_certified runs every back end from several compositions precisely because none of them works on every problem, and it keeps whichever answer the certificate proves. Left unguarded, a call that ends optimal = true still printed "returned MaxIters" and "did not certify optimality" from candidates along the way, which reads as a failed solve and is not one.

The verdict on the answer is untouched: equilibrate_certified warns or raises on its own certificate after the search, and verbose = true still reports every rejected start.

ChemistryLab._SOLVER_FACTORIES Constant
julia
_SOLVER_FACTORIES

Every back end an extension has registered, in load order. equilibrate uses the first as its default; equilibrate_certified uses all of them as starting points, because neither back end dominates the other — see solve_certified.

ChemistryLab.EquilibriumSolver Type
julia
struct EquilibriumSolver{F<:Function, S, V<:Val}

Encapsulates all fixed ingredients of a chemical equilibrium calculation: the potential function, the SciML solver, and the variable space.

Construct once, call repeatedly with different ChemicalState inputs.

Fields

  • μ: chemical potential closure μ(n, p) -> Vector{Float64}.

  • solver: any Optimization.jl-compatible solver (e.g. IpoptOptimizer()).

  • variable_space: variable space — Val(:linear) or Val(:log).

  • kwargs: solver keyword arguments forwarded to solve.

Examples

julia
julia> cs = ChemicalSystem([
           Species("H2O"; aggregate_state=AS_AQUEOUS, class=SC_AQSOLVENT),
           Species("H+";  aggregate_state=AS_AQUEOUS, class=SC_AQSOLUTE),
           Species("OH-"; aggregate_state=AS_AQUEOUS, class=SC_AQSOLUTE),
       ]);

julia> solver = EquilibriumSolver(cs, DiluteSolutionModel(), IpoptOptimizer());

julia> solver isa EquilibriumSolver
true
ChemistryLab.EquilibriumSolver Method
julia
EquilibriumSolver(cs, model, solver; variable_space=Val(:linear), kwargs...)

Construct an EquilibriumSolver from a ChemicalSystem, an activity model, and a SciML solver.

The potential function is built once at construction time from cs and model. Repeated calls to solve with different ChemicalState inputs reuse it.

Arguments

  • cs: the ChemicalSystem defining species and conservation matrix.

  • model: an AbstractActivityModel (e.g. DiluteSolutionModel()).

  • solver: any Optimization.jl solver.

  • variable_space: Val(:linear) (default) or Val(:log).

  • kwargs...: forwarded to the underlying solve call (tolerances, verbosity...).

ChemistryLab._attach_sensitivity Method
julia
_attach_sensitivity(state, nstar, μ, ϵ; b = nothing) -> ChemicalState

Differentiate the equilibrium map implicitly at a composition already found, and return it carrying the dual parts.

Split out of _solve_dual so the same derivative can be attached to an answer obtained by any route — in particular to a certified one, which is solved in real arithmetic by construction. Pushing duals through the iteration itself would be wrong anyway: an active set has a discrete component, so the map b ↦ n*(b) is smooth only piecewise, and the derivative belongs at the solution with the active set frozen.

ChemistryLab._build_n0 Method
julia
_build_n0(state::ChemicalState) -> Vector

Extract the dimensionless mole vector from a ChemicalState. Type is inferred from the state — compatible with ForwardDiff dual numbers.

ChemistryLab._build_params Method
julia
_build_params(state::ChemicalState; ϵ=1e-16) -> NamedTuple

Extract dimensionless parameters from a ChemicalState. ΔₐG⁰overRT is evaluated at the current T and P of the state. Units are stripped — compatible with ForwardDiff dual numbers.

Returned fields

  • ΔₐG⁰overRT: vector of standard Gibbs energies divided by RT (dimensionless).

  • T: temperature in K (plain number, Dual-safe).

  • P: pressure in Pa (plain number, Dual-safe).

  • ϵ: regularization floor (default 1e-16).

T and P are included so that temperature-dependent activity models (e.g. HKFActivityModel with temperature_dependent=true) can recompute their parameters inside the potential closure.

ChemistryLab._check_converged Method
julia
_check_converged(sol, what) -> sol

Return sol, raising or warning according to STRICT_CONVERGENCE when the optimizer did not converge. Before this existed, neither extension looked at the return code and a non-converged iterate was written into the state as though it were the equilibrium.

ChemistryLab._equilibrium_sensitivity Method
julia
_equilibrium_sensitivity(A, H, gθ, bdot, nstar; maxpin = 8) -> Vector

Sensitivity of an equilibrium composition, from the optimality conditions.

Equilibrium is the Gibbs minimization of (Leal et al., 2017), the problem Reaktoro solves: min G(n) subject to A n = b, n >= 0, whose first-order conditions are grad G(n) - A' y - z = 0, A n = b, n_i z_i = 0, with y the element potentials and z >= 0 the stability multipliers.

Differentiating them gives, on the FREE set (species present, z_i = 0),

julia
    | H   A' | | ndot |   | -dgradG |
    | A   0  | | ydot | = |   bdot  |

with ndot = 0 on the complement and H = grad^2 G at the solution. One factorization serves every partial derivative, and the answer is exact — no finite difference, no step size.

The complementarity block is not optional, and dropping it fails loudly rather than subtly: on calcite + CO2 in water with a gas phase declared, the unreduced system puts the whole perturbation into the absent gas species — n(CO2,g) = 5.8e-9, held at its bound — returning ndot = e_CO2(g), which satisfies A ndot = bdot to 4e-16 and means nothing.

No back-end returns z, so the active set is recovered here: a species that is negligible on the scale of the system yet takes a leading share of the response is pinned, and the system re-solved. Each pass pins at least one species, so the loop terminates.

H is singular by construction, and correctly so — a pure phase has unit activity, hence a zero row. The saddle-point form handles that; any method inverting H does not.

ChemistryLab._exploring_starts Method
julia
_exploring_starts(f)

Run f with _EXPLORING_STARTS set, restoring it afterwards. Nested calls are safe: the previous value is saved rather than assumed false.

ChemistryLab._primal Method
julia
_primal(state::ChemicalState) -> ChemicalState

The same state with every dual number replaced by its value.

ChemistryLab._solve_dual Method
julia
_solve_dual(esolver, state, ϵ) -> ChemicalState

Equilibrium of a composition carrying dual numbers.

No optimization solver is asked to iterate on dual numbers — most cannot, and Ipopt never will, being a C library. The equilibrium is solved once at the primal values and the sensitivities come from _equilibrium_sensitivity, the implicit-function-theorem route on the optimality conditions.

Called from the back-end solve methods, which dispatch on the solver type; making this a method of solve dispatching on the state would be ambiguous with them.

ChemistryLab.activity_model Method
julia
activity_model(solver::EquilibriumSolver) -> AbstractActivityModel

The activity model the solver's potential function was built from.

ChemistryLab.equilibrate Method
julia
equilibrate(state::ChemicalState, solver; model=..., variable_space=..., ϵ=...) -> ChemicalState
equilibrate(state::ChemicalState; kwargs...) -> ChemicalState

Compute the chemical equilibrium state by minimizing the Gibbs free energy.

Two-argument form (solver explicit, always available once an extension is loaded):

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

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

One-argument form — solves by every available route and returns the answer optimality_certificate proves globally optimal:

julia
state_eq = equilibrate(state)                  # certified
state_eq = equilibrate(state; certify = false) # single back end, as before

The certified route is the default because a single back end is not reliable here: measured on calcite dissolving in water, the interior point returns a composition whose charge balance is wrong in the second digit (3 %), because the fraction-to-boundary rule caps its step and the residual stops moving. The dual Newton gets that case to 1e-12 but fails to admit a supersaturated phase on a low-water cement. Offering both and keeping a proved answer certifies all ten cases of the reference battery; either alone certifies at most nine.

Use equilibrate_certified when the certificate itself is wanted, and certify = false for the old single-back-end behavior. certify = true has no effect on a system without an aqueous phase or without H2O@, where the dual route does not apply.

When both extensions are loaded, OptimaSolverExt provides the default single back end.

Arguments

  • state: initial ChemicalState — defines the system, T, P, and composition.

  • solver: any SciML-compatible solver (e.g. IpoptOptimizer(), OptimaOptimizer()).

  • model: activity model (default: DiluteSolutionModel()).

  • variable_space: Val(:linear) (default) or Val(:log).

  • ϵ: regularization floor for mole amounts (default: 1e-16).

  • kwargs...: forwarded to the underlying solver.

ChemistryLab.register_solver_factory! Method
julia
register_solver_factory!(f)

Called from an extension's __init__ to make its back end available to equilibrate_certified. Idempotent.