04 - Enthalpy of Formation from NASA Polynomials
The standard enthalpy of formation \(\Delta_f H^\circ\) is the enthalpy change when one mole of a species is formed from its elements in their reference states at 298.15 K and 1 bar. By definition it is zero for an element in its reference state (O₂, N₂, H₂, graphite, …).
pyglenn offers two equivalent routes to \(\Delta_f H^\circ\):
calculate_formation_enthalpy(id)— reads the dedicatedheat_of_formation_298Kdatabase column (populated since pyglenn 0.2.0);calculate_properties(id, 298.15).h_relative— the standardized enthalpy, which equals \(\Delta_f H^\circ\) at 298.15 K because NASA polynomials fit the standardized scale.
This notebook validates both routes against literature and shows they agree.
[1]:
from pyglenn import ThermochemicalCalculator, R
print("Universal gas constant R =", R, "J/(mol.K)")
Universal gas constant R = 8.31446261815324 J/(mol.K)
[2]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
plt.rcParams["figure.figsize"] = (8, 4.5)
plt.rcParams["axes.grid"] = True
pd.set_option("display.float_format", lambda v: f"{v:,.3f}")
1. The dedicated method (populated since 0.2.0)
The heat_of_formation_298K column is now populated in the bundled database, so calculate_formation_enthalpy returns the tabulated value directly in J/mol. (Before 0.2.0 this column was empty and the method returned None.)
[3]:
with ThermochemicalCalculator() as calc:
for name in ["H2O", "CO2", "CH4", "O2", "N2"]:
sid = calc.get_available_species(name, exact_match=True)[0].id
dhf = calc.calculate_formation_enthalpy(sid)
print(f"{name:5s} calculate_formation_enthalpy -> {dhf/1000:9.2f} kJ/mol")
H2O calculate_formation_enthalpy -> -241.83 kJ/mol
CO2 calculate_formation_enthalpy -> -393.51 kJ/mol
CH4 calculate_formation_enthalpy -> -74.60 kJ/mol
O2 calculate_formation_enthalpy -> 0.00 kJ/mol
N2 calculate_formation_enthalpy -> 0.00 kJ/mol
2. The same value from the polynomial
The standardized enthalpy already contains the formation enthalpy, so at 298.15 K it is \(\Delta_f H^\circ\) — the same quantity the dedicated method returns. Reference-state elements come out at essentially zero (to within rounding of the fit).
[4]:
def formation_enthalpy(calc, name, T=298.15):
"""Delta_f H at T (default 298.15 K), in J/mol, from the standardized H."""
return calc.calculate_properties(calc.get_available_species(name, exact_match=True)[0].id, T).h_relative
with ThermochemicalCalculator() as calc:
print("Reference-state elements (should be ~0):")
for el in ["O2", "N2", "H2", "C(gr)"]:
print(f" {el:6s} {formation_enthalpy(calc, el)/1000: .4f} kJ/mol")
print("\nDedicated method vs polynomial recovery:")
for name in ["H2O", "CO2", "CH4", "NH3", "NO"]:
sid = calc.get_available_species(name, exact_match=True)[0].id
dedicated = calc.calculate_formation_enthalpy(sid) / 1000.0
poly = formation_enthalpy(calc, name) / 1000.0
print(f" {name:5s} dedicated={dedicated:9.2f} polynomial={poly:9.2f} diff={abs(dedicated-poly):.2e} kJ/mol")
Reference-state elements (should be ~0):
O2 -0.0000 kJ/mol
N2 0.0000 kJ/mol
H2 -0.0000 kJ/mol
C(gr) 0.0000 kJ/mol
Dedicated method vs polynomial recovery:
H2O dedicated= -241.83 polynomial= -241.83 diff=3.40e-07 kJ/mol
CO2 dedicated= -393.51 polynomial= -393.51 diff=3.21e-08 kJ/mol
CH4 dedicated= -74.60 polynomial= -74.60 diff=1.39e-08 kJ/mol
NH3 dedicated= -45.94 polynomial= -45.94 diff=2.52e-08 kJ/mol
NO dedicated= 91.27 polynomial= 91.27 diff=1.24e-08 kJ/mol
3. A validated table of formation enthalpies
We compute \(\Delta_f H^\circ\) for a set of common species and compare against accepted literature values (CODATA / JANAF, gas phase). The agreement is well within 1 kJ/mol — the polynomials reproduce the reference data they were fit to.
[5]:
# Literature standard enthalpies of formation at 298.15 K, kJ/mol (gas phase)
REFERENCE = {
"H2": 0.0, "O2": 0.0, "N2": 0.0, "C(gr)": 0.0,
"H2O": -241.83, "CO2": -393.52, "CO": -110.53, "CH4": -74.87,
"NH3": -45.90, "NO": 91.29, "NO2": 33.10, "C2H5OH": -235.0,
}
records = []
with ThermochemicalCalculator() as calc:
for name, ref in REFERENCE.items():
sid = calc.get_available_species(name, exact_match=True)[0].id
M = calc.db.get_species_data(sid).molecular_weight
dhf = formation_enthalpy(calc, name) / 1000.0
records.append({
"species": name, "M [g/mol]": M,
"pyglenn [kJ/mol]": dhf, "reference [kJ/mol]": ref,
"abs.err [kJ/mol]": abs(dhf - ref),
})
hf_df = pd.DataFrame(records).set_index("species")
print(hf_df.to_string())
print(f"\nLargest absolute error: {hf_df['abs.err [kJ/mol]'].max():.3f} kJ/mol")
M [g/mol] pyglenn [kJ/mol] reference [kJ/mol] abs.err [kJ/mol]
species
H2 2.016 -0.000 0.000 0.000
O2 31.999 -0.000 0.000 0.000
N2 28.013 0.000 0.000 0.000
C(gr) 12.011 0.000 0.000 0.000
H2O 18.015 -241.826 -241.830 0.004
CO2 44.010 -393.510 -393.520 0.010
CO 28.010 -110.535 -110.530 0.005
CH4 16.042 -74.600 -74.870 0.270
NH3 17.031 -45.940 -45.900 0.040
NO 30.006 91.271 91.290 0.019
NO2 46.005 34.193 33.100 1.093
C2H5OH 46.068 -234.950 -235.000 0.050
Largest absolute error: 1.093 kJ/mol
4. Visualising formation enthalpies
A horizontal bar chart makes the exothermic (negative \(\Delta_f H^\circ\), stable) versus endothermic (positive, energetic) split obvious. Elements sit exactly on the zero line.
[6]:
order = hf_df.sort_values("pyglenn [kJ/mol]")
vals = order["pyglenn [kJ/mol]"]
colors = ["#c0392b" if v > 0 else ("#7f8c8d" if abs(v) < 1e-3 else "#2471a3")
for v in vals]
fig, ax = plt.subplots(figsize=(8, 5))
ax.barh(order.index, vals, color=colors)
ax.axvline(0, color="k", lw=0.8)
ax.set_xlabel(r"$\Delta_f H^\circ$ at 298.15 K [kJ/mol]")
ax.set_title("Standard enthalpies of formation (from pyglenn)")
plt.show()
5. Only 298.15 K is “formation”
\(\Delta_f H^\circ\) is defined at 298.15 K. The standardized enthalpy \(H^\circ(T)\) keeps varying with temperature, but only its 298.15 K value equals the tabulated formation enthalpy. For CO₂ we plot \(H^\circ(T)\) and highlight that reference point; the vertical offset above it is the sensible enthalpy \(H^\circ(T)-H^\circ(298.15)\) from notebook 03.
[7]:
Tgrid = np.linspace(300, 2500, 60)
with ThermochemicalCalculator() as calc:
co2 = calc.get_available_species("CO2", exact_match=True)[0].id
H = np.array([calc.calculate_properties(co2, t).h_relative for t in Tgrid]) / 1000.0
hf298 = formation_enthalpy(calc, "CO2") / 1000.0
fig, ax = plt.subplots()
ax.plot(Tgrid, H, label=r"$H^\circ(T)$ of CO$_2$")
ax.scatter([298.15], [hf298], color="red", zorder=5,
label=fr"$\Delta_f H^\circ$ = {hf298:.1f} kJ/mol at 298.15 K")
ax.set_xlabel("Temperature [K]")
ax.set_ylabel("Standardized enthalpy [kJ/mol]")
ax.set_title("Formation enthalpy is the 298.15 K value of the standardized enthalpy")
ax.legend()
plt.show()
Summary
calculate_formation_enthalpynow returns the tabulated \(\Delta_f H^\circ\) (theheat_of_formation_298Kcolumn is populated since 0.2.0).calculate_properties(id, 298.15).h_relativegives the same value — NASA’s standardized enthalpy carries the formation enthalpy.Both routes reproduce literature values to < 1 kJ/mol.
Next: notebook 05 uses exactly this to compute reaction enthalpies and heats of combustion.