Text · original UTF-8 line locators · 12.6 KiB
#!/usr/bin/env python3
"""Derive the canonical OCx25 modelling targets from the bound reference inputs.
Standard library only. Run it on the exact bytes TfP serves for a Problem:
python3 prepare_targets.py --experiments experiments.csv --xrf xrf.csv --out-dir targets/
It writes:
targets/targets_her.csv one row per (XRF composition, source): HER full-cell
voltage at 50 mA/cm2 and its half-cell conversion
targets/targets_co2rr.csv one row per (XRF composition, source): CO2RR
Faradaic efficiencies and production rates at a
full-cell voltage of 3.3 V (log-linear interpolation)
The steps are a port of the processing in the OCx24 reference code
(fairchem, src/fairchem/applications/ocx/core/data_handling.py and voltage.py,
MIT licence), with one deliberate simplification so that the derivation needs
no 1.5 GB computational file: replicates are pooled by synthesis target
(the sample id without its replicate suffix, per source), and the pooled
row carries the element-wise mean of its replicates' XRF compositions,
rather than pooling by the nearest composition in the computational
screening set. Everything else (dropping the avg/std rows,
XRF lookup by source, abs(voltage), mean over replicates, the electron
counts, the log-linear interpolation to 3.3 V without extrapolation, the
half-cell conversion, the element-system label used for leave-one-
composition-out folds) follows the reference code. See reference.md.
"""
from __future__ import annotations
import argparse
import csv
import math
import os
import re
from collections import OrderedDict, defaultdict
FARADAY = 96485.0 # (s A)/mol
PRODUCTS = ["fe_h2", "fe_co", "fe_ch4", "fe_c2h4", "fe_liquid"]
# Electron counts exactly as the reference implementation uses them
# (get_all_production_rates). Note that it uses 2 for C2H4; the physical
# value is 12. Keep the reference convention so numbers stay comparable.
ELECTRONS = {"fe_h2": 2, "fe_co": 2, "fe_ch4": 8, "fe_liquid": 2, "fe_c2h4": 2}
RATE_NAME = {
"fe_h2": "h2_pr",
"fe_co": "co_pr",
"fe_ch4": "ch4_pr",
"fe_liquid": "liquid_pr",
"fe_c2h4": "c2h4_pr",
}
INTERPOLATE_V = 3.3
def read_csv(path: str) -> list[dict[str, str]]:
with open(path, newline="", encoding="utf-8") as f:
return list(csv.DictReader(f))
def strip_rep(sample_id: str) -> str:
return re.sub(r"_rep\d+", "", sample_id)
def xrf_mapping(xrf_rows: list[dict[str, str]]) -> dict[str, str]:
"""vsp samples use the XRF measured on the electrode (gde); uoft samples use
the wafer measurement keyed by the sample id without its replicate suffix."""
mapping: dict[str, str] = {}
for row in xrf_rows:
source = row["source"]
sid = row["sample id"]
if source == "vsp":
value = row.get("xrf-gde_formula", "")
if value:
mapping[sid] = value
elif source == "uoft":
value = row.get("xrf-wafer_formula", "")
if value:
mapping[strip_rep(sid)] = value
return mapping
def lookup_xrf(row: dict[str, str], mapping: dict[str, str]) -> str | None:
if row["source"] == "vsp":
return mapping.get(row["sample id"])
if row["source"] == "uoft":
return mapping.get(strip_rep(row["sample id"]))
return None
def elements_of(comp: str) -> str:
"""'Ag-0.5-Au-0.5' -> 'AgAu' (the leave-one-composition-out group label)."""
return "".join(comp.split("-")[::2])
def to_float(value: str) -> float:
value = (value or "").strip()
if value == "":
return math.nan
return float(value)
def mean(values: list[float]) -> float:
values = [v for v in values if not math.isnan(v)]
if not values:
return math.nan
return sum(values) / len(values)
# --- half-cell conversion (port of fairchem voltage.py) --------------------
def _she2rhe(ushe: float, ph: float, ref_pot: float) -> float:
return ushe + ref_pot + 0.059 * ph
def _rhe2she(urhe: float, ph: float, ref_pot: float) -> float:
return urhe - ref_pot - 0.059 * ph
def _correct_potential(pot: float, r: float, ph: float, j: float, area: float, ref_pot: float) -> float:
if pot < 0:
return _she2rhe(pot + j / 1000 * area * r, ph, ref_pot)
return _she2rhe(pot - j / 1000 * area * r, ph, ref_pot)
def _linear_fit(xs: list[float], ys: list[float]) -> tuple[float, float]:
n = len(xs)
mx = sum(xs) / n
my = sum(ys) / n
sxx = sum((x - mx) ** 2 for x in xs)
sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
a = sxy / sxx
return a, my - a * mx
def get_she(vcell: float) -> float:
"""Full-cell voltage -> cathode potential vs SHE, using the calibration the
reference code embeds (neutral CO2RR, 4 cm2 cell, sputtered Cu, 0.1 M
bicarbonate, Ag/AgCl 3 M KCl at 0.23 V vs SHE). Affine in vcell."""
ref_pot = 0.23
cathode_ph = 12.5
area = 4.0
j = [50.0, 100.0, 200.0]
cathode_pot = [-1.62, -2.0, -2.3]
cathode_r = [0.48, 0.34, 0.3]
fullcell_pot = [3.0, 3.4, 3.7]
cathode_pot_corr = [
_correct_potential(cathode_pot[i], cathode_r[i], cathode_ph, j[i], area, ref_pot) for i in range(3)
]
a, b = _linear_fit(fullcell_pot, cathode_pot_corr)
urhe = a * vcell + b
return _rhe2she(urhe, cathode_ph, ref_pot)
# --- processing ------------------------------------------------------------
def base_rows(exp_rows: list[dict[str, str]], mapping: dict[str, str], reaction: str) -> list[dict]:
out = []
for row in exp_rows:
if row["reaction"] != reaction:
continue
if row["rep"] in ("avg", "std"):
continue
xrf = lookup_xrf(row, mapping)
if not xrf:
continue
out.append(
{
"sample id": row["sample id"],
"source": row["source"],
"composition": row["composition"],
"xrf comp": xrf,
"current density": float(row["current density"]),
"voltage": abs(to_float(row["voltage"])),
**{p: to_float(row.get(p, "")) for p in PRODUCTS},
}
)
return out
def parse_comp(comp: str) -> dict[str, float]:
parts = comp.split("-")
return {parts[i]: float(parts[i + 1]) for i in range(0, len(parts) - 1, 2)}
def format_comp(fractions: dict[str, float]) -> str:
return "-".join(f"{el}-{fractions[el]:.3f}" for el in sorted(fractions))
def mean_comp(comps: list[str]) -> str:
"""Element-wise mean of replicate XRF compositions, renormalised to 1."""
totals: dict[str, float] = defaultdict(float)
for c in comps:
for el, x in parse_comp(c).items():
totals[el] += x
n = len(comps)
fractions = {el: v / n for el, v in totals.items()}
s = sum(fractions.values()) or 1.0
return format_comp({el: v / s for el, v in fractions.items()})