Thermodynamic Functions

ChemistryLab represents temperature-dependent thermodynamic properties (Cp°, ΔₐH°, S°, ΔₐG°, log K°, …) as callable function objects rather than plain numbers. This page explains the three concrete types — SymbolicFunc, NumericFunc, and ThermoFactory — and shows how to build, combine, and extend them.

Overview

All three types share the same calling convention:

f(; T = 298.15, P = 1e5)          # raw numeric result (SI units)
f(; T = 350.0u"K", unit = true)   # result as a Quantity with units
f()                                # uses defaults stored in refs

Simple function objects

Constant

using ChemistryLab
using DynamicQuantities

# Numeric constant (unitless)
f0 = SymbolicFunc(42.0)
f0()

42.0

# Constant with physical unit
Cp_ref = SymbolicFunc(37.14u"J/mol/K")
Cp_ref()

37.14

Cp_ref(unit = true)    # returns a Quantity

37.14 m² kg s⁻² K⁻¹ mol⁻¹

Single variable

# Just T — useful as a building block for arithmetic
f_T = SymbolicFunc(:T)
f_T(T = 350.0)

350.0

Expression with parameters

Provide the expression as a Julia Expr literal. Any symbol that is not a declared variable (default list: [:T, :P, :t, :x, :y, :z]) is treated as a parameter and must be given as a keyword argument.

# Linear: Cp = a₀ + a₁·T
f_lin = SymbolicFunc(:(a₀ + a₁ * T); a₀ = 30.0, a₁ = 5e-3)
f_lin(T = 298.15)

31.49075

# Quadratic: Cp = a + b·T + c/T²
f_quad = SymbolicFunc(:(a + b * T + c / T^2); a = 25.0, b = 8e-3, c = -1.5e5)
f_quad(T = 500.0)

28.4

Calling conventions

# Pass T as a plain number (K)
f_lin(T = 400.0)

32.0

# Pass T as a Quantity — unit-aware evaluation
f_lin(T = 400.0u"K")

32.0

# Return result with units attached
f_lin(T = 400.0u"K", unit = true)

32.0 

Arithmetic on function objects

SymbolicFunc objects support all standard arithmetic operations. The result is always another SymbolicFunc that can be evaluated the same way.

using ChemistryLab
using DynamicQuantities

f1 = SymbolicFunc(:(a₀ + a₁ * T); a₀ = 30.0, a₁ = 2e-3)
f2 = SymbolicFunc(:(b₀ + b₁ * T); b₀ = 10.0, b₁ = 1e-3)

f_sum  = f1 + f2      # Cp of a reaction = Σ νᵢ Cpᵢ
f_diff = f1 - f2
f_scal = 2.5 * f1     # scaling by a stoichiometric coefficient

f_sum(T = 298.15)

40.89445

f_diff(T = 298.15)

20.29815

f_scal(T = 500.0)

77.5

Built-in thermodynamic models

ChemistryLab ships three built-in models, registered in THERMO_MODELS and pre-compiled in THERMO_FACTORIES.

:cp_ft_equation — Maier-Kelley polynomial (default)

The standard heat-capacity polynomial for crystalline and gas-phase species:

Cp°(T) = a₀ + a₁T + a₂/T² + a₃/√T + a₄T² + a₅T³ + a₆T⁴ + a₇/T³ + a₈/T + a₉√T + a₁₀log(T)

S°, ΔₐH°, and ΔₐG° are the analytical integrals, adjusted so that the reference-state values S°(Tᵣ), ΔₐH°(Tᵣ), ΔₐG°(Tᵣ) match the provided data.

Parameters (all keyword arguments, in SI):

Coefficient units:

Example — CO₂ (gas):

using ChemistryLab
using DynamicQuantities

params_CO2 = Dict(
    :S⁰   => 213.785u"J/K/mol",
    :ΔₐH⁰ => -393510.0u"J/mol",
    :ΔₐG⁰ => -394373.0u"J/mol",
    :a₀   => 33.98u"J/K/mol",
    :a₁   => 23.88e-3u"J/(mol*K^2)",
    :a₂   => 0.0u"J*K/mol",
    :a₃   => 0.0u"J/(mol*K^0.5)",
    :T    => 298.15u"K",
)

dtf = build_thermo_functions(:cp_ft_equation, params_CO2)

OrderedCollections.OrderedDict{Symbol, SymbolicFunc{1, @NamedTuple{T::DynamicQuantities.Quantity{Float64, DynamicQuantities.Dimensions{DynamicQuantities.FRInt32}}}, Float64, DynamicQuantities.Dimensions{DynamicQuantities.FRInt32}}} with 4 entries:
  :Cp⁰  => 33.98 + 0.02388T [m² kg s⁻² K⁻¹ mol⁻¹] ◆ vars=(T) ◆ T=298.15 K
  :ΔₐH⁰ => -4.04703e5 + 33.98T + 0.01194(T^2) [m² kg s⁻² mol⁻¹] ◆ vars=(T) ◆ T=298.15 K
  :S⁰   => 13.0608 + 0.02388T + 33.98log(T) [m² kg s⁻² K⁻¹ mol⁻¹] ◆ vars=(T) ◆ T=298.15 K
  :ΔₐG⁰ => -3.41826e5 + 20.9192T - 0.01194(T^2) - 33.98T*log(T) [m² kg s⁻² mol⁻¹] ◆ vars=(T) ◆ T=298.15 K

The returned OrderedDict has four entries:

keys(dtf)

KeySet for a OrderedCollections.OrderedDict{Symbol, SymbolicFunc{1, @NamedTuple{T::DynamicQuantities.Quantity{Float64, DynamicQuantities.Dimensions{DynamicQuantities.FRInt32}}}, Float64, DynamicQuantities.Dimensions{DynamicQuantities.FRInt32}}} with 4 entries. Keys:
  :Cp⁰
  :ΔₐH⁰
  :S⁰
  :ΔₐG⁰

Evaluate at different temperatures:

dtf[:Cp⁰](T = 298.15)   # J/mol/K at 25 °C

41.099821999999996

dtf[:Cp⁰](T = 500.0)    # J/mol/K at 227 °C

45.919999999999995

dtf[:ΔₐG⁰](T = 500.0)   # J/mol at 227 °C

-439937.13910932385

dtf[:ΔₐG⁰](T = 500.0u"K", unit = true)   # with units

-439937.13910932385 m² kg s⁻² mol⁻¹

:logk_fpt_function — van't Hoff log K fit

Fits the equilibrium constant as a polynomial in T:

log₁₀K(T) = A₀ + A₁T + A₂/T + A₃ log T + A₄/T² + A₅T² + A₆/√T

Parameters: A₀…A₆ (dimensionless / appropriate T-powers), T (reference, K).

This model produces a single function :logKr (the fit directly) rather than the full {Cp, H, S, G} set. It is used internally when loading database reactions.

:solute_hkf88_reaktoro — HKF aqueous model

The Helgeson-Kirkham-Flowers (1981/1988) model for aqueous solutes. This is a numeric model (returns NumericFunc objects) because it depends on the equation of state of water and cannot be written as a simple closed-form polynomial.

Parameters (all in SI): a1–a4, c1, c2, wref, z (formal charge), S⁰, ΔfH⁰, ΔfG⁰.

Attaching properties to a species

After building the thermodynamic function dictionary, assign each entry to the species properties dict:

using ChemistryLab
using DynamicQuantities

params_CO2 = Dict(
    :S⁰   => 213.785u"J/K/mol",
    :ΔₐH⁰ => -393510.0u"J/mol",
    :ΔₐG⁰ => -394373.0u"J/mol",
    :a₀   => 33.98u"J/K/mol",
    :a₁   => 23.88e-3u"J/(mol*K^2)",
    :a₂   => 0.0u"J*K/mol",
    :a₃   => 0.0u"J/(mol*K^0.5)",
    :T    => 298.15u"K",
)

CO2 = Species("CO2"; aggregate_state = AS_GAS, class = SC_GASFLUID)
dtf = build_thermo_functions(:cp_ft_equation, params_CO2)
for (k, v) in dtf
    CO2[k] = v
end

CO2[:Cp⁰](T = 400.0)

43.532

CO2[:ΔₐG⁰](T = 400.0u"K", unit = true)

-416804.2494809263 m² kg s⁻² mol⁻¹

Defining a custom model

From a Cp expression (automatic integration)

The simplest way to add a new model is to provide a Julia expression for Cp(T). ChemistryLab will automatically integrate it to obtain H, S, and G.

using ChemistryLab
using DynamicQuantities

# Linear Cp = a + b·T — register with explicit units so build_thermo_functions works
add_thermo_model(:linear_Cp, :(a + b * T), [:T => "K", :a => "J/mol/K", :b => "J/(mol*K^2)"])

# Build functions for a hypothetical species
params = Dict(
    :a    => 28.0u"J/mol/K",
    :b    => 4e-3u"J/(mol*K^2)",
    :S⁰   => 130.0u"J/mol/K",
    :ΔₐH⁰ => -50e3u"J/mol",
    :ΔₐG⁰ => -35e3u"J/mol",
    :T    => 298.15u"K",
)

dtf = build_thermo_functions(:linear_Cp, params)
dtf[:Cp⁰](T = 400.0)

29.6

dtf[:ΔₐG⁰](T = 400.0)

-48700.766558235315

Constant Cp model

using ChemistryLab
using DynamicQuantities

# Cp = Cp0 — constant, independent of T
# Use :(Cp0 * T^0) so the expression is an Expr (not a bare Symbol) and T stays in the vars list
add_thermo_model(:const_Cp, :(Cp0 * T^0), [:T => "K", :Cp0 => "J/mol/K"])

params = Dict(
    :Cp0  => 75.3u"J/mol/K",       # water
    :S⁰   => 69.95u"J/mol/K",
    :ΔₐH⁰ => -285830.0u"J/mol",
    :ΔₐG⁰ => -237140.0u"J/mol",
    :T    => 298.15u"K",
)

dtf = build_thermo_functions(:const_Cp, params)
# Cp is constant → H is linear → S = Cp·ln(T) → G from H and S
dtf[:Cp⁰](T = 373.15)

75.3

dtf[:ΔₐG⁰](T = 373.15)

-243043.52886833143

From a full dictionary

For maximum control, supply all four expressions directly:

using ChemistryLab

add_thermo_model(:my_model, Dict(
    :Cp => :(a + b * T + c / T^2),
    :H  => :(a * T + b * T^2 / 2 - c / T),
    :S  => :(a * log(T) + b * T - c / (2 * T^2)),
    :G  => :(-(a * T * log(T) - a * T) - b * T^2 / 2 - c / (2 * T)),
    :units => [:a => "J/mol/K", :b => "J/(mol*K^2)", :c => "J*K/mol", :T => "K"],
))

ThermoFactory directly

If you need to evaluate a symbolic expression at many parameter values without going through a registered model, use ThermoFactory directly:

using ChemistryLab
using DynamicQuantities

factory = ThermoFactory(:(a + b * T + c / T^2), [:T])

# Stamp out SymbolicFunc instances for different species
f_CO2 = factory(; a = 33.98, b = 23.88e-3, c = 0.0)
f_H2O = factory(; a = 30.54, b = 10.29e-3, c = 0.0)

f_CO2(T = 500.0), f_H2O(T = 500.0)

(45.919999999999995, 35.685)

NumericFunc

NumericFunc wraps a plain Julia closure. The calling convention is identical to SymbolicFunc.

using ChemistryLab
using DynamicQuantities

# Build a NumericFunc manually (e.g. for a piecewise model)
f_piece = NumericFunc(
    (T, P) -> T < 500.0 ? 30.0 + 5e-3 * T : 32.5 + 2e-3 * T,
    (:T, :P),
    (T = 298.15u"K", P = 1e5u"Pa"),
    1.0u"J/mol/K",
)

f_piece(T = 300.0)

31.5

f_piece(T = 600.0)

33.7

NumericFunc objects support the same arithmetic operations as SymbolicFunc. Cross-type arithmetic (e.g. SymbolicFunc + NumericFunc) returns a NumericFunc.

Summary