Chapter 4: Hands-On: Screening in Python

Five Fictitious Catalysts, Four Steps Each, and a Ranking That the Volcano Reproduces Exactly

📖 Reading Time: 20-25 minutes 📊 Difficulty: Beginner 💻 Code Examples: 0 📝 Exercises: 0

Video Lecture

This video covers the same content as the text below. Choose your preferred learning format.

🌐 EN | 🇯🇵 JP | Last sync: 2026-08-18

Materials Informatics Dojo > Computational Chemistry of OER > Chapter 4

Three chapters of machinery deserve a machine. This one builds the smallest complete screening pipeline that still contains every logical move a real one makes: define candidates, compute their intermediates, construct free-energy diagrams, extract a performance number, rank, and then check the ranking against the cheap descriptor model from Chapter 3.

It runs in NumPy, in under a second, on five catalysts that do not exist.

That last point is the design of the chapter, not a limitation of it. The adsorption free energies below are illustrative teaching values on fictitious surfaces called Catalyst A through Catalyst E. They were chosen so that the five land at five different, instructive positions on the volcano — one binding oxygen too strongly, one too weakly, one close to the apex, two in between. They are not DFT results, they are not measurements, and they must not be attributed to any real material. What is real is everything downstream: every step energy, every diagram, every overpotential and every ranking in this chapter is arithmetic performed by the code on those stated inputs plus the two constants of Chapter 3, and every printed number is the code's actual output.

A screening pipeline whose inputs are honest teaching values and whose logic is exact is a better teacher than one with real numbers and hidden steps. Swap in DFT energies and nothing about the code changes.

These are not Chapter 2's catalysts. Chapter 2 used a separate set of fictitious surfaces, also lettered, defined directly by their four step energies. The five here are defined by their adsorption energies and are an unrelated cast. Names reset at the chapter boundary; nothing carries over but the constants.

4.1 The Plan

Five stages, each consuming what the last produced.

Stage Input Output
1. Define the dataset five stated \((\Delta G_{\ast\text{OH}}, \Delta G_{\ast\text{O}})\) pairs descriptors \(x = \Delta G_{\ast\text{O}} - \Delta G_{\ast\text{OH}}\)
2. Build the steps the dataset + scaling relation \(\Delta G_{\ast\text{OOH}}\) and all four \(\Delta G_i\), sum-checked
3. Free-energy diagrams the step energies cumulative profiles at \(U = 0\) and \(U = 1.23\) V
4. Overpotentials the step energies \(U_L\), \(\eta\), the potential-determining step, a ranking
5. Volcano placement the descriptors alone predicted \(\eta\), checked against stage 4

The point of stage 5 is not to compute anything new. It is to confirm that a model using one number per catalyst reproduces, exactly, a calculation that used four. If the two disagree the descriptor is not doing its job; if they agree, Chapter 3's collapse was real.

4.2 Stage 1: The Teaching Dataset

A screening study begins with a table of candidates and the quantities you were able to compute for each. Ours has five rows and two computed columns.

The two constants come with different licences and the code says so in its comments. \(1.23\) V is a definition — the standard equilibrium potential of the four-electron reaction — and \(4 \times 1.23 = 4.92\) eV is arithmetic on it. The scaling offset of \(3.2\) eV is an empirical, approximate regularity, introduced and hedged in Chapter 3, and it is the only quantity here that could be wrong in an interesting way.

import numpy as np

# ------------------------------------------------------------------
# CONSTANTS. Two are definitions/arithmetic, one is empirical.
# ------------------------------------------------------------------
E_EQ = 1.23                    # V, equilibrium potential of the OER (definition)
N_STEPS = 4
SUM_TOTAL = N_STEPS * E_EQ     # eV, arithmetic: 4 x 1.23
SCALING = 3.2                  # eV, dG_OOH - dG_OH; empirical and approximate

# ------------------------------------------------------------------
# THE TEACHING DATASET.
#
# Catalysts A-E are FICTITIOUS. The adsorption free energies below are
# ILLUSTRATIVE values invented for this chapter to place five surfaces
# at five different positions on the volcano. They are not measured,
# not computed from DFT, and must not be attributed to any real
# material. Everything downstream of them is arithmetic on these
# numbers plus the two constants above.
# ------------------------------------------------------------------
CATALYSTS = [
    # name, dG_OH (eV), dG_O (eV), intended role
    ("Catalyst A", 0.60, 1.60, "binds oxygen too strongly"),
    ("Catalyst B", 0.80, 2.20, "mid, strong-binding side"),
    ("Catalyst C", 1.00, 2.65, "near the apex"),
    ("Catalyst D", 1.10, 3.00, "mid, weak-binding side"),
    ("Catalyst E", 1.30, 3.60, "binds oxygen too weakly"),
]

names = [c[0] for c in CATALYSTS]
dG_OH = np.array([c[1] for c in CATALYSTS])
dG_O = np.array([c[2] for c in CATALYSTS])
roles = [c[3] for c in CATALYSTS]
descriptor = dG_O - dG_OH

print("Stage 1: the teaching dataset (ILLUSTRATIVE values, fictitious surfaces)")
print(f"{'catalyst':<12} {'dG_OH':>7} {'dG_O':>7} {'x = dG_O - dG_OH':>18}   role")
for k, name in enumerate(names):
    print(f"{name:<12} {dG_OH[k]:7.2f} {dG_O[k]:7.2f} {descriptor[k]:18.2f}   {roles[k]}")
print()
print(f"descriptor range : {descriptor.min():.2f} to {descriptor.max():.2f} eV")
print(f"all distinct     : {len(set(np.round(descriptor, 6))) == len(names)}")
print(f"dG_O > dG_OH for every surface (O is the later intermediate): "
      f"{bool(np.all(dG_O > dG_OH))}")

Output:

Stage 1: the teaching dataset (ILLUSTRATIVE values, fictitious surfaces)
catalyst       dG_OH    dG_O   x = dG_O - dG_OH   role
Catalyst A      0.60    1.60               1.00   binds oxygen too strongly
Catalyst B      0.80    2.20               1.40   mid, strong-binding side
Catalyst C      1.00    2.65               1.65   near the apex
Catalyst D      1.10    3.00               1.90   mid, weak-binding side
Catalyst E      1.30    3.60               2.30   binds oxygen too weakly

descriptor range : 1.00 to 2.30 eV
all distinct     : True
dG_O > dG_OH for every surface (O is the later intermediate): True

Reading the result. The descriptors span \(1.00\) to \(2.30\) eV, straddling the apex position of \(1.60\) eV that Chapter 3 derived — A and B to the strong-binding side, C just past the apex, D and E to the weak-binding side. That spread is a deliberate act of dataset design, and it is worth being conscious of it: a screening set that happens to sample only one leg of a volcano will produce a monotonic-looking relationship between descriptor and activity, and someone reading the plot will conclude that stronger binding is always better. It is not; it is better only until the apex.

The two sanity lines at the bottom cost nothing and catch typos. Distinct descriptors mean no ties will complicate the ranking later, and \(\Delta G_{\ast\text{O}} > \Delta G_{\ast\text{OH}}\) confirms the intermediates are ordered the way the mechanism requires — \(\ast\text{O}\) sits further along the reaction coordinate than \(\ast\text{OH}\), so a negative difference would signal an input error rather than an exotic surface.

4.3 Stage 2: From Two Numbers to Four Steps

Now the scaling relation earns its keep. We know \(\Delta G_{\ast\text{OH}}\) and \(\Delta G_{\ast\text{O}}\); we do not know \(\Delta G_{\ast\text{OOH}}\). Rather than compute the third intermediate, we predict it — and that prediction is the whole reason a one-descriptor volcano exists.

The four steps follow as differences of cumulative energies:

\[ \Delta G_1 = \Delta G_{\ast\text{OH}}, \quad \Delta G_2 = \Delta G_{\ast\text{O}} - \Delta G_{\ast\text{OH}}, \quad \Delta G_3 = \Delta G_{\ast\text{OOH}} - \Delta G_{\ast\text{O}}, \quad \Delta G_4 = 4.92 - \Delta G_{\ast\text{OOH}} \]

The last one is fixed by the sum rule, which makes the obvious check — do the four add to \(4.92\)? — circular by construction. So the code does a second, non-circular check: it rebuilds \(\Delta G_4\) along an independent algebraic route and compares.

# ------------------------------------------------------------------
# Stage 2: the scaling relation supplies dG_OOH, and the four steps
# follow as differences. The final step is fixed by the sum rule,
# which is exactly why the sum check below can only ever be trivially
# satisfied -- so we ALSO rebuild step 4 independently and compare.
# ------------------------------------------------------------------
dG_OOH = dG_OH + SCALING

dG1 = dG_OH - 0.0                 # H2O + *   -> *OH  + H+ + e-
dG2 = dG_O - dG_OH                # *OH       -> *O   + H+ + e-
dG3 = dG_OOH - dG_O               # *O + H2O  -> *OOH + H+ + e-
dG4 = SUM_TOTAL - dG_OOH          # *OOH      -> * + O2 + H+ + e-

steps = np.vstack([dG1, dG2, dG3, dG4]).T          # shape (5 catalysts, 4 steps)
STEP_LABELS = ["1  H2O -> *OH", "2  *OH -> *O",
               "3  *O -> *OOH", "4  *OOH -> O2"]

print("Stage 2: intermediates and the four step energies (eV)")
print(f"{'catalyst':<12} {'dG_OH':>7} {'dG_O':>7} {'dG_OOH':>8} "
      f"{'dG1':>7} {'dG2':>7} {'dG3':>7} {'dG4':>7} {'sum':>8}")
for k, name in enumerate(names):
    print(f"{name:<12} {dG_OH[k]:7.2f} {dG_O[k]:7.2f} {dG_OOH[k]:8.2f} "
          f"{dG1[k]:7.2f} {dG2[k]:7.2f} {dG3[k]:7.2f} {dG4[k]:7.2f} "
          f"{steps[k].sum():8.4f}")
print()

# CROSS-CHECK 1: the sum constraint, to full floating-point tolerance.
sums = steps.sum(axis=1)
print(f"sum constraint: every catalyst sums to {SUM_TOTAL:.2f} eV : "
      f"{bool(np.allclose(sums, SUM_TOTAL, atol=1e-12))}")
print(f"  largest deviation from {SUM_TOTAL:.2f} eV : "
      f"{np.max(np.abs(sums - SUM_TOTAL)):.2e} eV")

# CROSS-CHECK 2: rebuild dG4 from the intermediates alone, without
# using the sum rule, and confirm the two routes agree.
dG4_independent = SUM_TOTAL - SCALING - dG_OH
print(f"  dG4 rebuilt as SUM_TOTAL - SCALING - dG_OH agrees : "
      f"{bool(np.allclose(dG4, dG4_independent, atol=1e-12))}")

# CROSS-CHECK 3: the scaling relation collapses the middle pair.
print(f"  dG2 + dG3 equals the scaling constant {SCALING:.2f} eV for all five : "
      f"{bool(np.allclose(dG2 + dG3, SCALING, atol=1e-12))}")
print(f"  dG1 + dG4 equals {SUM_TOTAL - SCALING:.2f} eV for all five        : "
      f"{bool(np.allclose(dG1 + dG4, SUM_TOTAL - SCALING, atol=1e-12))}")
print()
print("Every step is positive (no spontaneous step at U = 0) : "
      f"{bool(np.all(steps > 0.0))}")

Output:

Stage 2: intermediates and the four step energies (eV)
catalyst       dG_OH    dG_O   dG_OOH     dG1     dG2     dG3     dG4      sum
Catalyst A      0.60    1.60     3.80    0.60    1.00    2.20    1.12   4.9200
Catalyst B      0.80    2.20     4.00    0.80    1.40    1.80    0.92   4.9200
Catalyst C      1.00    2.65     4.20    1.00    1.65    1.55    0.72   4.9200
Catalyst D      1.10    3.00     4.30    1.10    1.90    1.30    0.62   4.9200
Catalyst E      1.30    3.60     4.50    1.30    2.30    0.90    0.42   4.9200

sum constraint: every catalyst sums to 4.92 eV : True
  largest deviation from 4.92 eV : 0.00e+00 eV
  dG4 rebuilt as SUM_TOTAL - SCALING - dG_OH agrees : True
  dG2 + dG3 equals the scaling constant 3.20 eV for all five : True
  dG1 + dG4 equals 1.72 eV for all five        : True

Every step is positive (no spontaneous step at U = 0) : True

Reading the result. Look down the \(\Delta G_2\) and \(\Delta G_3\) columns together and the weld is visible to the naked eye. As \(\Delta G_2\) climbs from \(1.00\) to \(2.30\) eV across the five surfaces, \(\Delta G_3\) falls from \(2.20\) to \(0.90\) — every increase in one is an exactly equal decrease in the other, because their sum is nailed to \(3.20\) eV. Five different materials, and this particular pair of numbers has only one degree of freedom among them.

The same thing happens in the outer pair, and it is easy to miss. \(\Delta G_1 + \Delta G_4 = 4.92 - 3.20 = 1.72\) eV for every catalyst, again because the scaling relation ties \(\Delta G_{\ast\text{OOH}}\) to \(\Delta G_{\ast\text{OH}}\). Four steps, two welded pairs, two free parameters. That is the entire design space these surfaces live in.

The three cross-checks deserve a word, because their asymmetry is the lesson. Check 1 — the sum — is guaranteed by how \(\Delta G_4\) was constructed, so it can only ever print True and it verifies nothing about the physics. Its value is narrow but real: it catches a broken array shape or a vstack transposed the wrong way. Check 2 is the meaningful one, arriving at \(\Delta G_4\) by a different algebraic path. In a pipeline this small the distinction is pedantic; in a pipeline with a thousand candidates and six modules it is the difference between a test suite and a decoration.

4.4 Stage 3: Free-Energy Diagrams

The free-energy diagram is the standard picture of an electrocatalytic mechanism, and it is more informative than any single number extracted from it.

Plot the cumulative free energy against reaction progress: start at zero for the clean surface, add \(\Delta G_1\) to reach \(\ast\text{OH}\), add \(\Delta G_2\) to reach \(\ast\text{O}\), and so on. At an applied potential \(U\), every electron already transferred lowers the current level by \(eU\), so the level after \(i\) steps is

\[ G_i(U) = \sum_{j \le i} \Delta G_j - i\,eU \]

The applied potential tilts the whole staircase downward, uniformly, one tread at a time. Catalysis begins when no tread points upward — that is the definition of the limiting potential we extract in stage 4.

We draw each catalyst twice: at \(U = 0\), where the diagram is pure chemistry, and at \(U = 1.23\) V, where the thermodynamic minimum has been supplied and whatever is still uphill is the catalyst's own fault.

# ------------------------------------------------------------------
# Stage 3: free-energy diagrams.
#
# At applied potential U every electron transferred lowers the state
# by eU, so the level after i steps is
#     G_i(U) = sum(dG_1..dG_i) - i * e * U
# The diagram is the list of levels; the catalyst is "all downhill"
# when every consecutive difference is <= 0.
# ------------------------------------------------------------------
def profile(step_row, U):
    """Cumulative free-energy levels at potential U, starting at 0."""
    return np.concatenate([[0.0], np.cumsum(step_row) - U * np.arange(1, 5)])

LEVEL_LABELS = ["*", "*OH", "*O", "*OOH", "* + O2"]

for U in (0.0, E_EQ):
    print(f"Stage 3: free-energy profile at U = {U:.2f} V  (levels in eV)")
    print(f"{'catalyst':<12} " + " ".join(f"{lab:>9}" for lab in LEVEL_LABELS)
          + "   largest uphill step")
    for k, name in enumerate(names):
        lv = profile(steps[k], U)
        rises = np.diff(lv)
        j = int(np.argmax(rises))
        print(f"{name:<12} " + " ".join(f"{v:9.2f}" for v in lv)
              + f"   {STEP_LABELS[j]}  (+{rises[j]:.2f} eV)")
    print()

# The end point at U = E_EQ is a check on the sum rule, seen from the
# other side: four electrons at 1.23 V exactly pay off 4.92 eV.
final_at_eq = np.array([profile(steps[k], E_EQ)[-1] for k in range(len(names))])
print(f"at U = {E_EQ:.2f} V every profile ends at 0.00 eV : "
      f"{bool(np.allclose(final_at_eq, 0.0, atol=1e-12))}")
print(f"  largest final level : {np.max(np.abs(final_at_eq)):.2e} eV")

# Which step is still uphill at U = E_EQ? That step IS the bottleneck.
print()
print(f"{'catalyst':<12} " + " ".join(f"{lab:>16}" for lab in STEP_LABELS))
for k, name in enumerate(names):
    rises = np.diff(profile(steps[k], E_EQ))
    print(f"{name:<12} " + " ".join(f"{r:16.2f}" for r in rises))
print()
n_uphill = [int(np.sum(np.diff(profile(steps[k], E_EQ)) > 1e-12))
            for k in range(len(names))]
print(f"steps still uphill at U = {E_EQ:.2f} V, per catalyst : {n_uphill}")
print("every catalyst has at least one uphill step at the equilibrium "
      f"potential : {all(n > 0 for n in n_uphill)}")

Output:

Stage 3: free-energy profile at U = 0.00 V  (levels in eV)
catalyst             *       *OH        *O      *OOH    * + O2   largest uphill step
Catalyst A        0.00      0.60      1.60      3.80      4.92   3  *O -> *OOH  (+2.20 eV)
Catalyst B        0.00      0.80      2.20      4.00      4.92   3  *O -> *OOH  (+1.80 eV)
Catalyst C        0.00      1.00      2.65      4.20      4.92   2  *OH -> *O  (+1.65 eV)
Catalyst D        0.00      1.10      3.00      4.30      4.92   2  *OH -> *O  (+1.90 eV)
Catalyst E        0.00      1.30      3.60      4.50      4.92   2  *OH -> *O  (+2.30 eV)

Stage 3: free-energy profile at U = 1.23 V  (levels in eV)
catalyst             *       *OH        *O      *OOH    * + O2   largest uphill step
Catalyst A        0.00     -0.63     -0.86      0.11      0.00   3  *O -> *OOH  (+0.97 eV)
Catalyst B        0.00     -0.43     -0.26      0.31      0.00   3  *O -> *OOH  (+0.57 eV)
Catalyst C        0.00     -0.23      0.19      0.51      0.00   2  *OH -> *O  (+0.42 eV)
Catalyst D        0.00     -0.13      0.54      0.61      0.00   2  *OH -> *O  (+0.67 eV)
Catalyst E        0.00      0.07      1.14      0.81      0.00   2  *OH -> *O  (+1.07 eV)

at U = 1.23 V every profile ends at 0.00 eV : True
  largest final level : 0.00e+00 eV

catalyst        1  H2O -> *OH     2  *OH -> *O    3  *O -> *OOH    4  *OOH -> O2
Catalyst A              -0.63            -0.23             0.97            -0.11
Catalyst B              -0.43             0.17             0.57            -0.31
Catalyst C              -0.23             0.42             0.32            -0.51
Catalyst D              -0.13             0.67             0.07            -0.61
Catalyst E               0.07             1.07            -0.33            -0.81

steps still uphill at U = 1.23 V, per catalyst : [1, 2, 2, 2, 2]
every catalyst has at least one uphill step at the equilibrium potential : True

Reading the result. Four things, and the third is the one to carry forward.

4.5 Stage 4: Overpotentials and the Ranking

Now compress each diagram into one number.

The limiting potential \(U_L\) is the smallest potential at which no step is uphill. Since every step is tilted by the same \(eU\), that is simply the largest step energy expressed as a potential:

\[ U_L = \max_i \Delta G_i / e, \qquad \eta = U_L - 1.23\ \text{V} \]

The step attaining the maximum is the potential-determining step. Reporting it alongside \(\eta\) costs one line and turns a score into a diagnosis: it tells you which direction to push the material.

# ------------------------------------------------------------------
# Stage 4: the theoretical overpotential.
#
#     U_L  = max_i dG_i / e        (the limiting potential: the
#                                   smallest U at which no step is
#                                   uphill)
#     eta  = U_L - E_EQ
# ------------------------------------------------------------------
U_L = steps.max(axis=1)
eta = U_L - E_EQ
pds = [STEP_LABELS[int(np.argmax(steps[k]))] for k in range(len(names))]

print("Stage 4: limiting potential and overpotential")
print(f"{'catalyst':<12} {'U_L (V)':>9} {'eta (V)':>9}   potential-determining step")
for k, name in enumerate(names):
    print(f"{name:<12} {U_L[k]:9.3f} {eta[k]:9.3f}   {pds[k]}")
print()

order = list(np.argsort(eta, kind="stable"))
print("Ranking, best (lowest overpotential) first")
for rank, k in enumerate(order, start=1):
    print(f"  {rank}. {names[k]:<12} eta = {eta[k]:.3f} V   "
          f"limited by {pds[k]:<16} x = {descriptor[k]:.2f} eV   ({roles[k]})")
print()

eta_floor = SCALING / 2.0 - E_EQ
print(f"floor implied by the scaling relation : {eta_floor:.3f} V")
print(f"best catalyst in this set             : {names[order[0]]}, "
      f"eta = {eta[order[0]]:.3f} V")
print(f"gap between the two                   : "
      f"{eta[order[0]] - eta_floor:.3f} V")
print(f"no catalyst beats the floor           : "
      f"{bool(np.all(eta >= eta_floor - 1e-12))}")
print(f"spread across the set                 : "
      f"{eta.max() - eta.min():.3f} V")

Output:

Stage 4: limiting potential and overpotential
catalyst       U_L (V)   eta (V)   potential-determining step
Catalyst A       2.200     0.970   3  *O -> *OOH
Catalyst B       1.800     0.570   3  *O -> *OOH
Catalyst C       1.650     0.420   2  *OH -> *O
Catalyst D       1.900     0.670   2  *OH -> *O
Catalyst E       2.300     1.070   2  *OH -> *O

Ranking, best (lowest overpotential) first
  1. Catalyst C   eta = 0.420 V   limited by 2  *OH -> *O     x = 1.65 eV   (near the apex)
  2. Catalyst B   eta = 0.570 V   limited by 3  *O -> *OOH    x = 1.40 eV   (mid, strong-binding side)
  3. Catalyst D   eta = 0.670 V   limited by 2  *OH -> *O     x = 1.90 eV   (mid, weak-binding side)
  4. Catalyst A   eta = 0.970 V   limited by 3  *O -> *OOH    x = 1.00 eV   (binds oxygen too strongly)
  5. Catalyst E   eta = 1.070 V   limited by 2  *OH -> *O     x = 2.30 eV   (binds oxygen too weakly)

floor implied by the scaling relation : 0.370 V
best catalyst in this set             : Catalyst C, eta = 0.420 V
gap between the two                   : 0.050 V
no catalyst beats the floor           : True
spread across the set                 : 0.650 V

Reading the result. The near-apex catalyst wins, and the reason it wins is printed next to it.

Catalyst C's limiting step is \(1.650\) eV — \(\Delta G_2\), by a margin of only \(0.10\) eV over its \(\Delta G_3\) of \(1.55\). The two welded steps are nearly equal, which is precisely the condition Chapter 3 identified as optimal: with their sum pinned at \(3.20\) eV, the maximum is smallest when they are balanced. C is not better because it binds oxygen strongly, or weakly; it is better because it binds oxygen evenly, splitting the welded pair almost down the middle. Its \(\eta = 0.420\) V sits just \(0.050\) V above the floor of \(0.370\) V, and that \(0.050\) V is exactly its imbalance — half of the \(0.10\) eV difference between its two middle steps.

The ranking is also a warning about intuition. Catalyst B (\(x = 1.40\)) and Catalyst D (\(x = 1.90\)) sit on opposite legs, and B wins — not for any chemical reason, but because \(|1.40 - 1.60| = 0.20\) is smaller than \(|1.90 - 1.60| = 0.30\). Distance from the apex is the only thing that matters, and it does not care which side you approach from. Meanwhile Catalyst A, the strong binder, and Catalyst E, the weak binder, finish last with \(0.970\) and \(1.070\) V — a \(0.650\) V spread from best to worst across a set whose adsorption energies differ by well under an electronvolt. Small differences in binding, large differences in performance. That sensitivity is what makes screening worth doing, and it is also what makes screening errors expensive.

4.6 Stage 5: Onto the Volcano

Stage 4 used all four step energies for each catalyst. Chapter 3 claimed that one number would do. Stage 5 tests that claim, with no new physics: it maps each descriptor onto the volcano curve, reads off a predicted overpotential, and compares against stage 4 element by element.

If this fails, either the volcano is wrong or the pipeline is. If it succeeds, the four-dimensional problem really was one-dimensional all along.

# ------------------------------------------------------------------
# Stage 5: place the five catalysts on the Chapter 3 volcano.
#
# The volcano knows ONE number per catalyst, the descriptor
# x = dG_O - dG_OH, and predicts
#     eta_volcano(x) = max(x, SCALING - x) - E_EQ
# Stage 4 used all four step energies. If the two agree, the
# descriptor really has absorbed everything that mattered.
# ------------------------------------------------------------------
X_APEX = SCALING / 2.0

def eta_volcano(x):
    return np.maximum(x, SCALING - x) - E_EQ

eta_pred = eta_volcano(descriptor)
leg = ["strong-binding (left)" if x < X_APEX else
       ("weak-binding (right)" if x > X_APEX else "apex") for x in descriptor]

print("Stage 5: volcano placement")
print(f"apex at x = {X_APEX:.2f} eV, floor eta = {eta_floor:.3f} V")
print(f"{'catalyst':<12} {'x':>6} {'|x - apex|':>11} {'eta_volcano':>12} "
      f"{'eta_stage4':>11} {'activity -eta':>14}   leg")
for k, name in enumerate(names):
    print(f"{name:<12} {descriptor[k]:6.2f} {abs(descriptor[k] - X_APEX):11.2f} "
          f"{eta_pred[k]:12.3f} {eta[k]:11.3f} {-eta[k]:14.3f}   {leg[k]}")
print()

# CROSS-CHECK 1: the descriptor prediction reproduces stage 4 exactly.
print(f"eta_volcano == eta_stage4 for all five : "
      f"{bool(np.allclose(eta_pred, eta, atol=1e-12))}")
print(f"  largest discrepancy : {np.max(np.abs(eta_pred - eta)):.2e} V")

# CROSS-CHECK 2: the rankings are identical, element by element.
order_volcano = list(np.argsort(eta_pred, kind="stable"))
print(f"  stage 4 ranking : {[names[k] for k in order]}")
print(f"  stage 5 ranking : {[names[k] for k in order_volcano]}")
print(f"  identical       : {order == order_volcano}")

# CROSS-CHECK 3: distance from the apex alone predicts the loss.
print(f"  eta == floor + |x - apex| for all five : "
      f"{bool(np.allclose(eta, eta_floor + np.abs(descriptor - X_APEX), atol=1e-12))}")
print()

# Why the agreement is not automatic: steps 1 and 4 had to stay out of
# the way. Show the margin by which they did.
margin = U_L - np.maximum(dG1, dG4)
print("Margin by which steps 1 and 4 stayed below the limiting step (eV)")
for k, name in enumerate(names):
    print(f"  {name:<12} max(dG1, dG4) = {max(dG1[k], dG4[k]):.2f}   "
          f"dG_max = {U_L[k]:.2f}   margin = {margin[k]:.2f}")
print(f"  all margins positive : {bool(np.all(margin > 0.0))}")
print("  (had any margin gone negative, the one-descriptor volcano would "
      "have mis-ranked that catalyst)")

Output:

Stage 5: volcano placement
apex at x = 1.60 eV, floor eta = 0.370 V
catalyst          x  |x - apex|  eta_volcano  eta_stage4  activity -eta   leg
Catalyst A     1.00        0.60        0.970       0.970         -0.970   strong-binding (left)
Catalyst B     1.40        0.20        0.570       0.570         -0.570   strong-binding (left)
Catalyst C     1.65        0.05        0.420       0.420         -0.420   weak-binding (right)
Catalyst D     1.90        0.30        0.670       0.670         -0.670   weak-binding (right)
Catalyst E     2.30        0.70        1.070       1.070         -1.070   weak-binding (right)

eta_volcano == eta_stage4 for all five : True
  largest discrepancy : 2.22e-16 V
  stage 4 ranking : ['Catalyst C', 'Catalyst B', 'Catalyst D', 'Catalyst A', 'Catalyst E']
  stage 5 ranking : ['Catalyst C', 'Catalyst B', 'Catalyst D', 'Catalyst A', 'Catalyst E']
  identical       : True
  eta == floor + |x - apex| for all five : True

Margin by which steps 1 and 4 stayed below the limiting step (eV)
  Catalyst A   max(dG1, dG4) = 1.12   dG_max = 2.20   margin = 1.08
  Catalyst B   max(dG1, dG4) = 0.92   dG_max = 1.80   margin = 0.88
  Catalyst C   max(dG1, dG4) = 1.00   dG_max = 1.65   margin = 0.65
  Catalyst D   max(dG1, dG4) = 1.10   dG_max = 1.90   margin = 0.80
  Catalyst E   max(dG1, dG4) = 1.30   dG_max = 2.30   margin = 1.00
  all margins positive : True
  (had any margin gone negative, the one-descriptor volcano would have mis-ranked that catalyst)

Reading the result. Three points.

4.7 What a Real Screening Adds

This pipeline is complete in its logic and hollow in its inputs. Everything that would make it a research tool sits in the gap between "illustrative values" and "computed energies", so it is worth being specific about what fills that gap.

DFT energies instead of teaching values. Each \(\Delta G_{\ast\text{OH}}\) and \(\Delta G_{\ast\text{O}}\) becomes a slab calculation: build a surface, choose a facet and a termination, decide a coverage, relax the geometry, and add the zero-point and entropic corrections that Chapter 2 introduced to turn an electronic energy into a free energy. That is hours of compute and a dozen human decisions per candidate, and the decisions — which facet, which termination, which coverage — often move the result more than the functional does.

Error bars, and the honesty to propagate them. Section 4.5 showed a \(0.650\) V spread across catalysts whose binding energies differ by less than an electronvolt. Now note that DFT adsorption energies routinely carry uncertainties of a few tenths of an electronvolt, and that the scaling relation itself has scatter around its \(3.2\) eV offset. Put those together with the one-to-one exchange rate of Section 4.6 and the conclusion is uncomfortable: the uncertainty on \(\eta\) is comparable to the differences you are ranking by. Catalyst C beat Catalyst B by \(0.15\) V here. In a real study that margin might not survive contact with the error bars, and a screening that reports a ranking without them is reporting a ranking it cannot defend. The correct output of a real pipeline is not "C is best" but "C, B and D are indistinguishable; A and E are clearly worse".

Stability filters, which are usually the real bottleneck. Nothing in this chapter asked whether any of these surfaces survives the potentials and acidity of actual operation. Many of the best computational OER candidates dissolve, oxidize past the intended state, or reconstruct into something else entirely under operating conditions. Activity screening produces a list; stability screening — Pourbaix analysis, dissolution potentials, surface reconstruction — is what shortens it, and it eliminates far more candidates than activity does. A pipeline that ranks only by \(\eta\) is answering half the question.

Kinetics, which the whole framework has been quietly ignoring. The theoretical overpotential is thermodynamic: it asks whether each step is downhill, never how fast it goes. Real activation barriers, coverage effects, the electrolyte, and mass transport are all outside the model. Two catalysts with identical \(\eta\) can differ by orders of magnitude in measured current.

And then there is the assumption underneath all of it. Every number in this chapter is downstream of the scaling relation, which we hedged carefully as an approximate empirical regularity and then used as though it were exact — because that is what a screening pipeline does. The five catalysts obeyed it perfectly, by construction. Real surfaces scatter around it, and the interesting ones are the outliers.

That is where Chapter 5 begins. It asks what happens when you hand this problem to a machine learning model: what a learned descriptor buys over a derived one, why a model trained on scaling-obedient surfaces will confidently reproduce a volcano it has no way to question, and how the limitations of the computational hydrogen electrode — the ones we accepted in Chapter 2 and have been compounding ever since — propagate into every prediction a model makes from data built this way. It is the chapter that audits the other four.

🎯 Exercise Problems

  1. A sixth candidate. Add Catalyst F with \(\Delta G_{\ast\text{OH}} = 0.95\) eV and \(\Delta G_{\ast\text{O}} = 2.55\) eV. Predict its overpotential and its rank before running the code, then check. How close does it come to the floor, and what would \(\Delta G_{\ast\text{O}}\) have to be for it to sit exactly at the apex?
  2. Break the descriptor. Construct a candidate with \(x = 1.60\) eV — apex-perfect — but \(\Delta G_{\ast\text{OH}} = 1.85\) eV. Run stages 2 through 5 on it. Which cross-check in stage 5 fails first, and what does the margin calculation report?
  3. Sensitivity to the empirical constant. Re-run the whole pipeline with the scaling offset set to \(3.0\) and \(3.4\) eV. Which of the five catalysts changes rank, and why does the identity of the potential-determining step change for some but not others?
  4. Error bars, crudely. Add independent Gaussian noise of standard deviation \(0.15\) eV to each \(\Delta G_{\ast\text{OH}}\) and \(\Delta G_{\ast\text{O}}\), repeat the ranking 10,000 times, and report how often each catalyst comes first. Restate the chapter's conclusion as a probability rather than an ordering.
  5. The check that was not run. Stage 3 reported that Catalyst E has an uphill first step at \(U = 1.23\) V. Write a diagnostic that flags any candidate whose \(\Delta G_1\) or \(\Delta G_4\) comes within \(0.20\) eV of its limiting step, and explain why such a candidate deserves a full four-step treatment even when a descriptor model is available.

Summary

This chapter built a five-stage OER screening pipeline in NumPy and ran it end to end on five fictitious catalysts with illustrative adsorption free energies, chosen to occupy five positions on the Chapter 3 volcano. Stage 1 defined the dataset and its descriptors \(x = \Delta G_{\ast\text{O}} - \Delta G_{\ast\text{OH}}\), spanning \(1.00\) to \(2.30\) eV around the apex at \(1.60\) eV. Stage 2 applied the scaling relation to predict \(\Delta G_{\ast\text{OOH}}\) and built all four step energies, verifying the \(4.92\) eV sum to zero floating-point error, rebuilding \(\Delta G_4\) by an independent route, and exhibiting both welds directly: \(\Delta G_2 + \Delta G_3 = 3.20\) eV and \(\Delta G_1 + \Delta G_4 = 1.72\) eV for every surface. Stage 3 drew cumulative free-energy profiles at \(U = 0\) and \(U = 1.23\) V — all ending at exactly \(0.00\) eV at the equilibrium potential, a check on the sum rule from the electrochemical side — and showed the two volcano legs as two distinct bottlenecks: step 3 for the strong binders A and B, step 2 for C, D and E. Stage 4 extracted limiting potentials and overpotentials, ranking C (0.420 V) < B (0.570) < D (0.670) < A (0.970) < E (1.070); the near-apex catalyst C won because its two welded steps were nearly equal (\(1.65\) against \(1.55\) eV), placing it only \(0.050\) V above the \(0.370\) V floor. Stage 5 mapped the five onto the volcano using the descriptor alone and reproduced stage 4's overpotentials to \(2 \times 10^{-16}\) V and its ranking exactly, confirmed \(\eta = \eta_{\text{floor}} + |x - x_{\text{apex}}|\) for all five, and printed the per-catalyst margin (\(0.65\) to \(1.08\) eV) by which steps 1 and 4 stayed clear of the limiting step — the condition that made the one-descriptor shortcut valid rather than lucky. Finally we named what a real screening adds: DFT energies with their facet, termination and coverage decisions; error bars comparable in size to the differences being ranked; stability filters that eliminate more candidates than activity does; and kinetics, which this framework never modelled at all.

Chapter 5 turns the audit on the whole construction. It asks what machine learning adds to a problem this constrained, and what it silently inherits — from the computational hydrogen electrode, from the scaling relation, and from the fact that a model trained on well-behaved surfaces has no way to recognize the outliers that would actually matter.

← Chapter 3: Scaling Relations and the Volcano Chapter 5: ML Screening and the Limits of CHE →

Disclaimer

Back to Series Index Next Chapter →