Kinetics API

Chemical kinetics module: rate models, kinetic reactions, ODE problem setup, solvers, and calorimetry.

Rate constants and models

ChemistryLab.KINETICS_RATE_FACTORIES Constant

KINETICS_RATE_FACTORIES

Compiled ThermoFactory objects for each kinetic rate model. Populated by __init__() from KINETICS_RATE_MODELS.

Keys are model name symbols (e.g. :arrhenius). Values are ThermoFactory callables that return SymbolicFunc{1} instances.

Usage

factory = KINETICS_RATE_FACTORIES[:arrhenius]
k = factory(; k₀=1e-5, Ea=50000.0, T_ref=298.15, R_gas=R_GAS)
k(; T = 298.15)   # → 1e-5  (rate constant at reference temperature)

ChemistryLab.KINETICS_RATE_MODELS Constant

KINETICS_RATE_MODELS

Dictionary of raw kinetic rate-constant model expressions, analogous to THERMO_MODELS.

Each entry maps a model name (:arrhenius, …) to a Dict containing:

• :k — symbolic Expr for the rate constant as a function of variables.

• :vars — list of variable symbols (e.g. [:T]).

• :units — list of Symbol => Quantity pairs for parameters and variables.

• :output_unit — Quantity representing the output unit.

At package initialization, every entry is compiled into a ThermoFactory stored in KINETICS_RATE_FACTORIES.

Example

# The acid mechanism of calcite dissolution, Palandri & Kharaka (2004), Table 33
calcite = literature_row("PalandriKharaka2004", "carbonate_rates", "calcite")
k_acid = KINETICS_RATE_FACTORIES[:arrhenius](;
    k₀    = 10.0^calcite.acid_log_k,           # mol/(m² s) at T_ref
    Ea    = ustrip(us"J/mol", calcite.acid_E),  # J/mol
    T_ref = 298.15,                            # K
)
k_acid(; T = 310.0)   # → Float64 rate constant

ChemistryLab.PK84_PARAMS_C2S Constant

PK84_PARAMS_C2S :: NamedTuple

Canonical parameters of Parrott and Killoh (1984) for belite (C₂S = Ca₂SiO₄).

With n₁ = 1 the Avrami branch reduces to k₁(1-ξ), which never limits the rate: belite hydration is governed by the power law throughout.

See PK84_PARAMS_C3S.

ChemistryLab.PK84_PARAMS_C3A Constant

PK84_PARAMS_C3A :: NamedTuple

Canonical parameters of Parrott and Killoh (1984) for tricalcium aluminate (C₃A = Ca₃Al₂O₆). See PK84_PARAMS_C3S.

ChemistryLab.PK84_PARAMS_C3S Constant

PK84_PARAMS_C3S :: NamedTuple

Canonical parameters of Parrott and Killoh (1984) for alite (C₃S = Ca₃SiO₅), with keys k₁, n₁, k₂, k₃, n₃, Ea, T_ref. They are valid for the Blaine fineness PK_BLAINE_REF and a reference temperature of 20 °C.

The rate constants and exponents are those of Table 3 of Lavergne et al. (2018), the activation energy that of their Table 4, both read from data/literature/Lavergne2018.json.

Pass to parrott_killoh_avrami, not to parrott_killoh — the two use different functional forms and their parameters are not transferable.

ChemistryLab.PK84_PARAMS_C4AF Constant

PK84_PARAMS_C4AF :: NamedTuple

Canonical parameters of Parrott and Killoh (1984) for tetracalcium aluminoferrite (C₄AF = Ca₄Al₂Fe₂O₁₀). See PK84_PARAMS_C3S.

ChemistryLab.PK_AVRAMI_SEED Constant

PK_AVRAMI_SEED

Lower bound imposed on the normalized degree of hydration inside the Avrami branch of parrott_killoh_avrami, so that the rate is strictly positive at α = 0 and the ODE leaves its degenerate initial point. See the discussion in parrott_killoh_avrami.

One visible consequence: a phase whose rate is governed by the Avrami branch near α = 0 (C₃S, C₃A, C₄AF) starts more slowly than one governed by the power law (C₂S), so alite only overtakes belite after a few minutes. The crossover falls inside the induction period, which this model does not describe anyway.

ChemistryLab.PK_BLAINE_REF Constant

PK_BLAINE_REF :: Quantity

The Blaine fineness the constants of Parrott and Killoh (1984) were calibrated at, B₀ of Lavergne et al. (2018), p. 39, read from data/literature/Lavergne2018.json.

It is the default reference of blaine_factor and the one parrott_killoh_avrami applies, written once so the two cannot drift apart — a rate constant and the fineness it was measured against are one datum in two places.

ChemistryLab.PK_PARAMS_C2S Constant

PK_PARAMS_C2S :: NamedTuple

Parameters of the deprecated smoothed variant for belite (C₂S = Ca₂SiO₄). See PK_PARAMS_C3S.

ChemistryLab.PK_PARAMS_C3A Constant

PK_PARAMS_C3A :: NamedTuple

Parameters of the deprecated smoothed variant for tricalcium aluminate (C₃A = Ca₃Al₂O₆) in the presence of sulfate (gypsum), corresponding to ettringite formation. See PK_PARAMS_C3S.

ChemistryLab.PK_PARAMS_C3S Constant

PK_PARAMS_C3S :: NamedTuple

Parameters of the deprecated smoothed variant parrott_killoh for alite (C₃S = Ca₃SiO₅), with keys K₁, N₁, K₂, N₂, K₃, N₃, B, Ea, T_ref.

Their provenance is unestablished. They are read from the table smoothed_variant of data/literature/ParrottKilloh1984.json, which records them as such.

pk = parrott_killoh(PK_PARAMS_C3S, "C3S")
# or with the water availability cap of Powers (1948):
pk = parrott_killoh(PK_PARAMS_C3S, "C3S"; α_max = powers_alpha_max(w_c))

ChemistryLab.PK_PARAMS_C4AF Constant

PK_PARAMS_C4AF :: NamedTuple

Parameters of the deprecated smoothed variant for tetracalcium aluminoferrite (C₄AF = Ca₄Al₂Fe₂O₁₀). See PK_PARAMS_C3S.

ChemistryLab.WALLER_PARAMS_FLY_ASH Constant

WALLER_PARAMS_FLY_ASH :: NamedTuple

Parameters of Waller (1999) for class-F fly ash, with keys τ, n, blaine_ref, Ea, T_ref, as adjusted by Lavergne et al. (2018), p. 42 to SEM image analyses of a fly ash of assumed pozzolanic activity. The values are read from data/literature/Lavergne2018.json.

Pass to waller.

ChemistryLab.WALLER_PARAMS_SILICA_FUME Constant

WALLER_PARAMS_SILICA_FUME :: NamedTuple

Parameters of Waller (1999) applied to silica fume — identical kinetics to WALLER_PARAMS_FLY_ASH, the higher reactivity being carried by the fineness. Pass to waller the effective Blaine fineness recommended by Lavergne et al. (2018), blaine = literature_value("Lavergne2018", "blaine_silica_fume"); the BET surface of silica fume is not a Blaine fineness and must not be used here.

ChemistryLab.KineticFunc Type

KineticFunc{F, R <: NamedTuple, Q}

Compiled kinetic rate function, analogous to NumericFunc for thermodynamics.

Calling convention (positional, not keyword):

kf(T, P, t, n, lna, n_initial) -> Real   # [mol/s]

where:

• T [K], P [Pa]: temperature and pressure (plain Real or ForwardDiff.Dual).

• t [s]: current time.

• n::StateView: moles of all species (named access: n["C3S"]).

• lna::StateView: log-activities of all species.

• n_initial::StateView: initial moles, in the number type of the initial state.

• return: net dissolution rate [mol/s], positive = dissolution.

AD-compatible when the compiled closure is AD-compatible.

Examples

julia> idx = Dict("C3S" => 1);

julia> n_sv  = StateView([1.0], idx);

julia> lna_sv = StateView([0.0], idx);

julia> pk = parrott_killoh_avrami(PK84_PARAMS_C3S, "C3S");

julia> pk(293.15, 1e5, 0.0, n_sv, lna_sv, n_sv) > 0
true

ChemistryLab.PoreHumidity Type

PoreHumidity(retention, system; reference, T = temperature(reference))

The internal relative humidity of a sealed paste, computed from the composition it currently has.

Pass it as the humidity keyword of parrott_killoh_avrami or waller and the rate law stops reading a humidity imposed from outside and starts reading the one the material makes for itself. That is what closes the loop: hydration consumes water, the pore space empties, the humidity falls, and humidity_factor throttles the reaction — self-desiccation, which is what the α_max = w/c / 0.42 of Powers (1948) describes empirically and what powers_alpha_max otherwise supplies as an input.

The humidity is the water activity the retention law returns at the current degree of saturation of the pore space, S = V_liquid / (V_ref − V_solid), both volumes recomputed from the composition the rate law is handed. V_ref is the fresh paste's total volume — the same reference the two-argument porosity uses, and the same sealed-curing convention: the volume the reactions destroy stays inside as empty porosity.

This is where the arrest comes from, and it is kinetic

The Kelvin term is far too small to arrest hydration thermodynamically. Measured on a CEM I paste, imposing a water activity anywhere from 0.95 down to 0.05 leaves the equilibrium assemblage unchanged: the shift is RT ln a_w = −553 J/mol of water at a_w = 0.80, worth about 1.8 kJ per mole of alite against a hydration Gibbs energy of order −100 kJ/mol. Nulling that would need a_w ≈ 5e-6, a Kelvin radius smaller than a water molecule.

A real paste stops at 75–80 % RH because transport and nucleation stop, not because the reaction has become unfavorable. So the humidity belongs in the rate law, through humidity_factor, and CapillaryWater is what makes the water activity of the equilibrium state mean the same thing.

Examples

h = PoreHumidity(retention, cs; reference = fresh)
rxn[:rate] = parrott_killoh_avrami(PK84_PARAMS_C3S, "C3S"; humidity = h)

See also: humidity_factor, WaterRetention, CapillaryWater, powers_alpha_max.

ChemistryLab.RateMechanism Type

RateMechanism(k::AbstractFunc, p::Real, q::Real, catalysts,
              inhibitors = RateModelInhibitor[]) -> RateMechanism

Construct a RateMechanism with catalysts and, optionally, inhibitors, their numbers promoted to one type.

ChemistryLab.RateMechanism Type

struct RateMechanism{F<:AbstractFunc, T<:Real}

A single kinetic mechanism (acid/neutral/base/…) contributing to the overall mineral dissolution or precipitation rate.

The mechanism rate is:

r_mech = k(T) × [Π_catalysts aᵢ^nᵢ] × [Π_inhibitors (1 + Kⱼ aⱼ)^(-mⱼ)] × sign(1 - Ω) × |1 - Ω^p|^q

Fields

• k: rate constant as AbstractFunc (typically SymbolicFunc{1} from arrhenius_rate_constant). Called as k(; T=...).

• p: saturation exponent p in (1 - Ω^p)^q. Default 1.0.

• q: outer exponent q. Default 1.0.

• catalysts: vector of RateModelCatalyst (may be empty).

• inhibitors: vector of RateModelInhibitor (may be empty, which it is unless given).

Examples

calcite = literature_row("PalandriKharaka2004", "carbonate_rates", "calcite")
k_acid = arrhenius_rate_constant(10.0^calcite.acid_log_k, calcite.acid_E)
mech   = RateMechanism(k_acid, 1.0, 1.0, [RateModelCatalyst("H+", calcite.acid_n_H)])

ChemistryLab.RateMechanism Method

RateMechanism(k::AbstractFunc, p::Real, q::Real) -> RateMechanism

Construct a RateMechanism with no catalyst contributions.

ChemistryLab.RateModelCatalyst Type

struct RateModelCatalyst{T<:Real}

Describes the contribution of a catalyst species to a reaction mechanism rate.

The catalyst multiplies the base rate by exp(n * ln aᵢ) = aᵢ^n, where aᵢ is the activity of the catalyst species.

Fields

• species: PHREEQC-format (Parkhurst and Appelo, 2013) formula string of the catalyst species (e.g. "H+", "OH-").

• n: power exponent (dimensionless).

Examples

acid_catalyst   = RateModelCatalyst("H+",  0.5)    # ∝ a(H+)^0.5
base_catalyst   = RateModelCatalyst("OH-", 0.5)    # ∝ a(OH-)^0.5
co2_catalyst    = RateModelCatalyst("CO2", 1.0)    # ∝ a(CO2)

ChemistryLab.RateModelInhibitor Type

RateModelInhibitor(species, K; m = 1)

A species that slows a mechanism down: the mechanism's rate is multiplied by

the activity of species. The factor is one where the species is absent, and falls as where  : unlike a catalyst with a negative order, with  , which grows without bound as the species vanishes, it leaves the rate measured without the species unchanged.

Fields

• species: the species, by its symbol in the system.

• K: the inverse of the activity at which the factor is .

• m: the order where  .

Examples

inh = RateModelInhibitor("Ca+2", 2000.0)       # halves the rate at a(Ca²⁺) = 5e-4

ChemistryLab.StateView Type

StateView{T, V <: AbstractVector{T}, I <: AbstractDict}

Thin wrapper giving O(1) named access to a species data vector.

sv["C3S"] === sv.data[sv.index["C3S"]]

The index dict is built once at KineticsProblem construction; data is a plain vector (mutated in-place or re-wrapped each ODE step) — no dict allocation in the hot path.

Examples

julia> idx = Dict("Ca++" => 1, "C3S" => 2);

julia> sv = StateView([0.5, 1.0], idx);

julia> sv["C3S"]
1.0

julia> haskey(sv, "Ca++")
true

ChemistryLab._fineness_at Method

_fineness_at(β₀, p, fraction) -> Real

The fineness factor at a remaining fraction n/n₀. p === nothing returns β₀ identically — same object, no arithmetic — so a law with a frozen factor computes exactly the number it computed before this existed.

ChemistryLab._fineness_parts Method

_fineness_parts(blaine, blaine_ref) -> (β₀, p)

The fineness factor split into the part that can be computed once and the part that has to be read from the state at every evaluation.

p === nothing means the factor is frozen, which is what a bare fineness has always meant and what every rate law shipped in this package was calibrated with. A ShrinkingCoreArea instead returns its exponent, and the caller multiplies β₀ by g(n/n₀) at each evaluation.

The split exists so that the frozen path keeps costing nothing: a constant fineness is still one area_ratio for the whole integration, evaluated here.

ChemistryLab.add_kinetics_rate_model Method

add_kinetics_rate_model(name::Symbol, dict_model::Dict)

Register a new kinetic rate-constant model in KINETICS_RATE_MODELS and compile it into KINETICS_RATE_FACTORIES.

dict_model must contain at minimum :k (expression), :vars (variable list), :units (parameter units), and :output_unit.

Example

add_kinetics_rate_model(:power_law, Dict(
    :k    => :(k₀ * (T / T_ref)^n),
    :vars => [:T],
    :units => [:T => u"K", :T_ref => u"K", :k₀ => u"mol/(m^2*s)", :n => u"1"],
    :output_unit => u"mol/(m^2*s)",
))

ChemistryLab.arrhenius_rate_constant Method

arrhenius_rate_constant(k₀, Ea; T_ref=298.15, R_gas=R_GAS) -> NumericFunc

Build a temperature-dependent Arrhenius rate constant as a NumericFunc:

k(T) = k₀ × exp(-Eₐ / R × (1/T - 1/T_ref))

The returned object is callable as k(; T=...) and fully AD-compatible (ForwardDiff-safe: the closure captures k₀, Ea, T_ref, R_gas directly, so dual numbers propagate correctly through all parameters).

Arithmetic between SymbolicFunc/NumericFunc objects is supported, so rate constants can be composed with activity or surface-area functions.

Arguments

• k₀: pre-exponential factor at T_ref. Plain Real → SI [mol/(m² s)]; Quantity → automatically converted (e.g. 5e-4u"mol/(m^2*s)").

• Ea: activation energy. Plain Real → SI [J/mol]; Quantity → converted (e.g. 62.0u"kJ/mol").

• T_ref: reference temperature. Plain Real → SI [K]; Quantity → converted (e.g. 298.15u"K"). Default 298.15.

• R_gas: gas constant [J/(mol K)] (plain Real only; default R_GAS, the CODATA value taken from DynamicQuantities.Constants).

Returns

A NumericFunc with variable T (in K) and refs = (T = T_ref * u"K",).

Examples

julia> k = arrhenius_rate_constant(5.0e-4, 62000.0);

julia> isapprox(k(; T = 298.15), 5.0e-4; rtol = 1e-10)
true

julia> k(; T = 350.0) > k(; T = 298.15)   # higher T → higher k
true

Unit-aware: k₀ in mmol/(m²·s), Ea in kJ/mol, T_ref in K — all converted to SI:

k = arrhenius_rate_constant(0.5u"mmol/(m^2*s)", 62.0u"kJ/mol"; T_ref = 298.15u"K")

AD-compatible through all parameters:

ForwardDiff.derivative(T  -> arrhenius_rate_constant(5e-4, 62000.0)(; T = T),  298.15)
ForwardDiff.derivative(Ea -> arrhenius_rate_constant(5e-4, Ea)(; T = 350.0),   62000.0)
ForwardDiff.derivative(k₀ -> arrhenius_rate_constant(k₀,   62000.0)(; T = 298.15), 5e-4)

ChemistryLab.blaine_factor Method

blaine_factor(blaine; blaine_ref = PK_BLAINE_REF) -> Real

Fineness correction of the hydration rate: the parameters of Parrott and Killoh (1984) were adjusted for a cement of Blaine fineness blaine_ref, and the rate scales as blaine / blaine_ref.

Both arguments accept a DynamicQuantities.Quantity, a plain Real in m²/kg, or a BlaineSurfaceArea. The default reference is PK_BLAINE_REF for clinker phases; pass blaine_ref = WALLER_PARAMS_FLY_ASH.blaine_ref for the Waller kinetics of additions.

A BET area is not a Blaine fineness

Passing a BETSurfaceArea raises instead of returning a number. Silica fume is about 20 000 m²/kg by BET and about 2 000 m²/kg by effective Blaine, so the substitution would multiply its hydration rate by ten. A bare number is still read as a Blaine fineness, which is the historical contract; type the argument to get the check.

Examples

julia> blaine_factor(PK_BLAINE_REF)
1.0

julia> round(blaine_factor(1.2 * PK_BLAINE_REF); digits = 4)
1.2

ChemistryLab.humidity_factor Method

humidity_factor(h) -> Real

Reduction coefficient β_h applied to the hydration rate at internal relative humidity h ∈ [0, 1], as Lavergne et al. (2018), Eq. 10 write it:

Hydration is taken to stop below 80 % relative humidity, on thermodynamic grounds. The cut is a genuine discontinuity: the one-sided limit from above is β_h(0.80⁺) ≈ 0.0953 while the value at and below 0.80 is exactly 0. The jump is mild in practice, the rate having already fallen by an order of magnitude from β_h(0.99) ≈ 0.914.

Examples

julia> humidity_factor(0.75)
0.0

julia> humidity_factor(0.80)
0.0

julia> round(humidity_factor(0.801); digits = 4)
0.0968

ChemistryLab.palandri_kharaka Method

palandri_kharaka(mineral; mechanisms = (:acid, :neutral, :base), pco2 = nothing,
                 fe3 = "Fe+3", o2 = "O2@", assume_Ea = nothing) -> Vector{RateMechanism}

The dissolution mechanisms of mineral as Palandri and Kharaka (2004) tabulate them, ready for transition_state: of the mechanisms named in mechanisms, each one the report gives a rate constant for, with its Arrhenius rate constant (log k at 25 °C and the activation energy, through arrhenius_rate_constant) and its catalyst.

mechanism	rate per unit area, far from equilibrium	catalyst
:acid		"H+", n > 0
:neutral		none
:base		"H+", n < 0
:carbonate		the species pco2

Each catalyst is the activity of a species raised to its order, the product of the report's Eq. (3a). The sulfides of its Table 35 carry two more: an acid mechanism in H⁺ and Fe³⁺ (pyrite:  , its Eq. 3b), and a neutral mechanism in dissolved O₂, the species named by fe3 and o2; the order of pyrite in H⁺ is negative although measured at low pH, the report reading it as a hydroxide catalysis (pp. 3–4).

The report writes the base mechanism as an order in H⁺ that is negative, which is how it is applied here rather than converted to an order in OH⁻. Each mechanism takes the saturation term 1 − Ω (p = q = 1): the report gives other exponents for a few minerals in its text, not in the tables.

mineral is a row of the transcribed tables ("calcite", "quartz, BET surface area", "gibbsite", "anhydrite", …), listed by palandri_kharaka_minerals. Where the report gives several sets for one mineral (quartz normalized by its geometric or by its BET area, amorphous silica from two data sets), the row says which: the rate per unit area has to be used with the area it was normalized by.

The carbonate mechanism reads the partial pressure of CO₂, which a rate law sees as the activity of a gas species, P/P° (P_STANDARD): asking for :carbonate requires pco2, the symbol of that species in the system.

An activation energy the report does not give (that of the neutral mechanism of gypsum, which its text says the data could not determine) is refused, unless assume_Ea supplies one; the rate constant is then that value's, an assumption of the caller, and exact at 25 °C whatever it is. So is one it prints negative, the acid mechanism of kyanite (Table 23), whose sign the transcription's notes argue is a misprint: it is stored as printed, and used only as the caller decides.

Examples

mech = palandri_kharaka("calcite")             # acid and neutral
tst = transition_state(mech, cs, rxn, BETSurfaceArea(90.0))

See also: palandri_kharaka_minerals, RateMechanism.

ChemistryLab.palandri_kharaka_minerals Method

palandri_kharaka_minerals() -> Vector{String}

The minerals whose dissolution mechanisms palandri_kharaka builds: the rows of the tables of Palandri and Kharaka (2004) transcribed in data/literature/PalandriKharaka2004.json.

ChemistryLab.parrott_killoh Method

parrott_killoh(params::NamedTuple, mineral_name::AbstractString; α_max=1.0) -> KineticFunc

Build a smoothed three-mechanism clinker hydration rate as a KineticFunc.

Deprecated, and no longer attributed to Parrott & Killoh

Use parrott_killoh_avrami with PK84_PARAMS_C3S and siblings instead. This function is kept so that existing scripts keep running, and it warns once per session.

Why the attribution is withdrawn. The formulas below are not those of Parrott and Killoh (1984): the nucleation–growth term carries no Avrami logarithm, and K₃ — a shell-formation coefficient — sits in the diffusion expression where the canonical formulation uses K₂. Nor do the shipped parameters match any published set: N₁ = 3.3 is the canonical n₃, and the canonical k₃ = 1.1 has no counterpart at all. The primary source, Parrott and Killoh (1984), is a conference proceedings without a DOI that could not be consulted, so the attribution is retracted rather than repaired by an invented calibration.

Why it matters. With PK_PARAMS_* the diffusion branch takes over very early — measured at α/α_max = 0.003 for C₂S, 0.013 for C₃S, 0.057 for C₃A, and over the whole range for C₄AF — and governs throughout the interval a seven-day run traverses, because its prefactor 3K₃/N₃ = 0.0018 d⁻¹ is 28 times smaller than the canonical k₂ = 0.05 d⁻¹. The rate then integrates in closed form,

α(t) = α_max · [1 − (1 − √(2·K₃·t / N₃))³]

which is independent of the phase, because K₃ = 0.0024 d⁻¹ and N₃ = 4 are identical in all four parameter sets. Measured on a CEM I at w/c = 0.40 over seven days, C₃S, C₂S and C₃A all land on α ≈ 0.239, and C₄AF lower still at 0.193 — there its own nucleation-growth branch is slower than diffusion and limits instead. The weighted mean comes to 0.234. The signature is unmistakable: K₁ spans a factor of 18 across the four phases and changes almost nothing.

params must be a NamedTuple with keys K₁, N₁, K₂, N₂, K₃, N₃, B, Ea, T_ref. All dimensional values accept plain Real (SI) or DynamicQuantities.Quantity.

mineral_name is the PHREEQC formula string (e.g. "C3S") used to look up the mineral moles in the n and n_initial StateViews.

Three competing mechanisms determine the rate (Parrott and Killoh, 1984):

Mechanism	Formula
Nucleation–growth	r_NG = (K₁/N₁)(1-ξ)^N₁ / (1 + B·ξ^N₃)
Interaction	r_I = K₂(1-ξ)^N₂
Diffusion	r_D = 3K₃(1-ξ)^(2/3) / (N₃·(1-(1-ξ)^(1/3)))

The rate [mol/s] is n_initial × Aₜ × min(max(r_NG, r_I), r_D) where ξ = α / α_max is the normalized degree of hydration and Aₜ = exp(-Ea/R × (1/T - 1/T_ref)) is the Arrhenius factor.

α_max can be set to apply the water/cement ratio limit of Powers (1948): α_max = powers_alpha_max(w_c).

Returns

A KineticFunc — callable as pk(T, P, t, n::StateView, lna::StateView, n_initial::StateView) -> Real [mol/s]. AD-compatible (ForwardDiff-safe): no Float64 casts in the evaluation path.

Examples

Shown rather than run: the constructor warns by design, and a deprecation warning on stderr is not something a doctest should have to reproduce verbatim. The numerical behavior of this variant is pinned by the test suite instead — see the branch oracle in test/kinetics/test_rate_models.jl, which asserts that PK_PARAMS_* puts C3S, C2S and C3A on the diffusion branch within a few percent of hydration and that C4AF is limited by its own nucleation branch.

pk = parrott_killoh(PK_PARAMS_C3S, "C3S")
idx = Dict("C3S" => 1)
n0 = StateView([1.0], idx)
lna = StateView([0.0], idx)
pk(293.15, 1e5, 0.0, n0, lna, n0) > 0      # true

See also: PK_PARAMS_C3S, PK_PARAMS_C2S, PK_PARAMS_C3A, PK_PARAMS_C4AF.

ChemistryLab.parrott_killoh_avrami Method

parrott_killoh_avrami(params::NamedTuple, mineral_name::AbstractString;
                     α_max = 1.0, blaine = nothing, humidity = nothing,
                     w_c = nothing, H = nothing) -> KineticFunc

Build the clinker hydration rate of Parrott and Killoh (1984) in its canonical formulation, as reported by Lothenbach et al. (2008) and used by Lavergne et al. (2018).

params must be a NamedTuple with keys k₁, n₁, k₂, k₃, n₃, Ea, T_ref — see PK84_PARAMS_C3S and siblings. Dimensional values accept plain Real (SI) or DynamicQuantities.Quantity.

Three competing mechanisms limit the rate, and the slowest one wins:

Mechanism	Formula
Nucleation–growth (Avrami)	α̇₁ = (k₁/n₁)(1-ξ)(-ln(1-ξ))^(1-n₁)
Diffusion (Jander)	α̇₂ = k₂(1-ξ)^(2/3) / (1-(1-ξ)^(1/3))
Shell formation (power law)	α̇₃ = k₃(1-ξ)^n₃

so that α̇ = min(α̇₁, α̇₂, α̇₃), with ξ = α/α_max the normalized degree of hydration. The returned rate [mol/s] is n_initial × Aₜ × β_B × β_h × f × α̇, where Aₜ = exp(-Ea/R × (1/T - 1/T_ref)) is the Arrhenius factor, β_B the Blaine fineness factor (blaine_factor), β_h the relative-humidity reduction (humidity_factor) and f the water/cement factor (pk_wc_factor). Each is 1 when its keyword is nothing.

Two Parrott–Killoh variants ship with ChemistryLab

parrott_killoh implements a different, smoothed variant (min(max(r_NG, r_I), r_D) with a B-damped nucleation term) together with the parameter set of PK_PARAMS_C3S and siblings. The two are not interchangeable: their parameters are not transferable, and only parrott_killoh_avrami with PK84_PARAMS_C3S reproduces the α(t) curves published in the cement literature cited above.

With the canonical parameters, C₂S has no nucleation–growth stage and C₃S has no diffusion-controlled stage — an artifact of the 1984 fit that Parrott and Killoh (1984) acknowledged, and a convenient signature to check an implementation against.

Keyword arguments

• α_max: the water availability cap of Powers (1948) — see powers_alpha_max.

• blaine: Blaine fineness of the binder, as a Quantity, a plain Real in m²/kg, a BlaineSurfaceArea — all three frozen for the whole integration — or a ShrinkingCoreArea, which makes the factor follow the grains as they are consumed. nothing (default) means no correction. See the warning below before using the last one.

• humidity: internal relative humidity, either a constant in [0, 1] or a callable t -> h(t). nothing (default) means no correction.

• w_c: the water/cement ratio, which slows the hydration once the degree exceeds a critical fraction of it (pk_wc_factor); H the critical degree of the phase in place of that of Parrott and Killoh (1984), as Lothenbach et al. (2008) fit one per clinker phase. nothing (default) means no correction.

An evolving fineness is not a free improvement

The constants of Parrott and Killoh (1984) were fitted with β_B constant, so passing a ShrinkingCoreArea leaves this law outside the calibration it came with and it has to be recalibrated — the machinery is scripts/hydration_calibration.jl.

Worse, the extra freedom largely already exists in the law. With α_max = 1 the remaining fraction is n/n₀ = 1 - ξ, so multiplying by (1-ξ)^p turns the shell-formation branch k₃(1-ξ)^{n₃} into k₃(1-ξ)^{n₃+p}: wherever that branch is the active one, p and n₃ are the same parameter written twice, and fitting both is fitting a sum. The Jander branch has no such exponent, so p is a genuine degree of freedom only while diffusion controls. Which of the two holds over a given dataset is a measurement, and identifiability is what makes it — do not report p and n₃ from one curve without it.

Returns

A KineticFunc — callable as pk(T, P, t, n::StateView, lna::StateView, n_initial::StateView) -> Real [mol/s]. AD-compatible (ForwardDiff-safe): no Float64 casts in the evaluation path.

Examples

julia> pk = parrott_killoh_avrami(PK84_PARAMS_C3S, "C3S"; blaine = 380u"m^2/kg");

julia> idx = Dict("C3S" => 1);

julia> n0 = StateView([1.0], idx);

julia> lna = StateView([0.0], idx);

julia> pk(293.15, 1e5, 3600.0, StateView([0.9], idx), lna, n0) > 0
true

See also: PK84_PARAMS_C3S, waller, blaine_factor, humidity_factor, powers_alpha_max.

ChemistryLab.pk_wc_factor Method

pk_wc_factor(α, w_c; H = nothing) -> Real

The water/cement factor of the law of Parrott and Killoh (1984): 1 while the degree of hydration α of the phase is below a critical fraction of w_c, and beyond it

as Lothenbach and Winnefeld (2006), Section 4.1 state it, which stops the phase at    ; the printed constants being rounded (  ), it steps up by    at the threshold. With H, the critical degree of the phase that Lothenbach et al. (2008), Section 3.2 fit for each clinker phase,

The constants are read from data/literature/. The water/cement ratio is a proxy for the room left to grow in, so the factor is the long-term ceiling of the law; powers_alpha_max is another statement of the same limit.

ChemistryLab.pore_saturation Method

pore_saturation(h::PoreHumidity, n) -> Real

Degree of saturation of the pore space at composition n, V_liquid / V_pore with V_pore = V_ref − V_solid. Clamped to [0, 1].

ChemistryLab.powers_alpha_max Method

powers_alpha_max(w_c; curing = :sealed) -> Real

The upper bound of Powers (1948) on the degree of hydration set by the availability of water, α_max = min(1, w/c / k): a paste below w/c = k cannot hydrate completely, k being 0.42 sealed or 0.36 water-cured, according to curing.

The 0.42 is not a stoichiometric demand, and reading it as one leads to the wrong conclusion about what a Gibbs minimization should return. It is about 0.23 g of non-evaporable water per gram of cement — the water written into the hydrate formulae, which is a mass balance — plus about 0.19 g of gel water held in the C-S-H gel pores, which is physically present and chemically unavailable. In a sealed paste hydration stops by self-desiccation with water still in the specimen, so this bound is a statement about transport and access, not about thermodynamics: an equilibrium calculation on the same mix consumes all the clinker well below 0.42, and only runs out of water near the stoichiometric demand.

The two curing conventions

curing = :sealed is a specimen that exchanges nothing with its surroundings — the convention of porosity and of the w/c example. curing = :saturated is a specimen kept under water after setting, free to draw in what the chemical shrinkage empties; the bound is then 0.36 and a mix that would arrest sealed can go on reacting. Neither is a property of the cement: they are two boundary conditions on the same paste, and which one applies is the caller's to state.

Where water is abundant — w/c above the coefficient, or a cure that keeps supplying it — the bound is 1 and this function stops doing anything, which is the correct answer rather than a degenerate case: nothing about water is then limiting the reaction.

Pass the result as the α_max keyword of parrott_killoh, parrott_killoh_avrami or waller.

It is a ceiling, not a schedule

α_max says how far the reaction can go, never how far it has got. At an early age the degree of reaction is set by the kinetics and is far below this bound; the bound binds only at long times, and only for the constituents whose kinetics would otherwise have taken them past it. For a constituent that reacts slowly — a slag, and a fly ash still more — the binding limit at 28 days is its own dissolution rate, not the water. Take the smaller of the two.

Examples

julia> powers_alpha_max(0.5)
1.0

julia> round(powers_alpha_max(0.32); digits = 4)
0.7619

julia> round(powers_alpha_max(0.32; curing = :saturated); digits = 4)
0.8889

ChemistryLab.saturation_ratio Method

saturation_ratio(stoich::AbstractVector, lna::AbstractVector,
                 ΔₐG⁰overRT::AbstractVector; ϵ=1e-16) -> Real

Compute the saturation ratio Ω = IAP / K for a kinetic reaction.

ln Ω = Σᵢ νᵢ ln aᵢ − ln K
     = Σᵢ νᵢ ln aᵢ + ΔᵣG⁰/(RT)   (note: ln K = −ΔᵣG⁰/RT = −Σᵢ νᵢ ΔₐG⁰ᵢ/RT)

where stoich[i] is the stoichiometric coefficient (positive for products, negative for reactants), lna[i] is the log-activity of species i, and ΔₐG⁰overRT[i] is the dimensionless standard Gibbs energy of formation ΔₐG⁰ᵢ / RT for species i.

Arguments

• stoich: stoichiometric coefficient vector for this reaction (length = number of species).

• lna: log-activity vector (same indexing as species in system).

• ΔₐG⁰overRT: dimensionless standard Gibbs energies ΔₐG⁰ᵢ/RT.

• ϵ: floor to avoid exp overflow when Ω → ∞.

Returns

Ω = exp(ln_IAP - ln_K) where ln_K = -ΔᵣG⁰/RT.

AD-compatible (ForwardDiff-safe).

ChemistryLab.snellings2013_glass Method

snellings2013_glass(oxides; Ea, T_ref = 293.15, modifiers = ("CaO",),
                    inhibitors = RateModelInhibitor{Float64}[],
                    extrapolate = false) -> RateMechanism

The dissolution of a calcium aluminosilicate glass far from equilibrium at pH 13, as Snellings (2013) measured it on six synthetic glasses from the composition of a slag to silica: per unit BET area and per mole of the glass's cations, at 20 °C,

the moles of each cation in the analysis oxides (mass fractions over the molar masses of the library), and  ,    the regression printed on (Snellings, 2013; Fig. 8), read from data/literature/Snellings2013.json. The abscissa is that of the figure, the ratio above with aluminum counted twice, as the compositions of the six glasses of (Snellings, 2013; Table I) give it; the regression holds with that ratio.

What Snellings (2013) covers, and what it does not:

• Composition. The glasses hold CaO, Al₂O₃ and SiO₂ only, so only the oxides named in modifiers enter the numerator: CaO by default; the MgO, Na₂O or K₂O of a real slag are left out unless named, as an assumption of the caller. runs from 0 (silica) to 0.556 (the most calcic glass); outside, the constant is refused unless extrapolate = true.

• Temperature. The rates were measured at 20 °C only: the activation energy Ea (J/mol, or a quantity) has no default, and the constant is arrhenius_rate_constant referred to T_ref (K), 20 °C.

• pH. At 13 only; Snellings (2013) gives no dependence on pH.

• Solution. In NaOH alone. Dissolved calcium slows every glass and aluminum the glasses whose calcium only balances their aluminum ((Snellings, 2013; Table II)); the paper gives the measurements, not a law, so inhibitors (RateModelInhibitor) are the caller's. Dissolved silicon up to 11 mM, and the affinity above  , change nothing measurable.

The rate is per mole of cations: the glass species is glass_species with M =cation_molar_mass(oxides), and its area a BETSurfaceArea with the material's BET surface.

Examples

slag = Dict("CaO" => 0.43, "Al2O3" => 0.19, "SiO2" => 0.38)    # glass G1
mech = snellings2013_glass(slag; Ea = 60.0e3)                 # Ea: the caller's
log10(mech.k(; T = 293.15))                                   # 2.74 x - 8.47

See also: transition_state, glass_species.

ChemistryLab.waller Method

waller(params::NamedTuple, mineral_name::AbstractString;
       α_max = 1.0, blaine = nothing, humidity = nothing) -> KineticFunc

Build the reaction rate of a pozzolanic or latent-hydraulic addition (fly ash, silica fume, ground granulated slag) after Waller (1999), as a KineticFunc.

The degree of reaction follows a sigmoid in log-time,

whose rate, written as a function of the current degree so that it composes with temperature, fineness and humidity corrections, is

params must be a NamedTuple with keys τ, n, Ea, T_ref and, optionally, blaine_ref — see WALLER_PARAMS_FLY_ASH.

Pozzolanic reactions are markedly more temperature-sensitive than the hydraulic reactions of clinker: the shipped activation energy is higher than that of any clinker phase in PK84_PARAMS_C3S and its siblings.

Silica fume is far finer than the cement, by an order of magnitude when its surface is measured by BET. Its reactivity is nonetheless represented here through an effective Blaine fineness (WALLER_PARAMS_SILICA_FUME) — the two measurements probe different physical phenomena and are not interchangeable.

Keyword arguments

Identical to parrott_killoh_avrami. The Blaine correction is taken relative to params.blaine_ref, the fineness of the fly ash the kinetics were adjusted to, not to the clinker reference of PK_BLAINE_REF.

Here an evolving area is a real degree of freedom

The degeneracy warned about on parrott_killoh_avrami does not carry over in the same form. The sigmoid rate is (n/τ)(1-ξ)^{1+1/n} ξ^{1-1/n}, and n sets both exponents at once, in opposite directions; a ShrinkingCoreArea exponent p shifts only the (1-ξ) one. So p is not a rewriting of n — it is a separate shape, correlated with n but distinguishable. The calibration caveat still applies in full: these constants were fitted with β_B frozen.

Examples

julia> fa = waller(WALLER_PARAMS_FLY_ASH, "FlyAsh");

julia> idx = Dict("FlyAsh" => 1);

julia> n0 = StateView([1.0], idx);

julia> lna = StateView([0.0], idx);

julia> fa(293.15, 1e5, 86400.0, StateView([0.9], idx), lna, n0) > 0
true

No parameter set is shipped for a slag: give its characteristic time from a source of your own, waller(merge(WALLER_PARAMS_FLY_ASH, (τ = τ_slag,)), "GGBS").

See also: WALLER_PARAMS_FLY_ASH, WALLER_PARAMS_SILICA_FUME, parrott_killoh_avrami.

Kinetic reactions

The area models the rate factories consume live on their own page, since they are shared with surface chemistry: see Surfaces API.

ChemistryLab.KineticReaction Type

struct KineticReaction{R<:AbstractReaction, F, H}

Associates a chemical Reaction with a compiled KineticFunc.

Following Leal et al. (2017), reactions — not individual species — carry kinetics. A single mineral can therefore appear as a reactant in multiple KineticReaction objects (e.g. C₃A → ettringite and C₃A → monosulphate for multi-pathway cement hydration). The ODE state is indexed by unique mineral species, and contributions from all reactions that consume the same mineral are accumulated.

Fields

• reaction: the underlying Reaction / CemReaction.

• rate_fn: a KineticFunc (or any callable matching the six-argument signature (T, P, t, n, lna, n_initial) -> Real) computing r [mol/s].

• idx_mineral: index of the primary (controlling) mineral species in the parent ChemicalSystem. Determined automatically as the first solid (AS_CRYSTAL) reactant.

• stoich: stoichiometric coefficient vector for all species in the system. Sign convention: positive for products, negative for reactants.

• heat_per_mol: enthalpy of reaction [J/mol], positive = exothermic (heat released). When nothing (default), the enthalpy is derived from the stoichiometric sum of species :ΔₐH⁰ values.

Constructors

From a species name (convenience, builds a minimal dissolution Reaction):

pk = parrott_killoh(PK_PARAMS_C3S, "C3S")
kr = KineticReaction(cs, "C3S", pk)
kr = KineticReaction(cs, "C3S", pk; heat_per_mol = 114_634.0)

From an explicit Reaction (multi-pathway):

pk_c3a = parrott_killoh(PK_PARAMS_C3A, "C3A")
kr_ett  = KineticReaction(cs, rxn_C3A_ettringite,   pk_c3a)
kr_mono = KineticReaction(cs, rxn_C3A_monosulphate, pk_c3a)

Reaction-centric (rate stored in rxn.properties[:rate]):

rxn[:rate] = parrott_killoh(PK_PARAMS_C3S, "C3S")
kr = KineticReaction(cs, rxn)

ChemistryLab.KineticReaction Method

KineticReaction(cs::ChemicalSystem, rxn::AbstractReaction) -> KineticReaction

Reaction-centric constructor: build a KineticReaction from a Reaction that carries its kinetics in reaction.properties.

Required property:

• rxn[:rate] — a KineticFunc or any callable matching (T, P, t, n, lna, n_initial) -> Real. Non-KineticFunc callables are wrapped automatically in a KineticFunc with empty refs.

Optional property:

• rxn[:heat_per_mol] — a Number giving the molar enthalpy [J/mol] for calorimetry.

Examples

pk = parrott_killoh(PK_PARAMS_C3S, "C3S")
rxn[:rate]         = pk
rxn[:heat_per_mol] = 114_634.0
kr = KineticReaction(cs, rxn)

# Build problem directly from a list of annotated Reaction objects:
kp = KineticsProblem(cs, [rxn_C3S, rxn_C3A, rxn_C2S], state0, tspan)

ChemistryLab.KineticReaction Method

KineticReaction(cs::ChemicalSystem, species_name::AbstractString, rate_fn;
                heat_per_mol=nothing) -> KineticReaction

Convenience constructor: look up species_name in cs and build a minimal dissolution Reaction (species as sole reactant, no products) automatically.

The default stoichiometry places -1.0 at the mineral index and 0.0 everywhere else.

ChemistryLab.KineticReaction Method

KineticReaction(rxn, rate_fn, idx_mineral, stoich; heat_per_mol=nothing)

Low-level constructor: explicit Reaction, rate callable, index, and stoichiometry.

ChemistryLab.KineticReaction Method

KineticReaction(cs::ChemicalSystem, rxn::AbstractReaction, rate_fn;
                heat_per_mol=nothing) -> KineticReaction

Construct a KineticReaction from an explicit Reaction object.

The controlling mineral index (idx_mineral) is determined automatically as the index of the first solid (AS_CRYSTAL) reactant found in rxn.reactants that is present in cs. The stoichiometric vector is derived from the reaction stoichiometry.

This constructor is the recommended entry point for multi-pathway kinetics:

pk_c3a = parrott_killoh(PK_PARAMS_C3A, "C3A")
kr_ett  = KineticReaction(cs, cs.dict_reactions["C3A_ettringite"],   pk_c3a)
kr_mono = KineticReaction(cs, cs.dict_reactions["C3A_monosulphate"], pk_c3a)
kp = KineticsProblem(cs, [kr_C3S, kr_ett, kr_mono], state0, tspan)

ChemistryLab.first_order_rate Method

first_order_rate(k, cs, rxn, surface_model; p=1.0, q=1.0, ϵ=1e-16) -> KineticFunc

Build a single-mechanism first-order TST rate as a KineticFunc.

r = A(n) × k(T) × sign(1 - Ω) × |1 - Ω^p|^q

This is a convenience wrapper around transition_state with one no-catalyst mechanism. Useful as a minimal test case or for empirical fits.

Arguments

• k: rate constant as an AbstractFunc (e.g. from arrhenius_rate_constant).

• cs, rxn, surface: same as transition_state.

• p, q: saturation exponents (defaults 1.0).

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

Examples

k = arrhenius_rate_constant(1e-7, 40000.0)
rf = first_order_rate(k, cs, rxn, SurfaceSupport("calcite", "Cal", BETSurfaceArea(90.0)))
kr = KineticReaction(cs, rxn, rf)

ChemistryLab.molar_mass Method

molar_mass(kr::KineticReaction) -> Float64

Return the molar mass of the mineral species [kg/mol], used by every specific-area model to turn an amount into an area.

Searches kr.reaction.reactants for a species with an :M property, and raises when none has one: a default here is worse than a refusal, because it scales the reactive area and therefore the whole rate.

ChemistryLab.saturation_ratio Method

saturation_ratio(cs::ChemicalSystem, rxn::AbstractReaction) -> Function

The saturation ratio of rxn, as a function to call inside a rate law: (T, P, lna) -> Ω, with lna the log-activities the law receives (a StateView, or any collection indexed by species symbol), T in K and P in Pa:

over every participant of rxn, products counted positive. Ω = 1 where the equilibrium solver puts the reaction, since both read the same activities and the same standard Gibbs energies: a law written as F (1 − Ω) with this Ω vanishes there (What a rate law may be).

Every participant must be a species of cs and carry a standard Gibbs energy; one missing either is refused here rather than left out of the sum.

Example

sat = saturation_ratio(cs, rxn)
law = (T, P, t, n, lna, n0) -> k * n["Cal"] * (1 - sat(T, P, lna))

See also the method on vectors, for a computation over the whole system.

ChemistryLab.transition_state Method

transition_state(mechanisms, cs, rxn, surface; ϵ=1e-16) -> KineticFunc

Build a Transition-State Theory (TST) dissolution/precipitation rate function from a list of RateMechanism objects, returning a KineticFunc.

The compiled closure captures:

• the host name and its molar mass, from the SurfaceSupport when one is given, otherwise rediscovered from rxn + cs

• the area model, evaluated at every step from the current and initial amounts, so an area that follows the microstructure needs no change here

• stoichiometry and ΔₐG⁰ callables for all aqueous species (T-dependent Ω)

The net rate [mol/s] is:

r = A(n) × Σ_m [ k_m(T) × Π_cat(aᵢ^nᵢ) × Π_inh(1 + Kⱼ aⱼ)^(-mⱼ) × (1 - Ω^p) × |1 - Ω^p|^(q-1) ]

where Ω(T) = exp(Σ νᵢ ln aᵢ + Σ νᵢ ΔₐG°ᵢ(T)/(RT)) is re-evaluated at every ODE step — correct for variable-temperature semi-adiabatic calorimetry.

Arguments

• mechanisms: vector of RateMechanism (acid, neutral, base, …).

• cs: ChemicalSystem supplying ΔₐG⁰ callables for aqueous species.

• rxn: AbstractReaction defining stoichiometry and the mineral species.

• surface: a SurfaceSupport, which names the host solid and carries its area model, or an AbstractSurfaceModel alone, in which case the host is the first solid reactant of rxn.

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

The molar mass is no longer guessed

A host species without an :M property used to fall back to 0.1 kg/mol, silently, which is wrong by up to an order of magnitude and scales the whole rate. It now raises, naming the species.

Returns

A KineticFunc callable as f(T, P, t, n::StateView, lna::StateView, n_initial::StateView) -> Real [mol/s].

AD-compatible: all operations use generic Julia arithmetic; no Float64 casts.

References

• Palandri and Kharaka (2004).

• Leal et al. (2017).

A slow surface

A reaction moving a site from one state to another, at a rate that vanishes where the equilibrium puts the surface. See A slow surface.

ChemistryLab.sorption_rate Method

sorption_rate(k, cs, rxn) -> KineticFunc

The rate of a reaction that takes one surface site from one state to another, written as an elementary step whose reverse is fixed by the equilibrium constant of the reaction:

in mol/s. n_s is the amount of the site species among the reactants, the product runs over the other reactants with their coefficients ν_i, and Ω is the saturation ratio of the reaction as written, over every participant, products counted positive. k is in s⁻¹ (the activities are dimensionless): a number, or an AbstractFunc of the temperature such as arrhenius_rate_constant.

Why this form

With site fractions for the activities of the two site states, x_s = n_s/N and x_c = n_c/N over the family's N sites, the second term is (k/K) n_c Π a_p^{ν_p}, the products other than the site: the law is the difference of a forward and a backward mass-action rate whose ratio is the equilibrium constant K of the database. Three properties follow, whatever k:

• the rate vanishes where the equilibrium solver puts the surface, Ω = 1, because both read the same activities and the same standard energies;

• it has the sign of the affinity 𝒜 = −RT ln Ω, so 𝒜 r ≥ 0;

• two reactions closing a cycle on one family have rate constants whose ratios multiply to the product of their K, which holds by construction.

A non-ideal site family (an exchange convention, a surface potential) puts its activity coefficients in Ω, hence in the backward term only. Thermodynamics fixes the ratio of the two terms, not how a non-ideality is shared between them; this is the choice made here.

Requirements, each checked

• Exactly one reactant and one product occupy a surface site, each with a coefficient of one, and both belong to the same site family.

• The site reactant is the species the reaction controls, which KineticReaction takes as the first reactant when the reaction has no crystal: write it first.

• Every participant is a species of cs and carries a standard Gibbs energy.

Used in a KineticsProblem, the reaction makes the whole site family kinetic, the free site included; each of its states then needs a reaction that forms it.

Example

# ≡XOH + Ca+2 → ≡XOHCa+2, the free site first
rxn = Reaction(OrderedDict(free => 1, ca => 1), OrderedDict(bound => 1);
               symbol = "Ca on X", equal_sign = '→')
kr = KineticReaction(cs, rxn, sorption_rate(1.0e-3, cs, rxn))

See also transition_state for minerals.

Laws that run one way

A law restricted to dissolution or to growth, switched off smoothly at saturation. See What a rate law may be.

ChemistryLab.dissolution_only Method

dissolution_only(rate, cs, rxn; width = 1e-3) -> KineticFunc

The rate law rate, of the reaction rxn written as the dissolution of its controlled species, restricted to one direction: dissolution only, while the solution is undersaturated and the phase is present, or growth only, while it is supersaturated.

with w = width, Ω the saturation ratio of rxn (saturation_ratio), n and n₀ the current and initial amounts of the controlled species, and ε = SHRINK_FLOOR. The factor s is zero on the forbidden side and continuous with its slope at Ω = 1, whatever rate is, so neither the rate nor the Jacobian of a stiff integrator jumps there. Above a few widths of ln Ω the law is rate itself, to within (w/ln Ω)².

This is the law of a phase that dissolves and never reprecipitates — an anhydrous clinker phase, or a calcium sulfate whose solution stays supersaturated without it growing back — and, with growth_only, of a phase that grows on its seeds and is never redissolved. The Langmuir factor of an inhibitor is already a mechanism of transition_state, RateModelInhibitor.

Example

law = first_order_rate(k, cs, rxn, ShrinkingCoreArea(BlaineSurfaceArea(blaine)))
kr = KineticReaction(cs, rxn, dissolution_only(law, cs, rxn))

Not symbolic: the factor switches on the sign of ln Ω.

ChemistryLab.growth_only Method

growth_only(rate, cs, rxn; width = 1e-3) -> KineticFunc

dissolution_only in the other direction: rate while the solution is supersaturated in the controlled species of rxn, zero otherwise, switched on across a width width of ln Ω, with no depletion factor.

Seeds, nucleation and growth

A phase that precipitates from its own nuclei: a reservoir of seeds fed by a nucleation rate, and growth on the surface those seeds offer.

ChemistryLab.AbstractNucleationLaw Type

AbstractNucleationLaw

The frequency of nucleation as a function of the saturation ratio, per unit of the extent nuclei form on (nucleation_rate): ThresholdNucleation or ClassicalNucleation.

ChemistryLab.AbstractSeedGeometry Type

AbstractSeedGeometry

The shape of the seeds of a SeededGrowthArea: SphericalSeeds or NeedleSeeds. A geometry gives the volume of one seed and the area a population grows on.

ChemistryLab.ClassicalNucleation Type

ClassicalNucleation(J₂, B) -> ClassicalNucleation

The frequency of classical nucleation theory, written through its value J₂ at Ω = 2,

B the barrier, 16πγ³v²/(3(kT)³) for spherical nuclei. Every derivative of J vanishes as Ω → 1⁺, so the frequency switches on smoothly.

ChemistryLab.FixedExtent Type

FixedExtent(x)

Nuclei form at a frequency per unit of a fixed extent x — an area in m², a volume of solution, whatever J is given per.

ChemistryLab.NeedleSeeds Type

NeedleSeeds(; length, aspect, max_length, width = 1e-3) -> NeedleSeeds

Seeds that are needles: cylinders of length l and radius aspect·l, growing homothetically from the length length up to max_length, where they meet one another. With N = V_res/v₀ seeds, v₀ = π aspect² l₀³, holding the volume V, the current length is l = l₀ (V/V_res)^{1/3}, and the area the growth law multiplies is the rate of change of the volume with the length,

so that a growth law k(Ω) per unit of that area is a lengthening rate dl/dt = V_m k(Ω). c stops the growth at l_max, continuously with its slope, across a fraction width of the length: c(z) = z²/(z + width) above zero and zero below. Lengths as Quantity or Real in m.

ChemistryLab.SeededGrowthArea Type

SeededGrowthArea(reservoir, geometry; T = 298.15, P = 1e5) -> SeededGrowthArea

The area a phase grows on when it grows on its own seeds: those of reservoir, a copy of the phase (with_symbol) that a nucleation_rate feeds, of the shape geometry. Give it as the surface of a transition_state law of the phase, typically restricted to growth (growth_only):

seeds = SeededGrowthArea(cs["Portlandite_nuc"], SphericalSeeds(50u"nm"))
grow = growth_only(transition_state([RateMechanism(k, 1, 1)], cs, rxn, seeds), cs, rxn)

The volume of the population is V = V_m (n_res + n), the molar volume V_m of the reservoir (the phase's, the reservoir being a copy of it) at T and P. The area is that of the geometry (SphericalSeeds, NeedleSeeds), zero without seeds.

ChemistryLab.SphericalSeeds Type

SphericalSeeds(radius) -> SphericalSeeds

Seeds that are spheres of radius radius (a length, Quantity or Real in m), growing as N equal spheres. With N = V_res/v₀ seeds, v₀ = 4πr₀³/3, holding the volume V of the phase and its seeds,

the area of the seeds themselves while nothing has grown on them.

ChemistryLab.SubstrateAmount Type

SubstrateAmount(species)

Nuclei form on another solid in proportion to its amount, J per mole of it.

ChemistryLab.SubstrateSurface Type

SubstrateSurface(species, area)

Nuclei form on the surface of another solid, J per m² of it: the area model area of species (total_area) at its current and initial amounts.

ChemistryLab.ThresholdNucleation Type

ThresholdNucleation(J, Ω_c; width = 1e-3) -> ThresholdNucleation

Nucleation above a critical saturation ratio Ω_c, at a frequency proportional to the excess, J c(Ω/Ω_c − 1) with c(z) = z²/(z + width) above zero and zero below: nothing forms below Ω_c, and the frequency is continuous with its slope there.

ChemistryLab.nucleation_rate Method

nucleation_rate(law, cs, rxn, geometry; on, saturation = nothing,
                T = 298.15, P = 1e5) -> KineticFunc

The rate at which the reservoir of seeds of a phase fills by nucleation: rxn is the dissolution of the reservoir, a copy of the phase, and the rate is

negative, the reaction being written as a dissolution: J from law (ThresholdNucleation, ClassicalNucleation), E the extent nuclei form on (on: a FixedExtent, the amount of another solid, SubstrateAmount, or its surface, SubstrateSurface), v₀ the volume of one seed (seed_volume) and V_m the molar volume of the reservoir at T and P. Ω is the saturation ratio of rxn, the phase's, unless saturation, a function (T, P, lna) -> Ω, says otherwise — the sum of the ratios of the end-members of a solid solution, for instance.

The growth on these seeds is a separate law, of the phase, on a SeededGrowthArea.

ChemistryLab.seed_volume Method

seed_volume(geometry) -> Real

The volume of one seed, m³.

A layer covering a surface

ChemistryLab.CoveringLayer Type

CoveringLayer(layer, substrate, area, Γ; blockable = 1, floor = 1e-6) -> CoveringLayer

A layer, the species named layer, covering the surface of the solid substrate, whose area model is area, at Γ blockable sites per m² (mol/m², a Quantity or a plain Real). The amount of layer counts the sites it blocks, so that the sites of the substrate and the fraction left free are

φ floored at floor (continuously with its slope). At most a fraction blockable of the sites can be blocked, the rest — the defects — staying open. covered multiplies a rate law of the substrate by φ, and covering_rate is the law of the layer itself.

ChemistryLab.covered Method

covered(rate, cs, layer::CoveringLayer) -> KineticFunc

The rate law rate of the substrate of layer, multiplied by the fraction of its sites the layer leaves free, φ (CoveringLayer): a dissolution that only the open sites carry.

ChemistryLab.covering_rate Method

covering_rate(layer::CoveringLayer, cs, rxn; k_form, τ_open, opening = nothing) -> KineticFunc

The rate law of the layer, rxn being the dissolution of the species layer counts, Ω its saturation ratio. While the phase of the layer is stable it forms on the blockable sites still open; once it is not, it opens from the sites already free, a logistic recession:

negative while it forms, the reaction being written as a dissolution; b is blockable, N the sites of the substrate, k_form in s⁻¹, τ_open in s. At a fixed Ω < 1 the blocked fraction θ = n_ℓ/N follows dθ/dt = −θ(1 − θ)(1 − Ω)/τ. The two branches meet at zero at Ω = 1, continuously but not with their slopes. opening, a function (T, P, t, n, lna, n₀) -> g of the run, slows the opening — a species that holds the layer in place, say — and is one when not given.

Kinetics problem

ChemistryLab.EQ_RESIDUAL_TOL Constant

EQ_RESIDUAL_TOL

Largest per-element balance violation, relative to that element's own total, that still counts as a converged speciation. Above it the point is not handed on as a warm start and the solve is retried from a budget-clipped guess.

ChemistryLab.RESTORE_MAXIT Constant

RESTORE_MAXIT

Alternating-projection sweeps allowed when restoring the feasibility of an in-run guess. Exposed because the right value is a trade: the projection converges linearly and a cement can need tens of thousands of sweeps, while this runs at every re-speciation of a run.

ChemistryLab._CONTINUATION_STEPS Constant

_CONTINUATION_STEPS

How many bisection steps a replay may take between the last certified instant and one it cannot certify directly. Each step halves the jump in the component totals, so eight of them reduce it by a factor 256.

ChemistryLab._EQ_GUESS_FLOOR Constant

_EQ_GUESS_FLOOR

Lower floor applied to the starting guess of the equilibrium sub-solve, chosen strictly above the 1e-16 lower bound that EquilibriumProblem imposes.

An interior-point method started on its own bound stalls short of its tolerance and returns MaxIters. Loosening that tolerance is not the fix: measured against Reaktoro on the calcite reference case, the worst species error grows from 4.3 % at tol = 1e-10 to 38 % at 1e-8 and 252 % at 1e-7. Moving the guess inside keeps the tight tolerance and removes the stalls.

ChemistryLab.KineticsProblem Type

struct KineticsProblem{CS, CAL, ES, AM}

Encapsulates a kinetics simulation following Leal et al. (2017).

The ODE state vector u is structured as:

• Without re-speciation: u = [nₖ₁, …, nₖ_K, ξ₁, …, ξ_M, [T | Q]]

• With re-speciation: u = [bₑ₁, …, bₑ_C, nₖ₁, …, nₖ_K, ξ₁, …, ξ_M, [T | Q]]

where bₑ are the element amounts in the equilibrium partition, nₖ the moles of kinetic species and ξ the extents of the kinetic reactions. The trailing slot is present only with a calorimeter: the temperature for SemiAdiabaticCalorimeter, the accumulated heat for IsothermalCalorimeter.

Fields

• system: ChemicalSystem.

• kinetic_reactions: vector of KineticReaction objects.

• initial_state: ChemicalState providing initial moles, T, P.

• tspan: (t_start, t_end) time interval [s].

• calorimeter: nothing, IsothermalCalorimeter, or SemiAdiabaticCalorimeter.

• activity_model: AbstractActivityModel for log-activities.

• equilibrium_solver: solver for re-speciation, or nothing.

• idx_kinetic: indices of kinetic species in system.species.

• idx_equilibrium: indices of equilibrium species.

• ν: stoichiometric matrix (M × N) = SM.N' restricted to kinetic reactions.

• νe, νk: partitions of ν for equilibrium / kinetic species.

• Ae: formula matrix restricted to equilibrium species (C × Nₑ).

Threads

A problem is not meant to be integrated on two threads at once: the runs would share its equilibrium solver, whose optimizer may keep its last answer as the next start. Independent trajectories run in parallel with one problem each.

See also: integrate, KineticsSolver.

ChemistryLab._ReadRecorder Type

_ReadRecorder(data)

A vector that records which of its entries are read, for _rates_read_speciation.

ChemistryLab._abs_residual Method

_abs_residual(Ae, n_e, be) -> Float64

Largest element-balance violation in moles. Reported alongside the relative measure because it is the one a chemist can judge: 1e-10 mol is machine precision whatever the system, and 7e-2 mol is not.

ChemistryLab._budget_clip! Method

_budget_clip!(n_eq, Ae, be)

Clip a starting guess to the element budget: no species may exceed what the totals bₑ can supply, nⱼ ≤ minᵢ bᵢ/Aᵢⱼ over the rows it consumes.

This does not change the feasible set — it only moves the guess into it. It is what unblocks the OPC case: once the sulfate is spent the warm start still carried ettringite at the aluminum budget, three times the sulfur available, and the interior-point solve could not walk back from there.

ChemistryLab._build_kinetics_problem Method

KineticsProblem(cs, kinetic_reactions, initial_state, tspan; ...) -> KineticsProblem

Construct a KineticsProblem from an explicit list of reactions.

Each element of kinetic_reactions must be either a KineticReaction or a Reaction with a :rate entry in its properties. Reaction objects are automatically wrapped via KineticReaction(cs, rxn).

KineticsProblem(cs, initial_state, tspan; ...) -> KineticsProblem

Construct from a ChemicalSystem that has kinetic_species declared (reactions and rates auto-generated via the kinetic_species keyword).

Arguments

• cs: ChemicalSystem.

• kinetic_reactions: AbstractVector of KineticReaction or Reaction objects carrying a :rate property.

• initial_state: ChemicalState providing initial moles, T, P.

• tspan: (t0, tf) time interval. Plain Real → [s]; Quantity → converted.

• calorimeter: nothing (no thermal coupling), IsothermalCalorimeter, or SemiAdiabaticCalorimeter.

• activity_model: activity model for log-activity computation (default: dilute).

• equilibrium_solver: nothing (no re-speciation) or an EquilibriumSolver.

Examples

# From explicit reactions
rxn = Reaction(OrderedDict(sp("C3S") => 1.0, sp("H2O@") => 3.33),
               OrderedDict(sp("Jennite") => 0.167, sp("Portlandite") => 1.5))
rxn[:rate] = parrott_killoh(PK_PARAMS_C3S, "C3S"; α_max)
kp = KineticsProblem(cs, [rxn], state0, (0.0, 7 * 86400.0))

# From kinetic_species in ChemicalSystem
cs = ChemicalSystem(species, primaries;
    kinetic_species = Dict("C3S" => pk_C3S, "C2S" => pk_C2S))
kp = KineticsProblem(cs, state0, (0.0, 7 * 86400.0))

ChemistryLab._cell_point Method

_cell_point(p, u) -> (T, H)

The temperature of the cell and the enthalpy of the paste at the state u of a run under partial equilibrium, on plain numbers, the partition left in p as the warm start of the next instant.

ChemistryLab._cell_residual Method

_cell_residual(p, bv, nkv, n_e, Tv, ΔHv) -> (R, C)

The residual R(T) = H(φ(bₑ, T), nₖ, T) − H₀ + C_v (T − T₀) − ΔH of the energy balance of the cell at the temperature Tv, and its derivative along the equilibrium,

C = R'(T) = C_v + Σᵢ nᵢ c°ₚ,ᵢ(T) + Σₑ ΔₐH⁰ₑ(T) ∂φₑ/∂T ,

the heat capacity of the cell at equilibrium: the partition n_e, solved at Tv, lifted in temperature by the implicit-function theorem at the answer. On plain numbers.

ChemistryLab._cell_temperature Method

_cell_temperature(p, u) -> Union{Nothing, Real}

The temperature of a semi-adiabatic cell under partial equilibrium at the state u, the root of its energy balance (_cell_temperature_value). On dual numbers the root is lifted by the implicit-function theorem: the residual is evaluated on the state's own numbers at the root T*, its partition lifted in bₑ, and

T = T* − (R(u, T*) − R*) / C ,

which carries ∂T/∂u = −(∂R/∂u)/C exactly and keeps the value T* the partition was solved at. At a reaction that holds the cell at its temperature (C infinite, _cell_temperature_bracketed) that derivative is zero: the temperature of a reaction between phases of fixed composition does not depend on how much of them there is. nothing when the root cannot be found.

ChemistryLab._cell_temperature_bracketed Method

_cell_temperature_bracketed(p, bv, nkv, ΔHv, (lo, Rlo), (hi, Rhi))
    -> Union{Nothing, Tuple{Float64, Float64}}

The root of the balance of the cell between lo, where its residual R is Rlo < 0, and hi, where it is Rhi ≥ 0, on the values of R alone: the missing side, if one is (an infinite end), found by steps doubling from _CELL_T_STEP, then the bracket bisected. Returns its lower end, where the partition is that of the assemblage below the root, and the capacity there.

The bracket is bisected to _CELL_T_JUMP first. When R jumps across it, by more than ten times what the capacity at its lower end accounts for (or that capacity cannot be read), the root is a reaction between phases of fixed composition, which takes its heat at one temperature: the cell stays at it while the heat of the paste goes into the reaction, and its temperature does not depend on the state there, which the capacity returned, infinite, says (∂T/∂u = 0, _cell_temperature); the bracket is then bisected on to _CELL_T_STEP, the temperature of the reaction. Judged on that width, the threshold is several joules for a paste, well above the noise of a residual that differences enthalpies of megajoules and well below the heat of such a reaction (7 kJ for the gypsum of 100 g of binder turning to anhydrite, 5.8 kJ for the ettringite and monocarbonate of a CEM I mortar turning to monosulfate). A regular root is bisected on to _CELL_T_TOL. nothing when the partition cannot be solved, or no bracket is found.

ChemistryLab._cell_temperature_value Method

_cell_temperature_value(p, bv, nkv, ΔHv) -> Union{Nothing, Tuple{Float64, Float64}}

The temperature of the cell at the plain state (bv, nkv, ΔHv) and its heat capacity at equilibrium there (_cell_residual): Newton on R, from the last root found (p.cell_root), the nearest state evaluated, or else from the last accepted temperature, until its step falls below _CELL_T_TOL. The temperature returned is the last one the partition was solved at, which the evaluation that follows finds in the cache of _rhs_values.

R grows with T, and the temperatures where it was found negative and positive bracket the root. Where Newton cannot go on, the root is bisected in that bracket, widened first if one side is missing (_cell_temperature_bracketed): when the capacity it reads is not positive, which the partition lifted in temperature gives where its assemblage is about to change, and when its iterations run out, which they do across the jump of R at a reaction between phases of fixed composition (ettringite and monocarbonate giving way to monosulfate in a paste near 53 °C). While the cell is held at such a reaction (the last root found had an infinite capacity), the root is bracketed from it at once. Where Newton converges, nothing else is done. nothing when the partition cannot be solved.

ChemistryLab._equilibrium_subsystem Method

_equilibrium_subsystem(system, idx_equilibrium) -> ChemicalSystem

The chemical system restricted to the equilibrium partition, as the partitioned formulation of Leal et al. (2017) requires.

Its formula matrix is exactly system.CSM.A[:, idx_equilibrium], in the same species order, so the element amounts bₑ carried by the ODE state are handed to it unchanged. The primaries are those of the parent system that survive the restriction — the kinetic minerals never do, not being in the partition.

ChemistryLab._exhausted_coreactant Method

_exhausted_coreactant(p, u; rtol = _FEASIBILITY_RTOL) -> Union{Nothing, Tuple{Int, Float64}}

The first species the extents of the state u take below zero, as its index in the system and its amount, when the species other than the kinetic ones follow the stoichiometry, n = n(0) + νᵀξ, that is without an equilibrium partition; nothing otherwise. Each amount is judged against the most of it the element totals allow, times rtol: the check of a trajectory passes the relative tolerance of its integration, below which the extents that make the amount are not known. A phase consumed in the first seconds slightly ahead of the reaction that forms it, by a few parts in 10⁸ of its calcium, is then not reported.

Only the species whose amount is an amount of substance are judged: a phase, a gas, the solvent. Without a partition a solute stands for the total of its component, whose sign carries the acidity: a clinker phase dissolving into the primaries consumes H+ below zero, which is hydroxide, not a shortfall.

A rate law that reads only its own phase does not see a co-reactant run out: under a Parrott–Killoh rate, C3A + 3 Gp + 26 H2O → ettringite goes on after the gypsum is gone. The right-hand side floors the amount it reads at zero, so the run went on, reported a success, and created the sulfate the extent demanded beyond what the system held.

ChemistryLab._initialize_cell! Method

_initialize_cell!(p, u0)

Set the reference of the energy balance of a calorimeter under partial equilibrium: H₀, the enthalpy of the paste at the partition of the initial element amounts at the initial temperature, from the solve the right-hand side uses, so that the balance holds at the start with ΔH = 0 and the first equilibrium releases nothing. That partition is the warm start of the run.

ChemistryLab._kinetic_partition Method

_kinetic_partition(system, kin_rxns) -> Vector{Int}

The kinetic species of a run: the species each reaction controls and, when one of them occupies a surface site, every member of its site family.

A family's members share one site budget, so they cannot be split between the two partitions. A family stays on the equilibrium side when no reaction controls one of its members, which is the fast surface the partial-equilibrium formulation assumes. When a reaction controls a member, the adsorption itself is the slow step: the free site goes with the complexes, and the site row is then conserved by the stoichiometry of the reactions instead of by the minimization.

A slow family whose budget follows its host (SITES_FOLLOW_HOST) needs a kinetic host: its sites then appear and vanish with an amount the run integrates, and the right-hand side moves the members with them (_kinetic_site_hosts).

Two declarations are refused, each by name:

• a slow family whose budget follows a host at equilibrium: its sites would be created or destroyed by the host's equilibrium while their occupancy is integrated, which nothing here accounts for;

• an electrostatic family sharing its support with a family left at equilibrium: the two would see one potential, computed from a charge the equilibrium solve does not hold.

ChemistryLab._kinetic_site_hosts Method

_kinetic_site_hosts(system, idx_kinetic, idx_equilibrium, Ae) -> Vector{<:_AbstractSiteHost}

The site families of system whose support is SITES_FOLLOW_HOST and whose host is a kinetic species, each with what the right-hand side needs to move its budget: the family, the host's position j in the kinetic partition and its molar mass M; for a fast family, col, the column of Aₑ of its free site; for a slow one, the positions of its members in the kinetic partition and what each holds in the rows of bₑ.

The budget is a function of the host, not a coefficient

At equilibrium a budget that follows its host has to be ν n_host, linear, for the minimization to solve the two together (sites_per_host). A kinetic host is not solved for: its amount is a variable of the ODE, so the budget is known before each solve, and any capacity will do. Over a ShrinkingCoreArea, N(n) = Γ 𝒜₀ (n/n₀)^p: a grain dissolving from the outside loses its sites as the square of its radius, not as its volume. The run carries N in the site row of bₑ,

so that every reader of bₑ — the re-speciation, speciated_states, the certificate — sees the budget of the host's current amount.

Why the free site has to be bare

A site that disappears with its grain is added to bₑ as one free site, with a minus sign, whatever it held. What it held stays in the partition: the calcium of a ≡SOCa⁺ goes back to the solution, which is the point. That conserves matter and charge only if the free site carries nothing but the site — no atom, no charge — since otherwise its atoms would leave with it. A family whose free site is XsOH is refused here; written relative to a bare Xs, the same model reads Xs + H₂O = XsOH⁻ + H⁺ and Xs + Ca²⁺ + H₂O = XsOHCa⁺ + H⁺, with the constants unchanged (the water is at unit activity in the convention the constants are given in).

A slow family

When a rate law controls the occupancy (sorption_rate), the members are kinetic and the budget is their sum, N = Σₖ nₖ. Sites that vanish with the host take every state in proportion to its amount, and what each held returns to the partition; sites the host brings are born free:

with Ṅ = (dN/dn_host)(dn_host/dt) and 𝐜ₖ the column of member k in the rows of bₑ. The members sum to N(n_host) and stay positive; the switch at Ṅ = 0 is continuous, not differentiable.

ChemistryLab._kinetic_state_infeasible Method

_kinetic_state_infeasible(p, u) -> Bool

Whether the kinetic amounts of the state u are outside what the system can hold: one of them negative beyond rounding, or the kinetic species together holding more of an element than the system was given. Each amount is judged against the most of it the element totals allow, so that a trace mineral is held to its own scale. Without an equilibrium partition, also whether the extents take a species they consume below zero (_exhausted_coreactant, to rtol).

ChemistryLab._kinetics_number_type Method

_kinetics_number_type(kp) -> Type

The number type of a run of kp: Float64, or the dual type covering the initial state, the calorimeter and every rate law.

ChemistryLab._lna_params Method

_lna_params(p, T) -> NamedTuple

The parameters the activity model of a run is called with at the temperature T of the cell: p, at that temperature when the cell has one of its own, and with the standard Gibbs energies over RT at that temperature when a phase of the system reads them, as the activities of a phase under the compound energy formalism do (CompoundEnergyModel). Without them that model refuses to compute, since no default would be right.

ChemistryLab._one_speciation Method

_one_speciation(p, guess, be) -> (ok, n_e, abs_res)

One equilibrium solve of φ(bₑ) from guess. Returns ok = false only if the solve threw; a solve that returns a poor composition still returns true, with its element-balance violation in moles, so the caller can compare attempts.

ChemistryLab._paste_enthalpy Method

_paste_enthalpy(p, n_e, nk, T) -> Real

H = Σᵢ nᵢ ΔₐH⁰ᵢ(T) [J] of the composition whose equilibrium partition is n_e, in the order of p.idx_equilibrium, and whose kinetic amounts are nk: what enthalpy gives a state of that composition. Generic in the number types of all three.

ChemistryLab._rates_read_speciation Method

_rates_read_speciation(p) -> Bool

Whether a rate law of the run reads the equilibrium partition: an amount or a log activity of an equilibrium species, or anything the activity model derives from them. Two probes, on the composition the run starts from: the entries of n and ln a each law reads, recorded, and, for a run on plain numbers, the derivative of each law along the amounts of the partition, through the activity model.

A law that reads the partition makes the right-hand side a function of the speciation, which build_kinetics_ode then solves at every evaluation; one that does not, a Parrott–Killoh or Avrami law on its own degree of reaction, leaves it to the step's re-speciation, which is then exact.

ChemistryLab._reconstruction_guess! Method

_reconstruction_guess!(buf, p, be) -> buf

The independent guess: the composition the specimen was cast with, carried onto the current element budget. It holds no active set at all, which is exactly what recommends it where the assemblage is switching.

ChemistryLab._refuse_missing_enthalpy Method

_refuse_missing_enthalpy(system, h_fns)

Refuse a calorimeter under partial equilibrium when a species carries no enthalpy of formation: its term would drop out of −dH/dt, and the heat with it.

ChemistryLab._respeciate_solve! Method

_respeciate_solve!(p, n_eq, be) -> Bool

Solve φ(bₑ) from the guess n_eq, write the result into p.n_full, and record the element-balance residual. Returns false if the solve threw.

Two guesses are available and the order between them is chosen, not fixed. The warm one is the previous speciation; the reconstruction is the cast composition carried onto the current budget. Whichever runs first, the other is tried when the first leaves too much matter unaccounted for, and the better of the two is kept — so the answer does not depend on the order, only the cost does.

WHY THE ORDER IS WORTH CHOOSING. The warm start carries the previous ACTIVE SET, and where the assemblage switches it is the wrong one: an interior-point method started inside a set of phases that no longer exists does not cross over, it exhausts its iterations. Measured on a six-hour hydration, 143 respeciations:

	first solve	then
the 98 calls the warm guess handled	1247 ms	—
the 45 calls where the assemblage switched	8087 ms, rejected	1467 ms from the reconstruction

Eight seconds spent learning that a guess is wrong, then one and a half to get the answer without it. That was 66 % of the whole integration.

An assemblage switch is not a single-step event — a phase takes several steps to exhaust — so the previous call's outcome says which guess to try first, and eq_switching carries it. This is information the problem already has, not a tuning parameter: no threshold is introduced, and the tolerance that decides "too much matter unaccounted for" is that of a converged speciation, EQ_RESIDUAL_TOL of each element row (_unbalanced), whichever guess runs first.

ChemistryLab._restore_feasibility! Method

_restore_feasibility!(n_eq, Ae, be; maxit, tol) -> n_eq

Move a starting guess into {Aₑn = bₑ, n ≥ 0} by alternating projection: clamp to the box, then project onto the affine set through the 7×7 system AₑAₑᵀ.

The Gibbs minimization is posed with A n = b as a hard equality, so a guess that violates it starts the interior-point method outside its own feasible set. On a full OPC this was not a detail: the solve stopped on a point demanding 0.732 mol of sulfate against the 0.267 mol available — 174 % over — and no tolerance, iteration budget, barrier setting or bound changed it, because none of them addresses an infeasible start. A projected-gradient phase-1 proved the set is non-empty (residual 3·10⁻¹²); this is the cheap way to land in it.

Ending on the box clamp rather than the affine step matters: the affine projection alone leaves small negative amounts, which the barrier cannot accept.

maxit deserves attention. Alternating projection converges linearly, at a rate set by the angle between the box and the affine set, and on a cement that angle can be small: at the six-hour instant of an ordinary Portland cement — where the iron row carries 0.013 mol across thirteen species — 200 sweeps left a residual of 8.4e-1, 2000 left 6.7e-2, and 20 000 were needed to reach 6.7e-9.

The default is small on purpose, and measured: it runs at every re-speciation of a run, and on a full OPC the worst in-run balance is 1.1 mol at 200 sweeps against 8.5 at 2000 and 41 at 100 000. A better guess producing a worse answer is a behavior of the back end not yet understood; until it is understood the in-run budget stays where it measures best. Note that the ranking depends on the back-end and should be re-measured if it changes.

A replay (speciated_states) runs a handful of times and buys accuracy instead, asking for a much larger budget.

ChemistryLab._rhs_partition Method

_rhs_partition(p, be, T) -> Union{Nothing, Vector}

The partition at the element amounts be and the temperature T of a state the right-hand side is evaluated at, carrying their dual numbers when they carry some: the Jacobian of a stiff method, a derivative with respect to the parameters of the run. Solved on the values (_rhs_values) and lifted by the implicit-function theorem at the answer, so that the Jacobian of the right-hand side holds ∂r/∂nₑ ⋅ ∂nₑ/∂bₑ exactly.

ChemistryLab._rhs_values Method

_rhs_values(p, bv, Tv) -> Union{Nothing, Vector{Float64}}

The partition at the element amounts bv and the temperature Tv, on plain numbers, by the certified solve warm-started from the last accepted partition, then from the cast composition carried onto bv, and, for an accessor walking a finished run, from the answer of the interior point; nothing when none certifies and the best leaves some element row off by more than a converged speciation does (_unbalanced). The last answer is cached, so that the evaluations of one point by the integrator and by the step's re-speciation solve it once.

ChemistryLab._row_residual Method

_row_residual(Ae, n_e, be) -> Float64

Largest element-balance violation |Aₑnₑ − bₑ| relative to that element's own budget. Unlike a single global scale it cannot hide a small element behind a large one, which is what let a 0.465 mol sulfur violation report as 1.4e-2.

ChemistryLab._with_equilibrium_solver Method

_with_equilibrium_solver(kp::KineticsProblem, es) -> KineticsProblem

Return kp with its equilibrium solver replaced by es, or kp itself when es is nothing. Used by integrate so that a solver passed on the KineticsSolver reaches the ODE; a solver already set on the problem wins when both are given, and the mismatch is reported.

ChemistryLab.build_kinetics_ode Method

build_kinetics_ode(kp::KineticsProblem) -> Function

Build the ODE right-hand-side f!(du, u, p, t) implementing Leal et al. (2017).

State layout:

• u[1:n_be] = bₑ (element amounts in equilibrium partition)

• u[n_be+1 : n_be+n_nk] = nₖ (moles of kinetic species)

• u[n_be+n_nk+1 : n_be+n_nk+M] = ξ (extents of the M kinetic reactions)

• u[end] = with a calorimeter: under partial equilibrium, the change ΔH of the enthalpy of the cell; otherwise T with a SemiAdiabaticCalorimeter, Q with an IsothermalCalorimeter; absent without one

ODE equations (Leal et al., 2017; Eq. 66):

• dnₖ/dt = νₖᵀ r

• dbₑ/dt = Aₑ νₑᵀ r

• dξ/dt = r

• under partial equilibrium, dΔH/dt = −φ(T − T_env) (semi-adiabatic) or 0 (isothermal), the temperature the root of the cell's energy balance (_cell_temperature);

• otherwise dT/dt = (q̇ − φ(ΔT)) / Cp_total (semi-adiabatic) or dQ/dt = q̇ (isothermal), q̇ the heat of the kinetic reactions,

where nₑ = φ(bₑ, T) is the equilibrium partition, solved at the element amounts and the temperature of the state.

ChemistryLab.build_kinetics_params Method

build_kinetics_params(kp::KineticsProblem; ϵ=1e-30) -> NamedTuple

Build the immutable parameter tuple p passed to the ODE function.

Key fields: T, P, ϵ, lna_fn, kin_rxns, species_index, n_initial_full, n_full, cp_fns, rates_buf, index ranges n_be, n_nk, idx_kinetic, idx_equilibrium, νe, νk, Ae.

ChemistryLab.build_u0 Method

build_u0(kp::KineticsProblem; R = _kinetics_number_type(kp)) -> Vector{R}

Build the initial ODE state vector, in the number type R of the run.

Structure of u:

• Without re-speciation: u = [nₖ₁, …, nₖ_K, ξ₁, …, ξ_M]

• With re-speciation: u = [bₑ₁, …, bₑ_C, nₖ₁, …, nₖ_K, ξ₁, …, ξ_M]

• Semi-adiabatic adds T at the end: u = [..., T₀]

The extents of reaction ξ are carried alongside the kinetic moles, integrating dξ/dt = r. They are redundant with nₖ — the two satisfy nₖ = nₖ(0) + νₖᵀ ξ — but they are what makes the non-kinetic amounts available inside the residual, through n = n(0) + νᵀ ξ. Without them, a rate law gating on a species that is not itself kinetic would read a value frozen at t = 0. They also make reaction_extents and state_at exact rather than quadrature-limited.

The calorimeter's slot stays last and is addressed from the end of the vector, so it is unaffected by the presence of ξ.

ChemistryLab.respeciate! Method

respeciate!(p, u) -> Bool

Solve the equilibrium sub-problem once, and write the result into the running composition p.n_full. Returns true when a solve actually happened.

This is the second half of the operator-splitting step: the ODE advances the kinetic minerals with the speciation held frozen, then this function re-equilibrates the equilibrium partition under the element amounts the ODE has just produced.

The element amounts bₑ carried by the state vector are the constraint of that sub-problem (Leal et al., 2017; Eq. 54). solve conserves A·n, so what has to be handed to it is a composition whose element totals are exactly bₑ — here the previous speciation, projected onto bₑ through the pseudo-inverse of Aₑ. Handing over p.n_full unchanged would discard bₑ entirely and leave the element balance to drift.

ChemistryLab.system_enthalpy Method

system_enthalpy(p, u, T) -> Float64

Σᵢ nᵢ ΔₐH⁰ᵢ(T) [J] over the composition of the ODE state u: the kinetic minerals read from u itself, the equilibrium partition from p.n_full.

Why the heat cannot come from the kinetic reactions here

heat_rate sums rᵢ (−ΔᵣH⁰ᵢ) over the KINETIC reactions, which is right when those reactions produce the hydrates. Under partial equilibrium they do not: they dissolve the anhydrous phases into ions, and the hydrates are precipitated by the Gibbs minimization, whose heat that sum cannot see. Measured on an ordinary Portland cement, counting only the dissolution put the semi-adiabatic temperature rise at 207 K where the test gives some tens of kelvin.

The enthalpy of the whole system has no such blind spot. It is a state function, so the heat released between two states at the same temperature is their difference — reactants, ions and hydrates all counted once, with no reaction stoichiometry to write down. This is Lavergne et al. (2018), Eqs. 17–21.

Kinetics solver

ChemistryLab._DEFAULT_KINETICS_SOLVER_FACTORY Constant

_DEFAULT_KINETICS_SOLVER_FACTORY

Ref holding a zero-argument factory () -> KineticsSolver for the default ODE solver used by integrate(kp) (no explicit solver argument).

Priority:

• KineticsOrdinaryDiffEqExt.__init__ registers Rodas5P() on load.

• Only one registration; extensions override it as needed.

Users can override with:

ChemistryLab._DEFAULT_KINETICS_SOLVER_FACTORY[] = () -> KineticsSolver(ode_solver=MyAlg())

ChemistryLab.KineticsSolver Type

struct KineticsSolver{ODE_S, ES}

Bundles the ODE algorithm and optional equilibrium solver used by integrate. Construct once, reuse across multiple KineticsProblem instances.

Fields

• ode_solver: any SciMLBase.AbstractODEAlgorithm (e.g. Rodas5P() from OrdinaryDiffEq), or :auto to let OrdinaryDiffEq choose — its default polyalgorithm detects stiffness at run time and switches. nothing keeps Rodas5P(), which stays the default because on the problems here :auto lands on the same stiff method and returns the same answers: switching would change nothing while removing a reproducible choice. Use nothing before OrdinaryDiffEq is loaded; an error will be raised at solve time.

• equilibrium_solver: optional EquilibriumSolver to re-equilibrate the equilibrium partition once per accepted step of the integrator (a discrete callback), and, when a rate law reads the partition, at every evaluation of the right-hand side (see the speciation keyword of integrate). When nothing, the kinetic minerals evolve without re-speciation (faster, less accurate).

• kwargs: keyword arguments forwarded to DifferentialEquations.solve (e.g. reltol, abstol, saveat, maxiters).

Examples

using OrdinaryDiffEq          # activates KineticsOrdinaryDiffEqExt
using OptimaSolver            # the equilibrium back end and its certified solver

es = EquilibriumSolver(cs, HKFActivityModel(), OptimaOptimizer())
ks = KineticsSolver(; ode_solver=Rodas5P(), equilibrium_solver=es,
                     reltol=1e-8, abstol=1e-10)
sol = integrate(kp, ks)

ChemistryLab.KineticsSolver Method

KineticsSolver(; ode_solver=nothing, equilibrium_solver=nothing, kwargs...) -> KineticsSolver

Construct a KineticsSolver.

kwargs are forwarded to DifferentialEquations.solve (e.g. reltol=1e-8, abstol=1e-10, saveat=0:60:3600).

ChemistryLab.integrate Method

integrate(kp::KineticsProblem; kwargs...) -> ODESolution

Shortcut that uses the default solver registered by _DEFAULT_KINETICS_SOLVER_FACTORY.

Requires OrdinaryDiffEq to be loaded (which sets the default to Rodas5P()).

using OrdinaryDiffEq
sol = integrate(kp)              # uses Rodas5P() by default
sol = integrate(kp; reltol=1e-6) # forward kwargs to the ODE solver

speciation = :auto (the default), :rhs or :frozen says where the right-hand side reads the equilibrium partition: solved at every evaluation, with its derivative in the Jacobian, when a rate law reads it (:rhs, which :auto takes then), or as the last accepted step left it (:frozen, exact for laws that read only the kinetic amounts). A trajectory that reaches kinetic amounts the system cannot hold is returned with retcode = Unstable.

Calorimetry

ChemistryLab.AbstractCalorimeter Type

abstract type AbstractCalorimeter end

Base type for calorimeter models that can be coupled to a kinetics simulation.

Concrete subtypes:

• IsothermalCalorimeter: T = constant, the heat Q(t) released to the bath.

• SemiAdiabaticCalorimeter: variable-T cell (Lavergne et al., 2018).

Under partial equilibrium both balance the enthalpy of the whole composition, the state carrying the change of the enthalpy of the cell; see the theory page Kinetics under partial equilibrium.

ChemistryLab.IsothermalCalorimeter Type

struct IsothermalCalorimeter{T} <: AbstractCalorimeter

Isothermal calorimeter: temperature held constant at T [K]; the heat Q(t) [J] the paste releases to the bath. When the kinetic reactions produce the hydrates, Q = ∫₀ᵗ q̇ dτ with q̇ their heat, integrated as the trailing ODE state. Under partial equilibrium the trailing state is the change of the enthalpy of the cell, zero, and Q = H₀ − H is the enthalpy the whole composition has lost — see cumulative_heat.

Examples

cal = IsothermalCalorimeter(298.15u"K")
kp = KineticsProblem(cs, reactions, state0, tspan; calorimeter = cal)
sol = integrate(kp, ks)
t, Q    = cumulative_heat(sol, cal)
t, qdot = heat_flow(sol, cal)

ChemistryLab.IsothermalCalorimeter Method

IsothermalCalorimeter(T) -> IsothermalCalorimeter

Plain Real → assumed SI [K]; Quantity → converted to K.

ChemistryLab.SemiAdiabaticCalorimeter Type

struct SemiAdiabaticCalorimeter{C, T, F} <: AbstractCalorimeter

Semi-adiabatic calorimeter following the energy balance of Lavergne et al. (2018). When the kinetic reactions produce the hydrates, the temperature is integrated,

where:

• q̇(t) [W] is the instantaneous heat-generation rate,

• φ(ΔT) [W] is the heat-loss function (e.g. linear L·ΔT or quadratic a·ΔT + b·ΔT²),

• Cp [J/K] is the fixed calorimeter heat capacity,

• Σᵢ nᵢ Cp°ᵢ(T) is the temperature- and mole-dependent sample heat capacity (computed from p.cp_fns at every ODE step when available).

Under partial equilibrium the same balance is written on the enthalpy of the cell: the state carries its change ΔH, dΔH/dt = −φ(T − T_env), and the temperature is the root of H(φ(bₑ, T), nₖ, T) − H₀ + Cp (T − T₀) = ΔH, the enthalpy H of the paste taken over its whole composition at the partition solved at that temperature (see the theory page Kinetics under partial equilibrium). temperature_profile returns it.

Fields

• Cp: heat capacity of everything the ODE does not compute for itself — the vessel, the flask, an inert filler — in J/K, stored as a Quantity.

• heat_loss: callable φ(ΔT::Real) -> Real [W].

• T_env: ambient temperature [K] (stored as Quantity).

• T0: initial temperature [K] (stored as Quantity).

Cp must NOT include the sample

The denominator of the balance above is Cp + Σᵢ nᵢ Cp°ᵢ(T), and the second term is recomputed from the database at every step. Adding the sample's own heat capacity to Cp therefore counts it twice and understates the temperature rise — by a factor of about 1.75 on a cement paste at w/c = 0.4, where the paste contributes some 2.5 kJ/K against 0.9 kJ/K for the flask. Three scripts shipped with this package made exactly that mistake, and the field description above used to invite it by saying "calorimeter + sample". Pass the vessel, and let the solver add the paste.

Examples

# Linear heat loss — Newton cooling
cal = SemiAdiabaticCalorimeter(; Cp=4000.0u"J/K", T_env=293.15u"K", L=0.5u"W/K", T0=293.15u"K")

# Quadratic heat loss, the form of Lavergne et al. (2018), Eq. (23), with
# illustrative coefficients. `Cp` is the vessel alone; the paste's own heat
# capacity is added by the solver from the database.
cal = SemiAdiabaticCalorimeter(;
    Cp        = 900.0u"J/K",
    T_env     = 293.15u"K",
    heat_loss = ΔT -> 0.3*ΔT + 0.003*ΔT^2,
    T0        = 293.15u"K",
)

kp = KineticsProblem(cs, reactions, state0, tspan; calorimeter = cal)
sol = integrate(kp, ks)
t, T_vec = temperature_profile(sol, cal)
t, qdot  = heat_flow(sol, cal)

References

• Lavergne et al. (2018).

ChemistryLab.SemiAdiabaticCalorimeter Method

SemiAdiabaticCalorimeter(; Cp, T_env, T0, heat_loss=nothing, L=nothing)

Keyword constructor for SemiAdiabaticCalorimeter.

Exactly one of heat_loss or L must be provided:

• heat_loss: callable ΔT -> [W] (e.g. quadratic ΔT -> a*ΔT + b*ΔT^2).

• L: linear Newton cooling coefficient [W/K]. Sets heat_loss = ΔT -> L * ΔT.

All scalar fields accept plain Real (assumed SI) or Quantity:

• Cp → J/K; T_env → K; T0 → K; L → W/K.

ChemistryLab.cumulative_heat Method

cumulative_heat(sol, cal::IsothermalCalorimeter) -> (t, Q)

Cumulative heat Q(t) [J] the paste has released to the bath, at the instants the solution saved.

In the stoichiometric formulation Q = ∫₀ᵗ q̇ dτ is the trailing ODE state, q̇ being heat_rate, Σᵢ rᵢ(−ΔᵣH⁰ᵢ) over the kinetic reactions, which produce the hydrates. Under partial equilibrium those reactions only dissolve the anhydrous phases into ions and the hydrates are precipitated by the Gibbs minimization, so the heat is the enthalpy the paste has lost, Q = H₀ − H, H = Σᵢ nᵢ ΔₐH⁰ᵢ(T) over the whole composition, its partition solved at the state of each instant and H₀ that of the first equilibrium. Every species then needs an enthalpy of formation, and a system where one lacks it is refused. heat_release computes the same difference from certified speciations.

ChemistryLab.cumulative_heat Method

cumulative_heat(sol, cal::SemiAdiabaticCalorimeter) -> (t, Q)

The heat the paste has released [J] at the instants the solution saved. Under partial equilibrium it is the enthalpy it has lost, which went to warm the vessel or out through its walls, Q = C_v (T − T₀) − ΔH, exactly, ΔH the change of the enthalpy of the cell the state carries. In the stoichiometric formulation, the integral of heat_flow by the rectangle rule.

Under partial equilibrium, as heat_flow and temperature_profile, it is computed on the values of the run: on a run carrying dual numbers it carries no derivative, and a warning says so once.

ChemistryLab.extend_ode! Method

extend_ode!(du, u, p, n_kin, cal::SemiAdiabaticCalorimeter)

Append dT/dt = (q̇ − φ(ΔT)) / Cp_total(T, n) to the ODE right-hand side.

Cp_total = Cp + Σᵢ nᵢ Cp°ᵢ(T) is recomputed at every ODE step from p.cp_fns and p.n_full (Lavergne et al., 2018).

ChemistryLab.heat_flow Method

heat_flow(sol, cal::IsothermalCalorimeter) -> (t, qdot)

Instantaneous heat-generation rate q̇(t) [W] at the instants the solution saved. In the stoichiometric formulation it is the time derivative of the accumulated heat the solution carries, read from the solver's own interpolant (sol(t, Val{1})). Under partial equilibrium it is −dH/dt of the paste at the temperature of the bath, the partition lifted by the implicit-function theorem in the direction the run moves its element amounts.

ChemistryLab.heat_flow Method

heat_flow(sol, cal::SemiAdiabaticCalorimeter) -> (t, qdot)

The heat the paste releases, q̇ = −dH/dt [W], at the instants the solution saved: what warms the vessel and what leaves through its walls, q̇ = C_v dT/dt + φ(T − T_env), C_v the heat capacity of the vessel (cal.Cp). dT/dt is the time derivative of the solver's interpolant in the stoichiometric formulation, and under partial equilibrium the derivative of the root of the cell's energy balance along the rate of the state.

ChemistryLab.heat_rate Method

heat_rate(kinetic_reactions, rates, T_K) -> Real

Compute the instantaneous heat generation rate [W = J/s]:

q̇ = Σᵢ rᵢ(t) × ΔHᵣ,ᵢ(T)

where rᵢ [mol/s] is the net rate of the i-th kinetic reaction (positive = dissolution/forward) and ΔHᵣ,ᵢ(T) [J/mol] is the enthalpy of reaction.

AD-compatible: the ΔₐH⁰ callables accept ForwardDiff.Dual T inputs.

ChemistryLab.temperature_profile Method

temperature_profile(sol, cal::SemiAdiabaticCalorimeter; times = sol.t) -> (t, T)

The temperature of the cell [K] at times, by default the instants the solution saved, the state at another instant read from the solver's interpolant: the state's last entry in the stoichiometric formulation, and under partial equilibrium the root of the cell's energy balance at that state, each instant warm-started from the one before.

Post-processing a solution

ChemistryLab._cell_temperature_rate Method

_cell_temperature_rate(p, u, du) -> Float64

dT/dt of a semi-adiabatic cell under partial equilibrium at the state u moving at du: the derivative of the root of its energy balance along du, from the implicit-function theorem (_cell_temperature).

ChemistryLab._heat_rate_of_states Method

_heat_rate_of_states(sol, kp, states, times) -> Vector

q̇ = −dH/dt [W] at each of states, the composition of the run at times, H = Σᵢ nᵢ ΔₐH⁰ᵢ(T) differentiated along the trajectory: the kinetic amounts, the extents and, in a stoichiometric semi-adiabatic cell, the temperature move at the rates of the solver's interpolant (sol(t, Val{1})). Under partial equilibrium the partition is lifted from the certified state by the implicit-function theorem in the direction the run moves its element amounts and, in a semi-adiabatic cell, its temperature, whose rate is that of the root of the cell's balance (_cell_temperature): an exact derivative, at the certified partition rather than the in-run one.

ChemistryLab._paste_heat_rate Method

_paste_heat_rate(p, n_e, u, du, T) -> Float64

q̇ = −dH/dt [W] of the paste at the plain state u moving at du, its partition n_e at that state and the temperature T: one derivative along the move, the partition lifted from n_e by the implicit-function theorem in the direction of the element amounts and, in a semi-adiabatic cell under partial equilibrium, of the temperature, at the rate of the root of the cell's balance.

ChemistryLab._replay_continuation Method

_replay_continuation(sol, kp, p, des, sub, certified, n_eq, eq, t_prev, t, Tt, Pv, be)
    -> (proved, n_eq, eq, certified)

The continuation of speciated_states from the composition certified, proved at t_prev, to the instant t, whose budget is be and temperature Tt: a homotopy in the component totals, walked forward by an adaptive step. Returns whether t was proved, its composition and state (n_eq, eq as given when it was not), and the last composition certified on the way.

ChemistryLab._replay_full_search Method

_replay_full_search(proved, n_eq, eq, certified, sub, start, T, P, model, b)
    -> (proved, n_eq, eq, certified)

The last resort of speciated_states for an instant neither its starts nor the continuation proved: the full certified search of equilibrate_certified, its start from the linear program, its ideal starts and its restarts, on the instant's budget b, from start in the partition system sub. The first instant has no certified neighbor to walk from, and on a slag paste (Gruyaert et al., 2010), at 2 and 7 days, neither start certified it: the interior-point composition the replay then fell back to had a pH of 15.3, where the full search certifies 12.82 at once from the linear program.

Silent and never raising, as the fallback of equilibrate_certified runs it: the replay says which instants stay unproved. Returns its arguments unchanged when the instant is already proved or the search does not certify it.

ChemistryLab.degrees_of_hydration Method

degrees_of_hydration(sol, kp::KineticsProblem; times = sol.t) -> OrderedDict{String, Vector{Float64}}

Degree of reaction    of every kinetic species, keyed by species symbol.

Species with a zero initial amount are omitted rather than returned as NaN: a phase that was never present has no degree of reaction.

Examples

α = degrees_of_hydration(sol, kp)
α["C3S"][end]                       # final degree of hydration of alite

See also: mean_degree_of_hydration, state_at.

ChemistryLab.extent_residual Method

extent_residual(sol, kp::KineticsProblem, ξ; times = sol.t) -> Float64

Largest absolute discrepancy, in mol, between the kinetic-species moles carried by the integrator and those implied by the extents ξ, that is   .

The two are redundant by construction — nₖ integrates νₖᵀ r while ξ integrates r — so this measures how far the integrator has let them drift, and should sit at the solver's tolerance. A large value means the solve is too loose, not that the post-processing is inaccurate.

ChemistryLab.heat_release Method

heat_release(sol, kp; times = sol.t, reference = nothing, states = nothing)
    -> (times, Q, q̇)

Cumulative heat Q [J] and heat rate q̇ [W] released along the trajectory, from the certified speciations at times.

Q(t) = H(t₀) − H(t) with H = Σᵢ nᵢ ΔₐH⁰ᵢ(T), which is Eq. (17)–(21) of Lavergne et al. (2018): enthalpy is a state function, so its drop between two states at the same temperature is the heat given off, with reactants, ions and hydrates each counted once and no reaction stoichiometry to write down. Pass reference to measure from a state other than the first. Pass states when the certified replay at those instants is already in hand — it is the expensive part, and a caller that needs the compositions anyway should not pay for it twice.

q̇ is computed on the values of the run: on a run carrying dual numbers it carries no derivative, and a warning says so once.

Why this does not read the running composition

The obvious implementation — accumulate the enthalpy the integrator already has in hand at each accepted step — is wrong here, and not subtly. Under partial equilibrium the assemblage at each step comes from an in-run Gibbs minimization that is warm-started and not certified; the enthalpy is worth hundreds of kilojoules per hydrate, so an assemblage off by one phase moves the curve by more than the entire heat of hydration. Measured on an ordinary Portland cement, the in-run figures ran 12.7, 145, 1174, 936 and 631 J/g at 1 h, 6 h, 12 h, 1 d and 2 d — heat that increases, then decreases, which no calorimeter has ever seen.

speciated_states re-solves each instant and checks it against the KKT conditions, and that is what this reads.

Why the heat rate cannot come from the kinetic reactions either

heat_rate sums rᵢ(−ΔᵣH⁰ᵢ) over the kinetic reactions, which is right when those reactions produce the hydrates. Under partial equilibrium they only dissolve the anhydrous phases into ions; the hydrates are precipitated by the minimization, whose heat that sum cannot see. On the same cement it put a semi-adiabatic temperature rise at 207 K. The calorimeters therefore take −dH/dt under partial equilibrium, from the in-run partition and its sensitivity to the element amounts (see cumulative_heat); this function reads certified states and is the reference to check them against.

ChemistryLab.mean_degree_of_hydration Method

mean_degree_of_hydration(sol, kp::KineticsProblem; times = sol.t, weights = :mass,
                         species = nothing) -> Vector{<:Real}

Degree of reaction of the binder as a whole, averaged over the kinetic species.

weights selects the averaging:

• :mass (default) — weighted by the initial mass of each species, the convention used when a single ᾱ is quoted for a cement.

• :mole — weighted by initial moles.

species, a collection of symbols, restricts the average to those kinetic species (the clinker phases of a blend without its slag, for instance); a symbol that is not a kinetic species of kp is refused.

Only species with a non-zero initial amount take part, consistently with degrees_of_hydration.

ChemistryLab.reaction_extents Method

reaction_extents(sol, kp::KineticsProblem; times = sol.t) -> Matrix{Float64}

Extent of reaction    for every kinetic reaction of kp, along the solution sol.

Returns a length(times) × M matrix, M being the number of kinetic reactions, in mol.

The extents are carried by the ODE state — the integrator advances dξ/dt = r alongside the kinetic moles — so this is a read of the solution, not a quadrature. It is exact to the solver's own tolerance at any instant, and costs nothing.

Arguments

• sol: the ODESolution returned by integrate.

• kp: the KineticsProblem that produced it.

• times: instants at which the extents are wanted. Any instants within tspan, in any order; the dense output is interpolated.

See also: state_at, degrees_of_hydration, extent_residual.

ChemistryLab.speciated_states Method

speciated_states(sol, kp::KineticsProblem; times = sol.t) -> Vector{ChemicalState}

The speciated compositions along a solution: the kinetic species read from the ODE state, and the equilibrium partition recovered by re-solving φ(bₑ) at each instant with the element totals the run carried.

This is what state_at deliberately does not do. That function returns the purely kinetic reconstruction n(0) + νᵀξ, because the redistribution performed by the equilibrium solve is not recoverable from the stoichiometry — for a cement that reconstruction is meaningless, putting every dissolved element in solution with not one hydrate. Nor is the composition left in the solver's own buffers usable: that buffer is rewritten at every right-hand-side evaluation, Jacobian differences and rejected steps included, so it is not the accepted composition at any particular time.

Without an equilibrium solver on kp there is nothing to replay, and this returns state_at for each instant.

The instants must be ascending

Each solve is warm-started from the previous one, and the guess is first capped at the element budget and projected back into {Aₑn = bₑ, n ≥ 0}. Both matter, and neither is optional:

• the warm start is the equilibrium of the previous bₑ, so once an element has been spent — the sulfate of an ordinary Portland cement, after the gypsum is gone — it demands more of it than now exists and the interior-point solve begins outside its own feasible set;

• solved instead from a cold guess, φ(bₑ) does not converge at all for a cement: the H⁺ component of bₑ reaches −14 mol, and a solve started from pure water returns no hydrates and a pore solution at pH 6 while the run itself computed 2.2 mol of C-S-H.

The chain is also walked up to the first instant requested, through a few earlier times of the trajectory whose compositions are discarded. That instant is the one with no predecessor of its own, and without the run-up it starts from the cast composition, which carries no active set at all: on a reference OPC its interior-point answer held 56 interior species where the answer has 25 — every candidate hydrate present, four of them at 1e-5 to 1e-6 mol — and neither the certifying solve nor its continuation recovers from that, both inheriting the start. With the run-up that instant is certified, at stationarity 9.1e-13 and element balance 6.1e-15 mol.

Examples

states = speciated_states(sol, kp; times = [1, 7, 28] .* 86400.0)
pH(states[end])                      # pore solution at 28 days
volume_fractions(states[end], groups; reference = state0)

See also: state_at, reaction_extents, volume_fractions.

ChemistryLab.state_at Method

state_at(sol, kp::KineticsProblem, t; ξ = nothing) -> ChemicalState

Reconstruct the full ChemicalState of the system at time t from a kinetics solution.

Every species — kinetic and non-kinetic alike — is obtained from the same stoichiometric replay   , so the returned state satisfies element conservation exactly. Temperature is taken from the solution when a SemiAdiabaticCalorimeter is in use, and from the problem otherwise.

The extents are read from the ODE state, so this is exact to the solver's tolerance. Pass a precomputed extent matrix as ξ — any matrix whose last row holds ξ(t) — to avoid re-interpolating the solution when sweeping many instants.

Re-speciation is not replayed

When the run used an equilibrium solver, the aqueous partition was re-speciated at every accepted step and that redistribution is not recoverable from the stoichiometry alone. The returned state is then the purely kinetic reconstruction; call equilibrate on it to recover the speciated composition.

Examples

ξ = reaction_extents(sol, kp)                      # once
states = [state_at(sol, kp, t; ξ = ξ[i:i, :]) for (i, t) in enumerate(sol.t)]

See also: reaction_extents, degrees_of_hydration, volume_fractions.

The second law along a run

The affinity of each kinetic reaction and the power it dissipates, on the certified compositions of a run. See What a rate law may be.

ChemistryLab.dissipation Method

dissipation(sol, kp::KineticsProblem; times = sol.t, states = nothing,
            lnΩ_tol = 1e-6, rate_rtol = 1e-6) -> NamedTuple

How each kinetic reaction of a run stands with the second law: at each instant, its affinity 𝒜_j (reaction_affinities), its rate r_j evaluated by its own law on the certified composition, and the product 𝒜_j r_j, the power it dissipates, which must not be negative.

The fields are times, affinity (J/mol), rate (mol/s), power (𝒜 r, W), entropy_production (Σ_j 𝒜_j r_j / T, W/K, one value per instant) and violations, one named tuple (time, reaction, affinity, rate) per instant and reaction where 𝒜_j r_j < 0 while |ln Ω_j| > lnΩ_tol and |r_j| > rate_rtol · maxₜ |r_j|.

What it tells

A law written as k F (1 − Ω) with F ≥ 0, such as transition_state or sorption_rate, cannot appear in violations: its rate has the sign of the affinity by construction. A law that does not read the composition, a Parrott–Killoh or Waller degree of reaction, can: it goes on dissolving a phase the solution has become supersaturated with, and the run then dissipates a negative power. So can a law whose reverse rate is not tied to the equilibrium constant of the database: it stops where Ω ≠ 1, and between Ω = 1 and its own stopping point it runs against its affinity.

lnΩ_tol keeps out of violations the instants where the reaction sits at equilibrium to the accuracy of the speciation, whose sign of 1 − Ω is then noise; rate_rtol keeps out those where the law has stopped, its rate a rounding of zero, at an Ω that is not one. Stopping there is allowed: a rate may vanish while the affinity does not, which is a metastable state. Such a stop, at an Ω other than one, is visible in affinity at the last instants.

ChemistryLab.reaction_affinities Method

reaction_affinities(sol, kp::KineticsProblem; times = sol.t, states = nothing)
    -> Matrix{Float64}

The affinity of each kinetic reaction of kp (columns) at each instant of times (rows), in J/mol:

with ν_ij the coefficients of the reaction as the problem holds it, the species it controls at −1: a positive affinity drives the reaction the way a positive rate runs it.

The composition is the certified one, speciated_states, passed as states when already computed; the activities are those the rate laws of the run read, at the temperature of each state. A reaction with a participant that carries no standard Gibbs energy, a glass given a formula for its rate law, has no affinity: its column is NaN.

See also dissipation.

The implicit kinetic step

One fully implicit problem per step, with the reaction extents as unknowns of the same Gibbs minimization. See Which route: two ways to advance in time.

ChemistryLab.KineticStepSolver Type

KineticStepSolver(system, model, reactions; opts...)

A solver for one fully implicit kinetic step, following Leal et al. (2017).

reactions are KineticReaction objects. Every species not carried by a declared reaction is at equilibrium; every declared reaction is kinetic and must have a rate law, which KineticReaction already enforces.

Requires OptimaSolver.

Example

kss = KineticStepSolver(cs, DiluteSolutionModel(), [kr_calcite])
st1 = kinetic_step(kss, st0, 60.0u"s")

See also kinetic_step, integrate for the ODEProblem route.

ChemistryLab._certificate_error Method

_certificate_error(c) -> Float64

One number ranking how far a step is from satisfying its own conditions: stationarity, element balance, the worst supersaturation among absent phases, and the residual of the rate equation. Inf when no certificate was produced.

Used to choose between the pinned and free formulations of a step. A ranking is needed rather than the boolean optimal, because two answers may straddle the tolerance and then the choice is decided by rounding — the same code choosing differently on another machine. Ranked, the gap is fifteen orders of magnitude.

ChemistryLab._kinetic_step_eliminated Method

_kinetic_step_eliminated(kss, state, Δt; ...)

The :species step: the pinned species are computed from the extents rather than solved for, and the extents come from a Newton on nr unknowns around a plain equilibrium over the species that remain.

F(Δξ) = Δξ − Δt·r(n(Δξ)), where n(Δξ) puts the imposed amounts where the stoichiometry says and lets the minimization place the rest under the budget those amounts leave. Backward Euler, as kinetic_step is, since r is read at the composition the step produces.

ChemistryLab._normalize_certificate Method

_normalize_certificate(c) -> NamedTuple

One shape for the certificate whichever route produced it.

kinetic_step fills its certificate reference from OptimaSolver's kkt_certificate on the :reactions route and from optimality_certificate on the :species one, and those name the element-balance field differently — feasibility against balance. Returning two shapes from one keyword is a trap for the caller, so both are mapped onto the names optimality_certificate uses.

ChemistryLab._step_rate_closure Method

_step_rate_closure(kss, des, T, P, t, n0) -> (x, params) -> Vector

The reaction rates at composition x, in mol/s, one per declared reaction.

Built once per step so the StateViews and the initial amounts are not rebuilt at every residual evaluation — the residual and its Jacobian are evaluated once per Newton iteration, on dual numbers.

ChemistryLab._warm_x0 Method

_warm_x0(model, state, n0, ϵ) -> Vector{Float64}

The equilibrated starting guess for an implicit step, or n0 when it cannot be computed.

Why the guess is equilibrated at all

The dual solver admits a mixing phase by a tangent-plane test inside its active-set loop, and from a composition where the phase is entirely absent that admission fails. Measured on C₃S dissolving into a C-S-H solid solution: a cold start left one end-member at 2.7e-9 with a stationarity residual of 6.5, while the same problem started from an equilibrated guess certified at 3.6e-13 with all four end-members present. The failure belongs to the cold start, not to the kinetics — a plain equilibrium fails the same way.

autostart = false, because the continuation fallback of equilibrate_certified is for a caller who has no starting point. Here state is one — the previous instant of the integration — and a handful of extra solves would be paid at every implicit step. Same principle as Reaktoro's coupled use, where the solver is reused across instants.

Why a failure here is not a failure of the step

What this computes is a STARTING POINT, not a result, so strict convergence is suspended for it — the same distinction the continuation makes. Left set, the strict flag would turn an uncertified guess into a raise, the catch would swallow it, and the step would fall back on n0 anyway: a caller asking for strict results would silently get a worse start than a caller who did not.

And whatever else goes wrong, the fallback is the cold start. A guess that cannot be computed must not take the step down with it: the step is judged on its own result, not on the quality of where it began. Both suspensions are scoped to this task, so a trajectory integrated beside another cannot take its flags away — nor leave them behind.

ChemistryLab.kinetic_step Method

kinetic_step(kss, state, Δt; t = 0.0, ϵ = 1e-16, parameters = nothing)
    -> ChemicalState

One implicit kinetic step of duration Δt from state.

t is the absolute time handed to rate laws that depend on it. Pass a Ref as parameters to receive the reaction extents Δξ the step found.

warm_start equilibrates the starting GUESS when the system carries solid solutions, leaving the component totals untouched. Keep it on: a mixing phase is admitted by a tangent-plane test, and from a composition where the phase is absent the step does not reach its answer. Measured on C₃S dissolving into a C-S-H solution, a cold start did not certify: with OptimaSolver 0.7.4 it left one end-member at 2.7e-9 with a stationarity residual of 6.5, and with the exact outer Jacobian since it admits every end-member and stops at a KKT error of 1.75.

The rate laws are evaluated at the end-of-step composition, which is what makes the step implicit and stable on a stiff system. Δt is the caller's choice; kinetic_step_adaptive chooses it from a local error estimate instead.

One step far beyond the relaxation time can land on the other root

The step is a nonlinear problem, and for a rate law that vanishes at equilibrium it has a second solution: the composition with the mineral wholly dissolved. That one satisfies the element balance and the reactivity row while violating Δξ − Δt·M·r(n) = 0 by the entire extent, so the certificate refuses it — pass certificate and check optimal before trusting a step much larger than 1/k.

Which root the Newton finds depends on its path, not on the chemistry. Measured on calcite under r = k(1 − Ω) with k = 10⁻⁵ mol/s: with the outer Jacobian of OptimaSolver 0.7.4, formed by differences, steps of 10⁴ s and 10⁶ s converged to the right root on Julia 1.12 and to the other one on 1.13, and so did 10⁵ s on 1.13.1. With the exact outer Jacobian all three reach the equilibrium root and certify, on 1.13.1. Steps well inside the relaxation time certify on either, and the certificate stays the arbiter.

kinetic_step_adaptive is the guarantee whatever the path: it refuses an uncertified step and halves until one certifies. With 0.7.4 it came back in seven steps on the same case; now it takes the interval in one.

ChemistryLab.kinetic_step_adaptive Method

kinetic_step_adaptive(kss, state, Δt; reltol, abstol, ...) -> (state, Δt_used, Δt_next)

One kinetic step with its length chosen from a local error estimate.

Takes one step of Δt and two of Δt/2, compares the extents, and accepts the finer pair when the estimated error is within tolerance — halving and retrying otherwise. Returns the accepted state, the step actually taken, and a suggestion for the next one.

The estimate is Richardson's for a first-order method: ξ_fine − ξ_coarse is the error of the coarse step to leading order. The error controlled is therefore on the reaction extents, which is what a kinetic step advances; the composition follows from them through an exact equilibrium solve, so it carries no separate error of its own.

Δt is the step to attempt, max_halvings bounds the retries, and the tolerance is abstol + reltol·sⱼ per reaction, where sⱼ is the amount the extent acts on — Σᵢ |Kᵢⱼ| nᵢ, the kinetic minerals of reaction j at the start of the step.

Scaling by the extent instead, reltol·|Δξⱼ|, does not work, and the reason is worth stating because it is the obvious thing to write. Δξ ∝ Δt, so that tolerance vanishes with the step while the equilibrium solve's own noise does not: the measured error then behaves as noise / (reltol·Δt) and grows as the step shrinks. Measured, the controller halved to the floor without advancing and the final answer got worse as the tolerance was tightened — 2.5e-5 at 1e-3 against 9.2e-5 at 1e-5. An ODE integrator's reltol multiplies the solution, not the increment, for exactly this reason.

Example

st, dt_used, dt_next = kinetic_step_adaptive(kss, st0, 3600.0u"s"; reltol = 1e-4)

See also kinetic_step for a step of a length you choose, and integrate for the ODEProblem route, which brings SciML's own step-size control and its choice of stiff solvers.