Skip to content

Thermochemistry: the quantities and how they connect

This page is the shortest path from "what is a chemical potential" to the expressions this package actually evaluates, with the notation of the code rather than of a textbook. Everything here is either a definition or an identity; nothing is fitted.

1. One potential, and the two pieces it is made of

The Gibbs energy of a system holding nᵢ moles of each species is  , and the chemical potential of a species is what one more mole of it costs:

Two pieces, and they come from two different places in the package:

  • , the standard potential, is thermodynamic data — it belongs to the species alone and comes from the database (§3);

  • , the activity, is where the mixture enters, and it comes from the activity model.

Divide by and the code appears verbatim. build_potentials returns the closure

julia
μ(n, p) = p.ΔₐG⁰overRT .+ lna(n, p)      #  μᵢ/RT  =  ΔₐG⁰ᵢ/RT  +  ln aᵢ

so throughout this documentation

Working in is not cosmetic: it is dimensionless, it is what the minimizer differentiates, and it makes and the standard potential directly comparable in size.

2. Activities, standard states, and the conventions in use

An activity is always a ratio to a standard state, and the choice of standard state is a convention that must be stated or the numbers are meaningless. The three in use here:

classactivitystandard stateconvention
aqueous solute    mol/kg of solventmolality
aqueous solute (ideal model)   mol/Lmolarity
solvent (water)pure watermole fraction
pure solid, pure phase  the pure substance
gas in an ideal mixture the pure gas at mole fraction
solid-solution end-member the pure end-membermole fraction (§5)

concentration_scale is how a model declares which of the first two it uses, and the aqueous accessors read it rather than guessing:

julia
using ChemistryLab
concentration_scale(DiluteSolutionModel()), concentration_scale(HKFActivityModel())
(:molarity, :molality)

An activity coefficient is therefore not a property of a species but of a species in a solution, and   in the standard state's own limit — infinite dilution for a solute, purity for a solid.

3. Where comes from

The database stores, per species, the coefficients of a heat-capacity polynomial, and the package integrates it. That is the whole of the temperature dependence, and it is exact rather than a series truncation: THERMO_MODELS (src/thermodynamics/thermo_models.jl) holds and the analytic integrals of the same coefficients for , and , so the three identities

hold by construction rather than numerically. Aqueous solutes use a second model, :solute_hkf88_reaktoro, the Helgeson-Kirkham-Flowers equation of state, whose electrostatic part is the Born term evaluated with the dielectric constant of water — the same water model that supplies and to the activity models.

Two consequences worth carrying:

  • a reaction's temperature dependence is not an input; it follows from the of its reactants and products, which is the Gibbs-Helmholtz relation    (van 't Hoff, when is taken constant);

  • a species whose database entry has no model cannot be moved off at all, and the package says so rather than extrapolating.

4. Equilibrium is a constrained minimization, and its dual is the useful part

At fixed and , equilibrium is the composition that minimizes subject to the matter available. "The matter available" is a linear constraint, because atoms and charge are conserved whatever the reactions do:

where is the conservation matrix cs.SM.A — one row per component, one column per species — and the element (and charge) budget. This is the formulation of (Leal et al., 2017), and it is why no reaction list is needed: reactions are whatever moves inside the null space of .

The Lagrangian multipliers of the equality constraints are the interesting output. Writing for the multiplier of row , the first-order conditions are

The are component potentials (element potentials, when the components are elements): one number per component, from which the chemical potential of any species is a dot product. The two lines above are the complementarity conditions of the KKT system, and they are exactly what optimality_certificate checks — see Proving that an answer is the answer for the certificate itself and why convexity makes it sufficient rather than merely necessary.

5. Reaction quotients, , and the saturation index

Take a species and write its formation from the system's primary species,    — the column of is that reaction. For any reaction,

so  : the sign of says which way the reaction runs, and   is equilibrium. For a solid dissolving into its ions, is the ion activity product IAP and is the solubility product, and the logarithm of their ratio is the saturation index.

Because the rows of are labeled by primary species, a component potential is that primary's , and the index becomes a difference of potentials with no equilibrium constant to look up:

Negative means undersaturated, zero means in equilibrium with the solution, positive means the phase should have precipitated. This is what saturation_indices returns.

Two consequences of the identity are worth carrying, because they are what make the index trustworthy rather than merely conventional:

  • no equilibrium constant is looked up. never appears in the computation; it is implied by the standard potentials, so an index cannot disagree with the the rest of the calculation used;

  • every phase present at an equilibrium must come out at exactly zero. That is not a property of the answer, it is the first-order condition of §4 restated, so it is a check on the solve: measured on a CEM I paste, the twelve present solids land within  . If they do not, the state is not an equilibrium and no other index in the result means anything.

The identity itself — that the difference of component potentials above equals what saturation_indices computes — is asserted in the test suite (test/aqueous_properties.jl), where it is recomputed by hand from this formula and compared to the function.

6. Volume, porosity and chemical shrinkage

Volumes are treated as ideal: the volume of a phase is the sum of its species' standard molar volumes,

with no excess volume of mixing. That is an assumption, and it is the one behind every porosity this package reports. volume returns the split by aggregate state, porosity the void fraction relative to a reference state, and chemical_shrinkage the volume the reaction itself consumes — hydrates occupying less than the water and clinker they were made from, which is the mechanism behind self-desiccation (Self-desiccation).

A species carrying no contributes zero and would corrupt a porosity in silence, which is why the volume machinery reports missing_molar_volumes(state) and why CapillaryWater refuses to be constructed when that list is not empty.

7. What is not assumed, and what is

The two lists, side by side:

assumednot assumed
one well-mixed phase per aggregate statea reaction list, a reaction path, a sequence
ideal molar volumes, no excess volumeideal activities — that is the model's business
the activity model's own domain of validitythat the answer is a local minimum: the certificate proves global (see below)
that   is the whole of the conservationwhich phases appear: they are a result

Convexity is what makes the last cell true. is convex in for the built-in activity models, so a KKT point is the minimum and the certificate is sufficient rather than necessary — the one documented exception being CapillaryWater, where a composition-dependent shift of breaks the convexity argument and the certificate falls back to proving a KKT point.

See also

certificate