Part II: Thermodynamics | Chapter 3

Phase Equilibria

A rigorous treatment of phase transitions, coexistence curves, and the thermodynamic conditions governing equilibrium between phases of matter

1. Introduction: Phases, Phase Transitions, and Equilibrium

A phase is a macroscopically homogeneous region of matter with uniform chemical composition and physical properties. The three familiar phases — solid, liquid, and gas — are distinguished by their molecular-level organization: solids possess long-range translational order, liquids have short-range order only, and gases are essentially disordered with molecules separated by distances much larger than their size.

A phase transition occurs when a system changes from one phase to another in response to changes in thermodynamic variables such as temperature, pressure, or composition. At the macroscopic level, phase transitions are governed by thermodynamics: the stable phase at any given temperature and pressure is the one that minimizes the Gibbs free energy $G(T,P)$.

Conditions for Phase Equilibrium

When two phases $\alpha$ and $\beta$ coexist in equilibrium, three conditions must be satisfied simultaneously:

$$T^\alpha = T^\beta \quad \text{(thermal equilibrium)}$$

$$P^\alpha = P^\beta \quad \text{(mechanical equilibrium)}$$

$$\mu^\alpha = \mu^\beta \quad \text{(chemical equilibrium)}$$

The last condition — equality of chemical potentials — is the most powerful, as it determines which phase is thermodynamically favored. For a single-component system, the chemical potential equals the molar Gibbs energy: $\mu = G_m$.

Phase diagrams map the regions of stability for each phase as a function of thermodynamic variables. The boundaries between single-phase regions are coexistence curves along which two phases are in equilibrium. Special points include the triple point, where three phases coexist ($F = 0$), and the critical point, where the distinction between liquid and gas vanishes.

Key Thermodynamic Relations

The fundamental relation for the Gibbs energy of a single-component system is:

$$dG = -SdT + VdP$$

From this we extract the Maxwell relations: $\left(\frac{\partial G}{\partial T}\right)_P = -S$ and $\left(\frac{\partial G}{\partial P}\right)_T = V$.

2. Derivation: The Clausius-Clapeyron Equation

The Clausius-Clapeyron equation gives the slope of any coexistence curve in the P-T phase diagram. It follows directly from the condition of phase equilibrium.

Step 1: Equilibrium Along the Coexistence Curve

Along the coexistence curve between phases 1 and 2, the chemical potentials are equal at every point:

$$\mu_1(T, P) = \mu_2(T, P)$$

If we move along the coexistence curve by $dT$ and $dP$, both chemical potentials must change by the same amount to maintain equilibrium:

$$d\mu_1 = d\mu_2$$

Step 2: Apply the Gibbs-Duhem Relation

For a single-component system, the molar Gibbs energy satisfies $dG_m = -S_m dT + V_m dP$, where $S_m$ and $V_m$ are the molar entropy and volume. Since $\mu = G_m$:

$$-S_{m,1}\,dT + V_{m,1}\,dP = -S_{m,2}\,dT + V_{m,2}\,dP$$

Step 3: Solve for dP/dT

Rearranging:

$$(V_{m,2} - V_{m,1})\,dP = (S_{m,2} - S_{m,1})\,dT$$

$$\boxed{\frac{dP}{dT} = \frac{\Delta S_m}{\Delta V_m} = \frac{\Delta H_m}{T\,\Delta V_m}}$$

In the last step we used $\Delta S = \Delta H / T$ at the transition temperature, since the phase transition is reversible at equilibrium. This is the Clapeyron equation, exact for any first-order phase transition.

Step 4: Clausius-Clapeyron Approximation for Liquid-Vapor Equilibrium

For the liquid-vapor transition, we make two approximations: (i) the molar volume of the liquid is negligible compared to the vapor ($V_{m,\ell} \ll V_{m,g}$), and (ii) the vapor behaves as an ideal gas ($V_{m,g} = RT/P$). Then:

$$\frac{dP}{dT} = \frac{\Delta H_{\text{vap}}}{T \cdot (RT/P)} = \frac{P\,\Delta H_{\text{vap}}}{RT^2}$$

$$\frac{1}{P}\frac{dP}{dT} = \frac{d(\ln P)}{dT} = \frac{\Delta H_{\text{vap}}}{RT^2}$$

Step 5: Integration (Constant Enthalpy of Vaporization)

If $\Delta H_{\text{vap}}$ is approximately constant over the temperature range of interest, we integrate between states $(T_1, P_1)$ and $(T_2, P_2)$:

$$\int_{P_1}^{P_2} \frac{dP}{P} = \frac{\Delta H_{\text{vap}}}{R}\int_{T_1}^{T_2}\frac{dT}{T^2}$$

$$\boxed{\ln\frac{P_2}{P_1} = -\frac{\Delta H_{\text{vap}}}{R}\left(\frac{1}{T_2} - \frac{1}{T_1}\right)}$$

This is the Clausius-Clapeyron equation in its integrated form. It predicts that a plot of $\ln P$ versus $1/T$ yields a straight line with slope $-\Delta H_{\text{vap}}/R$.

Physical Interpretation

The Clausius-Clapeyron equation explains why boiling points increase with pressure: since$\Delta H_{\text{vap}} > 0$ and $\Delta V_m > 0$ for liquid-to-vapor, we have $dP/dT > 0$. The anomalous case of water's ice-liquid boundary ($dP/dT < 0$) arises because ice has a larger molar volume than liquid water ($\Delta V_m < 0$), so the melting curve has a negative slope.

3. Derivation: The Gibbs Phase Rule

The Gibbs phase rule is one of the most powerful results in classical thermodynamics. It determines the number of degrees of freedom (independently variable intensive properties) for a system at equilibrium.

Step 1: Count the Intensive Variables

Consider a system with $C$ chemical components distributed among $P$ phases. The intensive state of each phase is specified by:

Temperature $T$ and pressure $P$ (the same for all phases at equilibrium)
The mole fractions of each component in each phase: $x_1^{(\alpha)}, x_2^{(\alpha)}, \ldots, x_C^{(\alpha)}$

Within each phase, the mole fractions satisfy the constraint $\sum_{i=1}^C x_i^{(\alpha)} = 1$, so each phase contributes $(C-1)$ independent composition variables. The total number of independent intensive variables is:

$$\text{Total variables} = 2 + P(C - 1)$$

The “2” accounts for temperature and pressure, shared across all phases.

Step 2: Count the Equilibrium Constraints

At equilibrium, the chemical potential of each component must be equal across all phases. For component $i$:

$$\mu_i^{(1)} = \mu_i^{(2)} = \cdots = \mu_i^{(P)}$$

This gives $(P - 1)$ independent equations for each of the $C$ components:

$$\text{Total constraints} = C(P - 1)$$

Step 3: The Phase Rule

The number of degrees of freedom is the difference between variables and constraints:

$$F = [2 + P(C - 1)] - [C(P - 1)]$$

$$F = 2 + PC - P - CP + C$$

$$\boxed{F = C - P + 2}$$

This is the Gibbs Phase Rule, first derived by J. Willard Gibbs in 1876.

Application: Water Phase Diagram

Region	C	P	F = C - P + 2	Interpretation
Single phase	1	1	2	T and P independently variable — 2D area on phase diagram
Two-phase coexistence	1	2	1	Only one variable free — coexistence curves (lines)
Triple point	1	3	0	No degrees of freedom — a unique point (T = 273.16 K, P = 611.73 Pa)

At the triple point, the system is completely determined: three phases of a single component leave zero degrees of freedom. This is why the triple point of water serves as a fundamental temperature reference in the SI system.

4. Derivation: Chemical Potential and Phase Stability

The fundamental criterion for phase stability is that the system adopts the phase with the lowest molar Gibbs energy (chemical potential) at the given temperature and pressure.

Step 1: Temperature Dependence of the Chemical Potential

From $dG = -SdT + VdP$, at constant pressure:

$$\left(\frac{\partial \mu}{\partial T}\right)_P = -S_m$$

Since $S_m > 0$ for all phases, the chemical potential always decreases with increasing temperature. The key insight is that the rate of decrease differs by phase:

$$S_m^{(\text{solid})} < S_m^{(\text{liquid})} < S_m^{(\text{gas})}$$

The gas has the steepest downward slope because it has the largest molar entropy (molecular disorder).

Step 2: Phase Stability Criterion

At any temperature and pressure, the thermodynamically stable phase is the one with the lowest chemical potential. Consider plotting $\mu(T)$ for each phase at fixed P:

At low T, $\mu^{(\text{solid})}$ is lowest — the solid is stable
At the melting point $T_m$: $\mu^{(\text{solid})} = \mu^{(\text{liquid})}$
Between $T_m$ and $T_b$: $\mu^{(\text{liquid})}$ is lowest
At the boiling point $T_b$: $\mu^{(\text{liquid})} = \mu^{(\text{gas})}$
Above $T_b$: $\mu^{(\text{gas})}$ is lowest — the gas is stable

Step 3: Pressure Dependence and Phase Boundaries

At constant temperature, the pressure dependence is:

$$\left(\frac{\partial \mu}{\partial P}\right)_T = V_m$$

Since $V_m^{(\text{gas})} \gg V_m^{(\text{liquid})} > V_m^{(\text{solid})}$, increasing pressure raises the chemical potential of the gas most rapidly. This explains why high pressure favors condensed phases: the gas chemical potential rises fastest, making it energetically unfavorable.

Step 4: Shift of Transition Temperatures with Pressure

Increasing pressure shifts the intersection point of the $\mu(T)$ curves. For the melting transition, the shift in melting point with pressure is given by the Clapeyron equation:

$$\Delta T_m = \frac{T_m \Delta V_m}{\Delta H_m} \Delta P$$

For most substances $\Delta V_m > 0$ on melting, so $T_m$ increases with pressure. Water is anomalous: $\Delta V_m < 0$, so its melting pointdecreases under pressure — ice can melt under the blade of an ice skate.

5. Derivation: Ehrenfest Classification of Phase Transitions

Paul Ehrenfest proposed a classification of phase transitions based on the behavior of the Gibbs energy and its derivatives at the transition point. This scheme distinguishes transitions by the order of the first derivative of G that shows a discontinuity.

First-Order Phase Transitions

In a first-order transition, the Gibbs energy $G$ itself is continuous at the transition point, but its first derivatives are discontinuous:

$$\Delta S = -\Delta\!\left(\frac{\partial G}{\partial T}\right)_P \neq 0 \quad \Rightarrow \quad \text{latent heat } \Delta H = T\Delta S$$

$$\Delta V = \Delta\!\left(\frac{\partial G}{\partial P}\right)_T \neq 0 \quad \Rightarrow \quad \text{volume change}$$

Examples: melting, boiling, sublimation. These transitions involve a latent heat and a discontinuous change in volume. The two phases can coexist at the transition temperature.

Second-Order Phase Transitions

In a second-order (continuous) transition, $G$ and its first derivatives ($S$ and $V$) are continuous, but the second derivatives show discontinuities:

$$\Delta C_P = -T\,\Delta\!\left(\frac{\partial^2 G}{\partial T^2}\right)_P \neq 0 \quad \text{(heat capacity jump)}$$

$$\Delta \alpha = \frac{1}{V}\Delta\!\left(\frac{\partial^2 G}{\partial T\,\partial P}\right) \neq 0 \quad \text{(thermal expansion jump)}$$

$$\Delta \kappa_T = -\frac{1}{V}\Delta\!\left(\frac{\partial^2 G}{\partial P^2}\right)_T \neq 0 \quad \text{(compressibility jump)}$$

Examples: superconducting transitions in metals (in zero field), superfluid transition in $^4\text{He}$, ferromagnetic-paramagnetic transition at the Curie point. There is no latent heat, no volume change, and no phase coexistence.

Ehrenfest Relations for Second-Order Transitions

For a second-order transition, the Clapeyron equation ($dP/dT = \Delta S / \Delta V$) gives the indeterminate form 0/0, since both $\Delta S = 0$ and $\Delta V = 0$. Ehrenfest resolved this by applying L'Hôpital's rule. The continuity of S along the transition line requires:

$$dS_1 = dS_2 \;\;\Rightarrow\;\; \left(\frac{\partial S_1}{\partial T}\right)_P dT + \left(\frac{\partial S_1}{\partial P}\right)_T dP = \left(\frac{\partial S_2}{\partial T}\right)_P dT + \left(\frac{\partial S_2}{\partial P}\right)_T dP$$

Using $(\partial S/\partial T)_P = C_P/T$ and $(\partial S/\partial P)_T = -V\alpha$(Maxwell relation), we obtain the first Ehrenfest relation:

$$\boxed{\frac{dP}{dT} = \frac{\Delta C_P}{TV\,\Delta\alpha}}$$

Second Ehrenfest Relation

Similarly, the continuity of V along the transition line requires $dV_1 = dV_2$. Using$(\partial V/\partial T)_P = V\alpha$ and $(\partial V/\partial P)_T = -V\kappa_T$:

$$V\alpha_1\,dT - V\kappa_{T,1}\,dP = V\alpha_2\,dT - V\kappa_{T,2}\,dP$$

$$\boxed{\frac{dP}{dT} = \frac{\Delta\alpha}{\Delta\kappa_T}}$$

The two Ehrenfest relations can be combined to give the Prigogine-Defay ratio:$\Pi = \Delta C_P \Delta \kappa_T / (TV(\Delta\alpha)^2)$. For a true second-order transition, $\Pi = 1$; deviations indicate more complex behavior (e.g., the glass transition typically has $\Pi > 1$).

6. Applications of Phase Equilibria

Water: The Anomalous Phase Diagram

Water exhibits several unusual features in its phase diagram:

Negative solid-liquid slope: The ice I$_h$-water boundary has $dP/dT < 0$ because ice is less dense than liquid water ($\Delta V_m < 0$ on melting). This is due to the open tetrahedral hydrogen-bonding network in ice.
Multiple solid phases: At high pressures, water forms at least 17 distinct crystalline ice polymorphs (Ice II through Ice XIX), each with different hydrogen-bonding arrangements.
Critical point: $T_c = 647.1$ K, $P_c = 220.6$ bar. Above this, water exists as a supercritical fluid with unique solvent properties.
Triple point: $T = 273.16$ K, $P = 611.73$ Pa. This serves as the defining fixed point for the kelvin in the SI system.

Carbon Dioxide: Supercritical Extraction

CO$_2$ has a critical point at $T_c = 304.1$ K and $P_c = 73.8$ bar, making it accessible with moderate equipment. Supercritical CO$_2$ (scCO$_2$) combines:

Gas-like diffusivity and low viscosity (rapid mass transfer)
Liquid-like density and solvating power (effective extraction)
Tunable solvent strength via pressure adjustment

Applications include decaffeination of coffee, extraction of essential oils and flavors, pharmaceutical processing, and green chemistry solvent replacement.

Helium: Quantum Phase Diagram

Helium-4 ($^4\text{He}$) has a unique phase diagram reflecting its quantum nature:

No triple point at ambient pressure — helium remains liquid down to absolute zero due to its large zero-point energy
Superfluid transition ($\lambda$-point) at $T_\lambda = 2.172$ K: a second-order phase transition where viscosity drops to zero due to Bose-Einstein condensation
Solidification requires pressures above ~25 bar, even at $T = 0$ K

Triple Points as Temperature Standards

Because triple points have $F = 0$ degrees of freedom, they occur at unique, reproducible temperatures and pressures. Several triple points serve as primary fixed points on the International Temperature Scale (ITS-90):

Substance	Triple Point (K)	Use
Hydrogen	13.8033	Cryogenic calibration
Neon	24.5561	Low-temperature thermometry
Oxygen	54.3584	Low-temperature thermometry
Argon	83.8058	Cryogenic calibration
Mercury	234.3156	Intermediate temperatures
Water	273.16	Primary SI reference
Gallium	302.9146	Near-ambient calibration

7. Historical Context

The Pioneers of Phase Equilibria

Benoît Paul Émile Clapeyron (1834) published the first quantitative relation between vapor pressure and temperature, building on the earlier qualitative insights of Sadi Carnot. Clapeyron's work introduced the P-V indicator diagram and established the thermodynamic foundations of steam engine theory.

Rudolf Clausius (1850) refined Clapeyron's relation using his new concept of entropy. The Clausius-Clapeyron equation in its modern form combines both contributions: Clapeyron's geometric insight about coexistence curves and Clausius's thermodynamic framework connecting heat and entropy.

Josiah Willard Gibbs (1876) revolutionized the field with his monumental paper “On the Equilibrium of Heterogeneous Substances.” Gibbs introduced the chemical potential, derived the phase rule ($F = C - P + 2$), and established the systematic framework for understanding multi-component, multi-phase systems. His work was largely ignored in the English-speaking world until translated and promoted by Wilhelm Ostwald and James Clerk Maxwell.

Paul Ehrenfest (1933) proposed his classification of phase transitions at the Solvay Conference. While the strict Ehrenfest classification has been superseded by modern renormalization group theory (which shows that most “second-order” transitions involve divergences rather than simple discontinuities), the scheme remains a valuable pedagogical tool and correctly distinguishes the fundamental categories of phase transition behavior.

Lev Landau (1937) developed a general theory of second-order phase transitions based on symmetry breaking and order parameters. Landau theory provides a mean-field description valid far from the critical point and introduced concepts (order parameter, broken symmetry, critical exponents) that remain central to modern condensed matter physics.

8. Python Simulation: Phase Diagram of Water

The following simulation uses the Clausius-Clapeyron equation with real thermodynamic data to construct the phase diagram of water, showing all three coexistence curves and the triple point.

Python

phase_diagram_water.py168 lines

#!/usr/bin/env python3
"""
phase_diagram_water.py
Compute and plot the phase diagram of water using the Clausius-Clapeyron
equation with real thermodynamic data. Shows solid-liquid, liquid-vapor,
and solid-vapor coexistence curves, plus the triple point and critical point.
Uses numpy only -- no scipy.
"""
import numpy as np

# ============================================================
# Physical constants and reference data
# ============================================================
R = 8.314462618        # J/(mol*K)

# Triple point of water
T_trip = 273.16        # K
P_trip = 611.73        # Pa

# Critical point of water
T_crit = 647.096       # K
P_crit = 22.064e6      # Pa (220.64 bar)

# Thermodynamic data at/near the triple point
DH_vap = 45054.0       # J/mol  (enthalpy of vaporization at 373 K, approx)
DH_sub = 51059.0       # J/mol  (enthalpy of sublimation at 273 K)
DH_fus = 6010.0        # J/mol  (enthalpy of fusion at 273.15 K)

# Molar volumes at the triple point
Vm_ice   = 19.65e-6    # m^3/mol (ice Ih)
Vm_water = 18.02e-6    # m^3/mol (liquid water)
DV_fus = Vm_water - Vm_ice  # negative! (water anomaly)

print("=" * 65)
print("   PHASE DIAGRAM OF WATER (Clausius-Clapeyron)")
print("=" * 65)
print(f"\nTriple point:   T = {T_trip} K, P = {P_trip:.2f} Pa")
print(f"Critical point: T = {T_crit:.1f} K, P = {P_crit/1e6:.3f} MPa")
print(f"\nDelta_H_vap = {DH_vap:.0f} J/mol")
print(f"Delta_H_sub = {DH_sub:.0f} J/mol")
print(f"Delta_H_fus = {DH_fus:.0f} J/mol")
print(f"Delta_V_fus = {DV_fus*1e6:.2f} cm^3/mol (negative -- anomalous!)")

# ============================================================
# 1. LIQUID-VAPOR coexistence (Clausius-Clapeyron)
# ============================================================
# ln(P/P_trip) = -(DH_vap/R)(1/T - 1/T_trip)
T_lv = np.linspace(T_trip, T_crit, 500)
P_lv = P_trip * np.exp(-(DH_vap / R) * (1.0 / T_lv - 1.0 / T_trip))

# Cap at critical point pressure
P_lv = np.clip(P_lv, 0, P_crit)

print("\n--- Liquid-Vapor Coexistence ---")
for T_check in [300, 350, 373.15, 400, 450, 500]:
    P_check = P_trip * np.exp(-(DH_vap / R) * (1.0 / T_check - 1.0 / T_trip))
    print(f"  T = {T_check:7.2f} K  =>  P = {P_check/1000:.2f} kPa ({P_check/101325:.4f} atm)")

# ============================================================
# 2. SOLID-VAPOR coexistence (sublimation curve)
# ============================================================
# ln(P/P_trip) = -(DH_sub/R)(1/T - 1/T_trip)
T_sv = np.linspace(200, T_trip, 300)
P_sv = P_trip * np.exp(-(DH_sub / R) * (1.0 / T_sv - 1.0 / T_trip))

print("\n--- Solid-Vapor Coexistence (Sublimation) ---")
for T_check in [200, 220, 240, 260, 273.16]:
    P_check = P_trip * np.exp(-(DH_sub / R) * (1.0 / T_check - 1.0 / T_trip))
    print(f"  T = {T_check:7.2f} K  =>  P = {P_check:.4f} Pa")

# ============================================================
# 3. SOLID-LIQUID coexistence (Clapeyron equation, exact)
# ============================================================
# dP/dT = DH_fus / (T * DV_fus)
# Integrate: P - P_trip = (DH_fus/DV_fus) * ln(T/T_trip)
T_sl = np.linspace(T_trip, 251, 300)  # goes to lower T at higher P (negative slope!)
P_sl = P_trip + (DH_fus / DV_fus) * np.log(T_sl / T_trip)

print("\n--- Solid-Liquid Coexistence (Melting) ---")
print(f"  Slope dP/dT = DH_fus/(T_trip*DV_fus) = {DH_fus/(T_trip*DV_fus)/1e6:.2f} MPa/K")
print(f"  NOTE: Negative slope (ice less dense than water)")
for P_check_MPa in [0.1, 100, 200, 500, 1000]:
    P_val = P_check_MPa * 1e6
    T_val = T_trip * np.exp((P_val - P_trip) * DV_fus / DH_fus)
    print(f"  P = {P_check_MPa:6.0f} MPa  =>  T = {T_val:.2f} K ({T_val - 273.15:.2f} C)")

# ============================================================
# 4. Verify the Clausius-Clapeyron integrated form
# ============================================================
print("\n--- Verification: Normal Boiling Point ---")
T1 = T_trip
P1 = P_trip
P2_atm = 101325.0  # 1 atm
# ln(P2/P1) = -(DH_vap/R)(1/T2 - 1/T1)
# 1/T2 = 1/T1 - (R/DH_vap)*ln(P2/P1)
inv_T2 = 1.0/T1 - (R/DH_vap)*np.log(P2_atm/P1)
T2_calc = 1.0 / inv_T2
print(f"  Calculated boiling point at 1 atm: T = {T2_calc:.2f} K ({T2_calc - 273.15:.2f} C)")
print(f"  Experimental value: 373.15 K (100.00 C)")
print(f"  Deviation: {abs(T2_calc - 373.15):.2f} K")

# ============================================================
# PLOTTING
# ============================================================
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

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

# --- Left panel: Full phase diagram (log scale) ---
ax1 = axes[0]
ax1.semilogy(T_lv, P_lv / 1e6, 'r-', lw=2.5, label='Liquid-Vapor')
ax1.semilogy(T_sv, P_sv / 1e6, 'b-', lw=2.5, label='Solid-Vapor')
ax1.semilogy(T_sl, P_sl / 1e6, 'g-', lw=2.5, label='Solid-Liquid')

# Triple point
ax1.plot(T_trip, P_trip/1e6, 'ko', ms=10, zorder=5)
ax1.annotate('Triple Point\n(273.16 K, 611.73 Pa)',
             xy=(T_trip, P_trip/1e6), xytext=(310, 1e-5),
             fontsize=9, color='black',
             arrowprops=dict(arrowstyle='->', color='black'),
             bbox=dict(boxstyle='round,pad=0.3', facecolor='yellow', alpha=0.8))

# Critical point
ax1.plot(T_crit, P_crit/1e6, 'r*', ms=15, zorder=5)
ax1.annotate('Critical Point\n(647.1 K, 22.06 MPa)',
             xy=(T_crit, P_crit/1e6), xytext=(500, 50),
             fontsize=9, color='red',
             arrowprops=dict(arrowstyle='->', color='red'),
             bbox=dict(boxstyle='round,pad=0.3', facecolor='lightyellow', alpha=0.8))

# Phase labels
ax1.text(220, 0.001, 'SOLID', fontsize=14, fontweight='bold', color='blue', alpha=0.6)
ax1.text(350, 0.0001, 'VAPOR', fontsize=14, fontweight='bold', color='red', alpha=0.6)
ax1.text(300, 100, 'LIQUID', fontsize=14, fontweight='bold', color='green', alpha=0.6)

ax1.set_xlabel('Temperature (K)', fontsize=13)
ax1.set_ylabel('Pressure (MPa)', fontsize=13)
ax1.set_title('Phase Diagram of Water', fontsize=14, fontweight='bold')
ax1.legend(loc='lower right', fontsize=10)
ax1.set_xlim(190, 680)
ax1.set_ylim(1e-7, 1e3)
ax1.grid(True, alpha=0.3)

# --- Right panel: Clausius-Clapeyron plot (ln P vs 1/T) ---
ax2 = axes[1]
inv_T_lv = 1000.0 / T_lv  # 1000/T for readable axis
inv_T_sv = 1000.0 / T_sv
ax2.plot(inv_T_lv, np.log(P_lv), 'r-', lw=2.5, label='Liquid-Vapor')
ax2.plot(inv_T_sv, np.log(P_sv), 'b-', lw=2.5, label='Solid-Vapor')

# Show the slopes
slope_vap = -DH_vap / (R * 1000)  # per (1000/T)
slope_sub = -DH_sub / (R * 1000)
ax2.text(2.0, 14, f'Slope (L-V) = {slope_vap:.2f}', fontsize=10, color='red')
ax2.text(3.8, 4, f'Slope (S-V) = {slope_sub:.2f}', fontsize=10, color='blue')

ax2.set_xlabel('1000/T (K$^{-1}$)', fontsize=13)
ax2.set_ylabel('ln(P / Pa)', fontsize=13)
ax2.set_title('Clausius-Clapeyron Plot', fontsize=14, fontweight='bold')
ax2.legend(loc='upper right', fontsize=10)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('water_phase_diagram.png', dpi=150, bbox_inches='tight')
print("\nPlot saved to water_phase_diagram.png")
print("=" * 65)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

9. Fortran Simulation: Vapor Pressures from the Antoine Equation

The Antoine equation is an empirical relation widely used for correlating vapor pressure data. It takes the form $\log_{10}(P) = A - B/(C + T)$, where A, B, C are substance-specific parameters. This Fortran program computes vapor pressures and normal boiling points for several common substances using tabulated Antoine parameters.

Antoine Equation Vapor Pressures and Boiling Points

Fortran

antoine_vapor_pressure.f90156 lines

program antoine_vapor_pressure
  implicit none

! Antoine equation: log10(P_mmHg) = A - B/(C + T_celsius)
  ! Parameters for common substances (valid ranges given)

integer, parameter :: n_sub = 6
  character(len=20) :: names(n_sub)
  double precision :: A(n_sub), B(n_sub), C(n_sub)
  double precision :: T_min(n_sub), T_max(n_sub)  ! validity range (C)
  double precision :: T, P_mmHg, P_Pa, P_atm
  double precision :: T_bp, T_lo, T_hi, T_mid, P_mid, f_mid
  integer :: i, j, iter
  double precision, parameter :: mmHg_to_Pa = 133.322d0

! --- Substance data (Antoine parameters, P in mmHg, T in C) ---
  names(1) = 'Water'
  A(1) = 8.07131d0;  B(1) = 1730.63d0;  C(1) = 233.426d0
  T_min(1) = 1.0d0;  T_max(1) = 100.0d0

names(2) = 'Ethanol'
  A(2) = 8.20417d0;  B(2) = 1642.89d0;  C(2) = 230.300d0
  T_min(2) = -57.0d0; T_max(2) = 80.0d0

names(3) = 'Acetone'
  A(3) = 7.02447d0;  B(3) = 1161.0d0;   C(3) = 224.0d0
  T_min(3) = -20.0d0; T_max(3) = 77.0d0

names(4) = 'Benzene'
  A(4) = 6.90565d0;  B(4) = 1211.033d0;  C(4) = 220.790d0
  T_min(4) = 8.0d0;   T_max(4) = 80.0d0

names(5) = 'Hexane'
  A(5) = 6.87776d0;  B(5) = 1171.53d0;   C(5) = 224.366d0
  T_min(5) = -25.0d0;  T_max(5) = 92.0d0

names(6) = 'Diethyl Ether'
  A(6) = 6.92032d0;  B(6) = 1064.07d0;   C(6) = 228.800d0
  T_min(6) = -60.0d0;  T_max(6) = 35.0d0

write(*,'(A)') '=================================================================='
  write(*,'(A)') '  ANTOINE EQUATION: VAPOR PRESSURES AND BOILING POINTS'
  write(*,'(A)') '=================================================================='
  write(*,'(A)') ''
  write(*,'(A)') 'Antoine equation: log10(P/mmHg) = A - B/(C + T_celsius)'
  write(*,'(A)') ''

! --- Print Antoine parameters ---
  write(*,'(A)') '--- Antoine Parameters ---'
  write(*,'(A6,A20,A12,A12,A12)') ' #', 'Substance', 'A', 'B', 'C'
  write(*,'(A)') '------------------------------------------------------------------'
  do i = 1, n_sub
    write(*,'(I3,3X,A20,3F12.4)') i, names(i), A(i), B(i), C(i)
  end do
  write(*,'(A)') ''

! --- Compute vapor pressures at selected temperatures ---
  write(*,'(A)') '--- Vapor Pressures at Selected Temperatures ---'
  write(*,'(A)') ''
  do i = 1, n_sub
    write(*,'(A,A)') 'Substance: ', trim(names(i))
    write(*,'(A8,A15,A15,A15)') 'T (C)', 'P (mmHg)', 'P (kPa)', 'P (atm)'
    write(*,'(A)') '------------------------------------------------------'
    do j = 0, 4
      T = T_min(i) + j * (T_max(i) - T_min(i)) / 4.0d0
      P_mmHg = 10.0d0**(A(i) - B(i)/(C(i) + T))
      P_Pa = P_mmHg * mmHg_to_Pa
      P_atm = P_Pa / 101325.0d0
      write(*,'(F8.1,F15.3,F15.4,F15.6)') T, P_mmHg, P_Pa/1000.0d0, P_atm
    end do
    write(*,'(A)') ''
  end do

! --- Compute normal boiling points (P = 760 mmHg) using bisection ---
  write(*,'(A)') '--- Normal Boiling Points (P = 760 mmHg = 1 atm) ---'
  write(*,'(A)') ''
  write(*,'(A20,A15,A15)') 'Substance', 'T_bp (C)', 'T_bp (K)'
  write(*,'(A)') '--------------------------------------------------'

do i = 1, n_sub
    ! Bisection to find T where P = 760 mmHg
    T_lo = -200.0d0
    T_hi = 500.0d0

do iter = 1, 100
      T_mid = 0.5d0 * (T_lo + T_hi)
      P_mid = 10.0d0**(A(i) - B(i)/(C(i) + T_mid))
      f_mid = P_mid - 760.0d0

if (abs(f_mid) < 1.0d-8) exit

if (f_mid > 0.0d0) then
        T_hi = T_mid
      else
        T_lo = T_mid
      end if
    end do

T_bp = T_mid
    write(*,'(A20,F15.3,F15.3)') trim(names(i)), T_bp, T_bp + 273.15d0
  end do

write(*,'(A)') ''
  write(*,'(A)') '--- Comparison with Experimental Boiling Points ---'
  write(*,'(A)') ''
  write(*,'(A)') 'Water:          calc vs exp 100.0 C'
  write(*,'(A)') 'Ethanol:        calc vs exp  78.4 C'
  write(*,'(A)') 'Acetone:        calc vs exp  56.1 C'
  write(*,'(A)') 'Benzene:        calc vs exp  80.1 C'
  write(*,'(A)') 'Hexane:         calc vs exp  68.7 C'
  write(*,'(A)') 'Diethyl Ether:  calc vs exp  34.6 C'

write(*,'(A)') ''

! --- Enthalpy of vaporization from Clausius-Clapeyron ---
  write(*,'(A)') '--- Estimated Delta_H_vap from Antoine Slope ---'
  write(*,'(A)') ''
  write(*,'(A)') 'Using: d(lnP)/d(1/T) = -Delta_H_vap/R at T_bp'
  write(*,'(A)') ''
  write(*,'(A20,A18)') 'Substance', 'Delta_H_vap (kJ/mol)'
  write(*,'(A)') '----------------------------------------------'

do i = 1, n_sub
    ! Find boiling point first
    T_lo = -200.0d0
    T_hi = 500.0d0
    do iter = 1, 100
      T_mid = 0.5d0 * (T_lo + T_hi)
      P_mid = 10.0d0**(A(i) - B(i)/(C(i) + T_mid))
      f_mid = P_mid - 760.0d0
      if (abs(f_mid) < 1.0d-8) exit
      if (f_mid > 0.0d0) then
        T_hi = T_mid
      else
        T_lo = T_mid
      end if
    end do
    T_bp = T_mid + 273.15d0  ! Convert to K

! Antoine: ln(P) = ln(10) * [A - B/(C + T)]
    ! d(lnP)/dT = ln(10) * B / (C + T)^2
    ! d(lnP)/d(1/T) = -T^2 * d(lnP)/dT = -ln(10)*B*T^2/(C+T-273.15)^2
    ! Wait: T in Antoine is Celsius. Let T_C = T_K - 273.15
    ! d(lnP)/d(1/T_K) = -T_K^2 * dln(P)/dT_K = -T_K^2 * ln(10)*B/(C+T_C)^2
    ! Delta_H_vap = R * ln(10) * B * T_K^2 / (C + T_C)^2

T = T_bp - 273.15d0  ! Celsius
    P_Pa = 8.314462618d0 * log(10.0d0) * B(i) * T_bp**2 &
           / (C(i) + T)**2
    write(*,'(A20,F18.2)') trim(names(i)), P_Pa / 1000.0d0
  end do

write(*,'(A)') ''
  write(*,'(A)') '=================================================================='

end program antoine_vapor_pressure

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

10. Summary and Key Results

Central Equations of Phase Equilibria

Clapeyron Equation (exact)

$$\frac{dP}{dT} = \frac{\Delta S_m}{\Delta V_m} = \frac{\Delta H_m}{T\,\Delta V_m}$$

Clausius-Clapeyron Equation (integrated, liquid-vapor)

$$\ln\frac{P_2}{P_1} = -\frac{\Delta H_{\text{vap}}}{R}\left(\frac{1}{T_2} - \frac{1}{T_1}\right)$$

Gibbs Phase Rule

$$F = C - P + 2$$

Ehrenfest Relations (second-order transitions)

$$\frac{dP}{dT} = \frac{\Delta C_P}{TV\,\Delta\alpha} = \frac{\Delta\alpha}{\Delta\kappa_T}$$