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 aChemicalSystem; 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)3cs.idx_solvent # index of the solvent1-element Vector{Int64}:
1cs.idx_solutes # indices of solutes2-element Vector{Int64}:
2
3pprint(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)
endH₂O = OH⁻ + H⁺
H⁺ + HCO₃⁻ = CO₂ + H₂OThe components have to span the species
A component list is not a preference. Every species must be writable as a combination of the components, because that combination is the conservation law the equilibrium solver enforces for it. A list that cannot express a species is refused, by name:
try
ChemicalSystem([H2O, Hp, OHm], [H2O]) # water alone
catch e
println(sprint(showerror, e))
endArgumentError: these species cannot be written over the chosen components, so no conservation law covers them: H+ (H+ ◆ H⁺), OH- (OH- ◆ OH⁻). The components are [H2O]. Add a component carrying the missing element, or drop the species: decomposing it anyway projects it onto the components and lets the solver create it out of nothing.Water alone cannot describe an acid-base system: those three species span a two-dimensional space, and one component cannot reach it. The refusal matters more than it looks, because the decomposition is a least-squares projection and a projection never fails — asked for this, it used to answer H+ = 0.4 H2O and OH- = 0.6 H2O, which balances arithmetically and lets a solver make H+ out of water with no OH- and no charge to pay for it.
The check is a rank comparison in exact rational arithmetic, so it is the same on every machine and has no tolerance to tune. When it fires, add a component carrying the missing element — or let the constructor choose for you by naming no components at all, in which case every species is a candidate.
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 stringSpecies{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 indexSpecies{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)2Getting 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)ChemicalState type
Construction
A ChemicalState is created from a ChemicalSystem. All mole amounts start at zero; temperature defaults to 298.15 K and pressure to 10⁵ Pa.
using ChemistryLab
using DynamicQuantities
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], [H2O, Hp])
state = ChemicalState(cs)ChemicalState{Species{Int64}, AbstractReaction, DynamicQuantities.Quantity{Float64, DynamicQuantities.SymbolicDimensions{DynamicQuantities.FRInt32}}, Float64}
┌──────────────────────────────────────────────────────┐
│ T : 298.15 K │
│ P : 1.0 bar │
╞══════════════════════════════════════════════════════╡
│ # liquid #│ n [mol]│ m [g]│
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ tot. liquid│ 0.0│ 0.0│
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ H2O│ 0.0│ 0.0│
│ H+│ 0.0│ 0.0│
│ OH-│ 0.0│ 0.0│
╞══════════════════════════════════════════════════════╡
│ # TOTAL #│ n [mol]│ m [g]│
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ │ 0.0│ 0.0│
╞══════════════════════════════════════════════════════╡
└──────────────────────────────────────────────────────┘A state displays itself, and that display is the normal way to look at one: the species grouped by phase with their amounts, masses and volumes, a total per phase, and the scalar diagnostics — pH, pOH, porosity, saturation — underneath. Fresh out of the constructor everything is zero, which is what an empty state looks like:
stateChemicalState{Species{Int64}, AbstractReaction, DynamicQuantities.Quantity{Float64, DynamicQuantities.SymbolicDimensions{DynamicQuantities.FRInt32}}, Float64}
┌──────────────────────────────────────────────────────┐
│ T : 298.15 K │
│ P : 1.0 bar │
╞══════════════════════════════════════════════════════╡
│ # liquid #│ n [mol]│ m [g]│
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ tot. liquid│ 0.0│ 0.0│
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ H2O│ 0.0│ 0.0│
│ H+│ 0.0│ 0.0│
│ OH-│ 0.0│ 0.0│
╞══════════════════════════════════════════════════════╡
│ # TOTAL #│ n [mol]│ m [g]│
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ │ 0.0│ 0.0│
╞══════════════════════════════════════════════════════╡
└──────────────────────────────────────────────────────┘The temperature and pressure are read from the fields:
ustrip(state.T[]) # temperature in K298.15ustrip(state.P[]) # pressure in Pa100000.0Override temperature and pressure at construction:
state2 = ChemicalState(cs; T = 350.0u"K", P = 2e5u"Pa")
ustrip(state2.T[])350.0Setting molar amounts
set_quantity! accepts any compatible unit — moles, kilograms, grams, or concentration (mol/L) multiplied by a volume:
using ChemistryLab
using DynamicQuantities
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)
Cal = Species("Cal"; aggregate_state = AS_CRYSTAL, class = SC_COMPONENT)
# Three components, not two: calcite carries elements that water and the proton
# do not, so a two-component basis cannot express it — and a decomposition that
# cannot express a species is refused rather than projected onto the ones that
# are there. `OH-` needs no component of its own, being `H2O - H+`.
cs = ChemicalSystem([H2O, Hp, OHm, Cal], [H2O, Hp, Cal])
state = ChemicalState(cs)
# 1 litre of water (by mass)
set_quantity!(state, "H2O", 1.0u"kg")
# pH 4 water: 10⁻⁴ mol/L × volume of liquid phase
V = volume(state)
set_quantity!(state, "H+", 1e-4u"mol/L" * V.liquid)
set_quantity!(state, "OH-", 1e-10u"mol/L" * V.liquid)
# 1 mmol of calcite
set_quantity!(state, "Cal", 1e-3u"mol")
# Inspect moles
ustrip.(state.n)4-element Vector{Float64}:
55.509297826565565
0.0
0.0
0.001The state itself shows the same amounts in context — which phase each species belongs to, and what the phase totals are:
stateChemicalState{Species{Int64}, AbstractReaction, DynamicQuantities.Quantity{Float64, DynamicQuantities.SymbolicDimensions{DynamicQuantities.FRInt32}}, Float64}
┌──────────────────────────────────────────────────────┐
│ T : 298.15 K │
│ P : 1.0 bar │
╞══════════════════════════════════════════════════════╡
│ # liquid #│ n [mol]│ m [g]│
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ tot. liquid│ 55.5093│ 1000.0│
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ H2O│ 55.5093│ 1000.0│
│ H+│ 0.0│ 0.0│
│ OH-│ 0.0│ 0.0│
╞══════════════════════════════════════════════════════╡
│ # solid #│ n [mol]│ m [g]│
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ tot. solid│ 0.001│ 0.040078│
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ Cal│ 0.001│ 0.040078│
╞══════════════════════════════════════════════════════╡
│ # TOTAL #│ n [mol]│ m [g]│
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ │ 55.5103│ 1000.04│
╞══════════════════════════════════════════════════════╡
└──────────────────────────────────────────────────────┘The volumes come out at zero and the porosity at NaN here, and that is not a defect: these species were built from formulas, which carry a composition and a molar mass but no molar volume. Volume and porosity below takes the same state from a database instead, and the same display then fills in.
Changing temperature and pressure
set_temperature!(state, 350.0u"K")
set_pressure!(state, 2e5u"Pa")
ustrip(state.T[])350.0Derived quantities
All derived quantities are recomputed automatically after any set_*! call:
using ChemistryLab
using DynamicQuantities
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], [H2O, Hp])
state = ChemicalState(cs)
set_quantity!(state, "H2O", 1.0u"kg")
# H⁺ and OH⁻ are auto-seeded at neutral pH — no manual seeding needed
pH(state)pOH(state)Automatic H⁺/OH⁻ seeding
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
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.
Volume and porosity
A volume needs a molar volume, and a species built from a formula does not have one: Species("H2O") knows its composition and its molar mass, both computed from the formula, but nothing about how much room a mole of it takes. So the state above — real enough for pH, which needs only amounts — reports 0.0 m³ for every volume and NaN for the porosity, a quotient of two zeros.
Volumes therefore come from a database, where the molar volume is a measured quantity:
using ChemistryLab
using DynamicQuantities
substances = build_species(datapath("cemdata18-thermofun.json"); verbose = false)
sp = speciation(substances, ["Portlandite"]; aggregate_state = [AS_AQUEOUS])
cs = ChemicalSystem(sp, CEMDATA_PRIMARIES)
wet = ChemicalState(cs)
set_quantity!(wet, "H2O@", 1.0u"kg")
set_quantity!(wet, "Portlandite", 5.0u"mol")
wetChemicalState{Species, AbstractReaction, DynamicQuantities.Quantity{Float64, DynamicQuantities.SymbolicDimensions{DynamicQuantities.FRInt32}}, Float64}
┌────────────────────────────────────────────────────────────────────────────────────────────────┐
│ T : 298.15 K │
│ P : 1.0 bar │
╞════════════════════════════════════════════════════════════════════════════════════════════════╡
│ # liquid #│ n [mol]│ m [g]│ V [cm³]│ c [mol/L]│
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ tot. liquid│ 55.5093│ 1000.0│ 1002.96│ │
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ H2O@│ 55.5093│ 1000.0│ 1002.96│ 55.3453│
│ H+│ 1.00282e-7│ 1.01084e-7│ 0.0│ 9.99856e-8│
│ OH-│ 1.00282e-7│ 1.7055e-6│ -4.72111e-7│ 9.99856e-8│
│ Ca+2│ 0.0│ 0.0│ -0.0│ 0.0│
│ CaOH+│ 0.0│ 0.0│ 0.0│ 0.0│
│ H2@│ 0.0│ 0.0│ 0.0│ 0.0│
│ O2@│ 0.0│ 0.0│ 0.0│ 0.0│
╞════════════════════════════════════════════════════════════════════════════════════════════════╡
│ # solid #│ n [mol]│ m [g]│ V [cm³]│ │
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ tot. solid│ 5.0│ 370.46│ 165.3│ │
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ Portlandite│ 5.0│ 370.46│ 165.3│ │
╞════════════════════════════════════════════════════════════════════════════════════════════════╡
│ # TOTAL #│ n [mol]│ m [g]│ V [cm³]│ │
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ │ 60.5093│ 1370.46│ 1168.26│ │
╞════════════════════════════════════════════════════════════════════════════════════════════════╡
│ pH : 7.0001 │
│ pOH : 7.0001 │
│ porosity : 0.858508 │
│ saturation : 1.0 │
└────────────────────────────────────────────────────────────────────────────────────────────────┘This is the same display as above, on species that carry molar volumes: every column is now filled, the liquid and the solid are totaled separately, and the porosity at the bottom is a number rather than NaN.
A negative volume in that table is physics, not a defect
OH- shows a volume of about -4.7e-7 cm³ for 1.0e-7 mol, which is a partial molar volume of −4.7 cm³/mol. Partial molar volumes of ions are genuinely negative: the charge pulls the surrounding water in tighter than bulk water is packed, so adding the ion makes the solution smaller. The figure is CEMDATA18's, not an artifact of the arithmetic, and the phase totals are right to add it with its sign.
v = volume(wet)
println("V liquid = ", v.liquid)
println("V solid = ", v.solid)
println("V total = ", v.total)V liquid = 0.00100296403235498 m³
V solid = 0.0001652999997139 m³
V total = 0.0011682640320688801 m³moles splits the amounts the same way:
m = moles(wet)
println("n liquid = ", m.liquid)
println("n solid = ", m.solid)n liquid = 55.509298027129454 mol
n solid = 5.0 molporosity is then the liquid share of the total volume — for a suspension of five moles of portlandite in a kilogram of water, most of it:
porosity(wet)0.8585080125926927A volume of zero means missing data, not an empty phase
Every accessor that divides by a volume — porosity, and the molarity conventions of pH — returns NaN when the species carry no molar volume. NaN here is the honest answer to 0/0 and a sign that the species came from formulas rather than from a database, not a defect in the state.
Rescaling a state
A recipe is usually written for a chosen basis — one kilogram of paste, one mole of binder, one cubic meter of concrete — and rescale! multiplies every amount by the one factor that puts the state on it. The composition does not change; only the size of the sample does, so every intensive quantity (pH, porosity, the mole fractions) is left exactly where it was.
The target's dimension chooses what is held: an amount rescales the total moles, a mass the total mass, a volume the total volume.
rescale!(wet, 1.0u"kg") # the same paste, weighed out to one kilogram
# `.total`, not `sum`: `mass` returns `(liquid, solid, gas, total)`, so summing
# the NamedTuple adds the phases AND their total and reports twice the mass.
println("total mass = ", mass(wet).total)
println("porosity = ", porosity(wet), " (unchanged: it is intensive)")total mass = 0.9999999999999998 kg
porosity = 0.8585080125926927 (unchanged: it is intensive)A volume target needs molar volumes for the same reason as above, and the rescaling is refused rather than silently wrong if the current total is zero.
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 molstateChemicalState{Species{Int64}, AbstractReaction, DynamicQuantities.Quantity{Float64, DynamicQuantities.SymbolicDimensions{DynamicQuantities.FRInt32}}, Float64}
┌──────────────────────────────────────────────────────┐
│ T : 298.15 K │
│ P : 1.0 bar │
╞══════════════════════════════════════════════════════╡
│ # liquid #│ n [mol]│ m [g]│
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ tot. liquid│ 55.5103│ 1000.04│
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ H2O│ 55.5093│ 1000.0│
│ CO2│ 0.001│ 0.044009│
│ H+│ 0.0│ 0.0│
│ OH-│ 0.0│ 0.0│
│ HCO3-│ 0.0│ 0.0│
│ CO3-2│ 0.0│ 0.0│
╞══════════════════════════════════════════════════════╡
│ # TOTAL #│ n [mol]│ m [g]│
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ │ 55.5103│ 1000.04│
╞══════════════════════════════════════════════════════╡
└──────────────────────────────────────────────────────┘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.