Solution Chemistry

struct Composition

A Stock-level analog acting on Species values rather than Reagents. It describes real quantities of species ("5g of Na⁺") without committing to which physical Reagent supplies them. See composition for the (one-directional) derivation from a Stock.

Stores molar amounts natively (not mass) – this is what dissociation math (_reagent_moles, dissolution(::Reagent)) naturally produces, and what total_concentration/pH need, with no dependency on any individual Species's own molecular_weight. mass is therefore the derived accessor, returning missing when a chemical is present but its mass can't be computed (unknown molecular_weight) – as opposed to 0u"g" when the chemical isn't present at all.

_reagent_moles(r::Reagent, quantity)

Convert a reagent's tracked quantity (mass for a Solid, volume for a Liquid) to moles, using molecular_weight(r) – and, for a volume, density(r) first (volume × density → mass). Returns missing if a required physical property isn't known.

atoms(s::Stock) -> Dict{Symbol,Unitful.Amount}
atoms(s::Stock,element::Symbol) -> Unitful.Amount
atoms(c::Composition) -> Dict{Symbol,Unitful.Amount}
atoms(c::Composition,element::Symbol) -> Unitful.Amount

Return the amount (in moles) of each element's atoms in a stock or Composition, found by summing every species' amount times its atom count from its elemental formula. With an element symbol such as :Fe, return just that element's amount (0u"mol" if it is absent). Multiply by Unitful.Na to get a number of atoms.

Only reagents contribute; organisms in a culture are not counted. A reagent contributes nothing, and the result is a lower bound, when it has no registered formula or when its amount cannot be computed (for example a liquid without a density). In that case atoms emits one warning naming the reagents. A reagent declared opaque contributes nothing without a warning.

julia> s = 18.015u"mL" * water + 0.1u"mol" * NaCl; # water defined as "H2O", NaCl as Na⁺ + Cl⁻

julia> atoms(s,:H)
2.0 mol

julia> atoms(s,:Cl)
0.1 mol

composition(s::Stock)

Derive the Composition (species-level, in real molar quantities) of s: for each Reagent in s, convert its tracked quantity to moles (via _reagent_moles) and distribute across its dissolution(::Reagent) Dissolution's species, summing across every reagent present.

There is deliberately no Stock-from-Composition conversion — going the other way is underdetermined (a target chemical profile doesn't uniquely determine which reagent(s) produce it).

mass(r::Composition,c::Species)

Return the mass of c tracked in Composition r: 0u"g" if c isn't present at all, missing if c is present but its molecular_weight is unknown (so a mass genuinely can't be computed), otherwise the mass derived from r's stored molar amount.

molar_amount(r::Composition,c::Species)

Return the molar amount of c tracked in Composition r (0u"mol" if absent) – Composition's native, directly-stored unit.

net_hydrogen_ion_concentration(s::Stock)

Return the net molar concentration of hydrogen ions in s: H⁺ minus OH⁻, derived from composition(::Stock). This nets the one acid/base ion pair explicitly by name rather than relying on signed stoichiometry in registered Dissolutions (which are otherwise always non-negative, real formulas) — see pH, its only consumer.

pH(s::Stock; ionic_strength_correction=true, water_correction=false, water_correction_source=nothing)

Estimate the pH of s. General case: partitions s's composition into any registered AcidBaseSystem families (see acid_base_system) plus every other charged species, and root-finds the generalized electroneutrality condition (see CHESSCore._charge_balance_residual) – optionally corrected for ionic strength via the Davies equation (see activity_coefficient), which self-consistently iterates since γ depends on the very speciation it corrects. Set ionic_strength_correction=false for the ideal (infinite-dilution) estimate instead.

Fast path: when no reagent in s has a registered AcidBaseSystem (the common complete/strong-electrolyte-only case – e.g. a stock built only from NaCl, HCl, NaOH, ...), this computes the exact original closed-form formula directly from net_hydrogen_ion_concentration (7.0 if net-neutral, -log10 of the net H⁺ concentration if acidic, 14 - (-log10(...)) if basic) – guaranteed bit-identical to CHESSCore's original strong-electrolyte-only pH, not merely expected to converge to the same place via the general solver. Strong bases should still register their real dissolution (e.g. Na⁺ + OH⁻ for NaOH) via set_dissolution! – mixing an acid and a base nets out through this same mechanism either way.

If s has no registered AcidBaseSystem and no charged species at all (e.g. a Mixture of pure glucose, or an Empty stock) – as opposed to charged species that happen to exactly cancel (e.g. equimolar HCl + NaOH, or plain NaCl) – the returned 7.0 isn't a computed result, just the neutral default with nothing behind it; this case emits a @warn so it isn't mistaken for a real answer.

Set water_correction=true to add a background open-system correction (e.g. dissolved atmospheric CO2) to s's equilibrium before estimating pH – most stocks don't need this level of detail and should leave it false (the default). Modeled as an OpenSystemSpecies (a fixed free-species concentration, not a fixed total dose) so its contribution grows correctly with pH – important for poorly-buffered basic solutions (e.g. a pure conjugate-base salt), which real atmospheric CO2 shifts much more than a fixed dose would. With no water_correction_source given, this uses default_water_correction (registered by a lab module, e.g. CHESSLabConstants's atmospheric-CO2 correction, via set_default_water_correction!); pass water_correction_source explicitly to use a different OpenSystemSpecies instead.

See also speciation for a per-family breakdown at the solved pH.

total_concentration(s::Stock,chem::Species)

Return the molar concentration of chem in s, derived from composition(::Stock) — e.g. the total potassium ion concentration in a stock containing dissolved KCl.

volume(r::Composition,c::Species,density::Unitful.Density)
volume(r::Composition,c::Species,source::Reagent)

Return the volume of c tracked in Composition r, using an explicit density or a reference Reagent to source it from. Species deliberately has no density field of its own — density is a property of physical state, not of the species (e.g. ice vs. liquid water), so it stays on Reagent. missing if mass can't be computed (unknown molecular_weight).

const Kw

The thermodynamic (activity-basis) autoionization constant of water at 25°C, 1.0e-14. Used by pH's general equilibrium solver; corrected for ionic strength (via activity_coefficient) exactly like any other equilibrium constant when ionic_strength_correction=true.

const acid_base_systems

Registry mapping a specific Reagent value to its AcidBaseSystem – mirrors dissolutions (a per-value, not per-type, registry, for the same reason: individual reagents are themselves named values of a small set of types). Prefer set_acid_base_system! over mutating this directly. Absence (see acid_base_system) means "inert with respect to proton transfer," not an error – the overwhelming majority of reagents (salts, sugars, antibiotics, ...) are exactly that.

const default_water_correction

Holds the OpenSystemSpecies used by pH(::Stock)'s water_correction=true when no water_correction_source is given – nothing until a lab module registers one via set_default_water_correction! (e.g. CHESSLabConstants's atmospheric-CO2 correction). Mirrors acid_base_systems in spirit, but there's exactly one relevant default here rather than a per-reagent registry.

const ion_parameters

Registry mapping a specific Species value to its Truesdell-Jones (extended Debye-Hückel) ion-size parameter å (Ångströms) and empirical b parameter, as (å,b). Prefer set_ion_parameters! over mutating this directly. Absence means "no citable per-species parameter" – _activity_coefficient falls back to activity_coefficient (Davies) for any Species not registered here, exactly like acid_base_systems's absence-means-inert convention.

abstract type AbstractAnalyticalSpecies

Common supertype for the equilibrium solver's per-family inputs – an AnalyticalSpecies (fixed total concentration, the mass-balance case; every family a Stock's own composition produces) or an OpenSystemSpecies (fixed free-species concentration, an open-system boundary condition like atmospheric CO2 – see pH's water_correction keyword). Both are consumed identically by _species_concentrations, which dispatches on the concrete type.

struct AcidBaseSystem

A linear "conjugate family" of Species protonation states, from fully-protonated to fully-deprotonated – e.g. phosphoric acid's [H3PO4, H2PO4⁻, HPO4²⁻, PO4³⁻], or a monoprotic acid's [HA, A⁻]. pKa[i] is the equilibrium constant (as a p-value) linking species[i] ⇌ species[i+1] + H⁺.

This same linear-chain shape represents zwitterions too (e.g. an amino acid's [cation, zwitterion, anion]) – nothing about the type assumes the fully-protonated reference state is neutral; only that each successive state carries exactly one fewer positive charge, which the constructor verifies against each Species's own registered charge.

Registered against a Reagent via set_acid_base_system! – independent of (and additional to) that reagent's Dissolution, which continues to answer the separate question of complete-dissociation mass bookkeeping (molecular_weight/composition). A Reagent may have both, neither, or just one.

struct AnalyticalSpecies

An AcidBaseSystem family together with its total analytical concentration in a particular Stock (the sum, across every protonation state already present in the stock's composition, divided by volume_estimate) – the mass-balance input to the equilibrium solver. Built by _analytical_species.

struct OpenSystemSpecies

An AcidBaseSystem family held at a fixed free-species concentration (system.species[1]) rather than a fixed total – the correct model for a species in equilibrium with an external reservoir at constant activity/partial pressure (e.g. dissolved atmospheric CO2, whose [CO2(aq)] is set by Henry's law regardless of how much HCO3⁻/CO3²⁻ the solution's pH happens to pull it into). Unlike AnalyticalSpecies, the family's total concentration is not fixed – it grows with pH as more of the conjugate base forms, which is exactly the physics a fixed-total dose gets wrong for poorly-buffered basic solutions (see pH's water_correction keyword docs). Registered as the default correction via set_default_water_correction!.

struct SpeciationResult

Per-family diagnostic output of speciation: fractions[i] and concentrations[i] correspond to system.species[i] at the solved pH.

_activity_coefficient(chem::Species, z::Integer, I::Real)

Per-species activity coefficient at ionic strength I (dimensionless mol/L): Truesdell-Jones (log10(γ) = -A z² √I/(1+B å √I) + b I, standard extended-Debye-Hückel constants A=0.5085, B=0.3281 at 25°C) using chem's registered (å,b) from ion_parameters if present, else activity_coefficient (Davies) as the documented-less-reliable fallback for any species without a citable literature parameter (most custom organic Species values – amino acid states, antibiotics, etc. – fall back this way; only the well-characterized inorganic ions a lab module registers get the upgrade).

_alpha_fractions(pKa::Vector{Float64}, pH::Real)

The equilibrium mole fraction of each species in an AcidBaseSystem (fully-protonated first) at the given pH – _beta_products normalized to sum to 1.0. Returns a Vector{Float64} of length length(pKa)+1.

_analytical_species(s::Stock)

Partition s's composition into the equilibrium solver's inputs: one AnalyticalSpecies per distinct AcidBaseSystem registered on any reagent actually present in s (deduplicated, since e.g. potassium and sodium phosphate salts share the same phosphate AcidBaseSystem), and a list of (chemical, concentration) pairs for every other charged Species in the composition – excluding H⁺/OH⁻ themselves (handled as the solver's own pH-dependent unknowns) and any chemical that's part of a family (handled via that family's speciation instead). Carrying the Species itself (not just its charge) lets the activity-correction machinery look up per-species Truesdell-Jones parameters (see _activity_coefficient) rather than only a shared per-charge-magnitude value.

_beta_products(pKa::Vector{Float64}, pH::Real)

The equilibrium concentration ratio [speciesⱼ]/[species₁] for each species in an AcidBaseSystem (beta[1]=1.0 for species₁ itself), computed in log10-space for numerical stability across wide pKa spans (e.g. phosphoric acid's pKa1≈2.1 to pKa3≈12.4). Unnormalized – see _alpha_fractions for the normalized (sum-to-1) mole-fraction version used when species₁'s concentration isn't itself fixed; this unnormalized form is what an OpenSystemSpecies needs, since there [species₁] (not the family total) is the fixed quantity.

_bisect_pH(f; lo=-2.0, hi=16.0, tol=1e-9, maxiter=200)

Small self-contained bisection root-finder for pH's charge-balance objective f (a monotonically decreasing function of pH, bounded on [lo,hi]) – deliberately not Newton's method (no derivative needed, guaranteed convergence given a valid bracket, no divergence risk with near-degenerate pKa's) and deliberately not a dependency on Roots.jl (this objective's shape doesn't need it). Emits a warning and returns a boundary estimate if f doesn't change sign across [lo,hi] (an extreme edge case), rather than erroring.

_bracket_titrant(f; n0, growth=2.0, max_expand=60)

Expand outward from f(0u"mol") until f changes sign, establishing a bracket (lo,hi) (molar amounts) for adjust_pH's bisection – mirrors _bisect_pH's own guaranteed-convergence-given-a-valid-bracket philosophy, but for "how much titrant" rather than "what pH." Throws an error (not a warning) if no sign change is found within max_expand doublings, since returning a stock that doesn't actually reach the target would be actively misleading, unlike _bisect_pH's own boundary-estimate fallback for an already-bounded pH range.

_charge_balance_residual(pH, families, pKa_eff, strong_net, Kw_eff)

The generalized electroneutrality condition pH is root-found against: [H⁺] - [OH⁻] + Σ_families Σⱼ [speciesⱼ](pH)·chargeⱼ + strong_net, all as dimensionless mol/L values – computed via _species_concentrations, so this works identically for a fixed-total AnalyticalSpecies family and a fixed-free-species OpenSystemSpecies one. Reduces algebraically to CHESSCore's original strong-electrolyte-only net_hydrogen_ion_concentration when families is empty (see pH's fast path, which uses that exact original formula directly rather than this general one).

_check_davies_reliability(pH::Real, I::Float64, families, strong_ions, pKa_eff)

Warn if I exceeds a charge-scaled threshold (0.5/z², following directly from Davies' own z² correction term – an approximation of where the correction becomes unreliable, not a literature-cited cutoff) for the highest charge magnitude materially and still-Davies-corrected present – i.e. a species contributing more than 1% of the total ionic strength I at the solved pH (via _species_concentrations, not just appearing somewhere in a family's chain definition) AND with no registered ion_parameters (real Truesdell-Jones species – e.g. CO3²⁻, Ca²⁺ once a lab module registers them – are exempt, since they no longer use Davies at all). This distinction matters for e.g. the carbonic acid family: its chain includes CO3²⁻, but at near-neutral pH that species' actual concentration is negligible, and checking chain membership alone would false-positive on every water_correction=true call regardless of whether CO3²⁻ is actually present in any meaningful amount. A no-op when everything materially present is monovalent, neutral, or already using real Truesdell-Jones parameters.

_conditional_pKa(sys::AcidBaseSystem, gammas::Dict{Species,Float64})

Convert sys's thermodynamic (activity-basis) pKas to concentration-basis "conditional" pKa's using the current per-species activity coefficients gammas (keyed by Species identity – see _activity_coefficient; real Truesdell-Jones values where registered, Davies fallback otherwise). The thermodynamic constant is Kaⱼ = {H⁺}{speciesⱼ₊₁}/{speciesⱼ} = γ(H⁺)[H⁺]·γ(speciesⱼ₊₁)[speciesⱼ₊₁] / (γ(speciesⱼ)[speciesⱼ]); solving for the concentration-basis Ka'ⱼ = [H⁺][speciesⱼ₊₁]/[speciesⱼ] gives Ka'ⱼ = Kaⱼ·γ(speciesⱼ)/(γ(H⁺)·γ(speciesⱼ₊₁)), i.e. pKa'ⱼ = pKaⱼ + log10(γ(H⁺)·γ(speciesⱼ₊₁)/γ(speciesⱼ)). Keying by species identity (not just charge magnitude) means two same-charge species with different registered ion-size parameters get correctly different shifts here, rather than being forced to share one value.

_dose(r::Reagent, n::Unitful.Amount)

Build the Stock contribution of n moles of r, regardless of whether r is a Solid (molar quantities are already supported directly, see Stocks.jl) or a Liquid (only volume-constructible there, so n is converted to a volume via molecular_weight/density first). Used by adjust_pH to titrate with either kind of reagent uniformly. Missing molecular_weight/density surface the same errors convert already raises elsewhere; Gas titrants aren't supported.

_ionic_strength(pH, families, strong_ions, gammas, pKa_eff)

Ionic strength I = 0.5·Σcᵢzᵢ² at the given pH, over every charged species actually present: [H⁺], [OH⁻] (using the activity-corrected Kw_eff = Kw/(γ(H⁺)·γ(OH⁻)) – each its own per-species γ from gammas, not one shared value, since H⁺ and OH⁻ are different species with different real ion-size parameters), every strong ion, and every family state (via pKa_eff, so this stays self-consistent with whatever pKa's were actually used to solve for pH).

_pH_strong_electrolyte_only(net_conc)

CHESSCore's original strong-electrolyte-only pH formula (net_conc = net H⁺-minus-OH⁻ molar concentration): 7.0 if zero, -log10(net_conc) if net H⁺, 14 - (-log10(-net_conc)) if net OH⁻. pH(::Stock)'s dedicated fast path when no acid/base family is present, so the common strong-electrolyte case is guaranteed bit-identical to the pre-equilibrium-solver behavior rather than merely "solved and hoped to land in the same place."

_relevant_chemicals(families, strong_ions)

Every distinct Species whose activity coefficient the solver needs each iteration: H⁺, OH⁻, every species across every family's chain (including neutral ones, which always resolve to γ=1.0 – see _activity_coefficient), and every strong-ion chemical.

_solve_ph(families, strong_ions; ionic_strength_correction=true)

Shared solver core behind pH and speciation: root-finds pH via _bisect_pH, then, if ionic_strength_correction, iterates the self-consistent activity-coefficient loop (recompute ionic strength from the current speciation, recompute each species' γ – Truesdell-Jones or Davies, per _activity_coefficient – recompute conditional pKa's/Kw, re-solve – typically converging within a handful of iterations) until γ stabilizes or a small iteration cap is hit (in which case a warning is emitted and the best estimate so far is returned, rather than erroring). Returns (pH, pKa_eff) – the second so callers needing per-species fractions (speciation) can report ones consistent with the pH actually solved for.

_species_concentrations(fam::AbstractAnalyticalSpecies, pH::Real, pKa_eff::Dict{AcidBaseSystem,Vector{Float64}})

The concentration of every species in fam.system at the given pH, as a plain Vector{Float64} (mol/L, dimensionless): total_concentration .* alphas for a fixed-total AnalyticalSpecies, or free_concentration .* betas for a fixed-free-species OpenSystemSpecies – the shared interface _charge_balance_residual/_ionic_strength/speciation consume so neither needs to branch on which kind of family it's looking at.

acid_base_system(r::Reagent)

Return the AcidBaseSystem registered for r, or nothing if r is inert with respect to proton transfer (the default for any reagent not registered via set_acid_base_system!).

activity_coefficient(z::Integer, ionic_strength::Real)

The Davies-equation activity coefficient for an ion of charge z at the given ionic_strength (dimensionless mol/L value, i.e. already ustripped): log10(γ) = -0.509 z² (√I/(1+√I) - 0.3I) at 25°C. γ is always 1.0 for a neutral species (z==0), matching the convention that activity corrections only apply to charged species.

Reliability is charge-dependent, not just concentration-dependent: because of the z² term, a divalent ion receives 4x the correction of a monovalent ion at the same ionic strength, and a trivalent ion 9x – so Davies degrades for |z|≥2 well within the ~0.1-0.5 mol/L range it's nominally reliable for monovalent ions (a well-documented limitation of Davies specifically, not an implementation quirk here). pH(::Stock) surfaces this via a warning (see _check_davies_reliability) rather than silently returning a number that may be off by a full pH unit or more for polyvalent systems at moderate-to-high concentration – confirmed empirically against real measurements (bicarbonate/succinate-heavy stocks in CHESSLabConstants's empirical measurement library, CHESSLabConstants/test/empirical_measurements.jl).

adjust_pH(s::Stock, target_ph::Real, acid::Reagent, base::Reagent;
          ionic_strength_correction=true, tol=1e-4, max_iter=60) -> Stock

Return a new Stock – s plus however much of acid or base is needed – whose pH is within tol of target_ph (must be in [1,14]). s itself is never modified: this falls out of Stock's existing immutable design (+/* always construct new values, see Stocks.jl), the same way every other stock-combining operation in CHESSCore already works.

Compares pH(s) to target_ph to decide direction – acid if s is currently more basic than the target, base if more acidic – then root-finds the required molar amount via an expanding bracket (see _bracket_titrant) followed by bisection, re-evaluating pH(s + n*titrant) at each trial n (through _dose, so acid/base can be either a Solid or a Liquid). This reuses pH(::Stock) exactly as-is for every evaluation, so it works for any registered acid/base chemistry – strong or weak titrant, buffered or unbuffered s – with no separate solver logic.

If s is already within tol of target_ph, returns s unchanged (no titrant needed).

ionic_strength_correction is forwarded to every internal pH call, matching pH(::Stock)'s own keyword. If max_iter bisection steps aren't enough to reach tol, emits a warning and returns the best candidate found, matching the warn-and-degrade style used throughout this solver (e.g. _bisect_pH, _solve_ph) – rather than erroring outright, since a slightly off-target result is still useful.

pH(families::Vector{<:AbstractAnalyticalSpecies}, strong_ions; ionic_strength_correction=true)

The general, low-level entry point behind pH(::Stock): root-finds the solution pH of an arbitrary mixture of acid/base families (each an AnalyticalSpecies or OpenSystemSpecies) plus strong_ions (a list of (chemical, concentration) pairs for fully-dissociated strong electrolytes – carrying the Species itself, not just its charge, so per-species activity parameters can be looked up; see _activity_coefficient). Set ionic_strength_correction=false to get the ideal/infinite-dilution pH instead (no activity correction).

set_acid_base_system!(r::Reagent,sys::AcidBaseSystem)

Register that r participates in acid/base equilibrium according to sys once dissolved – see acid_base_system, pH.

set_default_water_correction!(spec::OpenSystemSpecies)

Register spec as the default open-system correction pH(::Stock) applies when water_correction=true and no water_correction_source is passed explicitly.

set_ion_parameters!(chem::Species, a::Float64, b::Float64)

Register chem's Truesdell-Jones ion-size parameter a (å, Ångströms) and empirical b parameter – see ion_parameters, _activity_coefficient.

speciation(s::Stock; ionic_strength_correction=true)

Per-family diagnostic breakdown of s's acid/base equilibrium at its solved pH: a SpeciationResult for each distinct AcidBaseSystem present, giving the equilibrium fraction and concentration of every protonation state. Returns an empty vector if s has no registered acid/base families (pure strong-electrolyte stocks have nothing to speciate).