Lab 01: Thermodynamics · From G–X curves to equilibrium phase fractions¶

Thermodynamics and Kinetics of Materials · NPTEL course companion
Lab 01 · Thermodynamics — dual energy axes, explicit SER and three-phase invariant equilibrium.

Level: undergraduate materials science · Suggested time: 75–90 minutes

You have learned that, when multiple phases are in equilibrium, the chemical potential of each component must be the same in every phase in equilibrium. For component i, we write μᵢᵅ = μᵢᵝ = μᵢᵞ. Here, i denotes the component, and α, β and γ label the phases in equilibrium. We compare the chemical potential of the same component across the phases; we do not require different components to have the same chemical potential. Chemical potential is the change in Gibbs energy when we add one mole of a component, keeping temperature, pressure and the amounts of the other components fixed. Its unit is J/mol.

The phases can have different compositions even though these chemical potentials are equal. Each phase has its own relation between Gibbs energy and composition. We use these relations to find the compositions at which the chemical potentials of every component are equal across the phases. For a binary system, this means satisfying the equilibrium condition for both component A and component B.

We must also conserve each component. In this lab, the system is closed and there are no chemical reactions, so the total number of moles of component i remains unchanged. If nᵅ is the number of moles in phase α and xᵢᵅ is its mole fraction of component i, then nᵅxᵢᵅ + nᵝxᵢᵝ + nᵞxᵢᵞ = nᵢ,total. Each term is an amount in mol; a mole fraction is dimensionless. This balance tells us how the component is distributed among the phases.

Together, the equilibrium conditions and component balances determine the permitted equilibrium compositions and phase amounts. The common tangent identifies coexisting compositions on the molar Gibbs energy–composition (G–X) curves. Once two-phase compositions are known, the component balance gives the phase amounts. At a three-phase invariant equilibrium, overall composition alone does not uniquely determine all three phase amounts; an additional constraint is needed. In this lab, we show each step using G–X curves. Molar Gibbs energy has unit J/mol, and composition is expressed as a dimensionless mole fraction.

You will draw tangents, read chemical potentials from intercepts, compute a common tangent, and use conservation to obtain phase amounts. Predict each result before running its cell. Work through the seven exercises at your own pace. Record your predictions, plots and explanations, then compare them with the worked solutions and self-checks.

Prerequisites: mole fraction, molar Gibbs energy, derivatives and binary phase equilibria. Only basic Python is needed. This is an original, fictional A–B model, not an assessed alloy database or a prediction for a real material.

Learn to build the code: the final appendix explains the implementation line by line, including NumPy arrays, zip, itertools, solvers and Matplotlib. The core laboratory takes about 75–90 minutes; work through the coding guide separately at your own pace.

1. Setup and run¶

Use Python 3.10 or newer and JupyterLab. In a terminal, create an environment:

python -m venv .venv
# Linux/macOS:
source .venv/bin/activate
# Windows PowerShell instead:
# .venv\Scripts\Activate.ps1
python -m pip install numpy scipy matplotlib jupyterlab
python -m jupyter lab

Open this notebook and choose Kernel → Restart Kernel and Run All Cells. Alternatively, upload the notebook to Google Colab and run it; install missing packages there with %pip install numpy scipy matplotlib. No database, API key or internet connection is needed after installation. Saved outputs let you inspect the reference run without executing anything.

The parameter cell in §2 is the place to experiment. After a change, rerun from that cell downward. The self-check section restores the baseline parameters so you can compare your results with a reproducible reference calculation.

2. A self-contained two-phase binary¶

At fixed temperature $T$ and pressure, let $x=n_B/(n_A+n_B)$ and $g=G/(n_A+n_B)$, in J mol$^{-1}$. There are two ideal-solution phases, solid $\alpha$ and liquid $L$, with the same composition variable and consistent pure-component references:

$$g^\phi(x,T)=(1-x)g_A^{0,\phi}+xg_B^{0,\phi} +RT[(1-x)\ln(1-x)+x\ln x].$$

The reference energies are intentionally constant with temperature for this small teaching model. A prefers $\alpha$, B prefers $L$; configurational entropy favors mixing. We neglect interfaces, elasticity, finite-size effects, kinetics and all other phases. Phase fractions below are mole fractions of material in each phase, not mass or volume fractions.

Use $0\ln0=0$ for energy at pure endpoints. Chemical potentials have logarithmic limits there, so tangents are evaluated only for $0<x<1$. Each branch is convex: $g''=RT/[x(1-x)]>0$. Two-phase equilibrium here arises from two different branches, not a spinodal instability of one branch.

Reference convention in this first calculation: these four constants share one reference for A and one for B across both phases; they are not a table of SER data. The zeros at pure A in solid and pure B in liquid are a convenient choice at the demonstration temperature. They do not authorize separately zeroing every phase's endmembers. Section 6 constructs explicit SER functions and converts between references. The temperature sweep of these frozen constants is only an entropy sensitivity experiment, not a liquid–solid phase diagram.

3. The tangent tells us both chemical potentials¶

At fixed $T,P$, $G=ng(x)$ with $n=n_A+n_B$. Differentiating with respect to each component amount gives $$\mu_A=\left(\frac{\partial G}{\partial n_A}\right)_{T,P,n_B}=g-xg',\qquad \mu_B=\left(\frac{\partial G}{\partial n_B}\right)_{T,P,n_A}=g+(1-x)g'.$$ Thus $g'=\mu_B-\mu_A$, not either individual chemical potential. The tangent through $x_0$ is $$\ell(x)=g(x_0)+g'(x_0)(x-x_0)=(1-x)\mu_A(x_0)+x\mu_B(x_0).$$ Its extrapolated values at $x=0$ and $x=1$ give $\mu_A$ and $\mu_B$. These are not generally the pure-component energies.

Expected plot: one convex $\alpha$ curve and its tangent, touching at $x=0.35$. The left and right tangent intercepts are explicitly marked. Moving $x_0$ changes both intercepts.

Both vertical axes use exactly the same energy scale. The left axis is at $x_B=0$ and the right at $x_B=1$; the colored points are the tangent intercepts, not the endmember energies.

4. A common tangent enforces equilibrium¶

For coexisting compositions $x_\alpha,x_L$, $$\mu_A^\alpha(x_\alpha)=\mu_A^L(x_L),\qquad \mu_B^\alpha(x_\alpha)=\mu_B^L(x_L).$$ Equivalently, the tangents have the same slope and intercept: $$g_\alpha'(x_\alpha)=g_L'(x_L)=m,\quad g_\alpha(x_\alpha)-mx_\alpha=g_L(x_L)-mx_L=b.$$ Parallel tangents alone are insufficient. Also require a supporting line below all allowed energy branches; otherwise a stationary construction need not be globally stable.

For a trial slope $m$, our ideal-solution formula can be inverted: $$x_\phi(m)=\frac{1}{1+\exp[(g_B^{0,\phi}-g_A^{0,\phi}-m)/(RT)]}.$$ We vary $m$ until the two tangent intercepts coincide. brentq is a bracketed scalar root solver. expit and logaddexp evaluate the expressions without unnecessary exponential overflow. This inversion is specific to the ideal model, not a universal CALPHAD solver.

5. Equilibrium amounts: the lever rule¶

Let $z$ be the overall B mole fraction, with $f_\alpha+f_L=1$. Conservation requires $$z=f_\alpha x_\alpha+f_L x_L.$$ For $x_\alpha<z<x_L$, $$f_L=\frac{z-x_\alpha}{x_L-x_\alpha},\qquad f_\alpha=\frac{x_L-z}{x_L-x_\alpha},\quad g_{\rm eq}(z)=f_\alpha g_\alpha(x_\alpha)+f_L g_L(x_L).$$ Outside this interval the stable state is a single phase at composition $z$; do not extrapolate the lever rule to negative fractions. At a boundary the amount of the incipient phase is zero. The common tangent extended outside the interval is not an attainable mixture of these two endpoint compositions.

Expected plot: two convex curves, a dashed supporting line and a thick lower equilibrium envelope. Inside the tie line the envelope is straight. The marked energy saving at $z$ compares equilibrium to the lower homogeneous branch at the same overall composition. It is an energy per mole, not a kinetic rate or activation barrier.

Expected plots below: phase fractions vary linearly only inside the two-phase interval; equilibrium chemical potentials remain constant across that interval even while $z$ and phase amounts change. Outside it, chemical potentials follow the appropriate homogeneous phase. Endpoints are omitted from the chemical-potential panel because one component is infinitely dilute.

6. Explicit SER example: solid α and liquid for each species¶

The reference belongs to the species, not to the phase being evaluated. A and B may have different stable crystal structures in a real database. For this fictional example, both stable elemental phases at $T_r=298.15$ K and $P_r=1$ bar are called α.

SER (standard/stable element reference) conventionally uses the enthalpy $H_i^{\rm SER}(T_r,P_r)$ of each element in that state. Write the tabulated quantity as $$\bar g_i^{0,\phi}(T)=G_i^{0,\phi}(T,P_r)-H_i^{\rm SER}(T_r,P_r),\quad \phi\in\{\alpha,L\}.$$ The same $H_A^{\rm SER}$ is used for A in α and in liquid, and the same $H_B^{\rm SER}$ for B in both phases. The bars emphasize the energy reference; our earlier $g$ and $\mu$ are likewise always reference-dependent values.

SER does not mean $G_i^{\rm SER}(298.15)=0$. At the reference temperature, $$\bar g_i^{0,\alpha}(T_r)=-T_r S_i^\alpha(T_r),$$ which is generally nonzero. The entropy at 298.15 K is not arbitrarily erased. The reference state is fixed; it does not switch to liquid when the calculation temperature crosses the melting point.

Build original unary functions, rather than copy a database¶

Use constant heat capacities, with hypothetical parameters:

Integrating $C_p$ from $T_r$ and retaining the entropy term gives $$\bar g_i^{0,\alpha}(T)=C_{p,i}(T-T_r)-T\left[S_i^\alpha(T_r)+C_{p,i}\ln(T/T_r)\right].$$ With equal solid/liquid heat capacities, the fusion difference is $$\Delta g_{fus,i}(T)=\Delta H_{fus,i}\left(1-\frac{T}{T_{m,i}}\right),\qquad \bar g_i^{0,L}(T)=\bar g_i^{0,\alpha}(T)+\Delta g_{fus,i}(T).$$ Both phases of a species meet at its melting temperature. Above it, that species' liquid is lower in Gibbs energy. These functions are only an instructional constant-$C_p$ model, not a fit for real A or B or a low-temperature extrapolation to 0 K.

At 800 K the fusion differences are +8000 J/mol for A and −12000 J/mol for B. The four original constants can therefore be recovered by subtracting one shared reference per species: $r_A=\bar g_A^{0,\alpha}(800)$, $r_B=\bar g_B^{0,L}(800)$. This connects the earlier illustrative table to a consistent SER construction.

The same equilibrium in two reference conventions¶

For any chosen species references $r_A(T)$ and $r_B(T)$, applied identically to every phase, $$\widetilde g^\phi(x,T)=\bar g^\phi(x,T)-[(1-x)r_A(T)+xr_B(T)].$$ Then $\widetilde\mu_A^\phi=\bar\mu_A^\phi-r_A$ and $\widetilde\mu_B^\phi=\bar\mu_B^\phi-r_B$. Since each species receives the same shift in each phase, equality of chemical potentials and the equilibrium phase amounts are unchanged. A temperature-dependent reference is acceptable for this fixed-$T$ equilibrium comparison; derivatives used for entropy and enthalpy must also include the reference's temperature dependence.

In the right-hand plot choose $r_i(T)=\bar g_i^{0,\alpha}(T)$. This makes both pure-solid endmembers zero, while retaining the liquid–solid differences. This is a common solid-Gibbs reference at the calculation temperature, not the SER enthalpy convention.

Expected plots: the left uses SER-referenced energies; the right uses pure-solid Gibbs references. Their slopes and intercepts differ, but both tangency compositions and phase fractions agree. Both left/right vertical axes show the intercepts on identical scales within each panel. Numerical scales need not match between panels because the references differ.

If instead we subtract each phase's own pure endmembers, both ideal curves collapse to the same $RT[x\ln x+(1-x)\ln(1-x)]$. That is useful for displaying mixing energies, but discards the fusion differences and cannot be used alone to determine solid–liquid equilibrium.

7. Three-phase equilibrium: why a binary eutectic is invariant¶

At fixed pressure the Gibbs phase rule is $$F=C-P+1.$$ For a binary system with three phases ($\alpha+L+\beta$), $F=2-3+1=0$: the equilibrium temperature and the three phase compositions are fixed. A change in temperature destroys three-phase coexistence. Pressure is held at 1 bar here; without that constraint the full rule would give $F=C-P+2=1$.

Use a separate illustrative quadratic model, not an extrapolation of the earlier ideal/SER model. All phases share the same species reference convention. Let $T_E=1000$ K, $$g^\alpha=K(x-0.15)^2,\quad g^\beta=K(x-0.85)^2,\quad g^L=K(x-0.50)^2-\Delta s(T-T_E),$$ where $K=20000$ J/mol and $\Delta s=20$ J/(mol K). These are local teaching approximations, not realistic pure-element limits or assessed alloy data. No entropy of ideal mixing is added to these parabolas.

At $T_E$, one horizontal supporting tangent touches all three curves: $$x_\alpha=0.15,\quad x_L=0.50,\quad x_\beta=0.85,\qquad \mu_A^\alpha=\mu_A^L=\mu_A^\beta=0,\quad \mu_B^\alpha=\mu_B^L=\mu_B^\beta=0.$$ The zero values are a convenient shared energy reference, not a requirement for an invariant. Adding a common affine reference term would tilt this same three-contact tangent and leave the compositions unchanged.

Expected plots: below $T_E$, the stable common tangent connects α and β; at $T_E$ it also touches liquid; above $T_E$, separate α–L and L–β tangents occur, with a single-liquid interval between them. Each panel has identical energy scales on its left/right axes. At $T_E$ the colored endpoint markers identify the common chemical potentials.

For parabolas $g_i=K(x-c_i)^2+h_i$, the pair-tangent slope is $m=(h_j-h_i)/(c_j-c_i)$ and the contact compositions are $x_i=c_i+m/(2K)$, $x_j=c_j+m/(2K)$. A candidate is kept only if its contacts are in [0,1] and it lies below all three branches.

Invariant compositions do not uniquely fix the phase amounts¶

The phase rule counts intensive degrees of freedom, not the amounts of material in the phases. At $T_E$ with overall composition $z$, $$f_\alpha+f_L+f_\beta=1,\qquad 0.15f_\alpha+0.50f_L+0.85f_\beta=z.$$ There are only two independent equations for three unknown fractions. Choose $\lambda=f_L$; then $$f_\beta=\frac{z-0.15-0.35\lambda}{0.70},\qquad f_\alpha=1-\lambda-f_\beta.$$ All nonnegative solutions have the same equilibrium Gibbs energy on the common tangent. An additional condition, such as total enthalpy during a thermal arrest, can fix the amounts. $F=0$ does not remove this freedom in phase amounts.

At $z=0.40$, $0\le f_L\le5/7$. Interior values give all three phases; the limiting values have one absent phase. The endpoint states are limits of the three-phase assemblage, not states with three positive fractions.

The eutectic reaction is $L(0.50)\rightleftharpoons\alpha(0.15)+\beta(0.85)$. For one mole of liquid at its eutectic composition, complete freezing gives 0.5 mol α and 0.5 mol β. At $T_E$, partial freezing can coexist with the remaining liquid without changing any of these phase compositions.

8. Try it yourself¶

Write a prediction first, then use the code. Keep units and distinguish $z$ from $x_\alpha,x_L$.

1. Intercepts. At $T=800$ K, $x_0=0.35$ on $\alpha$, calculate both chemical potentials and the slope. Explain why the slope alone cannot give both potentials.
2. Equilibrium amounts. At 800 K, find $x_\alpha,x_L$ and the fractions for $z=0.20,0.40,0.70$. Which quantities change with $z$? Verify the B balance for one case. Then classify $z=0.05,0.90$.
3. The crossing trap. Find where $g_\alpha(x)=g_L(x)$. At that composition, are the slopes equal? Compare the homogeneous energy with the equilibrium energy. Explain why equal $g$ at the same composition is insufficient.
4. Temperature experiment. Repeat the common-tangent calculation at 600, 800 and 1200 K, holding $z=0.40$. Tabulate both phase compositions and fractions. Does the tie-line interval widen or narrow? Explain using mixing entropy, within this model's assumptions.
5. Reference-energy challenge. Add the same affine term $c+dx$ to both branches, using $c=1500$, $d=-2000$ J/mol. Predict changes in compositions, fractions, common-tangent slope and both chemical potentials. Explain why relative phase energies, rather than a common reference choice, determine equilibrium.

Check your understanding: annotate your common-tangent and phase-fraction plots, complete the temperature table, and explain each result in your own words. Try the exercises before consulting the worked solutions.

6. SER and phase references. At 800 K, reproduce the four SER-referenced pure-element energies. Why is the solid Gibbs energy not zero at 298.15 K? Convert both phases to the common pure-solid Gibbs reference. Explain why zeroing liquid endmembers independently would erase physical information. Identify the fusion-energy difference for each species and explain its sign at 800 K.

7. Binary invariant. Apply the phase rule at fixed pressure to α + L + β. For $z=0.40$, calculate all three fractions first with $f_L=0.20$ and then with $f_L=0.60$. Explain why both states are equilibrium states at the same invariant temperature and compositions. What additional measurement could determine the fractions uniquely? Use the neighboring-temperature plots to explain why three-phase coexistence is lost when temperature changes.

9. Worked solutions and self-checks¶

Compare your reasoning and calculations with the worked solutions below. The checks restore the original model parameters. If a result differs, revisit the relevant equation, units and reference convention, then rerun the calculation.

An independent analytic result follows directly from the ideal-solution chemical potentials. Define $$K_A=\exp[(g_A^{0,\alpha}-g_A^{0,L})/(RT)],\quad K_B=\exp[(g_B^{0,\alpha}-g_B^{0,L})/(RT)].$$ Equality of potentials gives $1-x_L=K_A(1-x_\alpha)$ and $x_L=K_Bx_\alpha$, hence $$x_\alpha=\frac{1-K_A}{K_B-K_A},\qquad x_L=K_Bx_\alpha.$$ This provides a check independent of the slope-root construction.

Interpretation of solutions.

1. The slope is $\mu_B-\mu_A$; the intercept supplies the missing information.

2. At fixed $T,P$, the two coexisting compositions and chemical potentials stay fixed as $z$ changes inside the tie line. Fractions change to conserve material. At $z=0.05$ only $\alpha$ is stable; at $z=0.90$ only $L$ is stable, each at composition $z$.

3. The curves cross at $x=0.4$, but their slopes differ by 20000 J/mol. Separating into different compositions lowers the total Gibbs energy. Equality of both chemical potentials is the equilibrium criterion.

4. The interval narrows with increasing $T$ over the stated range. Mixing becomes more favorable. The fixed reference energies are a teaching simplification: do not infer a real alloy phase diagram or a finite critical temperature from this model.

5. A common affine shift changes neither compositions nor fractions. It shifts $m$ by $d$, $\mu_A$ by $c$, and $\mu_B$ by $c+d$. Both homogeneous and equilibrium energy at $z$ shift by $c+dz$, so the energy saving is unchanged.

6. The SER example's printed table gives the four values. At 298.15 K, the solid values are −8944.5 J/mol (A) and −11926 J/mol (B), because $G-H^{SER}=-T_rS(T_r)$. Subtracting each species' solid Gibbs energy from both phases gives α: (0, 0), liquid: (+8000, −12000) J/mol at 800 K. A is below its 1400 K melting point and favors solid; B is above its 600 K melting point and favors liquid. Subtracting different phase-specific references would remove these differences and destroy the equilibrium comparison. The numerical check above confirms identical compositions/fractions and the correct chemical-potential shifts.

7. At fixed pressure, $F=2-3+1=0$. With $f_L=0.20$, $(f_\alpha,f_L,f_\beta)=(0.542857,0.20,0.257143)$; with $f_L=0.60$, they are $(0.342857,0.60,0.057143)$. Both sum to one and give $z=0.40$. Phase amounts are not counted by $F$. Total enthalpy (with known phase enthalpies), for example, supplies a third independent balance. Below $T_E$ liquid lies above the α–β tangent; above $T_E$ it lies below that line, so no single tangent can support all three phases.

10. Where this leads¶

The same equilibrium conditions underlie calculations with pycalphad, where phase models may contain assessed temperature-dependent reference energies, excess mixing and sublattices. This first notebook keeps the model visible and avoids redistributing any third-party thermodynamic database. A later laboratory can compare an openly licensed database calculation with this construction.

FiPy addresses the next question: how spatial fields evolve through diffusion and other PDEs. Equilibrium compositions alone do not tell us the time required, diffusion profiles, precipitate shapes or nucleation barriers. Those require additional physical models and parameters.

Course and next-step documentation:

• Thermodynamics and Kinetics of Materials — NPTEL
• Course teaching page and lecture notes
• pycalphad documentation
• FiPy documentation (NIST)

Model provenance and reuse¶

All numerical reference energies above were chosen for this original teaching example; none were extracted from a proprietary or published alloy database. No third-party figures or course text are reproduced. I permit free noncommercial teaching, learning and demonstrations with acknowledgment. Research (including academic research) and commercial use require my explicit written permission. Earlier MIT and GPL releases retain their existing permissions. Read the current code-use terms. My original teaching text and figures follow CC BY 4.0. Dependencies retain their own licenses.

Earlier MIT releases retain their existing permissions.

Reference-convention documentation (definitions only; no unary data copied):

• Thermo-Calc: stable element reference and reference enthalpy
• ESPEI: custom unary reference states, lattice stabilities and SER data

11. Build the code yourself: a line-by-line guide¶

This appendix starts from a new, empty notebook. Work through the small code blocks in order. They are displayed as examples so that reading this appendix does not overwrite the calculations above. To run one, copy it into a new Code cell and press Shift+Enter. Explanations and equations belong in Markdown cells. Save frequently, and use Restart Kernel and Run All Cells to check that your notebook does not depend on an old variable left in memory.

The sequence is deliberate: define the data → evaluate one energy → evaluate an array → derive chemical potentials → find equilibrium → check conservation → plot → extend to three phases. A plot is the final view of a calculation, not a substitute for checking it.

11.1 Imports: what is a library, a module and an alias?¶

import numpy as np                         # Numerical arrays and elementwise mathematics.
import matplotlib.pyplot as plt           # Figure creation and plotting functions.
from scipy.special import xlogy, expit     # Two specific numerical functions.
from scipy.optimize import brentq, linprog # Root finding and linear programming.
from itertools import combinations        # Unordered selections without repetition.

• import numpy as np imports the NumPy package and binds it to the short name np. np is a convention, not a reserved Python word. Thus np.log(...) means “call NumPy's logarithm.”
• matplotlib.pyplot is a module within Matplotlib. plt.subplots(...) is a function in that module. Importing matplotlib alone does not establish the name plt.
• from ... import ... brings the named functions into your current namespace, so write brentq(...), not scipy.brentq(...). Importing a function does not run it.
• itertools is part of Python's standard library: there is no separate pip install itertools. NumPy, SciPy and Matplotlib are external packages.
• # begins a comment. Python ignores the rest of that line; the comment explains the intention to the reader.

Try: type help(combinations) or help(plt.subplots) in a cell. In Jupyter, plt.subplots? also opens help. If an import fails, install the missing package in the environment used by the current kernel; %pip install numpy scipy matplotlib in a notebook targets that kernel more reliably than an unrelated terminal's pip.

11.2 Scalars, lists, tuples, dictionaries and arrays¶

R = 8.314462618                            # Gas constant, J/(mol K).
T_demo = 800.0                             # A floating-point scalar, in kelvin.
z_demo = 0.40                              # Overall mole fraction of B.
endmembers = {                            # Dictionary: phase name -> two energies.
    "alpha": (0.0, 12000.0),               # Tuple ordered as (g_A^0, g_B^0).
    "liquid": (8000.0, 0.0),
}
a, b = endmembers["alpha"]               # Unpack the two tuple entries into two names.
grid = np.linspace(0.0, 1.0, 1001)          # 1001 evenly spaced points, including endpoints.
print(a, b, grid.shape, grid[0], grid[-1])  # Inspect values, size and first/last elements.

= assigns a value; == compares two values. Python does not attach physical units to these numbers. Choose and document one unit system, then convert only when displaying results.

A list, such as [600., 800., 1200.], is an ordered, mutable collection. A tuple, such as (0., 12000.), is an ordered collection whose entries cannot be reassigned in place. A dictionary maps keys to values: endmembers["alpha"] retrieves a tuple using its string key. Names are case sensitive: "Alpha" is not "alpha".

a, b = ... is unpacking. The right-hand side must supply exactly two items. The names are local conventions: here a means the A endmember energy, not the α phase fraction. Later b may denote a tangent intercept inside a different function; use descriptive names when this risks confusion.

A NumPy array is designed for numerical operations. grid.shape is (1001,): a one-dimensional array with 1001 elements. Python indexing starts at zero; grid[-1] means the last element. np.linspace(0, 1, 1001) has 1000 intervals, so the spacing is 0.001.

python_list = [0.0, 0.5, 1.0]
numeric_array = np.array(python_list)
print(python_list * 2)                     # Repeats the list: six entries.
print(numeric_array * 2)                   # Multiplies each number: [0., 1., 2.].
print(numeric_array**2)                    # Elementwise squaring; ^ is NOT exponentiation.

In array expressions, multiplying an array by the scalar R*T_demo applies that scalar to every element. This is a simple example of broadcasting. There is no need to write a Python loop over all compositions.

11.3 Write one energy function and understand every line¶

def energy(xx, temp, phase):
    """Return molar Gibbs energy in J/mol, for a scalar or an array of x_B."""
    xx = np.asarray(xx, dtype=float)
    if temp <= 0 or np.any((xx < 0) | (xx > 1)):
        raise ValueError("Use T > 0 K and 0 <= x_B <= 1.")
    ga, gb = endmembers[phase]
    mixing = R*temp*(xlogy(xx, xx) + xlogy(1-xx, 1-xx))
    return (1-xx)*ga + xx*gb + mixing

print(energy(0.4, T_demo, "alpha"))
energy_curve = energy(grid, T_demo, "alpha")
print(energy_curve.shape)

The main notebook calls its equivalent function g; this appendix uses energy to make the role explicit. The initial constant-endmember model is the same. For a temperature-dependent model, obtain the endmembers at the requested temp before evaluating this expression.

11.4 Derivatives, chemical potentials and multiple return values¶

def energy_slope(xx, temp, phase):
    xx = np.asarray(xx, dtype=float)
    if temp <= 0 or np.any((xx <= 0) | (xx >= 1)):
        raise ValueError("Derivatives require T > 0 and 0 < x_B < 1.")
    ga, gb = endmembers[phase]
    return gb-ga + R*temp*np.log(xx/(1-xx))


def potentials(xx, temp, phase):
    slope = energy_slope(xx, temp, phase)
    gg = energy(xx, temp, phase)
    return gg-xx*slope, gg+(1-xx)*slope

mu_a, mu_b = potentials(0.35, T_demo, "alpha")
print(f"mu_A = {mu_a:.3f}; mu_B = {mu_b:.3f} J/mol")

The derivative uses the open interval $(0,1)$: unlike energy, the ideal-solution chemical potentials diverge at infinite dilution. We use the analytical derivative instead of estimating it from neighboring grid points; this avoids a numerical differentiation error in the tangent construction.

return value1, value2 returns a tuple. The assignment mu_a, mu_b = ... unpacks it. In the main notebook mu(...) plays this role. An f-string begins with f; {mu_a:.3f} inserts the value using three digits after the decimal point. Formatting does not change the numerical value held in memory.

Try: check that mu_b-mu_a equals energy_slope(...) at the same state. Neither potential alone is the slope.

11.5 Build the common-tangent solver, one operation at a time¶

def common_tangent(temp):
    def composition(m, phase):
        ga, gb = endmembers[phase]
        return expit((m-(gb-ga))/(R*temp))

    def intercept(m, phase):
        ga, gb = endmembers[phase]
        return ga-R*temp*np.logaddexp(0, (m-(gb-ga))/(R*temp))

    def difference(m):
        return intercept(m, "alpha")-intercept(m, "liquid")

    left, right = -200000.0, 200000.0
    if difference(left)*difference(right) >= 0:
        raise ValueError("This bracket does not enclose a sign change.")
    m = brentq(difference, left, right, xtol=1e-9)
    xa = composition(m, "alpha")
    xl = composition(m, "liquid")
    b = intercept(m, "alpha")
    if not 0 < xa < xl < 1:
        raise ValueError("Expected an A-rich solid and a B-rich liquid.")
    return xa, xl, m, b

xa, xl, m, b = common_tangent(T_demo)
print(f"x_alpha={xa:.6f}, x_liquid={xl:.6f}, slope={m:.3f}, intercept={b:.3f}")

1. common_tangent(temp) is the outer function. The nested functions can use temp and endmembers from the surrounding scope. This is why temp does not have to be passed separately into each nested function call.
2. composition(m, phase) solves $g'(x)=m$ analytically for this ideal model. expit(y)=1/(1+\exp(-y)) evaluates the logistic function stably; manually evaluating a very large exponential can overflow.
3. intercept(...) evaluates $g(x)-mx$ at the composition selected by that slope. The formula $g_A^0-RT\ln[1+\exp((m-(g_B^0-g_A^0))/(RT))]$ follows from substituting the logistic composition. np.logaddexp(0,y) computes $\ln(1+e^y)$ stably. It is not log(0+y).
4. difference(m) returns one scalar residual. At its zero, both tangent intercepts agree; equal slopes are already enforced by selecting both compositions using the same m.
5. The bracket endpoints are trial slopes, in J/mol because mole fraction is dimensionless. They are not compositions. For this example their residuals have opposite signs.
6. In brentq(difference, ...), pass the function itself: no parentheses after difference. The solver chooses trial slopes and calls that function repeatedly. Writing difference(m) here would pass a number instead of a callable.
7. brentq combines bracketing with interpolation to find the root. xtol controls an absolute tolerance on the root variable; it does not automatically guarantee every physical residual to that same numerical tolerance. We separately check chemical-potential equality.
8. The solver returns the slope m, not the phase compositions; the next two lines convert it into xa and xl. The tangent is then $b+mx$.
9. The ordering check detects a result incompatible with this particular phase ordering. Different phase models or reversed ordering require a different implementation.
10. Equal tangents are not enough in a more general model: check that the line supports every allowed branch. The main notebook includes that stability check and an independent minimization.

The main notebook's coexistence(T) implements this same algorithm and guards its supported temperature range. A failed bracket must be investigated; do not merely widen it indefinitely or silently accept a failed solver.

11.6 Lever-rule arrays, clip, where, and the underscore¶

def phase_amounts(zz, temp):
    zz = np.asarray(zz, dtype=float)
    if np.any((zz < 0) | (zz > 1)):
        raise ValueError("Overall compositions must lie in [0,1].")
    xa, xl, slope, intercept = common_tangent(temp)
    fl = np.clip((zz-xa)/(xl-xa), 0.0, 1.0)
    fa = 1-fl
    geq = np.where(zz <= xa, energy(zz,temp,"alpha"),
           np.where(zz >= xl, energy(zz,temp,"liquid"), intercept+slope*zz))
    return fa, fl, geq

fa, fl, geq = phase_amounts(z_demo, T_demo)
print(f"Solid={fa:.6f}, liquid={fl:.6f}, g_eq={geq:.3f} J/mol")
fa_array, fl_array, _ = phase_amounts(grid, T_demo)

np.clip(values, 0, 1) limits each value to the interval [0,1]. Here it implements the correct single-phase continuation of the fraction, after the two coexistence compositions are known. It must not be used to conceal a physically invalid equilibrium solution.

np.where(condition, a, b) selects a where the condition is true and b elsewhere. The nested expression selects the α curve on the left, liquid curve on the right, and common tangent only inside the tie line. Python evaluates its input expressions before np.where selects among them: it does not protect an invalid logarithm in a discarded branch. This is one reason our energy function handles the endpoints explicitly.

_ is an ordinary variable name conventionally used for a returned value that will not be used. Here the third returned array is computed, then bound to _; it is not a command to skip that computation.

Check: inside the tie line, verify fa*xa + fl*xl == z_demo within numerical tolerance. Outside it, the occupied single phase has composition z_demo, so do not use the tie-line endpoint composition in that balance.

11.7 Matplotlib from the beginning: Figure, Axes and Axis¶

fig, ax = plt.subplots(figsize=(8, 5), layout="constrained")
ax.plot(grid, energy(grid,T_demo,"alpha")/1000, color="#1665a7", label="alpha")
ax.plot(grid, energy(grid,T_demo,"liquid")/1000, color="#bb4b25", label="liquid")
ax.plot(grid, (b+m*grid)/1000, color="black", linestyle="--", label="Common tangent")
ax.scatter([xa,xl], [(b+m*xa)/1000,(b+m*xl)/1000], color="black", zorder=5)
ax.set(xlim=(0,1), xlabel="B mole fraction", ylabel="Molar Gibbs energy / kJ mol$^{-1}$")
ax.set_title("Two-phase equilibrium at 800 K")
ax.legend(loc="upper center")
fig.savefig("my-common-tangent.png", dpi=180, bbox_inches="tight")
plt.show()

The names are easy to confuse: Axes is one plotting panel, even though the word looks plural. An Axis is the x- or y-axis within that panel and manages ticks, tick labels and scale. ax.xaxis and ax.yaxis access those objects. Use fig for operations on the whole figure and ax for operations on a panel.

11.8 Both vertical axes, intercept labels and coordinate systems¶

fig, ax = plt.subplots(figsize=(8,5), layout="constrained")
ax.plot(grid, (b+m*grid)/1000, "k--", label="Common tangent")
ax.set_xlim(0,1)
ax.set_ylim(-3,0)
ax.spines["right"].set_visible(True)
right = ax.secondary_yaxis("right", functions=(lambda y: y, lambda y: y))
ax.set_ylabel("Energy / kJ mol$^{-1}$")
right.set_ylabel("Same energy scale / kJ mol$^{-1}$")
ax.scatter([0,1], [b/1000,(b+m)/1000], color=["#1665a7","#bb4b25"], clip_on=False)
ax.annotate(r"$\mu_A$", xy=(0,b/1000), xytext=(10,12),
            textcoords="offset points", ha="left")
ax.annotate(r"$\mu_B$", xy=(1,(b+m)/1000), xytext=(-10,12),
            textcoords="offset points", ha="right")
ax.text(0.5,0.95,"The same scale on both sides", transform=ax.transAxes, ha="center")
plt.show()

• ax.set_xlim(0,1) places the plot boundaries at the pure-component compositions. Without it, automatic x-margins can put the left/right spines outside $x_B=0,1$, making the intercept construction visually misleading.
• Spines are the border lines around a panel. ax.spines["right"] retrieves the right border. The main notebook also explicitly positions the left spine at data coordinate zero.
• secondary_yaxis("right", functions=(forward,inverse)) builds a linked right-hand scale. Both functions here are lambda y: y: the identity transformation. It therefore uses the same physical scale as the left axis. A lambda is a small anonymous function; lambda y: y means “take y and return y unchanged.”
• We deliberately avoid an independently autoscaled ax.twinx() for these intercepts. A second, unrelated y-scale could visually imply false chemical-potential values. The intercepts belong to the same energy coordinate system.
• b/1000 is the left intercept and (b+m)/1000 the right intercept. clip_on=False allows endpoint markers centered on the border to remain fully visible.
• For annotate, xy is the point being labeled, in data coordinates by default. xytext gives the label displacement. With textcoords="offset points", the displacement is in typographic points (1/72 inch), not composition or energy units. ha controls horizontal alignment.
• ax.text(..., transform=ax.transAxes) uses panel-relative coordinates: (0,0) is the lower-left corner and (1,1) the upper-right, regardless of the data limits. This is useful for labels that should remain in a corner as the numerical limits change.
• r"$\mu_A$" is a raw Python string containing Matplotlib mathematical text. The r preserves backslashes; the dollar signs request mathematical rendering. Plain Python strings, mathematical notation and Markdown are different layers.

The notebook's intercept_axes(...) packages this repeated formatting into one helper function. Its unit=1000.0 default keeps intercept labels consistent with kJ/mol plots. pad=32 on a title adds space in points, while fontsize, bbox, alpha and zorder respectively control text size, its background box, opacity and drawing order.

11.9 Multipanel figures: why use zip?¶

temperatures = [600.0, 800.0, 1200.0]
fig, axes = plt.subplots(1,3, figsize=(14,4), layout="constrained")
for ax, temp in zip(axes, temperatures):
    for phase in endmembers:
        ax.plot(grid, energy(grid,temp,phase)/1000, label=phase)
    ax.set_title(f"T = {temp:g} K")
    ax.set(xlabel="B mole fraction", ylabel="Energy / kJ mol$^{-1}$")
    ax.legend()
plt.show()

plt.subplots(1,3) makes one row of three panels. Here axes is a one-dimensional array of three Axes handles, not one panel. By contrast, plt.subplots() returns one Axes; calling ax.plot(...) on the whole array of panels will fail. For a 2×2 layout, use axes[row,column] or flatten it with axes.ravel() to iterate through all panels.

zip(axes, temperatures) walks through the two collections together. The first loop iteration receives (axes[0], 600.0), the second (axes[1], 800.0), and the third (axes[2], 1200.0). The loop header unpacks each pair into ax and temp. This is clearer than manually keeping an index synchronized with two lists.

print(list(zip(["left","middle","right"], [600,800,1200])))
print(list(zip(["left","middle","right"], [600,800])))  # Stops at the shorter input!

zip returns an iterator: items are produced as you request them. list(...) collects them so you can inspect them. An exhausted iterator does not automatically restart; call zip(...) again to make another one. Normally zip silently stops at the shortest input. In Python 3.10+, zip(axes, temperatures, strict=True) can catch an unintended length mismatch.

The SER plot uses zip(axes, labels, tables) with three collections. Each iteration receives one panel, its title and its endmember table. This keeps the title and data paired. It does not combine the energy models mathematically.

The inner loop for phase in endmembers visits the dictionary keys in insertion order. Thus each panel gets a line for each phase. Indentation determines which operations repeat for every phase and which run once per panel. The :g in {temp:g} uses general number formatting, so 800.0 is displayed as 800.

11.10 What itertools.combinations does in the three-phase calculation¶

phase_names = ["alpha", "liquid", "beta"]
pairs = list(combinations(phase_names,2))
print(pairs)

The output is [("alpha","liquid"), ("alpha","beta"), ("liquid","beta")]. Each unordered pair appears once. There is no (alpha,alpha) and no duplicate reversed pair (liquid,alpha). For three objects there are $\binom{3}{2}=3$ pairs. For four phases there would be six.

This is different from zip: zip pairs corresponding positions in separate sequences; combinations selects all distinct groups of a specified size from one sequence. itertools.product would produce Cartesian products, while permutations distinguishes order. Those would introduce unnecessary repeated pairs here.

centers_demo = {"alpha":.15, "liquid":.50, "beta":.85}
K_demo = 20000.0
TE_demo = 1000.0

def offset_demo(phase,temp):
    return -20.0*(temp-TE_demo) if phase == "liquid" else 0.0

for pi,pj in combinations(centers_demo,2):
    ci,cj = centers_demo[pi],centers_demo[pj]
    hi,hj = offset_demo(pi,TE_demo),offset_demo(pj,TE_demo)
    slope_pair = (hj-hi)/(cj-ci)
    xi,xj = ci+slope_pair/(2*K_demo),cj+slope_pair/(2*K_demo)
    intercept_pair = K_demo*(xi-ci)**2+hi-slope_pair*xi
    print(pi,pj,xi,xj,slope_pair,intercept_pair)

Line by line: construct the centers dictionary; set curvature and invariant temperature; define the liquid's temperature-dependent vertical offset; generate a pair; retrieve its centers and heights; calculate the shared tangent slope; recover the two contact compositions; calculate its intercept; print the candidate. The expression a if condition else b is a conditional expression returning one of two values. It is appropriate for the scalar phase-name comparison here.

At $T_E$ all three pairs yield slope zero and intercept zero. They are three descriptions of one tangent, not three distinct equilibrium lines.

11.11 Why every candidate tangent needs a stability check¶

For a trial line $b+mx$ and branch $k$, $$g_k(x)-(b+mx)=K(x-c_k)^2+h_k-b-mx.$$ Its unconstrained minimum is at $x=c_k+m/(2K)$. On [0,1], clip that position to the interval. If the minimum is negative, the line passes above the branch somewhere, so it cannot be the lower supporting tangent.

trial_temp = 1050.0
trial_m, trial_b = 0.0, 0.0                # The old alpha–beta horizontal tangent.
gaps = []                                 # Empty list to collect minimum gaps.
for phase,center in centers_demo.items():  # .items() yields (key,value) pairs.
    contact = np.clip(center+trial_m/(2*K_demo),0,1)
    value = K_demo*(contact-center)**2+offset_demo(phase,trial_temp)
    gaps.append(value-trial_m*contact-trial_b)
print(gaps)
print("Supporting line:", min(gaps) >= -1e-8)

At 1050 K the liquid minimum is −1000 J/mol, so this candidate fails. append adds one item to a list; min(gaps) selects the smallest scalar gap. The small negative tolerance permits floating-point roundoff, not a physically meaningful stability violation.

The main notebook writes the same loop as a list comprehension, [expression for pk,ck in centers.items()]. Read it as “evaluate this expression for each phase/center pair and collect the values.” Expand it into the loop above if you find that easier to debug.

Inside stable_three_tangents, continue skips a candidate whose contact compositions lie outside [0,1]. A retained record (pi,pj,xi,xj,m,b) contains phase names, contact compositions, slope and intercept. Later for pi,pj,xi,xj,m,b in ... unpacks each record. q[4:] selects its slope and intercept using a slice from index 4 to the end.

For plotting, seen records already drawn (m,b) pairs. any(np.allclose([m,b],old) for old in seen) asks whether any previous line is numerically the same. The parenthesized expression is a generator, evaluated lazily. any stops as soon as it finds a match. This prevents drawing the identical invariant tangent three times; it does not eliminate the three physical phases.

11.12 Reading and extending the SER implementation¶

The SER block has a nested dictionary: unary["A"] is itself a dictionary, and unary["A"]["Cp"] selects its heat capacity. Each property should be accompanied by an explicit unit. A short local name such as u = unary[element] avoids repeating that lookup.

reference_temperature = 298.15
unary_demo = {"A":{"Cp":25.,"S_ref":30.,"Tm":1400.,"H_fus":56000/3},
              "B":{"Cp":30.,"S_ref":40.,"Tm":600.,"H_fus":36000.}}
solid_values, liquid_values = [], []
for element in ["A","B"]:
    u = unary_demo[element]
    gs = u["Cp"]*(T_demo-reference_temperature) - T_demo*(
        u["S_ref"]+u["Cp"]*np.log(T_demo/reference_temperature))
    fusion = u["H_fus"]*(1-T_demo/u["Tm"])
    solid_values.append(gs)
    liquid_values.append(gs+fusion)
ser_table = {"alpha":tuple(solid_values), "liquid":tuple(liquid_values)}
solid_reference = np.array(ser_table["alpha"])
reduced_table = {phase:tuple(np.array(values)-solid_reference)
                 for phase,values in ser_table.items()}
print(ser_table)
print(reduced_table)

• The explicit for element in ["A","B"] fixes the ordering of the resulting tuples. It must match the energy function's (A,B) convention.
• The multiline expression is continued inside parentheses; no backslash continuation is needed. np.log is the natural logarithm, and its temperature ratio is dimensionless.
• Append the solid Gibbs value to one list and the liquid Gibbs value to the other. The latter is the same-reference solid value plus the fusion difference.
• Convert the lists into tuples when forming the phase table. This communicates that they are fixed ordered endmember pairs for this calculation.
• reduced_table = {... for ...} is a dictionary comprehension. For each phase it subtracts the same two species reference values, converts the array back to a tuple and stores the result under that phase's key.

The main comparison temporarily assigns different tables to refs. saved_refs = refs.copy() makes a shallow copy; that is sufficient here because the values are tuples of numbers. try: ... finally: refs = saved_refs restores the original table even if plotting or calculation raises an error. Without the restoration, later cells could silently use the wrong reference convention. In a larger project, explicitly passing the endmember table into every function is preferable to relying on a mutable global variable. The small notebook uses a shared table to keep the first lesson compact.

11.13 Numerical self-checks: what an assertion does¶

np.testing.assert_allclose(
    potentials(xa,T_demo,"alpha"), potentials(xl,T_demo,"liquid"), atol=1e-6)
assert 0 <= fa <= 1 and 0 <= fl <= 1
np.testing.assert_allclose(fa+fl,1.0,atol=1e-12)
np.testing.assert_allclose(fa*xa+fl*xl,z_demo,atol=1e-12)
print("Chemical-potential equality and material balance passed.")

assert condition stops with an AssertionError when its condition is false. A passed assertion is normally silent; it is not a displayed proof. np.testing.assert_allclose compares scalars or arrays with numerical tolerances and reports a useful difference on failure. Its acceptance uses both absolute and relative tolerances; set rtol=0 if you intend a purely absolute bound. Tolerances must be meaningful in the units of the quantity being checked.

Do not usually test floating-point results with exact ==: calculations involving logarithms and root solvers produce rounding error. Conversely, making a tolerance extremely loose merely to obtain a passing result defeats the check.

Useful independent checks in this notebook include the analytical coexistence formula, conservation, global supporting-line stability, invariance under a shared reference shift and pure-element equality at the melting temperature. They test physical consequences, not just whether the code ran without an exception.

11.14 The independent linprog check, explained line by line¶

check_grid = np.linspace(0,1,2001)
costs = np.r_[energy(check_grid,T_demo,"alpha"), energy(check_grid,T_demo,"liquid")]
state_compositions = np.tile(check_grid,2)
balances = np.vstack([np.ones_like(state_compositions), state_compositions])
result = linprog(costs, A_eq=balances, b_eq=[1,z_demo], bounds=(0,None), method="highs")
assert result.success, result.message
print("Discrete minimum:",result.fun,"J/mol")
print("Difference from tangent solution:",result.fun-float(geq),"J/mol")

Because the grid allows only finitely many compositions, its minimum should be slightly higher than the continuous tangent minimum (within solver roundoff). Increase resolution to investigate convergence. A small grid error is not evidence that the thermodynamic equality conditions are wrong.

11.15 Assemble, troubleshoot and extend your own notebook¶

Build the notebook in this order:

1. Purpose and assumptions: define components, phases, pressure, composition basis, energy reference and units in Markdown.
2. Imports and parameters: put editable controls in one obvious cell. Separate physical parameters from plot settings.
3. Pure functions: implement energy, derivative and chemical potentials. Check one scalar value and array shapes before plotting.
4. Equilibrium calculation: implement and inspect the residual, bracket the root and compute the phase compositions.
5. Physical checks: verify both chemical potentials, supporting-line stability and material balance.
6. Plots: create axes handles, label units, mark tangency points and display the two intercepts on consistent scales.
7. Parameter exploration: use loops and zip for comparisons; state which model parameters actually vary with temperature.
8. Extensions: introduce a reference conversion or extra phase only after the two-phase calculation is understood.
9. Reproducibility: restart and run all cells, inspect the saved figures, then save the notebook. A notebook with stale outputs is not a reproducible calculation.

Small coding challenges: replace zip with an index-based loop and explain the difference; use zip(..., strict=True) to detect a mismatch; expand a list comprehension into an ordinary loop; draw a 2×2 panel layout using axes.ravel(); change only the figure size, then only the numerical grid resolution, and explain why those changes have different meanings; pass the endmember table as an explicit function argument instead of a global variable. Keep the existing physical self-checks passing after each change.