ChemistryLab.jl

Appearance

Equilibrium ​

Activity models ​

ChemistryLab.ELECTROSTATIC_STIFFNESS_LIMIT Constant
julia
ELECTROSTATIC_STIFFNESS_LIMIT

The largest electrostatic_stiffness at which eliminating the surface potential — rather than carrying it as an unknown — still converges here, 5.0.

It is a measured threshold, not a derived one, and both sides of it were measured. A constant capacitance of 3 F/m² on ferrihydrite gives 4.7 and solves to a stationarity of 2e-16; 2 F/m² gives 7.0 and does not solve. Across eighteen diffuse-layer points spanning three ionic strengths and six pH values, everything at or below 3.4 converged and everything at or above 6.2 did not, with nothing in between. The bracket is therefore [3.4, 6.2] and 5.0 sits in it.

See electrostatic_stiffness for what the number means and why a threshold exists at all.

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

Fallback effective electrostatic radii r_e,j [Å] indexed by formal charge, for species absent from Table 3 of Helgeson et al. (1981): Table H.1-1 of the TOUGHREACT V2 user's guide (Xu et al., 2012), read from data/literature/Xu2012.json.

Turned into an ion size by HKFActivityModel exactly as the radii of REJ_HKF are, when that table has no entry for the ion.

See also: REJ_HKF, HKFActivityModel.

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

Effective electrostatic radii r_e,j [Å] of aqueous ions from Helgeson et al. (1981), Table 3, read from data/literature/Helgeson1981.json.

Keys are PHREEQC-format formula strings (e.g. "Na+", "Ca+2", "SO4-2"). A radius is not an ion size: HKFActivityModel turns it into one with Eq. (125) of Helgeson et al. (1981), for the electrolyte the ion forms with the NaCl background (see its docstring), at priority 3 of its lookup: model.å > 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, solid-solution end-members, and species occupying a surface site. 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 build _MixingTerms(cs) once and call _mixing_lna!, or the members of a solid solution, or of a surface site family, silently get ln a = 0 — unit activity, which is a plausible number and a wrong one;

  • concentration_scale has no fallback: a model declares which scale its solute standard state uses, and nothing guesses it.

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 the partner of these by the Gibbs–Duhem relation,

the integral of I f'(I) for the ionic term f(I) = √I/(1+√I) − b I.

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 comes from its own activity coefficients

Every ion carries the same function of I, times zᵢ², so the water activity that the Gibbs–Duhem relation pairs with it has the closed form above, and with it the ions and the solvent derive from one excess Gibbs energy. The neutral solutes are the exception: their salting-out term bₙ I has no partner in the coefficients of the ions, so with neutral species present and bₙ ≠ 0 the model is not the gradient of a Gibbs energy, and its certificate says so (:self_consistent). With bₙ = 0, or without neutral solutes, it is.

Inputs, their defaults, and where each default comes from

fielddefaultunitprovenance
A0.5114(kg/mol)^½the LLNL aqueous model at 25 °C (Parkhurst and Appelo, 2013; p. 118), 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/molPHREEQC's coefficient for an uncharged species, log γ = b I (Parkhurst and Appelo, 2013; p. 201)
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 for the water activity, which follows from them.

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, an activity is a concentration, and the water's activity is the one the solutes' imply.

phaselawexpression
solventideal diluteln a = −M_w Σ mᵢ, the partner of ln mᵢ by Gibbs–Duhem
aqueous solutesHenryln a = ln(cᵢ/c°), c° = 1 mol/L
pure crystals—ln a = 0
gasideal mixtureln a = ln xᵢ + ln(P/P°)
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.

The water's activity is not a separate choice. With the solutes at ln mᵢ, the Gibbs–Duhem relation Σⱼ nⱼ ∂ln aⱼ/∂nᵢ = 0 leaves the solvent exactly one partner, ln a_w = −M_w Σ mᵢ = −Σ nᵢ/n_w, so that the chemical potentials are the gradient of one Gibbs energy, homogeneous and convex, and every back end minimizing it solves the same equations. Raoult's mole fraction x_w agrees with it to first order, the two differing by (M_w Σm)²/2, but breaks that relation by 1 − x_w.

Its water activity counts solute molecules, and that is the real limit

ln a_w = −M_w Σ m counts the solutes and knows nothing about what they are. A cement pore solution at Σm = 0.5 mol/kg gets a_w = 0.991 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 solutes' log-activities 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)^½the LLNL aqueous model at 25 °C as Parkhurst and Appelo (2013) tabulate it (p. 118) — 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. Helgeson et al. (1981), Table 1, gives 0.5091 at 25 °C from the water properties of their time
B0.3288Å⁻¹(kg/mol)^½the same table; Helgeson et al. (1981) give 0.3283
Ḃ0.041kg/molthe same table: the B-dot of the LLNL model at 25 °C. What the term is is set out in Anderson and Crerar (1993) §17.7.1
Kₙ0.1kg/molthe coefficient b of log γ = b I that PHREEQC gives an uncharged species with no parameters of its own (Parkhurst and Appelo, 2013; p. 201); overridden per species by sp[:Kₙ], which is how CO₂(aq) gets its own
å_default3.72Åthe distance of closest approach of NaCl, (Helgeson et al., 1981) Table 2; a last resort, reached only for a charge no table covers (`
ånothingÅone common ion size 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., 2012; Table H.1-1); each is turned into an ion size as set out under Ion-size lookup below.

Every default has a source, and a source is not a fit

Each default is read from data/literature/, with the table and the page it was read from. That makes them traceable, not tailored: Ḃ is the value of a NaCl-dominated solution, å_default the ion size of NaCl, and a pore solution that is neither calls for its own, a keyword away. A and B are additionally derived here from the water model, and agree with the table to 1e-3.

What is approximated in the water activity

The γᵢ use a per-ion size; the osmotic coefficient uses one charge-weighted mean size å_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. In NaCl they do not — both ions have the size of the salt, 3.72 Å — but in a solution of CaCl₂, Na₂SO₄ or a cement pore solution they do, and the inconsistency is then measurable.

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 size at all.

Ion-size lookup

The ion-size parameter åᵢ is resolved in order:

  1. model.å — a common ion size 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 ion size in the species properties dict.

  3. The effective electrostatic radius r_e,i of REJ_HKF ((Helgeson et al., 1981; Table 3), keyed by PHREEQC formula), turned into an ion size by (Helgeson et al., 1981; Eq. 125) for the electrolyte the ion forms with the NaCl background — åᵢ = 2 (r_e,i + r_e,c |zᵢ|) / (1 + |zᵢ|), r_e,c that of Cl⁻ for a cation and of Na⁺ for an anion, as TOUGHREACT does (Xu et al., 2012; Eqs. H.4–H.5). It reproduces (Helgeson et al., 1981; Table 2): 3.72 Å for NaCl, 4.32 for CaCl₂, 4.65 for Na₂SO₄.

  4. REJ_CHARGE_DEFAULT — a radius by formal charge, turned into an ion size the same way.

  5. model.å_default — an ion size, reached only for a charge no table covers, i.e. |z| ≥ 5. It is not a way to impose a common ion size; 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. The provenance of every default is set out on HKFActivityModel.

å imposes one common ion size 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 (Kulik et al., 2013), 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 ion sizes from the radii of REJ_HKF, the B-dot of the LLNL
# model at 25 °C.
HKFActivityModel()

# One common ion size for every ion: that of NaCl.
HKFActivityModel(å = literature_value("Helgeson1981", "nacl_distance_of_closest_approach"))

# The model Cemdata18 prescribes, its Eq. (C.1): one ion size and one B-dot for
# every ion, the same B-dot for the neutral species.
cemdata18_activity_model(:KOH)

# The Debye-Hückel limiting law with the B-dot that one particular GEM-Selektor
# run of a Portland cement was configured with, Ḃ_gems ≈ 0.0976, recovered from
# the coefficients it printed (test/reference/gems_cemdata18_portland.json). That
# run is not the configuration Cemdata18 prescribes.
HKFActivityModel(å = 0.0, Ḃ = Ḃ_gems, Kₙ = 0.0)
ChemistryLab.LLNLActivityModel Type
julia
struct LLNLActivityModel{T<:Real} <: AbstractActivityModel

The aqueous model of Lawrence Livermore National Laboratory, as a PHREEQC database that carries LLNL_AQUEOUS_MODEL_PARAMETERS defines it (llnl.dat): the B-dot equation with Debye–Hückel parameters tabulated in temperature, and a formula of its own for aqueous carbon dioxide (Parkhurst and Appelo, 2013; p. 120).

Formulas

For an ion given a size å by -llnl_gamma:

A, B and Ḃ interpolated linearly in temperature between the points of the grid of the database. For a neutral species marked -co2_llnl_gamma, the temperature T in kelvin and C, F, G, E, H the coefficients of -co2_coefs:

Every other species has γ = 1: a neutral species given a size, and an ion given none. A species given (å, b) by -gamma has log₁₀ γ = b I, without the Debye–Hückel term. These are what PHREEQC 3.7.3 computes with such a database, measured against it: the manual does not state them.

The solvent follows Raoult's law, ln a_w = ln x_w, or, with water = :phreeqc, PHREEQC's a_w = 1 − 0.017 Σᵢ mᵢ (Parkhurst and Appelo, 1999; Eq. 25). The model is not the gradient of one Gibbs energy, and a certificate says so (:self_consistent).

Fields

  • temperatures: the grid, in °C; A, B, Bdot the parameters at its points. A temperature outside the grid is refused, as PHREEQC refuses it.

  • co2: (C, F, G, E, H).

  • sizes: species symbol ⇒ å, in Å.

  • co2_species: the symbols of the species that take the formula of carbon dioxide.

  • gamma: species symbol ⇒ (å, b).

  • water: :raoult or :phreeqc.

Built from a database by database_activity_model.

ChemistryLab.LLNLActivityModel Method
julia
LLNLActivityModel(; temperatures, A, B, Bdot, co2, sizes, co2_species = String[],
                    gamma = Dict(), water = :raoult)

Construct an LLNLActivityModel.

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

The activity model of a PHREEQC database: the WATEQ equation of Truesdell and Jones (1974) for every species its -gamma option gives parameters to, and PHREEQC's defaults for the others (Parkhurst and Appelo, 2013; description of SOLUTION_SPECIES).

Formulas

For a species with parameters (å, b):

which for an uncharged species is b_i I. For a species without parameters, the Davies equation for an ion and 0.1 I for a neutral species, the two defaults PHREEQC applies:

The solvent follows Raoult's law, ln a_w = ln x_w, or, with water = :phreeqc, the approximation PHREEQC applies to every database without a Pitzer or SIT block (Parkhurst and Appelo, 1999; Eq. 25),

the sum over the solutes. Neither is the partner the Gibbs–Duhem relation gives the per-species parameters, which have no closed-form one, so the model is not the gradient of one Gibbs energy, and a certificate says so (:self_consistent).

Fields

  • A, B: the Debye–Hückel constants, at 25 °C unless temperature_dependent, in which case they are computed from the water properties at every call.

  • parameters: species symbol ⇒ (å, b), as phreeqc_gamma_parameters reads them from a database file.

  • water: :raoult or :phreeqc, the activity of the solvent.

Example

julia
params = phreeqc_gamma_parameters(datapath("phreeqc.dat"))
model = TruesdellJonesActivityModel(; parameters = params)
ChemistryLab.TruesdellJonesActivityModel Method
julia
TruesdellJonesActivityModel(; parameters, A=0.5114, B=0.3288, temperature_dependent=false, water=:raoult)

Construct a TruesdellJonesActivityModel from a dictionary of species parameters, typically phreeqc_gamma_parameters(path).

ChemistryLab._DebyeHuckelLna Type
julia
_DebyeHuckelLna{K, G, X, M}

The log activities lna(n, p) of a system whose aqueous phase follows a model of the Debye–Hückel family: HKFActivityModel, DaviesActivityModel, TruesdellJonesActivityModel and LLNLActivityModel. The frame is the same for all four,

with nᵢ floored at the activity floor, and the solvent row the Gibbs–Duhem partner of the solutes' terms. The kernel K holds what is the model's own: its coefficients A and B (_coefficients), the terms of its solutes (_solute_lna!) and its solvent row (_solvent_row). The gas phase and the mixing phases follow (_gas_lna!, _mixing_lna!). M is the number type of the data the kernel and the mixing models capture.

ChemistryLab._MixingTerms Type
julia
_MixingTerms(cs::ChemicalSystem)

The solid solutions and the site families of cs, prepared once when an activity closure is built: their member indices and mixing models, the denticity and the charge of each site member, and what a diffuse layer needs to evaluate the ionic strength and the charge of each support. _mixing_lna! applies them.

Every activity model holds one and calls that function, so that a new model cannot leave either mixing term out. The fields are concretely typed through the type parameters, the closure capturing the struct as it is.

ChemistryLab._aqueous_ionic_strength Method
julia
_aqueous_ionic_strength(_n, idx_ions, z_ions, idx_solvent, M_w) -> Real

The molal ionic strength I = ½ Σ zᵢ² mᵢ of the aqueous solution, from the amounts _n and the solvent mass _n[idx_solvent] · M_w.

Written here, beside the site mixing, rather than reused from an activity model: three of the four models compute an ionic strength of their own and one does not, so taking it from them would make the surface answer depend on which aqueous model is in force through a path no one asked for. This is one definition, evaluated identically under all four.

The element type follows _n, so the whole path differentiates.

ChemistryLab._davies_osmotic_W Method
julia
_davies_osmotic_W(A, b, I, ϵ) -> Real

The excess part of the Davies solvent row, ln a_w = −M_w Σm + M_w W(I), with

f being the ionic term of the solutes, regularized as they regularize it, so that the Gibbs–Duhem relation between the solutes and the solvent holds to rounding. The antiderivative of s/(1+s) in u = t + ϵ is u − 2√u + 2 ln(1 + √u).

ChemistryLab._electrostatic_ln_a Method
julia
_electrostatic_ln_a(model, k, z, n, I, T) -> Real

The electrical work of adding one mole of member k to a charged surface, in units of RT, given the members' formal charges z and amounts n, the ionic strength I of the bulk solution, and the temperature T.

Zero for a model that does not describe a surface potential, which is every one of them but ConstantCapacitance, DiffuseLayer and ChargePlanes, whose work is a sum over its planes. The first two write it as z_k ψ̃ with ψ̃ = FΨ/RT; they differ only in the closure that gives Ψ from the surface charge density σ = F Σ_j z_j n_j / 𝒜.

  • Constant capacitance: σ = CΨ inverts to ψ̃ = F² Σ_j z_j n_j / (C 𝒜 RT), which is linear in the composition and ignores the solution entirely.

  • Diffuse layer: Gouy-Chapman's σ = κ√I sinh(ψ̃/2) inverts to ψ̃ = 2 asinh(σ / (κ√I)), with κ = √(8 ε_r ε₀ R T ρ). It is monotone in Ψ, which is why it inverts at all, and asinh keeps it smooth.

Neither needs an unknown of its own. What the second costs instead is written in is_gradient_consistent: I depends on the aqueous composition, and the aqueous activities do not depend in return on the surface, so this term is not the gradient of a Gibbs energy.

I is floored at eps() so that an unscreened surface gives a large finite potential rather than a NaN. That floor is a numerical guard, not a physical one: Gouy-Chapman genuinely diverges as the solution runs out of ions, and a diffuse layer in pure water is the model outside its range of validity.

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._gas_lna! Function
julia
_gas_lna!(out, n, idx_gas, p) -> out

The log activities of an ideal gas mixture, ln aᵢ = ln xᵢ + ln(P/P°): the fugacity xᵢ P over the standard-state pressure. Every activity model calls this one for its gas species, so that the pressure enters each of them in the same way, and it is what gives a gas ∂μᵢ/∂P = RT/P, the molar volume of an ideal gas.

ChemistryLab._gibbs_duhem_defect Method
julia
_gibbs_duhem_defect(J, n) -> (defect, column)

The worst relative violation of the Gibbs–Duhem relation Σᵢ nᵢ Jᵢⱼ = 0 over the columns of a Jacobian J of log activities (or of chemical potentials) at n, each column measured against Σᵢ |nᵢ Jᵢⱼ|, and the column where it is worst (0 when every column vanishes).

ChemistryLab._hkf_ion_size Method
julia
_hkf_ion_size(r_e, z) -> Real

The ion-size parameter of an ion of effective electrostatic radius r_e [Å] and charge z in a NaCl-dominated solution: Eq. (125) of Helgeson et al. (1981), å_k = 2 Σⱼ νⱼ,ₖ r_e,j / νₖ, for the electrolyte the ion forms with the counter-ion of the background, Cl⁻ for a cation and Na⁺ for an anion,

julia
å = 2 (r_e + r_e,c |z|) / (1 + |z|),

which is how TOUGHREACT forms it (Xu et al., 2012; Eqs. H.4–H.5). A radius is not an ion size: Na⁺ has r_e = 1.91 Å, and NaCl å = 3.72 Å, the value Helgeson et al. (1981), Table 2 tabulate.

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

Compute the σ function used in the osmotic coefficient formula (Helgeson et al., 1981; Eqs. 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._jacobian_asymmetry Method
julia
_jacobian_asymmetry(J) -> (worst, (i, j))

The largest relative asymmetry |Jᵢⱼ − Jⱼᵢ| / max(|Jᵢⱼ|, |Jⱼᵢ|) over the pairs of a square Jacobian, and the pair it is reached on ((0, 0) when every entry is negligible). Second derivatives of one energy commute, so a Jacobian of the log activities that is not symmetric says they are not the gradient of one Gibbs energy.

ChemistryLab._ln_pressure_ratio Method
julia
_ln_pressure_ratio(p) -> Real

ln(P/P°) for the parameters p an activity closure is called with, P in pascals and P° = P_STANDARD. Parameters built without a pressure, as a few internal probes build them, are at P°.

ChemistryLab._log10γ_ion Method
julia
_log10γ_ion(model, z, å, I, sqrtI, A, B) -> Real

log₁₀ γ for an ion of charge z and size å at ionic strength I, the compiled half of the pair log10_gamma_expression writes out.

Takes sqrtI rather than computing it, because the caller already regularizes it as sqrt(I + ϵ) to keep a ForwardDiff.Dual finite at I = 0, and takes A and B rather than reading them off model, because they depend on temperature and pressure through the water model and the caller has them.

This is the function every activity evaluation calls, and the accessor calls it too, so the two cannot drift. DaviesActivityModel fixes one ion size for everything and ignores å and B.

ChemistryLab._mixing_lna! Method
julia
_mixing_lna!(out, _n, mix::_MixingTerms, p, ϵ) -> out

Write into out the log activities of the members of every solid solution and every surface site family that mix describes, at the amounts _n and the parameters p. The temperature is read from p (298.15 K when it carries none), and a surface potential carried as an unknown of the solve arrives the way an adiabatic temperature does, through p.ψ_site.

ChemistryLab._site_excess_ln_gamma Method
julia
_site_excess_ln_gamma(model, k, x, T) -> Real

The departure from ideality of the k-th member of a site family, given the site fractions x of the whole family. Zero for IdealSiteMixing and for every model that adds no excess term to the site fractions.

The twin of _excess_ln_gamma for solid solutions, and deliberately a separate generic: a site fraction and a mole fraction obey different closures — a site family's total is pinned by a conservation row, a solid solution's is free — so a model written for one is not automatically valid for the other.

ChemistryLab._site_mixing_lna! Method
julia
_site_mixing_lna!(out, _n, site_groups, site_models, site_denticity,
                  site_charges, I, T, ϵ, ψ_site)

Fill out[i] with ln a_i = ln x_i + ln γ_i for every species occupying a surface site, x_i being its fraction of its family's site budget.

site_groups[k] lists the members of the k-th family, the free site first, site_models[k] is how they mix, and site_denticity[k] how many of the family's sites each member occupies.

Which fraction the activity is — mole or equivalent — is the convention's, and it enters through _site_weight. For a family of monodentate species the two coincide, which is why an oxide surface never has to choose.

Why this is not _solid_solution_lna! under another name

The arithmetic is the same and the meaning is not. A solid solution's total is an unknown the minimization is free to move; a site family's total is fixed by a conservation row, because sites are neither created nor destroyed while the support is fixed. The free site is a member like any other, and it is what makes saturation happen: as it runs out, ln x of every occupied state rises, and the cost of binding one more molecule with it. Langmuir is the consequence, not the premise.

ϵ regularizes the logarithm at exhaustion, exactly as it does for a solid solution, and the element type follows _n, so the whole path differentiates.

ChemistryLab._site_weight Method
julia
_site_weight(model, d::Int) -> Int

What one molecule of a family member counts as, in the fraction the convention takes for its activity.

One particle for a mole fraction — IdealSiteMixing and VanselowMixing — and d units of charge for an equivalent fraction, GainesThomasMixing, d being how many of the family's pseudo-elements the species carries.

That single number is the whole difference between the two exchange conventions. It is one for every monodentate species, which is why the two agree on a homovalent exchange and part company on a heterovalent one.

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

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). g is the vector of ΔₐG⁰/RT of every species, which only a model whose activities depend on the energies of its members reads (CompoundEnergyModel).

ForwardDiff-compatible.

ChemistryLab._ss_log_activities! Method
julia
_ss_log_activities!(out, grp, x, model, T, ϵ) -> out

Write ln aₖ into out[grp[k]] for the members of one solid solution at the mole fractions x. The default is ln xₖ + ln γₖ, with ln γₖ from _excess_ln_gamma; a model whose activities are not naturally written that way adds a method here rather than a special case in the loop above.

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

Return a closure lna(n, p) -> Vector computing log-activities for the model of Davies (1962). 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 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ᵢ) + ln(P/P°).

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.activity_model Method
julia
activity_model(cs::ChemicalSystem, model::LLNLActivityModel) -> Function

Return a closure lna(n, p) -> Vector computing log-activities for the LLNL aqueous model, at the temperature p.T (T_STANDARD when p carries none).

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

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

Return a closure lna(n, p) -> Vector computing log-activities for the Truesdell–Jones model of a PHREEQC database. The parameters of each species are resolved once, at construction.

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 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.cemdata18_activity_model Function
julia
cemdata18_activity_model(electrolyte = :KOH; temperature_dependent = true)
    -> HKFActivityModel

The extended Debye–Hückel model with which Cemdata18 was derived and which it prescribes (Lothenbach et al., 2019; Appendix C, Eq. C.1):

with one ion-size parameter a and one b_γ common to every ion, set by the electrolyte that dominates the pore solution:

electrolytea (Å)b_γ at 25 °C (kg/mol)
:KOH3.670.123
:NaOH3.310.098

The values are read from data/literature/Lothenbach2019.json. In terms of HKFActivityModel this is å = a, Ḃ = Kₙ = b_γ.

  • A and B depend on temperature and pressure, as Lothenbach et al. (2019) state; with temperature_dependent = true they are computed from the water properties at every call. The paper gives b_γ at 25 °C only, and it is kept at that value: Helgeson et al. (1981) tabulate b_γ against temperature for HCl, LiCl, NaCl, KCl, MgCl₂ and CaCl₂ (Table 26, p. 1457), and not for KOH or NaOH, the two electrolytes the paper gives parameters for.

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

  • A certified answer computed with it is :self_consistent, not a :global_minimum. The term b_γ I gives every ion i the contribution b_γ ∂I/∂n_j ∝ b_γ z_j² to the derivative of its log activity with respect to the amount of j, and for two ions of different charges that is not symmetric in i and j, so no Gibbs energy has these activities as its gradient.

Choose the electrolyte from the pore solution: in a Portland cement potassium usually dominates.

Examples

julia
model = cemdata18_activity_model(:KOH)
activity_model_range(model)      # 1.0 mol/kg
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.diffuse_layer_potential Method
julia
diffuse_layer_potential(model, z, n, I, T) -> Real

The dimensionless surface potential ψ̃ = FΨ/RT that Gouy-Chapman's closure gives for a surface carrying Σ_j z_j n_j moles of charge on model.area, screened by a solution of ionic strength I:

ρ = 1000 kg/m³ turns a molal ionic strength into the volumetric one the relation is written for. At 25 °C in water κ = 0.117215, against the 0.1174 PHREEQC carries in its source.

This is the eliminated form of the potential — a function of the composition. It is what the activity model uses when the solve does not carry Ψ as an unknown, and what the closure equation is written against when it does, which is why it lives in one place.

ChemistryLab.electrostatic_stiffness Method
julia
electrostatic_stiffness(model, z, n, I, T) -> Real

How strongly the surface potential reacts to the composition that raises it: the dimensionless sensitivity N |∂ψ̃/∂n|, with N = Σ_k n_k the family's site budget. Zero for a model that carries no potential.

What it predicts, and why a threshold exists

Both electrostatic models here are written as activity coefficients, which the dual solver reaches by successive substitution: it forms a composition, reads the activities it implies, and forms a composition again. That iteration contracts only while the activity moves less than the composition does — which is exactly this number being smaller than one, up to the damping the solver already applies. Above it the iteration walks off, and what comes back is not a second root but a point that violates mass action outright.

It is not a property of the physics. The feedback is always negative — charging the surface always opposes further charging, so the equilibrium is unique and stable. What fails is the elimination, and the answer to that is to carry Ψ as an unknown with its own equation, which linearizes the coupling instead of iterating it; the dual solver does so for a DiffuseLayer by default (surface_potential = :auto).

The shape of it

The first does not depend on the composition at all, so a capacitance can be judged once. The second does, through u, and that is what makes a diffuse layer hard in the middle and easy at the edges: far from the point of zero charge asinh flattens, the surface screens itself, and the same system that diverges at pH 7 converges at pH 4. It also scales as 1/√I, so a dilute background is the difficult case, not the easy one.

Use

Evaluate it at a composition you believe — a solve without the electrostatics, or another code's answer — and compare against ELECTROSTATIC_STIFFNESS_LIMIT before trusting a solve. It is a forecast, not a verdict: the verdict is the stationarity residual of optimality_certificate, which separates the two regimes by fourteen orders of magnitude and never lets a diverged answer pass silently.

ChemistryLab.electrostatic_stiffness Method
julia
electrostatic_stiffness(state; model = DiluteSolutionModel(), T = 298.15) -> Vector

electrostatic_stiffness for every site family of state's system, in the order site_families returns them, evaluated at the composition state holds and at the ionic strength that composition implies.

This is the form to reach for: it takes the state a solve produced, or the one a solve without electrostatics produced, and says whether the elimination will hold there.

ChemistryLab.excess_ln_gamma_expression Function
julia
excess_ln_gamma_expression(model, k::Int, n::Int = 2) -> Num

The symbolic ln γₖ of a solid-solution model, in the variables x₁ … xₙ and T, with the model's own parameter values substituted.

Why a package that already computes this also writes it down

thermo_factories.jl keeps a thermodynamic model in two forms: a symbolic expression, which makes the law and the unit of every parameter visible so that a wrong entry is caught by reading rather than by a wrong result, and a compiled function, type-stable and fit for a solver's inner loop. Solid solutions had only the second, with the formula transcribed into a docstring beside it – two copies of one law, free to drift.

This is the first form. The second stays exactly where it was, on the hot path, untouched: _excess_ln_gamma is still what every activity evaluation calls, and this function is never in that path. What the pair buys is a test that they agree, which is what makes the written formula true by construction rather than by proofreading.

Example

julia
julia> using Symbolics

julia> expr = excess_ln_gamma_expression(RedlichKisterModel(a0 = 4000.0), 1);

julia> expr isa Num
true

See also: _excess_ln_gamma, RedlichKisterModel, RegularSolutionModel.

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 (Haar et al., 1984)) and water_electro_props_jn (Johnson-Norton dielectric constant (Johnson and Norton, 1991)).

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

Examples

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

julia> isapprox(p.A, HKFActivityModel().A; rtol=1e-3)   # the tabulated default at 25 °C
true

julia> isapprox(p.B, HKFActivityModel().B; rtol=1e-3)
true
ChemistryLab.log10_gamma_expression Method
julia
log10_gamma_expression(model, z, å) -> Num

The symbolic log₁₀ γ of an ion of charge z and size å, in the ionic strength I, with the model's own parameters substituted.

Why a package that computes this also writes it down

thermo_factories.jl keeps a thermodynamic model in two forms: a symbolic expression, which makes the law and its parameters visible so that a wrong entry is caught by reading, and a compiled function for a solver's inner loop. The activity kernels had only the second.

This is the first form, and the compiled path is untouched: _log10γ_ion is still what every activity evaluation calls, and this is never in that path. A test asserts the two agree, which is what makes a formula quoted in prose true by construction rather than by proofreading.

What it is for, beyond inspection

Activity models shows that a model derived from a single excess Gibbs energy must satisfy ∂ln aᵢ/∂nⱼ = ∂ln aⱼ/∂nᵢ, and that for these kernels

so symmetry demands that f'(I; zᵢ, åᵢ)/zᵢ² not depend on i. With this function that f is an expression one can differentiate symbolically and look at, rather than a residual one can only measure.

Example

julia
julia> using Symbolics

julia> expr = log10_gamma_expression(DaviesActivityModel(), -2, 0.0);

julia> expr isa Num
true

See also: _log10γ_ion, HKFActivityModel, DaviesActivityModel.

ChemistryLab.site_gradient_asymmetry Method
julia
site_gradient_asymmetry(cs, model, n; T = 298.15, ϵ = 1.0e-30) -> NamedTuple

How far the activity map of cs under model is from being the gradient of a potential, at the composition n. Returns (; worst, pair, gd):

  • worst: the largest relative asymmetry over pairs, |J_ij − J_ji| / max(|J_ij|, |J_ji|) with J_ij = ∂ln a_i/∂n_j;

  • pair: the species symbols of the pair that attains it;

  • gd: the largest relative Gibbs-Duhem residual, |Σ_i n_i J_ij| / Σ_i |n_i J_ij|.

A twice-differentiable Gibbs energy has a symmetric Hessian, so an activity map that comes from one is symmetric here and both numbers are zero to rounding. This is the measurement behind is_gradient_consistent: rather than assert that DiffuseLayer breaks the certificate's premise, run it and read the number.

Why the scale is per pair and not the whole Jacobian

Because the Jacobian holds self-derivatives of trace species of order 1/n, normalizing by its largest entry makes every real inconsistency look like zero. The equivalent measurement in test/activities.jl records that this was done once and reported 1e-8 for a model off by 20 %. The per-pair scale is what makes the number mean what it says, and the two measurements deliberately agree.

Differentiates the whole activity closure with ForwardDiff, so it costs about length(n) evaluations of it. It is a diagnostic, not part of a solve.

Real gases ​

The equation of state of Peng and Robinson (1976) for a gas phase: the critical constants a gas carries, its fugacity coefficients and the compressibility factor of its phase. See Real gases.

ChemistryLab._PengRobinsonMixing Type
julia
_PengRobinsonMixing

The Peng-Robinson description of the gas phase of a system: per gas, in the order of idx_gas, the critical temperature (K) and pressure (Pa) and the coefficient , and the matrix of binary interaction parameters.

ChemistryLab._gas_mixing Method
julia
_gas_mixing(cs) -> Union{Nothing, _PengRobinsonMixing}

The mixing model of the gas phase of cs: nothing for an ideal mixture, a Peng-Robinson description when its gases carry critical constants. A phase in which some gases carry them and others do not is refused: there is no mixing rule between an equation of state and its absence.

ChemistryLab._pr_ln_phi Method
julia
_pr_ln_phi(mix, y, T, P) -> (lnφ, Z)

The fugacity coefficients of the gases at the mole fractions y, the temperature T (K) and the pressure P (Pa), and the compressibility factor of the phase, by the van der Waals mixing rules,    and   (Peng and Robinson (1976), their Eq. 19 and the mixture form of their fugacity coefficient).

ChemistryLab._pr_phase_volume Method
julia
_pr_phase_volume(mix, n, T, P) -> V

The volume of the gas phase,   in m³, for the amounts n of its gases (mol, in the order of the mixture), at T (K) and P (Pa): the volume of the ideal gas scaled by the compressibility factor of the mixture.

ChemistryLab._pr_residual_gibbs Method
julia
_pr_residual_gibbs(mix, n, T, P) -> G_res/(RT)

The residual Gibbs energy of the gas phase over , for the amounts n (in the order of the mixture), at T (K) and P (Pa). Its derivatives with respect to the amounts are the logarithms of the fugacity coefficients, which the tests check _pr_ln_phi against.

ChemistryLab._pr_roots Method
julia
_pr_roots(a, b) -> Vector{Float64}

The real roots above b of the cubic in Z at the reduced parameters a, b (plain numbers), in closed form, each sharpened by one Newton step: one root in a single-phase region, three where the equation has a vapor and a liquid root.

ChemistryLab._pr_z Method
julia
_pr_z(A, B) -> Z

The compressibility factor: of the real roots of the cubic above B, the one of least Gibbs energy, which is the stable phase where the equation has a vapor and a liquid root. The roots are found on the values, then lifted into the dual numbers of A and B by two Newton steps on the cubic: the value does not move and the derivatives of the root are exact to second order, as the implicit-function theorem gives them.

ChemistryLab.compressibility_factor Method
julia
compressibility_factor(state::ChemicalState) -> Float64

The compressibility factor   of the gas phase of state: 1 for an ideal gas phase, the root of the cubic of Peng and Robinson (1976) of least Gibbs energy when the gases carry critical constants (peng_robinson). The gas phase occupies .

See also: fugacity_coefficients.

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

The fugacity coefficient of every gas of state, at its temperature, pressure and gas composition: 1 in an ideal gas phase, the value of the equation of state of Peng and Robinson (1976) when the gases carry critical constants (peng_robinson). The fugacity of a gas is .

julia
co2 = peng_robinson(db["CO2"])
st = ChemicalState(ChemicalSystem([co2], ["CO2"]); T = 298.15u"K", P = 50u"bar", n = [1.0u"mol"])
fugacity_coefficients(st)["CO2"]             # 0.741

See also: compressibility_factor.

ChemistryLab.peng_robinson Method
julia
peng_robinson(s::Species; T_c, P_c, ω, kij = Pair{Symbol, Float64}[]) -> Species
peng_robinson(s::Species; kij = Pair{Symbol, Float64}[]) -> Species

The gas s described by the equation of state of Peng and Robinson (1976), through its critical temperature T_c, critical pressure P_c and acentric factor ω, with the binary interaction parameters kij, a list of Symbol(partner) => k, zero for a partner it does not name. Without the three constants, they are those PHREEQC ships in phreeqc.dat for the gas of the same symbol (Parkhurst and Appelo (2013)), which literature_table("ParkhurstAppelo2013", "gas_critical_constants") lists.

A gas phase is described as a whole: either every gas of a system carries critical constants and the phase follows the equation of state, or none does and the phase is ideal. The activity of a gas is then    , its fugacity coefficient given by the equation of state, and the standard state stays the ideal gas at .

T_c and P_c are quantities, or plain numbers in K and Pa. The copy shares the thermodynamic functions of s; only its properties gain the constants, in SI units, under :T_c (K), :P_c (Pa), :ω and :kij.

Examples

julia
co2 = peng_robinson(db["CO2"])                       # PHREEQC's constants
co2 = peng_robinson(db["CO2"]; T_c = 304.128u"K", P_c = 73.773u"bar", ω = 0.2239)
n2 = peng_robinson(db["N2"]; kij = [:CO2 => -0.02])

Specific ion Interaction Theory (SIT) ​

Debye-Hückel with a deviation term that is a sum over ion pairs rather than one coefficient times the ionic strength — so a sodium sulfate solution and a sodium chloride one of the same ionic strength get different activity coefficients. See Activity models §6a.

Nothing here carries a default compilation: ε is caller input, and an unlisted pair is taken as zero by the literature's convention, which missing_epsilon_pairs reports rather than leaves implied. The coefficients carry their own provenance, so a borrowed value never becomes indistinguishable from a measured one.

ChemistryLab.SITActivityModel Type
julia
struct SITActivityModel{T<:Real, P} <: AbstractActivityModel

The Specific ion Interaction Theory (SIT): the Debye-Hückel screening term and one interaction coefficient per pair of ions of opposite charge.

with the sum over ions k of charge opposite to i, m_k in molality and b = 1.5 kg^{1/2} mol^{-1/2}.

What it is, against its neighbors in this package

Read the three together and the family is one idea with three answers to the same question — what to do beyond Debye-Hückel's limiting law.

ion sizecomposition-dependent term
DaviesActivityModelnone, absorbed into 1 + √I−A z² b I, one parameter for every ion
HKFActivityModelper ion, åḂ I, one parameter for every ion
SITnone, absorbed into 1 + 1.5√IΣ ε(i,k) m_k, one parameter per pair
PitzerActivityModelnonebinary and ternary, several per pair

So SIT sits between Davies and Pitzer exactly where its parameter count does. It is the model several thermodynamic databases are calibrated in, so that a log K taken from one of them and used under Davies is not the constant that was fitted.

Validity, stated because it is narrower than the equation looks

SIT is a truncation: the pairwise term is the first order of a virial expansion, so at high ionic strength the terms it drops stop being small. No source transcribed in this package states the ionic strength where that happens, so activity_model_range returns nothing for it rather than a number of its own. It also says nothing about ion pairs that are better described as species — where a complex forms, SIT expects it in the speciation and not in ε.

Fields

  • A: the Debye-Hückel slope, 0.509 at 25 °C, or recomputed per temperature when temperature_dependent.

  • b: the denominator coefficient, 1.5 by convention. It is a convention, not a fit: changing it makes every published ε inconsistent with the model that produced it, so it is a field only to make that visible.

  • parameters: a SITParameters compilation. Empty by default, which reduces the model to its Debye-Hückel term — a legitimate limiting case, and never silently a full SIT calculation.

  • temperature_dependent: whether A follows temperature.

  • water: :raoult or :osmotic.

Water

The solvent activity is Raoult's mole fraction, as in TruesdellJonesActivityModel, unless water = :osmotic, which takes it from the osmotic coefficient of the same model, as PHREEQC does with a database that carries a SIT block:

the sum over the solutes and over the pairs of the compilation: the partner the Gibbs–Duhem relation gives the coefficients of the solutes. Raoult's mole fraction is a departure from a full SIT treatment; it does not enter a comparison made at prescribed proton activity.

See also: SITParameters, build_sit_parameters, missing_epsilon_pairs.

ChemistryLab.SITActivityModel Method
julia
SITActivityModel(; A = 0.509, b = 1.5, parameters = SITParameters(),
                   temperature_dependent = false, water = :raoult) -> SITActivityModel

Build a SITActivityModel.

ChemistryLab.SITParameters Type
julia
struct SITParameters{T<:Real}

A compilation of Traced coefficients, keyed by the unordered pair of species symbols, with the source of the compilation as a whole.

SIT sums ε over ions of opposite charge only, so the pair is unordered and the lookup is symmetric: sit_epsilon(p, "Na+", "Cl-") and sit_epsilon(p, "Cl-", "Na+") are the same entry.

Fields

  • epsilon: (a, b) => Traced, with a ≤ b as strings so one pair has one key.

  • source: what the compilation is, as a whole — a database name and version, a paper, or a note saying it was assembled here.

ChemistryLab.SITParameters Method
julia
SITParameters(pairs; source = "unstated", kind = PROV_UNSTATED) -> SITParameters

Build a compilation from (a, b) => ε pairs, where ε is a number or a Traced. A bare number takes kind and source; a Traced keeps its own, so a borrowed coefficient does not silently acquire the compilation's standing.

Keys are normalized so that a pair given either way round is one entry, and giving the same pair twice with different values is an error rather than a last-one-wins.

ChemistryLab._sit_epsilon_matrix Method
julia
_sit_epsilon_matrix(cs, model) -> Matrix

The ε a system needs, as a dense matrix over its species indices: zero for a like-charged pair, for two neutrals, and for a pair the compilation does not carry.

A neutral solute against an ion is kept when — and only when — the compilation supplies a coefficient for it. SIT leaves a neutral ideal by default, and that default is the absence of a coefficient rather than a rule applied here: excluding neutrals outright made the closure's neutral branch unable to do anything at all, so a coefficient somebody wrote down was dropped on the way in and the calculation silently ignored it.

Its element type is the compilation's, so a Dual coefficient stays a Dual all the way to the activity — which is what differentiating with respect to an interaction coefficient requires.

Built once when the closure is, so the inner loop is a matrix read rather than a dictionary lookup.

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

A closure lna(n, p) -> Vector of log-activities under SIT.

Unlike DaviesActivityModel and HKFActivityModel, SIT has no _log10γ_ion: its correction depends on the whole composition through Σ ε(i,k) m_k, not on the ion's charge and the ionic strength alone, so there is no per-ion formula to share with activity_coefficients. That is the same reason PitzerActivityModel has none, and it is a property of the model rather than an omission.

AD-compatible: every computation inside accepts ForwardDiff.Dual.

ChemistryLab.missing_epsilon_pairs Method
julia
missing_epsilon_pairs(cs::ChemicalSystem, model::SITActivityModel) -> Vector

Every oppositely charged pair of aqueous solutes in cs for which model has no ε, as (a, b) symbol tuples.

The SIT literature takes an unlisted coefficient as zero. That convention is usually harmless and is never a measurement, so this reports what a calculation is resting on rather than leaving it to be assumed. An empty result means every pair the system can form is carried by the compilation.

Neutral solutes are not reported. A neutral with no coefficient is ideal, which is SIT's own answer and not a gap in the compilation; listing every neutral-ion pair as missing would bury the pairs that are.

ChemistryLab.sit_epsilon Method
julia
sit_epsilon(p::SITParameters, a, b) -> Traced or nothing

The coefficient for the unordered pair (a, b), or nothing when the compilation does not carry it.

nothing and zero are different answers, and the distinction is the point: the SIT literature's convention is to take an unlisted ε as zero, which is a statement about what is customary and not about what was measured. See missing_epsilon_pairs for asking a system which of its pairs are resting on that convention.

ChemistryLab.build_sit_parameters Method
julia
build_sit_parameters(path; source = nothing) -> SITParameters

Read the SIT / -epsilon block of a PHREEQC-format database into a SITParameters compilation.

Reading rather than transcribing is deliberate. A compilation of several hundred interaction coefficients copied by hand is a transcription, which is the error class every generator in test/reference/ exists to remove, and this package has already found two standard energies that had drifted that way.

No database is stored in this package. The compilation PHREEQC distributes in sit.dat, the ANDRA/RWM ThermoChimie database, is obtained by datapath("sit.dat") from the PHREEQC distribution; any other file the caller has is read the same way. source defaults to the file's name and a truncated SHA-256 of its contents, so a result can always say which compilation it came from.

Format

julia
SIT
-epsilon
    Na+     Cl-     0.03
    ...

Three whitespace-separated fields per line: two species and the coefficient in kg/mol. The block ends at the next keyword. The species are named as in the database and keyed by the rule of every reader: SO4-- is SO4-2, a neutral species gets @. A line of another shape is not read, and a warning lists it.

Example

julia
p = build_sit_parameters("/path/to/sit.dat")
model = SITActivityModel(; parameters = p)
missing_epsilon_pairs(cs, model)     # what the compilation does not cover

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, etheta = true,
                      missing_pairs = :refuse, A = 0.5114)

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.

missing_pairs = :zero takes the coefficients of such a pair as zero instead, which is what PHREEQC does with a database whose PITZER block leaves pairs out (pitzer.dat describes neither H⁺–OH⁻ nor Ca²⁺–CO₃²⁻): it is the convention of the database, and database_activity_model builds the model of such a database with it.

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.

The higher-order electrostatic terms

Two ions of like sign and unsymmetrical charge (Na⁺ with Ca²⁺, Cl⁻ with SO₄²⁻) interact through besides , the terms of Pitzer (1975): with  ,

and its derivative with respect to . is evaluated by the Chebyshev approximation of Harvie, as PHREEQC evaluates it (Plummer et al., 1988). They vanish for a symmetrical pair, so a single electrolyte is unaffected, and they enter and the osmotic coefficient from the same excess Gibbs energy, so that the Gibbs–Duhem relation keeps holding exactly. etheta = false leaves them out, for a set whose were fitted without them (Pitzer and Mayorga's, for instance): and are fitted together, and a set is used with the convention it was fitted with.

Temperature

A is the Debye–Hückel constant of the decimal logarithm, from which  . It is used as given unless temperature_dependent = true, which takes from the water model at the temperature of the state, and evaluates the temperature terms of the coefficients when the set carries them (see PitzerParameters). A published set fitted at one temperature carries none, and is then used at its own values whatever the temperature, which is what the set can support.

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}(),
                   temperature = Dict{Symbol, Any}())

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.

Temperature

The tables give each coefficient at 25 °C, A₀. temperature may give, for any of them, the further coefficients of the form PHREEQC uses (Parkhurst and Appelo, 2013),

keyed by the table and then by the entry's key: temperature = Dict(:beta0 => Dict(("Na+", "Cl-") => (A₁, A₂, A₃, A₄, A₅))), fewer than five meaning the rest are zero. They are used by a PitzerActivityModel built with temperature_dependent = true, and vanish at exactly.

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._etheta_pairs Method
julia
_etheta_pairs(z, Aφ, sI, Ie) -> (Eθ, Eθ′)

The higher-order electrostatic terms of every pair of a set of like-sign ions of charges z, at the ionic strength Ie (regularized as the model regularizes it) with sI = √Ie: symmetric matrices, zero for two ions of the same charge.

ChemistryLab._pitzer_J Method
julia
_pitzer_J(x) -> (J, xJ′)

Pitzer's function and  , by Harvie's Chebyshev series in the variable z(x), evaluated by Clenshaw's recurrence with its derivative. The derivative is that of the approximation itself, so that is exactly the derivative of and the model keeps deriving from one energy. x > 0; a ForwardDiff.Dual passes through.

ChemistryLab._pitzer_tables_at Method
julia
_pitzer_tables_at(T, tables0, tablesT) -> tables

The coefficient tables at the temperature T: each entry A₀ of tables0 moved by its temperature terms in tablesT (see PitzerParameters). Each term vanishes at T_r exactly, so the tables are there their own values.

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, unless the model takes it as zero (missing_pairs = :zero).

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 — the cement set of Reardon (1990) 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(path; format = :auto) -> PitzerParameters

Read a Pitzer interaction-parameter set from a TOML file (format = :toml) or from the PITZER block of a PHREEQC-format database (format = :phreeqc); :auto takes the TOML reader for a .toml file and the PHREEQC one otherwise.

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"

Temperature terms (see PitzerParameters) are given as temperature = { beta0 = [A₁, A₂, A₃, A₄, A₅], Cphi = […] } in a [[binary]] entry and as temperature = [A₁, …] in a [[theta]], [[psi]] or [[lambda]] one, fewer than five meaning the rest are zero.

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.

The PHREEQC format

julia
PITZER
-B0
    Na+    Cl-    0.0765   A₁  A₂  A₃  A₄  A₅
-C0
    ...
-THETA
    Ca+2   Na+    0.07
-LAMDA
    CO2    Na+    0.1
-PSI
    Ca+2   Na+    Cl-     -0.007

-B0, -B1, -B2 and -C0 (Cφ) take a cation and an anion in either order, -THETA two ions of like sign, -LAMDA a neutral species and an ion, -PSI a triplet, each followed by A₀ and up to five temperature coefficients (Parkhurst and Appelo, 2013). Species are named as PHREEQC writes them, with @ appended to a neutral species, which is how this package names them. The identifiers this model has no counterpart for (-ZETA, -ETA, -MU, -ALPHAS, -APHI, which the water model provides) are skipped, and a warning names them; -use_etheta is the model's own keyword, etheta. No database is stored in this package: this reads the file the caller has (datapath("pitzer.dat") obtains PHREEQC's own), and every origin records its name and a truncated SHA-256 of its contents.

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._fill_site_potentials! Method
julia
_fill_site_potentials!(yv, cs, g)

Give a site primary that is not among the species the potential it actually has, instead of zero.

Why zero is right for charge and wrong for a site

The rule above it — a primary absent from the species contributes nothing — was written for the charge row, and there it is exact: every phase whose index is reported is neutral, so its coefficient in that row is zero and the potential multiplying it never appears.

A site row is not like that. A coupled family declares its component as the bare site, which is a component and not a substance, so it is absent from the species by construction. Every surface species carries it with coefficient one, and the host carries it with −ν. Zeroing it therefore shifts every surface index by y_site/ln 10 and the host's by −ν y_site/ln 10 — and y_site is not small: measured at +95.7 on an amphoteric oxide, so the surface indices would be off by more than forty log units.

Where the value comes from

Not from a convention. The free site is a species, it is always present, so its stationarity is exact:

and every y_c but the site's own is known. Solving that one equation for the missing one is the whole function. It is the same relation the solver enforced, read backwards.

A missing primary that no site family owns is left at zero, which keeps the charge row behaving as it did.

ChemistryLab._homotopy_rung Method
julia
_homotopy_rung(cs, A, i_w, n0, model, λ, start, ϵ, verbose, atol, rtol, T, P)
    -> 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.

The first rung is built at the temperature T and pressure P of the state the walk continues from, and every later one starts from the last. Built at the default 25 °C instead, the whole walk solved a problem at the wrong temperature, and the certificate, which reads the temperature off the state it is given, certified that problem: a pore solution meant for 20 °C came back 0.19 low in pH, the shift of pKw between the two temperatures, whenever this route was the one that succeeded.

ChemistryLab._repair_clipped_potentials! Method
julia
_repair_clipped_potentials!(yv, cs, A, g, n, ϵ)

Replace the potential of any primary whose species has collapsed to the activity floor, where reading it off that species is reading the floor.

The defect, measured

A primary's potential is read as that primary species' own μ/RT, which is exact while the species is there. Below the floor it is not there: ln a is clipped at ϵ and μ/RT reports the clip.

Measured on amorphous ferric hydroxide carrying a surface. Without the surface the solution sits at pH 7, Fe³⁺ is a 10⁻¹⁶ species, the clip barely bites and the solid reports log SI ≈ 0 as it must. Add the surface and its protolysis takes the pH to 8.1, where Fe³⁺ is nearer 10⁻²⁵: the clip is then worth nine decades, and the solid — present, and at an equilibrium the certificate accepts — reported log SI = 3.3. The certificate was right; it reads the solver's own multipliers. The index was reading a floor.

What replaces it, and why it cannot disturb anything else

Every aqueous species clear of the floor is a free variable, so at an equilibrium each one states μ_i/RT = Σ_c A_{c,i} y_c exactly. Over the set of them that is an overdetermined system for y, and the repair is its minimum-norm least-squares correction. Minimum norm matters: it puts nothing in the directions those species do not determine, so a primary supplied by _fill_site_potentials! keeps the value that function gave it.

This one assumes the state is an equilibrium

The rest of saturation_indices does not: it forms each species from the primaries' actual potentials, which is defined at any composition. This repair is different, because a primary that is not there has no potential to read and the only way to recover one is the stationarity of the species that are. It therefore fires only when a primary is at or below the floor — never otherwise, so an arbitrary state returns exactly what it always did — and on such a state the repaired index is a statement about the nearest equilibrium rather than about the state.

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.

For a model with a formula per ion, 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 same holds for DaviesActivityModel, TruesdellJonesActivityModel and LLNLActivityModel.

Not computed as a ratio of activity to concentration for those models. That ratio agrees for an abundant solute — and the tests check that it does — but it diverged for a species parked at the solver's lower bound while their closures added an ϵ inside the logarithm, returning 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.

A model whose coefficients couple the solutes, PitzerActivityModel and SITActivityModel, has no formula per ion to share. Its coefficient is read from its own closure, on the molality scale its concentration_scale declares,

with the closure's own floor ϵₐ, so that ln mᵢ cancels the term the closure added and γ is exact for a trace species too. Any model on the molality scale whose closure writes ln aᵢ = ln γᵢ + ln mᵢ that way is served by the same path; one on another scale is refused.

The solvent is reported as γ_w = a_w / x_w, which for the dilute model differs from 1 at second order in M_w Σ m, its water activity being exp(−M_w Σ m) rather than x_w.

Throws if the system has no aqueous phase.

Examples

julia
γ = activity_coefficients(eq, HKFActivityModel(å = 0.0, Ḃ = Ḃ_gems, Kₙ = 0.0))   # see HKFActivityModel
γ["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.activity_model_range Method
julia
activity_model_range(model) -> Union{Float64, Nothing}

The ionic strength, in mol/kg, up to which the manual states model is valid, or nothing when it states none.

The values are those of the table "Choosing a model" in the solving manual: the extended Debye–Hückel of HKFActivityModel up to about 1 mol/kg, the Davies equation of DaviesActivityModel up to about 0.5. Neither formula announces that it has left its range, and the consequence can be worse than an inaccurate coefficient: with å = 0, the Debye–Hückel limiting law plus a B-dot term, log γ keeps falling with I up to several mol/kg, so a pore solution can run away to an ionic strength no paste has, and the certified search not conclude. Comparing ionic_strength(eq) with this value is how a caller knows the answer is inside the model's range.

For HKFActivityModel the bound is that of the B-dot family as a whole. With a common ion size å = 0 the Debye–Hückel term is the limiting law, which the manual says is valid only far below it: read the bound as the most the model can be asked, not as what it delivers.

julia
ionic_strength(eq) <= something(activity_model_range(model), Inf)
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.saturation_indices Method
julia
saturation_indices(state, model, species; ϵ = 1e-16) -> OrderedDict{String, <:Real}

LogSI at state of each of species, solids that need not be species of the system: each is taken as a pure phase, of activity one, formed from the system's primaries by its primary_decomposition, with its standard Gibbs energy at the state's temperature and pressure,

julia
LogSI_s = [Σ_c a_cs y_c − ΔₐG⁰_s/RT] / ln 10 ,

y_c the potentials of the primaries read as the method without species reads them. For a pure phase of the system it is that method's index, to rounding. A species outside the span of the primaries is refused by name.

Examples

julia
si = saturation_indices(eq, model, [db["Gp"], db["ettringite"]])
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.MulticomponentRedlichKisterModel Type
julia
struct MulticomponentRedlichKisterModel{T<:Real} <: AbstractSolidSolutionModel

Redlich-Kister mixing of any number of end-members, as Redlich and Kister (1948) write it (their Eqs. 17 and 18): a series in the difference of the two mole fractions for each pair, and a term in the product of three for each triple,

every pair and triple evaluated at the mole fractions of the solution as they are, without renormalization. pairs maps an ordered pair (i, j) to its coefficients [L₀, L₁, L₂, …] in J/mol, the difference being xᵢ − xⱼ in that order; ternary maps an ordered triple (i, j, l) to [C, D₁, D₂] (or fewer) in J/mol. A pair or triple left out is zero. The coefficients B, C, D of Redlich and Kister (1948) are in units of 2.303 RT (decimal logarithms): multiply them by log(10) R T.

Two end-members with pairs = Dict((1, 2) => [a₀, a₁, a₂]) are RedlichKisterModel(; a0, a1, a2), and the first order is SubregularSolutionModel with a₀ = (W₁₂ + W₂₁)/2, a₁ = (W₂₁ − W₁₂)/2; the test suite checks both, and the activity coefficients against the worked ternary of (Redlich and Kister, 1948; Eqs. 22–23).

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

Examples

julia
julia> m = MulticomponentRedlichKisterModel(Dict((1, 2) => [4000.0, 500.0], (2, 3) => [2500.0]);
                                            ternary = Dict((1, 2, 3) => [1000.0]));

julia> m.pairs[1], m.ternary[1]
((1, 2, [4000.0, 500.0]), (1, 2, 3, [1000.0, 0.0, 0.0]))
ChemistryLab.MulticomponentRedlichKisterModel Method
julia
MulticomponentRedlichKisterModel(pairs; ternary = Dict(), n = nothing)
    -> MulticomponentRedlichKisterModel

Construct a MulticomponentRedlichKisterModel from the series of each ordered pair and the terms of each ordered triple, in J/mol (bare numbers) or as energies per mole (Quantity). n, the number of end-members, defaults to the largest index named.

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

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.

CEMDATA18 ships the AFm and AFt binaries under two names each, so that the database itself declares each binary twice; instances = 2 declares the second composition with a keyword.

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)

instances = :auto declares one composition and allows a second. The phase is solved with one, and equilibrate_certified adds the second only when its certificate finds the phase wanting to split, through with_instances: the answer is then a state of the enlarged system, whose solid_solutions hold the two instances. A caller whose composition stays outside the gap never pays for the second copy. :auto is refused for a convex model, as instances = 2 is.

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.SubregularSolutionModel Type
julia
struct SubregularSolutionModel{T<:Real} <: AbstractSolidSolutionModel

Asymmetric (subregular) Margules solid solution of any number of end-members, in the form Helffrich and Wood (1989) give it (their Eq. 5′):

W[i, j] is in J/mol; the diagonal is ignored. On the edge – the excess is  , so is at infinite dilution of in . ternary maps a triple of end-member indices to in J/mol; a triple left out is zero. The ternary coefficients are not determined by the binaries, as Helffrich and Wood (1989) show, and no coefficient of higher order exists in this model.

On the simplex the binary terms are those of RedlichKisterModel to first order,    with    and   , since     there. So two end-members are RedlichKisterModel(a0 = (W₁₂ + W₂₁)/2, a1 = (W₂₁ − W₁₂)/2), and a symmetric W without ternary terms is RegularSolutionModel; the test suite checks both, and the activity coefficients against (Helffrich and Wood, 1989; Eq. 6′).

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

Examples

julia
julia> m = SubregularSolutionModel([0.0 3000.0 0.0; 5000.0 0.0 0.0; 0.0 0.0 0.0];
                                   ternary = Dict((1, 2, 3) => 1500.0));

julia> m.W[2, 1], only(m.ternary)
(5000.0, (1, 2, 3, 1500.0))
ChemistryLab.SubregularSolutionModel Method
julia
SubregularSolutionModel(W::AbstractMatrix; ternary = Dict()) -> SubregularSolutionModel

Construct a SubregularSolutionModel from the matrix of the binary coefficients W[i, j] = and the ternary coefficients ternary, keyed by triples of end-member indices in any order, all in J/mol (a bare number) or as energies per mole (a Quantity).

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

The asymmetric formalism of Holland and Powell (2003), a macroscopic van Laar model of any number of end-members: each end-member carries a size parameter  , the mole fractions are weighted into  , and

W symmetric, in J/mol, with a zero diagonal (ignored), and α dimensionless: only their ratios matter. Equal sizes give the regular model (RegularSolutionModel), which Holland and Powell (2003) call the symmetric formalism, and two end-members the van Laar binary, asymmetric as soon as the sizes differ. W and α are those of the temperature and pressure of the calculation; the papers that fit them give their dependence, which the caller evaluates.

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

Examples

julia
julia> m = VanLaarModel([0.0 16000.0; 16000.0 0.0], [1.0, 0.643]);

julia> m.α
2-element Vector{Float64}:
 1.0
 0.643
ChemistryLab.VanLaarModel Method
julia
VanLaarModel(W::AbstractMatrix, α::AbstractVector) -> VanLaarModel

Construct a VanLaarModel from the symmetric interaction energies W (J/mol, or energies per mole) and the size parameters α of the end-members.

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

An interaction parameter as a bare number of joules per mole, checking the unit when one is given.

The error it exists to catch

RedlichKisterModel(a0 = 4000.0) and RedlichKisterModel(a0 = 0.188) are both accepted, and one of them is a published parameter set read in the wrong convention. The literature writes these coefficients two ways: CEMDATA18 gives the AFm sulfate/hydroxide binary as A₀ = 0.188, A₁ = 2.49 in RT units, while this package – and the Guggenheim form its ln γ is written in – takes J/mol and divides by RT internally. The factor between them is RT ≈ 2478 J/mol at 25 °C, and nothing in a bare Float64 says which one was meant.

A Quantity says. u"J/mol" passes and is stripped; u"kJ/mol" is converted; anything else raises rather than being silently reinterpreted. A bare number is still accepted and still means J/mol, so no existing call changes – the unit is an opportunity to be explicit, not a new obligation.

RT units have no unit to give, which is the point: a caller holding A₀ = 0.188 has to multiply by RT themselves, and that multiplication is where the convention becomes visible.

Zero is exempt, and has to be

[0.0 4.0u"kJ/mol"; 4.0u"kJ/mol" 0.0] is how one writes a binary W matrix, and Julia promotes that literal to a common element type before this function ever sees it: the bare 0.0 on the diagonal arrives as a dimensionless quantity. Refusing it would mean writing 0.0u"J/mol" on a diagonal the model documents as ignored. A zero also carries no convention to get wrong – 0 J/mol, 0 kJ/mol and 0 RT are the same number – so it passes whatever dimension it arrived with, and every nonzero value is still checked.

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, and SubregularSolutionModel the two-parameter one, a₀ = (W₁₂ + W₂₁)/2 with a₁ = (W₂₁ − W₁₂)/2, which is why they 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._warn_degenerate_end_members Method
julia
_warn_degenerate_end_members(name, end_members, model)

Warn when a member of a solid solution is, to within 0.1 RT per formula unit, a mixture of two others: its composition is λ times one member's plus 1 − λ times another's, 0 < λ < 1, and so is its standard Gibbs energy at 25 °C. The phase then holds one substance twice, once as a member and once as a mixture of two, and ideal mixing counts its configurations twice. The published MgAl-OH-LDH ternary of Cemdata18 is such a case (its M6A member is the average of M4A and M8A, 0.3 J/mol apart), and acknowledge_degenerate = true keeps it as published. A member that is such a mixture with a Gibbs energy well away from it, as the ordered T5C of CSH3T (−4.35 kJ/mol) or the siliceous hydrogarnet C3AS0.41H5.18 (−10.4 kJ/mol), is an ordering energy the model means, and is not reported.

For a sublattice model the question is asked of the occupancy instead: a direction along which no site fraction changes and the Gibbs energy does not either is refused, since the minimization then has no unique answer.

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. It does not compute x̄, the overall composition of the phase in a full aqueous equilibrium: the equilibrium does. When the element balance pins x̄ inside the gap, a phase declared twice (instances = 2) separates onto the pair by itself; when nothing pins it, two instances started at one composition stay there, the symmetric state being stationary, and equilibrate_split or instances = :auto seeds the split.

So read x̄ as the overall composition you have — from a 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.mixing_convexity Method
julia
mixing_convexity(model, n_members; T = 298.15) -> (; verdict, witness, how)

Whether the molar Gibbs energy of mixing of a solid solution is convex over its whole composition simplex: verdict is :convex, :nonconvex or :undecided, witness a composition (mole fractions) where it is concave when :nonconvex, nothing otherwise, and how the argument, in words.

Convexity is what makes a certified answer a global minimum and what lets one composition describe the phase; spinodal_interval answers the question for a binary, and this answers it for any number of end-members.

The arguments, by model

  • Ideal mixing and ideal mixing on sublattices (SublatticeModel) are convex by construction: is convex, and a sum of such terms in site fractions linear in is too.

  • A binary with a Redlich-Kister form is decided by the scan of spinodal_interval.

  • A regular model (RegularSolutionModel), with  : on the tangent space of the simplex,    for every unit with   (the ℓ1 norm of such a is at least , and  ), so the Hessian of is at least   there, an orthonormal basis of the tangent space. Hence convex when   , which for two end-members is the classical  . Not convex when some  : at the midpoint of that edge the second derivative along it is   . Between the two, the smallest eigenvalue of the projected Hessian is searched on a lattice of compositions: a negative one is a witness, and none found leaves the question :undecided.

  • A compound-energy model (CompoundEnergyModel) is judged in its site fractions, where its Gibbs energy lives, and needs g, the of its end-members at T (:undecided without them). On each site the configurational term gives a tangent curvature of at least (the bound above) and the site interactions add the smallest tangent eigenvalue of their matrix ; the reference surface adds only a coupling between two sites, constant when there are two: the matrix of the of the compounds projected on the two tangent spaces. With two sites the energy is therefore convex when    , both factors positive. Otherwise the Hessian is sampled on a lattice of site fractions: a negative eigenvalue is a witness, none found leaves :undecided. The CASH+ core model is convex by that bound; CASH+NK is not decided by it.

  • Any other model is sampled the same way, through ForwardDiff of its : a witness or :undecided, never a proof of convexity.

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 taken by forward-mode differentiation on a grid of compositions rather than in closed form, so that the same routine covers a₂ and the symmetric RegularSolutionModel without a separate derivation. The interval returned is the span of the grid points where it is negative. 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 mixing_convexity.

See also: RedlichKisterModel, RegularSolutionModel.

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

Mixing on sublattices with the two terms of the compound energy formalism that ideal site mixing (SublatticeModel) leaves out: the reference surface, on which the Gibbs energies of the end-members are interpolated in the site fractions, and regular interactions between the species of one site. It is the model of the CASH+ C-S-H of Kulik et al. (2022).

julia
CompoundEnergyModel(lattice::SublatticeModel; interactions = ())

interactions lists (site, species_1, species_2, W), the site and its two species named as lattice.sites and lattice.species name them, and W in J/mol of formula units.

The Gibbs energy

Per formula unit, with the fraction of site held by species and the standard Gibbs energy of end-member , which puts species on site :

It depends on the site fractions alone. The first term is the reference surface: it equals at the composition of end-member and is linear in the fractions of each site, so it differs from the mechanical mixture only where the end-members are not independent in the site fractions. That difference is the energy of the reciprocal reactions, such as TSvh + TCCh = TSCh + TCvh in CASH+, whose Gibbs energy is not zero.

Activities

The chemical potential of end-member is . With the reference surface and the number of sites,

The second term vanishes at each end-member and is unchanged when every is shifted by the energies of its elements, so the convention of the database does not matter. It needs those energies: the solver hands them to the model with the amounts, at the temperature of the solve.

The end-members must be every compound

The reference surface is spanned by the compounds, one per choice of a species on every site. The model requires the end-members to be exactly those compounds, each once; the constructor refuses any other set. The Gibbs energy and the element balance then depend on the site fractions alone, and the end-member amounts are not unique: six end-members for four independent site fractions in CASH+, twelve for six in CASH+NK. Every split of the same site fractions between the end-members is the same equilibrium.

The split the solver is given

To make the split unique, the energy the solver minimizes adds  , with

the divergence of from the product of its site fractions: it vanishes when  , and only then. At fixed site fractions the minimum is that product, where the added energy is zero and so is its gradient,  . The equilibrium, its site fractions and the chemical potentials of the members are those of the model; the member amounts are the product of the site fractions. With the null directions gone, the activity of each member goes as its mole fraction when it vanishes, and the phase is inverted as a mixture of end-members is. Without it, the answer is one of the equivalent splits, which one depending on where the search started; measured on twelve pastes of CASH+NK (Ca/Si 1 and 1.6, sodium, potassium or both), the two forms certify the same equilibria at the same cost, and only this one returns the same member amounts whatever the start.

Example

julia
julia> lat = SublatticeModel([1.0, 1.0], ["A" "A" "B" "B"; "X" "Y" "X" "Y"]; sites = ["s1", "s2"]);

julia> m = CompoundEnergyModel(lat; interactions = [("s1", "A", "B", -5000.0)]);

julia> length(m.interactions)
1

See also: compound_energy_model, site_fractions.

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

Ideal mixing on several distinct sites (sublattices) of one formula unit: the model of a solid solution whose end-members differ by what occupies a few structural positions, such as the bridging tetrahedra and interlayer sites of a C-S-H.

The model

A formula unit has sites s = 1 … S, site s counted mₛ times in the formula. Every end-member k puts exactly one species σₛ(k) on each site. The fraction of site s held by species i is the sum of the mole fractions of the end-members that put i there,

and the species mix at random on each site, independently of the other sites. The configurational Gibbs energy per mole of formula units is therefore

and the activity of an end-member is the product of its own site fractions:

That is exactly  : with    the amounts on each site, the derivative of is   , and the bracket loses its last two terms because every end-member fills every site once. The log activities are the gradient of one energy, so a certificate keeps its :global_minimum scope.

Two properties the solver has to know

Convexity. is a sum of convex functions of the site fractions, which are linear in , so it is convex. It is strictly convex exactly when the occupancy matrix — one row per (site, species), one column per end-member,    — has full column rank: a null vector of is a direction along which no site fraction changes. rank stores that rank.

The dilute limit. When  , the site fractions of the species that alone carries vanish with it, and the others do not, so

An end-member with   has no species of its own: its activity stays finite as it disappears, so it can be exactly absent from a phase that is present — its stationarity is then an inequality, as a pure phase's is. exponents stores the .

Fields

  • sites: the name of each site;

  • multiplicity: , the coefficient of each site in the formula unit of the end-member records the model is used with;

  • species: for each site, the names of the species found there;

  • occupancy: occupancy[s, k] is the index in species[s] of the species end-member k puts on site s;

  • exponents: ;

  • rank: the rank of .

A site that holds one species in every end-member contributes nothing and may be left out.

Construction

julia
SublatticeModel(multiplicity, occupancy; sites)

with occupancy a matrix of species names, one row per site and one column per end-member, in the order of the end-members of the phase. Two end-members with the same occupancy on every site are refused: the model cannot tell them apart. sublattice_model builds the model of a published phase from data/literature.

Example

The CSH3T model of Kulik (2011), Eq. 19, on its own formula unit: two sites holding Si or Ca, and the pentameric end-member T5C carrying Ca on the first and Si on the second.

julia
julia> m = SublatticeModel([1, 1], ["Si" "Ca" "Ca"; "Si" "Si" "Ca"]; sites = ["BTI1", "BTI2"]);

julia> m.exponents   # TobH and T2C each own a species on one site; T5C owns none
3-element Vector{Int64}:
 1
 0
 1

julia> m.rank
3

See also: site_fractions, sublattice_model, IdealSolidSolutionModel.

References

ChemistryLab.SublatticeModel Method
julia
SublatticeModel(multiplicity, occupancy; sites) -> SublatticeModel

Build the model from the coefficient of each site and a matrix of species names, occupancy[s, k] being what end-member k puts on site s. The species of each site are numbered in the order they first appear along the row.

ChemistryLab.compound_energy_model Method
julia
compound_energy_model(ref, end_members) -> CompoundEnergyModel

The compound-energy model a published source gives for a phase, read from data/literature/<key>.json: the sites and occupancy as sublattice_model reads them, and the site interactions from the table <model>_interactions (columns site, species_1, species_2, W).

Example

julia
subs = Dict(symbol(s) => s for s in build_species(datapath("cemdata18-cashplus.json")))
m = compound_energy_model("Kulik2022:cashplus", [subs[n] for n in ("TSvh", "TSCh", "Tvvh", "TCvh", "TvCh", "TCCh")])
ChemistryLab.site_fractions Method
julia
site_fractions(model::SublatticeModel, x) -> Vector{Vector}

The site fractions at the mole fractions x of the end-members: site_fractions(model, x)[s][i] is the fraction of site s held by model.species[s][i]. Each inner vector sums to sum(x). The element type is that of x, so the function carries dual and symbolic numbers.

ChemistryLab.sublattice_model Method
julia
sublattice_model(ref, end_members) -> SublatticeModel

The sublattice model a published source gives for a phase, read from data/literature/<key>.json, with its end-members in the order of end_members.

ref is "<key>:<model>", such as "Myers2014:cnash" or "Kulik2011:csh3t"; the file holds the tables <model>_sites (the coefficient of each site in the source's formula unit) and <model>_occupancy (the species each end-member puts on each site, one row per end-member and site, with the database name of the end-member in its first column).

Formula units

A database may store an end-member on another formula unit than the source's: the Cemdata18 records of CSH3T are half the units of Kulik (2011). The site coefficients belong to the source's unit, and a model applied to records f times as large has coefficients f times as large, because the configurational entropy is extensive. When end_members are species, their formulas are compared with the table <model>_end_members, the formula the source prints for each end-member: every record must be the same multiple f of its published formula, and the coefficients are multiplied by f. A record that is not such a multiple is refused, since the model then describes another substance. When end_members are names, f = 1 and no check is made.

Example

julia
subs = Dict(symbol(s) => s for s in build_species(datapath("cemdata18-thermofun.json")))
m = sublattice_model("Kulik2011:csh3t", [subs[n] for n in ("CSH3T-TobH", "CSH3T-T5C", "CSH3T-T2C")])
m.multiplicity   # [0.5, 0.5]: the records are half the published unit

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,
                      inner_tol, inner_maxit, inner_fall_bound,
                      lenient_line_search, 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.

The keywords are OptimaSolver's DualNewtonOptions: maxit bounds the Newton search, whose active set is reclassified at every iteration, and inner_tol and inner_maxit the inner fixed point that recovers the phase compositions. max_active_updates is no longer used, and deprecated. tol and si_tol are also the thresholds of the certificate this solver issues, so loosening them loosens the proof, not only the search. inner_fall_bound (30) and lenient_line_search (false) change the path of the search and not the conditions of its answer; the steps of a kinetic run set them to Inf and true, which is where they pay (see OptimaSolver's DualNewtonOptions).

See also: optimality_certificate, speciated_states.

ChemistryLab._activity_range_report Method
julia
_activity_range_report(model, state) -> (I, I_max, within)

The effective ionic strength of state (NaN without an aqueous phase), the range activity_model_range states for model, and whether the first lies within the second (nothing when either is unknown).

ChemistryLab._aqueous_convexity Method
julia
_aqueous_convexity(model, system, p) -> (proved, reason)

Whether the Gibbs energy of the aqueous phase under model is proved convex over its whole domain, every composition of the declared solutes at the temperature and pressure of p, and if not, why. It is a sufficient condition, and it is what a :global_minimum scope rests on; Gibbs–Duhem consistency, which the symmetry of the Jacobian measures, is checked before and is not repeated here.

  • The ideal dilute model: its energy, Σₛ nₛ (ln(nₛ/(n_w M_w)) − 1) plus a linear term, has the Hessian Σₛ nₛ (vₛ/nₛ − v_w/n_w)², positive.

  • A Debye–Hückel form with one function of I times zᵢ² on every ion, with the solvent row that integrates it: the HKF model with one ion size, Ḃ = 0 and no salting-out, and Davies with neutral species carrying none. Its excess energy adds a term of rank one to the ideal Hessian, which the Cauchy–Schwarz inequality bounds by it as long as ln(10)·A·z_max²·√I / (2(1 + κ√I)²) ≤ 1, κ = B å (one for Davies, whose b I term only adds convexity), and the supremum over I, at √I = 1/κ, gives ln(10)·A·z_max²/(8κ) ≤ 1. See What the certificate proves, and when for the derivation.

  • Any other model: not proved, whether or not it is convex.

ChemistryLab._certificate_scope Method
julia
_certificate_scope(des, p, n, constraint; Aq = nothing, floor = _CERTIFICATE_FLOOR)
    -> (scope, reasons, curvature)

What optimal = true proves for the problem at the composition n:

  • :global_minimum when the log activities are the gradient of one Gibbs energy, that energy is proved convex over the whole domain — the ideal phases, the aqueous models of _aqueous_convexity, convex solid solutions and ideal site mixing — and the constraint is variational, so that the optimality conditions are sufficient;

  • :local_minimum when the activities are such a gradient and the constraint variational, convexity is not proved, and the Hessian of the energy over the directions that conserve matter is positive definite at the answer: a strict local minimum, which a non-convex model can hold beside others;

  • :kkt_point when that Hessian is not positive definite — a saddle if a direction lowers the energy, a minimum that is not strict if one leaves it flat — or when the constraint is not variational;

  • :self_consistent when the activities are not the gradient of one energy, which is the case of the extended activity models in general use and of a diffuse layer: the conditions then state a composition consistent with its own activities, not the minimum of an energy.

reasons says which property decided, in words, and curvature is the smallest eigenvalue of _reduced_curvature (NaN where it is not computed: no energy, or a constraint that is not variational).

ChemistryLab._complete_floored_solutes Method
julia
_complete_floored_solutes(des, prob, res, floor) -> Vector{Float64}

The amounts of res, with each solute below the activity floor given the amount the potentials of the answer give it.

The activity models read a solute through max(n, floor) (see _ACTIVITY_FLOOR), so below the floor its log activity does not depend on its amount. The dual solve recovers a solute from its own stationarity, w ← w + r with r = uᵢ − ∇fᵢ, a step that assumes the log activity follows ln n: below the floor it gains r per sweep, whatever the distance to its answer, and it can be left there. The equation of the model is then solved exactly by n = floor·exp(r) when r > 0; when r ≤ 0 any amount below the floor satisfies it, and the one computed is kept. The element balance moves by the amounts added, of the order of the floor.

Measured with the floor at 1e-16, on a cement paste at pH 14 in 7.8 g of water, whose equilibrium holds 1.2e-16 mol of H⁺: a solve left it at 3e-100 mol, and pH(eq, model), which reads the amount, came out 0.09 high. With the floor at 1e-30 such a solute is far above it, and this is the net for one that is not.

ChemistryLab._dual_problem Function
julia
_dual_problem(des, p, n0, blocks = nothing) -> 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. blocks are the constraint blocks of the solve, those of a fixed temperature and pressure when nothing.

ChemistryLab._is_variational Method
julia
_is_variational(constraint) -> Bool

Whether a constraint keeps the equilibrium a minimization whose optimality conditions are sufficient. CapillaryWater shifts the activity of water by an amount that is not the gradient of a convex energy for an arbitrary retention law, and the constraints that make the temperature or the pressure an unknown turn the solve into a KKT system whose convexity is not established here.

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

The balance as the certificate judges it: balance_relative, the worst row relative to what it holds, as OptimaSolver's certificate judges it. Until 0.41.0 the larger of that and balance in moles, which held a large row to an amount of matter and made the verdict depend on the size of the system. A certificate without the relative figure gives balance.

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 (with its one-sided form on the members of a present phase below the floor). 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.

The balance enters as the certificate judges it (_judged_balance): relative to what each row holds below one mole, in moles above it, so that a candidate that lost a trace ranks below one that kept it. The others are log-activities. It ranks candidates of one problem; it is not a measure to compare two problems with.

ChemistryLab._lift_equilibrium Method
julia
_lift_equilibrium(des, state, eq, b; ϵ, constraint, q, surface_potential, strip_tag,
                  floor = _CERTIFICATE_FLOOR) -> (ChemicalState, q)

The equilibrium eq, found on the values of the data, carrying the derivatives the dual numbers of state (its temperature and pressure) and of the budget b imply. The problem is the one solve poses, constraint and surface unknowns included, built at the caller's duals; its answer is lifted by the implicit-function theorem at eq with the active set frozen (OptimaSolver's dual_newton_tangent), so the derivatives are those of the solution, not of the iteration that found it, and each level of a nested differentiation is taken in turn. q is what the solve found for the constraint's unknowns, returned with its derivatives. A pure phase holding floor or less is absent: the floor of the certificate for a certified answer, the solver's lower bound for an interior-point one, which leaves an absent phase there.

ChemistryLab._reduced_curvature Method
julia
_reduced_curvature(J, A, n, Aq, floor) -> Float64

The smallest eigenvalue of the Hessian of the Gibbs energy over the directions that conserve matter and move only the species present, J being the Jacobian of the log activities at n (the Hessian of G/RT, symmetric where this is called).

The species present are those above floor; the directions are the null space of their columns of the conservation matrix, with the columns Aq of a constraint's titrant, along which the energy is linear. The coordinates are scaled by √nᵢ, so that an ideal solute contributes exactly one whatever its amount, and the figure is dimensionless. Positive, the point is a strict local minimum on its active set; negative, a direction lowers the energy and the point is a saddle. Inf when no direction is free.

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, stationarity_floored, balance, balance_relative,
         worst_supersaturation, n_interior,
         n_absent_component, param_residual, worst_violation_split,
         split_phases, split_trials, optimal, scope, scope_reasons,
         reduced_curvature, stationarity_abs,
         ionic_strength, activity_range, within_activity_range)

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. Whether the problem is one is not always the case, and scope says what the proof covers here: :global_minimum, :local_minimum, :kkt_point or :self_consistent, with scope_reasons naming the property that decided and reduced_curvature the smallest eigenvalue of the Hessian over the directions that conserve matter (see _certificate_scope). The extended activity models in general use (B-dot, Davies with neutral species) are not the gradient of one Gibbs energy, and a certified equilibrium computed with them is a composition consistent with its own activities rather than the minimum of an energy. stationarity_abs is the stationarity residual in RT units, before the scaling stationarity applies. Use it to audit any solver — including EquilibriumSolver, whose interior-point back ends return their last iterate without stating whether it is the minimum.

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). The balance is reported twice: balance, the worst row in moles, and balance_relative, the worst row relative to what it holds. Each row is judged on the larger of the two below one mole and in moles above it, so that a trace is held to its own amount, as PHREEQC and GEMS hold a mass balance to its element total; judged in moles alone, a trace of 1e-9 mol could be 10 % wrong and certified. A component whose budget is below 1e-12 of the largest is one nobody supplies, its carriers held at the floor, and its row is judged in moles. stationarity_floored is the first, one-sided, on the members of a present phase held below floor: such a member may hold more than its exact amount, by truncation, but not less, and it fails when the search left it far below what the multipliers give it.

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.

ionic_strength is the effective ionic strength of the composition (NaN without an aqueous phase), activity_range the ionic strength up to which the activity model is stated valid (activity_model_range, nothing when no range is stated), and within_activity_range whether the first lies within the second (nothing when either is unknown). An answer can certify outside the range of its model: the certificate says so, and the caller decides.

ChemistryLab.solve_certified Method
julia
solve_certified(des, starts; b = nothing, ϵ = 1e-16, floor = 1e-25, constraint = FixedTP(),
                parameters = nothing, memo = nothing, stop = nothing, report = true,
                y0 = nothing, dual_state = nothing)
    -> (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.

Each candidate is judged on the optimality conditions alone; what the proof covers (scope, reduced_curvature, the range of the activity model) is computed once, for the answer returned. report = false returns the verdict without it, for a caller that only reads optimal and the residuals. stop, a function of a candidate's certificate, ends the search on the first candidate it accepts; parameters, a Ref, receives the unknowns of the constraint at the answer. y0 and dual_state are those of the solve (SciMLBase.solve(::DualEquilibriumSolver, …)): the multipliers each start is resumed from, for a start that is the answer of an earlier solve, and a Ref that receives the multipliers of the answer returned.

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.

Offering the same start twice

memo, when given, is an IdDict that records the answer and the certificate reached from each start, keyed on the start object itself, and a start met again is answered from it instead of being solved a second time. The dual solve and the certificate are deterministic functions of the start once des, b, ϵ, floor and constraint are fixed, so the record is the result, bit for bit; the caller owns the dictionary and must not share it between calls that differ in any of those. equilibrate_certified keeps one for the duration of a call, because its search offers the same cached starts again after the continuation and after each restart: measured on a CEM I paste that does not certify at once, twelve of the thirty-six dual solves of one call were repeats of earlier ones.

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

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

An answer of this solver holds an absent pure phase at exactly zero, and is a start as it stands. y0, the multipliers of that answer, resumes the search where it ended, as the next solve of a sequence; a Ref passed as dual_state receives those of this one, (; y, q, route, iterations), to hand to the next.

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._aqueous_form Method
julia
_aqueous_form(model, cs, members) -> Union{Nothing, NamedTuple}

What the inversion of the aqueous phase needs from an activity model: the charge and the log₁₀ γ(I) of each member, the molar mass of the solvent and the Debye–Hückel parameters at given conditions. nothing for a model whose coefficients depend on more than the ionic strength (Pitzer, SIT), which the sweeps of OptimaSolver then handle.

ChemistryLab._aqueous_inverter Function
julia
_aqueous_inverter(des, pq = nothing) -> Union{Nothing, Function}

The invert of OptimaSolver's SolutionPhase for the aqueous phase of des: (c, ref, w, q, params) -> log-amounts, or nothing from the call when no ionic strength solves the system. pq(q, params) gives the parameters the activity model sees when a constraint makes them unknowns (the temperature of an adiabatic solve); without it they are params. For SIT and Pitzer, whose coefficients depend on more than the ionic strength, it is Newton's method on the solutes (_newton_inverter). Returns nothing itself for a model neither covers.

ChemistryLab._branch_of_iterate Method
julia
_branch_of_iterate(F, s1, s_w, s_ceil; h0 = 0.05) -> Float64

The root of F(s) = ln S(eˢ) − s nearest s_w, the log ionic strength of the composition the solve holds, when s_w lies above the first root s1 and below s_ceil; s1 itself otherwise, and whenever no other root lies nearer.

A charged species of high valence makes F turn back up above its first root: its activity coefficient falls so fast with I that its amount, and I with it, grow faster than I, and F crosses zero three times. Under the balances the equilibrium can sit on the middle root, where the open solution would be unstable: a solution of amorphous silica in KOH at 80 °C, Cemdata18's extended Debye–Hückel and the tetramer Si4O10-4 (z = −4) carrying most of the dissolved silicon, has its root at 1.38 mol/kg between two others at 0.45 and 4.3. The first root alone never gives that composition, and the solve then never closes its balances; following the branch of the iterate does. Where F has one root, the search finds s1 and returns it unchanged.

The search widens a window about s_w, doubling it from h0, inside [s1, s_ceil], until F changes sign on one side, and returns the root there, bracketed. s_ceil keeps it below the stated range of the model times a factor (_BRANCH_CEILING_FACTOR).

ChemistryLab._ionic_strength_root Method
julia
_ionic_strength_root(F, dF, s_start; smax = log(1e4), hmax = 4) -> Union{Nothing, Float64}

The first root of F(s) = ln S(eˢ) − s above s_start, where F crosses from positive to negative: the self-consistent ionic strength eˢ of a solution whose solutes follow from it, on the branch connected to the dilute limit. nothing when that branch holds none.

The first one, and not the one nearest the composition the solve holds. The limiting law has up to three: past its range it can dip below zero over a narrow interval only, two roots close together, and the coefficient Ḃ I brings a third back far above, at thousands of mol/kg for a CEM I pore solution with sodium chloride. The physical answer is the dilute branch; taking the root nearest a start made the answer depend on the path, and a search stepping up from a start sitting on the root itself, with F at +2e-16, stepped over a dip 0.15 wide in ln I. The steps are Newton's from the left, which a convex F keeps short of its first root, and a dip a step spans shows as F′ changing sign there and is found whatever its width.

A dip that stays above zero ends the dilute branch without a root, and so does smax, 1e4 mol/kg, a molality no aqueous model is stated for: the potentials then hold no solution on it, which is the answer. Searching on past such a dip reached the third root, and from the potentials of a cold start the solve then began at a composition three orders of magnitude off, which it never left: the start of a chloride-loaded paste failed where recovering the solutes one by one from that start had certified it at once. s_start lies below the root (_invert_aqueous); a start already past it, with potentials that hold more than smax at zero ionic strength, is searched downwards, under the same ceiling.

ChemistryLab._newton_inverter Function
julia
_newton_inverter(des, pq = nothing) -> Union{Nothing, Function}

The invert of the aqueous phase for SIT and Pitzer: the solutes recovered from their potentials by Newton's method on their log-amounts, with the model's exact Jacobian. nothing from the call when the iteration finds no composition, which OptimaSolver reads as it reads the inversion through the ionic strength: an iterate is then swept, and a trial of its line search passed over. Returns nothing itself for a model it does not cover.

ChemistryLab._newton_predictor_form Method
julia
_newton_predictor_form(model, cs, members) -> Union{Nothing, NamedTuple}

The Debye–Hückel part of a model the Newton inversion covers, as an _aqueous_form, which gives that inversion its start; nothing for a model it does not cover.

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.UncertifiedEquilibrium Type
julia
UncertifiedEquilibrium <: Exception

Raised when no route produced a certified equilibrium and the caller asked for one: equilibrate(state; certify = true), or any solve under STRICT_CONVERGENCE. The answer with the smallest KKT error and its certificate travel with it, to be audited:

julia
try
    equilibrate(state; certify = true)
catch e
    e isa UncertifiedEquilibrium || rethrow()
    e.certificate.balance_relative, e.state
end

The message names the conservation row the answer misses most, relative to what it holds: an element present as a trace is where a wrong answer shows first.

Fields

  • msg: what failed, with the stationarity, the balance and that row;

  • state: the answer with the smallest KKT error, a ChemicalState;

  • certificate: its certificate (optimality_certificate).

ChemistryLab._DeferredStarts Type
julia
_DeferredStarts(make)

Starts computed only when the search reaches them, and once: make() returns an iterable of starting states, or nothing, and is called the first time the iterator is iterated; the search offers its starts again after a continuation or a restart, and the second pass reads what the first made. For a start that is a solve in its own right, such as _ideal_mixing_start, which a search that certifies from an earlier start never needs.

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 — unless offer_tail is false.

ChemistryLab._activity_range_hint Method
julia
_activity_range_hint(eq, model) -> String

A sentence for the refusal of equilibrate_certified when the answer it gives up on has an ionic strength past the range the manual states for model (activity_model_range), and an empty string otherwise.

Measured on blended cement pastes taken to full reaction, the refusal can be the activity model's rather than the solver's. The Debye–Hückel limiting law, HKFActivityModel(å = 0), has its log γ still falling at several mol/kg, and the search does not conclude; with an ion size per ion the same budgets certify. Without this sentence the message names a stationarity and a balance, and nothing in it points at the activity model.

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

nothing for real-valued data; the certified answer with its derivative attached when the amounts, the temperature, the pressure or the budget b carry ForwardDiff.Dual numbers. The derivative is that of the answer the search certified, lifted one level of duals at a time (_lift_equilibrium), so nested differentiations are each exact and kept apart by their tags.

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._auto_split Method
julia
_auto_split(eq, cert; model, b, kwargs...) -> (state, certificate)

What instances = :auto asks for. When the certificate finds a phase declared :auto wanting to split, the system is rebuilt with a second instance of it (with_instances), the answer is carried over with the copies at zero, and the passes of equilibrate_split run on the enlarged system, seeded from the composition the stability test gives. Their answer is kept only when its KKT error is smaller, so the result is never worse than the first answer; it then carries instances_added, the phases that were given a second instance.

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_mixing_start Method
julia
_ideal_mixing_start(state, model, bfix, ϵ, constraint, verbose; kwargs...)
    -> Union{ChemicalState, Nothing}

A certified answer to the same problem with every sublattice phase mixing its end-members ideally, as a state of state.system, to start the search from. nothing when the system has no sublattice phase, or when that solve does not certify.

The two problems differ by the site-mixing term alone, and the one without it is the solve the package has always done. Measured on a CEM I paste with the CNASH gel of Myers et al. (2014) (CEMDATA18, its eq. C.1 activity model): certified in 20.5 s from the state as given, in 0.67 s from this start, which costs 0.25 s itself.

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._linear_program Method
julia
_linear_program(des, state, bfix) -> Union{NamedTuple, Nothing}

The linear program over the pure phases of the problem equilibrate_certified is about to solve: minimize gᵀn subject to A n = b, n ≥ 0, with the standard potentials, the conservation matrix and the budget of the dual solve itself (OptimaSolver's lp_start). nothing when it could not be set up.

The program is what the equilibrium becomes when every activity is one, and it answers two questions at the cost of a few hundred pivots.

  • Can the budget be met at all? The minimization has no solution when no non-negative amounts of the declared species reproduce b, and the program says so with a Farkas vector z, Aᵀz ≥ 0 and bᵀz < 0, verified: a combination of the balances that every species raises and the budget lowers. The cascade would otherwise spend its restarts and continuation on a problem without an answer.

  • Where to start. Its vertex is the assemblage of pure phases of lowest standard energy, and its multipliers give every other species an amount; see _lp_lifted_state.

ChemistryLab._lp_lifted_state Method
julia
_lp_lifted_state(des, state, lp, ϵ) -> ChemicalState

The start the linear program gives: its vertex for the species it holds, and for the others the amount the multipliers y of the program give them in their own phase. With u = −Aᵀy, a species absent from the vertex has the activity aⱼ = exp(uⱼ − gⱼ), which is at most one, and its amount follows from what the activity of its phase is:

  • a solute: a molality, so aⱼ mol/kg of the water at the vertex;

  • a member of a solid solution present at the vertex: a mole fraction, so aⱼ times the amount of the phase there;

  • a pure phase, or a member of a solid solution absent from the vertex: zero. The phase is not part of the assemblage, and the solve decides whether it enters.

A dead species (one whose component the budget lacks) stays at zero.

No member of a phase starts below ϵ, the floor of the search, as none does in a state as given. The potentials of the program know nothing of a trace species: at the pH of a cement they give H⁺ an activity near 1e-100, far below anything the search resolves, and a logarithmic iteration climbs back from such an amount slowly if at all.

The amounts are those of an ideal phase at the potentials of the program, which is what the dual solve needs to start from; they do not meet the budget exactly, which it does not need either, b being fixed for the whole call.

The first form gave every absent species aⱼ mol whatever its phase, capped at one mole. On a cement that is dozens of species near a mole, and pure phases the program found undersaturated among them. Measured from the cast state of four cements, the start was then 11 to 66 mol off a 4-mol budget, and on two of them nothing certified from it, so the search lost 9 and 15 s before the state as given certified. Placed in its phase, the start is 1.3 to 1.8 mol off, and all four certify from it at the first attempt, in 0.1 to 0.2 s.

ChemistryLab._mixes_on_sites Method
julia
_ideal_mixing_system(cs) -> Union{ChemicalSystem, Nothing}

cs with every sublattice phase replaced by ideal mixing of the same end-members, the species, their order and every other phase unchanged; nothing when cs has no sublattice phase.

ChemistryLab._refuse_infeasible_budget Method
julia
_refuse_infeasible_budget(des, state, bfix, lp, ϵ) -> (state, certificate)

The answer to a budget no amounts of the declared species can meet: a copy of the state as given, and its certificate with optimal = false, budget_feasible = false, route = :infeasible and unplaceable, the reason in words, read from the Farkas vector of the linear program.

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.

The receiver is emptied first

The two instances are copies of the same end-members, so moving the receiver's whole content back into the donor leaves the budget as it was, and the receiver then holds the trial composition alone. Without it the trial material mixes with what the receiver already holds. Measured on two instances of a carbonate binary that both sat at x = 0.117: the receiver came out at 0.209 instead of the 0.973 asked for, no pass reached the pair, and which instance was the donor was decided by the rounding of two equal amounts. With the receiver emptied, the seed is the same whichever instance gives, and the first pass certifies.

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, lp_start, dual,
                      fallback_model, fallback_on) -> (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.

state is not modified: the answer is a new state, and the same state can be solved again, or under another model, without a copy.

The linear program over the pure phases comes first

With every activity set to one, the equilibrium is a linear program: minimize gᵀn subject to A n = b, n ≥ 0. It is solved before anything else (OptimaSolver's lp_start), and it does two things.

  • It refuses an impossible budget. When no non-negative amounts of the declared species meet b, the program proves it with a combination of the balances that every species raises and the budget lowers. The state is then returned as given, with optimal = false, budget_feasible = false, route = :infeasible and unplaceable, the reason in words; nothing is searched, and STRICT_CONVERGENCE raises.

  • It gives the first start. Its vertex, with every species it leaves out given the amount its multipliers give it in its own phase, is where the search begins: on two cold cements, 0.05 s and 0.26 s against 11.3 s and 10.1 s without it, to the same composition.

lp_start = false skips it, and so does autostart = false. The certificate reports budget_feasible, the start the answer came from as route (:lp_start, :state, :ideal_mixing, :ideal, :continuation or :restart) and the number of dual solves the search ran, n_dual_solves.

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.

Options of the solvers

dual is a NamedTuple of keywords for the DualEquilibriumSolver the route certifies with — maxit, max_active_updates, inner_tol, inner_maxit, tol, si_tol — for instance dual = (; maxit = 1000). tol and si_tol are also the thresholds of the certificate, so loosening them loosens the proof. The other keywords go to the interior-point starts: variable_space sets their formulation, Ipopt takes the common arguments of Optimization.jl (maxiters, reltol, maxtime, verbose), which OptimizationIpopt maps to its options, and OptimaOptimizer ignores the rest.

certificate.optimal == true is a proof of a global minimum when the log activities are the gradient of one Gibbs energy and that energy is proved convex over the whole domain: ideal mixing, convex solid solutions, and the aqueous models of a convexity bound (the ideal dilute model, Debye–Hückel with one ion size and Davies on ions of charge two at most). certificate.scope says which case holds: :global_minimum; :local_minimum when convexity is not proved (a concave mixing phase, Pitzer, a convexity bound above one) and the Hessian over the directions that conserve matter is positive definite at the answer, its smallest eigenvalue being reduced_curvature; :kkt_point when that curvature is not positive, or a constraint leaves the sufficiency unestablished; :self_consistent when the activities are not the gradient of one energy, as with the B-dot and Davies models at their default settings, and the answer is a speciation consistent with its own activities.

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.

A second activity model, when the first does not hold

fallback_model names an activity model to use when model does not give a certified answer (fallback_on = :refusal, the default), or also when its certified answer lies past the ionic strength model is stated valid for (fallback_on = :out_of_range, see activity_model_range). The fallback is solved from the first answer when it certified, from state otherwise. When it certifies, its answer is returned, a warning names both models and the reason, and the certificate carries fallback_used = true, fallback_reason and the first certificate as primary_certificate. When it does not, the first answer is returned with fallback_used = false. Without a fallback, nothing changes.

julia
eq, cert = equilibrate_certified(state; model = cemdata18_activity_model(),
                                 fallback_model = HKFActivityModel(),
                                 fallback_on = :out_of_range)

A phase that may unmix: instances = :auto

A solid solution declared SolidSolutionPhase(...; instances = :auto) is solved with one composition. When the certificate finds it wanting to split, the system is rebuilt with a second instance (with_instances), the answer is carried over, and the passes of equilibrate_split look for the pair. Their answer is kept when its KKT error is smaller: it is then a state of the enlarged system, and the certificate names the phases in instances_added. Under STRICT_CONVERGENCE the first solve is a search and the final answer is judged.

Example

julia
using ChemistryLab, OptimaSolver
eq, cert = equilibrate_certified(state)
cert.optimal          # true — proved globally optimal
cert.balance          # element balance residual, in mol
cert.balance_relative # the same, relative to what each row holds
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 or the phase given two instances. Then 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 (Michelsen, 1982), 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.

Under STRICT_CONVERGENCE the passes are searches, not results: they run with strictness suspended, and the answer they end on is judged strictly — an error if it does not certify. Honored by the first pass, which is by construction the one expected not to certify, the flag would make the function raise before it had seeded anything.

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.

ChemistryLab.with_instances Method
julia
with_instances(state, cs) -> ChemicalState

state carried over to cs, a system with_instances built from its own: every species keeps its amount, and the copies it did not have start at zero, so the element budget is unchanged. Temperature and pressure are kept.

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 the α_max = w/c / 0.42 of Powers (1948) 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
co    = literature_row("BaroghelBouny1999", "retention_fit", "CO")
r     = VanGenuchten(; a = co.a, m = 1 / co.b)      # 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 barely depend on pressure, and not at all in most records, 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

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 forward-mode differentiation.

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 the minerals do not depend on pressure at all — V⁰(1 bar) = V⁰(100 bar) exactly for Cal — water compresses by about 4.5e-10 per pascal, and 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 5e-5 at 1 bar, meaning some 200 bar to change the volume by one percent, beyond what a Newton step on that residual can take before the pressure leaves the domain of the equations of state.

Which is the physics, not a solver defect: the volume of a nearly incompressible condensed system is all but fixed by its composition, and pressure has little 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; a real gas phase counts for Z N R T / P instead (peng_robinson).

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. Baroghel-Bouny et al. (1999), Table 5, read from data/literature/BaroghelBouny1999.json with the water-to-cement ratios of their Table 1 (CO and CH are pastes, BO and BH concretes, CH and BH contain 10 % silica fume):

mixW/Ca (MPa)bm = 1/b
BO0.4818.62372.27480.4396
BH0.2646.93642.06010.48541
CO0.3437.54792.16840.46117
CH0.1996.28371.9540.51177

Examples

julia
# The ordinary cement paste of Baroghel-Bouny et al. (1999), their mix CO
co = literature_row("BaroghelBouny1999", "retention_fit", "CO")
r = VanGenuchten(; a = co.a, m = 1 / co.b)

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
fit = literature_row("BaroghelBouny1999", "retention_fit", "CO")
co = VanGenuchten(; a = fit.a, m = 1 / fit.b)
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 :: Threads.Atomic{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. The caller does this, not the package: nothing in src/ resets or reports the counter, so a figure read without resetting first covers the whole session and not the run. (This docstring claimed integrate did it; it never has.)

It is an Atomic rather than a Ref only so that increments cannot be lost when solves run on several tasks, which += 1 on a Ref does not guarantee. Reading and writing it are unchanged.

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 does not by itself indicate whether a point is the minimum 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.

This is the session default, set by the caller and read by every solve. Internally the package sometimes has to suspend it — while it is computing a starting point rather than an answer — and it does so through _relaxed_convergence, which is scoped to one task and leaves this Ref untouched. Read the effective value with _strict_convergence, never this Ref directly.

ChemistryLab._AUTO_SPLIT Constant
julia
_AUTO_SPLIT

Whether a phase declared instances = :auto may be given its second instance during the solves of this task: true, except inside a solve that builds a starting point. Such a solve (the ideal pre-solve, the ideal-mixing pre-solve) hands its answer to a search in the system it was given, so its answer must stay in that system. Given the second instance, it returned a state with two more species than the search it was meant for, and the dual solve built from both indexed past the end of the conservation matrix.

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.

Scoped to the task that set it, for the reason given under _STRICT_OVERRIDE — and here the leak was the quieter of the two: a flag left set silences the non-convergence warnings of every later solve in the session, so a genuine failure stops announcing itself.

ChemistryLab._POLISH Constant
julia
_POLISH

Whether the answer of a back end (Ipopt, the interior point of OptimaSolver, any optimizer behind an EquilibriumSolver) is polished by the dual Newton of the system before it is returned: true unless suspended.

A back end minimizes a scalar, n⋅μ(n), or steers on a gradient, and the dual Newton solves the conditions of equilibrium themselves, μ(n) = −Aᵀy on the species present, with the activities of the model whatever they derive from. The two coincide when the activities are the gradient of one Gibbs energy, homogeneous of degree one: then n⋅μ(n) is that energy and its gradient is μ. For the extended Debye–Hückel models in general use they are not, and the minimum of n⋅μ(n) is another composition than the equilibrium. Polished, every route returns the composition its certificate describes, and the derivatives lifted at it are those of the map it returns.

It is the default of the polish keyword of a back end's solve. The searches of this package, which ask a back end only for a starting point they polish themselves, pass polish = false rather than run the solve under _unpolished: a scoped value around a solve is inferred through, and on the first cement equilibrium of a session that cost seconds of compilation.

ChemistryLab._SOLVER_FACTORIES Constant
julia
_SOLVER_FACTORIES

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

The search tries them in this order, and the order is not neutral: a start that certifies ends the search, so which one comes first decides both the time and, on a hard case, the verdict. Putting the default back end first was measured and not kept. On the documentation's cements it saved about an eighth of the time overall, but it slowed two pages down and it lost the certificate of a composite binder at seven days (element balance 1.0 against 7.4e-13), which the search starting from Ipopt finds.

ChemistryLab._STRICT_OVERRIDE Constant
julia
_STRICT_OVERRIDE

A per-task override of STRICT_CONVERGENCE, nothing when none is in force.

It exists because the two readers of that setting want different things. A caller sets it once for a session; the package suspends it for the duration of a start search. Saving the Ref, writing it, and restoring it in a finally does the second correctly only while no two of them overlap — and two solves on two tasks do overlap. Measured, with A entering first and leaving first: B lost the relaxation inside its own region, and the flag was left set after both had finished, so every later solve in that session inherited it.

A ScopedValue has exactly the semantics the finally was imitating — dynamic extent, inherited by child tasks, invisible to siblings — and cannot be left behind, because there is nothing to restore.

ChemistryLab._STRIP_TAGS Constant
julia
_STRIP_TAGS

The tags of the dual numbers a certified search is run without, outermost first.

A derivative is taken by solving on the values and lifting the answer by the implicit-function theorem (_lift_equilibrium). The values have to be those of every input the duals of the differentiation reach: the state and the budget, which are stripped before the search, but also the thermodynamic data, the parameters of the activity and mixing models, and the targets of a constraint, which live in the species and in closures the search builds as it goes. While a search runs in this scope, _build_params, the activity closures and the blocks of a constraint return their values without these tags. A tuple and not one tag, because a differentiation nested in another strips its own level on top of the outer one's.

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, in the number type of n and p.

  • 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...).

A solver is not meant to be used by two threads at once: its optimizer may keep the last answer as the next start (OptimaOptimizer). Give each thread its own.

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

  • ϵa: the floor of the activities, min(ϵ, _activity_floor_amount(state)), the molality _ACTIVITY_FLOOR times the state's water.

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;
                         pinned = falses(n), pinnable = falses(n), 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.

Which species are held at zero

pinned names the species the caller knows to be absent, the pure phases at their bound. A member of a mixing phase is never absent while the phase exists, and an aqueous trace in particular keeps its conservation law however small its amount: pinning it would erase the perturbation of the component it carries.

No back-end returns z, so a pure phase near its bound may still be left free. Among the species marked pinnable — the pure phases — one that is negligible on the scale of the system yet takes a leading share of the response is pinned, and the system solved again; on calcite and CO₂ in water with a gas phase declared, the absent gas left free took the whole perturbation, satisfying A ndot = bdot to 4e-16 and meaning nothing. Each pass pins at least one species, so the loop terminates.

How the system is solved

A trace species has a curvature near 1/n, up to 1e300, beside a conservation block of order one, so the matrix is equilibrated by rows and columns before a rank-revealing solve; a pseudo-inverse of the unscaled matrix discards the conservation equations silently. H is singular by construction, a pure phase having unit activity and hence a zero row, which the saddle-point form handles and a method inverting H does not. Stationarity and the imposed component perturbation are checked on the answer, and a failure of either raises rather than returning a sensitivity that does not solve the problem.

ChemistryLab._exploring_starts Method
julia
_exploring_starts(f)

Run f with _EXPLORING_STARTS set, for this task only. Nesting is safe, and so is overlapping with another task's search.

ChemistryLab._finish_backend_solve Method
julia
_finish_backend_solve(esolver, state, eq; ϵ, b = nothing, certificate = nothing,
                      polish = _POLISH[]) -> ChemicalState

The answer eq a back end returned from state, polished by the dual Newton when polish holds (its default is _POLISH), OptimaSolver is loaded and the system has an aqueous phase with H2O@, and returned as it is otherwise. certificate, a Ref, receives the certificate of the answer returned, or nothing when none was computed.

A polish that does not certify keeps the better of the two answers, ranked as the certified search ranks them, and says so as a non-converged solve does: a warning, or an error under STRICT_CONVERGENCE.

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

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

ChemistryLab._pure_phase_indices Method
julia
_pure_phase_indices(cs) -> Vector{Int}

The species of cs the dual solver treats as pure phases: every species outside the aqueous phase that belongs to no solid solution and no site family. These are the ones a sensitivity may hold at zero when absent.

ChemistryLab._refuse_state_keywords Method
julia
_refuse_state_keywords(kwargs, caller)

Throw if kwargs names a temperature or a pressure.

The solves take the temperature and the pressure from the ChemicalState, and forward the keywords they do not know to the optimizer's options. So a T = 293.15u"K" handed to one of them was not an error: it reached the optimizer, was ignored there, and the solve ran at the state's own temperature without a word. On a C-S-H pore solution measured at 20 °C that is a 0.17 difference in pH — the change in pKw between 20 and 25 °C — reported as if it were the model's.

ChemistryLab._relaxed_convergence Method
julia
_relaxed_convergence(f)

Run f with strict convergence suspended, for this task only.

Used wherever the package is computing a starting point rather than a result: a back end that reports MaxIters on a candidate must not raise, or the catch around it silently loses that candidate and a caller asking for strict results gets a worse search than one who did not.

ChemistryLab._require_gibbs_duhem Method
julia
_require_gibbs_duhem(esolver, system, p)

Refuse a back end that minimizes n⋅μ(n) on a model whose activities do not satisfy the Gibbs–Duhem relation, when no dual Newton is there to polish its answer.

The gradient of n⋅μ(n) is μ + Jᵀn, with J = ∂μ/∂n, and Jᵀn = 0 is the Gibbs–Duhem relation, Σᵢ nᵢ ∂μᵢ/∂nⱼ = 0. Where it holds, the minimum of n⋅μ(n) is the equilibrium; where it fails, it is another composition, which nothing would then correct. It is checked at a composition holding every species, a kilogram of solvent with a tenth of a mole of each solute and one mole of everything else, where the terms that break it are all present.

ChemistryLab._scaled_kkt_solve Method
julia
_scaled_kkt_solve(A, H, gθ, bdot, free) -> Vector

The species part of the solution of the sensitivity system restricted to the species free, equilibrated by rows and columns and solved by a truncated singular value decomposition; raises when the answer does not satisfy the equations to 1e-8 of their scale.

ChemistryLab._solve_dual Method
julia
_solve_dual(esolver, state, ϵ; b = nothing) -> ChemicalState

Equilibrium of a problem carrying dual numbers: in the state (its amounts, its temperature, its pressure), in the budget b, in the standard potentials of its species or in the parameters of its activity model.

No optimization solver is asked to iterate on dual numbers — most cannot, and Ipopt, a C library, does not operate on them. The equilibrium is solved by esolver on the values of the outermost level of duals, polished by the dual Newton (see _POLISH), and the answer is lifted by the implicit-function theorem at it, as the certified route lifts its own (_lift_equilibrium): exact at every level of a nested differentiation, with the active set read off the answer. The conditions lifted are those of the dual Newton, μ(n) = −Aᵀy on the species present, and the polish is what makes them the conditions the returned answer satisfies; lifted at an unpolished answer of a back end that minimized n⋅μ(n), they described another map than the one it returned. A pure phase holding no more than the certificate's floor is absent (10ϵ at an answer the polish was suspended for).

Without the certified solver of the system (OptimaSolver not loaded, or no aqueous phase), only the amounts of the state may carry duals, one level, and the derivative is that of the optimality conditions of the unconstrained problem (_attach_sensitivity).

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._strict_convergence Method
julia
_strict_convergence() -> Bool

The effective strict-convergence setting: the innermost _STRICT_OVERRIDE in force, or the session's STRICT_CONVERGENCE when there is none.

ChemistryLab._unpolished Method
julia
_unpolished(f)

Run f with _POLISH off, for this task only: every back-end solve inside it returns its own answer. For one solve, pass polish = false instead.

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 accepts; its scope says what that proves (a global minimum, a KKT point, or a self-consistent speciation):

julia
state_eq = equilibrate(state)                  # certified where it can be
state_eq = equilibrate(state; certify = false) # single back end, polished
cert = Ref{Any}()
state_eq = equilibrate(state; certificate = cert)   # cert[] === nothing if none

The certified route is the default because a single back end does not certify every case here: measured on calcite dissolving in water, the interior point returns a composition whose charge balance is met to 3 % only, 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, or pass a Ref as certificate, which receives it, or nothing when no certificate was computed. certify has three values:

  • nothing (the default): the certified search where it applies — OptimaSolver loaded, an aqueous phase with H2O@ — and the single back end otherwise;

  • true: the certified search, an ArgumentError where it does not apply, and an UncertifiedEquilibrium where no route certifies, rather than an answer that was never certified;

  • false: the single back end, its answer polished by the dual Newton where the certified search would apply.

Every back end's answer is polished in that way, the two-argument form's included, so that all routes return a composition satisfying the same conditions of equilibrium (see ChemistryLab._POLISH). Without OptimaSolver, a back end that minimizes n⋅μ(n), such as Ipopt, refuses an activity model that breaks the Gibbs–Duhem relation, since its minimum is then not the equilibrium.

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

  • certificate: a Ref that receives the certificate of the answer returned, or nothing when none was computed.

  • kwargs...: forwarded to the underlying solver. The temperature and the pressure are not among them: they are the state's, and a T or a P given here is refused rather than passed on and ignored.

state is not modified: the answer is a new state.

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.