The file thermo_data.csv contains gas-phase thermodynamic data on several compounds. The comma-separated columns are: chemical formula, InChIKey, standard molar enthalpy of formation (\(\Delta_\mathrm{f}H^\ominus_\mathrm{m}\), in \(\mathrm{kJ\,mol^{-1}}\)), standard molar entropy (\(S^\ominus_\mathrm{m}\), in \(\mathrm{J\,K^{-1}\,mol^{-1}}\)), and the coefficients, \(A - H\) for the temperature dependence of \(\Delta_\mathrm{f}H^\ominus_\mathrm{m}\), \(S^\ominus_\mathrm{m}\) and the standard molar heat capacity, \(C^\ominus_{p, \mathrm{m}}\) according to the Shomate equations:
\begin{align*} C^\ominus_{p, \mathrm{m}} & = A + Bt + Ct^2 + Dt^3 + \frac{E}{t^2},\\ \Delta_\mathrm{f}H^\ominus_\mathrm{m} & = \Delta_\mathrm{f}H^\ominus_\mathrm{m}(298\;\mathrm{K}) + At + \frac{Bt^2}{2} + \frac{Ct^3}{3} + \frac{Dt^4}{4} - \frac{E}{t} + F - H,\\ S^\ominus_\mathrm{m} & = A\ln t + Bt + \frac{Ct^2}{2} + \frac{Dt^3}{3} - \frac{E}{2t^2} + G, \end{align*}
where \(t = (T \,/\,\mathrm{K})/1000\).
We will read in the data file with pandas' read_csv method. Comments are denoted by the '#' character and we would like to index the rows by the chemical formula:
import numpy as np
from scipy.constants import R # Gas constant, J.K-1.mol-1
import pandas as pd
df = pd.read_csv('thermo-data.csv', comment='#', index_col=0)
df
InChIKey DfHo_298 So_298 A ... H
formula
CH4 VNWKT...SA-N -74.87 188.660 -0.703029 ... -74.8731
CO UGFAI...SA-N -110.53 197.660 25.567590 ... -110.5271
... ... ... ... ... ... ...
O2 MYMOF...SA-N 0.00 205.152 30.032350 ... 0.0000
The column headers are taken from the first (non-comment) row of the file, but the spaces after the commas in this line have found their way into the column names:
df.columns
Index([' InChIKey', ' DfHo_298', ' So_298', ' A', ' B', ' C', ' D', ' E', ' F',
' G', ' H'],
dtype='object')
There is an easy way to strip these stray spaces: call the str.strip method on them:
df.columns = df.columns.str.strip()
df.columns
Index(['InChIKey', 'DfHo_298', 'So_298', 'A', 'B', 'C', 'D', 'E', 'F', 'G',
'H'],
dtype='object')
The data are now accessible using regular pandas syntax:
df.loc['CO'] # Thermodynamic data relating to CO
InChIKey UGFAIRIUMAVXCW-UHFFFAOYSA-N
DfHo_298 -110.53
So_298 197.66
A 25.56759
B 6.09613
C 4.054656
D -2.671301
E 0.131021
F -118.0089
G 227.3665
H -110.5271
Name: CO, dtype: object
It will be helpful to define a couple of functions to calculate \(\Delta_\mathrm{f}H^\ominus_\mathrm{m}\) and \(S^\ominus_\mathrm{m}\) for a gas-phase compound specified by its formula, at some arbitrary temperature \(T\), using the above Shomate equations:
def get_DfH(formula, T):
"""Return the standard molar enthalpy of formation of formula at T K."""
s = df.loc[formula]
t = T / 1000
return (s.DfHo_298 + s.A * t + s.B * t**2 / 2 + s.C * t**3 / 3
+ s.D * t**4 / 4 - s.E / t + s.F - s.H)
def get_S(formula, T):
"""Return the standard molar entropy of formula at T K."""
s = df.loc[formula]
t = T / 1000
return (s.A * np.log(t) + s.B * t + s.C * t**2 / 2 + s.D * t**3 / 3
- s.E / 2 / t**2 + s.G)
These formulas can be used to calculate the standard reaction Gibbs free energy change, \(\Delta_\mathrm{r}G^\ominus_\mathrm{m}\), and hence the equilibrium constant of reactions at different temperatures. For example,
$$ \ce{CO(g) + H2O(g) -> H2(g) + CO2(g)} $$
\begin{align*} \Delta_\mathrm{r}H^\ominus_\mathrm{m} &= \Delta_\mathrm{f}H^\ominus_\mathrm{m}(\mathrm{H_2(g)}) + \Delta_\mathrm{f}H^\ominus_\mathrm{m}(\mathrm{CO_2(g)}) - \Delta_\mathrm{f}H^\ominus_\mathrm{m}(\mathrm{CO(g)}) - \Delta_\mathrm{f}H^\ominus_\mathrm{m}(\mathrm{H_2)(g)}),\\ \Delta_\mathrm{r}S^\ominus_\mathrm{m} &= S^\ominus_\mathrm{m}(\mathrm{H_2(g)}) + S^\ominus_\mathrm{m}(\mathrm{CO_2(g)}) - S^\ominus_\mathrm{m}(\mathrm{CO(g)}) - S^\ominus_\mathrm{m}(\mathrm{H_2)(g)}),\\ \Delta_\mathrm{r}G^\ominus_\mathrm{m} &= \Delta_\mathrm{r}H^\ominus_\mathrm{m} - T \Delta_\mathrm{r}S^\ominus_\mathrm{m},\\ K &= \exp \left( -\frac{\Delta_\mathrm{r}G^\ominus_\mathrm{m}}{RT} \right), \end{align*}
where the relevant quantities are calculated at the temperature of interest.
Using the above functions, the molar enthalpy of formation and standard molar entropy of the reactants and products can be calculated at \(T = 850\;\mathrm{K}\):
# Temperature, K
T = 850
# Get the enthalpy and entropy of formation
DfH, DfS = {}, {}
species = ('CO', 'H2O', 'H2', 'CO2')
for formula in species:
DfH[formula] = get_DfH(formula, T)
DfS[formula] = get_DfS(formula, T)
From this, the standard reaction enthalpy, entropy and Gibbs free energy changes can be calculated:
# Standard reaction enthalpy and entropy:
DrH = DfH['H2'] + DfH['CO2'] - DfH['CO'] - DfH['H2O']
DrS = DfS['H2'] + DfS['CO2'] - DfS['CO'] - DfS['H2O']
# Standard Gibbs free energy; NB convert DrS to kJ.K-1.mol-1:
DrG = DrH - T * DrS / 1000
print(f'DrH({T} K) = {DrH:.1f} kJ.mol-1')
print(f'DrS({T} K) = {DrS:.1f} J.K-1.mol-1')
print(f'DrG({T} K) = {DrG:.1f} kJ.mol-1')
DrH(850 K) = -36.3 kJ.mol-1
DrS(850 K) = -33.4 J.K-1.mol-1
DrG(850 K) = -7.9 kJ.mol-1
Finally, the equilibrium constant is;
K = np.exp(-DrG * 1000 / R / T)
print(f'K = {K:.2f}')
K = 3.05
The above code is easily generalized to any reaction:
def get_K(species, T):
"""
Given a set of stoichiometric numbers (a, b, c, ...) and species formulas
(A, B, C, ...) specifying a gas-phase reaction: aA + bB + ... <-> cC + dD + ...
where reactant stoichiometric numbers are negative and products are positive,
return the equilibrium constant for this reaction at temperature T (in K).
e.g. for CO + H2O <-> H2 + CO2,
species = [(-1, 'CO'), (-1, 'H2O'), (1, 'H2'), (1, 'CO2')]
"""
DrH = sum(n * get_DfH(formula, T) for n, formula in species)
DrS = sum(n * get_S(formula, T) for n, formula in species)
DrG = DrH - T * DrS / 1000
K = np.exp(-DrG * 1000 / R / T)
return DrH, DrS, DrG, K
DrH, DrS, DrG, K = get_K([(-1, 'CO'), (-1, 'H2O'), (1, 'H2'), (1, 'CO2')], 850)
print(K)
3.0532234828745612
For a reaction in which the equilibrium lies far to the right (products), consider the combustion of methane:
$$ \ce{CH4(g) + 3/2O2(g) -> CO2(g) + 2H2O(g)} $$
DrH, DrS, DrG, K = get_K([(-1, 'CH4'), (-1.5, 'O2'), (1, 'CO2'), (2, 'H2O')], 850)
print(f'DrH({T} K) = {DrH:.1f} kJ.mol-1')
print(f'DrS({T} K) = {DrS:.1f} J.K-1.mol-1')
print(f'DrG({T} K) = {DrG:.1f} kJ.mol-1')
print(f'K = {K:.2e}')
DrH(850 K) = -791.2 kJ.mol-1
DrS(850 K) = 119.7 J.K-1.mol-1
DrG(850 K) = -892.9 kJ.mol-1
K = 7.38e+54
As expected for a highly exothermic reaction which also involves an increase in system entropy, the equilibrium constant, \(K\), is huge and this reaction proceeds to completion.