.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "benchmarks/gen_benchmarks/plot_starrett_saumon_2013_electronic.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_benchmarks_gen_benchmarks_plot_starrett_saumon_2013_electronic.py: 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 :cite:t:`StarrettSaumon2013` (doi: `10.1103/PhysRevE.87.013104 `__). It was added specifically to make pressure ionization auditable before interpreting the :math:`\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 :math:`E=-0.0125` Hartree (about -0.340 eV), not -0.0125 eV. For a bound level of energy :math:`E<0`, Otter uses Eq. (81), .. math:: M(E)=\operatorname{erf}\!\left[-2\sqrt{\ln 2}\,E/\gamma\right], where :math:`\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 :math:`E` and :math:`\gamma` entry uses Hartree, while :math:`M(E)` and both ionization definitions are dimensionless. .. raw:: html
Al electronic diagnostics at 2.7 g cm-3 and \(T=2\,\mathrm{eV}\)
State\(E\) [Ha] \(M(E)\)
Paper SCOtter-ISOtter-SC Paper SCOtter-ISOtter-SC
1s-54.6-54.6-54.6 1.001.001.00
2s-3.41-3.39-3.40 1.001.001.00
2p-2.04-2.02-2.02 1.001.001.00
3sunboundunboundunbound ———
Al electronic diagnostics at 2.7 g cm-3 and \(T=15\,\mathrm{eV}\)
State\(E\) [Ha] \(M(E)\)
Paper SCOtter-ISOtter-SC Paper SCOtter-ISOtter-SC
1s-54.9-54.8-54.8 1.001.001.00
2s-3.60-3.50-3.57 1.001.001.00
2p-2.23-2.14-2.21 1.001.001.00
3s-0.0125unbound-0.00763 0.134—0.0842
Scattering width \(\gamma\) [Ha]
\(T\) [eV]Paper SC Otter-ISOtter-SC
20.06980.03670.0617
150.1740.1090.170
As a formula-only audit, evaluating Otter's Eq. (81) with the *published* 15-eV 3s inputs :math:`E=-0.0125` Ha and :math:`\gamma=0.174` Ha gives :math:`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: .. math:: 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 :math:`\Gamma_{\rm OCP}` and :math:`\Gamma_{\rm TCP}` columns are transcribed paper parameters, not Otter observables; this benchmark does not manufacture nominally corresponding Otter values. .. raw:: html
Al at 2.7 g cm-3 (\(n_I^0=8.93\times10^{-3}\,\mathrm{Bohr}^{-3}\)): Table II versus Otter
\(T\) [eV] \(Z^*\) \(\bar Z\) Paper coupling
PaperOtter-ISOtter-SC PaperOtter-ISOtter-SC \(\Gamma_{\rm OCP}\)\(\Gamma_{\rm TCP}\)
21.982.092.09 3.003.003.0041.05.05
62.112.192.19 3.003.003.0013.62.04
102.242.332.33 3.003.033.028.141.60
152.512.562.56 3.183.243.206.121.34
Fe principal Hugoniot: Table III versus Otter
\(\rho\) [g cm\(^{-3}\)] \(n_I^0\) [Bohr\(^{-3}\)] \(T\) [eV] \(Z^*\) \(\bar Z\) Paper coupling
PaperOtter-ISOtter-SC PaperOtter-ISOtter-SC \(\Gamma_{\rm OCP}\)\(\Gamma_{\rm TCP}\)
22.50.036010 5.855.855.85 8.788.818.7411214.3
34.50.0551100 9.549.549.54 11.611.911.522.44.85
39.650.06341000 20.420.420.4 21.721.921.68.223.32
34.370.05495000 25.125.125.1 25.525.625.62.181.35
Otter-IS is the converged production ion-sphere workflow. Otter-SC is the experimental feedback path based on Sec. 2.4 of :cite:t:`StarrettSaumon2014`: 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. .. GENERATED FROM PYTHON SOURCE LINES 163-795 .. rst-class:: sphx-glr-script-out .. code-block:: none 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 | .. code-block:: Python 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) .. _sphx_glr_download_benchmarks_gen_benchmarks_plot_starrett_saumon_2013_electronic.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_starrett_saumon_2013_electronic.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_starrett_saumon_2013_electronic.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_starrett_saumon_2013_electronic.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_