Surfaces API

How much area a solid offers, and to what. One abstraction serves the reactive area of a dissolution rate law and the binding area that carries surface sites: the three areas of a problem — reactive, accessible, binding — are different numbers, but they are the same kind of object, and the measurement each comes from is carried by its type.

Areas

ChemistryLab.SHRINK_FLOOR Constant

julia

SHRINK_FLOOR

Fractional floor regularizing the power law of ShrinkingCoreArea near exhaustion, 1e-8.

What it is for, and why a floor rather than a branch

f^p with a fractional p has an infinite derivative at f = 0, so a solid that an integrator drives to exhaustion produces an infinite Jacobian entry and an automatic derivative that is Inf or NaN. The regularized form used here,

is exactly zero at f = 0, exactly one at f = 1, has a finite derivative everywhere, reduces to g(f) = f when p = 1, and differs from f^p by a relative f_c/(3f) at p = 2/3 — about 3e-9 at half depletion. It is written without a comparison so that nothing branches on a value that may be a Dual.

ChemistryLab.AbstractSpecificArea Type

julia

abstract type AbstractSpecificArea <: AbstractSurfaceModel end

A surface model defined by a specific area in m²/kg, together with the measurement it comes from.

The measurement is carried by the type, not by a field, and that is the whole point: area_ratio is defined only between two models of the same type, so comparing a BET area with a Blaine fineness is a MethodError-free, named refusal rather than a plausible wrong number.

Concrete subtypes implement specific_area and area_method, and inherit total_area as specific_area(m) * max(n, 0) * M.

ChemistryLab.AbstractSurfaceModel Type

julia

abstract type AbstractSurfaceModel end

How much area a solid offers, in m², as a function of how much of it is left.

Concrete subtypes implement

julia

total_area(model, n::Real, n₀::Real, M::Real) -> Real

where n is the current molar amount [mol], n₀ the initial one [mol], and M the molar mass [kg/mol]. The result is a total area in m², never a specific one.

Why n₀ is in the signature from the start

Because an area that cannot see where it started cannot describe a microstructure that evolves. n alone expresses "area proportional to the remaining mass"; the ratio n/n₀ is what expresses a shrinking grain (ShrinkingCoreArea). Passing both costs nothing and means no rate law has to change its signature when a richer area model is added later.

surface_area is the three-argument form kept for callers that have no initial amount to give; it is total_area(model, n, n, M).

All methods are AD-compatible: no Float64 conversion, and the element type of the result follows its arguments.

ChemistryLab.BETSurfaceArea Type

julia

struct BETSurfaceArea{T<:Real} <: AbstractSpecificArea

A specific area in m²/kg from a BET (Brunauer-Emmett-Teller) gas-adsorption measurement.

julia

A = A_specific × n × M        [m²]

the standard route in reactive-transport modeling (Palandri and Kharaka, 2004). Using it presumes the area accessible to the dissolving solute is the one the adsorbed gas saw, which is an assumption about the pore structure, not a measurement of it.

Fields

A_specific: specific BET area [m²/kg].

Examples

julia

julia> BETSurfaceArea(90.0)
BETSurfaceArea{Float64}(90.0)

julia> BETSurfaceArea(0.09u"m^2/g")
BETSurfaceArea{Float64}(90.0)

ChemistryLab.BlaineSurfaceArea Type

julia

struct BlaineSurfaceArea{T<:Real} <: AbstractSpecificArea

A specific area in m²/kg from a Blaine air-permeability fineness.

Why this is a separate type from BETSurfaceArea

Because the two measure different things, and the package already says so in prose three times over. Silica fume is about 20 000 m²/kg by BET, while the effective Blaine fineness recommended for its Waller kinetics (Waller, 1999) is about 2 000 m²/kg — a factor of ten. Passing the first where the second is expected produces a hydration rate ten times too fast, with nothing to catch it.

Making the measurement a type is what turns that warning into a refusal: area_ratio between a Blaine and a BET area raises, because no method covers the pair.

Fields

A_specific: Blaine fineness [m²/kg].

Examples

julia

julia> BlaineSurfaceArea(385.0)
BlaineSurfaceArea{Float64}(385.0)

julia> BlaineSurfaceArea(380u"m^2/kg")
BlaineSurfaceArea{Float64}(380.0)

See also: blaine_factor, area_ratio.

ChemistryLab.FixedSurfaceArea Type

julia

struct FixedSurfaceArea{T<:Real} <: AbstractSurfaceModel

A total reactive area in m², independent of how much solid is left.

Suitable for a short run, or when the area is controlled externally — a measured powder in a reactor whose area is not expected to move over the experiment.

Fields

A: total area [m²].

Examples

julia

julia> FixedSurfaceArea(0.5)
FixedSurfaceArea{Float64}(0.5)

julia> FixedSurfaceArea(500.0u"cm^2")
FixedSurfaceArea{Float64}(0.05)

See also: BETSurfaceArea, total_area.

ChemistryLab.FixedSurfaceArea Method

julia

FixedSurfaceArea(A) -> FixedSurfaceArea

Build a FixedSurfaceArea from a total area.

A is a plain Real in m², or a Quantity convertible to m². The element type is kept, so a ForwardDiff.Dual area differentiates through.

ChemistryLab.GeometricSurfaceArea Type

julia

struct GeometricSurfaceArea{T<:Real} <: AbstractSpecificArea

A specific area in m²/kg computed from a geometric idealization — a particle size distribution and an assumed shape — rather than measured by adsorption or permeability.

Kept distinct from the two measured kinds for the same reason they are distinct from each other: a geometric area is typically orders of magnitude below a BET area on the same powder, because it counts no internal roughness.

Fields

A_specific: geometric specific area [m²/kg].

ChemistryLab.ShrinkingCoreArea Type

julia

struct ShrinkingCoreArea{A<:AbstractSpecificArea, T<:Real} <: AbstractSurfaceModel

An area that follows a population of grains being consumed, rather than the mass that remains.

with 𝒜₀ the initial area computed from initial, and g the regularized power law described under SHRINK_FLOOR.

Choosing the exponent

p = 2/3 is the geometric answer for spheres consumed from the outside: area goes as the square of a radius, amount as its cube. p = 1 recovers an area proportional to the remaining mass, that is, exactly what a plain BETSurfaceArea gives — so the two are one family, and the exponent is the modeling choice that separates them. Values in between are common fits and should be labeled as fits.

What this does not do

It does not describe an area blocked by something growing on it — a hydrate layer covering a clinker grain reduces the accessible area without consuming the grain, so it depends on the whole assemblage and not on n alone. That is a different model, and the precedent for writing it is PoreHumidity, which already reads the current composition to produce a geometric quantity.

Fields

initial: the specific-area model giving the area at n = n₀, carrying its own measurement method.
exponent: p, dimensionless.

Examples

julia

ShrinkingCoreArea(BETSurfaceArea(90.0), 2 // 3)
ShrinkingCoreArea(GeometricSurfaceArea(0.4u"m^2/g"))      # p = 2/3 by default

Size classes of a powder

ChemistryLab.AbstractPSD Type

julia

AbstractPSD

A particle size distribution by volume: the fraction F(d) of the volume of the powder held by particles of diameter below d. Its subtypes are LogNormalPSD and TabulatedPSD, and size_classes cuts one into classes.

ChemistryLab.LogNormalPSD Type

julia

LogNormalPSD(; d50, d90) -> LogNormalPSD
LogNormalPSD(; d50, σ) -> LogNormalPSD

A size distribution by volume whose diameter is lognormal,

with Φ the standard normal distribution function. d50 and d90 are lengths, a Quantity or a plain Real in m. Differentiable in both, and generic in their number type.

ChemistryLab.SizeClass Type

julia

SizeClass

One class of a ground mineral, as size_classes returns it.

Fields

species: the species of the class — the mineral itself for the first, a copy under its own symbol (with_symbol) for the others;
amount: its initial amount, mol;
surface: its area model, ShrinkingCoreArea(GeometricSurfaceArea(s));
diameter: its Sauter mean diameter, the diameter of the sphere with the class's ratio of volume to surface, m;
fraction: the fraction of the volume of the powder it holds.

ChemistryLab.TabulatedPSD Type

julia

TabulatedPSD(d, F) -> TabulatedPSD

A measured size distribution by volume: diameters d in ascending order (a vector of lengths, Quantity or plain Real in m) and the cumulative volume fraction F passing at each, from zero at the first diameter to its largest value at the last, in percent or as a fraction (it is divided by its last value). Between two points F is taken linear in ln d, which is how a laser granulometer reports it.

ChemistryLab.inverse_diameter_moment Method

julia

inverse_diameter_moment(psd, Fa, Fb) -> Real

∫ (1/d) dF over the particles between the volume fractions Fa and Fb, in 1/m: times 6/ρ, the specific surface those particles add as spheres. Exact for both distributions — a closed form for LogNormalPSD,

and, for TabulatedPSD, the integral of e^{-u} over F linear in u = ln d between points — so that the classes of a distribution add up to its whole surface to rounding.

ChemistryLab.size_classes Method

julia

size_classes(species, psd::AbstractPSD, K; amount | mass, density = nothing,
             T = 298.15, P = 1e5, shape = 1, names = k -> …) -> Vector{SizeClass}

Cut a ground mineral into K classes of equal volume of its size distribution psd, each ready to receive the same rate law.

Class k holds the volume fractions ((k−1)/K, k/K), a K-th of the powder, and the particles in it, taken as spheres scaled by shape, have the specific surface

d_k the Sauter diameter of the class (inverse_diameter_moment). The classes add up to the amount given and to the specific surface of the whole distribution, both to rounding; one class is the mineral itself with that whole surface. Each class is a SizeClass: its species, its amount, and the area model ShrinkingCoreArea(GeometricSurfaceArea(s_k)), a grain shrinking from the outside.

Keywords

amount or mass: the whole powder, in mol or in kg (a Quantity or a plain Real in SI);
density: in kg/m³; by default the molar mass over the molar volume of the species, M/V⁰(T, P), which it must then carry;
shape: the ratio of the particles' surface to that of spheres of the same volume, one by default;
names: the symbol of class k, by default the mineral's own for the first and "$(symbol)_$k" for the others.

Example

julia

classes = size_classes(subs["C3S"], LogNormalPSD(; d50 = 10u"μm", d90 = 30u"μm"), 4;
                       mass = 1u"kg")
species = vcat(aqueous, [c.species for c in classes])
# one reaction per class, the same law on each class's own surface:
laws = [transition_state(mechs, cs, dissolution(c.species), c.surface) for c in classes]

Differentiable in the parameters of the distribution and in the amount.

Site families

A site is a conserved quantity that is not a chemical element, and it is carried the way electric charge already was: as a pseudo-element in a species' formula, turned into a conservation row by the ordinary matrix assembly.

ChemistryLab.AbstractSiteCapacity Type

julia

abstract type AbstractSiteCapacity end

How many moles of sites a family offers.

Concrete subtypes implement

julia

site_moles(capacity, support, n_host::Real, n_host₀::Real, M_host::Real) -> Real

returning moles of sites.

Why three of them, and not one

Because the published data comes in three shapes, and converting between them needs a quantity that is usually not measured. A specific area with a site density per square meter (AreaSiteDensity) is the oxide literature's form; a capacity per kilogram of dry solid (MassSiteDensity) is the clay literature's, and turning the second into the first would mean inventing a BET area to divide by. A prescribed total (TotalSiteAmount) is what a fixed sorbent in a batch experiment actually gives.

See also: SiteFamily, site_moles.

ChemistryLab.AbstractSiteMixingModel Type

julia

abstract type AbstractSiteMixingModel end

How the species of a SiteFamily mix on their shared budget of sites.

Concrete subtypes implement

julia

_site_excess_ln_gamma(model, k, x, T) -> Real

the departure from ideality of the k-th member, given the site fractions x of the whole family. The ideal part, ln x_k, is added by the caller, exactly as it is for a solid solution.

The models provided are IdealSiteMixing, the two exchange conventions VanselowMixing and GainesThomasMixing, and the electrostatic wrappers ConstantCapacitance and DiffuseLayer. The others the literature uses — Frumkin's interaction term, the quasi-chemical approximation for a multidentate adsorbate — are separate models with their own parameters and their own validation, and this is where they would attach.

See also: IdealSiteMixing, SiteFamily.

ChemistryLab.AreaSiteDensity Type

julia

struct AreaSiteDensity{T<:Real} <: AbstractSiteCapacity

Sites per unit area, in mol/m². The area comes from the family's SurfaceSupport, so the measurement it was made with travels with it.

Fields

Γ_C: site density [mol/m²].

ChemistryLab.ChargePlanes Type

julia

struct ChargePlanes{M<:AbstractSiteMixingModel, T<:Real} <: AbstractSiteMixingModel

A charged surface described by three planes of charge and a diffuse layer: the surface plane 0, where the protons and the inner-sphere complexes sit, a plane 1 at the head of the Stern layer, and a plane 2 where the diffuse layer begins. Each surface species puts a charge on each plane, and the planes are separated by two capacitors,

with , the charge density of plane p and the charge species k places on it, and as in DiffuseLayer. The electrical work of species k is .

One description covers the models of the literature, by the charges the species are given and the capacitances:

the triple-layer model of Davis et al. (1978): protons and inner-sphere complexes on plane 0, the ion pairs of the background electrolyte on plane 1, two finite capacitances;
the basic Stern model (Westall and Hohl, 1980): the ion pairs on the plane where the diffuse layer begins, C2 = Inf;
the charge distribution model of Hiemstra and Van Riemsdijk (1996): the charge of an inner-sphere complex shared between planes 0 and 1 in fractions, the singly and triply coordinated surface groups of a mineral carrying fractional charges of their own.

The charges a species carries are its :plane_charges property, three numbers, given with with_plane_charges; a species without it puts its formal charge on plane 0. The intrinsic charge of a MUSIC site, −1/2 on a singly coordinated group of goethite, belongs there too: it raises the potentials, while the formal charges, whole numbers, keep the conservation of charge, which the constant difference between the two does not change.

The potentials are carried as unknowns of the solve, one per plane and per surface, closed by the three equations above; a diffuse layer makes the electrical term depend on the ionic strength of a solution that does not depend on the surface in return, so, as with DiffuseLayer, what a solve returns is a self-consistent speciation rather than a certified minimum.

Fields

base: the site mixing this decorates, IdealSiteMixing by default.
area: the charged area [m²].
C1, C2: the capacitances of the inner and outer Stern layers [F/m²]; C2 = Inf merges planes 1 and 2.
ε_r, water, pressure: the relative permittivity of the solvent, for the diffuse layer, as for DiffuseLayer: water's at the temperature of each solve unless a value is given.

ChemistryLab.ChargePlanes Method

julia

ChargePlanes(; area, C1, C2 = Inf, temperature = 298.15, pressure = 1.0e5,
               ε_r = nothing, base = IdealSiteMixing()) -> ChargePlanes

Build a ChargePlanes surface: area in m², C1 and C2 in F/m², each a plain Real in SI or a Quantity. Without ε_r, the diffuse layer is screened with the permittivity of water at the temperature of each solve, as for DiffuseLayer.

ChemistryLab.ConstantCapacitance Type

julia

struct ConstantCapacitance{M<:AbstractSiteMixingModel, T<:Real} <: AbstractSiteMixingModel

A charged surface, in the constant-capacitance model: the potential is proportional to the charge the surface carries.

with C the capacitance in F/m², 𝒜 the surface area in m², and z_k the formal charge of each member. It decorates another mixing model rather than replacing one: the site fractions are still whatever base says they are, and this adds the electrical work of putting a charge on a charged surface.

Why this needs no unknown of its own

The literature presents an electrostatic surface model as one extra unknown per surface, Ψ, with one extra equation to close it. It is not needed here: σ = CΨ makes Ψ an explicit function of the composition,

so the whole model is a composition-dependent term in the chemical potential — which is what an activity coefficient is. It belongs with the mixing, not with the constraints, and putting it there is what keeps the solver unchanged.

The convexity, written rather than assumed

The electrical work of charging the surface is G_el = ∫₀^σ Ψ(s)\,ds \cdot 𝒜, and with Ψ = σ/C that integrates to

a quadratic form in n with Hessian \frac{F^2}{C\mathcal{A}} z z^{\mathsf T}, positive semi-definite for any positive capacitance. So this term is convex, the site mixing it decorates is convex, and the equilibrium certificate covers the sum unchanged. That is a proof, not an expectation, and it is the reason this model could be added without reopening the certificate.

Its gradient is ∂G_el/∂n_j = z_j F Ψ, the electrochemical work the physics asks for — which is the check that the energy above is the right one.

How far the solve reaches, measured

Convexity makes the minimum unique, so any failure to find it is numerical and not a second answer. And there is one: eliminating Ψ puts the whole electrostatic stiffness into the composition dependence of an activity, and the Newton loses it when that stiffness grows.

The scale of it is dimensionless and worth computing before a run:

N being the site budget — the potential the surface would reach with every site charged. Measured on hydrous ferric oxide, N = 2·10⁻⁴ mol on 53.4 m²:

`C` [F/m²]	`ψ̃_max`	stationarity of the solve
10	1.4	3e-16
5	2.8	2e-16
3	4.7	2e-16
2	7.0	0.04 — lost
1.2	11.7	0.07 — lost

So ψ̃_max ≲ 5 is reached directly, and stepping C down from a large value while reusing the previous composition extends the range but does not remove the limit. Oxide capacitances of 1–3 F/m² therefore sit at the edge of it.

What lifts it is the formulation, not the tolerance. Carrying Ψ as an unknown of the solve with σ = CΨ as its closing equation is mathematically the same problem — that is what eliminating it proved — but the Newton then controls the potential directly instead of meeting it through a stiff exponential. The literature's extra unknown is a preconditioner, and it is what this model is missing rather than a correction to it.

Fields

base: the site mixing this decorates, usually IdealSiteMixing.
C: capacitance [F/m²]. Values of 1–3 F/m² are the usual range for an oxide.
area: the charged area [m²].

What it does to a titration

A surface that has already taken protons resists taking more, because the work of adding a charge to an object that is already charged grows with the charge. The visible effect is a flattened titration curve: the transitions spread over more pH units than the constants alone would give. A set of constants fitted with an electrostatic term and used without one — or the reverse — therefore describes a different surface.

See also: IdealSiteMixing, SiteFamily.

ChemistryLab.ConstantCapacitance Method

julia

ConstantCapacitance(base, C, area) -> ConstantCapacitance
ConstantCapacitance(; C, area, base = IdealSiteMixing()) -> ConstantCapacitance

Build a ConstantCapacitance. C is a capacitance in F/m² and area an area in m², each a plain Real in SI or a Quantity.

ChemistryLab.DiffuseLayer Type

julia

DiffuseLayer(base, area, ε_r, scale = 1.0) -> DiffuseLayer
DiffuseLayer(; area, temperature = 298.15, pressure = 1.0e5, ε_r = nothing,
               scale = 1.0, base = IdealSiteMixing()) -> DiffuseLayer
DiffuseLayer(; specific_area, ...) -> DiffuseLayer

Build a DiffuseLayer. area is in m², a plain Real in SI or a Quantity. specific_area, given instead, is an area per mole or per kilogram of the host of the family's support — a Quantity in m²/mol or m²/kg, or a plain Real in m²/kg — and the area then follows the host (see the type). Without ε_r, the permittivity is that of water at the temperature of each solve and at pressure, from this package's own model, as PHREEQC takes it; temperature only sets the value shown. A given ε_r (and the positional form) is held at every temperature. scale multiplies the electrostatic term, as with_electrostatic_scale sets it.

The permittivity of water costs 7 µs to evaluate (measured), so following the temperature costs a solve nothing worth counting.

ChemistryLab.DiffuseLayer Type

julia

struct DiffuseLayer{M<:AbstractSiteMixingModel, T<:Real, A} <: AbstractSiteMixingModel

A charged surface, in the diffuse-layer (Gouy-Chapman) model: the potential is raised by the charge the surface carries and screened by the ions in solution, so it depends on the ionic strength as well as on the charge.

with I the molal ionic strength of the aqueous solution, ρ = 1000 kg/m³ the factor that turns it into a volumetric concentration, 𝒜 the surface area in m², and ε_r the relative permittivity of the solvent. Like ConstantCapacitance it decorates another mixing model: the site fractions stay whatever base says, and this adds the electrical work.

This is the model behind the calibration of Dzombak and Morel (1990) and the default of PHREEQC's SURFACE block — which is what makes it checkable against another code rather than only against itself.

It needs no unknown of its own either, and that is not obvious

A diffuse layer is usually presented as one extra unknown Ψ per surface with one extra equation to close it, on the grounds that the relation above is transcendental. It is transcendental in Ψ, but it is monotone in Ψ, so it inverts in closed form:

asinh is smooth and bounded in its derivative everywhere, so the term costs the solver nothing and differentiates cleanly. κ is computed from ε_r, and ε_r from this package's own Johnson-Norton model of water rather than typed in: at 25 °C that gives κ = 0.117215, against the 0.1174 PHREEQC writes into its source. Agreeing with PHREEQC while using a different dielectric constant is a stronger statement than agreeing with its arithmetic.

What this model costs the certificate

Unlike ConstantCapacitance, this term is not the gradient of any Gibbs energy, and no implementation can make it one. At fixed I it is: the electrical work integrates to

whose derivative in n_k is exactly RT z_k \tilde\psi. But I is itself a function of the aqueous composition, and the aqueous activity of an ion does not depend in return on how much is bound to the surface. The Jacobian of the activity map is therefore asymmetric, and an asymmetric Jacobian is not the Hessian of anything.

That is a property of the Dzombak-Morel model, not of this code: it treats the bulk solution as a reservoir whose ionic strength is a parameter, which is precisely the approximation that lets the diffuse layer be written without carrying its ion inventory. PHREEQC's default SURFACE makes the same one.

Here the consequence is reported. is_gradient_consistent returns false for this model, a system that uses it says so, and what comes back from a solve is a self-consistent speciation — mass action and conservation satisfied together — not a certified minimum. Everything the other models certify, they still certify; this one buys agreement with the published calibrations at that price, and the price is written on it.

Fields

base: the mixing model this decorates, IdealSiteMixing() by default.
area: the surface area carrying the charge — a number of m², or, for a layer given a specific_area, the area per mole of sites of its support that the system resolves it to (see below).
ε_r: the relative permittivity of the solvent, dimensionless. Unless a value is given, the permittivity of water at the temperature of each solve, as PHREEQC takes it (at 50 °C, 69.8 against 78.2 at 25 °C); ε_r then holds its value at the construction temperature, for display. A given value is held at every temperature.
water: whether the permittivity is that of water at the temperature of the solve, and pressure the pressure it is evaluated at, in Pa.

A surface whose amount changes

A surface whose host precipitates or dissolves during a calculation carries a charge per unit area that is a property of the phase, so its area has to follow the host as its sites do. Give the layer a specific_area, per mole or per kilogram of the host of its support, instead of an area. The support must follow its host (SITES_FOLLOW_HOST) with a budget proportional to it, ν moles of sites per mole of host (sites_per_host), and every member must occupy one site. When the system is built, the specific area a becomes an area per mole of sites of the support, s = a/ν, and the area at every solve is

the sums running over the members of every family on the support. At the solution Σₖ nₖ = ν n_host, so 𝒜 = a·n_host, and the charge per unit area is a function of the site fractions alone: it does not change when the surface grows or shrinks at a fixed occupancy. The area is read from the site amounts rather than from the host's, so it holds whether the host is a species of the equilibrium, a kinetic species of a KineticsProblem, or both.

ChemistryLab.GainesThomasMixing Type

julia

struct GainesThomasMixing <: AbstractSiteMixingModel

Cation exchange in the Gaines-Thomas convention: the activity of an exchanger species is its equivalent fraction, counting the charge it compensates.

with z_i read from the formula as the number of the family's pseudo-elements the species carries — which on a permanent-charge exchanger is the charge it neutralizes.

The two conventions are not interchangeable, and the conversion is not a factor

For a homovalent exchange, Na⁺/K⁺, the two fractions are proportional and the selectivity coefficients agree. For a heterovalent one, Na⁺/Ca²⁺, they do not: one calcium and one sodium are one particle each but one and two charges, so x_i and E_i diverge, and so do the constants fitted under each.

The literature sometimes quotes conversion factors of 2, 3 or 4. Those are trace-composition limits, not constant offsets: the exact relation depends on the exchanger's composition, which is what the calculation is solving for. This package therefore declares the convention and converts nothing implicitly. A constant fitted under one convention used under the other is a different model.

ChemistryLab.IdealSiteMixing Type

julia

struct IdealSiteMixing <: AbstractSiteMixingModel

Occupied and free sites mix ideally: the activity of a member is its site fraction, a_j = n_j / N_t.

What this is, and what it is not

It is the production form of the Langmuir model, and Langmuir falls out of it rather than being imposed: eliminating the free site from the equilibrium of n_free + Σ n_j = N_t gives n_j/N_t = β_j /(1 + Σ_k β_k) with β_j = K_j a_j, saturation and competition included, with no isotherm written anywhere.

It is not a place to also apply a surface activity coefficient of the 1/(1 − θ) kind. That factor is what an eliminated-free-site formulation needs in order to recover what this one already produces; applying both counts the same physics twice. The ratio is available as a diagnostic — see the tests — and the identity between the two forms is asserted rather than assumed.

ChemistryLab.MassSiteDensity Type

julia

struct MassSiteDensity{T<:Real} <: AbstractSiteCapacity

Sites per kilogram of dry support, in mol/kg — the form a clay exchange capacity is published in, and the one that needs no area at all.

Fields

q: site capacity [mol/kg].

ChemistryLab.SiteFamily Type

julia

SiteFamily(name, free_site, complexes; capacity, support, model = IdealSiteMixing()) -> SiteFamily

Build and validate a SiteFamily.

Every member is requalified to AS_SURFACE / SC_SURFCOMPLEX, so a species built from a database record can be passed as it is.

Refused at construction, each because the alternative is a wrong number rather than an error:

a free site whose formula carries no site symbol, or two members carrying different ones — the family would have no budget to share;
a free site occupying no site, or more than one under a mixing model that does not describe multidentate occupancy;
a complex of denticity zero (it is not a member), or above one under such a model (ideal site mixing does not describe it; see the note on denticity above);
a duplicate member.

ChemistryLab.SiteFamily Type

julia

struct SiteFamily{S<:AbstractSpecies, C<:AbstractSiteCapacity, U<:SurfaceSupport}

One family of surface sites: a group of species that share a finite site budget and mix on it.

A family is to a surface what a SolidSolutionPhase is to a solid: a named group of species with its own mixing. The one thing it has that a solid solution does not is a conservation row — its total is fixed by a budget rather than free — and that row is produced by the ordinary matrix assembly, because every member carries the family's pseudo-element in its formula.

The free site is a species, and it is the reference

XsOH is as much a member as XsOH2+ is. It occupies a site, it takes part in the mixing, and its amount is what makes saturation happen: as free sites run out, adding another bound molecule costs more. Eliminating it gives the familiar Langmuir expression, but the production form keeps it explicit — one species per state of a site, with reaction-consistent standard potentials.

It is also the primary the system should be built on, and that is not a detail of taste. A bare Species("Xs") declared as the site primary is inconsistent in charge with a neutral free site, which flips the exact rank test in StoichMatrix and silently adds a spurious charge component; and it never enters cs.species, so saturation_indices would zero its potential and report a wrong index for every surface species without a word.

Denticity is read from the formula, not declared beside it

How many sites a molecule occupies is the coefficient of the family's symbol in its formula — Xs2OCa occupies two. Declaring it separately would create a way for the two to disagree. Ideal mixing of occupied and free sites is exact only for a monodentate species, so a family under it refuses anything but one; a mixing model that describes multidentate occupancy (supports_multidentate, the Vanselow and Gaines–Thomas conventions) accepts more.

Fields

name: the family's label, free-form ("Hfo_s", "Hfo_w").
site: its pseudo-element, one of SITE_SYMBOLS.
free_site: the unoccupied site species.
complexes: the occupied ones.
capacity: how many moles of sites, as an AbstractSiteCapacity.
support: the SurfaceSupport carrying it.
model: how the members mix on the budget, an AbstractSiteMixingModel.

Examples

julia

ssa = literature_value("DzombakMorel1990", "specific_surface_area")   # 600 m²/g
support = SurfaceSupport("hydrous ferric oxide", "Fe(OH)3", BETSurfaceArea(ssa))
free    = Species("XsOH"; aggregate_state = AS_SURFACE, class = SC_SURFCOMPLEX)
bound   = Species("XsOH2+"; aggregate_state = AS_SURFACE, class = SC_SURFCOMPLEX)
SiteFamily("Hfo_s", free, [bound]; capacity = TotalSiteAmount(5.0e-6u"mol"), support)

ChemistryLab.TotalSiteAmount Type

julia

struct TotalSiteAmount{T<:Real} <: AbstractSiteCapacity

A prescribed total number of moles of sites, independent of how much support is present. The honest description of a batch experiment on a fixed sorbent, and the form the first milestone's fixed-support scope uses.

Fields

N: total sites [mol].

ChemistryLab.VanselowMixing Type

julia

struct VanselowMixing <: AbstractSiteMixingModel

Cation exchange in the Vanselow convention: the activity of an exchanger species is its mole fraction, counting molecules.

A calcium and a sodium on the exchanger count as one particle each, whatever charge they neutralize.

See also: GainesThomasMixing, and the note on conversion there.

ChemistryLab._refuse_stacked_electrostatics Method

julia

_refuse_stacked_electrostatics(base, outer)

Refuse an electrostatic decorator applied to another one.

Stacking a diffuse layer on a constant capacitance is how a Stern or a triple-layer model is drawn, and adding the two potentials is not how it works: the two capacitances belong to different charge planes, and each surface species sits on one plane or the other. Summing them puts every species on both. The plane-resolved model is ChargePlanes; refusing here is what keeps someone from assembling a wrong one out of right parts.

ChemistryLab._support_members Method

julia

_support_members(cs) -> (indices, charges)

Per family, the species indices and the formal charges of every member of every family on its support, in one flat list each and in the same order — what a surface charge density has to be summed over.

ChemistryLab.charge_planes_potentials Method

julia

charge_planes_potentials(model, c, n, I, T) -> (ψ̃₀, ψ̃₁, ψ̃₂)

The dimensionless potentials of the three planes of a ChargePlanes surface carrying the members of amounts n, of plane charges c (one triplet per member), screened by a solution of ionic strength I at T: the diffuse layer first, then the two capacitors inward.

ChemistryLab.denticity Method

julia

denticity(family::SiteFamily, sp::AbstractSpecies) -> Int

How many sites of family one molecule of sp occupies, read from the coefficient of the family's pseudo-element in its formula. Zero when sp does not belong to the family.

ChemistryLab.is_electrostatic Method

julia

is_electrostatic(model) -> Bool

Whether a site mixing model adds the work of charging a surface, on top of whatever mixing it decorates.

true for ConstantCapacitance, DiffuseLayer and ChargePlanes, false otherwise — including for a model that merely decorates one, which is why the predicate exists rather than an isa test at each use.

ChemistryLab.is_gradient_consistent Method

julia

is_gradient_consistent(model) -> Bool

Whether the activity contribution of a site mixing model is the gradient of a Gibbs energy — which is what the equilibrium certificate assumes about every term it certifies.

true for every model here but DiffuseLayer and ChargePlanes, whose potentials depend on the ionic strength of a bulk solution that does not depend in return on the surface, making the activity Jacobian asymmetric. See that model's docstring for why this is the Dzombak-Morel approximation itself rather than a defect of the implementation, and site_gradient_asymmetry for the measurement.

ChemistryLab.name Method

julia

name(family::SiteFamily) -> String
site_members(family::SiteFamily) -> Vector{<:AbstractSpecies}
site_capacity(family::SiteFamily) -> AbstractSiteCapacity
surface_support(family::SiteFamily) -> SurfaceSupport

Accessors of a SiteFamily. members lists the free site first, which is the order the site mixing and the solver's reference member both expect.

ChemistryLab.needs_ionic_strength Method

julia

needs_ionic_strength(model) -> Bool

Whether evaluating a site mixing model requires the ionic strength of the aqueous solution.

DiffuseLayer and ChargePlanes do. The activity closures test this once, when they are built, and skip the ionic-strength sum entirely when no family asks for it — so a system without a diffuse layer pays nothing for the possibility of one.

ChemistryLab.reference_member Method

julia

reference_member(family::SiteFamily) -> AbstractSpecies

The member the mixing is written against, and the one the solver carries instead of inverting — site_members(family) puts it first for that reason.

On an oxide it is the unoccupied site, and family.free_site names it literally: as the surface fills, its amount falls and every occupied state costs more to form, which is where saturation comes from.

On a permanent-charge exchanger there is no unoccupied site at all: every charge is compensated by some cation, and what this returns is the form chosen as the reference of the exchange — usually the abundant monovalent one, Na-X. The field keeps the name free_site from the case it was written for; this accessor exists to say what it means in the case it was not.

ChemistryLab.site_mixing_model Method

julia

site_mixing_model(family::SiteFamily) -> AbstractSiteMixingModel

How the members of family mix on their shared site budget.

ChemistryLab.site_moles Function

julia

site_moles(capacity, support, n_host, n_host₀, M_host) -> Real

Moles of sites the family offers, given how much host solid is present.

n_host and n_host₀ are the current and initial molar amounts of the support species [mol] and M_host its molar mass [kg/mol]. A TotalSiteAmount ignores all three; the other two do not, which is what will let a support that precipitates carry its sites with it.

Who calls this

You do, or the checks of a state against its families that call it for you (declared_site_moles, check_site_budget, host_consistent_state). The solve itself does not, and that is exactly what "the support is fixed" means here: the site budget reaches the calculation through the initial amounts of the site-bearing species, so this is the helper that computes the number you put there, evaluated once, before the solve.

Writing it as a function of the host's amount rather than as a constant is what makes an evolving support a change of when it is called, not of the data model — see the theory chapter §10 for the two things that still stand in the way.

ChemistryLab.site_moles Method

julia

site_moles(family::SiteFamily, n_host, n_host₀, M_host) -> Real

The family's site budget, from its capacity and its support.

ChemistryLab.support_group Method

julia

support_group(cs) -> Vector{Vector{Int}}

For each site family of cs, the indices of every family sharing its support — itself included — identified by the support's name.

A surface carries one potential, not one per family. Ferrihydrite is the case that makes this concrete: the strong and weak sites of Dzombak and Morel (1990) are two families on one oxide, so a proton bound to a weak site charges the same surface a proton bound to a strong site does, and both feel the same Ψ. Computing a potential from one family's members alone would make the two sites electrostatically invisible to each other, which is neither the published model nor the physics.

Empty when cs declares no surface.

ChemistryLab.supports_multidentate Method

julia

supports_multidentate(model::AbstractSiteMixingModel) -> Bool

Whether model describes a species occupying more than one site.

false by default, and deliberately: the site balance holds for any denticity — a species carrying two of the family's pseudo-elements contributes twice to it — but a mixing law does not follow from the balance. Counting the ways a molecule can straddle two neighboring sites is a combinatorial problem with its own answer, and returning a number without having solved it would be worse than refusing.

The ion-exchange conventions are the exception, and not because they solved that problem: on a permanent-charge exchanger the "sites" being counted are units of charge, which a divalent cation neutralizes two of without straddling anything. Their multidentate case is bookkeeping, not combinatorics.

ChemistryLab.water_relative_permittivity Function

julia

water_relative_permittivity(T_K, P_Pa = 1.0e5) -> Real

The relative permittivity (dielectric constant) of liquid water at T_K kelvin and P_Pa pascal, from the model of Johnson and Norton (1991) this package already carries for the HKF activity model (Helgeson et al., 1981).

It is 78.245 at 25 °C and 1 bar, and falls to 66.68 at 60 °C — which is why a surface electrostatic model calibrated at room temperature is not transferable to a hydrating paste without saying so.

It evaluates the water equation of state, in about 7 µs (measured), and a DiffuseLayer that follows the temperature calls it at every evaluation of its potential.

ChemistryLab.with_electrostatic_scale Method

julia

with_electrostatic_scale(model, λ) -> model

The same model with its electrostatic term scaled by λ ∈ [0, 1], λ = 0 switching it off entirely and λ = 1 being the model itself.

This exists for one reason, and it is not a physical one: near the point of zero charge the diffuse-layer term is stiff, and a cold Newton solve walks off it. See ChemistryLab.electrostatic_stiffness for the number, and the manual for the continuation that uses this.

Returns the model unchanged for one that carries no electrostatics, so a continuation loop can be written without asking what it is solving.

ChemistryLab.with_plane_charges Method

julia

with_plane_charges(s::Species, c0, c1 = 0, c2 = 0) -> Species

A copy of the surface species s placing the charges c0, c1 and c2 on the planes 0, 1 and 2 of a ChargePlanes surface: the intrinsic charge of its site plus the change of charge of its complexation on each plane, as PHREEQC's -cd_music Δz0 Δz1 Δz2 gives the latter. Fractions are allowed, and so are dual numbers: a charge distribution being fitted carries its derivative.

A published surface model is a list of reactions and constants, written the way PHREEQC (Parkhurst and Appelo, 2013) writes them. site_family turns such a list into a family, including the reactions read_sorption_model reads from a database.

ChemistryLab.site_family Method

julia

site_family(name, reactions, aqueous; master, site, capacity, support,
            model = IdealSiteMixing(), T = 298.15u"K", P = 1.0e5u"Pa")
    -> SiteFamily

A site family built from its surface reactions, written as PHREEQC writes them, and their constants.

reactions holds equation => log_K pairs, or the SorptionReactions read_sorption_model returns. Each reaction forms one surface complex from the free site and aqueous species, with a coefficient of one on both. master is the prefix naming the surface species in the equations ("Hfo_w"), and it is replaced by site, one of the SITE_SYMBOLS ("Xw"), so that Hfo_wOZn+ becomes the species XwOZn+.

aqueous supplies the aqueous species the reactions name, as a vector of species or a dictionary by symbol. H2O in an equation is the solvent, looked up as H2O@ when the collection has no H2O.

The standard energies are built so that every reaction has its constant, with the free site at zero and the aqueous species at their own ΔₐG⁰:

The free site at zero is a gauge while the site budget is fixed, and it is the required value under SITES_FOLLOW_HOST, where the free sites are counted as part of the host; host_coupling_bias explains why.

A published surface constant carries no reaction enthalpy, and log K is then held at every temperature, as PHREEQC holds a constant given without one: the energy of each complex follows the aqueous species of its reaction (ASSUMED: ΔᵣH = 0). The keywords T and P are accepted for compatibility and do not change the answer.

Examples

julia

dat = read_sorption_model(datapath("phreeqc.dat"))
weak = [r for r in dat.surfaces["Hfo_w"].reactions if !haskey(r.stoichiometry, "Hfo_sOH")]
family = site_family("Hfo_w", weak, aqueous; master = "Hfo_w", site = "Xw",
                     capacity = TotalSiteAmount(2.0e-4u"mol"), support)

The surface potential as an unknown

A DiffuseLayer written as an activity coefficient is exact and not always reachable: the solver's inner loop recovers a mixing phase with lnγ read at the previous iterate, which contracts only while electrostatic_stiffness stays below ELECTROSTATIC_STIFFNESS_LIMIT. Carrying the potential as an unknown of the outer Newton removes that fixed point rather than taming it.

ChemistryLab._compose_blocks Method

julia

_compose_blocks(a, b) -> NamedTuple

One parameter block from two: the constraint's unknowns first, the system's after. Each side keeps its own indices into q, so neither has to know the other exists.

The one combination that is refused

Both sides may declare hq, the callback that rebuilds the activity model at the current unknowns — an adiabatic solve does, because the Debye-Hückel coefficients are functions of temperature, and a diffuse layer does, because it is told its potential. Composing two of them would mean threading one's parameter merge through the other's, and the shipped hq closures are opaque to that.

Rather than compose them wrongly, this refuses and says so. A diffuse layer under an adiabatic constraint is not a case anyone has yet; when someone does, the fix is to have each block contribute the params it merges instead of a whole hq, and this error message is where they will start.

ChemistryLab._potential_supports Method

julia

_potential_supports(cs) -> Vector{Vector{Int}}

The distinct supports of cs that need a potential unknown, each as the list of family indices on it.

One unknown per surface, not per family: ferrihydrite's strong and weak sites are two families on one oxide and share one Ψ. See support_group.

Refuses a support whose families disagree — one asking for a potential and another not, or two declaring different areas for the same surface. Both are incoherent rather than merely unusual, and both would otherwise produce a number.

ChemistryLab._scatter Function

julia

_scatter(cs, supports, q, slots = eachindex(supports)) -> Vector

Each surface's potential placed at every family standing on it, nothing elsewhere, so the activity kernel can index by family without knowing which families share a surface: a number for a single plane, a triplet for a surface of three planes.

ChemistryLab._surface_potential_blocks Method

julia

_surface_potential_blocks(des, state, p, n0) -> NamedTuple or nothing

The parameter block the system contributes, as opposed to the one its constraint does: one unknown ψ̃ = FΨ/RT per surface that needs it, with the Gouy-Chapman closure as its equation, or three, one per plane, for a ChargePlanes surface closed by its two capacitors and its diffuse layer.

nothing when no surface needs one, which is every system without a diffuse layer — so nothing pays for the possibility.

The closure, and why its residual is already dimensionless

σ_k is summed over every member of every family on the k-th surface. ψ̃ is in units of RT by construction, which is what the stationarity rows are in, so no scaling is needed — unlike the enthalpy residual of an adiabatic solve, which is 10⁵ J and has to be divided by RT before it can sit in the same Newton system.

q0 = 0 starts the solve on an uncharged surface, which is the composition a solve without electrostatics would return and therefore the natural cold start.

ChemistryLab.needs_potential_unknown Method

julia

needs_potential_unknown(model) -> Bool

Whether a site mixing model's electrostatic term must be carried as an unknown of the solve rather than evaluated from the composition.

true for DiffuseLayer and ChargePlanes, one unknown per plane, and false for everything else, including ConstantCapacitance — whose potential is linear in the composition, so the fixed point contracts and eliminating it costs nothing.

The ions of the diffuse layer

A DiffuseLayer leaves the ions that screen the surface in the solution, which then carries the counter-charge. A DonnanLayer counts them in a layer of water on the surface, as PHREEQC's SURFACE -Donnan does, and equilibrate_donnan solves with them withdrawn from the solution.

ChemistryLab.DonnanLayer Type

julia

DonnanLayer(; thickness = 1.0e-8u"m")

The ions of the diffuse layer that balances a charged surface, averaged over a layer of water of fixed thickness: the Donnan approach that (Appelo and Wersin, 2007) added to PHREEQC as SURFACE -Donnan, which equilibrate_donnan solves with and diffuse_layer_contents reads.

On each support carrying a DiffuseLayer family, of area 𝒜, the layer holds W_D = ρ 𝒜 t of water (ρ = 1000 kg/m³, t the thickness) and each solute i of the solution at

with m_i the molality of the free solution and ψ̃_D = FΨ_D/RT the average potential of the layer, the one whose charge balances the surface charge σ. It is not the surface potential: the surface complexes still see the potential Gouy-Chapman gives, and ψ̃_D is the smaller average over the layer. The left-hand side decreases in ψ̃_D, so the root is unique.

The solutes of the layer come out of the system's totals, and the free solution holds the rest. Where its water comes from is a convention, and equilibrate_donnan offers PHREEQC's two: added to the solution's, which PHREEQC does when a solution and a surface are reacted together, or taken from it, which PHREEQC does for a surface equilibrated beforehand and which is the one a closed pore solution obeys. The excess of a solute, n_i^D − m_i W_D, is what the layer holds beyond its water at the molality of the free solution: positive for a counter-ion, negative for a co-ion.

thickness accepts a length or a number of meters. PHREEQC's default is 10 nm; a layer thicker than the pores it lines is not a physical statement, and the fixed point of equilibrate_donnan stops converging well before it.

ChemistryLab.diffuse_layer_contents Method

julia

diffuse_layer_contents(state, layer::DonnanLayer) -> NamedTuple

What the Donnan layers of state hold, one per support carrying a DiffuseLayer family (see DonnanLayer):

amounts: moles of each species of the system in the layers, zero but for the aqueous solutes and the solvent, whose entry is the layers' water;
excess: the same, less what the layers' water holds at the molality of the free solution — the ions the layers add to the solution's totals;
water: the layers' water, kg;
potential: ψ̃_D = FΨ_D/RT of each layer, by support name;
charge: the surface charge each layer balances, moles of charge.

state is a solution with its surfaces: the molalities are those of its free water, so the layers' water is not in the solvent of state.

ChemistryLab.equilibrate_donnan Method

julia

equilibrate_donnan(state, layer::DonnanLayer; model = DiluteSolutionModel(),
                   water = :added, b = nothing, maxiter = 200, rtol = 1e-8,
                   kwargs...)
    -> NamedTuple

The equilibrium of state with the ions of its diffuse layers counted, as PHREEQC's SURFACE -Donnan counts them (see DonnanLayer).

The layers take their solutes out of the solution's totals, and what they take depends on the solution they leave. Each step solves the solution and its surfaces with equilibrate_certified on the totals less what the layers held, from the previous answer, until the layers' contents stop moving by more than rtol relative. A layer that holds r times what the free solution holds of a solute would make a plain substitution oscillate with a factor r, and diverge beyond one; each step therefore moves that solute only the fraction 1/(1 + r) of the way, which cancels the factor. b is the system's total budget, the layers included; by default the one state holds.

water says where the layers' water comes from. :added, the default, adds it to the solution's: PHREEQC's convention for a solution and a surface reacted together, and the one its -Donnan results are compared against. :taken takes it out of the solution's water, so that the free water is the total less the layers': the convention of a closed system, and the one to use when the layers hold a noticeable share of the water, which :added would otherwise add to the system and dilute it by. Keywords other than these are passed to equilibrate_certified.

Returns (state, certificate, layer, iterations): the free solution with its surfaces and solids, the certificate of its last solve, and diffuse_layer_contents of that state. An iteration that does not settle within maxiter steps is an error that gives the last relative change. The more of the water the layers hold, the slower they settle: a layer holding 70 % of it takes about a hundred steps.

The site budget, and the declaration it has to match

A family declares a capacity; a state carries amounts. Until these functions existed nothing compared the two, so a capacity could be changed with no effect on the answer and a certificate could report optimal = true over a budget the declaration contradicted.

ChemistryLab.REFERENCE_SITE_DENSITY Constant

julia

REFERENCE_SITE_DENSITY

REFERENCE_SITE_DENSITY_NM2 in mol/m², derived through AVOGADRO rather than written again — 2.0009e-5 mol/m².

That it comes out so close to a round 2 × 10⁻⁵ mol/m² is not a coincidence: that is the number the convention was chosen to be, and the value in nm⁻² is how it reads in the units the surface literature uses.

ChemistryLab.REFERENCE_SITE_DENSITY_NM2 Constant

julia

REFERENCE_SITE_DENSITY_NM2

The conventional reference site density Γ° = 12.05 nm⁻² of (Kulik, 2002), used to define the standard state of a monodentate surface species, read from data/literature/Kulik2002.json.

What it is for

An intrinsic adsorption constant is not a property of a surface alone: it is fitted at some total site density Γ_C, and the value depends on that choice. Dzombak and Morel (1990) fitted their hydrous-ferric-oxide constants at two site densities, one for the weak sites and one for the strong, a factor of forty apart, so two of their own constants are not directly comparable with each other, let alone with a constant from another compilation.

Fixing one conventional Γ° for every sorbent and every surface makes them comparable, and convert_logk_site_density is the conversion. This is a choice of concentration scale, not a measurement: it is one number, it is stated here rather than assumed, and it can be changed.

ChemistryLab.InconsistentSiteBudget Type

julia

InconsistentSiteBudget <: Exception

A state whose site amounts contradict the capacity its family declares.

Carries the family, the residual in mol, and the scale it was judged against, so the message says how far off it is rather than only that it is off.

ChemistryLab._constraint_matrix Method

julia

_constraint_matrix(cs::ChemicalSystem) -> AbstractMatrix

The matrix every solver path constrains the equilibrium with: cs.SM.A itself when no family follows its host, and conservation_matrix otherwise. Returning SM.A untouched keeps a system without a coupled family bit-identical to what it was, and routing every path through here keeps a coupled one from being solved on the uncoupled matrix by one of them.

ChemistryLab._is_bare_site Method

julia

_is_bare_site(sp, site::Symbol) -> Bool

Whether sp carries the site pseudo-element site and nothing else — no real atom, possibly a charge.

Such a species is a component rather than a substance: the site with no matter attached. A coupled family may take one as the primary of its site row, or take its free site; conservation_matrix holds in either basis. The test matters for _refuse_unidentifiable_site, which concerns the bare basis only.

ChemistryLab._refuse_biased_coupling Method

julia

_refuse_biased_coupling(cs, family)

Refuse a coupled family whose free site is not at zero energy; see host_coupling_bias for why a coupled family needs it and a fixed one does not.

ChemistryLab._refuse_charged_free_site Method

julia

_refuse_charged_free_site(family, free)

Refuse a coupled family whose free site carries a charge. The coupling moves ν free sites with each mole of host, so a charged one would carry charge in and out of the system as the host dissolves or grows.

ChemistryLab._refuse_host_short_of_site_matter Method

julia

_refuse_host_short_of_site_matter(family, host, free, ν)

Refuse a coupled family whose host formula does not contain ν free sites' worth of every element: the free sites are counted as part of the host, so their atoms have to be among the host's own.

ChemistryLab._refuse_unidentifiable_site Method

julia

_refuse_unidentifiable_site(cs, family, A, r)

Refuse a coupled family whose site potential the basis cannot determine.

What goes wrong, measured

Every species carrying a site symbol carries it together with a fixed amount of charge: XsOH decomposes as H₂O − H⁺ + Xs + Zz and XsOCa⁺ as Ca²⁺ + H₂O − 2H⁺ + Xs + Zz. Where that is the only way Zz enters, the site row and the charge row of the coupled matrix are the same row up to the host entry, so only y_Xs + y_Zz is identifiable.

The solver finds that out the hard way. Measured on portlandite carrying sites, with a neutral bare component: the two multipliers ran to +234 650 and −234 710 — four decades past any chemical potential — while their sum stayed at −60; the dual Newton stalled at a KKT error of 2.3e-9; and the host came out at 0.0999 mol where the uncoupled answer is 0.0883.

Giving the component the charge the free site carries with its site symbol removes the Zz primary altogether. The same run then converges at 4.8e-11, max|y| is 227, the host lands on 0.08826 against 0.08827 uncoupled, and the coupling holds to 1.8e-7.

Why this is measured and not inferred from Zz being a primary

Refusing on the presence of a charge component is what this did first, and it refuses systems that are perfectly well posed. Zz survives as a primary whenever charge is genuinely independent of the element rows — which is the ordinary situation in a redox system. Measured, on hydrous ferric oxide over Fe(OH)₃(am) with chloride and a mixed-valence iron speciation, singular values of the coupled matrix:

system	component	`σ` (smallest three)	identifiable rank
portlandite, no `Zz` row possible	`Xs+`	`1.19, 0.872`	4 of 4
portlandite	`Xs`	`1.21, 0.957, 1.7e-5`	4 of 5
HFO, Fe(III) only	`Xs`	`2.16, 1.39, 0.105`	4 of 5
HFO, Fe(II) and Fe(III), chloride	`Xs+`	`1.90, 1.35, 0.916`	6 of 6
HFO, Fe(II) and Fe(III), chloride	`Xs`	`2.14, 1.42, 0.656`	6 of 6

The last two carry a Zz primary and are not degenerate at all — the charge row is nonzero on an iron chloride complex, so it is not the site row. The old test refused both of them, and the charge it then suggested was itself refused on the next call: the two suggestions pointed at each other and the user went in a circle.

identifiable_rank on the singular values separates the two groups with the spectral gap it is built for — a factor of 13 at the tightest refusal against 2.2 at the tightest acceptance. The charge suggested on a refusal is the coefficient of Zz in the free site's own decomposition, which is correct exactly where the refusal now fires.

ChemistryLab._standard_gibbs Method

julia

_standard_gibbs(sp) -> Union{Float64, Nothing}

ΔₐG⁰ at 298.15 K and 1 bar, in J/mol, or nothing when the species carries none. A missing property reads back as an Int zero rather than as an error (species.jl), so the test is on the type and not on the value: a species whose energy is genuinely zero must not be confused with one that has none.

ChemistryLab.check_site_budget Method

julia

check_site_budget(state::ChemicalState; rtol = 1e-6)

Throw InconsistentSiteBudget when a family's state amounts and its declared capacity disagree by more than rtol of the declaration.

The tolerance is relative to the declaration, not absolute: a site budget is 1e-6 mol on an oxide and 1e-1 mol on a clay, so an absolute threshold would be meaningless on one of them. A declaration of exactly zero is compared absolutely against the same rtol, since nothing else is available.

Where 1e-6 comes from

From the floor, not from a round number. A state is normally initialized with the occupied sites at the solver's 1e-12 floor rather than at zero, so a family with two complexes and a 1e-3 mol capacity is correctly initialized and still 2e-9 off in relative terms — measured. A default of 1e-9 rejects that, which is a check that fires on the ordinary case.

This exists to catch a declaration error, and those are orders of magnitude: the case it was written for declared 1e-6 mol against a state carrying 1e-3, a factor of a thousand. 1e-6 catches that with nine decades to spare and leaves the floor alone.

ChemistryLab.conservation_matrix Method

julia

conservation_matrix(cs::ChemicalSystem) -> Matrix

The matrix the equilibrium is constrained with: SM.A when no family follows its host, and otherwise SM.A with, for each coupled family, ν times the column of its free site subtracted from the column of its host,

The host's database formula already contains the surface groups its sites are made of, since a hydroxide carries the hydroxyls its surface exposes, so the free sites are counted as part of the host rather than on top of it. On the site row the subtraction reads Σ dₖ nₖ − ν n_host = 0, the statement of site_coupling_rows; on the element rows it removes from the host the atoms the free sites carry, so that every atom is counted once; and, the free site being neutral, it leaves the charge conserved while the host dissolves or grows. This is the bookkeeping of Kulik (2002) for a sorbent whose surface groups belong to it, and it holds in any basis of primaries: the free site may be the primary of the site row, or a bare site component (Species("Xs+")) may be declared in its place.

Refused, each by name: a free site that carries a charge, a host whose formula does not contain ν free sites' worth of every element, a capacity that is not proportional to the host's amount (sites_per_host), and a free site whose energy is away from zero (host_coupling_bias).

ChemistryLab.convert_logk_site_density Method

julia

convert_logk_site_density(logK, Γ_C; Γ0 = REFERENCE_SITE_DENSITY_NM2,
                          free_site_side = :reactant) -> Real

An intrinsic adsorption constant moved from the total site density it was fitted at to a reference one — equation (21) of Kulik (2002):

with the sign of the last term inverted when the neutral functional group sits on the product side of the reaction, which is what free_site_side selects. Written ≡OH + H⁺ = ≡OH₂⁺ the free site is a reactant, the default.

Only the ratio of the two densities enters, so they may be given in any unit as long as it is the same one. The default Γ0 is in nm⁻², which is how the surface literature quotes a site density.

Example

The two densities of Dzombak and Morel (1990), the case Kulik (2002) uses to make the point that their weak and strong constants are not on one scale:

julia

julia> using ChemistryLab

julia> weak = literature_value("Kulik2002", "dzombak_morel_weak_site_density_nm2");

julia> strong = literature_value("Kulik2002", "dzombak_morel_strong_site_density_nm2");

julia> round(convert_logk_site_density(0.0, weak), digits = 2)
-0.73

julia> round(convert_logk_site_density(0.0, strong), digits = 2)
-2.33

A constant fitted at a density below the reference moves down, and by different amounts for the two site types — 1.6 log units apart here. Correlating one against the other without this conversion compares two different scales.

ChemistryLab.declared_site_moles Method

julia

declared_site_moles(state::ChemicalState, family::SiteFamily) -> Real

The moles of sites family declares, evaluated on state — that is, site_moles supplied with the host amount and molar mass the state carries.

A TotalSiteAmount ignores all three and returns its number. The other two do not, so a family that measures its capacity per unit mass or per unit area needs a host, and one declared without naming a species is refused here rather than quietly evaluated at zero — a capacity of zero is a different model, not a missing input.

ChemistryLab.host_consistent_state Method

julia

host_consistent_state(state::ChemicalState) -> ChemicalState

state with each family's free-site amount derived from its declared capacity instead of taken from the caller:

so site_budget_residual comes back zero and the conservation row the solver forms says what the family declared.

The free site is the one adjusted, and deliberately: it is the family's reference member, the species whose amount carries no chemistry of its own. An occupied site holds a sorbate whose element budget the rest of the system accounts for, so moving it would silently move that element too.

Refused when the complexes already hold more sites than the capacity allows — that is not an initialization to repair but a declaration to correct, and scaling the complexes down to fit would invent a composition nobody asked for.

Everything outside the site families is copied unchanged.

ChemistryLab.host_coupling_bias Method

julia

host_coupling_bias(cs::ChemicalSystem) -> OrderedDict{String, <:Real}

Per coupled family, the shift in log SI that the reference energy of the free site imposes on the host, in the units the answer is read in:

Why the free site of a coupled family has zero energy

With a fixed budget the free site's ΔₐG⁰ cancels out of every surface reaction, both sides carrying one site, and it is a gauge. With a budget that follows its host it no longer cancels. conservation_matrix counts the free sites as part of the host, so that an intact sorbent, every site free, has the composition of the host's database formula; it has the database energy only when the free site's energy is zero. That is the convention of (Kulik, 2002) for surface groups that belong to their sorbent, and any other value moves the host's solubility by the amount above. A family whose bias exceeds 0.05 log units is refused, and site_family builds families that satisfy it, its free site at zero and every complex written relative to it.

The energies are read at 298.15 K, where the surface constants are given.

ChemistryLab.present_site_moles Method

julia

present_site_moles(state::ChemicalState, family::SiteFamily) -> Real

The moles of sites the state actually carries for family: Σ dᵢ nᵢ over its members, with dᵢ the denticity read from each member's formula.

This is exactly the left-hand side of the site conservation row, so comparing it with declared_site_moles compares the declaration against what the solver will conserve.

ChemistryLab.site_budget_residual Method

julia

site_budget_residual(state::ChemicalState) -> OrderedDict{String, <:Real}

Per declared family, present − declared in mol.

Zero means the state's site amounts and the family's capacity say the same thing. Anything else means the calculation will conserve the amounts and ignore the declaration, which is what it did unconditionally before this existed.

Returns an empty dictionary for a system with no surface, so it is safe to call on anything.

ChemistryLab.site_coupling_rows Method

julia

site_coupling_rows(cs::ChemicalSystem) -> (Matrix, Vector{String})

One row per family whose support is SITES_FOLLOW_HOST, stating that the sites in use equal ν times the host's amount:

with d_k the denticity of each member. Returns the rows and the family names that label them; both are empty when nothing is coupled, so a system without a coupled surface gets back exactly what it had.

These rows state the constraint so that it can be checked, and conservation_matrix is how it is imposed: on the site row the two coincide. Appending them to SM.A instead would not work, since SM.A already carries a site row and two equations on one sum forbid the host to move.

ChemistryLab.sites_per_host Method

julia

sites_per_host(family::SiteFamily, M_host) -> Real

ν, the moles of sites one mole of the host carries, for a family whose support is SITES_FOLLOW_HOST.

This is Kulik (2002), Eq. 30 read as a coefficient: the moles of sites of a surface type are ψ · A_s · M · Γ times the moles of sorbent, so the whole dependence on the host is one number multiplying its amount.

The capacity has to be homogeneous of degree one, and that is measured

ν only exists if the capacity really is proportional to the host's amount, and which capacities are is not a list to maintain but a property to check. This evaluates site_moles at two scaled amounts and requires the answer to scale with them:

with n₀ held fixed. Scaling n₀ along with n would make a ShrinkingCoreArea's ratio n/n₀ equal one at every λ and hide exactly the nonlinearity this exists to find.

What the probe accepts and refuses follows without being enumerated: MassSiteDensity and AreaSiteDensity over a specific area pass, being q M n and Γ a M n; a TotalSiteAmount and an AreaSiteDensity over a FixedSurfaceArea are refused, being constants that do not follow anything; and a ShrinkingCoreArea passes only at p = 1, where it genuinely is linear. A capacity type written later is judged on the same evidence rather than on whether somebody remembered to add it here.

A capacity of zero is refused too: it passes the proportionality test for the empty reason that 0 = λ·0, and it is a family with no sites, which is a declaration to correct.

Surfaces API ​

Areas ​

Size classes of a powder ​

Site families ​

The surface potential as an unknown ​

The ions of the diffuse layer ​

The site budget, and the declaration it has to match ​

Surfaces API

Areas

Size classes of a powder

Site families

The surface potential as an unknown

The ions of the diffuse layer

The site budget, and the declaration it has to match