Gas-Phase Thermochemistry Skill

Goal

Compute temperature-dependent thermodynamic quantities — enthalpy ($H$), entropy ($S$), and Gibbs free energy ($G$) — for gas-phase molecules, and reaction thermochemistry ($\Delta H$, $\Delta S$, $\Delta G$) for balanced chemical reactions, using Machine Learning Interatomic Potentials (MLIPs) with ASE's ideal-gas / rigid-rotor / harmonic-oscillator (IGRRHO) framework.

[!IMPORTANT] This skill extends chem-vibration. It reuses the same MLIP-powered Hessian calculation and adds thermodynamic post-processing via ase.thermochemistry.IdealGasThermo.

Background

The IGRRHO partition function decomposes into independent translational, rotational, vibrational, and electronic contributions:

$$q = q_\text{trans} \cdot q_\text{rot} \cdot q_\text{vib} \cdot q_\text{elec}$$

From the partition function, thermodynamic quantities at temperature $T$ and pressure $P$ are:

Enthalpy: $H(T) = E_\text{pot} + E_\text{ZPE} + \sum_i \frac{h\nu_i}{e^{h\nu_i/k_BT}-1} + k_BT$ (+ rotational and translational terms)
Entropy: $S(T,P) = S_\text{trans} + S_\text{rot} + S_\text{vib} + S_\text{elec}$
Gibbs free energy: $G(T,P) = H(T) - TS(T,P)$

For a balanced reaction, $\Delta X = \sum_\text{products} n_i X_i - \sum_\text{reactants} n_j X_j$ where $X = H, S, G$.

1. Prerequisites

An MLIP wrapper must be available (MACEWrapper, MatGLWrapper, or FAIRCHEMWrapper).
ASE must be installed (provides ase.thermochemistry.IdealGasThermo).
Structures must be molecules or gas-phase species (non-periodic).

2. Choosing a Foundation Potential

Refer to the foundation-potentials skill for model selection.

[!IMPORTANT] For molecular thermochemistry, use models with good force accuracy (e.g., MACE-OMAT-0-small, MACE-MH-1 with omol head). Accurate forces are critical for reliable vibrational frequencies and thus thermodynamic quantities.

3. Calculation Workflow

Single-Molecule Mode

Compute absolute H(T), S(T,P), G(T,P) for one species:

# Env: mace-agent
python .agents/skills/chem-thermochemistry/scripts/calculate_thermochemistry.py \
    --molecule H2O \
    --temperature 298.15 \
    --pressure 101325 \
    --model_type mace \
    --model_name MACE-OMAT-0-small \
    --output_dir research/my_folder/thermo

Reaction Mode

Compute ΔH, ΔS, ΔG for a balanced gas-phase reaction:

# Env: mace-agent
python .agents/skills/chem-thermochemistry/scripts/calculate_thermochemistry.py \
    --reaction "2H2 + O2 -> 2H2O" \
    --temperature 298.15 \
    --model_type mace \
    --model_name MACE-OMAT-0-small \
    --output_dir research/my_folder/reaction_thermo

Key Parameters

| Argument | Description | |:---|:---| | --molecule | ASE built-in molecule name (single species mode) | | --structure | Path to structure file (single species mode) | | --reaction | Balanced reaction string, e.g., "2H2 + O2 -> 2H2O" | | --temperature | Temperature in K (default: 298.15) | | --pressure | Pressure in Pa (default: 101325) | | --model_type | MLIP backend: mace, matgl, fairchem | | --model_name | Model name (e.g., MACE-OMAT-0-small) | | --spin | Spin multiplicity override as "name:value" pairs (e.g., "O2:1") | | --symmetry_number | Symmetry number override as "name:value" pairs (e.g., "H2:2") | | --fmax | Force convergence for relaxation (default: 0.01 eV/Å) | | --output_dir | Output directory |

4. Output Files

thermochemistry_results.json: Full results including:
- Per species: potential energy, ZPE, vibrational frequencies, H(T), S(T,P), G(T,P)
- Reaction (if applicable): ΔH, ΔS, ΔG in both eV and kJ/mol
- Metadata: temperature, pressure, model, spin, symmetry numbers

5. Examples

H₂O Formation Reaction

See examples/H2O_formation/ for the water formation reaction:

# Env: mace-agent
python .agents/skills/chem-thermochemistry/scripts/calculate_thermochemistry.py \
    --reaction "2H2 + O2 -> 2H2O" \
    --temperature 298.15 \
    --model_type mace \
    --model_name MACE-OMAT-0-small \
    --output_dir .agents/skills/chem-thermochemistry/examples/H2O_formation

NIST reference values for 2H₂(g) + O₂(g) → 2H₂O(g) at 298.15 K: | Quantity | NIST Value | |:---|:---| | ΔH | −483.65 kJ/mol | | ΔG | −457.22 kJ/mol |

6. Constraints

Gas-phase only: This skill uses the ideal-gas approximation. Not applicable to condensed-phase or surface reactions.
Harmonic approximation: Vibrational contributions assume harmonic potentials. Accuracy degrades for floppy modes and near dissociation.
Spin and symmetry: The script includes a lookup table for common molecules, but exotic species may need manual --spin and --symmetry_number overrides.
Environments: Scripts require conda environments with MLIP packages:
- mace-agent for MACE models
- matgl-agent for MatGL/CHGNet models
- fairchem-agent for FairChem/UMA models

References

ASE Thermochemistry module: ASE Documentation
NIST-JANAF Thermochemical Tables: NIST
McQuarrie, D.A. "Statistical Mechanics", University Science Books, 2000.

Author: Bowen Deng Contact: GitHub @learningmatter-mit