ChemicalSystem and ChemicalState

These two types are the bridge between species definitions and equilibrium calculations:

ChemicalSystem — immutable container that groups species, stoichiometric matrices, and index maps. Built once; shared across many states.
ChemicalState — mutable snapshot of a system at a given temperature, pressure, and molar composition. Built from a ChemicalSystem; mutated in place.

ChemicalSystem type

What it holds

ChemicalSystem
├── species           → AbstractVector{<:AbstractSpecies}    (ordered list)
├── CSM               → CanonicalStoichMatrix                (elements × species)
├── SM                → StoichMatrix                         (primaries × species)
├── reactions         → Vector{<:AbstractReaction}           (derived or provided)
├── idx_aqueous       → Vector{Int}  (aqueous species)
├── idx_crystal       → Vector{Int}  (crystalline species)
├── idx_gas           → Vector{Int}  (gas-phase species)
├── idx_solvent       → Vector{Int}  (solvent, i.e. H₂O@)
├── idx_solutes       → Vector{Int}  (SC_AQSOLUTE)
├── idx_components    → Vector{Int}  (SC_COMPONENT)
└── solid_solutions   → Nothing | Vector{SolidSolutionPhase}

All derived fields (stoichiometric matrices, index maps) are computed once at construction and remain consistent for the lifetime of the object.

Minimal construction

using ChemistryLab

# Three aqueous species — no primaries specified → all species are primaries
H2O  = Species("H2O";  aggregate_state = AS_AQUEOUS, class = SC_AQSOLVENT)
Hp   = Species("H+";   aggregate_state = AS_AQUEOUS, class = SC_AQSOLUTE)
OHm  = Species("OH-";  aggregate_state = AS_AQUEOUS, class = SC_AQSOLUTE)

cs = ChemicalSystem([H2O, Hp, OHm])
length(cs)

cs.idx_solvent      # index of the solvent

1-element Vector{Int64}:
 1

cs.idx_solutes      # indices of solutes

2-element Vector{Int64}:
 2
 3

pprint(cs.CSM; label = :symbol)

┌────┬─────┬────┬─────┐
│    │ H2O │ H+ │ OH- │
├────┼─────┼────┼─────┤
│  H │   2 │  1 │   1 │
│  O │   1 │    │   1 │
│ Zz │     │  1 │  -1 │
└────┴─────┴────┴─────┘

Specifying primary species

Primary species (independent components) determine which species are "dependent" and how the stoichiometric matrix SM is built. Balanced reactions are then extracted from the null space of SM.

using ChemistryLab

H2O  = Species("H2O";  aggregate_state = AS_AQUEOUS, class = SC_AQSOLVENT)
Hp   = Species("H+";   aggregate_state = AS_AQUEOUS, class = SC_AQSOLUTE)
OHm  = Species("OH-";  aggregate_state = AS_AQUEOUS, class = SC_AQSOLUTE)
CO2  = Species("CO2";  aggregate_state = AS_AQUEOUS, class = SC_AQSOLUTE)
HCO3 = Species("HCO3-"; aggregate_state = AS_AQUEOUS, class = SC_AQSOLUTE)

# Primaries: H₂O, H⁺, HCO₃⁻  (a chemist's choice for the carbonate system)
primaries = [H2O, Hp, HCO3]

cs = ChemicalSystem([H2O, Hp, OHm, CO2, HCO3], primaries)
pprint(cs.SM; label = :symbol)

┌───────┬─────┬────┬─────┬─────┬───────┐
│       │ H2O │ H+ │ OH- │ CO2 │ HCO3- │
├───────┼─────┼────┼─────┼─────┼───────┤
│   H2O │   1 │    │   1 │  -1 │       │
│    H+ │     │  1 │  -1 │   1 │       │
│ HCO3- │     │    │     │   1 │     1 │
└───────┴─────┴────┴─────┴─────┴───────┘

# Independent reactions derived from the null space
rxns = reactions(cs.SM)
for r in rxns
    println(r.equation)
end

H₂O = OH⁻ + H⁺
H⁺ + HCO₃⁻ = CO₂ + H₂O

Filtered views

ChemicalSystem is an AbstractVector{<:AbstractSpecies}. Filtered views return sub-vectors without copying data:

using ChemistryLab

H2O   = Species("H2O";   aggregate_state = AS_AQUEOUS, class = SC_AQSOLVENT)
Hp    = Species("H+";    aggregate_state = AS_AQUEOUS, class = SC_AQSOLUTE)
Cal   = Species("Cal";   aggregate_state = AS_CRYSTAL, class = SC_COMPONENT)
CO2g  = Species("CO2";   aggregate_state = AS_GAS,     class = SC_GASFLUID)

cs = ChemicalSystem([H2O, Hp, Cal, CO2g])

println("aqueous:  ", symbol.(aqueous(cs)))
println("crystal:  ", symbol.(crystal(cs)))
println("gas:      ", symbol.(gas(cs)))
println("solvent:  ", symbol(solvent(cs)))
println("solutes:  ", symbol.(solutes(cs)))

aqueous:  ["H2O", "H+"]
crystal:  ["Cal"]
gas:      ["CO2"]
solvent:  H2O
solutes:  ["H+"]

Species lookup

using ChemistryLab

H2O  = Species("H2O"; aggregate_state = AS_AQUEOUS, class = SC_AQSOLVENT)
Hp   = Species("H+";  aggregate_state = AS_AQUEOUS, class = SC_AQSOLUTE)
cs   = ChemicalSystem([H2O, Hp])

cs["H2O"]       # lookup by symbol string

Species{Int64}
           name: H2O
         symbol: H2O
        formula: H2O ◆ H₂O
          atoms: H => 2, O => 1
         charge: 0
aggregate_state: AS_AQUEOUS
          class: SC_AQSOLVENT
     properties: M = 0.0180149999937744 kg mol⁻¹

cs[2]            # lookup by index

Species{Int64}
           name: H+
         symbol: H+
        formula: H+ ◆ H⁺
          atoms: H => 1
         charge: 1
aggregate_state: AS_AQUEOUS
          class: SC_AQSOLUTE
     properties: M = 0.001007999999651657 kg mol⁻¹

# Find the index of a species
findfirst(s -> symbol(s) == "H+", cs.species)

Getting reactions

When the system was built from a database (with thermodynamic data attached to species), each dependent species has a corresponding dissolution/formation reaction:

# From a database workflow
cs = ChemicalSystem(species, primaries)
r_cal = get_reaction(cs, "Cal")   # reaction for calcite
r_cal.logK⁰(T = 298.15)

When water (the aqueous solvent, SC_AQSOLVENT) is added to a ChemicalState that contains H⁺ and OH⁻ species, ChemistryLab automatically seeds them at the neutral pH concentration $c = 10^{-pK_w/2}$ based on the water autoprotolysis constant $K_w(T, P)$. This ensures a chemically consistent initial state without manual intervention.

The auto-seeding only triggers when both H⁺ and OH⁻ are at zero — if you set them explicitly (e.g. to impose a specific initial pH), your values are preserved.

v = volume(state)
println("V liquid = ", v.liquid)
println("V solid  = ", v.solid)
println("V total  = ", v.total)

V liquid = 0.0 m³
V solid  = 0.0 m³
V total  = 0.0 m³

m = moles(state)
println("n liquid = ", m.liquid)
println("n solid  = ", m.solid)

n liquid = 55.509297826565565 mol
n solid  = 0.0 mol

porosity(state)

NaN

Copying a state

copy creates a new ChemicalState that shares the underlying ChemicalSystem (no duplication) but has its own independent mole vector, temperature, and pressure:

using ChemistryLab
using DynamicQuantities

H2O = Species("H2O"; aggregate_state = AS_AQUEOUS, class = SC_AQSOLVENT)
Hp  = Species("H+";  aggregate_state = AS_AQUEOUS, class = SC_AQSOLUTE)
cs  = ChemicalSystem([H2O, Hp])

s1 = ChemicalState(cs)
set_quantity!(s1, "H2O", 1.0u"kg")

s2 = copy(s1)
set_temperature!(s2, 350.0u"K")   # only s2 is modified

ustrip(s1.T[]), ustrip(s2.T[])

(298.15, 350.0)

Full workflow example

Building a minimal carbonate system from scratch — no database required:

using ChemistryLab
using DynamicQuantities

# 1. Declare species
H2O  = Species("H2O";   aggregate_state = AS_AQUEOUS, class = SC_AQSOLVENT)
Hp   = Species("H+";    aggregate_state = AS_AQUEOUS, class = SC_AQSOLUTE)
OHm  = Species("OH-";   aggregate_state = AS_AQUEOUS, class = SC_AQSOLUTE)
CO2  = Species("CO2";   aggregate_state = AS_AQUEOUS, class = SC_AQSOLUTE)
HCO3 = Species("HCO3-"; aggregate_state = AS_AQUEOUS, class = SC_AQSOLUTE)
CO3  = Species("CO3-2";  aggregate_state = AS_AQUEOUS, class = SC_AQSOLUTE)

# 2. Build the system (H₂O, H⁺, CO₃²⁻ as primaries)
cs = ChemicalSystem([H2O, Hp, OHm, CO2, HCO3, CO3], [H2O, Hp, CO3])

println("Species: ", join(symbol.(cs.species), ", "))
println("Aqueous: ", join(symbol.(aqueous(cs)), ", "))

Species: H2O, H+, OH-, CO2, HCO3-, CO3-2
Aqueous: H2O, H+, OH-, CO2, HCO3-, CO3-2

# 3. Build an initial state
state = ChemicalState(cs)
set_quantity!(state, "H2O",   1.0u"kg")
set_quantity!(state, "H+",    1e-7u"mol/L" * volume(state).liquid)
set_quantity!(state, "OH-",   1e-7u"mol/L" * volume(state).liquid)
set_quantity!(state, "CO2",   1e-3u"mol")

println("pH       = ", pH(state))
println("n liquid = ", moles(state).liquid)

pH       = nothing
n liquid = 55.51029782656556 mol

Next step: equilibrium

Once you have a ChemicalSystem and a ChemicalState, pass the state to equilibrate to find the thermodynamic equilibrium. See the Equilibrium page for details.