Custom Extractants: Defining and Using Novel REE Extractants#
This notebook demonstrates how to create and use custom extractants with the difflow_ree module.
Why Custom Extractants?#
While difflow_ree ships 5 extractant systems (D2EHPA, PC88A, Cyanex272, TBP and naphthenic acid), you may need to:
Model novel extractants from research literature
Calibrate extractant mixtures or modified systems
Test hypothetical extractants for process optimization
Fit empirical data from lab experiments
What You’ll Learn#
Create custom extractants with pH-dependent properties
Register custom extractants with the database
Use custom extractants in simulations
Calibrate extractant parameters from data — and recognise a fit that has failed
Compare custom vs. standard extractants
import jax
import jax.numpy as jnp
from jax import grad, jit
jax.config.update("jax_enable_x64", True)
from difflow_ree import (
# Custom extractant creation
create_custom_extractant,
get_extractant_database,
# Standard functions
get_extractant,
list_extractants,
list_ree_elements,
REEDistribution,
REEExtractor,
REEExtractorParams,
)
from difflow.streams import make_stream, get_flows
1. Understanding the Distribution Model#
Custom extractants use the same pH-dependent distribution model:
\[\log_{10}(D) = a + b \cdot pH + c \cdot pH^2 + d\left(\frac{1}{T} - \frac{1}{T_{ref}}\right)\]
Required parameters for each REE element:
a,b,c: pH dependence coefficientsd: Temperature coefficient (K, optional)
Choosing the coefficients#
b is not a free parameter. Cation exchange of a trivalent ion on a dimeric
acidic extractant,
\[\mathrm{RE}^{3+} + 3\,\overline{(\mathrm{HA})_2} \rightleftharpoons \overline{\mathrm{RE}(\mathrm{HA}_2)_3} + 3\,\mathrm{H}^+,\]
releases three protons, so mass action fixes the slope at \(b = 3\) exactly —
a factor of 1000 in \(D\) per pH unit. Every element on a given extractant shares
it. All four acidic records shipped with difflow_ree are pinned there since the
#270 refit; three of them used to carry b between 2.35 and 2.55, hand-tuned,
which is where the artefacts below came from. If you fit b freely and get 2.4,
that is a statement about physics your model is missing — loading,
aggregation, activity coefficients, a third phase — not about the
stoichiometry.
Selectivity lives in a, not in b. With one shared slope,
\[\log_{10}\beta_{ij} = \log_{10}\frac{D_i}{D_j} = a_i - a_j,\]
so the separation factor is independent of pH, temperature and extractant
concentration. That is a real and useful property: it means the operating pH
sets recovery, and the reagent choice sets selectivity, and the two decisions
come apart. Staggering b between elements instead makes \(\beta\) drift with pH,
which is the artefact #265 and #270 removed from this package — and with it
the “optimal pH for separation” that an optimiser would otherwise climb straight
out of the fitted window.
a: the baseline. \(a_i = \log_{10} K_{e,i}\) referred toreference_concentration. \(D_i = 1\) at \(\mathrm{pH} = -a_i / b\). Adjacent lanthanides typically differ by \(\Delta a \approx 0.2\)–\(0.5\) (\(\beta \approx 1.5\)–\(3\)), rising across the series.b: 3 for a trivalent cation exchanger. Leave it there.c: 0 unless you have data that demands curvature. Section 7 shows what happens when you fit it from seven points: it comes back with the wrong sign and a standard error six times its own size, and it inflates the error bars onaandbfour- to nine-fold while it does so.d: enthalpy of extraction, negative for exothermic (\(-1000\) to \(-2500\) K).
2. Creating a Simple Custom Extractant#
Let’s create a hypothetical extractant with moderate selectivity.
# Create custom extractant
#
# b = 3.0 for every element -- the stoichiometric slope for a trivalent cation
# exchange (section 1). Selectivity is set entirely by the spread in `a`:
#
# beta(Nd/La) = 10**(a_Nd - a_La) = 10**1.18 = 15.1
# beta(Dy/Nd) = 10**(a_Dy - a_Nd) = 10**1.78 = 60.3
#
# and, because the slope is shared, both are the same at every pH.
# a_Nd = -6.00 puts D(Nd) = 1 at pH 2.00, which is where this reagent would be
# run; the validity window is set around that.
my_extractant = create_custom_extractant(
name="MyExtractant",
full_name="My Novel Phosphoric Acid Extractant",
formula="C10H20O4P",
molecular_weight=250.0,
# pH coefficients for each element
# Format: {element: {"a": value, "b": value, "c": value}}
ph_coefficients={
"La": {"a": -7.18, "b": 3.0, "c": 0.0},
"Nd": {"a": -6.00, "b": 3.0, "c": 0.0},
"Dy": {"a": -4.22, "b": 3.0, "c": 0.0},
},
# Temperature corrections (K)
temperature_coefficients={
"La": -1550.0,
"Nd": -1750.0,
"Dy": -2200.0,
},
# Physical properties
density=0.98, # g/mL
pKa=3.4, # Acid dissociation constant
extractant_type="acidic_phosphoric",
# Operational parameters
typical_concentration=0.5, # M
stoichiometry_protons=3, # H+ released per extraction
stoichiometry_extractant=3, # Extractant molecules per complex
valid_ph_range=(0.5, 3.0), # the range these coefficients are meant for
valid_temp_range=(283.0, 333.0), # K
reference_concentration=0.5, # M (the basis of `a`)
concentration_exponent=3.0, # D ∝ [HA]^3
cost_usd_kg=15.0,
)
print("✓ Custom extractant created successfully!")
print(f" Name: {my_extractant.name}")
print(f" Formula: {my_extractant.formula}")
print(f" MW: {my_extractant.molecular_weight} g/mol")
print(f" Elements: {list(my_extractant.ph_coefficients.keys())}")
print(f" Valid pH range: {my_extractant.valid_ph_range}")
✓ Custom extractant created successfully!
Name: MyExtractant
Formula: C10H20O4P
MW: 250.0 g/mol
Elements: ['La', 'Nd', 'Dy']
Valid pH range: (0.5, 3.0)
3. Registering and Using Custom Extractants#
Register the extractant with the database to use it in simulations.
# Get the global extractant database
db = get_extractant_database()
# Check what's already registered
print("Before registration:")
print(f" Available: {list_extractants()}")
# Register custom extractant
db.add_extractant("MyExtractant", my_extractant)
print("\nAfter registration:")
print(f" Available: {list_extractants()}")
# Verify it can be retrieved
retrieved = get_extractant("MyExtractant")
print(f"\n✓ Custom extractant registered and retrievable")
print(f" Full name: {retrieved.full_name}")
Before registration:
Available: ['D2EHPA', 'PC88A', 'Cyanex272', 'TBP', 'naphthenic_acid']
After registration:
Available: ['D2EHPA', 'PC88A', 'Cyanex272', 'TBP', 'naphthenic_acid', 'MyExtractant']
✓ Custom extractant registered and retrievable
Full name: My Novel Phosphoric Acid Extractant
4. Using Custom Extractants in Distribution Models#
# Create distribution model with custom extractant
dist_custom = REEDistribution(
extractant="MyExtractant",
elements=("La", "Nd", "Dy"),
concentration=0.5,
)
# pH 2.0 -- the middle of this reagent's window, and where D(Nd) = 1.
pH_op = 2.0
D_values = dist_custom.get_D_all(pH=pH_op, T=298.15)
print(f"Custom Extractant Performance at pH {pH_op}:")
print("="*50)
for elem, D in D_values.items():
print(f" D({elem}) = {float(D):10.4f}")
# Separation factors. With one shared slope these do not depend on pH; check it.
print(f"\nSeparation Factors:")
print(f"{'pH':<8} {'SF(Nd/La)':<14} {'SF(Dy/Nd)':<14}")
print("-"*36)
for pH in [0.5, 1.0, 2.0, 3.0]:
Dv = dist_custom.get_D_all(pH=pH, T=298.15)
print(f"{pH:<8.1f} {float(Dv['Nd']/Dv['La']):<14.4f} {float(Dv['Dy']/Dv['Nd']):<14.4f}")
print("\n\u2191 identical rows: b is shared, so log10(beta) = a_i - a_j and the")
print(" separation factor is a property of the reagent, not of the operating point.")
Custom Extractant Performance at pH 2.0:
==================================================
D(La) = 0.0661
D(Nd) = 1.0000
D(Dy) = 60.2560
Separation Factors:
pH SF(Nd/La) SF(Dy/Nd)
------------------------------------
0.5 15.1356 60.2560
1.0 15.1356 60.2560
2.0 15.1356 60.2560
3.0 15.1356 60.2560
↑ identical rows: b is shared, so log10(beta) = a_i - a_j and the
separation factor is a property of the reagent, not of the operating point.
5. Comparing Custom vs. Standard Extractants#
# Compare with D2EHPA. Pick a pH inside BOTH validity windows: MyExtractant is
# [0.5, 3.0] and D2EHPA, since the #270 refit, is [0.0, 2.0], so pH 2.0 is the
# one point legal for both -- and legal only because the bound test is
# inclusive. (No two records share a window any more, and comparing two
# correlations at a pH outside one of them compares a fit with an
# extrapolation.)
dist_d2ehpa = REEDistribution(
extractant="D2EHPA",
elements=("La", "Nd", "Dy"),
concentration=0.5,
)
D_d2ehpa = dist_d2ehpa.get_D_all(pH=pH_op, T=298.15)
print(f"Comparison: Custom vs. D2EHPA at pH {pH_op}")
print("="*70)
print(f"{'Element':<10} {'Custom D':<15} {'D2EHPA D':<15} {'Ratio':<12}")
print("-"*70)
for elem in ["La", "Nd", "Dy"]:
D_cust = float(D_values[elem])
D_std = float(D_d2ehpa[elem])
print(f"{elem:<10} {D_cust:<15.4f} {D_std:<15.4f} {D_cust/D_std:<12.2f}")
SF_custom = float(D_values["Nd"] / D_values["La"])
SF_d2ehpa = float(D_d2ehpa["Nd"] / D_d2ehpa["La"])
print(f"\nSeparation Factor (Nd/La) at pH {pH_op}:")
print(f" Custom: {SF_custom:.2f}")
print(f" D2EHPA: {SF_d2ehpa:.2f} (pH-independent too, since #270)")
print(f" Change: {((SF_custom/SF_d2ehpa - 1)*100):+.1f}%")
print("""
The ratio column is small because a_Nd was chosen to put D(Nd) = 1 at pH 2.0,
which is the TOP of D2EHPA's refitted window -- five decades into its range,
where it extracts everything. That is a statement about where each correlation
was anchored, not about which reagent is stronger.
The separation factor is the comparison worth making, and neither of these two
depends on pH: both records pin b = 3 across their elements, so each beta is a
single number per pair. Compare those, at whatever pH each reagent is actually
operated at.""")
Comparison: Custom vs. D2EHPA at pH 2.0
======================================================================
Element Custom D D2EHPA D Ratio
----------------------------------------------------------------------
La 0.0661 10120.4541 0.00
Nd 1.0000 115611.2242 0.00
Dy 60.2560 13128043.2857 0.00
Separation Factor (Nd/La) at pH 2.0:
Custom: 15.14
D2EHPA: 11.42 (pH-independent too, since #270)
Change: +32.5%
The ratio column is small because a_Nd was chosen to put D(Nd) = 1 at pH 2.0,
which is the TOP of D2EHPA's refitted window -- five decades into its range,
where it extracts everything. That is a statement about where each correlation
was anchored, not about which reagent is stronger.
The separation factor is the comparison worth making, and neither of these two
depends on pH: both records pin b = 3 across their elements, so each beta is a
single number per pair. Compare those, at whatever pH each reagent is actually
operated at.
6. Using Custom Extractants in Multi-Stage Units#
# Create extraction unit with custom extractant
# Note: include_loading=False -- a custom extractant has no loading isotherm
# data, and #268 derives the Langmuir constant from a declared reference state
# this record does not have.
params_custom = REEExtractorParams(
n_stages=5,
extractant="MyExtractant",
elements=("La", "Nd", "Dy"),
pH=pH_op,
include_loading=False,
)
extractor_custom = REEExtractor(params_custom)
# Define streams
feed = make_stream(
flows={"H2O": 10.0, "La": 0.01, "Nd": 0.02, "Dy": 0.01},
T=298.15, P=101325.0,
)
# Organic solvent stream.
#
# The stream must name the *extractant* and the *diluent* as species. A
# solvent whose carrier matches neither now raises instead of silently
# defaulting the organic flow to 1.0 (#192), and when loading is enabled the
# extractant molar flow is what sets the capacity of the organic phase
# (capacity = F_extractant / m, #191).
#
# 0.5 M MyExtractant in kerosene is roughly 10 mol% extractant (kerosene is
# ~0.75 g/mL and ~170 g/mol, so ~4.4 mol/L of diluent against 0.5 mol/L of
# extractant). The total organic flow is unchanged at 8.0 mol/s; it is just
# split 0.8 MyExtractant / 7.2 kerosene.
solvent = make_stream(
flows={"MyExtractant": 0.8, "kerosene": 7.2, "La": 0.0, "Nd": 0.0, "Dy": 0.0},
T=298.15, P=101325.0,
)
# Run extraction
raffinate, extract, info = extractor_custom(feed, solvent)
# Analyze results
feed_flows = get_flows(feed)
ext_flows = get_flows(extract)
print(f"Custom Extractant Performance (5 stages, pH {pH_op}):")
print("="*60)
print(f"{'Element':<10} {'Feed':<12} {'Extract':<12} {'Recovery %':<12}")
print("-"*60)
for elem in ["La", "Nd", "Dy"]:
feed_val = float(feed_flows[elem])
ext_val = float(ext_flows[elem])
recovery = (ext_val / feed_val) * 100
print(f"{elem:<10} {feed_val:<12.4f} {ext_val:<12.4f} {recovery:<12.1f}")
# Calculate Nd purity
total_REE = sum(float(ext_flows[e]) for e in ["La", "Nd", "Dy"])
nd_purity = float(ext_flows["Nd"]) / total_REE * 100
print(f"\nExtract Quality:")
print(f" Nd purity: {nd_purity:.1f}%")
print(" (Dy reports almost entirely to the extract at this pH: beta(Dy/Nd) = 60,")
print(" so Nd purity here is set by the feed's Dy content, not by the cascade.)")
Custom Extractant Performance (5 stages, pH 2.0):
============================================================
Element Feed Extract Recovery %
------------------------------------------------------------
La 0.0100 0.0005 5.3
Nd 0.0200 0.0145 72.7
Dy 0.0100 0.0100 100.0
Extract Quality:
Nd purity: 58.0%
(Dy reports almost entirely to the extract at this pH: beta(Dy/Nd) = 60,
so Nd purity here is set by the feed's Dy content, not by the cascade.)
7. Advanced: Calibrating Extractant Parameters from Data#
Fit pH coefficients to measured D values — and check the fit before believing it.
# ---------------------------------------------------------------------------
# The data.
#
# These points are SYNTHETIC, generated from a known correlation plus 4% (in
# log10 D) scatter, so that at the end we can ask whether the fit found the
# right answer. Earlier versions of this notebook fitted five made-up points
# whose slope was about 0.9 decades per pH unit; no trivalent cation exchange
# produces that curve, so there was nothing physical for a fit to find.
# ---------------------------------------------------------------------------
TRUE_a, TRUE_b, TRUE_c = -6.00, 3.00, 0.00
pH_data = jnp.array([1.00, 1.25, 1.50, 1.75, 2.00, 2.25, 2.50])
noise = 0.04 * jax.random.normal(jax.random.PRNGKey(0), pH_data.shape)
logD_data = TRUE_a + TRUE_b * pH_data + TRUE_c * pH_data**2 + noise
D_data = 10.0 ** logD_data
print("Synthetic 'measurements' (Nd):")
print(f"{'pH':<8} {'D':<14}")
print("-"*22)
for p, d in zip(pH_data, D_data):
print(f"{float(p):<8.2f} {float(d):<14.6g}")
# ---------------------------------------------------------------------------
# 1. How NOT to fit it: least squares on D itself, by gradient descent.
# ---------------------------------------------------------------------------
def loss_on_D(params):
a, b, c = params
D_pred = jnp.power(10.0, a + b * pH_data + c * pH_data**2)
return jnp.sum((D_pred - D_data)**2)
print("\n" + "="*70)
print("1. Least squares on D, gradient descent, learning rate 0.1")
print("="*70)
print(f"{'Iter':<6} {'a':<14} {'b':<14} {'c':<14} {'Loss':<12}")
print("-"*70)
p_bad = jnp.array([-8.0, 2.5, 0.01])
for i in range(5):
print(f"{i:<6} {float(p_bad[0]):<14.4g} {float(p_bad[1]):<14.4g} "
f"{float(p_bad[2]):<14.4g} {float(loss_on_D(p_bad)):<12.4g}")
p_bad = p_bad - 0.1 * grad(loss_on_D)(p_bad)
print(f"""
That is not a fit, it is a divergence, and it does not announce itself. One
step of 0.1 throws the parameters to order 1e28; every predicted D then
underflows to zero, the residual becomes -D_measured for every point, and the
loss parks on sum(D_measured**2) = {float(jnp.sum(D_data**2)):.0f} and STOPS
MOVING. A run that prints a frozen loss and a table of 100.00% errors has
failed. It looks converged.
Two things went wrong, and both are structural:
* SCALE. D spans four decades across these seven points, so squared error in
D is a fit to the single largest point and ignores the rest.
* CURVATURE. dD/da = D ln(10), so the gradient is proportional to D itself.
A step size that is reasonable at pH 1 overshoots by decades at pH 2.5.
The correlation is written in log10(D). Fit it there.""")
# ---------------------------------------------------------------------------
# 2. Least squares on log10(D) -- which is LINEAR in (a, b, c).
# ---------------------------------------------------------------------------
def fit_log10(columns):
"""Exact linear least squares of log10(D) on the given design columns."""
X = jnp.stack(columns, axis=1)
theta, *_ = jnp.linalg.lstsq(X, logD_data, rcond=None)
resid = logD_data - X @ theta
dof = len(pH_data) - X.shape[1]
s2 = float(resid @ resid) / dof
cov = s2 * jnp.linalg.inv(X.T @ X)
return theta, jnp.sqrt(jnp.diag(cov)), float(jnp.sqrt(jnp.mean(resid**2))), \
float(jnp.linalg.cond(X.T @ X))
ones = jnp.ones_like(pH_data)
print("\n" + "="*70)
print("2. Least squares on log10(D): a + b*pH + c*pH**2, all three free")
print("="*70)
theta3, se3, rms3, cond3 = fit_log10([ones, pH_data, pH_data**2])
for name, v, s, truth in zip("abc", theta3, se3, (TRUE_a, TRUE_b, TRUE_c)):
print(f" {name} = {float(v):>9.4f} +- {float(s):.4f} (true {truth})")
print(f" rms residual: {rms3:.4f} log units cond(X'X) = {cond3:,.0f}")
print("""
c comes back with the WRONG SIGN and a standard error six times its own size
-- which is to say, indistinguishable from zero -- and carrying it costs a
factor of four to nine on the error bars of a and b. pH and pH**2 are nearly
collinear over this window; the design is ill-conditioned, and seven points
cannot resolve a quadratic term. Drop it.""")
print("\n" + "="*70)
print("3. Same fit with c fixed at 0")
print("="*70)
theta2, se2, rms2, cond2 = fit_log10([ones, pH_data])
for name, v, s, truth in zip("ab", theta2, se2, (TRUE_a, TRUE_b)):
print(f" {name} = {float(v):>9.4f} +- {float(s):.4f} (true {truth})")
print(f" rms residual: {rms2:.4f} log units cond(X'X) = {cond2:,.0f}")
b_fit, b_se = float(theta2[1]), float(se2[1])
print(f"\n b = {b_fit:.3f} +- {b_se:.3f}: the stoichiometric slope of 3 sits "
f"{abs(b_fit - 3.0)/b_se:.1f} standard errors away.")
print(" The rms residual is unchanged, so the quadratic term was buying nothing.")
# ---------------------------------------------------------------------------
# 4. The same log-space objective by gradient descent, for comparison.
# The model is linear, so the Hessian is constant and the largest stable
# step is 1/lambda_max -- read it off rather than guessing 0.1.
# ---------------------------------------------------------------------------
X3 = jnp.stack([ones, pH_data, pH_data**2], axis=1)
lr = float(1.0 / jnp.linalg.eigvalsh(2.0 * X3.T @ X3).max())
def loss_on_logD(params):
a, b, c = params
return jnp.sum((a + b * pH_data + c * pH_data**2 - logD_data)**2)
@jit
def descend(p, n):
return jax.lax.fori_loop(0, n, lambda _, q: q - lr * grad(loss_on_logD)(q), p)
print("\n" + "="*70)
print(f"4. Gradient descent on the SAME objective, step 1/lambda_max = {lr:.2e}")
print("="*70)
print(f"{'Iter':<10} {'a':<12} {'b':<12} {'c':<12} {'Loss':<12}")
print("-"*70)
p_gd = jnp.array([-8.0, 2.5, 0.01])
seen = 0
for n in (100, 1_000, 10_000, 50_000):
p_gd = descend(p_gd, n - seen); seen = n
print(f"{n:<10} {float(p_gd[0]):<12.4f} {float(p_gd[1]):<12.4f} "
f"{float(p_gd[2]):<12.4f} {float(loss_on_logD(p_gd)):<12.6f}")
print(f"{'lstsq':<10} {float(theta3[0]):<12.4f} {float(theta3[1]):<12.4f} "
f"{float(theta3[2]):<12.4f} {float(loss_on_logD(theta3)):<12.6f}")
print(f"""
50,000 descent steps to approach what one call to lstsq returns exactly. The
condition number above, {cond3:,.0f}, is why: gradient descent converges at a
rate set by it. When a model is linear in its parameters, solve it; keep the
gradients for the parts that are not (the cascade, the flowsheet, the
economics), which is what the rest of difflow is for.""")
# ---------------------------------------------------------------------------
# 5. Fit quality of the answer we are keeping.
# ---------------------------------------------------------------------------
params = jnp.array([float(theta2[0]), float(theta2[1]), 0.0])
print("\n" + "="*70)
print("Fit quality (c = 0 fit):")
print("="*70)
print(f"{'pH':<8} {'D_measured':<15} {'D_predicted':<15} {'Error %':<12}")
print("-"*55)
for p, D_meas in zip(pH_data, D_data):
a, b, c = params
D_pred = float(jnp.power(10.0, a + b * p + c * p**2))
print(f"{float(p):<8.2f} {float(D_meas):<15.6g} {D_pred:<15.6g} "
f"{abs(D_pred - float(D_meas)) / float(D_meas) * 100:<12.2f}")
Synthetic 'measurements' (Nd):
pH D
----------------------
1.00 0.00098122
1.25 0.0052313
1.50 0.0373804
1.75 0.180931
2.00 1.00748
2.25 5.43395
2.50 35.2863
======================================================================
1. Least squares on D, gradient descent, learning rate 0.1
======================================================================
Iter a b c Loss
----------------------------------------------------------------------
0 -8 2.5 0.01 1274
1 -7.654 3.361 2.156 2.851e+28
2 -1.313e+28 -3.282e+28 -8.205e+28 1276
3 -1.313e+28 -3.282e+28 -8.205e+28 1276
4 -1.313e+28 -3.282e+28 -8.205e+28 1276
That is not a fit, it is a divergence, and it does not announce itself. One
step of 0.1 throws the parameters to order 1e28; every predicted D then
underflows to zero, the residual becomes -D_measured for every point, and the
loss parks on sum(D_measured**2) = 1276 and STOPS
MOVING. A run that prints a frozen loss and a table of 100.00% errors has
failed. It looks converged.
Two things went wrong, and both are structural:
* SCALE. D spans four decades across these seven points, so squared error in
D is a fit to the single largest point and ignores the rest.
* CURVATURE. dD/da = D ln(10), so the gradient is proportional to D itself.
A step size that is reasonable at pH 1 overshoots by decades at pH 2.5.
The correlation is written in log10(D). Fit it there.
======================================================================
2. Least squares on log10(D): a + b*pH + c*pH**2, all three free
======================================================================
a = -6.0544 +- 0.2192 (true -6.0)
b = 3.0593 +- 0.2647 (true 3.0)
c = -0.0116 +- 0.0750 (true 0.0)
rms residual: 0.0325 log units cond(X'X) = 8,418
c comes back with the WRONG SIGN and a standard error six times its own size
-- which is to say, indistinguishable from zero -- and carrying it costs a
factor of four to nine on the error bars of a and b. pH and pH**2 are nearly
collinear over this window; the design is ill-conditioned, and seven points
cannot resolve a quadratic term. Drop it.
======================================================================
3. Same fit with c fixed at 0
======================================================================
a = -6.0219 +- 0.0531 (true -6.0)
b = 3.0187 +- 0.0292 (true 3.0)
rms residual: 0.0326 log units cond(X'X) = 72
b = 3.019 +- 0.029: the stoichiometric slope of 3 sits 0.6 standard errors away.
The rms residual is unchanged, so the quadratic term was buying nothing.
======================================================================
4. Gradient descent on the SAME objective, step 1/lambda_max = 3.95e-03
======================================================================
Iter a b c Loss
----------------------------------------------------------------------
100 -6.7830 3.4791 -0.0362 0.242239
1000 -6.4877 3.5843 -0.1584 0.014687
10000 -6.2032 3.2395 -0.0620 0.008251
50000 -6.0557 3.0608 -0.0120 0.007391
lstsq -6.0544 3.0593 -0.0116 0.007391
50,000 descent steps to approach what one call to lstsq returns exactly. The
condition number above, 8,418, is why: gradient descent converges at a
rate set by it. When a model is linear in its parameters, solve it; keep the
gradients for the parts that are not (the cascade, the flowsheet, the
economics), which is what the rest of difflow is for.
======================================================================
Fit quality (c = 0 fit):
======================================================================
pH D_measured D_predicted Error %
-------------------------------------------------------
1.00 0.00098122 0.000992837 1.18
1.25 0.0052313 0.00564366 7.88
1.50 0.0373804 0.0320807 14.18
1.75 0.180931 0.182359 0.79
2.00 1.00748 1.0366 2.89
2.25 5.43395 5.89241 8.44
2.50 35.2863 33.4947 5.08
8. Creating an Extractant with Full Element Coverage#
For a whole separation train, define coefficients for every element you intend to carry. The shipped database has 15; this one covers 10 of them.
# Create an extractant covering ten of the fifteen database elements.
#
# One shared slope b = 3, and selectivity written entirely into `a`: adjacent
# lanthanides are 0.4 apart (beta = 2.5), with the Nd->Sm gap widened to 0.6
# because Pm is skipped and the Eu->Gd gap narrowed to 0.3 for the gadolinium
# break. a_Nd = -6.00 again puts D(Nd) = 1 at pH 2.0.
#
# Y carries no f electrons, so where it sits is a property of the EXTRACTANT,
# not of yttrium: on a phosphonic acid like PC88A it falls between Dy and Ho
# (beta(Y/Dy) = 3.14 from Torres et al. 2021), while on Cyanex 272 it drops down
# among the light-middle lanthanides, which is why that reagent is the one used
# to make yttrium. This record follows the phosphonic placement: a_Y is 0.5
# above a_Dy, for beta(Y/Dy) = 3.16.
comprehensive_extractant = create_custom_extractant(
name="ComprehensiveExtractant",
full_name="Comprehensive Custom Extractant",
formula="C12H24O4P",
molecular_weight=280.0,
ph_coefficients={
# Light REEs
"La": {"a": -7.20, "b": 3.0, "c": 0.0},
"Ce": {"a": -6.80, "b": 3.0, "c": 0.0},
"Pr": {"a": -6.40, "b": 3.0, "c": 0.0},
"Nd": {"a": -6.00, "b": 3.0, "c": 0.0},
# Middle REEs
"Sm": {"a": -5.40, "b": 3.0, "c": 0.0},
"Eu": {"a": -5.00, "b": 3.0, "c": 0.0},
"Gd": {"a": -4.70, "b": 3.0, "c": 0.0},
# Heavy REEs
"Tb": {"a": -4.30, "b": 3.0, "c": 0.0},
"Dy": {"a": -3.90, "b": 3.0, "c": 0.0},
"Y": {"a": -3.40, "b": 3.0, "c": 0.0},
},
temperature_coefficients={
"La": -1500, "Ce": -1600, "Pr": -1700, "Nd": -1800,
"Sm": -2000, "Eu": -2100, "Gd": -2200,
"Tb": -2300, "Dy": -2400, "Y": -2200,
},
pKa=3.3,
valid_ph_range=(0.5, 3.0),
cost_usd_kg=20.0,
)
# Register it
db.add_extractant("ComprehensiveExtractant", comprehensive_extractant)
print("✓ Comprehensive extractant created")
print(f" Covers {len(comprehensive_extractant.ph_coefficients)} of "
f"{len(list_ree_elements())} database elements")
print(f" Elements: {', '.join(comprehensive_extractant.ph_coefficients.keys())}")
# Test with larger element set
series = ("La", "Ce", "Pr", "Nd", "Sm", "Eu", "Gd", "Tb", "Dy", "Y")
dist_comp = REEDistribution(
extractant="ComprehensiveExtractant",
elements=series,
concentration=0.5,
)
D_comp = dist_comp.get_D_all(pH=pH_op, T=298.15)
print(f"\nDistribution coefficients and adjacent separation factors at pH {pH_op}:")
print(f"{'Element':<10} {'D':<14} {'beta vs previous':<18}")
print("-"*44)
prev = None
for elem in series:
D = float(D_comp[elem])
beta = "" if prev is None else f"{D / prev:.2f}"
print(f"{elem:<10} {D:<14.4f} {beta:<18}")
prev = D
✓ Comprehensive extractant created
Covers 10 of 15 database elements
Elements: La, Ce, Pr, Nd, Sm, Eu, Gd, Tb, Dy, Y
Distribution coefficients and adjacent separation factors at pH 2.0:
Element D beta vs previous
--------------------------------------------
La 0.0631
Ce 0.1585 2.51
Pr 0.3981 2.51
Nd 1.0000 2.51
Sm 3.9811 3.98
Eu 10.0000 2.51
Gd 19.9526 2.00
Tb 50.1187 2.51
Dy 125.8925 2.51
Y 398.1072 3.16
9. Cleanup: Removing Custom Extractants#
Remove custom extractants from the database when done.
# Show current extractants
print("Before cleanup:")
print(f" {list_extractants()}")
# Remove custom extractants
db.remove_extractant("MyExtractant")
db.remove_extractant("ComprehensiveExtractant")
print("\nAfter cleanup:")
print(f" {list_extractants()}")
print("\n✓ Database restored to standard extractants")
Before cleanup:
['D2EHPA', 'PC88A', 'Cyanex272', 'TBP', 'naphthenic_acid', 'MyExtractant', 'ComprehensiveExtractant']
After cleanup:
['D2EHPA', 'PC88A', 'Cyanex272', 'TBP', 'naphthenic_acid']
✓ Database restored to standard extractants
Summary#
This notebook demonstrated:
Creating custom extractants with
create_custom_extractant()Registering extractants with the global database
Using custom extractants in distribution models and unit operations
Comparing performance vs. standard extractants, inside both validity windows
Calibrating parameters from data — and diagnosing a fit that has failed
Best practices for coefficient selection
Guidelines for Custom Extractants#
pH coefficients:
b = 3for a trivalent cation exchanger. It is stoichiometry, not a knob. Share it across every element on the record.Put the selectivity in
a: \(\log_{10}\beta_{ij} = a_i - a_j\), so a shared slope makes \(\beta\) independent of pH, temperature and concentration.Adjacent lanthanides differ by \(\Delta a \approx 0.2\)–\(0.5\) (\(\beta \approx 1.5\)–\(3\)), rising across the series. Widen the gap where a member is skipped (Pm) and narrow it at the gadolinium break.
Keep
c = 0unless data over a wide pH range demands curvature. Seven points cannot resolve it (section 7).Declare
valid_ph_rangeas the range your coefficients are actually good for, andreference_concentrationas the basis ofa. Both are load-bearing:get_Dwarns outside the window, and at \(b = 3\) one pH unit outside it is three decades of error.
Where Y goes is a property of the extractant, not of yttrium — phosphoric and phosphonic acids place it in the heavy group, Cyanex 272 does not.
Validation:
Fit
log10(D), neverD. Squared error inDis a fit to your largest point.The correlation is linear in its parameters — use
lstsq, and report standard errors. A coefficient smaller than its own error bar is not a result.Check the recovered
bagainst 3. A slope well below it means the model is absorbing physics it does not contain.
Next Steps:
Use custom extractants in flowsheet optimization
Combine with economic analysis for cost comparisons
Explore extractant mixtures (weighted average of coefficients)