Note
Go to the end to download the full example code.
Starrett–Saumon electronic levels and ionization
This benchmark compares 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 covers bound energies, pressure-ionization weights and mean ionization.
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. The tables give Otter values to three significant digits and published values at the article’s precision. 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.03 | 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.51 | -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.00835 | 0.134 | — | 0.0916 |
| \(T\) [eV] | Paper SC | Otter-IS | Otter-SC |
|---|---|---|---|
| 2 | 0.0698 | 0.0378 | 0.0630 |
| 15 | 0.174 | 0.111 | 0.171 |
Evaluating 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. Saved numerical results retain full precision. Tables II and III use two 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 parameters from the paper.
| \(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.06 | 2.06 | 3.00 | 3.00 | 3.00 | 41.0 | 5.05 |
| 6 | 2.11 | 2.18 | 2.18 | 3.00 | 3.00 | 3.00 | 13.6 | 2.04 |
| 10 | 2.24 | 2.31 | 2.31 | 3.00 | 3.03 | 3.02 | 8.14 | 1.60 |
| 15 | 2.51 | 2.54 | 2.54 | 3.18 | 3.24 | 3.19 | 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 fixed-chemical-potential iteration differs from the simultaneous 2013 QTCP/TFTCP solver. Al uses quantum orbitals; Fe uses Thomas–Fermi electrons, matching the model named for Table III.
QM SC feedback uses zero-tail shallow-state matching by default, adaptive full-AA precision and the unmixed correlation-potential residual. Inner SCF and outer SC acceptance are distinct; the numerical output records both the outer residual and whether QM precision refinement was performed.
Reproduction
From the root of the complete Otter checkout, using Poetry, run:
poetry run python benchmarks/examples/plot_starrett_saumon_2013_electronic.py
Downloads are optional: .ipynb launches this repository script; .zip
contains both formats. See Reproducing examples and benchmarks for setup.
The script calculates the states from their input parameters and prints the numerical comparisons. No bundled Otter NPZ is required. Numerical outputs are written locally; literature reference tables remain inputs to the comparison.
Recorded results
The figures and output below are from the recorded validation run; running the source recalculates them with the installed Otter version.
Recorded terminal output
[IS start] al_rho2p7_te2_qm
[IS complete; SC start] al_rho2p7_te2_qm: 36.7 s
[IS start] al_rho2p7_te6_qm
[accepted] al_rho2p7_te2_qm: 253.0 s, SC outer=8
[IS complete; SC start] al_rho2p7_te6_qm: 35.0 s
[accepted] al_rho2p7_te6_qm: 198.4 s, SC outer=7
[IS start] al_rho2p7_te10_qm
[IS complete; SC start] al_rho2p7_te10_qm: 32.4 s
[accepted] al_rho2p7_te10_qm: 162.6 s, SC outer=7
[IS start] al_rho2p7_te15_qm
[IS complete; SC start] al_rho2p7_te15_qm: 27.7 s
[accepted] al_rho2p7_te15_qm: 321.4 s, SC outer=7
[IS start] fe_rho22.5_te10_tf
[IS complete; SC start] fe_rho22.5_te10_tf: 1.1 s
[accepted] fe_rho22.5_te10_tf: 10.9 s, SC outer=10
[IS start] fe_rho34.5_te100_tf
[IS complete; SC start] fe_rho34.5_te100_tf: 1.0 s
[accepted] fe_rho34.5_te100_tf: 10.0 s, SC outer=10
[IS start] fe_rho39.65_te1000_tf
[IS complete; SC start] fe_rho39.65_te1000_tf: 1.0 s
[IS start] fe_rho34.37_te5000_tf
[accepted] fe_rho39.65_te1000_tf: 12.9 s, SC outer=13
[IS complete; SC start] fe_rho34.37_te5000_tf: 1.0 s
[accepted] fe_rho34.37_te5000_tf: 12.5 s, SC outer=13
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.03 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.51 -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.00837 0.134 -- 0.0919
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.0379 0.0629
15 0.174 0.111 0.171
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.06 2.06 3.00 3.00 3.00
al_rho2p7_te6_qm 2.11 2.18 2.18 3.00 3.00 3.00
al_rho2p7_te10_qm 2.24 2.31 2.31 3.00 3.03 3.02
al_rho2p7_te15_qm 2.51 2.54 2.54 3.18 3.24 3.19
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
import os
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 = False
if os.environ.get("OTTER_RECOMPUTE_STARRETT_SAUMON_ELECTRONIC", "0") == "1":
USE_PRECOMPUTED_DATA = False
MAX_STATE_WORKERS = 1
# =============================================================================
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),
)
# Common IS/SC comparison protocol. The validated hot-Fe controls needed
# 13 outer iterations, beyond the library's default budget of 10.
SC_CONTROLS = SCFeedbackConfig(
max_outer=16,
v_corr_mix=0.5,
)
def repository_root() -> Path:
"""Locate source and reference inputs, independently of numerical outputs."""
candidates = [Path.cwd().resolve(), *Path.cwd().resolve().parents]
source_file = globals().get("__file__")
if source_file is not None:
candidates.extend(Path(source_file).resolve().parents)
for candidate in candidates:
if (candidate / "pyproject.toml").is_file() and (candidate / "src/otter").is_dir():
return candidate
raise FileNotFoundError("Run from an Otter source checkout with its dependencies installed.")
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 producer_source() -> dict[str, Any]:
"""Identify the actual script, without inventing a notebook source hash."""
source = globals().get("__file__")
path = Path(source).resolve() if source is not None else None
return {
"execution_mode": "script" if path is not None else "interactive",
"script_filename": path.name if path is not None else None,
"script_relative_path": path.relative_to(ROOT).as_posix()
if path is not None and path.is_relative_to(ROOT) else None,
"script_sha256_current": sha256_file(path) if path is not None else None,
}
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"
model_override = {"electronic_model": "tf"} if is_tf else {}
aa_options: dict[str, Any] = {}
if is_tf:
# The nondegenerate 5000-eV Fe chemical potential lies below the
# ordinary warm-dense bracket.
aa_options["mu_bounds"] = (-2000.0, 200.0)
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"]),
**model_override,
aa_overrides=aa_options,
hnc_closure_transform_tol=2.5e-3,
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)
print(f"[IS start] {state['state_id']}", flush=True)
started = time.perf_counter()
is_workflow = solve_plasma_workflow(config)
is_elapsed_s = time.perf_counter() - started
print(f"[IS complete; SC start] {state['state_id']}: {is_elapsed_s:.1f} s", flush=True)
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"]),
"sc_max_g_change": float(feedback["history"][-1]["max_g_change"]),
"sc_v_corr_residual_ha": float(feedback["history"][-1]["max_v_corr_residual_ha"]),
"sc_inner_full_refined": bool(feedback["history"][-1]["inner_full_refined"]),
}
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),
"storage_profile": np.asarray("electronic_summary"),
"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",
"sc_max_g_change",
"sc_v_corr_residual_ha",
"sc_inner_full_refined",
):
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",
"aa_input_overrides_by_state": {
state["state_id"]: workflow_config(state).aa_overrides for state in STATES
},
"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,
"sc_inner_full_tol_scale": SC_CONTROLS.inner_full_tol_scale,
"sc_potential_convergence_metric": "unmixed_current_output_minus_used_input",
},
"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] = {}
checkpoint_dir = OUTPUT_DIR / "completed_pairs"
checkpoint_dir.mkdir(parents=True, exist_ok=True)
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]}", flush=True)
else:
rows.append(row)
# Compact scalar/level rows, not a restartable AA density dump.
(checkpoint_dir / f"{state_id}.json").write_text(
json.dumps(row, indent=2, default=lambda value: np.asarray(value).tolist())
)
total = row["is_elapsed_s"] + row["sc_extension_elapsed_s"]
print(
f"[accepted] {state_id}: {total:.1f} s, "
f"SC outer={row['sc_iterations']}", flush=True,
)
(OUTPUT_DIR / "failures.json").write_text(json.dumps(failures, indent=2))
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",
"external_average_atom": True,
"qoz_hnc": True,
"bound_energy_cut_mode": "zero",
"aa_input_overrides_by_state": {
state["state_id"]: workflow_config(state).aa_overrides for state in STATES
},
"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,
"inner_full_tol_scale": SC_CONTROLS.inner_full_tol_scale,
"potential_convergence_metric": "unmixed_current_output_minus_used_input",
},
},
"producer": {
**producer_source(),
"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)}"
)
def main() -> None:
"""Calculate and print all level/ionization comparisons."""
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)
if __name__ == "__main__":
main()