← Part IV/Partition Functions

Partition Functions

Reading time: ~45 minutes | Derivations: 5 | Simulations: 2

1. Introduction: The Bridge Between Worlds

The partition function is the single most important quantity in statistical mechanics. It encodes all thermodynamic information about a system in one mathematical object — a sum over Boltzmann-weighted states that connects the microscopic energy spectrum to macroscopic observables.

The Canonical Partition Function

For a system with discrete energy levels $\{E_i\}$, the canonical partition function is:

$$Z = \sum_i g_i \, e^{-E_i / k_B T} = \sum_i g_i \, e^{-\beta E_i}$$

where $\beta = 1/(k_B T)$ is the thermodynamic beta, $g_i$ is the degeneracy of level $i$, and the sum runs over all distinct energy levels.

The German word Zustandssumme ("sum over states") captures the essence: we sum over every accessible quantum state, each weighted by its Boltzmann factor. At low temperature, only the ground state contributes significantly. As $T \to \infty$, all states participate equally, and $Z$ grows without bound.

For a single molecule, we write the molecular partition function as $q$. For $N$ indistinguishable, non-interacting molecules, the system partition function is:

$$Q = \frac{q^N}{N!}$$

The $N!$ accounts for the indistinguishability of identical particles (Gibbs's correction), preventing the Gibbs paradox in which the entropy of mixing two identical gases would be nonzero.

Why the Partition Function Matters

Every thermodynamic quantity can be derived from $Z$ by differentiation
It provides the normalization for the Boltzmann distribution
The logarithm of $Z$ is directly proportional to the Helmholtz free energy
Chemical equilibrium constants follow from ratios of partition functions
It naturally incorporates quantum effects through the energy level structure

2. Derivation: Translational Partition Function

A molecule of mass $m$ confined to a cubic box of side $L$(volume $V = L^3$) has particle-in-a-box energy levels in three dimensions:

$$E_{n_x, n_y, n_z} = \frac{h^2}{8mL^2}\left(n_x^2 + n_y^2 + n_z^2\right), \quad n_x, n_y, n_z = 1, 2, 3, \ldots$$

The translational partition function is therefore:

$$q_{\text{trans}} = \sum_{n_x=1}^{\infty} \sum_{n_y=1}^{\infty} \sum_{n_z=1}^{\infty} \exp\!\left[-\frac{h^2(n_x^2 + n_y^2 + n_z^2)}{8mL^2 k_B T}\right]$$

Because the energy is a sum of independent terms, the triple sum factors:

$$q_{\text{trans}} = \left(\sum_{n=1}^{\infty} e^{-n^2 h^2 / (8mL^2 k_B T)}\right)^{\!3}$$

High-Temperature (Continuum) Limit

Define $\alpha = h^2/(8mL^2 k_B T)$. For typical molecules at room temperature,$\alpha \sim 10^{-18}$, so the spacing between levels is vastly smaller than $k_B T$. We replace the sum by an integral:

$$\sum_{n=1}^{\infty} e^{-\alpha n^2} \approx \int_0^{\infty} e^{-\alpha n^2}\, dn = \frac{1}{2}\sqrt{\frac{\pi}{\alpha}}$$

This is the Gaussian integral $\int_0^\infty e^{-ax^2}\,dx = \frac{1}{2}\sqrt{\pi/a}$.

Substituting back and cubing:

$$q_{\text{trans}} = \left(\frac{1}{2}\sqrt{\frac{\pi}{\alpha}}\right)^3 = \left(\frac{2\pi m k_B T}{h^2}\right)^{3/2} L^3 = \left(\frac{2\pi m k_B T}{h^2}\right)^{3/2} V$$

Thermal de Broglie Wavelength

Define the thermal de Broglie wavelength:

$$\Lambda = \frac{h}{\sqrt{2\pi m k_B T}} = h\left(\frac{\beta}{2\pi m}\right)^{1/2}$$

Then the translational partition function takes the elegant form:

$$\boxed{q_{\text{trans}} = \frac{V}{\Lambda^3}}$$

Physical interpretation: $q_{\text{trans}}$ counts the number of thermal de Broglie volumes that fit inside the container. For an ideal gas at STP,$q_{\text{trans}} \sim 10^{30}$ — an enormous number, confirming that the classical (high-temperature) limit is valid.

3. Derivation: Rotational Partition Function

For a rigid rotor (linear molecule), the rotational energy levels are:

$$E_J = \frac{\hbar^2}{2I} J(J+1) = k_B \Theta_{\text{rot}} \, J(J+1), \quad J = 0, 1, 2, \ldots$$

where $I$ is the moment of inertia and we define the characteristic rotational temperature:

$$\Theta_{\text{rot}} = \frac{\hbar^2}{2I k_B} = \frac{hB}{k_B}$$

where $B$ is the rotational constant in frequency units. For N$_2$,$\Theta_{\text{rot}} \approx 2.88$ K; for H$_2$,$\Theta_{\text{rot}} \approx 87.6$ K.

Each level $J$ has degeneracy $g_J = 2J + 1$ (from the magnetic quantum number $m_J = -J, \ldots, +J$). The rotational partition function is:

$$q_{\text{rot}} = \frac{1}{\sigma} \sum_{J=0}^{\infty} (2J+1) \, e^{-J(J+1)\Theta_{\text{rot}}/T}$$

The symmetry number $\sigma$ corrects for overcounting of indistinguishable orientations: $\sigma = 1$ for heteronuclear diatomics (CO, HCl),$\sigma = 2$ for homonuclear diatomics (N$_2$, O$_2$).

High-Temperature Limit

When $T \gg \Theta_{\text{rot}}$, many rotational levels are populated and we may replace the sum by an integral. Let $x = J(J+1)$, so$dx = (2J+1)\,dJ$:

$$q_{\text{rot}} \approx \frac{1}{\sigma}\int_0^{\infty} e^{-x\,\Theta_{\text{rot}}/T}\,dx = \frac{1}{\sigma}\cdot\frac{T}{\Theta_{\text{rot}}}$$

$$\boxed{q_{\text{rot}} = \frac{T}{\sigma\,\Theta_{\text{rot}}}} \qquad (T \gg \Theta_{\text{rot}})$$

At 300 K for N$_2$: $q_{\text{rot}} \approx 300/(2 \times 2.88) \approx 52$.

Extension to Polyatomic Molecules

Symmetric Top

For a symmetric top with moments of inertia $I_A = I_B \neq I_C$:

$$q_{\text{rot}}^{\text{sym}} = \frac{1}{\sigma}\sqrt{\frac{\pi T^3}{\Theta_A^2 \Theta_C}}$$

where $\Theta_A = \hbar^2/(2I_A k_B)$ and similarly for $\Theta_C$.

Asymmetric Top

For a general asymmetric top with three distinct moments $I_A, I_B, I_C$:

$$q_{\text{rot}}^{\text{asym}} = \frac{1}{\sigma}\sqrt{\frac{\pi T^3}{\Theta_A \Theta_B \Theta_C}}$$

Water is an asymmetric top with $\sigma = 2$: $\Theta_A = 40.1$ K,$\Theta_B = 20.9$ K, $\Theta_C = 13.4$ K. At 300 K,$q_{\text{rot}} \approx 43$.

Note: Rotational contributions give $C_V = k_B$ per molecule for linear molecules (2 rotational degrees of freedom) and $C_V = \frac{3}{2}k_B$for nonlinear molecules (3 rotational degrees of freedom) in the classical limit.

4. Derivation: Vibrational Partition Function

In the harmonic approximation, a diatomic molecule has vibrational energy levels:

$$E_v = \left(v + \frac{1}{2}\right)h\nu, \quad v = 0, 1, 2, \ldots$$

Measuring energies from the zero-point level (setting $E_0 = 0$), the vibrational partition function becomes:

$$q_{\text{vib}} = \sum_{v=0}^{\infty} e^{-v h\nu/(k_B T)} = \sum_{v=0}^{\infty} \left(e^{-\Theta_{\text{vib}}/T}\right)^v$$

where we define the characteristic vibrational temperature:

$$\Theta_{\text{vib}} = \frac{h\nu}{k_B}$$

For N$_2$: $\Theta_{\text{vib}} \approx 3393$ K; for HCl: $\Theta_{\text{vib}} \approx 4303$ K.

This is a geometric series with ratio $r = e^{-\Theta_{\text{vib}}/T} < 1$. Using $\sum_{v=0}^{\infty} r^v = 1/(1-r)$:

$$\boxed{q_{\text{vib}} = \frac{1}{1 - e^{-\Theta_{\text{vib}}/T}}}$$

Vibrational Energy and Heat Capacity

The mean vibrational energy (per molecule, above zero-point) is obtained from$\langle E \rangle = -\partial \ln q_{\text{vib}} / \partial \beta$:

$$\langle E_{\text{vib}} \rangle = \frac{k_B \Theta_{\text{vib}}}{e^{\Theta_{\text{vib}}/T} - 1}$$

The vibrational contribution to the heat capacity is:

$$\frac{C_{V,\text{vib}}}{k_B} = \left(\frac{\Theta_{\text{vib}}}{T}\right)^2 \frac{e^{\Theta_{\text{vib}}/T}}{\left(e^{\Theta_{\text{vib}}/T} - 1\right)^2}$$

This is the Einstein function. Key limiting behaviors:

$T \gg \Theta_{\text{vib}}$: $C_{V,\text{vib}} \to k_B$ (classical equipartition)
$T \ll \Theta_{\text{vib}}$: $C_{V,\text{vib}} \to k_B (\Theta_{\text{vib}}/T)^2 e^{-\Theta_{\text{vib}}/T} \to 0$ (mode frozen out)

Polyatomic molecules: A nonlinear molecule with $N$ atoms has$3N - 6$ vibrational modes (or $3N - 5$ for linear). Each mode has its own $\Theta_{\text{vib},i}$, and the total vibrational partition function is $q_{\text{vib}} = \prod_{i} q_{\text{vib},i}$.

5. Derivation: Electronic Partition Function

The electronic partition function sums over electronic energy levels, each with its own degeneracy:

$$q_{\text{elec}} = g_0 + g_1 e^{-\varepsilon_1/(k_B T)} + g_2 e^{-\varepsilon_2/(k_B T)} + \cdots$$

where $g_i$ is the degeneracy of the $i$-th electronic state and$\varepsilon_i$ is its energy above the ground state.

The Usual Simplification

For most molecules, the first excited electronic state lies tens of thousands of cm$^{-1}$above the ground state, corresponding to $\Theta_{\text{elec}} \sim 10^4$–$10^5$ K. At ordinary temperatures, all excited-state Boltzmann factors are negligible:

$$\boxed{q_{\text{elec}} \approx g_0}$$

Examples: N$_2$ ($g_0 = 1$, $^1\Sigma_g^+$), O$_2$ ($g_0 = 3$, $^3\Sigma_g^-$), Cl ($g_0 = 4$, $^2P_{3/2}$).

Exception: Nitric Oxide (NO)

NO is a remarkable exception. Its ground state $^2\Pi_{1/2}$ ($g_0 = 2$) has a low-lying excited state $^2\Pi_{3/2}$ ($g_1 = 2$) at only$\varepsilon_1 = 121.1$ cm$^{-1}$, corresponding to$\Theta_{\text{elec}} \approx 174$ K. This is comparable to room temperature!

$$q_{\text{elec}}^{\text{NO}} = 2 + 2\,e^{-174/T}$$

At 300 K: $q_{\text{elec}}^{\text{NO}} \approx 2 + 2(0.562) = 3.12$. The excited state contributes significantly to the partition function and therefore to the entropy and heat capacity of NO gas.

Electronic Contribution to Heat Capacity

For a two-level electronic system with degeneracies $g_0, g_1$ and gap $\varepsilon_1$:

$$\frac{C_{V,\text{elec}}}{k_B} = \frac{g_1}{g_0}\left(\frac{\varepsilon_1}{k_B T}\right)^2 \frac{e^{-\varepsilon_1/(k_BT)}}{\left(1 + \frac{g_1}{g_0}e^{-\varepsilon_1/(k_BT)}\right)^2}$$

This is a Schottky anomaly: the heat capacity has a maximum at$k_BT \approx 0.42\,\varepsilon_1$ and vanishes at both low and high temperatures.

6. Derivation: Total Molecular Partition Function

The key assumption enabling the factorization of the molecular partition function is the separability of degrees of freedom. If the total molecular energy can be written as:

$$\varepsilon = \varepsilon_{\text{trans}} + \varepsilon_{\text{rot}} + \varepsilon_{\text{vib}} + \varepsilon_{\text{elec}}$$

then because the exponential of a sum is a product of exponentials, the molecular partition function factors:

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

This factorization is approximate: rotation-vibration coupling, centrifugal distortion, and electronic-vibrational interactions break strict separability. These corrections are typically small for low-lying states.

Thermodynamic Quantities from the Partition Function

For $N$ indistinguishable molecules, $Q = q^N/N!$. Using Stirling's approximation ($\ln N! \approx N\ln N - N$):

$$\ln Q = N\ln q - N\ln N + N$$

Internal Energy

$$E = -\frac{\partial \ln Q}{\partial \beta}\bigg|_{V} = -N\frac{\partial \ln q}{\partial \beta}\bigg|_{V}$$

Equivalently, $E = k_B T^2 \frac{\partial \ln Q}{\partial T}\Big|_V$.

Helmholtz Free Energy

$$\boxed{A = -k_B T \ln Q}$$

This is the fundamental bridge equation of statistical thermodynamics. All other thermodynamic potentials follow from $A$.

Entropy

From $S = -(\partial A/\partial T)_V$ or equivalently $S = (E - A)/T$:

$$\boxed{S = k_B \ln Q + \frac{E}{T}}$$

For an ideal gas, this yields the Sackur-Tetrode equation for translational entropy:

$$S_{\text{trans}} = Nk_B\left[\ln\left(\frac{V}{N\Lambda^3}\right) + \frac{5}{2}\right]$$

Pressure

$$P = -\left(\frac{\partial A}{\partial V}\right)_T = k_B T \left(\frac{\partial \ln Q}{\partial V}\right)_T$$

Since only $q_{\text{trans}} \propto V$, this immediately gives the ideal gas law:$P = Nk_BT/V$.

Chemical Potential

$$\mu = -k_B T \ln\left(\frac{q}{N}\right) = -k_B T \ln\left(\frac{q}{N}\right)$$

This connects to chemical equilibrium: for the reaction $aA + bB \rightleftharpoons cC + dD$, the equilibrium constant is:

$$K = \frac{(q_C/V)^c (q_D/V)^d}{(q_A/V)^a (q_B/V)^b} \, e^{-\Delta E_0/(k_BT)}$$

7. Applications

Heat Capacities of Gases

The heat capacity of an ideal gas reveals which degrees of freedom are active at a given temperature. For a diatomic molecule:

Low T ($T \ll \Theta_{\text{rot}}$): Only translation is active. $C_V = \frac{3}{2}Nk_B$, giving$\gamma = C_P/C_V = 5/3$ (monatomic-like).
Intermediate T ($\Theta_{\text{rot}} \ll T \ll \Theta_{\text{vib}}$): Translation + rotation active. $C_V = \frac{5}{2}Nk_B$, giving$\gamma = 7/5 = 1.4$. This is the regime for N$_2$ and O$_2$ at room temperature.
High T ($T \gg \Theta_{\text{vib}}$): All modes active. $C_V = \frac{7}{2}Nk_B$, giving $\gamma = 9/7 \approx 1.29$.

Chemical Equilibrium from Partition Functions

The equilibrium constant for any gas-phase reaction can be computed from molecular partition functions without any empirical input beyond spectroscopic data. For example, the dissociation of I$_2$:

$$\text{I}_2(g) \rightleftharpoons 2\text{I}(g)$$

$$K_P = \frac{k_BT(q_{\text{I}}/V)^2}{q_{\text{I}_2}/V}\,e^{-D_0/(k_BT)}$$

where $D_0$ is the dissociation energy from the zero-point level. Computing each $q$ from spectroscopic constants gives excellent agreement with experimental equilibrium constants.

Spectroscopic Determination of q

The molecular parameters entering partition functions — rotational constants ($B$), vibrational frequencies ($\nu$), electronic term values — are all measured spectroscopically. High-resolution infrared and microwave spectroscopy provide these constants to extraordinary precision (often 8+ significant figures for $B$), making statistical mechanical calculations of thermodynamic properties more accurate than direct calorimetric measurements for many small molecules.

8. Historical Context

Ludwig Boltzmann (1844–1906)

Boltzmann laid the foundations of statistical mechanics in the 1870s. His expression $S = k_B \ln W$ (carved on his tombstone) connects entropy to the number of microstates. He introduced the concept of the partition function as the normalization of the probability distribution over states, though the term "Zustandssumme" came later from Planck.

Max Planck (1858–1947)

In his derivation of the blackbody radiation law (1900), Planck introduced energy quantization $E = nh\nu$ and computed the partition function for quantum oscillators. This was the first quantum partition function, though Planck did not initially appreciate its revolutionary implications. His result for the average energy of an oscillator, $\langle E \rangle = h\nu/(e^{h\nu/k_BT} - 1)$, is precisely the vibrational contribution we derived above.

Albert Einstein (1879–1955)

Einstein's 1907 paper on the specific heat of solids was a triumph of partition function methods. By modeling a solid as $3N$ independent quantum oscillators all at the same frequency, he derived the Einstein heat capacity function that correctly predicts $C_V \to 0$ as $T \to 0$, resolving the classical equipartition catastrophe. Though quantitatively imperfect at very low temperatures (due to the single-frequency assumption), it demonstrated the power of quantum statistical mechanics.

Peter Debye (1884–1966)

Debye improved upon Einstein's model in 1912 by recognizing that a solid has a continuous spectrum of vibrational frequencies up to a cutoff $\nu_D$(the Debye frequency). His model gives $C_V \propto T^3$ at low temperatures, in excellent agreement with experiment. The Debye model remains the standard starting point for understanding the thermal properties of solids.

9. Simulation: Heat Capacity Contributions (Python)

This simulation computes and plots the individual contributions of translational, rotational, vibrational, and electronic modes to the constant-volume heat capacity $C_V/Nk_B$ as a function of temperature for two molecules: N$_2$ (linear diatomic) and H$_2$O (nonlinear triatomic). The rotational partition function for N$_2$is computed by exact summation over $J$ levels, capturing the quantum freeze-out at low temperature. Characteristic temperatures $\Theta_{\text{rot}}$ and$\Theta_{\text{vib}}$ are printed for comparison.

Heat Capacity from Partition Functions: N₂ and H₂O

Python

Computes translational, rotational, vibrational, and electronic contributions to Cv/NkB vs temperature for N₂ and H₂O. Characteristic temperatures and room-temperature values are printed.

partition_heat_capacity.py192 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ──────────────────────────────────────────────────────────────
# Partition Function Contributions to Heat Capacity
# Compute C_v/Nk_B from translational, rotational, vibrational,
# and electronic modes for N2 and H2O vs temperature
# ──────────────────────────────────────────────────────────────

k_B = 1.380649e-23    # J/K
h   = 6.62607015e-34  # J s
c   = 2.99792458e10   # cm/s
N_A = 6.02214076e23   # 1/mol

# Temperature range
T = np.linspace(50, 6000, 800)

# ── N2 (diatomic, linear) ──
# Molecular parameters
m_N2      = 28.014e-3 / N_A        # kg per molecule
B_N2      = 1.99824              # cm^-1 (rotational constant)
nu_N2     = 2358.57              # cm^-1 (vibrational frequency)
sigma_N2  = 2                    # symmetry number
g0_N2     = 1                    # ground electronic degeneracy
g1_N2     = 3                    # first excited state degeneracy (a 1Pi_g)
eps1_N2   = 69283.0              # cm^-1 first excited electronic state

Theta_rot_N2  = h * c * B_N2 / k_B          # rotational temperature
Theta_vib_N2  = h * c * nu_N2 / k_B         # vibrational temperature
Theta_elec_N2 = h * c * eps1_N2 / k_B       # electronic excitation temperature

# ── H2O (nonlinear, polyatomic) ──
m_H2O     = 18.015e-3 / N_A
# Three rotational constants for asymmetric top
A_H2O = 27.877   # cm^-1
B_H2O = 14.512   # cm^-1
C_H2O = 9.285    # cm^-1
sigma_H2O = 2    # symmetry number

Theta_A = h * c * A_H2O / k_B
Theta_B = h * c * B_H2O / k_B
Theta_C = h * c * C_H2O / k_B

# Three vibrational modes: symmetric stretch, bend, asymmetric stretch
nu_H2O = np.array([3657.05, 1594.78, 3755.93])  # cm^-1
Theta_vib_H2O = h * c * nu_H2O / k_B

g0_H2O = 1

# ──────────────────────────────────────
# Heat capacity contributions (C_v / k_B per molecule)
# ──────────────────────────────────────

# Translational: always 3/2 k_B
Cv_trans = np.full_like(T, 1.5)

# ── N2 rotational (high-T limit for linear molecule) ──
# Exact: sum over J with degeneracy (2J+1)
# C_rot/k_B from numerical differentiation of <E_rot>/k_BT
J_max = 200
J_arr = np.arange(0, J_max + 1)
# Energy: E_J = J(J+1) Theta_rot
# Degeneracy: (2J+1)

def rotational_Cv_linear(T_arr, Theta_rot, sigma):
    Cv = np.zeros_like(T_arr)
    for i, temp in enumerate(T_arr):
        if temp < 1.0:
            Cv[i] = 0.0
            continue
        x = Theta_rot / temp
        weights = (2 * J_arr + 1) * np.exp(-J_arr * (J_arr + 1) * x)
        E1 = np.sum(J_arr * (J_arr + 1) * x * weights) / np.sum(weights)
        E2 = np.sum((J_arr * (J_arr + 1) * x)**2 * weights) / np.sum(weights)
        Cv[i] = E2 - E1**2
    return Cv

Cv_rot_N2 = rotational_Cv_linear(T, Theta_rot_N2, sigma_N2)

# H2O rotational (nonlinear, high-T limit: 3/2 k_B)
# For asymmetric top, C_rot -> 3/2 at high T
# Use interpolation: frozen at low T, classical at high T
Theta_rot_avg_H2O = (Theta_A * Theta_B * Theta_C)**(1.0/3.0)
Cv_rot_H2O = 1.5 * (1.0 - np.exp(-T / (3.0 * Theta_rot_avg_H2O)))

# ── Vibrational heat capacity ──
def vibrational_Cv(T_arr, Theta_vib_arr):
    Cv = np.zeros_like(T_arr)
    for Theta_v in np.atleast_1d(Theta_vib_arr):
        x = Theta_v / T_arr
        # Avoid overflow
        mask = x < 500
        contrib = np.zeros_like(T_arr)
        contrib[mask] = x[mask]**2 * np.exp(x[mask]) / (np.exp(x[mask]) - 1.0)**2
        Cv += contrib
    return Cv

Cv_vib_N2  = vibrational_Cv(T, Theta_vib_N2)
Cv_vib_H2O = vibrational_Cv(T, Theta_vib_H2O)

# ── Electronic heat capacity ──
def electronic_Cv(T_arr, g0, g1, Theta_elec):
    x = Theta_elec / T_arr
    mask = x < 500
    Cv = np.zeros_like(T_arr)
    ratio = np.zeros_like(T_arr)
    ratio[mask] = (g1 / g0) * np.exp(-x[mask])
    Z = 1.0 + ratio
    frac = ratio / Z
    Cv = frac * x**2 * (1.0 - frac)
    return Cv

Cv_elec_N2 = electronic_Cv(T, g0_N2, g1_N2, Theta_elec_N2)
Cv_elec_H2O = np.zeros_like(T)  # negligible for H2O

# ── Total ──
Cv_total_N2  = Cv_trans + Cv_rot_N2 + Cv_vib_N2 + Cv_elec_N2
Cv_total_H2O = Cv_trans + Cv_rot_H2O + Cv_vib_H2O + Cv_elec_H2O

# ──────────────────────────────────────
# Plotting
# ──────────────────────────────────────
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# N2 plot
ax = axes[0]
ax.plot(T, Cv_trans,    '--', color='#60a5fa', linewidth=1.5, label='Translational (3/2)')
ax.plot(T, Cv_rot_N2,   '--', color='#34d399', linewidth=1.5, label='Rotational')
ax.plot(T, Cv_vib_N2,   '--', color='#f87171', linewidth=1.5, label='Vibrational')
ax.plot(T, Cv_elec_N2,  '--', color='#fbbf24', linewidth=1.5, label='Electronic')
ax.plot(T, Cv_total_N2, '-',  color='white',   linewidth=2.5, label='Total')
ax.set_xlabel('Temperature (K)', fontsize=12, color='white')
ax.set_ylabel(r'$C_V / Nk_B$', fontsize=12, color='white')
ax.set_title(r'$\mathrm{N_2}$: Heat Capacity Contributions', fontsize=14, color='white')
ax.legend(fontsize=9, facecolor='#1e293b', edgecolor='#6366f1', labelcolor='white')
ax.set_xlim(0, 6000)
ax.set_ylim(0, 4.5)
ax.axhline(y=2.5, color='gray', linestyle=':', alpha=0.5)
ax.axhline(y=3.5, color='gray', linestyle=':', alpha=0.5)
ax.set_facecolor('#0f172a')
ax.tick_params(colors='white')
ax.spines['bottom'].set_color('#6366f1')
ax.spines['left'].set_color('#6366f1')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.grid(True, alpha=0.15, color='#6366f1')

# H2O plot
ax = axes[1]
ax.plot(T, Cv_trans,     '--', color='#60a5fa', linewidth=1.5, label='Translational (3/2)')
ax.plot(T, Cv_rot_H2O,  '--', color='#34d399', linewidth=1.5, label='Rotational')
ax.plot(T, Cv_vib_H2O,  '--', color='#f87171', linewidth=1.5, label='Vibrational')
ax.plot(T, Cv_total_H2O, '-', color='white',   linewidth=2.5, label='Total')
ax.set_xlabel('Temperature (K)', fontsize=12, color='white')
ax.set_ylabel(r'$C_V / Nk_B$', fontsize=12, color='white')
ax.set_title(r'$\mathrm{H_2O}$: Heat Capacity Contributions', fontsize=14, color='white')
ax.legend(fontsize=9, facecolor='#1e293b', edgecolor='#6366f1', labelcolor='white')
ax.set_xlim(0, 6000)
ax.set_ylim(0, 7.0)
ax.axhline(y=3.0, color='gray', linestyle=':', alpha=0.5)
ax.axhline(y=6.0, color='gray', linestyle=':', alpha=0.5)
ax.set_facecolor('#0f172a')
ax.tick_params(colors='white')
ax.spines['bottom'].set_color('#6366f1')
ax.spines['left'].set_color('#6366f1')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.grid(True, alpha=0.15, color='#6366f1')

fig.patch.set_facecolor('#0f172a')
plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()
print("Heat capacity plot saved.")

# Print characteristic temperatures
print(f"\nN2 characteristic temperatures:")
print(f"  Theta_rot  = {Theta_rot_N2:.2f} K")
print(f"  Theta_vib  = {Theta_vib_N2:.1f} K")
print(f"  Theta_elec = {Theta_elec_N2:.0f} K")
print(f"\nH2O characteristic temperatures:")
print(f"  Theta_A = {Theta_A:.2f} K, Theta_B = {Theta_B:.2f} K, Theta_C = {Theta_C:.2f} K")
print(f"  Theta_vib = {Theta_vib_H2O[0]:.1f}, {Theta_vib_H2O[1]:.1f}, {Theta_vib_H2O[2]:.1f} K")
print(f"\nAt T = 300 K:")
print(f"  N2:  Cv/Nk_B = {np.interp(300, T, Cv_total_N2):.3f} (classical diatomic limit = 2.5)")
print(f"  H2O: Cv/Nk_B = {np.interp(300, T, Cv_total_H2O):.3f} (classical nonlinear limit = 3.0)")
print(f"\nAt T = 3000 K:")
print(f"  N2:  Cv/Nk_B = {np.interp(3000, T, Cv_total_N2):.3f}")
print(f"  H2O: Cv/Nk_B = {np.interp(3000, T, Cv_total_H2O):.3f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

10. Computation: Molecular Partition Functions (Fortran)

This Fortran program computes the translational, rotational, vibrational, and electronic partition functions for five diatomic molecules (N$_2$, O$_2$, CO, HCl, NO) at temperatures from 100 K to 5000 K. It uses exact summation for the rotational partition function at low temperatures and the high-temperature approximation when appropriate. Thermodynamic properties (energy, entropy, free energy) are computed at 298.15 K. Note the difference for NO, which has a doubly degenerate ground electronic state.

Diatomic Partition Functions & Thermodynamics

Fortran

Computes molecular partition functions (translational, rotational, vibrational, electronic) and thermodynamic properties for N₂, O₂, CO, HCl, and NO across a range of temperatures.

partition_functions.f90167 lines

program partition_functions
  implicit none

! ──────────────────────────────────────────────────────────────
  ! Molecular Partition Functions & Thermodynamic Properties
  ! for Diatomic Molecules: N2, O2, CO, HCl, NO
  ! ──────────────────────────────────────────────────────────────

integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: n_mol = 5
  integer, parameter :: n_temp = 8

real(dp), parameter :: k_B   = 1.380649e-23_dp    ! J/K
  real(dp), parameter :: h_bar = 1.054571817e-34_dp  ! J s
  real(dp), parameter :: h_pl  = 6.62607015e-34_dp   ! J s
  real(dp), parameter :: c_cm  = 2.99792458e10_dp    ! cm/s
  real(dp), parameter :: N_A   = 6.02214076e23_dp    ! 1/mol
  real(dp), parameter :: pi    = 3.14159265358979324_dp

character(len=4) :: names(n_mol)
  real(dp) :: mass(n_mol)        ! kg per molecule
  real(dp) :: B_rot(n_mol)       ! cm^-1
  real(dp) :: nu_vib(n_mol)      ! cm^-1
  integer  :: sigma(n_mol)       ! symmetry number
  integer  :: g_elec(n_mol)      ! ground state degeneracy
  real(dp) :: temps(n_temp)

real(dp) :: Theta_rot, Theta_vib
  real(dp) :: q_trans, q_rot, q_vib, q_elec, q_total
  real(dp) :: Lambda             ! thermal de Broglie wavelength
  real(dp) :: E_avg, Cv_total, S_total, A_helmholtz
  real(dp) :: x_vib, V_box
  real(dp) :: m_kg, T_curr
  integer  :: i, j

! Molecular data
  names = (/ 'N2  ', 'O2  ', 'CO  ', 'HCl ', 'NO  ' /)
  mass  = (/ 28.014_dp, 31.998_dp, 28.010_dp, 36.461_dp, 30.006_dp /) * 1.0e-3_dp / N_A

B_rot  = (/ 1.99824_dp, 1.43768_dp, 1.93128_dp, 10.5934_dp, 1.67195_dp /)
  nu_vib = (/ 2358.57_dp, 1580.19_dp, 2169.81_dp, 2990.95_dp, 1904.20_dp /)
  sigma  = (/ 2, 2, 1, 1, 1 /)
  g_elec = (/ 1, 3, 1, 1, 2 /)   ! NO has doubly degenerate ground state

temps = (/ 100.0_dp, 300.0_dp, 500.0_dp, 1000.0_dp, &
             2000.0_dp, 3000.0_dp, 4000.0_dp, 5000.0_dp /)

V_box = 1.0e-3_dp / N_A  ! Volume per molecule at ~1 mol/L (m^3)

! Header
  write(*,'(A)') '=================================================================='
  write(*,'(A)') '  Molecular Partition Functions & Thermodynamic Properties'
  write(*,'(A)') '  for Diatomic Molecules'
  write(*,'(A)') '=================================================================='
  write(*,*)

do i = 1, n_mol
    m_kg = mass(i)

write(*,'(A)') '------------------------------------------------------------------'
    write(*,'(A,A4)') '  Molecule: ', names(i)
    write(*,'(A,F10.4,A)') '  B_rot   = ', B_rot(i), ' cm^-1'
    write(*,'(A,F10.2,A)') '  nu_vib  = ', nu_vib(i), ' cm^-1'
    write(*,'(A,I4)')      '  sigma   = ', sigma(i)
    write(*,'(A,I4)')      '  g_elec  = ', g_elec(i)

Theta_rot = h_pl * c_cm * B_rot(i) / k_B
    Theta_vib = h_pl * c_cm * nu_vib(i) / k_B

write(*,'(A,F10.3,A)') '  Theta_rot = ', Theta_rot, ' K'
    write(*,'(A,F10.1,A)') '  Theta_vib = ', Theta_vib, ' K'
    write(*,'(A)') '------------------------------------------------------------------'
    write(*,'(A6, A14, A14, A14, A14, A14)') &
      'T(K)', 'q_trans', 'q_rot', 'q_vib', 'q_total', 'Cv/Nk_B'
    write(*,'(A)') '------------------------------------------------------------------'

do j = 1, n_temp
      T_curr = temps(j)

! Translational partition function
      Lambda = h_pl / sqrt(2.0_dp * pi * m_kg * k_B * T_curr)
      q_trans = V_box / Lambda**3

! Rotational partition function (high-T approximation)
      if (T_curr > 2.0_dp * Theta_rot) then
        q_rot = T_curr / (real(sigma(i), dp) * Theta_rot)
      else
        ! Exact sum for low T
        q_rot = 0.0_dp
        block
          integer :: J_val
          real(dp) :: E_J
          do J_val = 0, 300
            E_J = real(J_val * (J_val + 1), dp) * Theta_rot / T_curr
            if (E_J > 500.0_dp) exit
            q_rot = q_rot + real(2*J_val + 1, dp) * exp(-E_J)
          end do
        end block
        q_rot = q_rot / real(sigma(i), dp)
      end if

! Vibrational partition function
      x_vib = Theta_vib / T_curr
      if (x_vib < 500.0_dp) then
        q_vib = 1.0_dp / (1.0_dp - exp(-x_vib))
      else
        q_vib = 1.0_dp
      end if

! Electronic partition function
      q_elec = real(g_elec(i), dp)

! Total
      q_total = q_trans * q_rot * q_vib * q_elec

! Heat capacity Cv/Nk_B
      ! trans: 3/2, rot: 1 (linear), vib: x^2 e^x/(e^x - 1)^2
      Cv_total = 1.5_dp + 1.0_dp   ! trans + rot (high-T)
      if (x_vib < 500.0_dp) then
        Cv_total = Cv_total + x_vib**2 * exp(x_vib) / (exp(x_vib) - 1.0_dp)**2
      end if

write(*,'(F6.0, ES14.4, ES14.4, ES14.4, ES14.4, F14.4)') &
        T_curr, q_trans, q_rot, q_vib, q_total, Cv_total
    end do

write(*,*)

! Print thermodynamic properties at 298.15 K
    T_curr = 298.15_dp
    Lambda = h_pl / sqrt(2.0_dp * pi * m_kg * k_B * T_curr)
    q_trans = V_box / Lambda**3
    q_rot = T_curr / (real(sigma(i), dp) * Theta_rot)
    x_vib = Theta_vib / T_curr
    q_vib = 1.0_dp / (1.0_dp - exp(-x_vib))
    q_elec = real(g_elec(i), dp)
    q_total = q_trans * q_rot * q_vib * q_elec

! Thermodynamic properties per molecule (in units of k_B T)
    E_avg = (1.5_dp + 1.0_dp) * k_B * T_curr  ! trans + rot
    if (x_vib < 500.0_dp) then
      E_avg = E_avg + k_B * Theta_vib * (0.5_dp + 1.0_dp/(exp(x_vib) - 1.0_dp))
    end if

S_total = k_B * (log(q_total) + 2.5_dp + 1.0_dp)  ! Sackur-Tetrode + rot + vib
    A_helmholtz = -k_B * T_curr * log(q_total)

write(*,'(A,F8.2,A)') '  At T = 298.15 K:'
    write(*,'(A,ES14.4)')     '    Lambda (m)      = ', Lambda
    write(*,'(A,ES14.4)')     '    q_total         = ', q_total
    write(*,'(A,F14.6,A)')    '    <E>/k_BT        = ', E_avg / (k_B * T_curr)
    write(*,'(A,F14.6,A)')    '    S/k_B           = ', S_total / k_B
    write(*,'(A,F14.6,A)')    '    A/k_BT          = ', A_helmholtz / (k_B * T_curr)
    write(*,*)

end do

write(*,'(A)') '=================================================================='
  write(*,'(A)') '  Key observations:'
  write(*,'(A)') '  - q_trans >> q_rot >> q_vib at room temperature'
  write(*,'(A)') '  - Vibrational modes freeze out below Theta_vib'
  write(*,'(A)') '  - NO has g_elec=2 (doublet ground state Pi_{1/2})'
  write(*,'(A)') '  - HCl has large B_rot => low Theta_rot => early classical limit'
  write(*,'(A)') '=================================================================='

end program partition_functions

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Summary of Partition Functions

Mode	Partition Function	$C_V/k_B$ (classical)	Characteristic T
Translation	$V/\Lambda^3$	$3/2$	$\sim 10^{-15}$ K
Rotation (linear)	$T/(\sigma\Theta_{\text{rot}})$	$1$	$\sim 1$–$100$ K
Rotation (nonlinear)	$\frac{1}{\sigma}\sqrt{\pi T^3/(\Theta_A\Theta_B\Theta_C)}$	$3/2$	$\sim 10$–$50$ K
Vibration	$1/(1 - e^{-\Theta_{\text{vib}}/T})$	$1$ per mode	$\sim 1000$–$6000$ K
Electronic	$g_0 + g_1 e^{-\varepsilon_1/k_BT} + \cdots$	Schottky	$\sim 10^4$–$10^5$ K

Key Takeaways

The partition function $q$ is a sum over Boltzmann-weighted states that encodes all equilibrium thermodynamic information about a molecule.

Under the separability assumption, $q = q_{\text{trans}} \cdot q_{\text{rot}} \cdot q_{\text{vib}} \cdot q_{\text{elec}}$, allowing each contribution to be computed independently.

Characteristic temperatures ($\Theta_{\text{rot}}$, $\Theta_{\text{vib}}$) determine when each mode "thaws" — transitions from quantum freeze-out to classical equipartition.

The thermal de Broglie wavelength $\Lambda$ sets the quantum length scale for translation: when $\Lambda \ll$ the mean intermolecular spacing, classical behavior prevails.

All thermodynamic quantities follow from derivatives of $\ln Q$: energy, entropy, free energy, pressure, chemical potential, and equilibrium constants.

←Part IV Overview Back to Statistical Thermodynamics→