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_surface: indices of species in AS_SURFACE, i.e. bound to a site.

  • site_families: Nothing when no surface is declared, or a concrete Vector{<:SiteFamily}. Populated through the site_families keyword.

  • site_groups: for each family, the indices of its members in species, the free site first — the order the site mixing and the solver's reference member both rely on.

  • 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._resolve_site_families Method
julia
_resolve_site_families(site_families, species, idx_surface) -> (families, groups)

Resolve declared SiteFamily objects against the species list.

Returns (nothing, Vector{Int}[]) when none is declared, so a system without a surface is byte-identical to what it was before surfaces existed.

Everything refused here is refused rather than discovered later, because the alternative is a wrong number rather than an error:

  • a member that is not in the species list — the family would share a budget with a species the system cannot see;

  • an AS_SURFACE species belonging to no family — it would carry a site pseudo-element into the conservation matrix, and so a row, with nothing mixing on it. This is the surface counterpart of the check that already refuses two solid solutions sharing a composition;

  • two families sharing one pseudo-element — one symbol, one budget, one family, or the site balance silently merges them;

  • a member whose identity does not survive the lookup: the system's species of that name must have the family's formula, must carry the family's site symbol, and must be AS_SURFACE;

  • a species claimed by two families.

Why the last two are not redundant with the third

They were once believed to be. The reasoning was that two families can only share a member if they share a site symbol, that a species carries exactly one, and that a shared symbol is refused above — so no input could reach the case.

The reasoning skips a step. Members are matched to the system by symbol, and a symbol is a label: the species the lookup lands on need not be the species the family validated. Two families with different site symbols and the same member labels therefore resolved to the same indices, one family's conservation row disappeared, and nothing said so. The same gap let a family hold AS_SURFACE copies while the system kept the caller's unqualified originals, leaving idx_surface empty with the site mixing still running.

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.site_families Method
julia
site_families(cs::ChemicalSystem) -> Union{Nothing, Vector{<:SiteFamily}}

The surface site families declared on cs, or nothing when it has no surface.

Mirrors solid_solutions, and for the same reason: a family is a named group of species with its own mixing, and the system has to carry the declaration because the species alone do not say which budget they share.

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
ChemistryLab.surface Method
julia
surface(cs::ChemicalSystem) -> SubArray

A view of the species bound to a surface site, i.e. those in AS_SURFACE.

They are deliberately not in crystal(cs): idx_crystal is read by the solver and by the start repair as "a pure mineral phase", which a site occupancy is not. They are nonetheless counted in the solid compartment of a ChemicalState, because that is where their matter is.

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._solid_indices Method
julia
_solid_indices(system) -> Vector{Int}

The species counted in the solid compartment: the crystalline ones, and the surface complexes.

Written once rather than three times because it is a decision, not a detail. A species bound to a site is matter, and it sits on a mineral, so its amount, its mass and its volume belong with the solid. It is kept out of idx_crystal itself because that index means "a pure mineral phase" to the solver and to the start repair, which a site occupancy is not.

A surface complex has no standard molar volume, so it contributes nothing to the solid volume — the volume it occupies is the host's, already counted.

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.

Loss on ignition ​

What a solid assemblage would lose on heating, from the formulas its phases carry — the total a thermogram integrates to. See The water budget of a hydrating paste for what this gives and what a full thermogravimetric operator would still need.

ChemistryLab._ignition_molar_mass Method
julia
_ignition_molar_mass(system, symbols..., fallback) -> Quantity

The molar mass of a released gas, taken from the system's own species when it carries one under any of symbols, and from fallback otherwise.

Preferring the system's own is the rule this package keeps everywhere: a molar mass written here would drift from the database the rest of the calculation uses. The fallback exists because an assemblage may release a species it never declared — a paste with no gas phase still loses water on ignition.

ChemistryLab.bound_water Method
julia
bound_water(state::ChemicalState) -> Quantity

The water a solid assemblage would release on ignition, in kilograms, from the hydrogen written into its phases.

Every hydrogen in a solid leaves as water — Ca(OH)₂ → CaO + H₂O, CaSO₄·2H₂O → CaSO₄ + 2H₂O — so the amount is H/2 per formula unit, whatever form the hydrogen is written in. That is why this counts hydrogen and not formula water: Ca(OH)₂ has no H₂O in its formula and loses one on ignition, and a rule that looked for H₂O would report zero.

The aqueous phase is not included. Pore solution is water, and evaporable water is not what a hydrate bound — the distinction the water budget is about. Only species in the solid buckets are counted.

See also: ignition_loss, bound_water_per_phase.

ChemistryLab.bound_water_per_phase Method
julia
bound_water_per_phase(state::ChemicalState) -> Vector{Pair{String,Quantity}}

The ignition water each solid species carries, largest first.

What a thermogram would resolve if the decomposition windows were known — so this is the quantity a TGA observation operator needs per phase, and the place those windows would attach.

ChemistryLab.ignition_loss Method
julia
ignition_loss(state::ChemicalState) -> NamedTuple

What a solid assemblage would lose on ignition, as (; water, carbon_dioxide, total) in kilograms.

Hydrogen leaves as water and carbon as carbon dioxide, which is the convention a loss-on-ignition measurement reports and which is exact for the carbonates and hydrates a cement or a clay assemblage is made of.

Carbon is assumed to be carbonate

C in a solid is counted as CO₂. For a cement, a clay or a soil mineral assemblage that is what it is. For organic carbon it is not — the mass is right and the species is not, and a sample with organic matter needs the split made before this is read.

What this is and is not

It is the total a thermogram integrates to, computed from the formulas the database gives. It is not a thermogram: which phase releases what between which temperatures is a set of decomposition windows, and those are literature values rather than consequences of a formula. This package does not carry them, and inventing them would be fabricating the part of the measurement that does the identifying.

Example

julia
loss = ignition_loss(state)
ustrip(us"g", loss.water), ustrip(us"g", loss.total)

Thermogravimetry ​

Turning that total into a curve: a decomposition window per phase, carrying its own provenance, and the operator that makes them identifiable from a measured thermogram rather than only supplied. See A thermogram, and the windows it takes to have one.

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

When a phase releases what it releases, as a logistic step in temperature.

f is the fraction already released at temperature T, so the phase's contribution to a thermogram is m_i f(T) and to its derivative m_i f'(T). A logistic rather than a step because a decomposition is not instantaneous, and rather than something with more shape parameters because two — a midpoint and a width — are already at the edge of what one peak in a thermogram determines.

Fields

  • phase: the species symbol the window belongs to, as the system spells it.

  • midpoint: T₁/₂ in kelvin, where half of the release has happened.

  • width: w in kelvin. The release runs from roughly T₁/₂ − 3w to T₁/₂ + 3w, so a peak 100 K wide has w ≈ 17 K.

  • releases: :water or :carbon_dioxide.

  • fraction: how much of that phase's release this window accounts for, 1 by default. A phase can go in stages — gypsum loses its two waters in two steps, CaSO₄·2H₂O → CaSO₄·½H₂O → CaSO₄ — and one window per stage with fractions summing to one is how that is written. The sum is checked, in thermogram, against each phase and product it is given for.

Both temperature parameters are Traced, which is the point rather than a decoration: a window read from a paper and a window guessed to get a picture on the screen are the same two numbers and are not the same claim.

See also: thermogram, window_parameters.

ChemistryLab.DecompositionWindow Method
julia
DecompositionWindow(phase, midpoint, width; releases = :water, fraction = 1,
                    kind = PROV_UNSTATED, source = "") -> DecompositionWindow

Build a window. midpoint and width are in kelvin, as plain numbers or as Traced values that keep their own provenance.

Passing bare numbers with no kind leaves them PROV_UNSTATED, which is the weakest claim there is — deliberately, so a window nobody sourced never strengthens a result.

ChemistryLab._check_fractions Method
julia
_check_fractions(windows)

Refuse a set of windows whose fractions do not sum to one for some phase and product.

Two windows on one phase with the default fraction would release its mass twice, and the only symptom is a curve that integrates to more than ignition_loss — a silent doubling rather than an error. The check is here rather than at construction because a window does not know what it will be used with.

ChemistryLab._co2_per_phase Method
julia
_co2_per_phase(state) -> Dict{String, Quantity}

The carbon dioxide each solid species would release, by symbol.

ChemistryLab.phases_without_windows Method
julia
phases_without_windows(state, windows) -> Vector{Pair{String,Symbol}}

What state would release on heating and has no window saying when, as phase => product pairs.

A thermogram computed with one of these present integrates to less than ignition_loss says, and the difference is silent. This is how to see it before reading a curve.

Why the product is in the answer

A phase can need two windows. A carbonated hydrate carries hydrogen and carbon, releases water and carbon dioxide, and does so at different temperatures — a hemicarboaluminate is not an exotic case in a cement. Coverage counted per phase rather than per phase => product reports such a phase as covered when only half of it is, and the missing half is precisely the silent shortfall this exists to catch.

ChemistryLab.released_fraction Method
julia
released_fraction(w::DecompositionWindow, T) -> Real

The fraction of w's phase already released at temperature T in kelvin.

ChemistryLab.released_rate Method
julia
released_rate(w::DecompositionWindow, T) -> Real

df/dT — the shape one peak of a DTG curve has, per kelvin.

ChemistryLab.thermogram Method
julia
thermogram(state, windows; temperatures) -> NamedTuple

A thermogravimetric curve for state under windows, as (; temperature, mass, loss, dtg, by_phase):

  • temperature: the grid, in kelvin, as given.

  • mass: the sample mass remaining, in kilograms — the solid mass of state, residue included, because that is what sits on the pan.

  • loss: what has been released, in kilograms, counted from nothing rather than from the first point of the grid — so mass = mass₀ - loss with mass₀ the solid mass, and loss[1] is the tail already gone before the grid begins rather than zero by definition.

  • dtg: -dm/dT in kilograms per kelvin, which is the curve a thermogravimetric analysis actually resolves peaks in.

  • by_phase: the same loss, per phase, so a peak can be attributed.

A window whose phase state has nothing to release from contributes nothing and raises nothing — windows_without_phases is how a typo in a phase name is caught rather than read off a missing peak.

Every phase of state that carries hydrogen or carbon and has no window contributes its mass to the starting point and never leaves — which is the honest behavior, because a phase nobody said when to decompose has not been modeled. phases_without_windows is how to find out that this happened rather than to discover it from a curve that integrates to the wrong total.

A logistic has infinite tails

The curve does not start at exactly zero: a window centered at 400 K with a width of 15 K is already 1.3 ‰ through at 300 K, and that mass is counted as released before the grid begins. It is the model rather than an error — no renormalization happens here, because silently rescaling a curve to make it start at zero would put the discrepancy somewhere a reader cannot see.

Start the grid at least 6w below the lowest midpoint and the leak is under 0.3 %; 10w puts it under 5 × 10⁻⁵. thermogram(...).loss[1] is what it actually is.

What this is

A deconvolution model, not a kinetic one: it says what fraction of a phase has gone at a temperature, not how fast it goes at a heating rate. That is the form published TGA interpretations take, and it is what makes the windows identifiable from one curve — a rate model would need several heating rates before its parameters meant anything.

ChemistryLab.window_parameters Method
julia
window_parameters(windows) -> (θ, names)

The windows' parameters as one vector, [T₁/₂, w] per window, with names — what an optimizer and identifiability take.

The two are returned together because a parameter vector whose entries are not named is a parameter vector nobody can report.

ChemistryLab.windows_without_phases Method
julia
windows_without_phases(state, windows) -> Vector{String}

Which windows name a phase state has nothing to release from.

The mirror of phases_without_windows, and the one that catches a typo: a window on "Portlandit" contributes nothing, raises nothing, and leaves a curve that is simply missing a peak. Both directions are silent by themselves, which is why both are reported rather than one.

A window on a phase that is genuinely absent — a paste with no calcite — is not an error and shows up here too; the two are indistinguishable from inside, and naming them is all this can honestly do.

ChemistryLab.with_window_parameters Method
julia
with_window_parameters(windows, θ; kind = PROV_FITTED, source = "") -> Vector

windows with their parameters replaced by θ, in the order window_parameters produced.

The replacements are marked PROV_FITTED by default, because that is what a number coming out of an optimizer is. Pass kind explicitly when it is not — stepping a parameter to differentiate against it, for instance, does not make the result a fit.