Equilibrium

abstract type AbstractActivityModel end

Base type for all activity models. Each concrete subtype must implement activity_model(cs::ChemicalSystem, model::AbstractActivityModel) which returns a closure a(n, p) -> Vector{Float64} of log-activities.

struct DiluteSolutionModel <: AbstractActivityModel

Ideal dilute solution model:

• Solvent: Raoult's law — ln a = ln(x_solvent)
• Solutes: Henry's law — ln a = ln(c_i / c°) where c° = 1 mol/L
• Crystals: pure solid — ln a = 0
• Gas: ideal mixture — ln a = ln(x_i)

activity_model(cs::ChemicalSystem, ::DiluteSolutionModel) -> Function

Return a closure lna(n, p) -> Vector{Float64} computing the vector of log-activities for the dilute ideal solution model.

The returned function has signature lna(n, p) where:

• n: dimensionless mole vector (same indexing as cs.species)
• p: NamedTuple containing at least ϵ (floor value to avoid log(0))

All quantities are dimensionless — units are stripped at construction time.

build_potentials(cs::ChemicalSystem, model::AbstractActivityModel) -> Function

Return a closure μ(n, p) -> Vector{Float64} computing dimensionless chemical potentials μ_i / RT for all species.

$\mu_i / RT = \Delta_a G_i^0 / RT + \ln a_i$

The returned function is compatible with SciML solvers:

• n: dimensionless mole vector
• p: NamedTuple containing:
  • ΔₐG⁰overT: vector of standard Gibbs energies of formation divided by RT
  • ϵ: regularization floor (e.g. 1e-30)

All quantities are dimensionless — caller is responsible for stripping units from ΔₐG⁰overT before passing them in p.

Examples

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

julia> μ = build_potentials(cs, DiluteSolutionModel());

julia> n = [55.5, 0.1];

julia> p = (ΔₐG⁰overT = [-95.6, -105.6], ϵ = 1e-30);

julia> length(μ(n, p)) == 2
true

EquilibriumProblem

Definition of a chemical equilibrium problem.

Fields

• b: conservation vector (elemental abundances).
• A: stoichiometric matrix (conservation matrix).
• μ: chemical potential function μ(n, p).
• u0: initial guess for species amounts.
• p: coefficients for the potential function (default: NullParameters).
• lb: lower bounds for species amounts.
• ub: upper bounds for species amounts.

The problem solves for the species distribution that minimizes the Gibbs energy subject to mass conservation constraints A * n = b.

struct EquilibriumSolver{F<:Function, S, V<:Val}

Encapsulates all fixed ingredients of a chemical equilibrium calculation: the potential function, the SciML solver, and the variable space.

Construct once, call repeatedly with different ChemicalState inputs.

Fields

• μ: chemical potential closure μ(n, p) -> Vector{Float64}.
• solver: any Optimization.jl-compatible solver (e.g. IpoptOptimizer()).
• variable_space: variable space — Val(:linear) or Val(:log).
• kwargs: solver keyword arguments forwarded to solve.

Examples

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

julia> solver = EquilibriumSolver(cs, DiluteSolutionModel(), IpoptOptimizer());

julia> solver isa EquilibriumSolver
true

equilibrate(state::ChemicalState;
            model    = DiluteSolutionModel(),
            solver   = IpoptOptimizer(...),
            variable_space  = Val(:log),
            ϵ        = 1e-16,
            kwargs...) -> ChemicalState

Compute the chemical equilibrium state from an initial ChemicalState.

This is a high-level convenience function that wraps EquilibriumSolver and solve with sensible defaults. It is intended for users who do not need to fine-tune the solver.

Arguments

• state: initial ChemicalState — defines the system, T, P, and composition.
• model: activity model (default: DiluteSolutionModel()).
• solver: SciML-compatible solver (default: IpoptOptimizer with tight tolerances).
• variable_space: variable space — Val(:linear) or Val(:log) (default). Val(:log) is more robust for systems spanning many orders of magnitude.
• ϵ: regularization floor for mole amounts (default: 1e-16).
• kwargs...: additional keyword arguments forwarded to the solver.

Returns

A new ChemicalState at thermodynamic equilibrium, with all derived quantities (pH, pOH, volumes, porosity, saturation) already computed.

Examples

# Minimal usage — all defaults
state_eq = equilibrate(state)

# Custom temperature-dependent run
set_temperature!(state, 350.0u"K")
state_eq = equilibrate(state)

# Custom solver tolerance
state_eq = equilibrate(state; abstol=1e-12, reltol=1e-12)

# Custom activity model (future)
state_eq = equilibrate(state; model=DebyeHuckelModel(A=0.51, B=3.28))