Custom REE Elements: Adding New Elements to the Database

Custom REE Elements: Adding New Elements to the Database#

The difflow_ree element database covers all fifteen rare earths (La, Ce, Pr, Nd, Sm, Eu, Gd, Tb, Dy, Y, Ho, Er, Tm, Yb, Lu) and five extractant systems (D2EHPA, PC88A, Cyanex 272, TBP, naphthenic acid). What it does not cover uniformly is the extraction chemistry: four of the five extractants carry distribution coefficients for only ten of those fifteen.

So “missing” means two different things, and they take two different APIs:

What is missing	Symptom	Fix
The element itself — Sc, Pm, or a non-REE impurity like Fe or Th	`KeyError` from `ree_db.get(...)`	`create_custom_element` + `ree_db.add_element`
The element’s coefficients for one extractant — Ho, Er, Tm, Yb, Lu with PC88A	`ValueError` from `REEDistribution`, naming the gap	`ext_db.add_element_to_extractant`

Physical properties alone do not make an element usable. Nothing can compute a D for it until some extractant has coefficients for it.

What you will learn#

Read the coverage gap off ext_db.coverage() before a run, rather than mid-solve
Register a genuinely new element with create_custom_element and add_element (scandium)
Fit an existing element into an existing correlation with add_element_to_extractant (holmium into PC88A) — and why b is not yours to choose
Add separation factors with add_pair, and why you usually should not
Use the new element in distribution, extraction and gradient calculations
Clean up custom data when done

import math
import textwrap

import jax
import jax.numpy as jnp

jax.config.update("jax_enable_x64", True)

from difflow_ree import (
    # Custom element creation
    create_custom_element,
    # Database singletons
    get_ree_database,
    get_extractant_database,
    get_sf_database,
    # Equilibrium
    REEDistribution,
    # Units
    REEExtractor,
    REEExtractorParams,
)
from difflow.streams import make_stream, get_flows

1. What the database already covers#

coverage() answers “which elements can this extractant give me a D for?” before a flowsheet finds out mid-solve from a KeyError raised while iterating stages (#269).

ree_db = get_ree_database()
ext_db = get_extractant_database()
sf_db = get_sf_database()

elements = ree_db.list_elements()
extractants = ext_db.list_extractants()
print(f"Built-in elements ({len(elements)}):", elements)
print(f"Built-in extractants ({len(extractants)}):", extractants)
print()

for group in ("light", "middle", "heavy"):
    print(f"  {group}: {ree_db.list_by_group(group)}")
print()

# The gap that matters is in the correlations, not the property table.
print("Distribution-coefficient coverage:")
print(ext_db.coverage().as_text())

Built-in elements (15): ['La', 'Ce', 'Pr', 'Nd', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Y', 'Ho', 'Er', 'Tm', 'Yb', 'Lu']
Built-in extractants (5): ['D2EHPA', 'PC88A', 'Cyanex272', 'TBP', 'naphthenic_acid']

  light: ['La', 'Ce', 'Pr', 'Nd']
  middle: ['Sm', 'Eu', 'Gd']
  heavy: ['Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Y']

Distribution-coefficient coverage:
  extractant       covered  missing
  D2EHPA            10/15   Ho, Er, Tm, Yb, Lu
  PC88A             10/15   Ho, Er, Tm, Yb, Lu
  Cyanex272         10/15   Ho, Er, Tm, Yb, Lu
  TBP               10/15   Ho, Er, Tm, Yb, Lu
  naphthenic_acid   15/15   -

2. An element the database does not have: scandium#

Scandium is a rare earth by every classification that counts, and it is not one of the fifteen above. It has no 4f shell, and Sc(III) has an ionic radius of 74.5 pm — well inside Lu(III)’s 86.1 pm — so it sits past the end of the lanthanide contraction rather than on it. It is also a real feed component: most scandium on the market is a by-product of titanium-dioxide residues, bauxite residue and uranium leach liquors, which is to say it turns up in exactly the streams this package models.

The physical properties below are standard reference values. The group argument is a three-way enumeration (light / middle / heavy) that scandium does not honestly belong to; heavy is the least-wrong slot, and the choice affects nothing but list_by_group.

# Sources: CRC Handbook (physical properties), Shannon 1976 (ionic radius)
sc = create_custom_element(
    symbol="Sc",
    name="Scandium",
    atomic_number=21,
    atomic_weight=44.956,     # g/mol
    ionic_radius_pm=74.5,     # pm, Sc3+ CN=6; Lu3+ is 86.1, so Sc is smaller still
    density=2.985,            # g/cm3
    melting_point=1814,       # K
    group="heavy",            # see above -- Sc fits none of the three
    oxide_formula="Sc2O3",
    oxide_mw=137.91,          # g/mol
    price_usd_kg=3700.0,      # Sc2O3 99.99%; order of magnitude only, and volatile
)

ree_db.add_element("Sc", sc)

print(f"Added {sc.name} ({sc.symbol}), Z={sc.atomic_number}")
print(f"  Ionic radius: {sc.ionic_radius_pm} pm")
print(f"  Elements now: {len(ree_db.list_elements())}")
print()

# Properties are not enough: nothing can compute a D for Sc yet.
try:
    REEDistribution(extractant="PC88A", elements=("Dy", "Sc"), concentration=0.5)
except ValueError as err:
    print("REEDistribution refuses it, and says why:")
    print()
    print(textwrap.fill(str(err), 78))

Added Scandium (Sc), Z=21
  Ionic radius: 74.5 pm
  Elements now: 16

REEDistribution refuses it, and says why:

Extractant 'PC88A' has no coefficients for Sc (mechanism='cation_exchange'),
so no D can be computed for it. It covers: La, Ce, Pr, Nd, Sm, Eu, Gd, Tb, Dy,
Y. Either drop that element from the list, pick an extractant that covers them
(difflow_ree.get_extractant_database().coverage() reports the gaps across the
database), or add coefficients with add_element_to_extractant() -- with a
source, since extending a fitted correlation to new elements is a refit rather
than an interpolation (#269).

3. Extraction coefficients: holmium into PC88A#

Holmium is in the element database. What PC88A has no entry for is its distribution coefficient. The PC88A block was refit in #270 against Tanaka (2021), which measured La, Nd, Sm, Dy and Y; five more entries are interpolated across atomic number, and Ho, Er, Tm, Yb and Lu were left out rather than extrapolated past the end of the measured range.

The model is log10(D) = a + b*pH + c*pH^2, with a temperature correction d*(1/T - 1/T_ref). Two of those four coefficients are not free parameters:

b = 3 is stoichiometry, not a fit. Cation exchange runs

\[\mathrm{RE}^{3+} + 3\,\overline{(\mathrm{HA})_2} \rightleftharpoons \overline{\mathrm{RE}(\mathrm{HA}_2)_3} + 3\,\mathrm{H}^+\]

so mass action puts exactly three protons on the right and fixes \(d\log_{10}D/d\mathrm{pH} = 3\) — a thousandfold per pH unit. Every element in the refit PC88A block has b = 3.0.
c = 0 unless the data force otherwise. A quadratic term with no support bends the correlation outside the fitted window, and it does something worse besides: with one shared b, \(\log_{10}\beta_{ij} = a_i - a_j\) exactly, so every separation factor is a constant. Give Ho a different b or a nonzero c and that identity is gone — Ho’s factors start drifting with pH while everyone else’s stay put.

That leaves exactly one number to choose: a. Anchor it on measured neighbours. T21 gives a(Dy) = -0.3905 and a(Y) = +0.1058, i.e. beta(Y/Dy) = 3.14 with Y above Dy — in this system Y sits between Ho and Er, and the Er/Y factor that implies (about 1.5) matches Li (2019)’s measured 1.4–1.6. Splitting the Dy–Y interval so beta(Y/Ho) also lands near 1.5 puts Ho at a = -0.07, giving beta(Ho/Dy) = 2.09 — inside the 2.0–2.2 band reported for PC88A adjacent heavy pairs.

That is an interpolation between two measured points, not a measurement. Good enough to size a cascade; not good enough to publish.

pc88a = ext_db.get("PC88A")
a_Dy = pc88a.ph_coefficients["Dy"].a
a_Y = pc88a.ph_coefficients["Y"].a
a_Gd = pc88a.ph_coefficients["Gd"].a

A_HO = -0.07  # interpolated between the two measured anchors -- see above

ext_db.add_element_to_extractant(
    "PC88A",
    "Ho",
    ph_coefficients={
        "a": A_HO,   # the only chosen number here
        "b": 3.0,    # stoichiometric, and shared with every other element
        "c": 0.0,    # no quadratic term: none of the measured elements needs one
    },
    # Still hand-tuned, like every temperature coefficient in this block: T21 is
    # isothermal at 298 K and carries no information about d. This continues the
    # existing table's -100 per element step (Dy is -2300).
    temperature_coefficient=-2400,
)

print(f"a(Dy) = {a_Dy:+.4f}   a(Ho) = {A_HO:+.4f}   a(Y) = {a_Y:+.4f}")
print(f"  beta(Ho/Dy) = 10**(a_Ho - a_Dy) = {10 ** (A_HO - a_Dy):.2f}")
print(f"  beta(Y/Ho)  = 10**(a_Y - a_Ho)  = {10 ** (a_Y - A_HO):.2f}")
print()
print("Ho in PC88A: ", "Ho" in pc88a.ph_coefficients)
print("Ho in D2EHPA:", "Ho" in ext_db.get("D2EHPA").ph_coefficients)
print()
print(ext_db.coverage().as_text())

a(Dy) = -0.3905   a(Ho) = -0.0700   a(Y) = +0.1058
  beta(Ho/Dy) = 10**(a_Ho - a_Dy) = 2.09
  beta(Y/Ho)  = 10**(a_Y - a_Ho)  = 1.50

Ho in PC88A:  True
Ho in D2EHPA: False

  extractant       covered  missing
  D2EHPA            10/16   Ho, Er, Tm, Yb, Lu, Sc
  PC88A             11/16   Er, Tm, Yb, Lu, Sc
  Cyanex272         10/16   Ho, Er, Tm, Yb, Lu, Sc
  TBP               10/16   Ho, Er, Tm, Yb, Lu, Sc
  naphthenic_acid   15/16   Sc

4. Separation factors: derived by default, authored by exception#

Since #265 the shipped separation factors are computed from the distribution correlations, not tabulated. separation_factors.yaml lists which pairs an extractant is characterised for; the numbers come from the coefficients, so they cannot drift out of step with them. SeparationFactorData.derived names the pairs that were computed.

add_pair is the escape hatch, and it takes the value as given. Use it for a measured factor the correlation cannot reproduce. For anything the correlation already describes — which, now that Ho is in PC88A, includes every Ho pair — putting the element in the extractant is the better move: one description of the physics instead of two that can disagree. So the cell below reads its factors off the coefficients rather than inventing them.

The same applies to stage counts. get_stages_needed derives the Fenske minimum

\[N_{\min} = \frac{2\ln 99}{\left|\ln \beta\right|}\]

from the pair’s factor — the equimolar binary feed taken to 99% purity in both products, at total reflux. This was an 18-entry hand-authored table until #270, and it matched no beta in the package. Passing stages_99= to add_pair still overrides the derivation and silently wins, so pass it only when you mean to (§8 shows what that looks like).

# Convention: "heavier_lighter", SF = D_heavier / D_lighter.
# Read the factors off the coefficients added in section 3.
beta_Ho_Dy = 10 ** (A_HO - a_Dy)
beta_Y_Ho = 10 ** (a_Y - A_HO)
beta_Ho_Gd = 10 ** (A_HO - a_Gd)

sf_db.add_pair("PC88A", "Ho_Dy", beta_Ho_Dy, adjacent=True)
sf_db.add_pair("PC88A", "Y_Ho", beta_Y_Ho, adjacent=True)
sf_db.add_pair("PC88A", "Ho_Gd", beta_Ho_Gd, adjacent=False)  # non-adjacent group pair

print(f"SF(Ho/Dy) with PC88A: {sf_db.get_sf('PC88A', 'Ho_Dy'):.3f}")
print(f"SF(Y/Ho)  with PC88A: {sf_db.get_sf('PC88A', 'Y_Ho'):.3f}")
print(f"SF(Ho/Gd) with PC88A: {sf_db.get_sf('PC88A', 'Ho_Gd'):.3f}")
print()

data = sf_db.get("PC88A")
print("Derived from the correlation, or authored here?")
for pair in ("Dy_Tb", "Y_Dy", "Ho_Dy", "Ho_Gd"):
    print(f"  {pair:<7} {'derived' if pair in data.derived else 'AUTHORED'}")
print()

n_min = sf_db.get_stages_needed("PC88A", "Ho_Dy")
print(f"Fenske minimum for Ho/Dy at beta = {beta_Ho_Dy:.3f}: {n_min} stages")
print(f"  check: 2*ln(99)/ln(beta) = {2 * math.log(99) / math.log(beta_Ho_Dy):.2f}, rounded up")
print("  It is a thermodynamic floor at total reflux. A real cascade at a")
print("  finite solvent ratio needs several times more.")

SF(Ho/Dy) with PC88A: 2.092
SF(Y/Ho)  with PC88A: 1.499
SF(Ho/Gd) with PC88A: 26.953

Derived from the correlation, or authored here?
  Dy_Tb   derived
  Y_Dy    derived
  Ho_Dy   AUTHORED
  Ho_Gd   AUTHORED

Fenske minimum for Ho/Dy at beta = 2.092: 13 stages
  check: 2*ln(99)/ln(beta) = 12.45, rounded up
  It is a thermodynamic floor at total reflux. A real cascade at a
  finite solvent ratio needs several times more.

5. Using Ho in distribution calculations#

Ho now behaves like any built-in element. Watch the pH: the refit PC88A block is valid over pH 0.1–2.5, and for the heavy rare earths the interesting part of that is the acid end. At pH 1 the correlation already gives D(Ho) ~ 850 at trace loading.

dist = REEDistribution(
    extractant="PC88A",
    elements=("Gd", "Dy", "Ho", "Y"),
    concentration=0.5,
)

print("Distribution coefficients, PC88A 0.5 M, 25 C, trace loading")
print("-" * 56)
print(f"{'pH':>5} {'D(Gd)':>11} {'D(Dy)':>11} {'D(Ho)':>11} {'D(Y)':>11}")
for pH in (0.1, 0.5, 1.0, 2.0):
    D = dist.get_D_all(pH=pH, T=298.15)
    row = " ".join(f"{float(D[e]):11.4g}" for e in ("Gd", "Dy", "Ho", "Y"))
    print(f"{pH:>5.1f} {row}")
print()

print("Separation factors -- the extraction order is Y > Ho > Dy > Gd")
print("-" * 56)
print(f"{'pH':>5} {'Ho/Dy':>9} {'Y/Ho':>9} {'Ho/Gd':>9}")
for pH in (0.1, 0.5, 1.0, 2.0):
    D = dist.get_D_all(pH=pH, T=298.15)
    print(
        f"{pH:>5.1f} {float(D['Ho'] / D['Dy']):9.3f}"
        f" {float(D['Y'] / D['Ho']):9.3f} {float(D['Ho'] / D['Gd']):9.2f}"
    )
print()
print("Not one of them moves with pH. That is the shared b = 3.0, not a")
print("coincidence: log10(beta_ij) = a_i - a_j exactly, so beta is a constant.")

Distribution coefficients, PC88A 0.5 M, 25 C, trace loading
--------------------------------------------------------
   pH       D(Gd)       D(Dy)       D(Ho)        D(Y)
  0.1     0.06301      0.8119       1.698       2.546
  0.5      0.9986       12.87       26.92       40.35
  1.0       31.58       406.9       851.1        1276
  2.0   3.158e+04   4.069e+05   8.511e+05   1.276e+06

Separation factors -- the extraction order is Y > Ho > Dy > Gd
--------------------------------------------------------
   pH     Ho/Dy      Y/Ho     Ho/Gd
  0.1     2.092     1.499     26.95
  0.5     2.092     1.499     26.95
  1.0     2.092     1.499     26.95
  2.0     2.092     1.499     26.95

Not one of them moves with pH. That is the shared b = 3.0, not a
coincidence: log10(beta_ij) = a_i - a_j exactly, so beta is a constant.

6. Using Ho in a multi-stage extraction unit#

Custom elements work with the unit operation models too. Here we run a 5-stage extraction of a Gd/Ho/Y mixture with PC88A.

The operating point is pH 0.1 at S/F = 0.4 — lower on both counts than this notebook used before #270. That is a consequence of the refit rather than a modelling preference: at trace loading PC88A takes the heavy rare earths essentially completely above pH 0.5, so a calculation there answers “everything extracts” and displays no selectivity at all. Two things bring a real circuit back into range: the acid end of the measured window, and organic loading capacity (include_loading, or solve_free_extractant). This cell uses the first. The feed here is far too dilute for the second to bite — 0.045 mol/s of rare earth against 0.4 mol/s of extractant.

# Extract at the acid end of PC88A's measured window
pH_op = 0.1

params = REEExtractorParams(
    n_stages=5,
    extractant="PC88A",
    elements=("Gd", "Ho", "Y"),
    pH=pH_op,
    include_loading=False,  # feed is dilute enough that it changes nothing here
)
extractor = REEExtractor(params)

# Feed stream: mixed Gd/Ho/Y solution
feed = make_stream(
    flows={"H2O": 10.0, "Gd": 0.02, "Ho": 0.01, "Y": 0.015},
    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 PC88A 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). Total organic flow 4.0 mol/s against a 10.0 mol/s feed is
# S/F = 0.4, split 0.4 PC88A / 3.6 kerosene.
solvent = make_stream(
    flows={"PC88A": 0.4, "kerosene": 3.6, "Gd": 0.0, "Ho": 0.0, "Y": 0.0},
    T=298.15, P=101325.0,
)

# Run extraction
raffinate, extract, info = extractor(feed, solvent)

# Results
feed_flows = get_flows(feed)
ext_flows = get_flows(extract)

print(f"Extraction results (5 stages, PC88A, pH {pH_op}, S/F = 0.4):")
print("-" * 55)
print(f"{'Element':<10} {'Feed (mol/s)':<15} {'Extract':<15} {'Recovery %'}")
print("-" * 55)
for elem in ["Gd", "Ho", "Y"]:
    f_val = float(feed_flows[elem])
    e_val = float(ext_flows[elem])
    rec = (e_val / f_val) * 100 if f_val > 0 else 0
    print(f"{elem:<10} {f_val:<15.4f} {e_val:<15.4f} {rec:.1f}")

Extraction results (5 stages, PC88A, pH 0.1, S/F = 0.4):
-------------------------------------------------------
Element    Feed (mol/s)    Extract         Recovery %
-------------------------------------------------------
Gd         0.0200          0.0005          2.5
Ho         0.0100          0.0064          64.2
Y          0.0150          0.0126          83.9

7. Gradients work through custom elements#

The custom element goes through the same differentiable model, so JAX automatic differentiation works on it unchanged. Below is the gradient of Ho recovery with respect to extraction pH — followed by a check on what that number does and does not license you to say.

from jax import grad


def ho_recovery(pH):
    """Compute Ho recovery as a function of extraction pH."""
    params = REEExtractorParams(
        n_stages=5,
        extractant="PC88A",
        elements=("Gd", "Ho", "Y"),
        pH=pH,
        include_loading=False,
    )
    extractor = REEExtractor(params)

    feed = make_stream(
        flows={"H2O": 10.0, "Gd": 0.02, "Ho": 0.01, "Y": 0.015},
        T=298.15, P=101325.0,
    )
    # Same solvent as above: PC88A + kerosene, total organic flow 4.0 (#191/#192)
    solvent = make_stream(
        flows={"PC88A": 0.4, "kerosene": 3.6, "Gd": 0.0, "Ho": 0.0, "Y": 0.0},
        T=298.15, P=101325.0,
    )

    _, extract, _ = extractor(feed, solvent)
    ext_flows = get_flows(extract)
    return ext_flows["Ho"] / 0.01  # recovery fraction


d_recovery_d_pH = grad(ho_recovery)

rec = float(ho_recovery(pH_op))
drec = float(d_recovery_d_pH(pH_op))

print(f"At pH {pH_op}:")
print(f"  Ho recovery:     {rec:.4f} ({rec * 100:.1f}%)")
print(f"  d(recovery)/dpH: {drec:.4f} per pH unit")
print()

# A derivative is a statement about an infinitesimal step. Recovery is bounded
# above by 1, so the tangent line has to fail somewhere; the only question is
# how soon. Measure that instead of quoting the slope times a finite step.
print("Tangent line against the model it was taken from:")
print(f"{'delta pH':>9} {'linear':>10} {'actual':>10}")
for d_pH in (0.02, 0.05, 0.10):
    linear = rec + drec * d_pH
    actual = float(ho_recovery(pH_op + d_pH))
    print(f"{d_pH:>9.2f} {linear * 100:>9.1f}% {actual * 100:>9.1f}%")
print()
print("By delta pH = 0.10 the tangent is predicting more than 100% recovery,")
print("which no separation achieves. Use the gradient to point an optimizer")
print("downhill -- not to quote a finite change.")

At pH 0.1:
  Ho recovery:     0.6420 (64.2%)
  d(recovery)/dpH: 3.5963 per pH unit

Tangent line against the model it was taken from:
 delta pH     linear     actual
     0.02      71.4%      71.4%
     0.05      82.2%      81.4%
     0.10     100.2%      93.1%

By delta pH = 0.10 the tangent is predicting more than 100% recovery,
which no separation achieves. Use the gradient to point an optimizer
downhill -- not to quote a finite change.

8. Correcting or removing custom data#

If you need to update coefficients — after getting better literature data, say — remove and re-add. All custom data modifications are runtime-only and do not touch the YAML files on disk.

Note what does not change in a revision: b stays at 3.0 and c stays at 0.0. A better fit to Ho is a better a. If a refit genuinely wants a different slope, the thing to check first is the loading, the aqueous complexation or the diluent — not the stoichiometry.

# Update extraction coefficients: remove and re-add
A_HO_REVISED = -0.10

ext_db.remove_element_from_extractant("PC88A", "Ho")
ext_db.add_element_to_extractant(
    "PC88A", "Ho",
    ph_coefficients={"a": A_HO_REVISED, "b": 3.0, "c": 0.0},  # revised a; b, c pinned
    temperature_coefficient=-2400,
)
print(f"Revised a(Ho): {ext_db.get('PC88A').ph_coefficients['Ho'].a:+.4f}")

# Update element properties (e.g., a corrected price)
updated_sc = create_custom_element(
    symbol="Sc", name="Scandium", atomic_number=21,
    atomic_weight=44.956, ionic_radius_pm=74.5, density=2.985,
    melting_point=1814, group="heavy", oxide_formula="Sc2O3",
    oxide_mw=137.91, price_usd_kg=4200.0,  # updated price
)
ree_db.update_element("Sc", updated_sc)
print(f"Updated Sc price: ${ree_db.get('Sc').price_usd_kg}/kg")

# Keep the separation factor in step with the revised coefficients
beta_revised = 10 ** (A_HO_REVISED - a_Dy)
sf_db.remove_pair("PC88A", "Ho_Dy")
sf_db.add_pair("PC88A", "Ho_Dy", beta_revised)
print(f"Revised SF(Ho/Dy): {sf_db.get_sf('PC88A', 'Ho_Dy'):.3f}")
print(f"  Fenske minimum, derived: {sf_db.get_stages_needed('PC88A', 'Ho_Dy')} stages")

# And what an authored override looks like. Nothing checks it against beta.
sf_db.remove_pair("PC88A", "Ho_Dy")
sf_db.add_pair("PC88A", "Ho_Dy", beta_revised, stages_99=18)
print(f"  with stages_99=18:       {sf_db.get_stages_needed('PC88A', 'Ho_Dy')} stages"
      "  <- the override wins, silently")

Revised a(Ho): -0.1000
Updated Sc price: $4200.0/kg
Revised SF(Ho/Dy): 1.952
  Fenske minimum, derived: 14 stages
  with stages_99=18:       18 stages  <- the override wins, silently

9. Cleanup#

Remove Ho and Sc from all databases to restore the original state. Since the databases are singletons, this matters if other code in your session relies on the default element list.

# Remove Ho from extractant coefficients
ext_db.remove_element_from_extractant("PC88A", "Ho")

# Remove Ho separation factor pairs
for pair in ["Ho_Dy", "Y_Ho", "Ho_Gd"]:
    try:
        sf_db.remove_pair("PC88A", pair)
    except KeyError:
        pass

# Remove the custom element (Ho was never a custom element -- it ships with the
# property database. Only Sc was added here.)
ree_db.remove_element("Sc")

print("Restored databases:")
print("  Elements:  ", ree_db.list_elements())
print("  Heavy REEs:", ree_db.list_by_group("heavy"))
print()
print(ext_db.coverage().as_text())

Restored databases:
  Elements:   ['La', 'Ce', 'Pr', 'Nd', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Y', 'Ho', 'Er', 'Tm', 'Yb', 'Lu']
  Heavy REEs: ['Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Y']

  extractant       covered  missing
  D2EHPA            10/15   Ho, Er, Tm, Yb, Lu
  PC88A             10/15   Ho, Er, Tm, Yb, Lu
  Cyanex272         10/15   Ho, Er, Tm, Yb, Lu
  TBP               10/15   Ho, Er, Tm, Yb, Lu
  naphthenic_acid   15/15   -

Summary#

This notebook demonstrated the full workflow for extending the REE database:

Step	Function	Purpose
Check the gap	`ext_db.coverage()`	Which elements each extractant can give a `D` for, before a run
Create element	`create_custom_element(...)`	Build an `REEElement` from physical properties
Register element	`ree_db.add_element(symbol, elem)`	Make it available to all database queries
Add extraction data	`ext_db.add_element_to_extractant(...)`	pH/temperature coefficients for one extractant
Add SF data	`sf_db.add_pair(...)`	An authored separation factor, for what the correlation cannot reproduce
Update data	`remove_*` then re-add	Correct values when better literature data arrives

All modifications are runtime-only: they do not alter the YAML files on disk, so restarting Python restores the defaults.

Three things worth carrying out of here:

A property record is not an extraction model. Ho ships with the database and still cannot be simulated with PC88A. The coverage() report is the thing to check, not list_elements().
b and c are not free parameters. b = 3 is the mass-action stoichiometry of cation exchange and c = 0 is the default; the only number to fit for a new element is a. Break the shared slope and every separation factor involving that element starts drifting with pH while the rest stay constant.
Derived beats authored. Separation factors and Fenske stage counts are computed from the correlations (#265, #270) precisely so they cannot go stale. add_pair(..., stages_99=...) overrides the derivation and wins silently — an authored number is a commitment to maintain it by hand.

For creating entirely new extractants rather than adding elements to existing ones, see notebook 21_custom_extractants.ipynb. For fitting a to real D-vs-pH data — including why to fit log10 D and never D — see its section 8.