Part IV, Topic 4 | Lectures 25–28

Molecular Thermodynamics

Connecting molecular partition functions to macroscopic thermodynamic properties — entropy, free energy, heat capacity, and chemical equilibrium from first principles

1. Introduction: From Molecules to Bulk Thermodynamics

Statistical mechanics provides the fundamental bridge between the microscopic world of individual molecules — their energy levels, quantum states, and partition functions — and the macroscopic thermodynamic quantities we measure in the laboratory: entropy, enthalpy, free energy, and equilibrium constants.

The central quantity in this connection is the canonical partition function, which encodes all thermodynamic information about a system. For a system of $N$ independent, indistinguishable particles, the canonical partition function is:

$$Q(N,V,T) = \frac{q(V,T)^N}{N!}$$

where $q$ is the molecular partition function, and $N!$ corrects for indistinguishability

From $Q$, every thermodynamic function follows by differentiation. The Helmholtz free energy provides the most direct route:

$$A = -k_BT \ln Q$$

From this single expression we derive all other thermodynamic potentials:

$S = -\left(\frac{\partial A}{\partial T}\right)_V = k_B \ln Q + k_BT\left(\frac{\partial \ln Q}{\partial T}\right)_V$

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

$U = k_BT^2\left(\frac{\partial \ln Q}{\partial T}\right)_V$

$\mu = -k_BT\left(\frac{\partial \ln Q}{\partial N}\right)_{V,T}$

In this chapter, we apply these fundamental connections to derive explicit formulas for entropy, chemical equilibrium, heat capacities of solids, and residual entropy — all from molecular-level information.

2. Derivation 1: Entropy from the Partition Function

General Expression for Entropy

We begin from the Helmholtz free energy for $N$ indistinguishable particles:

$$A = -k_BT \ln Q = -k_BT \ln\left(\frac{q^N}{N!}\right) = -Nk_BT\ln q + k_BT\ln(N!)$$

Using Stirling's approximation $\ln(N!) \approx N\ln N - N$:

$$A = -Nk_BT\left[\ln\left(\frac{q}{N}\right) + 1\right]$$

The entropy is obtained from $S = -({\partial A}/{\partial T})_V$. Taking the derivative:

$$S = Nk_B\left[\ln\left(\frac{q}{N}\right) + 1\right] + Nk_BT\left(\frac{\partial \ln q}{\partial T}\right)_V$$

General entropy formula for $N$ indistinguishable particles

The Sackur-Tetrode Equation: Translational Entropy

For a monatomic ideal gas, the only contribution is translational. The translational partition function is:

$$q_{\text{trans}} = \frac{V}{\Lambda^3}, \qquad \Lambda = \frac{h}{\sqrt{2\pi mk_BT}}$$

where $\Lambda$ is the thermal de Broglie wavelength

We need $T(\partial \ln q/\partial T)_V$. Since $q_{\text{trans}} = V/\Lambda^3$ and $\Lambda \propto T^{-1/2}$:

$$\ln q_{\text{trans}} = \ln V - 3\ln\Lambda = \ln V + \frac{3}{2}\ln T + \text{const}$$

$$\therefore \quad T\left(\frac{\partial \ln q_{\text{trans}}}{\partial T}\right)_V = T \cdot \frac{3}{2T} = \frac{3}{2}$$

Substituting into our general entropy formula:

$$S_{\text{trans}} = Nk_B\left[\ln\left(\frac{q_{\text{trans}}}{N}\right) + 1 + \frac{3}{2}\right]$$

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

The Sackur-Tetrode Equation

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

where $\Lambda = h/\sqrt{2\pi mk_BT}$

For an ideal gas at pressure $P$, using $V/N = k_BT/P$:

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

Numerical Example: Argon at 298.15 K

For Ar ($m = 39.948$ u = $6.634 \times 10^{-26}$ kg) at 298.15 K, 1 atm:

$\Lambda = \frac{6.626 \times 10^{-34}}{\sqrt{2\pi(6.634 \times 10^{-26})(1.381 \times 10^{-23})(298.15)}} = 1.60 \times 10^{-11}$ m

$V/N = k_BT/P = (1.381 \times 10^{-23})(298.15)/(101325) = 4.065 \times 10^{-26}$ m$^3$

$S = Nk_B[5/2 + \ln(4.065 \times 10^{-26}/(1.60 \times 10^{-11})^3)]$

$= R[5/2 + \ln(9.92 \times 10^6)] = R[5/2 + 16.11] = R \times 18.61$

$S_{\text{trans}}^{\circ}(\text{Ar}) = 154.8$ J/(mol·K) (Experimental: 154.8 J/(mol·K))

Complete Molecular Entropy

For a polyatomic molecule, the total entropy includes all degrees of freedom. With factored partition functions$q = q_{\text{trans}} \cdot q_{\text{rot}} \cdot q_{\text{vib}} \cdot q_{\text{elec}}$:

$$S = S_{\text{trans}} + S_{\text{rot}} + S_{\text{vib}} + S_{\text{elec}}$$

Rotational (linear molecule, high-T limit): $S_{\text{rot}} = Nk_B\left[1 + \ln\left(\frac{T}{\sigma \Theta_{\text{rot}}}\right)\right]$

Vibrational (per mode): $S_{\text{vib}} = Nk_B\left[\frac{\Theta_{\text{vib}}/T}{e^{\Theta_{\text{vib}}/T}-1} - \ln\left(1 - e^{-\Theta_{\text{vib}}/T}\right)\right]$

Electronic: $S_{\text{elec}} = Nk_B\ln g_0$ (when only the ground state is occupied)

3. Derivation 2: Chemical Equilibrium from Partition Functions

Deriving the Equilibrium Constant

Consider a general gas-phase reaction $\sum_i \nu_i A_i = 0$ (with $\nu_i > 0$ for products,$\nu_i < 0$ for reactants). At equilibrium, the Gibbs free energy is minimized, which requires $\sum_i \nu_i \mu_i = 0$.

The chemical potential for species $i$ in an ideal gas mixture is obtained from the partition function. Starting from $A = -k_BT\ln Q$ for species $i$:

$$\mu_i = -k_BT\ln\left(\frac{q_i}{N_i}\right)$$

It is conventional to measure energies relative to the separated atoms at rest. Define $q_i = q_i^{\prime} \cdot e^{-\epsilon_{0,i}/k_BT}$ where $q_i^{\prime}$ is the partition function with the ground state energy set to zero, and $\epsilon_{0,i}$ is the zero-point energy of species $i$. Then:

$$\mu_i = -k_BT\ln\left(\frac{q_i^{\prime}}{N_i}\right) + \epsilon_{0,i}$$

The equilibrium condition $\sum_i \nu_i \mu_i = 0$ becomes:

$$\sum_i \nu_i \left[-k_BT\ln\left(\frac{q_i^{\prime}}{N_i}\right) + \epsilon_{0,i}\right] = 0$$

$$\sum_i \nu_i \ln\left(\frac{q_i^{\prime}}{N_i}\right) = \frac{1}{k_BT}\sum_i \nu_i \epsilon_{0,i} = \frac{\Delta E_0}{k_BT}$$

Writing the concentration $N_i/V = c_i$ and separating the volume dependence of the translational partition function via $q_i^{\prime} = (q_i^{\prime}/V) \cdot V$:

Statistical Mechanical Equilibrium Constant

$$\boxed{K = \prod_i \left(\frac{q_i^{\prime}}{V}\right)^{\nu_i} \cdot e^{-\Delta E_0 / k_BT}}$$

where $\Delta E_0 = \sum_i \nu_i \epsilon_{0,i}$ is the difference in zero-point energies between products and reactants

Application: Dissociation of I$_2$

Consider the dissociation equilibrium $\text{I}_2(g) \rightleftharpoons 2\text{I}(g)$. The equilibrium constant is:

$$K_c = \frac{(q_{\text{I}}^{\prime}/V)^2}{q_{\text{I}_2}^{\prime}/V} \cdot e^{-\Delta E_0/k_BT}$$

For the iodine atom: $q_{\text{I}}^{\prime}/V = g_{\text{elec}} \cdot \Lambda_{\text{I}}^{-3}$, where $g_{\text{elec}} = 4$ (ground state $^2P_{3/2}$).

For the I$_2$ molecule: $q_{\text{I}_2}^{\prime}/V = \Lambda_{\text{I}_2}^{-3} \cdot q_{\text{rot}} \cdot q_{\text{vib}} \cdot g_{\text{elec}}$, with spectroscopic constants $\tilde{\nu} = 214.5$ cm$^{-1}$, $B_e = 0.03737$ cm$^{-1}$, $\sigma = 2$, and $D_0 = 1.5422$ eV.

At T = 1000 K:

$\Theta_{\text{vib}} = hc\tilde{\nu}/k_B = 308.7$ K, $\Theta_{\text{rot}} = hcB_e/k_B = 0.0538$ K

$q_{\text{rot}} = T/(2\Theta_{\text{rot}}) = 9294$

$q_{\text{vib}} = 1/(1 - e^{-308.7/1000}) = 3.762$

$\Delta E_0/k_BT = D_0/k_BT = 1.5422 \times 1.602 \times 10^{-19}/(1.381 \times 10^{-23} \times 1000) = 17.89$

$K_c \approx 3.1 \times 10^{-5}$ mol/L — very small, as expected (I$_2$ is mostly undissociated at 1000 K)

4. Derivation 3: Einstein and Debye Models of Solids

The Einstein Model (1907)

Einstein proposed the first quantum model of heat capacity by treating each atom in a crystal as an independent quantum harmonic oscillator, all vibrating at the same frequency $\nu_E$. For $N$ atoms in three dimensions, there are $3N$ independent oscillators.

The partition function for a single harmonic oscillator is:

$$q_{\text{osc}} = \frac{e^{-\Theta_E/2T}}{1 - e^{-\Theta_E/T}}, \qquad \Theta_E = \frac{h\nu_E}{k_B}$$

$\Theta_E$ is the Einstein temperature characteristic of the material

The average energy of one oscillator (above the zero-point energy) is:

$$\langle \epsilon \rangle = k_BT^2 \frac{\partial \ln q_{\text{osc}}}{\partial T} = \frac{k_B\Theta_E}{e^{\Theta_E/T} - 1}$$

The total internal energy for $3N$ oscillators is $U = 3N\langle\epsilon\rangle$, and the heat capacity follows from $C_V = (\partial U/\partial T)_V$:

$$C_V = 3Nk_B \frac{\partial}{\partial T}\left[\frac{\Theta_E}{e^{\Theta_E/T} - 1}\right]$$

$$= 3Nk_B \left(\frac{\Theta_E}{T}\right)^2 \frac{e^{\Theta_E/T}}{\left(e^{\Theta_E/T} - 1\right)^2}$$

Einstein Heat Capacity

$$\boxed{C_V^{\text{Einstein}} = 3Nk_B\left(\frac{\Theta_E}{T}\right)^2 \frac{e^{\Theta_E/T}}{\left(e^{\Theta_E/T} - 1\right)^2}}$$

High-T Limit ($T \gg \Theta_E$)

When $\Theta_E/T \to 0$: $e^{\Theta_E/T} \approx 1 + \Theta_E/T$

$C_V \to 3Nk_B$ (Dulong–Petit law, the classical result)

Low-T Limit ($T \ll \Theta_E$)

When $\Theta_E/T \to \infty$: $C_V \to 3Nk_B(\Theta_E/T)^2 e^{-\Theta_E/T}$

Exponential decay — too fast compared to experiment ($C_V \propto T^3$)

The Debye Model (1912)

Debye improved upon Einstein's model by recognizing that atoms in a crystal do not vibrate independently — they are coupled, forming collective vibrational modes (phonons) with a continuous spectrum of frequencies from 0 up to a maximum cutoff frequency $\omega_D$.

The Debye approximation treats the crystal as an elastic continuum with a density of states:

$$g(\omega) = \frac{9N}{\omega_D^3}\omega^2 \qquad \text{for } 0 \leq \omega \leq \omega_D$$

The $\omega^2$ dependence comes from the density of states of an elastic continuum in 3D

The total energy is obtained by integrating over all modes, each weighted by the Planck distribution:

$$U = \int_0^{\omega_D} \frac{\hbar\omega}{e^{\hbar\omega/k_BT} - 1} g(\omega)\,d\omega = \frac{9N\hbar}{\omega_D^3}\int_0^{\omega_D}\frac{\omega^3}{e^{\hbar\omega/k_BT} - 1}\,d\omega$$

Substituting $x = \hbar\omega/k_BT$ and defining $\Theta_D = \hbar\omega_D/k_B$:

$$U = 9Nk_BT\left(\frac{T}{\Theta_D}\right)^3 \int_0^{\Theta_D/T} \frac{x^3}{e^x - 1}\,dx$$

Taking $C_V = \partial U/\partial T$ and differentiating under the integral sign:

Debye Heat Capacity

$$\boxed{C_V^{\text{Debye}} = 9Nk_B\left(\frac{T}{\Theta_D}\right)^3 \int_0^{\Theta_D/T} \frac{x^4 e^x}{(e^x - 1)^2}\,dx}$$

Low-Temperature Limit: The Debye $T^3$ Law

At very low temperatures ($T \ll \Theta_D$), the upper limit $\Theta_D/T \to \infty$, and the integral becomes a standard definite integral:

$$\int_0^{\infty} \frac{x^4 e^x}{(e^x - 1)^2}\,dx = \frac{4\pi^4}{15}$$

$$\therefore \quad C_V \xrightarrow{T \ll \Theta_D} \frac{12\pi^4}{5}Nk_B\left(\frac{T}{\Theta_D}\right)^3 \propto T^3$$

This $T^3$ dependence is in excellent agreement with experiments at low temperatures — a major triumph of the Debye model.

Feature	Einstein Model	Debye Model
Frequency spectrum	Single frequency $\nu_E$	Continuous: $g(\omega) \propto \omega^2$
High-T limit	$3Nk_B$ (correct)	$3Nk_B$ (correct)
Low-T behavior	$\sim e^{-\Theta_E/T}$ (too fast)	$\sim T^3$ (correct)
Free parameter	$\Theta_E$	$\Theta_D$
Physical picture	Independent oscillators	Coupled oscillators (phonons)

5. Derivation 4: Residual Entropy and the Third Law

The Third Law of Thermodynamics

The third law states that the entropy of a perfect crystal approaches zero as $T \to 0$:

$$\lim_{T \to 0} S = 0 \qquad \text{(perfect crystal)}$$

Statistical interpretation: at $T = 0$, only the ground state is occupied, so $W = 1$ and $S = k_B\ln 1 = 0$

Residual Entropy: When $S_0 > 0$

Some substances retain a nonzero entropy at $T = 0$ — this is called residual entropy. It arises when molecular disorder is frozen in at low temperatures, meaning the system cannot reach its unique ground state on experimental timescales.

The residual entropy is calculated from the Boltzmann formula:

Residual Entropy

$$\boxed{S_0 = k_B \ln W}$$

where $W$ is the number of energetically equivalent orientational configurations at $T = 0$

Classic Examples of Residual Entropy

Carbon Monoxide (CO)

CO and OC are nearly isoenergetic due to the similar sizes of C and O. In the crystal, each molecule can point in either direction with nearly equal probability.

For $N$ molecules, each with 2 orientations: $W = 2^N$

$S_0 = k_B\ln(2^N) = Nk_B\ln 2 = R\ln 2 = 5.76$ J/(mol·K)

Experimental: 4.6 J/(mol·K) — close to $R\ln 2$

Nitrous Oxide (N$_2$O)

N$_2$O is a linear molecule (NNO) that can orient as NNO or ONN in the crystal. Like CO, the two orientations are nearly degenerate.

$W = 2^N$, so $S_0 = R\ln 2 = 5.76$ J/(mol·K)

Experimental: 5.2 J/(mol·K)

Ice (H$_2$O)

In ice Ih, each oxygen is tetrahedrally coordinated with four neighbors via hydrogen bonds. Each O–H···O bond has the H atom closer to one oxygen (covalent) or the other (hydrogen bond). The "ice rules" (Bernal–Fowler rules) require two H atoms close to each O, but many configurations satisfy this.

Pauling's calculation: $W = (3/2)^N$

$S_0 = R\ln(3/2) = 3.37$ J/(mol·K)

Experimental: 3.4 J/(mol·K) — excellent agreement with Pauling's estimate

Conditions for Residual Entropy

Residual entropy occurs when:

• Multiple molecular orientations are nearly degenerate in energy ($\Delta E \ll k_BT$ at the freezing point)
• The reorientation barrier is high, preventing the system from annealing to the ordered ground state
• The disorder involves local configurations that are constrained but not uniquely determined by the crystal structure

Residual entropy is a direct macroscopic manifestation of microscopic disorder. It provides a powerful experimental test of the statistical mechanical interpretation of entropy, and the agreement between$S_0 = k_B\ln W$ and calorimetric measurements is one of the great confirmations of Boltzmann's entropy formula.

6. Applications of Molecular Thermodynamics

JANAF Thermochemical Tables

The Joint Army-Navy-Air Force (JANAF) tables compile standard thermodynamic functions ($S^{\circ}$,$H^{\circ} - H_0^{\circ}$, $\Delta_f H^{\circ}$, $\Delta_f G^{\circ}$) computed largely from spectroscopic data via the partition function approach developed in this chapter.

These tables are the foundation of combustion chemistry, rocket propulsion calculations, and high-temperature chemical engineering.

Computational Thermochemistry

Modern quantum chemistry programs (Gaussian, ORCA, etc.) routinely compute thermodynamic properties using the partition function formalism. From a frequency calculation on an optimized geometry, the software computes:

• $S^{\circ}$ from Sackur-Tetrode + rigid rotor + harmonic oscillator
• $C_p^{\circ}$ and $H^{\circ}$ from energy derivatives
• $G^{\circ}$ and equilibrium constants from $\Delta G^{\circ} = -RT\ln K$

Molecular Design of Materials

Understanding how molecular-level properties determine bulk thermodynamics enables rational material design:

• Low Debye temperature materials for thermoelectric applications
• High-entropy alloys with controlled configurational entropy
• Phase-change materials where entropy drives the transition

Atmospheric Chemistry

Equilibrium constants for high-temperature atmospheric reactions (combustion, reentry, lightning) are calculated from partition functions when direct measurements are unavailable:

• N$_2$ + O$_2$ $\rightleftharpoons$ 2NO (Zeldovich mechanism)
• Dissociation equilibria of diatomics at extreme temperatures
• Formation of exotic species in shock waves and plasmas

7. Historical Context

The development of molecular thermodynamics was driven by critical failures of classical physics and the need to reconcile the atomic hypothesis with macroscopic observations.

1907

Albert Einstein — Quantum Theory of Heat Capacity

Einstein proposed that each atom in a solid vibrates as a quantum harmonic oscillator, explaining the decrease of $C_V$ below the Dulong–Petit value at low temperatures. This was the first application of quantum theory to solid-state physics.

1912

Peter Debye — Phonon Model of Solids

Debye improved the Einstein model by introducing a continuous spectrum of vibrational frequencies, correctly predicting the $T^3$ law for heat capacity at low temperatures. His model treated the solid as an elastic continuum with a frequency cutoff.

1912

Otto Sackur & Hugo Tetrode — Translational Entropy

Working independently, Sackur and Tetrode derived the quantum expression for the translational entropy of an ideal gas. Their formula, which gives absolute entropy values in perfect agreement with experiment, was one of the first quantitative successes of quantum statistical mechanics.

1930s

William F. Giauque — Third Law and Residual Entropy

Giauque performed precise calorimetric measurements that revealed discrepancies between spectroscopic and calorimetric entropies for certain substances (CO, N$_2$O, H$_2$O ice). These "residual entropies" were explained by frozen-in molecular disorder. Giauque received the 1949 Nobel Prize in Chemistry for his work on chemical thermodynamics at very low temperatures.

1935

Linus Pauling — Residual Entropy of Ice

Pauling calculated the residual entropy of ice using a combinatorial argument based on the ice rules, obtaining $S_0 = R\ln(3/2)$ in remarkable agreement with Giauque's experimental value. This calculation remains a classic example of statistical mechanics applied to disordered systems.

8. Python Simulation: Einstein & Debye Heat Capacities

This simulation computes and plots the Einstein and Debye heat capacity models for diamond and aluminum, comparing both models with experimental data. The third panel focuses on the low-temperature regime where the Debye $T^3$ law dominates and the Einstein model fails.

Python

script.py165 lines

import numpy as np
import matplotlib.pyplot as plt

# ================================================================
# Einstein and Debye Heat Capacity Models vs Experimental Data
# ================================================================

fig, axes = plt.subplots(1, 3, figsize=(18, 6))
fig.patch.set_facecolor('#0f172a')

# ---- Model Functions ----

def einstein_cv(T, theta_E):
    """Einstein model: C_V / 3Nk = (theta_E/T)^2 * exp(theta_E/T) / (exp(theta_E/T) - 1)^2"""
    x = theta_E / T
    # Avoid overflow for very large x
    mask = x < 500
    result = np.zeros_like(T, dtype=float)
    result[mask] = x[mask]**2 * np.exp(x[mask]) / (np.exp(x[mask]) - 1.0)**2
    return result

def debye_cv(T, theta_D):
    """Debye model: C_V / 3Nk = 9(T/theta_D)^3 * integral_0^{theta_D/T} x^4 e^x / (e^x - 1)^2 dx"""
    results = np.zeros_like(T, dtype=float)
    for i, temp in enumerate(T):
        if temp < 1e-10:
            results[i] = 0.0
            continue
        x_max = theta_D / temp
        n_pts = 2000
        x = np.linspace(1e-10, x_max, n_pts)
        # Integrand: x^4 * e^x / (e^x - 1)^2
        mask = x < 500
        integrand = np.zeros_like(x)
        integrand[mask] = x[mask]**4 * np.exp(x[mask]) / (np.exp(x[mask]) - 1.0)**2
        integral = np.trapezoid(integrand, x)
        results[i] = 9.0 * (temp / theta_D)**3 * integral
    return results

# ---- Material Parameters ----
# Theta_E and Theta_D values (in K)
materials = {
    'Diamond': {
        'theta_E': 1320,
        'theta_D': 2230,
        'color_exp': '#f472b6',
        'color_E': '#fb923c',
        'color_D': '#a78bfa',
        # Experimental C_V data for diamond (C_V / 3R) at various T
        'T_exp': np.array([100, 200, 300, 400, 500, 600, 800, 1000, 1200, 1500]),
        'cv_exp': np.array([0.006, 0.085, 0.252, 0.428, 0.575, 0.683, 0.822, 0.899, 0.941, 0.972]),
    },
    'Aluminum': {
        'theta_E': 290,
        'theta_D': 428,
        'color_exp': '#4ade80',
        'color_E': '#facc15',
        'color_D': '#38bdf8',
        # Experimental C_V data for aluminum (C_V / 3R)
        'T_exp': np.array([20, 40, 80, 100, 150, 200, 250, 300, 400, 500]),
        'cv_exp': np.array([0.009, 0.078, 0.445, 0.585, 0.790, 0.885, 0.935, 0.960, 0.985, 0.995]),
    },
}

# ---- Panel 1: Diamond ----
ax1 = axes[0]
ax1.set_facecolor('#1e293b')
mat = materials['Diamond']
T = np.linspace(1, 2000, 500)
cv_einstein = einstein_cv(T, mat['theta_E'])
cv_debye = debye_cv(T, mat['theta_D'])

ax1.plot(T, cv_einstein, '-', color=mat['color_E'], linewidth=2, label=f"Einstein (\u0398_E={mat['theta_E']} K)")
ax1.plot(T, cv_debye, '--', color=mat['color_D'], linewidth=2, label=f"Debye (\u0398_D={mat['theta_D']} K)")
ax1.scatter(mat['T_exp'], mat['cv_exp'], color=mat['color_exp'], s=60, zorder=5,
            edgecolors='white', linewidths=0.5, label='Experimental')
ax1.set_xlabel('Temperature (K)', color='#cbd5e1', fontsize=11)
ax1.set_ylabel('C_V / 3Nk', color='#cbd5e1', fontsize=11)
ax1.set_title('Diamond', fontsize=14, fontweight='bold', color='#e2e8f0')
ax1.legend(fontsize=8, facecolor='#1e293b', edgecolor='#475569', labelcolor='#cbd5e1')
ax1.set_xlim(0, 2000)
ax1.set_ylim(0, 1.1)
ax1.tick_params(colors='#94a3b8')
for sp in ax1.spines.values():
    sp.set_color('#475569')
ax1.grid(True, alpha=0.15, color='#64748b')

# ---- Panel 2: Aluminum ----
ax2 = axes[1]
ax2.set_facecolor('#1e293b')
mat = materials['Aluminum']
T2 = np.linspace(1, 600, 500)
cv_einstein2 = einstein_cv(T2, mat['theta_E'])
cv_debye2 = debye_cv(T2, mat['theta_D'])

ax2.plot(T2, cv_einstein2, '-', color=mat['color_E'], linewidth=2, label=f"Einstein (\u0398_E={mat['theta_E']} K)")
ax2.plot(T2, cv_debye2, '--', color=mat['color_D'], linewidth=2, label=f"Debye (\u0398_D={mat['theta_D']} K)")
ax2.scatter(mat['T_exp'], mat['cv_exp'], color=mat['color_exp'], s=60, zorder=5,
            edgecolors='white', linewidths=0.5, label='Experimental')
ax2.set_xlabel('Temperature (K)', color='#cbd5e1', fontsize=11)
ax2.set_ylabel('C_V / 3Nk', color='#cbd5e1', fontsize=11)
ax2.set_title('Aluminum', fontsize=14, fontweight='bold', color='#e2e8f0')
ax2.legend(fontsize=8, facecolor='#1e293b', edgecolor='#475569', labelcolor='#cbd5e1')
ax2.set_xlim(0, 600)
ax2.set_ylim(0, 1.1)
ax2.tick_params(colors='#94a3b8')
for sp in ax2.spines.values():
    sp.set_color('#475569')
ax2.grid(True, alpha=0.15, color='#64748b')

# ---- Panel 3: Low-T Comparison (T^3 law) ----
ax3 = axes[2]
ax3.set_facecolor('#1e293b')
T_low = np.linspace(1, 100, 300)

# Diamond low-T
cv_debye_dia = debye_cv(T_low, 2230)
cv_einstein_dia = einstein_cv(T_low, 1320)

# T^3 law: C_V/3Nk = (4*pi^4/5) * (T/theta_D)^3
t3_dia = (4 * np.pi**4 / 5) * (T_low / 2230)**3

ax3.plot(T_low, cv_debye_dia, '-', color='#a78bfa', linewidth=2, label='Debye (Diamond)')
ax3.plot(T_low, cv_einstein_dia, '--', color='#fb923c', linewidth=2, label='Einstein (Diamond)')
ax3.plot(T_low, t3_dia, ':', color='#f472b6', linewidth=2, label='T\u00b3 law (Diamond)')

ax3.set_xlabel('Temperature (K)', color='#cbd5e1', fontsize=11)
ax3.set_ylabel('C_V / 3Nk', color='#cbd5e1', fontsize=11)
ax3.set_title('Low-T Behavior', fontsize=14, fontweight='bold', color='#e2e8f0')
ax3.legend(fontsize=8, facecolor='#1e293b', edgecolor='#475569', labelcolor='#cbd5e1')
ax3.set_xlim(0, 100)
ax3.set_ylim(0, 0.02)
ax3.tick_params(colors='#94a3b8')
for sp in ax3.spines.values():
    sp.set_color('#475569')
ax3.grid(True, alpha=0.15, color='#64748b')

plt.suptitle('Einstein vs Debye Heat Capacity Models',
             fontsize=16, fontweight='bold', color='#e2e8f0', y=1.02)
plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()

# Print summary table
print("Einstein and Debye Heat Capacity Models")
print("=" * 65)
print(f"{'Material':<12} {'Theta_E (K)':<14} {'Theta_D (K)':<14} {'C_V/3Nk at 300K'}")
print("-" * 65)
for name, m in materials.items():
    T_pt = np.array([300.0])
    cv_e = einstein_cv(T_pt, m['theta_E'])[0]
    cv_d = debye_cv(T_pt, m['theta_D'])[0]
    print(f"{name:<12} {m['theta_E']:<14} {m['theta_D']:<14} E={cv_e:.4f}  D={cv_d:.4f}")

print()
print("Key Observations:")
print("-" * 65)
print("1. Einstein model underestimates C_V at low T (exponential decay)")
print("2. Debye model correctly gives C_V ~ T^3 at low T")
print("3. Both models approach Dulong-Petit limit (C_V = 3Nk) at high T")
print("4. Diamond has very high Theta_D => classical limit reached only at ~2000 K")
print("5. Aluminum has lower Theta_D => reaches classical limit by ~400 K")
print()
print("Low-T Debye limit: C_V = (12/5)*pi^4*Nk*(T/Theta_D)^3")
print("High-T limit: C_V -> 3Nk (Dulong-Petit law)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

9. Fortran Simulation: Thermodynamic Functions from Spectroscopic Data

This Fortran program computes the standard molar entropy ($S^{\circ}$), integrated enthalpy ($H^{\circ} - H_0^{\circ}$), and their individual contributions (translational, rotational, vibrational, electronic) for diatomic molecules using spectroscopic constants. The Sackur-Tetrode equation is used for translational entropy, the rigid rotor approximation for rotational entropy, and the harmonic oscillator model for vibrational contributions.

Fortran

program.f90169 lines

program thermodynamic_functions
  implicit none
  ! Compute standard thermodynamic functions from spectroscopic data
  ! S^o, (H^o - H_0^o), and (G^o - H_0^o) for diatomic molecules

integer, parameter :: n_temps = 12
  double precision :: T_arr(n_temps)
  double precision :: R, kB, h_planck, c_light, Na
  double precision :: omega_e, omega_ex, B_e, sigma_sym, g0
  double precision :: D_0, mass1, mass2, mu_red
  double precision :: T, beta, q_trans, q_rot, q_vib, q_elec, q_total
  double precision :: S_trans, S_rot, S_vib, S_elec, S_total
  double precision :: H_trans, H_rot, H_vib, H_total
  double precision :: G_total
  double precision :: theta_rot, theta_vib, pi_val
  double precision :: x_vib, exp_vib
  integer :: i
  character(len=10) :: mol_name

R = 8.314472d0
  kB = 1.380649d-23
  h_planck = 6.62607015d-34
  c_light = 2.99792458d10
  Na = 6.02214076d23
  pi_val = 3.14159265358979d0

T_arr = (/200.0d0, 298.15d0, 400.0d0, 500.0d0, 600.0d0, &
            800.0d0, 1000.0d0, 1500.0d0, 2000.0d0, 3000.0d0, &
            4000.0d0, 5000.0d0/)

write(*,'(A)') '======================================================================'
  write(*,'(A)') '  Standard Thermodynamic Functions from Spectroscopic Data'
  write(*,'(A)') '======================================================================'

! ---- HCl molecule ----
  mol_name = 'HCl'
  omega_e = 2990.95d0    ! cm^-1, fundamental frequency
  omega_ex = 52.82d0     ! cm^-1, anharmonicity
  B_e = 10.5934d0        ! cm^-1, rotational constant
  sigma_sym = 1           ! symmetry number
  g0 = 1.0d0             ! electronic degeneracy
  mass1 = 1.00794d0      ! H mass in amu
  mass2 = 35.4527d0      ! Cl mass in amu
  mu_red = (mass1 * mass2) / (mass1 + mass2) * 1.66054d-27

theta_rot = h_planck * c_light * B_e / kB
  theta_vib = h_planck * c_light * omega_e / kB

write(*,*)
  write(*,'(A,A)') '  Molecule: ', mol_name
  write(*,'(A,F10.2,A)') '  omega_e  = ', omega_e, ' cm^-1'
  write(*,'(A,F10.2,A)') '  B_e      = ', B_e, ' cm^-1'
  write(*,'(A,F10.2,A)') '  Theta_rot= ', theta_rot, ' K'
  write(*,'(A,F10.2,A)') '  Theta_vib= ', theta_vib, ' K'
  write(*,'(A,I5)')      '  sigma    = ', sigma_sym
  write(*,*)
  write(*,'(A8,A14,A14,A14,A14,A14)') &
       'T(K)', 'S_trans', 'S_rot', 'S_vib', 'S_total', 'H-H0 (kJ)'
  write(*,'(A)') '----------------------------------------------------------------------'

do i = 1, n_temps
    T = T_arr(i)
    beta = 1.0d0 / (kB * T)

! Translational entropy: Sackur-Tetrode
    ! S_trans/R = 5/2 + ln[ (2*pi*m*kT/h^2)^(3/2) * kT/P ]
    S_trans = R * (2.5d0 + 1.5d0 * log(2.0d0 * pi_val * mu_red * kB * T &
              / (h_planck**2)) + log(kB * T / 101325.0d0))

! Rotational entropy (high-T limit): S_rot/R = 1 + ln(T/(sigma*theta_rot))
    S_rot = R * (1.0d0 + log(T / (dble(sigma_sym) * theta_rot)))

! Vibrational entropy
    x_vib = theta_vib / T
    exp_vib = exp(x_vib)
    if (x_vib < 500.0d0) then
      S_vib = R * (x_vib / (exp_vib - 1.0d0) - log(1.0d0 - 1.0d0/exp_vib))
    else
      S_vib = 0.0d0
    end if

! Electronic entropy
    S_elec = R * log(g0)

S_total = S_trans + S_rot + S_vib + S_elec

! Enthalpy contributions (H - H_0)
    H_trans = 2.5d0 * R * T
    H_rot = R * T
    if (x_vib < 500.0d0) then
      H_vib = R * theta_vib / (exp_vib - 1.0d0)
    else
      H_vib = 0.0d0
    end if
    H_total = (H_trans + H_rot + H_vib) / 1000.0d0  ! kJ/mol

write(*,'(F8.1, F14.3, F14.3, F14.3, F14.3, F14.3)') &
         T, S_trans, S_rot, S_vib, S_total, H_total
  end do

! ---- N2 molecule ----
  write(*,*)
  write(*,'(A)') '----------------------------------------------------------------------'
  mol_name = 'N2'
  omega_e = 2358.57d0
  omega_ex = 14.32d0
  B_e = 1.99824d0
  sigma_sym = 2
  g0 = 1.0d0
  mass1 = 14.0067d0
  mass2 = 14.0067d0
  mu_red = (mass1 * mass2) / (mass1 + mass2) * 1.66054d-27

theta_rot = h_planck * c_light * B_e / kB
  theta_vib = h_planck * c_light * omega_e / kB

write(*,'(A,A)') '  Molecule: ', mol_name
  write(*,'(A,F10.2,A)') '  omega_e  = ', omega_e, ' cm^-1'
  write(*,'(A,F10.2,A)') '  B_e      = ', B_e, ' cm^-1'
  write(*,'(A,F10.2,A)') '  Theta_rot= ', theta_rot, ' K'
  write(*,'(A,F10.2,A)') '  Theta_vib= ', theta_vib, ' K'
  write(*,'(A,I5)')      '  sigma    = ', sigma_sym
  write(*,*)
  write(*,'(A8,A14,A14,A14,A14,A14)') &
       'T(K)', 'S_trans', 'S_rot', 'S_vib', 'S_total', 'H-H0 (kJ)'
  write(*,'(A)') '----------------------------------------------------------------------'

do i = 1, n_temps
    T = T_arr(i)
    beta = 1.0d0 / (kB * T)

S_trans = R * (2.5d0 + 1.5d0 * log(2.0d0 * pi_val * mu_red * kB * T &
              / (h_planck**2)) + log(kB * T / 101325.0d0))
    S_rot = R * (1.0d0 + log(T / (dble(sigma_sym) * theta_rot)))

x_vib = theta_vib / T
    exp_vib = exp(x_vib)
    if (x_vib < 500.0d0) then
      S_vib = R * (x_vib / (exp_vib - 1.0d0) - log(1.0d0 - 1.0d0/exp_vib))
    else
      S_vib = 0.0d0
    end if
    S_elec = R * log(g0)
    S_total = S_trans + S_rot + S_vib + S_elec

H_trans = 2.5d0 * R * T
    H_rot = R * T
    if (x_vib < 500.0d0) then
      H_vib = R * theta_vib / (exp_vib - 1.0d0)
    else
      H_vib = 0.0d0
    end if
    H_total = (H_trans + H_rot + H_vib) / 1000.0d0

write(*,'(F8.1, F14.3, F14.3, F14.3, F14.3, F14.3)') &
         T, S_trans, S_rot, S_vib, S_total, H_total
  end do

write(*,*)
  write(*,'(A)') '======================================================================'
  write(*,'(A)') '  Notes on computation:'
  write(*,'(A)') '  - S_trans from Sackur-Tetrode equation at P = 1 atm'
  write(*,'(A)') '  - S_rot from high-T rigid rotor (T >> Theta_rot)'
  write(*,'(A)') '  - S_vib from harmonic oscillator partition function'
  write(*,'(A)') '  - All units: S in J/(mol*K), H-H0 in kJ/mol'
  write(*,'(A)') '  - Reference: NIST-JANAF Thermochemical Tables'
  write(*,'(A)') '======================================================================'

end program thermodynamic_functions

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Back to Part IV: Statistical Thermodynamics Continue to Part V →