Solving the Gibbs Free Energy Minimization Problem with Python

The Water-Gas Shift Equilibrium

Why Minimizing Gibbs Free Energy Matters

Chemical systems don’t settle into equilibrium by “reacting until reactants run out.” They settle where the total Gibbs free energy of the system is at its minimum at fixed temperature $T$ and pressure $P$. This is one of the most elegant ideas in thermodynamics: instead of solving equilibrium constant equations one reaction at a time, you can treat the whole problem as a constrained optimization problem.

For a mixture of $N$ species, the total Gibbs free energy is:

$$
G(n_1, n_2, \dots, n_N) = \sum_{i=1}^{N} n_i \left[ \mu_i^{\circ}(T) + RT \ln\left(\frac{y_i P}{P^{\circ}}\right) \right]
$$

where:

  • $n_i$ = moles of species $i$
  • $y_i = n_i / \sum_j n_j$ = mole fraction
  • $\mu_i^{\circ}(T)$ = standard chemical potential of species $i$ at temperature $T$
  • $R$ = universal gas constant
  • $P^{\circ}$ = reference pressure (1 atm)

The equilibrium composition is the set of $n_i$ that minimizes $G$, subject to atomic mass balance constraints:

$$
\sum_{i=1}^{N} a_{ki} , n_i = b_k \quad \text{for each atom } k
$$

where $a_{ki}$ is the number of atoms of element $k$ in species $i$, and $b_k$ is the total amount of element $k$ available from the feed. This formulation (the “direct minimization” method, originally due to White, Johnson & Dantzig) avoids ever writing down an equilibrium constant — the physics falls out entirely from the shape of $G$.

The Example: Water-Gas Shift Reaction

We’ll use a classic industrial reaction:

$$
\text{CO} + \text{H}_2\text{O} ;\rightleftharpoons; \text{CO}_2 + \text{H}_2
$$

This is a single-degree-of-freedom system: 4 species (CO, H₂O, CO₂, H₂), 3 atoms (C, H, O), so there’s exactly one independent reaction extent $\xi$. Starting from 1 mol CO + 1 mol H₂O, we can write:

$$
n_{CO} = 1-\xi,\quad n_{H_2O} = 1-\xi,\quad n_{CO_2} = \xi,\quad n_{H_2} = \xi
$$

We’ll solve the general $n$-dimensional problem with scipy.optimize.minimize (SLSQP), then also verify it against the simple 1-D $\xi$ picture, scan over temperature, and finally map out the equilibrium H₂ mole fraction as a 3D surface over $(T, P)$.

Full Source Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
# ============================================================
#  Gibbs Free Energy Minimization: Water-Gas Shift Equilibrium
# ============================================================
import numpy as np
from scipy.optimize import minimize, minimize_scalar
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 (enables 3D projection)

# ---------- 1. Thermodynamic data ----------
R = 8.314462618  # J/(mol*K), universal gas constant

species = ['CO', 'H2O', 'CO2', 'H2']

# Standard enthalpy of formation [J/mol]
dHf = np.array([-110.53e3, -241.82e3, -393.51e3, 0.0])

# Standard molar entropy [J/(mol*K)]
S0 = np.array([197.66, 188.83, 213.79, 130.68])

# Atom balance matrix (rows: C, H, O / columns: CO, H2O, CO2, H2)
A = np.array([
    [1, 0, 1, 0],   # Carbon
    [0, 2, 0, 2],   # Hydrogen
    [1, 1, 2, 0],   # Oxygen
], dtype=float)

# Feed composition [mol]: 1 mol CO + 1 mol H2O
feed = np.array([1.0, 1.0, 0.0, 0.0])
b = A.dot(feed)  # total atoms available (conserved quantities)


def mu0(T):
    """Approximate standard chemical potential of each species at temperature T [K]."""
    return dHf - T * S0


def gibbs_total(n, T, P, P0=1.0):
    """Total Gibbs free energy of the ideal-gas mixture [J]."""
    n = np.maximum(n, 1e-12)          # guard against log(0)
    n_tot = n.sum()
    y = n / n_tot                      # mole fractions
    mu = mu0(T) + R * T * np.log(y * P / P0)
    return float(np.sum(n * mu))


def solve_equilibrium(T, P, x0=None):
    """Minimize G(n) under atom-balance constraints -> equilibrium composition."""
    if x0 is None:
        x0 = feed + 0.3

    constraints = [
        {'type': 'eq', 'fun': (lambda n, i=i: A[i].dot(n) - b[i])}
        for i in range(A.shape[0])
    ]
    bounds = [(1e-9, None)] * len(species)

    result = minimize(
        gibbs_total, x0, args=(T, P),
        method='SLSQP', bounds=bounds, constraints=constraints,
        options={'maxiter': 300, 'ftol': 1e-12}
    )
    return result


# ---------- 2. Single-point equilibrium (T=1000K, P=1atm) ----------
T_ref, P_ref = 1000.0, 1.0
res_ref = solve_equilibrium(T_ref, P_ref)

print("=== Equilibrium at T=1000K, P=1atm ===")
print("Optimizer success:", res_ref.success, "|", res_ref.message)
for s, n in zip(species, res_ref.x):
    print(f"  n_{s:4s} = {n:.5f} mol   (y_{s} = {n/res_ref.x.sum():.5f})")
print(f"  G_min = {res_ref.fun:.3f} J\n")


# ---------- 3. 1-D view: G vs. reaction extent xi ----------
def n_from_xi(xi):
    return np.array([1 - xi, 1 - xi, xi, xi])

def G_of_xi(xi, T, P):
    return gibbs_total(n_from_xi(xi), T, P)

xi_grid = np.linspace(1e-4, 1 - 1e-4, 400)
G_grid = np.array([G_of_xi(x, T_ref, P_ref) for x in xi_grid])

opt_xi = minimize_scalar(G_of_xi, bounds=(1e-6, 1 - 1e-6),
                          method='bounded', args=(T_ref, P_ref))

fig1, ax1 = plt.subplots(figsize=(7, 5))
ax1.plot(xi_grid, G_grid / 1000, lw=2, color='tab:blue', label=r'$G(\xi)$')
ax1.scatter([opt_xi.x], [opt_xi.fun / 1000], color='red', zorder=5,
            label=f'Minimum ($\\xi$={opt_xi.x:.3f})')
ax1.set_xlabel(r'Reaction extent $\xi$')
ax1.set_ylabel('Total Gibbs free energy  G  [kJ]')
ax1.set_title('Gibbs Free Energy vs. Reaction Extent (T=1000K, P=1atm)')
ax1.legend()
ax1.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('gibbs_vs_xi.png', dpi=150)
plt.show()


# ---------- 4. Equilibrium composition vs. Temperature (P fixed) ----------
T_scan = np.linspace(400, 1600, 40)
comp_scan = np.zeros((len(T_scan), 4))
x0 = feed + 0.3
for i, T in enumerate(T_scan):
    r = solve_equilibrium(T, P_ref, x0=x0)
    comp_scan[i] = r.x
    x0 = r.x  # warm start: reuse previous solution as next guess (much faster)

y_scan = comp_scan / comp_scan.sum(axis=1, keepdims=True)

fig2, ax2 = plt.subplots(figsize=(7, 5))
for j, s in enumerate(species):
    ax2.plot(T_scan, y_scan[:, j], marker='o', ms=3, label=s)
ax2.set_xlabel('Temperature [K]')
ax2.set_ylabel('Equilibrium mole fraction  y')
ax2.set_title('Water-Gas Shift Equilibrium Composition vs. Temperature (P=1atm)')
ax2.legend()
ax2.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('composition_vs_T.png', dpi=150)
plt.show()


# ---------- 5. 3-D surface: y_H2 as a function of T and P (warm-start speed-up) ----------
T_range = np.linspace(400, 1600, 25)
P_range = np.linspace(0.5, 5.0, 25)
TT, PP = np.meshgrid(T_range, P_range)
YH2 = np.zeros_like(TT)

x0_seed = feed + 0.3
for j in range(len(P_range)):
    x0_row = x0_seed.copy()
    for i in range(len(T_range)):
        r = solve_equilibrium(TT[j, i], PP[j, i], x0=x0_row)
        y = r.x / r.x.sum()
        YH2[j, i] = y[3]   # H2 is species index 3
        x0_row = r.x        # warm start along the T-axis for this P row

fig3 = plt.figure(figsize=(9, 7))
ax3 = fig3.add_subplot(111, projection='3d')
surf = ax3.plot_surface(TT, PP, YH2, cmap='viridis', edgecolor='none', antialiased=True)
ax3.set_xlabel('Temperature [K]')
ax3.set_ylabel('Pressure [atm]')
ax3.set_zlabel('Equilibrium  y_H2')
ax3.set_title('Equilibrium H2 Mole Fraction Surface: f(T, P)')
fig3.colorbar(surf, shrink=0.6, aspect=12, label='y_H2')
plt.tight_layout()
plt.savefig('h2_surface_3d.png', dpi=150)
plt.show()

print("All computations finished.")

Detailed Code Walkthrough

Section 1 — Thermodynamic data. We store enthalpy of formation ($\Delta H_f^\circ$) and standard entropy ($S^\circ$) for each of the 4 species, and approximate the temperature-dependent chemical potential as $\mu_i^\circ(T) \approx \Delta H_{f,i}^\circ - T S_i^\circ$. This ignores heat-capacity ($C_p$) corrections, which is a standard simplification for a first-pass teaching example — the qualitative equilibrium behavior it produces is realistic.

The matrix A encodes how many atoms of C, H, and O are in each species (columns = CO, H₂O, CO₂, H₂). Multiplying A by the feed vector gives the conserved atom totals b — this is the constraint the optimizer must respect no matter how the mixture rearranges itself.

gibbs_total(n, T, P) implements the free-energy formula directly: it clips n away from zero (to avoid log(0)), computes mole fractions y, and sums $n_i \mu_i$ over all species.

solve_equilibrium(T, P, x0) is the core optimizer call. It builds 3 equality constraints (one per atom) as lambda closures — note the i=i trick, which is required so each lambda captures its own value of i rather than all sharing the last loop value. Bounds keep every $n_i$ strictly positive (no negative moles). SLSQP (Sequential Least Squares Programming) is used because it natively supports both bounds and nonlinear/linear equality constraints — ideal for this kind of problem.

Section 3 cross-checks the general optimizer against the simple 1-D reaction-extent formulation, where $\xi$ alone fully parametrizes the mixture. Scanning $\xi$ from 0 to 1 and plotting $G(\xi)$ should show a single smooth bowl-shaped curve with one minimum — exactly where minimize_scalar lands, and exactly where the SLSQP solution from Section 2 should agree.

Section 4 solves equilibrium repeatedly across a temperature sweep. The key performance trick is warm-starting: each solve’s result becomes the initial guess (x0) for the next temperature, since equilibrium composition changes only slightly between adjacent grid points. This cuts iteration counts dramatically compared to always starting from the same generic guess.

Section 5 extends the same warm-start idea to a full $(T, P)$ grid (25×25 = 625 optimizations), reusing the previous $T$-point’s solution as it sweeps across each pressure row. This keeps the whole grid scan fast even though it’s calling scipy.optimize.minimize hundreds of times.

=== Equilibrium at T=1000K, P=1atm ===
Optimizer success: True | Optimization terminated successfully
  n_CO   = 0.51293 mol   (y_CO = 0.25646)
  n_H2O  = 0.51293 mol   (y_H2O = 0.25646)
  n_CO2  = 0.48707 mol   (y_CO2 = 0.24354)
  n_H2   = 0.48707 mol   (y_H2 = 0.24354)
  G_min = -761468.144 J

All computations finished.

Reading the Graphs

Figure 1 — G vs. reaction extent ξ. This is the most direct visual proof of “free energy minimization” in action: the curve should form a single smooth bowl between $\xi=0$ and $\xi=1$, with the red marker sitting exactly at the bottom. There is no equilibrium constant, no algebra — just literally the lowest point of a curve.

Figure 2 — Composition vs. Temperature. Since $\Delta H_{rxn} = \Delta H_{f,CO_2} + \Delta H_{f,H_2} - \Delta H_{f,CO} - \Delta H_{f,H_2O} \approx -41\text{ kJ/mol}$, the water-gas shift reaction is mildly exothermic. By Le Chatelier’s principle, raising temperature should push the equilibrium back toward the reactants — so you should see $y_{H_2}$ and $y_{CO_2}$ decrease, while $y_{CO}$ and $y_{H_2O}$ increase, as $T$ rises across the scan. This is a genuine physical prediction falling directly out of the optimizer, not hard-coded.

Figure 3 — 3D surface of y_H₂ over (T, P). Look closely at the pressure axis: because this reaction has no change in total gas moles (2 moles of reactants → 2 moles of products), the equilibrium composition is theoretically independent of pressure — the surface should look essentially flat along the $P$ direction while curving strongly along $T$. That’s a subtle but important thermodynamic signature this graph should reveal at a glance: pressure only matters for reactions that change the mole count.

Wrap-Up

What makes this approach powerful is that it generalizes far beyond one reaction: add more species and more atom-balance rows, and the exact same solve_equilibrium function handles arbitrarily complex reacting systems — combustion, gasification, multi-phase equilibria — without ever writing a single equilibrium-constant expression by hand. The physics is entirely encoded in the shape of $G(n)$, and scipy.optimize.minimize does the rest.