Note
Go to the end to download the full example code.
Starrett–Saumon electronic levels and ionization
This benchmark compares current Otter ion-sphere (IS) and experimental self-consistent (SC) diagnostics with Tables I–III of Starrett and Saumon [2013] (doi: 10.1103/PhysRevE.87.013104). It was added specifically to make pressure ionization auditable before interpreting the \(\bar Z\) difference in Johnson et al. (2025).
The article’s reported results are not isolated ion-sphere calculations. Its average-atom and two-component-plasma (QTCP or TFTCP) equations are iterated self-consistently; Appendix B identifies the ion-sphere model as the initial guess. Table I is discussed as part of that coupled model and its 15-eV 3s entry is \(E=-0.0125\) Hartree (about -0.340 eV), not -0.0125 eV.
For a bound level of energy \(E<0\), Otter uses Eq. (81),
\[M(E)=\operatorname{erf}\!\left[-2\sqrt{\ln 2}\,E/\gamma\right],\]
where \(\gamma\) is obtained from the continuum phase shifts through the transport-cross-section, Ziman-conductivity, and Drude relaxation-time path. These sparse diagnostics are clearer as tables than as plots, so this gallery uses direct numerical comparisons throughout. Otter values are rounded to three significant digits; published values retain the precision of the article. Every \(E\) and \(\gamma\) entry uses Hartree, while \(M(E)\) and both ionization definitions are dimensionless.
| State | \(E\) [Ha] | \(M(E)\) | ||||
|---|---|---|---|---|---|---|
| Paper SC | Otter-IS | Otter-SC | Paper SC | Otter-IS | Otter-SC | |
| 1s | -54.6 | -54.6 | -54.6 | 1.00 | 1.00 | 1.00 |
| 2s | -3.41 | -3.39 | -3.40 | 1.00 | 1.00 | 1.00 |
| 2p | -2.04 | -2.02 | -2.02 | 1.00 | 1.00 | 1.00 |
| 3s | unbound | unbound | unbound | — | — | — |
| State | \(E\) [Ha] | \(M(E)\) | ||||
|---|---|---|---|---|---|---|
| Paper SC | Otter-IS | Otter-SC | Paper SC | Otter-IS | Otter-SC | |
| 1s | -54.9 | -54.8 | -54.8 | 1.00 | 1.00 | 1.00 |
| 2s | -3.60 | -3.50 | -3.57 | 1.00 | 1.00 | 1.00 |
| 2p | -2.23 | -2.14 | -2.21 | 1.00 | 1.00 | 1.00 |
| 3s | -0.0125 | unbound | -0.00763 | 0.134 | — | 0.0842 |
| \(T\) [eV] | Paper SC | Otter-IS | Otter-SC |
|---|---|---|---|
| 2 | 0.0698 | 0.0367 | 0.0617 |
| 15 | 0.174 | 0.109 | 0.170 |
As a formula-only audit, evaluating Otter’s Eq. (81) with the published 15-eV 3s inputs \(E=-0.0125\) Ha and \(\gamma=0.174\) Ha gives \(M=0.1343\), consistent with the published 0.134. This is not a fourth physical model, so it is reported as an audit rather than an extra table model. The accepted NPZ retains full precision; the executable summary follows the same three-significant-digit presentation. Tables II and III use two distinct ionization definitions:
\[Z^*=n_e^0/n_I^0, \qquad
\bar Z=Z-\int n_e^{\rm ion}(r)\,d^3r.\]
The following tables compare both definitions directly. The \(\Gamma_{\rm OCP}\) and \(\Gamma_{\rm TCP}\) columns are transcribed paper parameters, not Otter observables; this benchmark does not manufacture nominally corresponding Otter values.
| \(T\) [eV] | \(Z^*\) | \(\bar Z\) | Paper coupling | |||||
|---|---|---|---|---|---|---|---|---|
| Paper | Otter-IS | Otter-SC | Paper | Otter-IS | Otter-SC | \(\Gamma_{\rm OCP}\) | \(\Gamma_{\rm TCP}\) | |
| 2 | 1.98 | 2.09 | 2.09 | 3.00 | 3.00 | 3.00 | 41.0 | 5.05 |
| 6 | 2.11 | 2.19 | 2.19 | 3.00 | 3.00 | 3.00 | 13.6 | 2.04 |
| 10 | 2.24 | 2.33 | 2.33 | 3.00 | 3.03 | 3.02 | 8.14 | 1.60 |
| 15 | 2.51 | 2.56 | 2.56 | 3.18 | 3.24 | 3.20 | 6.12 | 1.34 |
| \(\rho\) [g cm\(^{-3}\)] | \(n_I^0\) [Bohr\(^{-3}\)] | \(T\) [eV] | \(Z^*\) | \(\bar Z\) | Paper coupling | |||||
|---|---|---|---|---|---|---|---|---|---|---|
| Paper | Otter-IS | Otter-SC | Paper | Otter-IS | Otter-SC | \(\Gamma_{\rm OCP}\) | \(\Gamma_{\rm TCP}\) | |||
| 22.5 | 0.0360 | 10 | 5.85 | 5.85 | 5.85 | 8.78 | 8.81 | 8.74 | 112 | 14.3 |
| 34.5 | 0.0551 | 100 | 9.54 | 9.54 | 9.54 | 11.6 | 11.9 | 11.5 | 22.4 | 4.85 |
| 39.65 | 0.0634 | 1000 | 20.4 | 20.4 | 20.4 | 21.7 | 21.9 | 21.6 | 8.22 | 3.32 |
| 34.37 | 0.0549 | 5000 | 25.1 | 25.1 | 25.1 | 25.5 | 25.6 | 25.6 | 2.18 | 1.35 |
Otter-IS is the converged production ion-sphere workflow. Otter-SC is the experimental feedback path based on Sec. 2.4 of Starrett and Saumon [2014]: it feeds QOZ/HNC ion structure back into the average atom while holding the converged IS chemical potential fixed. This is a useful controlled comparison, but it is not claimed to reproduce every detail of the simultaneous 2013 QTCP/TFTCP solver. Al uses quantum orbitals; Fe uses Thomas–Fermi electrons, matching the model named for Table III.
Set USE_PRECOMPUTED_DATA = False below to recompute all eight independent
IS/SC pairs. New calculations are staged under benchmarks/outputs and
never overwrite the accepted, checksummed baseline.
Using checksummed, accepted Otter IS/SC states.
Al bound levels: all energies E are in Hartree
T[eV] shell paper E[Ha] IS E[Ha] SC E[Ha] M paper M IS M SC
2 1s -54.6 -54.6 -54.6 1.00 1.00 1.00
2 2s -3.41 -3.39 -3.40 1.00 1.00 1.00
2 2p -2.04 -2.02 -2.02 1.00 1.00 1.00
2 3s unbound unbound unbound -- -- --
15 1s -54.9 -54.8 -54.8 1.00 1.00 1.00
15 2s -3.60 -3.50 -3.57 1.00 1.00 1.00
15 2p -2.23 -2.14 -2.21 1.00 1.00 1.00
15 3s -0.0125 unbound -0.00763 0.134 -- 0.0842
Eq. (81) formula audit at 15 eV 3s: paper inputs give M=0.134; published M=0.134.
Scattering width gamma: all values are in Hartree
T[eV] paper gamma[Ha] IS gamma[Ha] SC gamma[Ha]
2 0.0698 0.0367 0.0617
15 0.174 0.109 0.170
Mean ionization: all values are dimensionless
state Z* paper Z* IS Z* SC Zbar paper Zbar IS Zbar SC
al_rho2p7_te2_qm 1.98 2.09 2.09 3.00 3.00 3.00
al_rho2p7_te6_qm 2.11 2.19 2.19 3.00 3.00 3.00
al_rho2p7_te10_qm 2.24 2.33 2.33 3.00 3.03 3.02
al_rho2p7_te15_qm 2.51 2.56 2.56 3.18 3.24 3.20
fe_rho22.5_te10_tf 5.85 5.85 5.85 8.78 8.81 8.74
fe_rho34.5_te100_tf 9.54 9.54 9.54 11.6 11.9 11.5
fe_rho39.65_te1000_tf 20.4 20.4 20.4 21.7 21.9 21.6
fe_rho34.37_te5000_tf 25.1 25.1 25.1 25.5 25.6 25.6
from __future__ import annotations
from concurrent.futures import ProcessPoolExecutor, as_completed
import csv
import hashlib
import json
from pathlib import Path
import subprocess
import time
from typing import Any
import numpy as np
from otter import PlasmaWorkflowConfig, solve_plasma_workflow
from otter.data.helpers import ion_level_weight
from otter.experimental import SCFeedbackConfig, solve_sc_feedback_workflow
# =============================================================================
# User input
# =============================================================================
USE_PRECOMPUTED_DATA = True
MAX_STATE_WORKERS = 2
CONTINUUM_WORKERS_PER_STATE = 4
# =============================================================================
BENCHMARK_ID = "starrett_saumon_2013_electronic"
SCHEMA = "otter_starrett_saumon_2013_electronic_v2"
STRUCTURE_LABELS = ("is", "sc")
LEVEL_LABELS = ("1s", "2s", "2p", "3s")
AL_TEMPERATURES_EV = (2.0, 6.0, 10.0, 15.0)
FE_HUGONIOT = (
(22.5, 10.0),
(34.5, 100.0),
(39.65, 1000.0),
(34.37, 5000.0),
)
SC_CONTROLS = SCFeedbackConfig(
max_outer=16,
g_tol=5.0e-4,
v_corr_tol=5.0e-4,
v_corr_mix=0.5,
require_converged=True,
)
def repository_root() -> Path:
"""Locate the checkout when run directly or through Sphinx-Gallery."""
candidates = [Path.cwd().resolve(), *Path.cwd().resolve().parents]
source = Path(str(globals().get("__file__", Path.cwd()))).resolve()
candidates.extend([source.parent, *source.parents])
for candidate in candidates:
if (candidate / "pyproject.toml").is_file() and (
candidate / "src" / "otter"
).is_dir():
return candidate
raise FileNotFoundError("Cannot locate the Otter checkout.")
ROOT = repository_root()
REFERENCE_DIR = ROOT / "benchmarks" / "reference_data" / BENCHMARK_ID
BASELINE_DIR = ROOT / "benchmarks" / "baselines" / BENCHMARK_ID
OUTPUT_DIR = ROOT / "benchmarks" / "outputs" / BENCHMARK_ID
BASELINE_PATH = BASELINE_DIR / "electronic_states.npz"
BASELINE_MANIFEST = BASELINE_DIR / "manifest.json"
CANDIDATE_PATH = OUTPUT_DIR / "electronic_states.npz"
CANDIDATE_MANIFEST = OUTPUT_DIR / "manifest.json"
def sha256_file(path: Path) -> str:
"""Return the SHA-256 digest of one file."""
digest = hashlib.sha256()
with path.open("rb") as stream:
for block in iter(lambda: stream.read(1024 * 1024), b""):
digest.update(block)
return digest.hexdigest()
def state_definitions() -> tuple[dict[str, Any], ...]:
"""Return the four Al-QM and four Fe-TF benchmark states."""
aluminium = tuple(
{
"state_id": f"al_rho2p7_te{temperature:g}_qm",
"element": "Al",
"electronic_model": "qm",
"rho_g_cc": 2.7,
"temperature_ev": temperature,
}
for temperature in AL_TEMPERATURES_EV
)
iron = tuple(
{
"state_id": f"fe_rho{rho:g}_te{temperature:g}_tf",
"element": "Fe",
"electronic_model": "tf",
"rho_g_cc": rho,
"temperature_ev": temperature,
}
for rho, temperature in FE_HUGONIOT
)
return aluminium + iron
STATES = state_definitions()
def load_reference_tables() -> dict[str, Any]:
"""Verify and load the three published numerical tables."""
manifest = json.loads(
(REFERENCE_DIR / "manifest.json").read_text(encoding="utf-8")
)
if (
manifest.get("schema_version") != "otter_reference_manifest_v1"
or manifest.get("reference_id") != BENCHMARK_ID
or manifest["publication"].get("doi")
!= "10.1103/PhysRevE.87.013104"
):
raise ValueError("Unexpected Starrett--Saumon reference manifest.")
for record in manifest["files"]:
path = REFERENCE_DIR / str(record["path"])
if sha256_file(path) != str(record["sha256"]):
raise RuntimeError(f"Reference checksum mismatch: {path.name}.")
with (REFERENCE_DIR / "table_i_al_levels.csv").open(
encoding="utf-8", newline=""
) as stream:
level_rows = tuple(csv.DictReader(stream))
return {
"levels": level_rows,
"al_ionization": np.genfromtxt(
REFERENCE_DIR / "table_ii_al_ionization.csv",
delimiter=",",
names=True,
),
"fe_ionization": np.genfromtxt(
REFERENCE_DIR / "table_iii_fe_hugoniot.csv",
delimiter=",",
names=True,
),
}
def workflow_config(state: dict[str, Any]) -> PlasmaWorkflowConfig:
"""Return the complete IS workflow used to initialise SC feedback."""
is_tf = str(state["electronic_model"]) == "tf"
return PlasmaWorkflowConfig(
elements=[str(state["element"])],
temperature_ev=float(state["temperature_ev"]),
ion_temperature_ev=float(state["temperature_ev"]),
rho_g_cc=float(state["rho_g_cc"]),
electronic_model=str(state["electronic_model"]),
aa_overrides={
# The nondegenerate 5000-eV Fe chemical potential lies below the
# ordinary warm-dense bracket. Only the scalar search interval
# is enlarged; the converged TF equations are unchanged.
"mu_bounds": (-2000.0, 200.0) if is_tf else (-200.0, 200.0),
"cont_n_jobs": int(CONTINUUM_WORKERS_PER_STATE),
"cont_shards": int(2 * CONTINUUM_WORKERS_PER_STATE),
# SC feedback can move an s state through the E=0 threshold.
# Match that candidate to its all-space tail before allowing its
# finite-box sign to define the bound/free partition.
"bound_zero_tail_refine": not is_tf,
# The paper's 15-eV 3s level is -0.0125 Ha. Include that entire
# physical comparison window rather than the library's narrower
# near-zero diagnostic default (1e-3 Ha).
"bound_zero_tail_max_binding_ha": 3.0e-2,
},
hnc_closure_transform_tol=2.5e-3,
hnc_max_iter=500,
show_progress=False,
)
def level_label(l_value: int, radial_index: int) -> str:
"""Convert stored angular/radial indices to a spectroscopic label."""
letters = "spdfgh"
return f"{radial_index + l_value}{letters[l_value]}"
def reduce_electronic(
result: dict[str, Any], *, electronic_model: str
) -> dict[str, Any]:
"""Apply convergence gates and retain compact electronic diagnostics."""
if result.get("stage2_converged") is not True:
raise RuntimeError("full average-atom stage 2 did not converge")
threshold = str(result.get("threshold_state_status", "none")).lower()
if threshold == "unresolved":
raise RuntimeError("threshold-state representation is unresolved")
n_i = float(dict(result["meta"])["n_i_bohr3"])
zstar = float(result["n0"]) / n_i
zbar = float(result["zbar_partition"])
if not np.all(np.isfinite((n_i, zstar, zbar))):
raise RuntimeError("full AA returned a non-finite ionization diagnostic")
energies = np.full(len(LEVEL_LABELS), np.nan)
weights = np.full(len(LEVEL_LABELS), np.nan)
gamma = float(
result.get(
"ion_gamma",
dict(result.get("meta", {})).get("ion_gamma_final", np.nan),
)
)
if str(electronic_model) != "tf":
energy_table = np.asarray(result["bound_energy_ha"], dtype=float)
weight_table = np.asarray(result["bound_m"], dtype=float)
l_values = np.asarray(result["bound_l_list"], dtype=int)
n_indices = np.asarray(result["bound_n_index"], dtype=int)
for l_row, l_value in enumerate(l_values):
for column, radial_index in enumerate(n_indices):
label = level_label(int(l_value), int(radial_index))
if label not in LEVEL_LABELS:
continue
energy = float(energy_table[l_row, column])
if not np.isfinite(energy) or energy >= 0.0:
continue
index = LEVEL_LABELS.index(label)
energies[index] = energy
weights[index] = float(weight_table[l_row, column])
expected = ion_level_weight(energy, gamma)
if not np.isclose(weights[index], expected, atol=2.0e-13):
raise RuntimeError(f"stored M(E) is inconsistent for {label}")
history = list(result.get("history", ()))
final_error = float(history[-1].get("err", np.nan)) if history else np.nan
return {
"n_i_bohr3": n_i,
"r_ws_bohr": float(result["r_ws"]),
"zstar": zstar,
"zbar_partition": zbar,
"ion_gamma_ha": gamma,
"level_energy_ha": energies,
"level_m": weights,
"stage2_error": final_error,
"stage2_iters": int(result.get("stage2_iters", len(history))),
"threshold_state_status": threshold,
}
def electronic_result(workflow: dict[str, Any]) -> dict[str, Any]:
"""Return the single-species full+external result from one workflow."""
if str(workflow["electronic"]["kind"]) != "single_species":
raise ValueError("This benchmark expects one species per workflow.")
return dict(workflow["electronic"]["result"])
def solve_state(state: dict[str, Any]) -> dict[str, Any]:
"""Solve one IS workflow and continue it to experimental SC feedback."""
config = workflow_config(state)
started = time.perf_counter()
is_workflow = solve_plasma_workflow(config)
is_elapsed_s = time.perf_counter() - started
started = time.perf_counter()
sc_workflow = solve_sc_feedback_workflow(
config,
is_workflow,
feedback_cfg=SC_CONTROLS,
)
sc_extension_elapsed_s = time.perf_counter() - started
feedback = dict(sc_workflow["sc_feedback"])
return {
**state,
"is": reduce_electronic(
electronic_result(is_workflow),
electronic_model=str(state["electronic_model"]),
),
"sc": reduce_electronic(
electronic_result(sc_workflow),
electronic_model=str(state["electronic_model"]),
),
"is_elapsed_s": float(is_elapsed_s),
"sc_extension_elapsed_s": float(sc_extension_elapsed_s),
"sc_converged": bool(feedback["converged"]),
"sc_iterations": int(feedback["iterations"]),
"fixed_is_mu_ha": float(feedback["fixed_is_mu_ha"]),
}
def pack_rows(rows: list[dict[str, Any]]) -> dict[str, np.ndarray]:
"""Pack ordered IS/SC rows into one pickle-free portable archive."""
by_id = {str(row["state_id"]): row for row in rows}
ordered = [by_id[str(state["state_id"])] for state in STATES]
archive: dict[str, np.ndarray] = {
"schema_version": np.asarray(SCHEMA),
"structure_labels": np.asarray(STRUCTURE_LABELS),
"level_labels": np.asarray(LEVEL_LABELS),
}
for field in (
"state_id",
"element",
"electronic_model",
"rho_g_cc",
"temperature_ev",
"is_elapsed_s",
"sc_extension_elapsed_s",
"sc_converged",
"sc_iterations",
"fixed_is_mu_ha",
):
archive[field] = np.asarray([row[field] for row in ordered])
for field in (
"n_i_bohr3",
"r_ws_bohr",
"zstar",
"zbar_partition",
"ion_gamma_ha",
"stage2_error",
"stage2_iters",
"threshold_state_status",
):
archive[field] = np.asarray(
[[row[structure][field] for structure in STRUCTURE_LABELS]
for row in ordered]
)
for field in ("level_energy_ha", "level_m"):
archive[field] = np.asarray(
[[row[structure][field] for structure in STRUCTURE_LABELS]
for row in ordered],
dtype=float,
)
return add_archive_metadata(archive)
def add_archive_metadata(
arrays: dict[str, np.ndarray],
) -> dict[str, np.ndarray]:
"""Attach the compact provenance block required of accepted baselines."""
archive = {
key: np.asarray(value)
for key, value in arrays.items()
if key != "metadata_json"
}
state_count = int(np.asarray(archive["state_id"]).size)
metadata = {
"schema_version": "otter_compact_archive_metadata_v1",
"configuration": {
"scope": "is_and_experimental_sc_feedback",
"aa_n_points": 4096,
"bound_energy_cut_mode": "zero",
"bound_zero_tail_refine": True,
"bound_zero_tail_max_binding_ha": 3.0e-2,
"ion_gamma_mode": "scattering",
"al_electronic_model": "qm",
"fe_electronic_model": "tf",
"fe_mu_bounds_ha": [-2000.0, 200.0],
"sc_max_outer": int(SC_CONTROLS.max_outer),
"sc_g_tol": float(SC_CONTROLS.g_tol),
"sc_v_corr_tol_ha": float(SC_CONTROLS.v_corr_tol),
"sc_v_corr_mix": float(SC_CONTROLS.v_corr_mix),
"sc_fixed_is_mu": True,
},
"state": {"benchmark_id": BENCHMARK_ID, "count": state_count},
"producer": {
"project": "Otter",
"script_relative_path": (
"benchmarks/examples/"
"plot_starrett_saumon_2013_electronic.py"
),
},
"citation_keys": ["StarrettSaumon2013", "StarrettSaumon2014"],
"convergence": {
"accepted_state_pairs": state_count,
"failed_state_pairs": 0,
"is_stage2_converged": True,
"sc_feedback_converged": True,
},
"fields": sorted(archive),
"units": {
"rho_g_cc": "g cm^-3",
"temperature_ev": "eV",
"n_i_bohr3": "Bohr^-3",
"r_ws_bohr": "Bohr",
"ion_gamma_ha": "Hartree",
"level_energy_ha": "Hartree",
},
}
archive["metadata_json"] = np.asarray(
json.dumps(metadata, sort_keys=True, separators=(",", ":"))
)
return archive
def git_head() -> str:
"""Return the producing revision when Git metadata is available."""
completed = subprocess.run(
["git", "rev-parse", "HEAD"],
cwd=ROOT,
check=False,
capture_output=True,
text=True,
)
return completed.stdout.strip() if completed.returncode == 0 else "unknown"
def compute_candidate() -> dict[str, np.ndarray]:
"""Compute all IS/SC pairs in bounded parallelism and stage a candidate."""
rows: list[dict[str, Any]] = []
failures: dict[str, str] = {}
with ProcessPoolExecutor(max_workers=MAX_STATE_WORKERS) as pool:
futures = {pool.submit(solve_state, state): state for state in STATES}
for future in as_completed(futures):
state = futures[future]
state_id = str(state["state_id"])
try:
row = future.result()
except Exception as error:
failures[state_id] = f"{type(error).__name__}: {error}"
print(f"[rejected] {state_id}: {failures[state_id]}")
else:
rows.append(row)
total = row["is_elapsed_s"] + row["sc_extension_elapsed_s"]
print(
f"[accepted] {state_id}: {total:.1f} s, "
f"SC outer={row['sc_iterations']}"
)
if failures:
raise RuntimeError(f"IS/SC benchmark failures: {failures}")
archive = pack_rows(rows)
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
np.savez_compressed(CANDIDATE_PATH, **archive)
manifest = {
"schema_version": "otter_benchmark_manifest_v1",
"benchmark_id": BENCHMARK_ID,
"status": "candidate_not_accepted",
"configuration": {
"scope": "is_and_experimental_sc_feedback",
"bound_energy_cut_mode": "zero",
"bound_zero_tail_refine": True,
"bound_zero_tail_max_binding_ha": 3.0e-2,
"al_electronic_model": "qm",
"fe_electronic_model": "tf",
"fe_mu_bounds_ha": [-2000.0, 200.0],
"aa_n_points": 4096,
"sc_controls": {
"max_outer": int(SC_CONTROLS.max_outer),
"g_tol": float(SC_CONTROLS.g_tol),
"v_corr_tol_ha": float(SC_CONTROLS.v_corr_tol),
"v_corr_mix": float(SC_CONTROLS.v_corr_mix),
"fixed_is_mu": True,
},
},
"producer": {
"script_relative_path": str(Path(__file__).resolve().relative_to(ROOT)),
"script_sha256_current": sha256_file(Path(__file__).resolve()),
"git_commit": git_head(),
},
"state": {
"data_file": CANDIDATE_PATH.name,
"data_sha256": sha256_file(CANDIDATE_PATH),
"count": len(STATES),
},
}
CANDIDATE_MANIFEST.write_text(
json.dumps(manifest, indent=2, sort_keys=True) + "\n",
encoding="utf-8",
)
return archive
def validate_archive(state: dict[str, np.ndarray]) -> None:
"""Validate state order, IS/SC dimensions, and accepted diagnostics."""
if str(state["schema_version"].item()) != SCHEMA:
raise ValueError("Unsupported electronic benchmark schema.")
expected_ids = np.asarray([str(item["state_id"]) for item in STATES])
if not np.array_equal(np.asarray(state["state_id"], dtype=str), expected_ids):
raise ValueError("Electronic benchmark state order changed.")
if not np.array_equal(
np.asarray(state["structure_labels"], dtype=str), STRUCTURE_LABELS
):
raise ValueError("Electronic benchmark structure order changed.")
if np.asarray(state["level_energy_ha"]).shape != (len(STATES), 2, 4):
raise ValueError("Malformed bound-level table.")
for field in ("n_i_bohr3", "r_ws_bohr", "zstar", "zbar_partition"):
values = np.asarray(state[field], dtype=float)
if values.shape != (len(STATES), 2) or not np.all(np.isfinite(values)):
raise ValueError(f"Malformed or non-finite accepted field: {field}.")
if not np.all(np.asarray(state["sc_converged"], dtype=bool)):
raise ValueError("The accepted archive contains unconverged SC feedback.")
def load_precomputed() -> dict[str, np.ndarray]:
"""Verify and load the accepted Otter IS/SC baseline."""
manifest = json.loads(BASELINE_MANIFEST.read_text(encoding="utf-8"))
if (
manifest.get("schema_version") != "otter_benchmark_manifest_v1"
or manifest.get("benchmark_id") != BENCHMARK_ID
or manifest.get("status") != "accepted"
):
raise ValueError("Unexpected electronic benchmark manifest.")
record = dict(manifest["state"])
if record.get("data_file") != BASELINE_PATH.name:
raise ValueError("Electronic benchmark manifest names the wrong file.")
if sha256_file(BASELINE_PATH) != str(record.get("data_sha256")):
raise RuntimeError("Electronic benchmark baseline checksum mismatch.")
with np.load(BASELINE_PATH, allow_pickle=False) as archive:
state = {key: np.asarray(archive[key]) for key in archive.files}
validate_archive(state)
return state
def state_index(state: dict[str, np.ndarray], state_id: str) -> int:
"""Return the row index of one stable state ID."""
matches = np.flatnonzero(np.asarray(state["state_id"], dtype=str) == state_id)
if matches.size != 1:
raise KeyError(state_id)
return int(matches[0])
def print_comparison(state: dict[str, np.ndarray], reference: dict[str, Any]) -> None:
"""Print the table values using the gallery's three-digit convention."""
def cell(
value: float,
width: int,
*,
missing: str = "--",
) -> str:
"""Format three significant digits and label absent bound states."""
if not np.isfinite(value):
return f"{missing:>{width}}"
if value == 0.0:
rendered = "0.00"
else:
exponent = int(np.floor(np.log10(abs(value))))
decimals = max(0, 2 - exponent)
rendered = f"{value:.{decimals}f}"
return f"{rendered:>{width}}"
print("\nAl bound levels: all energies E are in Hartree")
print(
f"{'T[eV]':>6} {'shell':>5} {'paper E[Ha]':>12} "
f"{'IS E[Ha]':>11} {'SC E[Ha]':>11} "
f"{'M paper':>9} {'M IS':>9} {'M SC':>9}"
)
for row in reference["levels"]:
temperature = float(row["temperature_ev"])
shell = str(row["shell"])
index = state_index(state, f"al_rho2p7_te{temperature:g}_qm")
column = LEVEL_LABELS.index(shell)
paper_energy = float(row["energy_ha"])
values = (
temperature,
shell,
paper_energy,
float(state["level_energy_ha"][index, 0, column]),
float(state["level_energy_ha"][index, 1, column]),
float(row["m_weight"]),
float(state["level_m"][index, 0, column]),
float(state["level_m"][index, 1, column]),
)
print(
f"{values[0]:6g} {values[1]:>5} "
f"{cell(values[2], 12, missing='unbound')} "
f"{cell(values[3], 11, missing='unbound')} "
f"{cell(values[4], 11, missing='unbound')} "
f"{cell(values[5], 9)} {cell(values[6], 9)} "
f"{cell(values[7], 9)}"
)
paper_3s_m = ion_level_weight(-0.0125, 0.174)
print(
"Eq. (81) formula audit at 15 eV 3s: "
f"paper inputs give M={paper_3s_m:.3f}; published M=0.134."
)
print("\nScattering width gamma: all values are in Hartree")
print(
f"{'T[eV]':>6} {'paper gamma[Ha]':>16} "
f"{'IS gamma[Ha]':>13} {'SC gamma[Ha]':>13}"
)
for temperature in (2.0, 15.0):
index = state_index(state, f"al_rho2p7_te{temperature:g}_qm")
paper_gamma = float(
next(
row["gamma_ha"]
for row in reference["levels"]
if float(row["temperature_ev"]) == temperature
)
)
print(
f"{temperature:6g} {cell(paper_gamma, 16)} "
f"{cell(float(state['ion_gamma_ha'][index, 0]), 13)} "
f"{cell(float(state['ion_gamma_ha'][index, 1]), 13)}"
)
print("\nMean ionization: all values are dimensionless")
print(
f"{'state':>24} {'Z* paper':>9} {'Z* IS':>9} {'Z* SC':>9} "
f"{'Zbar paper':>11} {'Zbar IS':>10} {'Zbar SC':>10}"
)
for element, table in (
("Al", reference["al_ionization"]),
("Fe", reference["fe_ionization"]),
):
for row in np.atleast_1d(table):
temperature = float(row["temperature_ev"])
state_id = (
f"al_rho2p7_te{temperature:g}_qm"
if element == "Al"
else (
f"fe_rho{float(row['rho_g_cc']):g}_"
f"te{temperature:g}_tf"
)
)
index = state_index(state, state_id)
print(
f"{state_id:>24} {cell(float(row['zstar']), 9)} "
f"{cell(float(state['zstar'][index, 0]), 9)} "
f"{cell(float(state['zstar'][index, 1]), 9)} "
f"{cell(float(row['zbar']), 11)} "
f"{cell(float(state['zbar_partition'][index, 0]), 10)} "
f"{cell(float(state['zbar_partition'][index, 1]), 10)}"
)
# Load the same validated state table once for all gallery cells below.
if __name__ == "__main__":
reference = load_reference_tables()
state = compute_candidate() if not USE_PRECOMPUTED_DATA else load_precomputed()
validate_archive(state)
print(
"Using "
+ (
"checksummed, accepted Otter IS/SC states."
if USE_PRECOMPUTED_DATA
else "new candidate Otter IS/SC states."
)
)
print_comparison(state, reference)