Chemical Reactions

In ChemistryLab it is possible to build chemical reactions and manipulate them. A reaction is constructed as a structure, "a composite data type that allows you to store multiple values in a single object". The struct is organized as follows:

julia

struct Reaction{SR<:AbstractSpecies, TR<:Number, SP<:AbstractSpecies, TP<:Number}
    equation::String
    colored::String
    reactants::OrderedDict{SR, TR}
    products::OrderedDict{SP, TP}
    equal_sign::Char
    properties::OrderedDict{Symbol,PropertyType}
end

Parsing reactions

Reaction is a composite type struct which can be build from:

a string containing species

julia

equation = "13H⁺ + NO₃⁻ + CO₃²⁻ + 10e⁻ = 6H₂O@ + HCN@"
reac, prod, equal_sign = parse_equation(equation)

a string containing cement species

julia

eqC3S = "C₃S + 5.3H = 1.3CH + C₁.₇SH₄"
rC3S = CemReaction(eqC3S)

an operation on species

julia

using ChemistryLab
C3S = CemSpecies("C3S", symbol="C₃S", aggregate_state=AS_CRYSTAL, class=SC_COMPONENT)
H = CemSpecies("H", symbol="H₂O@", aggregate_state=AS_AQUEOUS, class=SC_AQSOLVENT)
CH = CemSpecies("CH", symbol="C₃S", aggregate_state=AS_CRYSTAL, class=SC_COMPONENT)
CSH = CemSpecies("C1.7SH4", aggregate_state=AS_CRYSTAL, class=SC_COMPONENT)
r = C3S + 5.3H ↔ 1.3CH + CSH
typeof(r)

Reaction{CemSpecies{Int64, Int64}, Float64, CemSpecies, Float64, Int64}

a balance calculation

julia

using ChemistryLab
C3S = CemSpecies("C3S")
H = CemSpecies("H")
CH = CemSpecies("CH")
CSH = CemSpecies("C1.7SH4")
r = Reaction([C3S, H, CH, CSH]; equal_sign='→')
pprint(r)

  equation: C₃S + 5.3H → 1.3CH + C₁.₇SH₄
 reactants: C₃S => 1, H => 5.3
  products: CH => 1.3, C₁.₇SH₄ => 1
    charge: 0

a balance calculation with symbolic numbers

julia

using ChemistryLab
using Symbolics
@variables a b g
CSH = CemSpecies(Dict(:C => a, :S => one(Num), :H => g))
C3S = CemSpecies("C3S")
H = CemSpecies("H")
CH = CemSpecies("CH")
r = Reaction([CSH, C3S, H, CH]; equal_sign='→')

  equation: CₐSHg → (3-a+g)H + C₃S + (-3+a)CH
 reactants: CₐSHg => 1
  products: H => 3 - a + g, C₃S => 1, CH => -3 + a
    charge: 0

julia

using ChemistryLab
using Symbolics
using PrettyTables
@variables a b g
CSH = CemSpecies(Dict(:C => a, :S => one(Num), :H => g))
C3S = CemSpecies("C3S")
H = CemSpecies("H")
CH = CemSpecies("CH")
# `apply`, not `map`: `map` over a reaction iterates its coefficients and hands
# back a plain vector, losing the reaction. `apply` rebuilds one.
r = apply(simplify, Reaction([C3S, H], [CH, CSH]; equal_sign = '→'))
r.equation

"C₃S + (3-a+g)H → (3-a)CH + CₐSHg"

The same decomposition as a matrix, with the symbols carried into it:

julia

using ChemistryLab
using Symbolics
using PrettyTables
@variables a b g
CSH = CemSpecies(Dict(:C => a, :S => one(Num), :H => g))
C3S = CemSpecies("C3S")
H = CemSpecies("H")
CH = CemSpecies("CH")
SM = StoichMatrix([C3S], [CSH, H, CH])
pprint(SM)

┌───────┬────────────┐
│       │        C3S │
├───────┼────────────┤
│     H │ -3 + a - g │
│ CₐSHg │          1 │
│    CH │      3 - a │
└───────┴────────────┘

Collection of arrow symbols

In ChemistryLab, there are collections of arrow symbols used in chemical reaction notation, such as:

>, →, ↣, ↦, ⇾, ⟶, ⟼, ⥟, ⥟, ⇀, ⇁, ⇒, ⟾ for reaction directionality from reactants to products;
<, ←, ↢, ↤, ⇽, ⟵, ⟻, ⥚, ⥞, ↼, ↽, ⇐, ⟽ for reaction directionality from products to reactants;
↔, ⟷, ⇄, ⇆, ⇌, ⇋, ⇔, ⟺ for reversible reactions and equilibrium states;
=, ≔, ⩴, ≕ to separate reactants from products in balanced equations.

Thermodynamic properties of reactions

When the species involved in a reaction carry thermodynamic data (loaded from a database), the reaction automatically exposes temperature-dependent thermodynamic functions. These are computed lazily on first access and stored in the reaction's properties dict:

Property	Description
`r.ΔᵣCp⁰`	Heat capacity of reaction (J mol⁻¹ K⁻¹)
`r.ΔᵣH⁰`	Enthalpy of reaction (J mol⁻¹)
`r.ΔᵣS⁰`	Entropy of reaction (J mol⁻¹ K⁻¹)
`r.ΔᵣG⁰`	Gibbs free energy of reaction (J mol⁻¹)
`r.logK⁰`	Decimal logarithm of the equilibrium constant

Each property is a SymbolicFunc callable with a keyword argument T (temperature in K):

julia

using ChemistryLab

# Load reactions from a database-built stoichiometric matrix
all_species = build_species(datapath("cemdata18-merged.json"); verbose = false)
species = speciation(all_species, split("Cal H2O@");
              aggregate_state=[AS_AQUEOUS], exclude_species=split("H2@ O2@ CH4@"))
dict_species = Dict(symbol(s) => s for s in species)
candidate_primaries = [dict_species[s] for s in CEMDATA_PRIMARIES if haskey(dict_species, s)]
cs = ChemicalSystem(species, candidate_primaries)
list_reactions = reactions(cs.SM)
dict_reactions = Dict(r.symbol => r for r in list_reactions)

r_cal = dict_reactions["Cal"]   # calcite dissolution reaction

# Evaluate at 25 °C
r_cal.logK⁰(T = 298.15)          # log₁₀ K at 25 °C
r_cal.ΔᵣG⁰(T = 298.15)           # ΔᵣG° at 25 °C  (J/mol)

Properties can also be set manually on a reaction:

julia

r = Reaction("H2 + O2 = H2O")
r[:ΔᵣH⁰] = -241800.0   # J/mol

Reactions from a stoichiometric matrix

The most common source of reactions in a database workflow is reactions(SM), which derives all independent reactions from a StoichMatrix:

julia

list_reactions = reactions(cs.SM)
dict_reactions = Dict(r.symbol => r for r in list_reactions)

Each reaction symbol matches the dependent species it describes. Temperature-dependent thermodynamic functions are computed automatically when the constituent species carry thermodynamic data.

Algebraic reaction balancing

The null space of the stoichiometric matrix automatically yields balanced chemical reactions — no manual coefficient guessing needed. This is the algebraic counterpart to the chemical insight that mass must be conserved.

Principle

Given a set of species, ChemicalSystem builds the stoichiometric matrix A (rows = elements, columns = species). Any vector ν in the null space of A satisfies Aν = 0, i.e. it represents a balanced reaction. reactions(cs.SM) extracts one independent balanced reaction per dependent species.

Example: combustion of methane

julia

using ChemistryLab

# Declare gas-phase species
CH4 = Species("CH4"; aggregate_state = AS_GAS, class = SC_GASFLUID)
O2  = Species("O2";  aggregate_state = AS_GAS, class = SC_GASFLUID)
CO2 = Species("CO2"; aggregate_state = AS_GAS, class = SC_GASFLUID)
H2O = Species("H2O"; aggregate_state = AS_GAS, class = SC_GASFLUID)

# CH4, O2, H2O are the primary (independent) components; CO2 is the derived species
primaries = [CH4, O2, H2O]
cs = ChemicalSystem([CH4, O2, CO2, H2O], primaries)

rxns = reactions(cs.SM)
pprint(rxns[1])

    symbol: CO2
  equation: CH₄ + 2O₂ = CO₂ + 2H₂O
 reactants: CH₄ => 1, O₂ => 2
  products: CO₂ => 1, H₂O => 2
    charge: 0

Reading the matrix behind it

The balance above is read off a matrix, and the matrix displays itself — there is no reason to transcribe one by hand into a table, and a transcription is one more thing that can be wrong. cs.CSM is the canonical form, atoms down the rows and species across the columns:

julia

pprint(cs.CSM; label = :symbol)

┌───┬─────┬────┬─────┬─────┐
│   │ CH4 │ O2 │ CO2 │ H2O │
├───┼─────┼────┼─────┼─────┤
│ C │   1 │    │   1 │     │
│ H │   4 │    │     │   2 │
│ O │     │  2 │   2 │   1 │
└───┴─────┴────┴─────┴─────┘

cs.SM is the other one, and the pair is worth separating. The canonical matrix above counts atoms. SM expresses every species over the primaries chosen for the system, which is what the balance actually uses:

julia

pprint(cs.SM; label = :symbol)

┌─────┬─────┬────┬─────┬─────┐
│     │ CH4 │ O2 │ CO2 │ H2O │
├─────┼─────┼────┼─────┼─────┤
│ H2O │     │    │  -2 │   1 │
│ CH4 │   1 │    │   1 │     │
│  O2 │     │  1 │   2 │     │
└─────┴─────┴────┴─────┴─────┘

With CH₄, O₂ and H₂O as primaries, CO₂ is the dependent species, and its column in SM is the reaction: CO₂ = CH₄ + 2 O₂ − 2 H₂O, that is CH₄ + 2 O₂ → CO₂ + 2 H₂O.

Example: combustion of alkanes (CₙH₂ₙ₊₂)

julia

using ChemistryLab

C2H6 = Species("C2H6"; aggregate_state = AS_GAS, class = SC_GASFLUID)
C3H8 = Species("C3H8"; aggregate_state = AS_GAS, class = SC_GASFLUID)
O2   = Species("O2";   aggregate_state = AS_GAS, class = SC_GASFLUID)
CO2  = Species("CO2";  aggregate_state = AS_GAS, class = SC_GASFLUID)
H2O  = Species("H2O";  aggregate_state = AS_GAS, class = SC_GASFLUID)

primaries = [CO2, H2O, O2]

# Ethane: C₂H₆ + 7/2 O₂ → 2 CO₂ + 3 H₂O
# reactions(SM) returns the reverse here; negate to get the forward combustion direction
cs_ethane = ChemicalSystem([C2H6, O2, CO2, H2O], primaries)
pprint(-reactions(cs_ethane.SM)[1])

    symbol: C2H6
  equation: C₂H₆ + 7//2O₂ = 3H₂O + 2CO₂
 reactants: C₂H₆ => 1, O₂ => 7//2
  products: H₂O => 3, CO₂ => 2
    charge: 0

julia

# Propane: C₃H₈ + 5 O₂ → 3 CO₂ + 4 H₂O
cs_propane = ChemicalSystem([C3H8, O2, CO2, H2O], primaries)
pprint(-reactions(cs_propane.SM)[1])

    symbol: C3H8
  equation: C₃H₈ + 5O₂ = 4H₂O + 3CO₂
 reactants: C₃H₈ => 1, O₂ => 5
  products: H₂O => 4, CO₂ => 3
    charge: 0

The general alkane, balanced once

The two blocks above each balance one alkane. The family can be balanced instead, by giving the species a symbolic carbon number and letting the same element balance carry it through — n sits in the formula exactly where an integer sits, and every operation downstream keeps it:

julia

using ChemistryLab
using Symbolics

@variables n
CnH = Species(Dict(:C => n, :H => 2n + 2))
O2, H2O, CO2 = Species.(split("O2 H2O CO2"))

r = apply(Symbolics.expand, Reaction([CnH, O2], [H2O, CO2]))
r.equation

"CₙH₂₊₂ₙ + ((1//2)+(3//2)n)O₂ = ((1//1)+n)H₂O + nCO₂"

That is the general law, derived rather than asserted: one half plus three halves of n moles of oxygen, n of carbon dioxide and n + 1 of water. The two numeric cases above are two of its values, and so is every other alkane:

julia

for vn in 1:6
    rn = apply(x -> Symbolics.value(substitute(x, Dict(n => vn))), r)
    println("  n = ", vn, " : ", rn.equation)
end

  n = 1 : CH₄ + 2//1O₂ = 2//1H₂O + CO₂
  n = 2 : C₂H₆ + 7//2O₂ = 3//1H₂O + 2CO₂
  n = 3 : C₃H₈ + 5//1O₂ = 4//1H₂O + 3CO₂
  n = 4 : C₄H₁₀ + 13//2O₂ = 5//1H₂O + 4CO₂
  n = 5 : C₅H₁₂ + 8//1O₂ = 6//1H₂O + 5CO₂
  n = 6 : C₆H₁₄ + 19//2O₂ = 7//1H₂O + 6CO₂

Do not convert the coefficients while substituting

Writing apply(x -> Int(substitute(x, d)), r) looks tidier and is wrong: half the alkanes need a half-integer of oxygen, Int(7//2) throws, and apply returns the original coefficient when the function throws rather than raising. The output then reads as a mixture of numbers and symbols that looks like a partial substitution instead of like a failure. Keep the rational.

Example: combustion of alkenes (CₙH₂ₙ)

julia

using ChemistryLab

C2H4 = Species("C2H4"; aggregate_state = AS_GAS, class = SC_GASFLUID)
C3H6 = Species("C3H6"; aggregate_state = AS_GAS, class = SC_GASFLUID)
O2   = Species("O2";   aggregate_state = AS_GAS, class = SC_GASFLUID)
CO2  = Species("CO2";  aggregate_state = AS_GAS, class = SC_GASFLUID)
H2O  = Species("H2O";  aggregate_state = AS_GAS, class = SC_GASFLUID)

primaries = [CO2, H2O, O2]

# Ethylene: C₂H₄ + 3 O₂ → 2 CO₂ + 2 H₂O
# reactions(SM) returns the reverse here; negate to get the forward combustion direction
cs_ethylene = ChemicalSystem([C2H4, O2, CO2, H2O], primaries)
pprint(-reactions(cs_ethylene.SM)[1])

    symbol: C2H4
  equation: C₂H₄ + 3O₂ = 2H₂O + 2CO₂
 reactants: C₂H₄ => 1, O₂ => 3
  products: H₂O => 2, CO₂ => 2
    charge: 0

julia

# Propylene: C₃H₆ + 9/2 O₂ → 3 CO₂ + 3 H₂O
cs_propylene = ChemicalSystem([C3H6, O2, CO2, H2O], primaries)
pprint(-reactions(cs_propylene.SM)[1])

    symbol: C3H6
  equation: C₃H₆ + 9//2O₂ = 3H₂O + 3CO₂
 reactants: C₃H₆ => 1, O₂ => 9//2
  products: H₂O => 3, CO₂ => 3
    charge: 0

Alkenes CₙH₂ₙ follow: CₙH₂ₙ + 3n/2 O₂ → n CO₂ + n H₂O.

Example: multiple simultaneous reactions

When more than one dependent species is present, reactions(cs.SM) returns one balanced equation per dependent species:

julia

using ChemistryLab

CH4  = Species("CH4";  aggregate_state = AS_GAS, class = SC_GASFLUID)
C2H6 = Species("C2H6"; aggregate_state = AS_GAS, class = SC_GASFLUID)
C2H4 = Species("C2H4"; aggregate_state = AS_GAS, class = SC_GASFLUID)
O2   = Species("O2";   aggregate_state = AS_GAS, class = SC_GASFLUID)
CO2  = Species("CO2";  aggregate_state = AS_GAS, class = SC_GASFLUID)
H2O  = Species("H2O";  aggregate_state = AS_GAS, class = SC_GASFLUID)

primaries = [CO2, H2O, O2]
cs = ChemicalSystem([CH4, C2H6, C2H4, O2, CO2, H2O], primaries)

rxns = reactions(cs.SM)
for r in rxns
    println((-r).equation)
end

CH₄ + 2O₂ = 2H₂O + CO₂
C₂H₆ + 7//2O₂ = 3H₂O + 2CO₂
C₂H₄ + 3O₂ = 2H₂O + 2CO₂

When to use algebraic balancing

Algebraic balancing via reactions(cs.SM) is most useful when:

The system has many species and manual balancing would be error-prone.
You load species from a database and want all independent reactions automatically.
You want to verify that a set of species is stoichiometrically consistent.

Chemical Reactions ​

Parsing reactions ​

Thermodynamic properties of reactions ​

Reactions from a stoichiometric matrix ​

Algebraic reaction balancing ​

Principle ​

Example: combustion of methane ​

Reading the matrix behind it ​

Example: combustion of alkanes (CₙH₂ₙ₊₂) ​

The general alkane, balanced once ​

Example: combustion of alkenes (CₙH₂ₙ) ​

Example: multiple simultaneous reactions ​

Chemical Reactions

Parsing reactions

Thermodynamic properties of reactions

Reactions from a stoichiometric matrix

Algebraic reaction balancing

Principle

Example: combustion of methane

Reading the matrix behind it

Example: combustion of alkanes (CₙH₂ₙ₊₂)

The general alkane, balanced once

Example: combustion of alkenes (CₙH₂ₙ)

Example: multiple simultaneous reactions