ChemistryLab.jl is a computational chemistry toolkit. Although initially dedicated to low-carbon cementitious materials and aqueous solutions and designed for researchers, engineers, and developers working with cement chemistry, its scope is actually wider. It provides formula handling, species management, stoichiometric matrix construction, and database interoperability (ThermoFun and Cemdata). Main features include chemical formula parsing, Unicode/Phreeqc notation conversion, reaction and equilibrium analysis, and data import/export.
- Chemical formula handling: Create, convert, and display formulas with charge management and Unicode/Phreeqc notation.
- Chemical species management:
SpeciesandCemSpeciestypes to represent solution and solid phase species;with_classto requalify a species without modifying the original. - Stoichiometric matrices: Automatic construction of matrices for reaction and equilibrium analysis.
- Database interoperability: Import and merge ThermoFun (.json) and Cemdata (.dat) data; load solid solution definitions from a TOML file with
build_solid_solutions. - Parsing tools: Convert chemical notations, extract charges, calculate molar mass, and more.
- Solid solutions: Define ideal (
IdealSolidSolutionModel) or non-ideal binary (RedlichKisterModel) mineral mixing phases viaSolidSolutionPhase; end-members are automatically requalified at construction time. - Activity models: Built-in aqueous activity models for equilibrium:
DiluteSolutionModel(ideal),HKFActivityModel(extended Debye-Hückel B-dot), andDaviesActivityModel. - Chemical equilibrium: Compute thermodynamic equilibrium compositions from initial states using Gibbs energy minimization (
equilibrate,ChemicalSystem,ChemicalState).
The package can be installed with the Julia package manager.
From the Julia REPL, type ] to enter the Pkg REPL mode and run:
pkg> add ChemistryLabOr, equivalently, via the Pkg API:
julia> import Pkg; Pkg.add("ChemistryLab")Let's imagine we want to study the equilibrium of calcite in water.
To do this, we can create a list of chemical species, retrieve the thermodynamic properties of these species from one of the databases integrated into ChemistryLab. We can then deduce the chemical species likely to appear in the reaction and calculate the associated stoichiometric matrix.
In this example, the database is cemdata. The .json file is included in ChemistryLab but is a copy of a file which can be found in ThermoHub.
using ChemistryLab
using DynamicQuantities
# From the repository root (or adapt the path for an installed package):
# all_species = build_species(joinpath(pkgdir(ChemistryLab), "data", "cemdata18-thermofun.json"))
all_species = build_species("data/cemdata18-thermofun.json")The chemical species likely to appear during calcite equilibrium in water are obtained in the following way:
species_calcite = speciation(all_species, split("Cal H2O@ CO2");
aggregate_state=[AS_AQUEOUS],
exclude_species=split("H2@ O2@ CH4@"))
dict_species_calcite = Dict(symbol(s) => s for s in species_calcite)The output of dict_species_calcite reads:
Dict{String, Species} with 12 entries:
"H+" => H+ {H+} [H+ ◆ H⁺]
"OH-" => OH- {OH-} [OH- ◆ OH⁻]
"CO2" => CO2 {CO2 g} [CO2 ◆ CO₂]
"Ca(HCO3)+" => Ca(HCO3)+ {CaHCO3+} [Ca(HCO3)+ ◆ Ca(HCO₃)⁺]
"Cal" => Cal {Calcite} [CaCO3 ◆ CaCO₃]
"CaOH+" => CaOH+ {CaOH+} [Ca(OH)+ ◆ Ca(OH)⁺]
"H2O@" => H2O@ {H2O l} [H2O@ ◆ H₂O@]
"Ca+2" => Ca+2 {Ca+2} [Ca+2 ◆ Ca²⁺]
"CO2@" => CO2@ {CO2 aq} [CO2@ ◆ CO₂@]
"HCO3-" => HCO3- {HCO3-} [HCO3- ◆ HCO₃⁻]
"CO3-2" => CO3-2 {CO3-2} [CO3-2 ◆ CO₃²⁻]
"Ca(CO3)@" => Ca(CO3)@ {CaCO3 aq} [CaCO3@ ◆ CaCO₃@]
During species creation, ChemistryLab calculates the molar mass of the species. It also constructs thermodynamic functions (heat capacity, entropy, enthalpy, and Gibbs free energy of formation) as a function of temperature.
dict_species_calcite["Cal"]Species{Int64}
name: Calcite
symbol: Cal
formula: CaCO3 ◆ CaCO₃
atoms: Ca => 1, C => 1, O => 3
charge: 0
aggregate_state: AS_CRYSTAL
class: SC_COMPONENT
properties: M = 0.10008599996541243 kg mol⁻¹
Tref = 298.15 K
Pref = 100000.0 m⁻¹ kg s⁻²
Cp⁰ = 104.5163192749 + 0.02192415855825T + -2.59408e6 / (T^2) [m² kg s⁻² K⁻¹ mol⁻¹] ◆ T=298.15 K
ΔₐH⁰ = -1.2482415842895252e6 + 104.5163192749T + 2.59408e6 / T + 0.010962079279125(T^2) [m² kg s⁻² mol⁻¹] ◆ T=298.15 K
S⁰ = -523.9438829693111 + 0.02192415855825T + 104.5163192749log(T) + 1.29704e6 / (T^2) [m² kg s⁻² K⁻¹ mol⁻¹] ◆ T=298.15 K
ΔₐG⁰ = -1.1423813547027335e6 + 628.460202244211T + 1.29704e6 / T - 0.010962079279125(T^2) - 104.5163192749T*log(T) [m² kg s⁻² mol⁻¹] ◆ T=298.15 K
V⁰ = 3.6933999061584004e-5 [m³ mol⁻¹]
Cp⁰_Tref = 81.87109375 m² kg s⁻² K⁻¹ mol⁻¹
ΔₐH⁰_Tref = -1.207405e6 m² kg s⁻² mol⁻¹
S⁰_Tref = 92.675598144531 m² kg s⁻² K⁻¹ mol⁻¹
ΔₐG⁰_Tref = -1.129176e6 m² kg s⁻² mol⁻¹
V⁰_Tref = 3.6933999061584004e-5 m³ mol⁻¹
The evolution of thermodynamic properties as a function of temperature, such as heat capacity, can thus be easily plotted.
using Plots
p1 = plot(xlabel="Temperature [K]", ylabel="Cp⁰ [J/mol/K]", title="Heat capacity of calcite \nas a function of temperature")
plot!(p1, θ -> dict_species_calcite["Cal"].Cp⁰(T = 273.15+θ), 0:0.1:100, label="Cp⁰")
Obtaining stoichiometric matrices requires the choice of a species-independent basis.
primaries = [dict_species_calcite[s] for s in split("H2O@ H+ CO3-2 Ca+2")]
SM = StoichMatrix(values(dict_species_calcite), primaries)┌───────┬────┬─────┬─────┬───────────┬─────┬───────┬──────┬──────┬──────┬───────┬───────┬──────────┐
│ │ H+ │ OH- │ CO2 │ Ca(HCO3)+ │ Cal │ CaOH+ │ H2O@ │ Ca+2 │ CO2@ │ HCO3- │ CO3-2 │ Ca(CO3)@ │
├───────┼────┼─────┼─────┼───────────┼─────┼───────┼──────┼──────┼──────┼───────┼───────┼──────────┤
│ H2O@ │ │ 1 │ -1 │ │ │ 1 │ 1 │ │ -1 │ │ │ │
│ H+ │ 1 │ -1 │ 2 │ 1 │ │ -1 │ │ │ 2 │ 1 │ │ │
│ CO3-2 │ │ │ 1 │ 1 │ 1 │ │ │ │ 1 │ 1 │ 1 │ 1 │
│ Ca+2 │ │ │ │ 1 │ 1 │ 1 │ │ 1 │ │ │ │ 1 │
└───────┴────┴─────┴─────┴───────────┴─────┴───────┴──────┴──────┴──────┴───────┴───────┴──────────┘
These stoichiometric matrices thus allow us to write the chemical reactions at work.
list_reactions = reactions(SM)
dict_reactions_calcite = Dict(r.symbol => r for r in list_reactions)Again, when constructing the reactions, the thermodynamic properties of the reactions as a function of temperature are deduced. It is thus possible to see, for example, the expression for the solubility product of calcite for the reaction under study and to plot its evolution.
dict_reactions_calcite["Cal"].logK⁰SymbolicFunc:
Expression: (-1.29704e6 + 125345.63212888106T - 2666.9195440882527(T^2) + 424.77184295654104(T^2)*log(T) + 0.010962079279125T*(T^2)) / (19.144757680815896(T^2)) [m² kg s⁻² mol⁻¹]
References: T=298.15 K
Variables: T
p1 = plot(xlabel="Temperature [K]", ylabel="pKs", title="Solubility product (pKs) of calcite \nas a function of temperature")
plot!(p1, θ -> dict_reactions_calcite["Cal"].logK⁰(T = 273.15+θ), 0:0.1:100, label="pKs")
Once the chemical system is defined, the equilibrium composition can be computed directly from an initial state. A ChemicalSystem groups all species and pre-computes the conservation matrix; a ChemicalState holds the mole amounts, temperature, and pressure.
# Build the chemical system from the species returned by speciation
cs = ChemicalSystem(collect(values(dict_species_calcite)))
# Define an initial state: dissolve 1e-3 mol of calcite in 1 kg of water (≈ 55.5 mol) at 25 °C
state0 = ChemicalState(cs; T = 298.15u"K", P = 1u"bar")
set_quantity!(state0, "H2O@", 55.5u"mol")
set_quantity!(state0, "Cal", 1e-3u"mol")
set_quantity!(state0, "H+", 1e-7u"mol")
set_quantity!(state0, "OH-", 1e-7u"mol")
# Compute thermodynamic equilibrium (Gibbs energy minimization)
state_eq = equilibrate(state0)Key results can then be inspected directly:
pH(state_eq) # equilibrium pH
moles(state_eq) # mole amounts by phase (liquid / solid / gas / total)
moles(state_eq, "Ca+2") # moles of a specific speciesThe equilibrate function uses IpoptOptimizer under the hood (via Optimization.jl) and accepts optional keyword arguments to tune the solver (abstol, reltol, variable_space, …), as well as an model keyword for the aqueous activity model:
state_eq = equilibrate(state0; model = HKFActivityModel()) # extended Debye-Hückel
state_eq = equilibrate(state0; model = DaviesActivityModel()) # Davies equationMineral solid solutions (e.g. C-S-H, AFm) are modelled by grouping end-member species into a SolidSolutionPhase. Database species with SC_COMPONENT class can be passed directly — requalification to SC_SSENDMEMBER is handled automatically:
substances = build_species("data/cemdata18-thermofun.json")
dict = Dict(symbol(s) => s for s in substances)
# Ideal solid solution (any number of end-members)
cshq = SolidSolutionPhase("CSHQ", [dict["CSHQ-TobD"], dict["CSHQ-TobH"],
dict["CSHQ-JenH"], dict["CSHQ-JenD"]])
# Non-ideal binary: Redlich-Kister (parameters in J/mol)
afm = SolidSolutionPhase("AFm", [dict["Ms"], dict["Mc"]];
model = RedlichKisterModel(a0 = 3000.0, a1 = 500.0))
# Or load all solid solution phases at once from a TOML file
ss_phases = build_solid_solutions("data/solid_solutions.toml", dict)
cs = ChemicalSystem(species_list, primaries; solid_solutions = ss_phases)A pre-built data/solid_solutions.toml (CSHQ, AFm, Hydrogarnet, Ettringite_ss, Hydrotalcite) is shipped with ChemistryLab for use with the cemdata18 database.
A ChemicalState supports scalar multiplication and in-place rescaling to a target total:
state2 = state_eq * 2.0 # double all amounts (non-mutating)
state_h = state_eq / 1000 # millimolar scale (non-mutating)
rescale!(state_eq, 1.0u"mol") # total moles → 1 mol (in-place)
rescale!(state_eq, 1.0u"kg") # total mass → 1 kg (in-place)
rescale!(state_eq, 1.0u"m^3") # total volume → 1 m³ (in-place)See the documentation and tutorials for examples on formula creation, species management, reaction parsing, and database merging.
MIT License. See LICENSE for details.
See CITATION.cff for citation details.
BibTeX entry:
@software{chemistrylab_jl,
authors = {Barthélémy, Jean-François and Soive, Anthony},
title = {ChemistryLab.jl: Numerical laboratory for computational chemistry},
doi = {10.5281/zenodo.17756074},
url = {https://github.com/ChemistryTools/ChemistryLab.jl}
}Developed by Jean-François Barthélémy and Anthony Soive, both researchers at Cerema in the research team UMR MCD.