ChemistryLab.jl

Appearance

Chemical systems and states ​

ChemicalSystem ​

ChemistryLab.ChemicalSystem Type
julia
struct ChemicalSystem{T<:AbstractSpecies, R<:AbstractReaction, C, S, SS} <: AbstractVector{T}

An immutable, fully typed collection of chemical species and reactions with derived index structures and stoichiometric matrices.

Immutability guarantees that all derived fields (dict_species, dict_reactions, index vectors, CSM, SM) remain consistent with species and reactions throughout the lifetime of the object. To modify the system, use merge to construct a new ChemicalSystem.

Fields

  • species: ordered list of all species.

  • dict_species: fast O(1) lookup by species symbol.

  • idx_aqueous, idx_crystal, idx_gas: indices by aggregate state.

  • idx_solutes, idx_solvent, idx_components, idx_gasfluid: indices by class.

  • reactions: ordered list of all reactions.

  • dict_reactions: fast O(1) lookup by reaction symbol.

  • CSM: canonical stoichiometric matrix.

  • SM: stoichiometric matrix with respect to primaries.

  • solid_solutions: Nothing when no solid solutions are present, or a concrete Vector{<:AbstractSolidSolutionPhase} describing each solid-solution phase and its end-members. Populated via the solid_solutions keyword constructor.

  • ss_groups: for each solid solution, the indices of its end-members in species.

  • idx_ssendmembers: union of all end-member indices (flattened ss_groups).

  • idx_kinetic: indices of kinetic species (empty when none declared).

ChemistryLab.ChemicalSystem Method
julia
ChemicalSystem(species, primaries=species; kinetic_species, solid_solutions) -> ChemicalSystem

Construct a fully typed ChemicalSystem from a vector of species, an optional vector of primary species, optional kinetic species with rates, and optional solid-solution phases.

All derived fields are computed once at construction time and remain consistent for the lifetime of the object.

Arguments

  • species: vector of AbstractSpecies.

  • primaries: subset used as independent components (default: all species).

  • kinetic_species: nothing (default) or a dictionary / vector of pairs mapping each kinetic species (by name String or Species object) to its rate function. Rate functions must be callable as (T, P, t, n, lna, n_initial) → Real [mol/s] (see KineticFunc). The rate is given per mole of kinetic species (stoichiometric coefficient = 1); the constructor corrects by 1/|νₖ| automatically. When provided, the nullspace N of the stoichiometric matrix is diagonalized so that each kinetic species appears in exactly one reaction. Those reactions are stored in the reactions field with their rate attached via rxn[:rate].

  • solid_solutions: vector of SolidSolutionPhase (default: nothing). When provided, end-members must already appear in species (matched by symbol) and must carry aggregate_state = AS_CRYSTAL and class = SC_SSENDMEMBER.

Examples

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

julia> cs = ChemicalSystem(sp);

julia> length(cs)
2

julia> cs["H2O"] == sp[1]
true
julia
julia> em1 = Species("AFm1"; aggregate_state=AS_CRYSTAL, class=SC_SSENDMEMBER);

julia> em2 = Species("AFm2"; aggregate_state=AS_CRYSTAL, class=SC_SSENDMEMBER);

julia> ss = SolidSolutionPhase("AFm", [em1, em2]);

julia> cs = ChemicalSystem([em1, em2]; solid_solutions=[ss]);

julia> cs.ss_groups
1-element Vector{Vector{Int64}}:
 [1, 2]

julia> cs.idx_ssendmembers
2-element Vector{Int64}:
 1
 2
ChemistryLab.ChemicalSystem Method
julia
ChemicalSystem(species, primaries::AbstractVector{<:AbstractString}; kinetic_species, solid_solutions) -> ChemicalSystem

Convenience constructor that resolves primary species from their symbol strings.

The components must span the species: every species is written as a combination of them, and that combination is the conservation law the equilibrium enforces for it. A list that cannot express a species is refused by name — see StoichMatrix.

Examples

julia
julia> sp = [
           Species("H2O";  aggregate_state=AS_AQUEOUS),
           Species("NaCl"; aggregate_state=AS_CRYSTAL),
       ];

julia> cs = ChemicalSystem(sp, ["H2O", "NaCl"]);

julia> symbol.(cs.SM.primaries)
2-element Vector{String}:
 "H2O"
 "NaCl"

Water alone would not do. Sodium and chlorine would have nowhere to be conserved, and ChemicalSystem(sp, ["H2O"]) raises an ArgumentError naming NaCl rather than projecting it onto H2O — a projection that would balance arithmetically and let the solver make salt out of water.

Base.getindex Method
julia
Base.getindex(cs::ChemicalSystem, i::AbstractString) -> AbstractSpecies

Return the species whose symbol matches i. Runs in O(1) via dict_species.

Examples

julia
julia> cs = ChemicalSystem([Species("H2O"; aggregate_state=AS_AQUEOUS)]);

julia> cs["H2O"] == Species("H2O"; aggregate_state=AS_AQUEOUS)
true
Base.getindex Method
julia
Base.getindex(cs::ChemicalSystem, i::Int) -> AbstractSpecies

Return the species at position i.

Examples

julia
julia> cs = ChemicalSystem([Species("H2O"; aggregate_state=AS_AQUEOUS)]);

julia> cs[1] == Species("H2O"; aggregate_state=AS_AQUEOUS)
true
Base.merge Method
julia
Base.merge(cs1::ChemicalSystem, cs2::ChemicalSystem) -> ChemicalSystem

Construct a new ChemicalSystem from the union of two systems.

Species and reactions are unioned by symbol — duplicates from cs2 are discarded. CSM and SM are built from scratch from the full species list. Primaries are taken as the union of both systems' primaries, filtered to those actually present in the merged species list.

In case of symbol conflict (species or reactions), cs1 takes priority over cs2. The return type is inferred from the merged collections and may differ from typeof(cs1) or typeof(cs2) if they contain different concrete types.

Examples

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

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

julia> cs = merge(cs1, cs2);

julia> length(cs)
3
Base.merge Method
julia
Base.merge(css::ChemicalSystem...) -> ChemicalSystem

Construct a new ChemicalSystem from the union of an arbitrary number of systems, processed left-to-right. Earlier systems take priority over later ones in case of symbol conflicts.

Examples

julia
julia> cs1 = ChemicalSystem([Species("H2O"; aggregate_state=AS_AQUEOUS, class=SC_AQSOLVENT)]);

julia> cs2 = ChemicalSystem([Species("H+"; aggregate_state=AS_AQUEOUS, class=SC_AQSOLUTE)]);

julia> cs3 = ChemicalSystem([Species("OH-"; aggregate_state=AS_AQUEOUS, class=SC_AQSOLUTE)]);

julia> cs = merge(cs1, cs2, cs3);

julia> length(cs)
3
Base.size Method
julia
Base.size(cs::ChemicalSystem) -> Tuple

Return the size of the underlying species vector.

Examples

julia
julia> cs = ChemicalSystem([Species("H2O"; aggregate_state=AS_AQUEOUS)]);

julia> size(cs)
(1,)
ChemistryLab._declared Method
julia
_declared(ss) -> String

The name of the declaration a solid-solution phase is an instance of.

Falls back to the phase's own name for any AbstractSolidSolutionPhase that does not carry the field, so a user-defined phase type keeps working.

ChemistryLab._expand_instances Method
julia
_expand_instances(species, solid_solutions) -> (species, solid_solutions)

Give every declaration asking for instances > 1 the extra species it needs.

A SolidSolutionPhase with instances = k becomes k phases: the declaration itself, then k-1 copies named "$name#2", "$name#3", … whose end-members are copies of the originals under "$symbol#2", "$symbol#3", … Each copy shares the whole property dictionary of the species it was made from, so the two are one substance under two labels and carry byte-identical thermodynamic data.

Why the copies are needed at all: the composition vector has one entry per species, so a species belongs to exactly one phase. Two coexisting compositions of one substance — which is what the Gibbs minimum is inside a spinodal — cannot be written down without the substance appearing twice.

Why it is done by duplicating the species rather than by letting two phases share them: it keeps ss_groups disjoint, so the activity assembly, the mole fraction fill and the optimality certificate need no change at all. Sharing an index would make the last group written win.

The duplicated column leaves the row rank of the conservation matrix unchanged — it is a copy of a column already there — so conservation is untouched, and the copies start at zero amount, so a budget computed as A n is unchanged too.

Returns the inputs unchanged, and without copying, when no declaration asks for more than one instance. That is every system the shipped data describes.

ChemistryLab._instances Method
julia
_instances(ss) -> Int

How many coexisting compositions a declaration asks for; 1 for a phase type that does not carry the field.

ChemistryLab._refuse_overlapping_solid_solutions Method
julia
_refuse_overlapping_solid_solutions(solid_solutions)

Refuse a system in which two declared solid solutions describe the same substance, naming the pair.

The case this exists for is the calcium silicate hydrate. CEMDATA18 carries three descriptions of it — CSHQ, CNASH_ss and the ECSH family — and they are three models of one gel, not three phases. Declaring two of them counts the same hydrate twice: the calcium, the silicon and the alkalis all enter the element balance once and come out distributed over two phases that are supposed to be alternatives.

The overlap is exact rather than approximate, which is what makes it detectable here: KSiOH (an end-member of CSHQ), ECSH1-KSH and ECSH2-KSH all carry the formula ((KOH)2.5SiO2H2O)0.2. Two end-members of two different declared phases with the same composition are therefore the signature, and the check is composition-based rather than name-based so that it does not depend on the database's naming.

Two end-members of the SAME phase may of course share nothing — that is a mixture — and a pure phase repeating a mixing phase's composition is a separate question the rank test upstream already refuses.

Instances of one declaration are exempt, and that exemption is the whole reason SolidSolutionPhase carries a declared field. A miscibility gap is represented by the same binary present twice, on purpose, as two coexisting compositions — which is composition overlap of exactly the kind this function refuses. Telling the deliberate case from the mistake cannot be done by composition, since they look identical; it is done by provenance, and two phases share a declared name only when ChemicalSystem itself made the second from the first.

ChemistryLab.aqueous Method
julia
aqueous(cs::ChemicalSystem) -> SubArray

Return a view of all aqueous species.

Examples

julia
julia> cs = ChemicalSystem([
           Species("H2O";  aggregate_state=AS_AQUEOUS),
           Species("NaCl"; aggregate_state=AS_CRYSTAL),
       ]);

julia> length(aqueous(cs))
1

julia> aggregate_state(aqueous(cs)[1]) == AS_AQUEOUS
true
ChemistryLab.components Method
julia
components(cs::ChemicalSystem) -> SubArray

Return a view of all component species.

Examples

julia
julia> cs = ChemicalSystem([Species("SiO2"; aggregate_state=AS_CRYSTAL, class=SC_COMPONENT)]);

julia> class(components(cs)[1]) == SC_COMPONENT
true
ChemistryLab.crystal Method
julia
crystal(cs::ChemicalSystem) -> SubArray

Return a view of all crystalline species.

Examples

julia
julia> cs = ChemicalSystem([
           Species("H2O";  aggregate_state=AS_AQUEOUS),
           Species("NaCl"; aggregate_state=AS_CRYSTAL),
       ]);

julia> aggregate_state(crystal(cs)[1]) == AS_CRYSTAL
true
ChemistryLab.gas Method
julia
gas(cs::ChemicalSystem) -> SubArray

Return a view of all gas-phase species.

Examples

julia
julia> cs = ChemicalSystem([Species("CO2"; aggregate_state=AS_GAS)]);

julia> aggregate_state(gas(cs)[1]) == AS_GAS
true
ChemistryLab.gasfluid Method
julia
gasfluid(cs::ChemicalSystem) -> SubArray

Return a view of all gas/fluid species.

Examples

julia
julia> cs = ChemicalSystem([Species("CO2"; aggregate_state=AS_GAS, class=SC_GASFLUID)]);

julia> class(gasfluid(cs)[1]) == SC_GASFLUID
true
ChemistryLab.get_reaction Method
julia
get_reaction(cs::ChemicalSystem, sym::AbstractString) -> AbstractReaction

Return the reaction identified by symbol sym. Runs in O(1) via dict_reactions.

Examples

julia
cs = ChemicalSystem(
    [Species("H2O"; aggregate_state=AS_AQUEOUS)];
);
get_reaction(cs, "some_rxn")  # returns the Reaction with that symbol
ChemistryLab.kinetic_species Method
julia
kinetic_species(cs::ChemicalSystem) -> SubArray

Return a view of the kinetic species declared at construction time. Empty when no kinetic species were declared.

ChemistryLab.solid_solutions Method
julia
solid_solutions(cs::ChemicalSystem) -> Nothing | Vector{<:AbstractSolidSolutionPhase}

Return the registered solid-solution phases, or nothing if none were declared.

Examples

julia
julia> em1 = Species("Em1"; aggregate_state=AS_CRYSTAL, class=SC_SSENDMEMBER);

julia> em2 = Species("Em2"; aggregate_state=AS_CRYSTAL, class=SC_SSENDMEMBER);

julia> cs = ChemicalSystem(
           [em1, em2];
           solid_solutions=[SolidSolutionPhase("SS", [em1, em2])],
       );

julia> solid_solutions(cs) isa Vector
true

julia> length(solid_solutions(cs))
1
ChemistryLab.solutes Method
julia
solutes(cs::ChemicalSystem) -> SubArray

Return a view of all aqueous solute species.

Examples

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

julia> class(solutes(cs)[1]) == SC_AQSOLUTE
true
ChemistryLab.solvent Method
julia
solvent(cs::ChemicalSystem) -> AbstractSpecies

Return the unique solvent species directly (not a view), since a chemical system contains at most one solvent.

Examples

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

julia> class(solvent(cs)) == SC_AQSOLVENT
true

ChemicalState ​

ChemistryLab.ChemicalState Type
julia
struct ChemicalState{C, S, Q<:AbstractQuantity, R<:Real}

Immutable container holding the thermodynamic state of a ChemicalSystem.

Molar amounts are always stored internally in mol regardless of the input unit. Each species can be provided independently as a molar amount (mol) or as a mass (g, kg, etc.) — the constructor converts each entry individually using the molar mass M stored in the corresponding species.

The struct itself is immutable — fields cannot be reassigned. However, n, T, and P are stored as Vector to allow in-place mutation via set_quantity!, set_temperature!, and set_pressure!.

system is a shared reference: cloning via Base.copy does not duplicate the underlying ChemicalSystem.

Fields

  • system: reference to the underlying ChemicalSystem.

  • n: molar amounts (mol), one per species — mutable in place.

  • T: temperature (K) — 1-element Vector, mutable in place.

  • P: pressure (Pa) — 1-element Vector, mutable in place.

  • n_phases: moles per phase (liquid, solid, gas, total) — PhaseQuantities{Q}.

  • m_phases: mass per phase (liquid, solid, gas, total) — PhaseQuantities{Q}.

  • V_phases: volume per phase (liquid, solid, gas, total) — PhaseQuantities{Q}.

  • pH: pH of the liquid phase, or nothing if H⁺ is absent — R | Nothing.

  • pOH: pOH of the liquid phase, or nothing if OH⁻ is absent — R | Nothing.

  • porosity: (V_liquid + V_gas) / V_total, or NaN if volumes unavailable — R.

  • saturation: V_liquid / (V_liquid + V_gas), or NaN if pore volume is zero — R.

Examples

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

julia> state = ChemicalState(cs; T=298.15u"K", P=1u"bar");

julia> length(state.n)
2

julia> ustrip(state.T[])
298.15
ChemistryLab.ChemicalState Method
julia
ChemicalState(system::ChemicalSystem, values::AbstractVector; T, P) -> ChemicalState

Construct a ChemicalState with explicit initial amounts or masses. Each entry is converted to moles independently — mixed units allowed.

Examples

julia
julia> cs = ChemicalSystem([
           Species("H2O"; aggregate_state=AS_AQUEOUS, class=SC_AQSOLVENT),
           Species("NaCl"; aggregate_state=AS_CRYSTAL),
       ]);

julia> state = ChemicalState(cs, [55.5u"mol", 5.844u"g"]);

julia> ustrip(moles(state, "H2O"))
55.5

julia> isapprox(ustrip(moles(state, "NaCl")), 0.1; rtol=1e-4)
true
ChemistryLab.ChemicalState Method
julia
ChemicalState(system::ChemicalSystem; T, P, n) -> ChemicalState

Construct a ChemicalState from a ChemicalSystem with optional initial temperature, pressure, and molar amounts (default: all zero).

Arguments

  • system: the ChemicalSystem describing the species.

  • T: temperature in K (default: 298.15u"K").

  • P: pressure (default: 1u"bar").

  • n: molar amounts in mol (default: zeros).

Examples

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

julia> state = ChemicalState(cs; T=298.15u"K", P=1u"bar");

julia> ustrip(state.T[])
298.15

julia> ustrip(state.P[]) ≈ 1e5
true

julia> all(iszero.(ustrip.(state.n)))
true
ChemistryLab.PhaseQuantities Type
julia
PhaseQuantities{Q}

Named tuple type alias (liquid::Q, solid::Q, gas::Q, total::Q) for phase-aggregated quantities (moles, mass, or volume).

Each field holds the total of the corresponding thermodynamic phase. Q is typically an AbstractQuantity carrying SI units (mol, kg, or m³).

Base.:* Method
julia
Base.:*(state::ChemicalState, α::Real) -> ChemicalState

Return a new ChemicalState with all molar amounts scaled by α. Temperature, pressure, and the underlying ChemicalSystem are unchanged. The operation is non-mutating — a copy is returned.

Examples

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

julia> state = ChemicalState(cs, [2.0u"mol"]);

julia> s2 = state * 3.0;

julia> ustrip(moles(s2, "H2O"))
6.0

julia> ustrip(moles(state, "H2O"))   # original unchanged
2.0
Base.:* Method
julia
Base.:*(α::Real, state::ChemicalState) -> ChemicalState

Equivalent to state * α.

Base.:+ Method
julia
Base.:+(s1::ChemicalState, s2::ChemicalState) -> ChemicalState

Combine two chemical states by adding their species amounts.

  • Same system (s1.system === s2.system): the result shares the system reference.

  • Different systems: a merged system is created via merge(s1.system, s2.system) (union of species, s1 takes priority for duplicates).

  • T, P: taken from s1. A warning is emitted if s2 has different T or P.

  • Derived quantities (pH, volumes, …) are recomputed from the summed moles.

Examples

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

julia> s1 = ChemicalState(cs, [2.0u"mol"]);

julia> s2 = ChemicalState(cs, [3.0u"mol"]);

julia> ustrip(moles(s1 + s2, "H2O"))
5.0
Base.:/ Method
julia
Base.:/(state::ChemicalState, α::Real) -> ChemicalState

Return a new ChemicalState with all molar amounts divided by α. The operation is non-mutating — a copy is returned.

Examples

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

julia> state = ChemicalState(cs, [4.0u"mol"]);

julia> s = state / 2.0;

julia> ustrip(moles(s, "H2O"))
2.0
Base.copy Method
julia
Base.copy(state::ChemicalState) -> ChemicalState

Create a clone of a ChemicalState that shares the same ChemicalSystem reference but owns independent copies of all mutable fields.

Modifying the clone does not affect the original, and vice versa. The underlying ChemicalSystem is not duplicated.

Examples

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

julia> state = ChemicalState(cs, [55.5u"mol", 0.1u"mol"]);

julia> clone = copy(state);

julia> set_quantity!(clone, "Na+", 0.5u"mol");

julia> ustrip(moles(state, "Na+"))
0.1

julia> ustrip(moles(clone, "Na+"))
0.5

julia> clone.system === state.system
true
Base.show Method
julia
Base.show(io::IO, state::ChemicalState)

Compact single-line representation of a ChemicalState.

Base.show Method
julia
Base.show(io::IO, ::MIME"text/plain", state::ChemicalState)

Detailed multi-line display for ChemicalState. Shows molar amounts, masses, and volumes for each species grouped by phase, with phase totals and scalar diagnostics (pH, pOH, porosity, saturation).

ChemistryLab._auto_seed_neutral_pH! Method
julia
_auto_seed_neutral_pH!(state::ChemicalState)

When the aqueous solvent (H₂O@) is present with non-zero amount and both H⁺ and OH⁻ are in the system at zero/negligible concentration, seed them at the neutral pH concentration   [mol/L] using the T- and P-dependent water autoprotolysis constant.

Does nothing if any of the three species is absent, if H⁺ or OH⁻ already have non-negligible amounts (i.e. the user set them explicitly), or if pKw cannot be computed.

ChemistryLab._compute_V_phases Method
julia
_compute_V_phases(system, n, T, P) -> NamedTuple

Compute volume per phase from n, T, P and standard molar volumes V⁰. Gas phase falls back to ideal gas law if V⁰ is not available for all gas species.

ChemistryLab._compute_m_phases Method
julia
_compute_m_phases(system, n) -> NamedTuple

Compute mass per phase from species vector n and molar masses.

ChemistryLab._compute_n_phases Method
julia
_compute_n_phases(system, n) -> NamedTuple

Compute moles per phase from species vector n.

ChemistryLab._compute_pH Method
julia
_compute_pH(system, n, V_liquid) -> Union{Real, Nothing}

Compute pH using the most reliable species between H⁺ and OH⁻.

  • If c_H⁺ ≥ c_OH⁻: pH = -log10(c_H⁺)

  • If c_OH⁻ > c_H⁺: pOH = -log10(c_OH⁻), then pH = pKw(T) - pOH

pKw(T) is retrieved from the reaction H2O = H+ + OH- in the system's reaction dictionary, evaluated at the current T and P. Returns nothing if neither H⁺ nor OH⁻ is present, or if liquid volume is zero.

ChemistryLab._compute_pKw Method
julia
_compute_pKw(system::ChemicalSystem, T, P) -> Union{Real, Nothing}

Compute pKw = -logK⁰(T, P) for the water dissociation reaction H2O@ = H+ + OH- reconstructed on the fly from the species present in system.

Returns nothing if any of H2O@, H+, or OH- is absent from the system.

ChemistryLab._compute_pOH Method
julia
_compute_pOH(system, n, T, P, V_liquid) -> Union{Real, Nothing}

Compute pOH symmetrically to _compute_pH:

  • If c_OH⁻ ≥ c_H⁺: pOH = -log10(c_OH⁻)

  • If c_H⁺ > c_OH⁻: pH = -log10(c_H⁺), then pOH = pKw(T) - pH

Returns nothing if neither species is present or volume is zero.

ChemistryLab._compute_porosity Method
julia
_compute_porosity(V_phases) -> Real

Compute porosity = (V_liquid + V_gas) / V_total. Returns NaN if total volume is zero.

ChemistryLab._compute_saturation Method
julia
_compute_saturation(V_phases) -> Real

Compute saturation = V_liquid / (V_liquid + V_gas). Returns NaN if pore volume is zero.

ChemistryLab._entry_to_moles Method
julia
_entry_to_moles(v::AbstractQuantity, s::AbstractSpecies) -> AbstractQuantity

Convert a single value v to moles for species s. If v has amount dimension (mol), it is returned as-is. If v has mass dimension, it is divided by the molar mass M of s. Otherwise an error is raised.

ChemistryLab._has_molar_volume Method
julia
_has_molar_volume(s::AbstractSpecies) -> Bool

Return true if species s has a standard molar volume V⁰ available.

ChemistryLab._molar_volume Method
julia
_molar_volume(s::AbstractSpecies) -> SymbolicFunc

Return the standard molar volume SymbolicFunc of species s. Must be called as _molar_volume(s)(T=T, P=P; unit=true) to get a quantity.

ChemistryLab._update_derived! Method
julia
_update_derived!(state::ChemicalState)

Recompute and update in place all derived quantities after any mutation of n, T, or P. Called automatically by set_quantity!, set_temperature!, and set_pressure!.

ChemistryLab.enthalpy Method
julia
enthalpy(state::ChemicalState, s) -> Union{AbstractQuantity, Nothing}

The contribution n × ΔₐH⁰(T, P) of one species, or nothing when it carries no enthalpy of formation. s may be a species or its symbol.

ChemistryLab.enthalpy Method
julia
enthalpy(state::ChemicalState) -> AbstractQuantity

Sum of nᵢ ΔₐH⁰ᵢ(T, P) over the species that carry an enthalpy of formation.

Only differences of this quantity are physical. ΔₐH⁰ is referred to the elements, so the absolute number has no meaning on its own; but two states built on the same element budget share that reference exactly, and it cancels in the difference. enthalpy(state₀) - enthalpy(state) is therefore the heat released in going from one to the other at fixed T and P — which is what a calorimeter measures, and what makes a hydration curve computable from the database alone, with no calibrated heat of reaction anywhere.

Species without ΔₐH⁰ contribute nothing; missing_enthalpy lists them, and a run where that list is non-empty is not a closed heat balance.

Examples

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

julia> state = ChemicalState(cs, [55.5u"mol"]);

julia> enthalpy(state) isa AbstractQuantity
true
ChemistryLab.heat_capacity Method
julia
heat_capacity(state::ChemicalState, s) -> Union{AbstractQuantity, Nothing}

The contribution n × Cp⁰(T, P) of one species, or nothing when it carries no heat capacity. s may be a species or its symbol.

ChemistryLab.heat_capacity Method
julia
heat_capacity(state::ChemicalState) -> AbstractQuantity

Sum of nᵢ Cp⁰ᵢ(T, P) over the species that carry a heat capacity.

This is the standard-state sum: it ignores excess contributions from mixing, which for a cement pore solution are small against the solids. It is what an adiabatic temperature rise ΔT = Q / C needs.

ChemistryLab.mass Method
julia
mass(state::ChemicalState, s::AbstractSpecies) -> AbstractQuantity

Return the mass of species s, computed as n × M.

Examples

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

julia> state = ChemicalState(cs, [55.5u"mol"]);

julia> ustrip(uconvert(us"g", mass(state, cs[1]))) ≈ 55.5 * 18.015
true
ChemistryLab.mass Method
julia
mass(state::ChemicalState, sym::AbstractString) -> AbstractQuantity

Return the mass of the species identified by symbol sym.

Examples

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

julia> state = ChemicalState(cs, [55.5u"mol"]);

julia> ustrip(uconvert(us"g", mass(state, "H2O"))) ≈ 55.5 * 18.015
true
ChemistryLab.mass Method
julia
mass(state::ChemicalState) -> NamedTuple

Return mass per phase (liquid, solid, gas, total).

Examples

julia
julia> cs = ChemicalSystem([
           Species("H2O";  aggregate_state=AS_AQUEOUS, class=SC_AQSOLVENT),
           Species("NaCl"; aggregate_state=AS_CRYSTAL),
       ]);

julia> state = ChemicalState(cs, [55.5u"mol", 0.05u"mol"]);

julia> mass(state).total isa AbstractQuantity
true
ChemistryLab.missing_enthalpy Method
julia
missing_enthalpy(state::ChemicalState) -> Vector{String}

Symbols of the species that carry no ΔₐH⁰, and whose amounts are therefore absent from enthalpy.

A heat balance is only closed when this is empty, and it is worth checking rather than assuming: a single missing hydrate silently removes its entire heat of formation from the curve.

ChemistryLab.moles Method
julia
moles(state::ChemicalState, s::AbstractSpecies) -> AbstractQuantity

Return the molar amount of species s in mol.

Examples

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

julia> state = ChemicalState(cs, [55.5u"mol"]);

julia> ustrip(moles(state, cs[1]))
55.5
ChemistryLab.moles Method
julia
moles(state::ChemicalState, sym::AbstractString) -> AbstractQuantity

Return the molar amount of the species identified by symbol sym.

Examples

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

julia> state = ChemicalState(cs, [55.5u"mol"]);

julia> ustrip(moles(state, "H2O"))
55.5
ChemistryLab.moles Method
julia
moles(state::ChemicalState) -> NamedTuple

Return moles per phase (liquid, solid, gas, total).

Examples

julia
julia> cs = ChemicalSystem([
           Species("H2O";  aggregate_state=AS_AQUEOUS, class=SC_AQSOLVENT),
           Species("NaCl"; aggregate_state=AS_CRYSTAL),
       ]);

julia> state = ChemicalState(cs, [55.5u"mol", 0.05u"mol"]);

julia> ustrip(moles(state).liquid)
55.5
ChemistryLab.pH Method
julia
pH(state::ChemicalState) -> Union{Real, Nothing}

Return the pH of the liquid phase, or nothing if H⁺ is absent.

This is −log₁₀ c(H⁺), a concentration in mol/L over the computed liquid volume, and when the solution is alkaline it is reconstructed from OH⁻ through pKw because OH⁻ is the better-resolved species. No activity coefficient enters.

Not the pH another geochemical code reports

GEM-Selektor, PHREEQC and Reaktoro report −log₁₀ a(H⁺), the activity on the molality scale. That is pH(state, model), a different quantity. On a Portland cement pore solution at I ≈ 0.2 mol/kg, with γ(H⁺) ≈ 0.61, the two differ by about 0.21 units — 13.31 here against 13.10 there. Comparing the wrong one against another code means chasing a convention rather than a result.

See also: pH(state, model), pOH, activity_coefficients.

Examples

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

julia> state = ChemicalState(cs, [55.5u"mol", 1e-7u"mol"]);

julia> pH(state) isa Union{Real, Nothing}
true
ChemistryLab.pOH Method
julia
pOH(state::ChemicalState) -> Union{Real, Nothing}

Return the pOH of the liquid phase, or nothing if OH⁻ is absent.

Like pH, this is a concentration in mol/L, not an activity. For the activity convention use pOH(state, model).

Examples

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

julia> state = ChemicalState(cs, [55.5u"mol", 1e-7u"mol"]);

julia> pOH(state) isa Union{Real, Nothing}
true
ChemistryLab.porosity Method
julia
porosity(state::ChemicalState, reference::ChemicalState) -> NamedTuple

Porosity of a setting binder, in the sealed-curing convention: referred to the volume of reference — the fresh material — held fixed, and accounting for the chemical shrinkage.

Returns (; liquid, void, total), all dimensionless fractions of the reference volume:

  • liquid — the water-filled porosity, V_liquid(state) / V_ref.

  • void — the empty porosity created by the reactions, (V_ref − V_total(state)) / V_ref. This is the Le Chatelier contraction: hydration products occupy less space than the reactants they consume, and in a sealed specimen the deficit is gas, not a reduction in size.

  • total — their sum, the quantity a transport or micromechanical model wants.

Consistent by construction with volume_fractions called with the same reference: total equals the sum of its "water"-like and "void" entries.

Examples

julia
p = porosity(state_28d, state_0)
p.total                     # 0.375 for a w/c = 0.5 paste
p.void                      # 0.072 — the chemical shrinkage
p.liquid / p.total          # degree of saturation

See also: saturation(state, reference), chemical_shrinkage, volume_fractions.

ChemistryLab.porosity Method
julia
porosity(state::ChemicalState) -> Real

Return (V_liquid + V_gas) / V_total, or NaN if the total volume is zero.

This is not the porosity of a setting binder

Both ends of that ratio are wrong for a hydrating cement:

  • the denominator is the current total volume, which shrinks as hydration proceeds, whereas a sealed specimen keeps the volume it was cast with;

  • the numerator has no gas term unless gas species were declared, so the empty porosity left by the chemical shrinkage — products occupying less space than the reactants they consume — is structurally invisible.

The two errors compound. On a w/c = 0.5 paste at 28 days this returns 0.327 where the porosity referred to the specimen is 0.375, the total volume having shrunk by 7.2 %.

Use porosity(state, reference) for a binder. This method remains the right one for a fixed-volume aqueous system, where nothing shrinks.

Examples

julia
julia> cs = ChemicalSystem([
           Species("H2O";  aggregate_state=AS_AQUEOUS, class=SC_AQSOLVENT),
           Species("NaCl"; aggregate_state=AS_CRYSTAL),
       ]);

julia> state = ChemicalState(cs, [55.5u"mol", 0.05u"mol"]);

julia> porosity(state) isa Real
true
ChemistryLab.pressure Method
julia
pressure(state::ChemicalState) -> AbstractQuantity

Return the current pressure.

Examples

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

julia> state = ChemicalState(cs; T=298.15u"K", P=1u"bar");

julia> isapprox(ustrip(pressure(state)), 1e5; rtol=1e-4)
true
ChemistryLab.rescale! Method
julia
rescale!(state::ChemicalState, target::AbstractQuantity) -> ChemicalState

Scale all molar amounts in-place so that the total quantity of the matching physical dimension equals target.

target dimensionQuantity brought to target
molmoles(state).total
kg (mass)mass(state).total
m³ (volume)volume(state).total

All derived quantities (pH, volume, porosity, …) are recomputed after scaling. Returns state for chaining.

Examples

julia
rescale!(state, 1.0u"mol")    # total moles  → 1 mol
rescale!(state, 1.0u"kg")     # total mass   → 1 kg
rescale!(state, 1.0u"m^3")    # total volume → 1 m³
rescale!(state, 500u"g")      # total mass   → 500 g
ChemistryLab.saturation Method
julia
saturation(state::ChemicalState, reference::ChemicalState) -> Real

Degree of saturation of a setting binder: the fraction of its porosity that is water-filled, liquid / total from porosity(state, reference).

Falls to NaN when the porosity is zero. Unlike the one-argument method, the empty porosity left by the chemical shrinkage is counted, so a sealed paste desaturates as it hydrates even though no water ever leaves it.

ChemistryLab.saturation Method
julia
saturation(state::ChemicalState) -> Real

Return the saturation V_liquid / (V_liquid + V_gas), or NaN if pore volume is zero.

Examples

julia
julia> cs = ChemicalSystem([
           Species("H2O";  aggregate_state=AS_AQUEOUS, class=SC_AQSOLVENT),
           Species("NaCl"; aggregate_state=AS_CRYSTAL),
       ]);

julia> state = ChemicalState(cs, [55.5u"mol", 0.05u"mol"]);

julia> saturation(state) isa Real
true
ChemistryLab.set_neutral_pH! Method
julia
set_neutral_pH!(state::ChemicalState) -> ChemicalState

Set H⁺ and OH⁻ concentrations to neutral pH at the current temperature and pressure, using the water autoprotolysis constant :

Requires the system to contain H2O@ (solvent), H+, and OH-. The liquid volume is estimated from the current water amount.

Unlike _auto_seed_neutral_pH! (which only triggers when H⁺/OH⁻ are zero), this function always overwrites the current values — useful inside loops where the state is reused across iterations.

Examples

julia
set_quantity!(s, "H2O@", 1.0u"kg")
set_neutral_pH!(s)   # H⁺ and OH⁻ at neutral, T/P-dependent
ChemistryLab.set_pressure! Method
julia
set_pressure!(state::ChemicalState, P::AbstractQuantity) -> ChemicalState

Set the pressure in place and update all derived quantities.

Examples

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

julia> state = ChemicalState(cs; T=298.15u"K", P=1u"bar");

julia> set_pressure!(state, 2u"bar");

julia> isapprox(ustrip(pressure(state)), 2e5; rtol=1e-4)
true
ChemistryLab.set_quantity! Method
julia
set_quantity!(state::ChemicalState, s::AbstractSpecies, n::AbstractQuantity) -> ChemicalState

Set the molar amount of species s in place and update all derived quantities. If n has mass dimension, it is automatically converted to moles using M.

Examples

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

julia> state = ChemicalState(cs, [55.5u"mol"]);

julia> set_quantity!(state, cs[1], 10.0u"mol");

julia> ustrip(moles(state, "H2O"))
10.0
ChemistryLab.set_quantity! Method
julia
set_quantity!(state::ChemicalState, sym::AbstractString, n::AbstractQuantity) -> ChemicalState

Set the molar amount of the species identified by symbol sym in place and update all derived quantities.

Examples

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

julia> state = ChemicalState(cs, [55.5u"mol"]);

julia> set_quantity!(state, "H2O", 10.0u"mol");

julia> ustrip(moles(state, "H2O"))
10.0
ChemistryLab.set_temperature! Method
julia
set_temperature!(state::ChemicalState, T::AbstractQuantity) -> ChemicalState

Set the temperature in place and update all derived quantities.

Examples

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

julia> state = ChemicalState(cs; T=298.15u"K", P=1u"bar");

julia> set_temperature!(state, 350.0u"K");

julia> ustrip(temperature(state))
350.0
ChemistryLab.temperature Method
julia
temperature(state::ChemicalState) -> AbstractQuantity

Return the current temperature.

Examples

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

julia> state = ChemicalState(cs; T=298.15u"K", P=1u"bar");

julia> ustrip(temperature(state))
298.15
ChemistryLab.volume Method
julia
volume(state::ChemicalState, s::AbstractSpecies) -> Union{AbstractQuantity, Nothing}

Return the volume contribution of species s as n × V⁰(T,P). Returns nothing if V⁰ is not available for s.

Examples

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

julia> state = ChemicalState(cs, [55.5u"mol"]);

julia> volume(state, cs[1]) isa Union{AbstractQuantity, Nothing}
true
ChemistryLab.volume Method
julia
volume(state::ChemicalState, sym::AbstractString) -> Union{AbstractQuantity, Nothing}

Return the volume contribution of the species identified by symbol sym. Returns nothing if V⁰ is not available.

Examples

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

julia> state = ChemicalState(cs, [55.5u"mol"]);

julia> volume(state, "H2O") isa Union{AbstractQuantity, Nothing}
true
ChemistryLab.volume Method
julia
volume(state::ChemicalState) -> NamedTuple

Return volume per phase (liquid, solid, gas, total).

Examples

julia
julia> cs = ChemicalSystem([
           Species("H2O";  aggregate_state=AS_AQUEOUS, class=SC_AQSOLVENT),
           Species("NaCl"; aggregate_state=AS_CRYSTAL),
       ]);

julia> state = ChemicalState(cs, [55.5u"mol", 0.05u"mol"]);

julia> volume(state).total isa AbstractQuantity
true

Volume fractions ​

ChemistryLab.chemical_shrinkage Method
julia
chemical_shrinkage(state::ChemicalState, state₀::ChemicalState) -> Quantity

Volume lost to the chemical reactions between the reference state state₀ and state,   , positive for the usual contraction (Le Chatelier).

Both states must share the same ChemicalSystem. Species without a standard molar volume V⁰ contribute nothing, exactly as in volume — see missing_molar_volumes to find out which those are before trusting the number.

Examples

julia
δV = chemical_shrinkage(state_28d, state_0)
ustrip(δV / volume(state_0).total)          # relative contraction, dimensionless

See also: volume_fractions, missing_molar_volumes.

ChemistryLab.missing_molar_volumes Method
julia
missing_molar_volumes(state::ChemicalState; atol = 0.0) -> Vector{String}

Symbols of the species that are present in state (amount in mol strictly above atol) but carry no standard molar volume V⁰, and are therefore invisible to volume, porosity and volume_fractions.

A non-empty result means the volume balance is incomplete. Databases differ in coverage — CEMDATA18 supplies V⁰ for the cement phases but a species added by hand, or one taken from an aqueous-only dataset, may not have it.

Examples

julia
isempty(missing_molar_volumes(state)) || @warn "incomplete volume balance" missing_molar_volumes(state)
ChemistryLab.volume_fractions Method
julia
volume_fractions(state::ChemicalState, groups; kwargs...)
    -> OrderedDict{String, Float64}

Volume fractions aggregated into named families.

groups is any iterable of name => symbols pairs, where symbols is a species symbol or a collection of them — the granularity a mean-field homogenization scheme consumes, where "C-S-H" is one phase rather than four solid-solution end members.

Keyword arguments are those of volume_fractions(state); when reference is given, the chemical-shrinkage void is appended under void_key as a group of its own.

Species that appear in no group are collected under other_key when their total fraction exceeds atol in magnitude, and dropped otherwise — silently discarding matter would defeat the purpose of a volume balance, and the leftover may be negative when only aqueous solutes are unassigned. A species listed in two groups raises an error rather than being counted twice.

Examples

julia
groups = [
    "anhydrous" => ["C3S", "C2S", "C3A", "C4AF"],
    "C-S-H"     => ["CSHQ-TobH", "CSHQ-TobD", "CSHQ-JenH", "CSHQ-JenD",
                    "KSiOH", "NaSiOH"],
    "CH"        => "Portlandite",
    "AFt"       => "ettringite",
    "water"     => "H2O@",
]
f = volume_fractions(state, groups; reference = state_0)

See also: volume_fractions(state).

ChemistryLab.volume_fractions Method
julia
volume_fractions(state::ChemicalState; reference = nothing, void_key = "void")
    -> OrderedDict{String, Float64}

Volume fraction of every species carrying a standard molar volume, keyed by species symbol. Species with a zero amount are omitted.

An individual fraction may be negative: aqueous solutes have negative standard partial molar volumes, electrostriction contracting the solvent around the ion (V⁰(OH⁻) ≈ -4.7 cm³/mol). Those contributions are kept, so that this function and volume always agree; they are negligible at the concentrations of a pore solution, and grouping folds them back into the liquid.

Two normalizations, selected by reference:

  • reference = nothing (default) — fractions are relative to the current total volume of the state, and sum to 1 by construction. This is the right choice for a closed, fully-saturated system.

  • reference::ChemicalState — fractions are relative to volume(reference).total, held fixed. They then sum to less than 1, and the deficit is returned under void_key as the chemical-shrinkage void. This is the sealed-curing convention of Lavergne et al. (2018): the specimen keeps its volume while the reactions consume some, and the resulting empty porosity is a phase of the microstructure, not a rounding error.

Passing reference = state itself is equivalent to the default.

Examples

julia
f = volume_fractions(state_28d; reference = state_0)
f["Portlandite"], f["void"]
sum(values(f))                       # 1.0

See also: volume_fractions(state, groups), chemical_shrinkage, missing_molar_volumes, porosity.

An oxide analysis as an element budget ​

ChemistryLab.glass_species Method
julia
glass_species(oxides; symbol, M = 100.0u"g/mol", name, ΔₐG⁰) -> Species

A pseudo-species for a material that has no formula unit — a blastfurnace slag, a fly ash, a calcined clay — built from its oxide analysis, so that a rate law can consume it and release every element the analysis reports.

oxide_budget solves this for an equilibrium: a glass enters as a contribution to b and never needs to be a species at all. A kinetic run cannot do that, because a rate law consumes a species. So the glass has to be given a formula, and where that formula comes from decides what the dissolution can produce.

Why not pick a representative mineral

Because the elements it leaves out cannot come back. A slag written as anorthite, CaAl₂Si₂O₈, carries no magnesium — so no hydrotalcite can form from it, which is the one phase a slag is certain to make. Writing the formula from the analysis puts every reported element into the budget in its reported proportion, and the question stops being which mineral the glass resembles.

The two masses, and why they differ

The formula unit carries the elements of the reported oxides. M is the mass of material it stands for. These are not the same number, and the difference is the part of the analysis that is not modeled — loss on ignition, and the minor oxides a datasheet omits. The analysis is not renormalized, for the reason oxide_budget gives: scaling the reported fractions up to one would invent material. modeled_mass_fraction in the returned species' properties records what fraction of the material the formula actually accounts for.

ΔₐG⁰ is optional and defaults to a placeholder far below anything the system contains, so that the dissolution is always favored. That is sound for a rate law of the waller or Parrott–Killoh kind, which never reads the saturation ratio; it is not sound for a mechanistic rate driven by Ω, and such a law needs a real Gibbs energy, which a glass does not have.

Examples

julia
slag = Dict("CaO" => 0.41, "SiO2" => 0.36, "Al2O3" => 0.11,
            "MgO" => 0.08, "SO3" => 0.02)
sp = glass_species(slag; symbol = "GGBS", M = 95.0u"g/mol")
atoms(sp)                        # Ca, Si, Al, Mg, S and O, in the reported ratio
rate = waller(WALLER_PARAMS_SLAG, "GGBS"; α_max = 0.9)

See also: oxide_budget, waller.

ChemistryLab.oxide_budget Method
julia
oxide_budget(oxides, primaries; mass = 100.0u"g") -> Vector{Float64}

The component totals b contributed by a material reported as an oxide analysis, for use as the right-hand side of A n = b.

oxides maps an oxide formula to its mass fraction — the shape of a cement or slag datasheet, "CaO" => 0.41 and so on. mass is how much of the material the budget is for. The analysis is not renormalized: if the fractions do not sum to one, what is missing is loss on ignition and minor oxides the sheet does not report, and silently scaling them up would invent material.

Each oxide's molar mass is computed from its formula by Species, never written down, and each is decomposed over primaries by primary_decomposition, which refuses an oxide the primaries cannot express.

julia
# A ground granulated blast-furnace slag, as a datasheet reports it
slag = Dict("CaO" => 0.41, "SiO2" => 0.36, "Al2O3" => 0.11, "MgO" => 0.08)
b = oxide_budget(slag, cs.SM.primaries; mass = 60.0u"g")

This says what a glass contains, not what it does

An oxide budget is a statement of composition and carries no information about reactivity. A slag and a quartz sand of the same analysis give the same b, and an equilibrium calculation will dissolve both completely. How much of the glass has actually reacted is a kinetic quantity, supplied from outside — by a degree of reaction, or by a rate law such as waller — and multiplying the budget by it is the caller's responsibility, not this function's.

ChemistryLab.primary_decomposition Method
julia
primary_decomposition(species, primaries) -> Vector{Float64}

The coefficients writing species as a linear combination of primaries.

Solved by column-pivoted QR rather than \, because the matrix is singular whenever the primaries are not all exercised, and refused above a residual of 1e-8: a species outside the span of the primaries has no decomposition, and returning a least-squares approximation of one would put elements into the budget that the species does not contain.