Note
Go to the end to download the full example code.
Starrett et al. CH1.36 mixture pair distributions
This benchmark covers the nine thermodynamic states in Figure 3 of Starrett et al. [2014].
The open markers are independent digitizations of the solid IS-QM curves in the published figure. The continuous curves are the September 2026 Otter recalculation of all nine states. Their configuration, producing-source provenance, diagnostics, and checksums are recorded in the manifest.
The orbital average atom, Appendix-A continuum, Appendix-B density tail,
pseudoatom/QOZ construction, and IS approximation follow
Starrett and Saumon [2014]; the mixture construction follows
Starrett et al. [2014]; and the finite-temperature jellium LFC follows
Chabrier [1990]. See the provenance, numerical metrics, model
configuration, and redistribution boundary for interpretation.
The comparison figure is exported as matching PNG and vector PDF files under
benchmarks/outputs/starrett_et_al_2014_mixtures_fig3/figures.
Reproduction
From the root of the complete Otter checkout, using Poetry, run:
poetry run python benchmarks/examples/plot_starrett_et_al_2014_mixtures_fig3.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 then plots the results. 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
Otter 0.4.1
[run] C:1 H:1.36 | rho=2.94 g/cm^3 | Te=1.72347 eV | Ti=1.72347 eV | QM -> QOZ/HNC
[mixture] iter= 1 max|dmu|=5.237e-01 Ha
[mixture] iter= 2 max|dmu|=6.484e-02 Ha
[mixture] iter= 3 max|dmu|=5.776e-03 Ha
[mixture] iter= 4 max|dmu|=4.172e-05 Ha
[AA] C (Z=6) | Rws=2.010544 a0 | full + external
[AA] converged=True | mu=0.54479406 Ha | Zbar=3.99874882 | time=27.183 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -9.08662 | 1.00000
2s | -0.00168 | 0.99982
[AA] H (Z=1) | Rws=1.438639 a0 | full + external
[AA] converged=True | mu=0.54475234 Ha | Zbar=1.00000000 | time=5.982 s
Bound levels
none
[mixture] mu=0.54477320 Ha
[mixture] C: Rws=2.010544 Bohr Zbar=3.99874882
[mixture] H: Rws=1.438639 Bohr Zbar=1.00000000
[ion] HNC iterations=4 QOZ=0.835 s HNC=0.620 s
[done] total time=263.379 s
Otter 0.4.1
[run] C:1 H:1.36 | rho=2.94 g/cm^3 | Te=4.30867 eV | Ti=4.30867 eV | QM -> QOZ/HNC
[computed] CH1.36, rho=2.94 g/cc, T=20 kK
[mixture] iter= 1 max|dmu|=5.069e-01 Ha
[mixture] iter= 2 max|dmu|=7.799e-02 Ha
[mixture] iter= 3 max|dmu|=5.069e-01 Ha
[mixture] iter= 4 max|dmu|=7.069e-03 Ha
[mixture] iter= 5 max|dmu|=5.922e-05 Ha
[AA] C (Z=6) | Rws=2.015627 a0 | full + external
[AA] converged=True | mu=0.55542440 Ha | Zbar=3.99877300 | time=27.313 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -9.06911 | 1.00000
2s | -0.00170 | 0.97121
[AA] H (Z=1) | Rws=1.431282 a0 | full + external
[AA] converged=True | mu=0.55536517 Ha | Zbar=1.00000000 | time=5.832 s
Bound levels
none
[mixture] mu=0.55539479 Ha
[mixture] C: Rws=2.015627 Bohr Zbar=3.99877300
[mixture] H: Rws=1.431282 Bohr Zbar=1.00000000
[ion] HNC iterations=4 QOZ=0.800 s HNC=0.566 s
[done] total time=247.366 s
Otter 0.4.1
[run] C:1 H:1.36 | rho=2.94 g/cm^3 | Te=8.61733 eV | Ti=8.61733 eV | QM -> QOZ/HNC
[computed] CH1.36, rho=2.94 g/cc, T=50 kK
[mixture] iter= 1 max|dmu|=5.088e-01 Ha
[mixture] iter= 2 max|dmu|=1.214e-01 Ha
[mixture] iter= 3 max|dmu|=8.568e-03 Ha
[mixture] iter= 4 max|dmu|=1.214e-01 Ha
[mixture] iter= 5 max|dmu|=5.088e-01 Ha
[mixture] iter= 6 max|dmu|=6.799e-04 Ha
[mixture] iter= 7 max|dmu|=4.190e-07 Ha
[AA] C (Z=6) | Rws=2.027399 a0 | full + external
[AA] converged=True | mu=0.52026735 Ha | Zbar=3.99951139 | time=19.861 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -9.05268 | 1.00000
2s | -0.00096 | 0.83834
[AA] H (Z=1) | Rws=1.413803 a0 | full + external
[AA] converged=True | mu=0.52026693 Ha | Zbar=1.00000000 | time=3.601 s
Bound levels
none
[mixture] mu=0.52026714 Ha
[mixture] C: Rws=2.027399 Bohr Zbar=3.99951139
[mixture] H: Rws=1.413803 Bohr Zbar=1.00000000
[ion] HNC iterations=4 QOZ=0.798 s HNC=0.517 s
[done] total time=237.349 s
[computed] CH1.36, rho=2.94 g/cc, T=100 kK
Otter 0.4.1
[run] C:1 H:1.36 | rho=5 g/cm^3 | Te=1.72347 eV | Ti=1.72347 eV | QM -> QOZ/HNC
[mixture] iter= 1 max|dmu|=6.407e-01 Ha
[mixture] iter= 2 max|dmu|=2.241e-01 Ha
[mixture] iter= 3 max|dmu|=1.717e-02 Ha
[mixture] iter= 4 max|dmu|=2.167e-04 Ha
[mixture] iter= 5 max|dmu|=1.754e-07 Ha
[AA] C (Z=6) | Rws=1.710000 a0 | full + external
[AA] converged=True | mu=0.97249293 Ha | Zbar=4.00003289 | time=7.575 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -8.67277 | 1.00000
[AA] H (Z=1) | Rws=1.166669 a0 | full + external
[AA] converged=True | mu=0.97249275 Ha | Zbar=1.00000000 | time=5.550 s
Bound levels
none
[mixture] mu=0.97249284 Ha
[mixture] C: Rws=1.710000 Bohr Zbar=4.00003289
[mixture] H: Rws=1.166669 Bohr Zbar=1.00000000
[ion] HNC iterations=5 QOZ=0.796 s HNC=0.678 s
[done] total time=53.164 s
[computed] CH1.36, rho=5 g/cc, T=20 kK
Otter 0.4.1
[run] C:1 H:1.36 | rho=5 g/cm^3 | Te=4.30867 eV | Ti=4.30867 eV | QM -> QOZ/HNC
[mixture] iter= 1 max|dmu|=6.231e-01 Ha
[mixture] iter= 2 max|dmu|=2.358e-01 Ha
[mixture] iter= 3 max|dmu|=1.860e-02 Ha
[mixture] iter= 4 max|dmu|=2.372e-04 Ha
[mixture] iter= 5 max|dmu|=1.764e-07 Ha
[AA] C (Z=6) | Rws=1.712653 a0 | full + external
[AA] converged=True | mu=0.97935107 Ha | Zbar=4.00003244 | time=8.284 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -8.66256 | 1.00000
[AA] H (Z=1) | Rws=1.162457 a0 | full + external
[AA] converged=True | mu=0.97935089 Ha | Zbar=1.00000000 | time=5.527 s
Bound levels
none
[mixture] mu=0.97935098 Ha
[mixture] C: Rws=1.712653 Bohr Zbar=4.00003244
[mixture] H: Rws=1.162457 Bohr Zbar=1.00000000
[ion] HNC iterations=4 QOZ=0.816 s HNC=0.575 s
[done] total time=52.132 s
Otter 0.4.1
[run] C:1 H:1.36 | rho=5 g/cm^3 | Te=8.61733 eV | Ti=8.61733 eV | QM -> QOZ/HNC
[computed] CH1.36, rho=5 g/cc, T=50 kK
[mixture] iter= 1 max|dmu|=6.127e-01 Ha
[mixture] iter= 2 max|dmu|=2.664e-01 Ha
[mixture] iter= 3 max|dmu|=1.885e-02 Ha
[mixture] iter= 4 max|dmu|=2.487e-04 Ha
[mixture] iter= 5 max|dmu|=2.445e-07 Ha
[AA] C (Z=6) | Rws=1.716904 a0 | full + external
[AA] converged=True | mu=0.96529366 Ha | Zbar=4.00003176 | time=4.784 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -8.64210 | 1.00000
[AA] H (Z=1) | Rws=1.155615 a0 | full + external
[AA] converged=True | mu=0.96529341 Ha | Zbar=1.00000000 | time=3.197 s
Bound levels
none
[mixture] mu=0.96529354 Ha
[mixture] C: Rws=1.716904 Bohr Zbar=4.00003176
[mixture] H: Rws=1.155615 Bohr Zbar=1.00000000
[ion] HNC iterations=4 QOZ=0.833 s HNC=0.530 s
[done] total time=45.668 s
[computed] CH1.36, rho=5 g/cc, T=100 kK
Otter 0.4.1
[run] C:1 H:1.36 | rho=15 g/cm^3 | Te=1.72347 eV | Ti=1.72347 eV | QM -> QOZ/HNC
[mixture] iter= 1 max|dmu|=9.246e-01 Ha
[mixture] iter= 2 max|dmu|=9.387e-01 Ha
[mixture] iter= 3 max|dmu|=5.438e-02 Ha
[mixture] iter= 4 max|dmu|=3.833e-03 Ha
[mixture] iter= 5 max|dmu|=1.837e-06 Ha
[AA] C (Z=6) | Rws=1.210867 a0 | full + external
[AA] converged=True | mu=2.64153191 Ha | Zbar=4.00158793 | time=7.613 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -7.35136 | 1.00000
[AA] H (Z=1) | Rws=0.765994 a0 | full + external
[AA] converged=True | mu=2.64153007 Ha | Zbar=1.00000000 | time=4.173 s
Bound levels
none
[mixture] mu=2.64153099 Ha
[mixture] C: Rws=1.210867 Bohr Zbar=4.00158793
[mixture] H: Rws=0.765994 Bohr Zbar=1.00000000
[ion] HNC iterations=5 QOZ=0.836 s HNC=3.352 s
[done] total time=45.866 s
Otter 0.4.1
[run] C:1 H:1.36 | rho=15 g/cm^3 | Te=4.30867 eV | Ti=4.30867 eV | QM -> QOZ/HNC
[computed] CH1.36, rho=15 g/cc, T=20 kK
[mixture] iter= 1 max|dmu|=9.209e-01 Ha
[mixture] iter= 2 max|dmu|=9.418e-01 Ha
[mixture] iter= 3 max|dmu|=5.455e-02 Ha
[mixture] iter= 4 max|dmu|=3.842e-03 Ha
[mixture] iter= 5 max|dmu|=1.815e-06 Ha
[AA] C (Z=6) | Rws=1.211047 a0 | full + external
[AA] converged=True | mu=2.64057330 Ha | Zbar=4.00158713 | time=6.849 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -7.34888 | 1.00000
[AA] H (Z=1) | Rws=0.765665 a0 | full + external
[AA] converged=True | mu=2.64057148 Ha | Zbar=1.00000000 | time=3.719 s
Bound levels
none
[mixture] mu=2.64057239 Ha
[mixture] C: Rws=1.211047 Bohr Zbar=4.00158713
[mixture] H: Rws=0.765665 Bohr Zbar=1.00000000
[ion] HNC iterations=7 QOZ=0.827 s HNC=0.767 s
[done] total time=41.771 s
[computed] CH1.36, rho=15 g/cc, T=50 kK
Otter 0.4.1
[run] C:1 H:1.36 | rho=15 g/cm^3 | Te=8.61733 eV | Ti=8.61733 eV | QM -> QOZ/HNC
[mixture] iter= 1 max|dmu|=9.142e-01 Ha
[mixture] iter= 2 max|dmu|=9.515e-01 Ha
[mixture] iter= 3 max|dmu|=5.473e-02 Ha
[mixture] iter= 4 max|dmu|=3.829e-03 Ha
[mixture] iter= 5 max|dmu|=1.802e-06 Ha
[AA] C (Z=6) | Rws=1.211455 a0 | full + external
[AA] converged=True | mu=2.63319264 Ha | Zbar=4.00158666 | time=3.423 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -7.34126 | 1.00000
[AA] H (Z=1) | Rws=0.764913 a0 | full + external
[AA] converged=True | mu=2.63319084 Ha | Zbar=1.00000000 | time=2.109 s
Bound levels
none
[mixture] mu=2.63319174 Ha
[mixture] C: Rws=1.211455 Bohr Zbar=4.00158666
[mixture] H: Rws=0.764913 Bohr Zbar=1.00000000
[ion] HNC iterations=4 QOZ=0.841 s HNC=0.527 s
[done] total time=36.215 s
[computed] CH1.36, rho=15 g/cc, T=100 kK
rho T[kK] pair RMSE MAE max
2.94 20 CC 1.0918e-02 7.0505e-03 6.3662e-02
2.94 20 CH 4.4822e-03 3.4401e-03 1.4188e-02
2.94 20 HH 6.4847e-03 4.8891e-03 1.8988e-02
2.94 50 CC 6.8528e-03 5.2132e-03 1.7566e-02
2.94 50 CH 7.2180e-03 4.9852e-03 2.1639e-02
2.94 50 HH 2.1868e-02 1.2860e-02 7.1116e-02
2.94 100 CC 4.9615e-03 4.0318e-03 1.0064e-02
2.94 100 CH 4.3664e-03 3.6604e-03 9.1943e-03
2.94 100 HH 3.8817e-03 3.0670e-03 9.7004e-03
5.00 20 CC 1.1302e-02 7.3957e-03 3.5907e-02
5.00 20 CH 1.1349e-02 8.3758e-03 3.6837e-02
5.00 20 HH 6.7852e-03 5.2309e-03 1.5502e-02
5.00 50 CC 6.8752e-03 5.1426e-03 1.9235e-02
5.00 50 CH 6.5198e-03 5.1063e-03 1.5333e-02
5.00 50 HH 5.7391e-03 4.6329e-03 1.1426e-02
5.00 100 CC 5.4052e-03 3.9032e-03 1.7649e-02
5.00 100 CH 6.5439e-03 5.0662e-03 1.5357e-02
5.00 100 HH 9.1775e-03 6.4749e-03 2.6109e-02
15.00 20 CC 1.7835e-02 1.0420e-02 7.0617e-02
15.00 20 CH 1.4540e-02 9.3449e-03 5.0283e-02
15.00 20 HH 1.3569e-02 7.3008e-03 5.2509e-02
15.00 50 CC 1.1677e-02 6.8037e-03 5.3533e-02
15.00 50 CH 1.0628e-02 7.4013e-03 3.3336e-02
15.00 50 HH 1.0681e-02 6.3658e-03 4.5084e-02
15.00 100 CC 7.6125e-03 4.1540e-03 3.6401e-02
15.00 100 CH 6.7544e-03 4.8768e-03 2.2534e-02
15.00 100 HH 5.2174e-03 3.8024e-03 1.6015e-02
from __future__ import annotations
from concurrent.futures import ProcessPoolExecutor, as_completed
from contextlib import ExitStack
import hashlib
import json
import os
from pathlib import Path
import time
from typing import Any
# Avoid hidden BLAS/OpenMP oversubscription in the outer state pool.
os.environ.setdefault("OMP_NUM_THREADS", "1")
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
os.environ.setdefault("MKL_NUM_THREADS", "1")
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.lines import Line2D
from otter import PlasmaWorkflowConfig, solve_plasma_workflow
from otter.numerics.constants import KELVIN_TO_EV
from otter.plotting import PAIR_COLORS, grid_figsize, save_figure, style_context
# =============================================================================
# User input
# =============================================================================
USE_PRECOMPUTED_DATA = False
if os.environ.get("OTTER_RECOMPUTE_STARRETT_FIG3", "0") == "1":
USE_PRECOMPUTED_DATA = False
DENSITIES_G_CC = (2.94, 5.0, 15.0)
TEMPERATURES_KK = (20, 50, 100)
STOICHIOMETRIC_COUNTS = (1.0, 1.36)
# Serial states; each AA inherits the single-worker default.
MAX_STATE_WORKERS = 1
SPECIES_PARALLEL_JOBS = 1
MU_E_TOL_HA = 1.0e-4
HNC_TOL = 1.0e-5
HNC_CLOSURE_TOL = 1.0e-4
R_RETAIN_MAX_BOHR = 20.0
# =============================================================================
PAIR_ORDER = ("CC", "CH", "HH")
PAIR_COLUMNS = {"CC": (0, 1), "CH": (4, 5), "HH": (2, 3)}
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()
PRECOMPUTED_DIR = (
ROOT
/ "benchmarks"
/ "baselines"
/ "starrett_et_al_2014_mixtures_fig3"
)
OUTPUT_DIR = (
ROOT
/ "benchmarks"
/ "outputs"
/ "starrett_et_al_2014_mixtures_fig3"
/ "gallery_recomputed"
)
FIGURE_DIR = (
ROOT
/ "benchmarks"
/ "outputs"
/ "starrett_et_al_2014_mixtures_fig3"
/ "figures"
)
def sha256_file(path: Path) -> str:
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 accepted_manifest() -> dict[str, Any]:
manifest = json.loads(
(PRECOMPUTED_DIR / "manifest.json").read_text(encoding="utf-8")
)
if manifest.get("benchmark_id") != (
"starrett_et_al_2014_mixtures_fig3_ch1p36"
):
raise ValueError("Unexpected Starrett Figure 3 manifest.")
if len(manifest.get("states", ())) != 9:
raise ValueError("The Figure 3 manifest must contain nine states.")
return manifest
def state_key(record: dict[str, Any]) -> tuple[float, int]:
return float(record["rho_g_cc"]), int(record["temperature_kk"])
def load_result(path: Path) -> dict[str, np.ndarray]:
"""Load one pickle-free accepted or newly computed result."""
with np.load(path, allow_pickle=False) as archive:
result = {key: np.asarray(archive[key]) for key in archive.files}
if any(value.dtype.hasobject for value in result.values()):
raise TypeError(f"Object arrays are forbidden in {path}.")
schema = str(result["schema_version"].item())
if schema not in {
"otter_starrett_mixtures_fig3_baseline_v1",
"otter_gallery_starrett_fig3_v1",
}:
raise ValueError(f"Unsupported numerical-result schema in {path}.")
if tuple(str(value) for value in result["pair_labels"]) != PAIR_ORDER:
raise ValueError(f"Unexpected pair order in {path}.")
return result
def load_precomputed_results(
manifest: dict[str, Any],
) -> dict[tuple[float, int], dict[str, np.ndarray]]:
"""Verify every accepted result and digitization checksum."""
loaded: dict[tuple[float, int], dict[str, np.ndarray]] = {}
for record in manifest["states"]:
result_path = PRECOMPUTED_DIR / str(record["baseline_file"])
reference_path = (
PRECOMPUTED_DIR / str(record["reference_file"])
).resolve()
if sha256_file(result_path) != str(record["baseline_sha256"]):
raise RuntimeError(f"Checksum mismatch for {result_path}.")
if sha256_file(reference_path) != str(record["reference_sha256"]):
raise RuntimeError(f"Checksum mismatch for {reference_path}.")
loaded[state_key(record)] = load_result(result_path)
return loaded
def workflow_config(
rho_g_cc: float,
temperature_kk: int,
) -> PlasmaWorkflowConfig:
"""Build one strict public Otter C--H mixture calculation."""
temperature_ev = 1000.0 * float(temperature_kk) * KELVIN_TO_EV
return PlasmaWorkflowConfig(
elements=["C", "H"],
counts=list(STOICHIOMETRIC_COUNTS),
temperature_ev=temperature_ev,
ion_temperature_ev=temperature_ev,
rho_g_cc=float(rho_g_cc),
species_parallel_jobs=int(SPECIES_PARALLEL_JOBS),
hnc_tol=float(HNC_TOL),
)
def strict_mixture_check(
workflow: dict[str, Any],
) -> tuple[dict[str, Any], dict[str, Any]]:
"""Reject unconverged common-mu, AA, external-AA, or HNC output."""
if str(workflow["electronic"]["kind"]) != "mixture":
raise RuntimeError("The Figure 3 benchmark requires a mixture result.")
electronic = dict(workflow["electronic"]["result"])
meta = dict(electronic.get("meta", {}))
if not bool(meta.get("root_success", False)):
raise RuntimeError("The common-mu root did not converge.")
if float(meta.get("mu_residual_max_ha", np.inf)) > MU_E_TOL_HA:
raise RuntimeError("The common-mu residual exceeds tolerance.")
if not bool(meta.get("final_mu_root_success", False)):
raise RuntimeError("The full+external rerun lost common-mu closure.")
if float(meta.get("final_mu_residual_max_ha", np.inf)) > MU_E_TOL_HA:
raise RuntimeError("The final common-mu residual exceeds tolerance.")
species = [dict(entry) for entry in electronic["species"]]
if tuple(str(entry["element"]) for entry in species) != ("C", "H"):
raise RuntimeError("Unexpected species order in the mixture result.")
for entry in species:
symbol = str(entry["element"])
result = dict(entry["result"])
if result.get("stage2_converged") is not True:
raise RuntimeError(f"{symbol}: full AA stage 2 did not converge.")
if str(result.get("threshold_state_status", "")).lower() == "unresolved":
raise RuntimeError(f"{symbol}: unresolved threshold state.")
if dict(result.get("ext_status", {})).get("converged") is not True:
raise RuntimeError(f"{symbol}: external AA did not converge.")
ion = dict(workflow["ion"])
if ion.get("hnc_converged") is not True:
raise RuntimeError("The mixture HNC did not reach a physical root.")
if float(ion["hnc_output_residual"]) > HNC_TOL:
raise RuntimeError("The mixture HNC residual exceeds tolerance.")
if float(ion["closure_transform_max_abs"]) > HNC_CLOSURE_TOL:
raise RuntimeError("The mixture g/S transform-closure audit failed.")
return electronic, ion
def solve_state(
rho_g_cc: float,
temperature_kk: int,
) -> tuple[tuple[float, int], dict[str, np.ndarray]]:
"""Run one complete AA -> pseudoatom -> mixture QOZ/HNC state."""
started = time.perf_counter()
workflow = solve_plasma_workflow(
workflow_config(rho_g_cc, temperature_kk)
)
elapsed_s = time.perf_counter() - started
electronic, ion = strict_mixture_check(workflow)
r = np.asarray(ion["r"], dtype=float)
gij = np.asarray(ion["gij_r"], dtype=float)
if gij.ndim != 3 or gij.shape[:2] != (2, 2):
raise ValueError(f"Unexpected gij array shape: {gij.shape}.")
mask = r <= R_RETAIN_MAX_BOHR
species = [dict(entry) for entry in electronic["species"]]
payload = {
"schema_version": np.asarray("otter_gallery_starrett_fig3_v1"),
"storage_profile": np.asarray("benchmark_analysis"),
"rho_g_cc": np.asarray(float(rho_g_cc)),
"temperature_kk": np.asarray(int(temperature_kk)),
"temperature_ev": np.asarray(
1000.0 * float(temperature_kk) * KELVIN_TO_EV
),
"species_symbols": np.asarray(("C", "H")),
"species_counts": np.asarray(STOICHIOMETRIC_COUNTS),
"pair_labels": np.asarray(PAIR_ORDER),
"r_bohr": r[mask],
"g_ab": np.asarray(
(gij[0, 0, mask], gij[0, 1, mask], gij[1, 1, mask])
),
"zbar_partition": np.asarray(ion["zbar_partition"], dtype=float),
"mu_ha": np.asarray(
[float(dict(entry["result"])["mu"]) for entry in species]
),
"root_residual_ha": np.asarray(
float(dict(electronic["meta"])["final_mu_residual_max_ha"])
),
"hnc_output_residual": np.asarray(
float(ion["hnc_output_residual"])
),
"hnc_closure_mismatch": np.asarray(
float(ion["closure_transform_max_abs"])
),
"producer_elapsed_s": np.asarray(elapsed_s),
}
return (float(rho_g_cc), int(temperature_kk)), payload
def solve_all_states() -> dict[tuple[float, int], dict[str, np.ndarray]]:
"""Calculate all nine independent states with a bounded outer pool."""
jobs = tuple(
(float(rho), int(temperature))
for rho in DENSITIES_G_CC
for temperature in TEMPERATURES_KK
)
loaded: dict[tuple[float, int], dict[str, np.ndarray]] = {}
with ProcessPoolExecutor(max_workers=MAX_STATE_WORKERS) as pool:
futures = {
pool.submit(solve_state, rho, temperature): (rho, temperature)
for rho, temperature in jobs
}
for future in as_completed(futures):
key, payload = future.result()
loaded[key] = payload
print(
f"[computed] CH1.36, rho={key[0]:g} g/cc, "
f"T={key[1]} kK"
)
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
for (rho, temperature), payload in loaded.items():
rho_token = f"{rho:.2f}".replace(".", "p")
path = OUTPUT_DIR / (
f"CH1p36_rho{rho_token}gcc_T{temperature}kK_otter.npz"
)
np.savez_compressed(path, **payload)
print(f"[saved] {path}")
return loaded
def load_reference(
path: Path,
) -> dict[str, tuple[np.ndarray, np.ndarray]]:
"""Read the two-header, six-column Figure 3 digitization."""
data = np.atleast_2d(
np.asarray(np.genfromtxt(path, delimiter=",", skip_header=2), dtype=float)
)
if data.shape[1] != 6:
raise ValueError(f"Expected six columns in {path}, got {data.shape}.")
curves: dict[str, tuple[np.ndarray, np.ndarray]] = {}
for pair in PAIR_ORDER:
r_column, g_column = PAIR_COLUMNS[pair]
r = data[:, r_column]
g = data[:, g_column]
mask = np.isfinite(r) & np.isfinite(g) & (r >= 0.0)
r = r[mask]
g = g[mask]
order = np.argsort(r)
r = r[order]
g = g[order]
unique_r = np.unique(r)
curves[pair] = (
unique_r,
np.asarray([np.mean(g[r == value]) for value in unique_r]),
)
return curves
def main() -> None:
manifest = accepted_manifest() if USE_PRECOMPUTED_DATA else None
results = (
load_precomputed_results(manifest)
if USE_PRECOMPUTED_DATA
else solve_all_states()
)
references = {
(2.94, 20): load_reference(ROOT / 'benchmarks/reference_data/starrett_et_al_2014_mixtures_fig3/fig3_gab_rho2p94gcc_T20kK.csv'),
(2.94, 50): load_reference(ROOT / 'benchmarks/reference_data/starrett_et_al_2014_mixtures_fig3/fig3_gab_rho2p94gcc_T50kK.csv'),
(2.94, 100): load_reference(ROOT / 'benchmarks/reference_data/starrett_et_al_2014_mixtures_fig3/fig3_gab_rho2p94gcc_T100kK.csv'),
(5.0, 20): load_reference(ROOT / 'benchmarks/reference_data/starrett_et_al_2014_mixtures_fig3/fig3_gab_rho5p00gcc_T20kK.csv'),
(5.0, 50): load_reference(ROOT / 'benchmarks/reference_data/starrett_et_al_2014_mixtures_fig3/fig3_gab_rho5p00gcc_T50kK.csv'),
(5.0, 100): load_reference(ROOT / 'benchmarks/reference_data/starrett_et_al_2014_mixtures_fig3/fig3_gab_rho5p00gcc_T100kK.csv'),
(15.0, 20): load_reference(ROOT / 'benchmarks/reference_data/starrett_et_al_2014_mixtures_fig3/fig3_gab_rho15p00gcc_T20kK.csv'),
(15.0, 50): load_reference(ROOT / 'benchmarks/reference_data/starrett_et_al_2014_mixtures_fig3/fig3_gab_rho15p00gcc_T50kK.csv'),
(15.0, 100): load_reference(ROOT / 'benchmarks/reference_data/starrett_et_al_2014_mixtures_fig3/fig3_gab_rho15p00gcc_T100kK.csv'),
}
model_legend_label = "Otter"
print(
"Using "
+ (
"checksummed accepted arrays."
if USE_PRECOMPUTED_DATA
else "new results calculated directly by this gallery script."
)
)
print(
f"{'rho':>6s} {'T[kK]':>6s} {'pair':>4s} "
f"{'RMSE':>10s} {'MAE':>10s} {'max':>10s}"
)
for key in sorted(results):
result = results[key]
r_model = np.asarray(result["r_bohr"], dtype=float)
g_ab = np.asarray(result["g_ab"], dtype=float)
for pair_index, pair in enumerate(PAIR_ORDER):
r_ref, g_ref = references[key][pair]
mask = (
(r_ref >= r_model[0])
& (r_ref <= min(6.0, float(r_model[-1])))
)
delta = np.interp(r_ref[mask], r_model, g_ab[pair_index]) - g_ref[mask]
print(
f"{key[0]:6.2f} {key[1]:6d} {pair:>4s} "
f"{np.sqrt(np.mean(delta**2)):10.4e} "
f"{np.mean(np.abs(delta)):10.4e} "
f"{np.max(np.abs(delta)):10.4e}"
)
# %%
# Figure 3 overlay
# ----------------
#
# Each open marker is digitized from the published solid IS-QM curve. Each
# continuous line is the selected numerical result: either the checksum-
# verified archived result or a new calculation performed above.
plot_style = ExitStack()
plot_style.enter_context(style_context("thesis", palette="bing"))
fig, axes = plt.subplots(
3,
3,
figsize=grid_figsize(3, 3),
sharex=True,
sharey=True,
)
for row, density in enumerate(DENSITIES_G_CC):
for column, temperature in enumerate(TEMPERATURES_KK):
axis = axes[row, column]
key = (float(density), int(temperature))
result = results[key]
r_model = np.asarray(result["r_bohr"], dtype=float)
g_model = np.asarray(result["g_ab"], dtype=float)
for pair_index, pair in enumerate(PAIR_ORDER):
r_ref, g_ref = references[key][pair]
axis.plot(
r_ref,
g_ref,
marker="o",
markersize=6.0,
markerfacecolor="none",
markeredgewidth=1.0,
color=PAIR_COLORS[pair],
linestyle="none",
alpha=0.75,
)
axis.plot(
r_model,
g_model[pair_index],
color=PAIR_COLORS[pair],
lw=1.4,
)
axis.set(xlim=(-0.5, 6.0), ylim=(-0.05, 2.0))
if row == 0:
axis.set_title(f"{temperature} kK")
if column == 0:
axis.set_ylabel(
rf"$\rho={density:g}$ g cm$^{{-3}}$"
+ "\n"
+ r"$g_{ab}(r)$"
)
if row == 2:
axis.set_xlabel(r"$r$ [$a_{\rm B}$]")
handles = [
*(
Line2D([], [], color=PAIR_COLORS[pair], lw=2.0, label=pair)
for pair in PAIR_ORDER
),
Line2D(
[],
[],
color="0.35",
marker="o",
markersize=7.5,
markeredgewidth=1.2,
markerfacecolor="none",
linestyle="none",
label="Starrett Fig. 3 digitized",
),
Line2D([], [], color="0.35", lw=2.0, label=model_legend_label),
]
fig.legend(
handles=handles,
loc="upper center",
bbox_to_anchor=(0.5, 0.945),
ncol=5,
frameon=False,
)
fig.suptitle(r"CH$_{1.36}$: Starrett et al. (2014), Figure 3", y=0.985)
fig.text(
0.5,
0.006,
"Reference data: Starrett et al. (2014), "
"doi:10.1103/PhysRevE.90.033110.",
ha="center",
va="bottom",
fontsize=7.5,
)
fig.tight_layout(
rect=(0.0, 0.025, 1.0, 0.90),
pad=0.45,
w_pad=0.25,
h_pad=0.25,
)
save_figure(
fig,
FIGURE_DIR / "starrett_et_al_2014_mixtures_fig3",
close=False,
)
plot_style.close()
plt.show()
if __name__ == "__main__":
main()