← Part II: Thermodynamics & Equations of State

Chemical Potential

Introduction: The Driving Force for Chemical Change

The chemical potential is arguably the most important intensive quantity in chemical thermodynamics. It governs every process involving the transfer of matter: chemical reactions, phase transitions, diffusion, osmosis, and the formation of mixtures. Just as temperature differences drive heat flow and pressure differences drive mechanical work, differences in chemical potential drive the transfer of matter from one phase, region, or chemical identity to another.

The chemical potential $\mu_i$ of species $i$ in a mixture answers the fundamental question: how much does the Gibbs energy of the system change when an infinitesimal amount of species $i$ is added, while temperature, pressure, and the amounts of all other species are held fixed?

At equilibrium, the chemical potential of each species must be the same in every phase and every region of the system. Any imbalance in chemical potential is a thermodynamic driving force that propels the system toward equilibrium. This principle unifies an enormous range of phenomena under a single mathematical framework.

Key Concept: Equilibrium Condition

For a system with multiple phases $\alpha, \beta, \gamma, \ldots$, thermodynamic equilibrium requires:

$$\mu_i^\alpha = \mu_i^\beta = \mu_i^\gamma = \cdots \quad \text{for every species } i$$

This condition, combined with thermal equilibrium ($T^\alpha = T^\beta$) and mechanical equilibrium ($P^\alpha = P^\beta$), completely specifies the state of a multi-phase, multi-component system at equilibrium.

Derivation 1: Chemical Potential from Thermodynamic Potentials

Definition via the Gibbs Energy

The chemical potential of species $i$ is defined as the partial molar Gibbs energy:

$$\mu_i \equiv \left(\frac{\partial G}{\partial n_i}\right)_{T,P,n_{j \neq i}}$$

This is the most natural definition because experiments are typically conducted at constant temperature and pressure. However, the chemical potential can equivalently be expressed as a partial derivative of any thermodynamic potential with respect to $n_i$, provided the appropriate natural variables are held constant.

Equivalence with Other Potentials

Starting from the fundamental equations for each thermodynamic potential, we can show that$\mu_i$ appears in all of them:

Internal Energy (natural variables S, V, n):

$$dU = TdS - PdV + \sum_i \mu_i \, dn_i$$

$\Rightarrow \mu_i = \left(\frac{\partial U}{\partial n_i}\right)_{S,V,n_{j \neq i}}$

Enthalpy (natural variables S, P, n):

$$dH = TdS + VdP + \sum_i \mu_i \, dn_i$$

$\Rightarrow \mu_i = \left(\frac{\partial H}{\partial n_i}\right)_{S,P,n_{j \neq i}}$

Helmholtz Free Energy (natural variables T, V, n):

$$dA = -SdT - PdV + \sum_i \mu_i \, dn_i$$

$\Rightarrow \mu_i = \left(\frac{\partial A}{\partial n_i}\right)_{T,V,n_{j \neq i}}$

Gibbs Free Energy (natural variables T, P, n):

$$dG = -SdT + VdP + \sum_i \mu_i \, dn_i$$

$\Rightarrow \mu_i = \left(\frac{\partial G}{\partial n_i}\right)_{T,P,n_{j \neq i}}$

Derivation of the Fundamental Equation for G

We start from the combined first and second laws for an open system. For a reversible process:

$$dU = TdS - PdV + \sum_i \mu_i \, dn_i$$

Since $G = U + PV - TS = H - TS$, we compute the total differential:

$$dG = dU + PdV + VdP - TdS - SdT$$

Substituting the expression for $dU$:

$$dG = (TdS - PdV + \sum_i \mu_i \, dn_i) + PdV + VdP - TdS - SdT$$

The $TdS$ and $PdV$ terms cancel, yielding:

$$\boxed{dG = VdP - SdT + \sum_i \mu_i \, dn_i}$$

This is the fundamental equation for the Gibbs energy of an open, multi-component system. At constant$T$ and $P$, it reduces to $dG = \sum_i \mu_i \, dn_i$, confirming that the chemical potential controls the change in Gibbs energy due to changes in composition.

Derivation 2: The Gibbs-Duhem Equation

Euler's Theorem and the Integrated Form of G

The Gibbs energy is an extensive quantity: if we scale all mole numbers by a factor $\lambda$while holding $T$ and $P$ fixed, then $G$ scales by the same factor. Mathematically, $G(T, P, \lambda n_1, \lambda n_2, \ldots) = \lambda \, G(T, P, n_1, n_2, \ldots)$.

By Euler's theorem for homogeneous functions of degree one, this implies:

$$G = \sum_i n_i \left(\frac{\partial G}{\partial n_i}\right)_{T,P,n_{j \neq i}} = \sum_i n_i \mu_i$$

This is the Euler equation for the Gibbs energy. It states that the total Gibbs energy is simply the sum of each species' mole number weighted by its chemical potential. This result is remarkable: it tells us that $G$ can be reconstructed entirely from the intensive quantities $\mu_i$.

Derivation of the Gibbs-Duhem Relation

Take the total differential of the Euler equation $G = \sum_i n_i \mu_i$:

$$dG = \sum_i n_i \, d\mu_i + \sum_i \mu_i \, dn_i$$

But we already know from the fundamental equation that:

$$dG = VdP - SdT + \sum_i \mu_i \, dn_i$$

Subtracting the fundamental equation from the total differential of the Euler equation:

$$\sum_i n_i \, d\mu_i + \sum_i \mu_i \, dn_i - VdP + SdT - \sum_i \mu_i \, dn_i = 0$$

The $\sum_i \mu_i \, dn_i$ terms cancel, giving the Gibbs-Duhem equation:

$$\boxed{\sum_i n_i \, d\mu_i = -SdT + VdP}$$

Gibbs-Duhem at Constant T and P

At constant temperature and pressure, the right-hand side vanishes:

$$\sum_i n_i \, d\mu_i = 0 \quad \text{(constant } T, P\text{)}$$

Dividing by the total number of moles $n = \sum_i n_i$ and introducing mole fractions $x_i = n_i / n$:

$$\boxed{\sum_i x_i \, d\mu_i = 0 \quad \text{(constant } T, P\text{)}}$$

This is a powerful constraint: in a mixture at constant $T$ and $P$, the chemical potentials of the components are not independent. If we know how $\mu$varies with composition for all but one component, the Gibbs-Duhem equation determines the remaining one. For a binary system, $x_1 \, d\mu_1 + x_2 \, d\mu_2 = 0$, so$d\mu_2 = -(x_1/x_2) \, d\mu_1$.

Derivation 3: Chemical Potential of Ideal Gas Mixtures

Pure Ideal Gas

For a pure ideal gas at constant temperature, the molar Gibbs energy satisfies:

$$dG_m = V_m \, dP = \frac{RT}{P} \, dP$$

Integrating from the standard pressure $P^\circ$ (typically 1 bar) to pressure $P$:

$$G_m(T,P) = G_m^\circ(T) + RT \ln\!\left(\frac{P}{P^\circ}\right)$$

Since for a pure substance the molar Gibbs energy equals the chemical potential ($G_m = \mu$):

$$\mu(T,P) = \mu^\circ(T) + RT \ln\!\left(\frac{P}{P^\circ}\right)$$

Ideal Gas Mixture: Partial Pressure Form

In an ideal gas mixture, each component behaves as if it alone occupied the entire volume at its partial pressure $P_i = x_i P$. The chemical potential of component $i$in the mixture is therefore:

$$\boxed{\mu_i(T, P, \{x\}) = \mu_i^\circ(T) + RT \ln\!\left(\frac{P_i}{P^\circ}\right)}$$

where $\mu_i^\circ(T)$ is the standard chemical potential of pure $i$at $P^\circ$ and temperature $T$.

Mole Fraction Form

Substituting $P_i = x_i P$ into the expression above:

$$\mu_i = \mu_i^\circ(T) + RT \ln\!\left(\frac{x_i P}{P^\circ}\right)$$

Using the properties of logarithms:

$$\boxed{\mu_i = \mu_i^\circ(T) + RT \ln x_i + RT \ln\!\left(\frac{P}{P^\circ}\right)}$$

The term $RT \ln x_i$ is always negative (since $0 < x_i < 1$), reflecting the fact that mixing always lowers the chemical potential in an ideal gas mixture. This is the thermodynamic origin of the entropy of mixing:$\Delta_{\text{mix}} G = nRT \sum_i x_i \ln x_i < 0$, confirming that mixing of ideal gases is always spontaneous.

Derivation 4: Equilibrium Constant from Chemical Potential

Equilibrium Condition for a Chemical Reaction

Consider a general gas-phase reaction $\sum_i \nu_i A_i = 0$ where $\nu_i$ are stoichiometric coefficients (positive for products, negative for reactants). At equilibrium at constant$T$ and $P$, the Gibbs energy is minimized, which requires:

$$\sum_i \nu_i \mu_i = 0 \quad \text{(at equilibrium)}$$

Substituting the ideal gas expression for each $\mu_i$:

$$\sum_i \nu_i \left[\mu_i^\circ(T) + RT \ln\!\left(\frac{P_i}{P^\circ}\right)\right] = 0$$

Separating the standard and non-standard parts:

$$\sum_i \nu_i \mu_i^\circ(T) + RT \sum_i \nu_i \ln\!\left(\frac{P_i}{P^\circ}\right) = 0$$

Deriving $\Delta G^\circ = -RT \ln K$

Identifying $\Delta G^\circ(T) = \sum_i \nu_i \mu_i^\circ(T)$ as the standard Gibbs energy of reaction, and using $\sum_i \nu_i \ln(P_i/P^\circ) = \ln \prod_i (P_i/P^\circ)^{\nu_i}$:

$$\Delta G^\circ + RT \ln \prod_i \left(\frac{P_i}{P^\circ}\right)^{\nu_i} = 0$$

We define the thermodynamic equilibrium constant as the product at equilibrium:

$$K(T) = \prod_i \left(\frac{P_i}{P^\circ}\right)^{\nu_i}\bigg|_{\text{eq}}$$

Yielding the central result:

$$\boxed{\Delta G^\circ = -RT \ln K}$$

This equation connects the measurable standard Gibbs energy of reaction (tabulated from calorimetric data) to the equilibrium composition. A large negative $\Delta G^\circ$ implies$K \gg 1$ (products favored), while a large positive $\Delta G^\circ$implies $K \ll 1$ (reactants favored).

The van't Hoff Equation

To find how $K$ depends on temperature, start from$\ln K = -\Delta G^\circ / (RT)$ and differentiate with respect to $T$. Using the Gibbs-Helmholtz equation:

$$\frac{\partial}{\partial T}\left(\frac{\Delta G^\circ}{T}\right)_P = -\frac{\Delta H^\circ}{T^2}$$

Since $\ln K = -\Delta G^\circ / (RT)$:

$$\boxed{\frac{d(\ln K)}{dT} = \frac{\Delta H^\circ}{RT^2}}$$

This is the van't Hoff equation. Its consequences are profound:

For an exothermic reaction ($\Delta H^\circ < 0$), increasing $T$ decreases $K$ — the equilibrium shifts toward reactants (Le Chatelier's principle).
For an endothermic reaction ($\Delta H^\circ > 0$), increasing $T$ increases $K$ — the equilibrium shifts toward products.
Integrating (assuming $\Delta H^\circ$ is approximately constant): $\ln\!\left(\frac{K(T_2)}{K(T_1)}\right) = -\frac{\Delta H^\circ}{R}\left(\frac{1}{T_2} - \frac{1}{T_1}\right)$

Applications of Chemical Potential

Chemical Equilibrium Calculations

The condition $\sum_i \nu_i \mu_i = 0$ combined with $\Delta G^\circ = -RT \ln K$allows quantitative prediction of equilibrium compositions from tabulated thermodynamic data.

Haber process: optimal T and P for NH$_3$ yield
Combustion equilibria at high temperatures
Acid-base equilibria in solution
Solubility product calculations

Reaction Spontaneity

The reaction Gibbs energy $\Delta_r G = \sum_i \nu_i \mu_i$ determines spontaneity under actual (non-standard) conditions:

$\Delta_r G = \Delta G^\circ + RT \ln Q$, where $Q$ is the reaction quotient. The reaction proceeds forward when $\Delta_r G < 0$(i.e., when $Q < K$).

Fugacity and Fugacity Coefficients

For real gases, the ideal gas expression is generalized by replacing pressure with fugacity $f_i$:

$\mu_i = \mu_i^\circ + RT \ln(f_i / P^\circ)$

The fugacity coefficient $\phi_i = f_i / (x_i P)$ measures the deviation from ideal behavior. For an ideal gas, $\phi_i = 1$. Fugacity coefficients can be computed from equations of state via $\ln \phi_i = \int_0^P (Z_i - 1) \, dP/P$.

Activity Coefficients

For non-ideal solutions (liquid or solid), the chemical potential is written in terms of activity$a_i$:

$\mu_i = \mu_i^\circ + RT \ln a_i = \mu_i^\circ + RT \ln(\gamma_i x_i)$

The activity coefficient $\gamma_i$ captures deviations from ideal solution behavior. Models such as Margules, van Laar, Wilson, NRTL, and UNIQUAC express $\gamma_i$ as functions of composition. The Gibbs-Duhem equation provides crucial consistency checks for these models.

Phase Equilibria and the Chemical Potential

Phase equilibrium requires $\mu_i^\alpha = \mu_i^\beta$ for each component across every pair of coexisting phases. This single condition, when combined with appropriate expressions for $\mu_i$ in each phase, produces all the classical results:

Clausius-Clapeyron equation for the slope of a one-component phase boundary: $dP/dT = \Delta_{trs} H / (T \Delta_{trs} V)$
Raoult's law for ideal solutions: $P_i = x_i P_i^*$ follows when $\gamma_i = 1$
Henry's law for dilute solutions: $P_i = K_{H,i} x_i$ with a composition-independent Henry's constant
Colligative properties (boiling point elevation, freezing point depression, osmotic pressure) all derive from the lowering of solvent chemical potential by solute

Historical Context

J. Willard Gibbs (1876)

Gibbs introduced the chemical potential in his monumental paper "On the Equilibrium of Heterogeneous Substances" (1876–1878). He recognized that the condition of equal chemical potential across phases provides the complete criterion for multi-component equilibrium. Gibbs' work was so far ahead of its time that it was largely unappreciated by chemists for decades, until translated into more accessible form by others. His formulation of the phase rule $F = C - P + 2$ is a direct consequence of the chemical potential framework.

Gilbert N. Lewis (1901–1923)

Lewis made Gibbs' abstract formalism practical by introducing the concepts of fugacity (1901) and activity (1907). Fugacity provides a "corrected pressure" that allows the ideal gas chemical potential formula to be applied to real gases. Activity extends this idea to condensed phases. Lewis and Randall's 1923 textbook "Thermodynamics" established these concepts as the standard working tools of chemical thermodynamics.

Pierre Duhem (1861–1916)

Duhem contributed to the mathematical foundations of chemical thermodynamics, particularly the constraint equation that bears his name (jointly with Gibbs). The Gibbs-Duhem equation establishes that the intensive properties of a system are not independently variable — a result with profound consequences for the thermodynamic consistency of experimental data and theoretical models. Duhem also contributed to the philosophy of science, arguing for the holistic nature of physical theories.

Simulation: Haber Process Equilibrium

The Haber process ($\text{N}_2 + 3\text{H}_2 \rightleftharpoons 2\text{NH}_3$) is the industrial synthesis of ammonia and a classic application of chemical potential and equilibrium thermodynamics. The reaction is exothermic ($\Delta H^\circ \approx -92$ kJ/mol), so the van't Hoff equation predicts that $K$ decreases with increasing temperature. However, higher temperatures are needed for reasonable reaction rates. The simulation below plots the equilibrium constant, NH$_3$ yield, and standard Gibbs energy as functions of temperature at several total pressures.

Haber Process: Equilibrium Composition vs Temperature

Python

Plots van't Hoff relation, K(T), NH3 mole fraction vs temperature at multiple pressures, and standard Gibbs energy for the Haber process N2 + 3H2 = 2NH3.

haber_equilibrium.py118 lines

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

# ──────────────────────────────────────────────────────────────
# Haber Process: Equilibrium Composition vs Temperature
# N2 + 3H2 <=> 2NH3
#
# Uses the van't Hoff equation to compute K(T), then solves
# for equilibrium mole fractions at a given total pressure.
# ──────────────────────────────────────────────────────────────

R = 8.314          # J/(mol K)
T_ref = 298.15     # K
DH0 = -92220.0     # J/mol  (standard enthalpy of reaction)
DS0 = -198.76      # J/(mol K) (standard entropy of reaction)
DG0_ref = DH0 - T_ref * DS0  # standard Gibbs energy at T_ref

# K at reference temperature
K_ref = np.exp(-DG0_ref / (R * T_ref))

# Temperature range
T_vals = np.linspace(300, 900, 500)

# van't Hoff: ln K(T) = ln K(T_ref) - (DH0/R)(1/T - 1/T_ref)
ln_K = np.log(K_ref) - (DH0 / R) * (1.0 / T_vals - 1.0 / T_ref)
K_vals = np.exp(ln_K)

# Solve equilibrium at several pressures
# Starting with stoichiometric mixture: 1 mol N2, 3 mol H2 => n_total_init = 4
# Let xi = extent of reaction (mol)
# n_N2 = 1 - xi, n_H2 = 3 - 3*xi, n_NH3 = 2*xi
# n_total = 4 - 2*xi
# Kp = (P_NH3)^2 / (P_N2 * P_H2^3)

pressures = [1.0, 10.0, 50.0, 200.0, 500.0]  # in bar
P_std = 1.0  # standard pressure in bar

fig, axes = plt.subplots(2, 2, figsize=(14, 11))
fig.patch.set_facecolor('#0a0a1a')
for ax in axes.flat:
    ax.set_facecolor('#0a0a1a')
    for spine in ax.spines.values():
        spine.set_edgecolor('#334155')
    ax.tick_params(colors='#94a3b8')

colors = ['#f97316', '#fb923c', '#fbbf24', '#a3e635', '#22d3ee']

# Panel 1: ln K vs 1/T
axes[0,0].plot(1000.0 / T_vals, ln_K, color='#f97316', linewidth=2.5)
axes[0,0].set_xlabel('1000/T  (1/K)', color='#94a3b8', fontsize=12)
axes[0,0].set_ylabel('ln K', color='#94a3b8', fontsize=12)
axes[0,0].set_title("van't Hoff Plot: ln K vs 1/T", color='white', fontsize=13)
axes[0,0].axhline(y=0, color='#64748b', linewidth=0.8, linestyle='--')

# Panel 2: K(T) on log scale
axes[0,1].semilogy(T_vals, K_vals, color='#fb923c', linewidth=2.5)
axes[0,1].set_xlabel('Temperature (K)', color='#94a3b8', fontsize=12)
axes[0,1].set_ylabel('K(T)', color='#94a3b8', fontsize=12)
axes[0,1].set_title('Equilibrium Constant vs Temperature', color='white', fontsize=13)
axes[0,1].axhline(y=1, color='#64748b', linewidth=0.8, linestyle='--')

# Panel 3: Mole fraction NH3 vs T at various pressures
for idx, P_total in enumerate(pressures):
    xi_eq = np.zeros_like(T_vals)
    for i, K in enumerate(K_vals):
        # Solve numerically: sweep xi from 0 to ~1
        xi_test = np.linspace(0.001, 0.999, 5000)
        n_N2 = 1.0 - xi_test
        n_H2 = 3.0 - 3.0 * xi_test
        n_NH3 = 2.0 * xi_test
        n_tot = 4.0 - 2.0 * xi_test
        P_N2 = (n_N2 / n_tot) * P_total
        P_H2 = (n_H2 / n_tot) * P_total
        P_NH3 = (n_NH3 / n_tot) * P_total
        Kp_calc = (P_NH3 / P_std)**2 / ((P_N2 / P_std) * (P_H2 / P_std)**3)
        diff = np.abs(np.log(Kp_calc + 1e-300) - np.log(K + 1e-300))
        best = np.argmin(diff)
        xi_eq[i] = xi_test[best]
    y_NH3 = 2.0 * xi_eq / (4.0 - 2.0 * xi_eq)
    axes[1,0].plot(T_vals, y_NH3, color=colors[idx], linewidth=2,
                   label=f'P = {P_total:.0f} bar')

axes[1,0].set_xlabel('Temperature (K)', color='#94a3b8', fontsize=12)
axes[1,0].set_ylabel('Mole Fraction NH3', color='#94a3b8', fontsize=12)
axes[1,0].set_title('NH3 Yield vs Temperature (Haber Process)', color='white', fontsize=13)
axes[1,0].legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[1,0].set_ylim(0, 1)

# Panel 4: DG_rxn vs T at P = 1 bar
DG_vals = DH0 - T_vals * DS0 + R * T_vals * np.log(1.0 / (K_vals + 1e-300)) * 0
# Actually DG0(T) approx = DH0 - T*DS0 (assuming constant DH, DS)
DG_approx = DH0 - T_vals * DS0
axes[1,1].plot(T_vals, DG_approx / 1000.0, color='#f97316', linewidth=2.5, label=r'$\Delta G^\circ(T)$')
axes[1,1].axhline(y=0, color='#ef4444', linewidth=1.0, linestyle='--', alpha=0.7)
axes[1,1].set_xlabel('Temperature (K)', color='#94a3b8', fontsize=12)
axes[1,1].set_ylabel(r'$\Delta G^\circ$ (kJ/mol)', color='#94a3b8', fontsize=12)
axes[1,1].set_title('Standard Gibbs Energy of Reaction', color='white', fontsize=13)
axes[1,1].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')

# Mark the crossover temperature
T_cross = DH0 / DS0
if 300 < T_cross < 900:
    axes[1,1].axvline(x=T_cross, color='#22d3ee', linewidth=1, linestyle=':', alpha=0.8)
    axes[1,1].annotate(f'T* = {T_cross:.0f} K', xy=(T_cross, 0),
                       xytext=(T_cross + 30, 10), fontsize=10, color='#22d3ee',
                       arrowprops=dict(arrowstyle='->', color='#22d3ee'))

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')
plt.close()
print("Haber process equilibrium analysis complete.")
print(f"K at 298 K: {K_ref:.2e}")
print(f"K at 500 K: {np.exp(np.log(K_ref) - (DH0/R)*(1/500 - 1/T_ref)):.2e}")
print(f"K at 700 K: {np.exp(np.log(K_ref) - (DH0/R)*(1/700 - 1/T_ref)):.2e}")
print(f"Crossover temperature (DG=0): {T_cross:.1f} K")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Computation: Equilibrium Constants for Multiple Reactions

The Fortran program below computes equilibrium constants and thermodynamic quantities for several industrially important reactions (Haber, hydrogen iodide, water-gas shift, sulfur trioxide formation) across a range of temperatures. It also performs a Gibbs-Duhem consistency check for a binary ideal mixture, verifying that $\sum x_i \, d\mu_i = 0$ at constant $T$and $P$.

Equilibrium Constants & Gibbs-Duhem Verification

Fortran

Computes K(T), DG0, and ln K for N2+3H2=2NH3, H2+I2=2HI, CO+H2O=CO2+H2, and 2SO2+O2=2SO3. Also verifies the Gibbs-Duhem equation for a binary ideal mixture.

chemical_potential_eq.f90128 lines

program chemical_potential_equilibrium
  implicit none
  ! ──────────────────────────────────────────────────────────────
  ! Compute equilibrium constants and compositions for several
  ! reactions at various temperatures using van't Hoff equation.
  ! ──────────────────────────────────────────────────────────────

integer, parameter :: dp = selected_real_kind(15, 307)
  real(dp), parameter :: R = 8.314_dp    ! J/(mol K)
  real(dp), parameter :: T_ref = 298.15_dp

integer, parameter :: n_reactions = 4
  integer, parameter :: n_temps = 8

character(len=40) :: rxn_names(n_reactions)
  real(dp) :: DH0(n_reactions), DS0(n_reactions)
  real(dp) :: T_vals(n_temps)
  real(dp) :: K_val, DG0, ln_K_ref, ln_K
  integer :: i, j

! Reaction data: DH0 (J/mol), DS0 (J/(mol K))
  rxn_names(1) = 'N2 + 3H2 <=> 2NH3'
  DH0(1) = -92220.0_dp;  DS0(1) = -198.76_dp

rxn_names(2) = 'H2 + I2 <=> 2HI'
  DH0(2) = -10170.0_dp;  DS0(2) = 21.81_dp

rxn_names(3) = 'CO + H2O <=> CO2 + H2'
  DH0(3) = -41166.0_dp;  DS0(3) = -42.08_dp

rxn_names(4) = '2SO2 + O2 <=> 2SO3'
  DH0(4) = -197780.0_dp; DS0(4) = -187.90_dp

! Temperature values
  T_vals(1) = 300.0_dp;  T_vals(2) = 400.0_dp
  T_vals(3) = 500.0_dp;  T_vals(4) = 600.0_dp
  T_vals(5) = 700.0_dp;  T_vals(6) = 800.0_dp
  T_vals(7) = 900.0_dp;  T_vals(8) = 1000.0_dp

write(*,'(A)') '================================================================='
  write(*,'(A)') '  Chemical Potential: Equilibrium Constants vs Temperature'
  write(*,'(A)') '================================================================='
  write(*,*)

do i = 1, n_reactions
    write(*,'(A,A)') 'Reaction: ', trim(rxn_names(i))
    write(*,'(A,F10.1,A,F8.2,A)') '  DH0 = ', DH0(i)/1000.0_dp, &
         ' kJ/mol,  DS0 = ', DS0(i), ' J/(mol K)'

! Crossover temperature
    if (abs(DS0(i)) > 1.0e-10_dp) then
      write(*,'(A,F8.1,A)') '  Crossover T (DG=0): ', DH0(i)/DS0(i), ' K'
    end if

write(*,'(A)') '  -------------------------------------------------'
    write(*,'(A)') '    T(K)       DG0(kJ/mol)      ln K          K'
    write(*,'(A)') '  -------------------------------------------------'

ln_K_ref = -(DH0(i) - T_ref * DS0(i)) / (R * T_ref)

do j = 1, n_temps
      ! van't Hoff equation
      ln_K = ln_K_ref - (DH0(i) / R) * (1.0_dp / T_vals(j) - 1.0_dp / T_ref)
      DG0 = DH0(i) - T_vals(j) * DS0(i)
      K_val = exp(ln_K)

if (abs(K_val) < 1.0e12_dp .and. abs(K_val) > 1.0e-12_dp) then
        write(*,'(F8.1,F16.3,F14.4,ES16.4)') &
             T_vals(j), DG0/1000.0_dp, ln_K, K_val
      else
        write(*,'(F8.1,F16.3,F14.4,A)') &
             T_vals(j), DG0/1000.0_dp, ln_K, '    (extreme)'
      end if
    end do

write(*,*)
    write(*,'(A)') '-----------------------------------------------------------------'
    write(*,*)
  end do

! ── Gibbs-Duhem consistency check for binary ideal mixture ──
  write(*,'(A)') '================================================================='
  write(*,'(A)') '  Gibbs-Duhem Check: Binary Ideal Mixture at 500 K, 1 bar'
  write(*,'(A)') '================================================================='
  write(*,*)

block
    real(dp) :: x1, x2, mu1, mu2, T, sum_xdmu
    real(dp) :: dx, mu1_prev, mu2_prev
    T = 500.0_dp
    dx = 0.01_dp

write(*,'(A)') '    x1       x2      mu1/RT     mu2/RT    x1*dmu1+x2*dmu2'
    write(*,'(A)') '  -------------------------------------------------------'

mu1_prev = 0.0_dp
    mu2_prev = 0.0_dp

do j = 1, 99
      x1 = j * dx
      x2 = 1.0_dp - x1
      ! For ideal mixture: mu_i/RT = mu_i^0/RT + ln(x_i)
      ! We set mu_i^0 = 0 for simplicity
      mu1 = log(x1)
      mu2 = log(x2)

if (j > 1) then
        sum_xdmu = x1 * (mu1 - mu1_prev) + x2 * (mu2 - mu2_prev)
      else
        sum_xdmu = 0.0_dp
      end if

if (mod(j, 10) == 0) then
        write(*,'(F8.3,F8.3,F11.4,F11.4,ES16.4)') &
             x1, x2, mu1, mu2, sum_xdmu
      end if

mu1_prev = mu1
      mu2_prev = mu2
    end do
  end block

write(*,*)
  write(*,'(A)') 'At constant T,P the Gibbs-Duhem relation requires'
  write(*,'(A)') 'sum(x_i * dmu_i) = 0. Values near zero confirm consistency.'

end program chemical_potential_equilibrium

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Summary of Key Results

Result	Equation	Significance
Definition of $\mu_i$	$\mu_i = (\partial G / \partial n_i)_{T,P,n_j}$	Partial molar Gibbs energy
Euler equation	$G = \sum_i n_i \mu_i$	G reconstructed from intensive quantities
Gibbs-Duhem	$\sum_i n_i d\mu_i = -SdT + VdP$	Constrains intensive variables
Ideal gas $\mu_i$	$\mu_i = \mu_i^\circ + RT\ln(P_i/P^\circ)$	Basis for equilibrium calculations
Equilibrium constant	$\Delta G^\circ = -RT\ln K$	Links thermodynamic data to composition
van't Hoff	$d(\ln K)/dT = \Delta H^\circ/(RT^2)$	Temperature dependence of K

← Part II: Thermodynamics Part II Overview →