Note
Go to the end to download the full example code.
Beryllium ionization: Döppner et al. (2023)
Beryllium ionization at \(T_e=50,100,150\) eV and densities 1–70 g cm-3 is compared with Fig. 3(a) of Döppner et al. [2023]. Otter uses the QM ion-sphere model with Dirac exchange. The two ionization definitions are
\[\bar Z_{\mathrm{AA,WS}}=4-4\pi\int_0^{R_{\rm WS}}n_{\rm ion}(r)r^2dr,
\qquad Z^*=n_e^0/n_i.\]
where \(n_{\rm ion}\) includes the pressure-ionization weight and radial cutoff of Starrett and Saumon [2014], Starrett and Saumon [2013]. Ionization definitions depend on the electron partition used by each model. The Stewart–Pyatt and OPAL reference curves in the 150 eV panel are at 160 eV, as labelled in the paper.
Reproduction
From the root of the complete Otter checkout, using Poetry, run:
poetry run python benchmarks/examples/plot_doppner_2023_be_ionization.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
[1/48] T=50 eV rho=1 g/cc Zbar=2.964777 Zstar=2.555507 resolved 23.3 s
[2/48] T=50 eV rho=3 g/cc Zbar=2.723907 Zstar=2.373393 resolved 17.5 s
[3/48] T=50 eV rho=6 g/cc Zbar=2.563412 Zstar=2.328884 resolved 12.4 s
[4/48] T=50 eV rho=10 g/cc Zbar=2.452816 Zstar=2.348229 resolved 11.2 s
[5/48] T=50 eV rho=15 g/cc Zbar=2.384223 Zstar=2.402004 resolved 9.9 s
[6/48] T=50 eV rho=20 g/cc Zbar=2.358773 Zstar=2.461425 resolved 8.8 s
[7/48] T=50 eV rho=25 g/cc Zbar=2.361489 Zstar=2.518648 resolved 8.6 s
[8/48] T=50 eV rho=30 g/cc Zbar=2.388167 Zstar=2.571559 resolved 8.0 s
[9/48] T=50 eV rho=35 g/cc Zbar=2.450928 Zstar=2.619794 resolved 8.3 s
[10/48] T=50 eV rho=40 g/cc Zbar=2.564377 Zstar=2.663580 resolved 8.1 s
[11/48] T=50 eV rho=45 g/cc Zbar=2.723122 Zstar=2.703326 resolved 7.5 s
[12/48] T=50 eV rho=50 g/cc Zbar=2.906421 Zstar=2.739482 resolved 7.4 s
[13/48] T=50 eV rho=55 g/cc Zbar=3.092768 Zstar=2.772465 resolved 7.5 s
[14/48] T=50 eV rho=60 g/cc Zbar=3.267179 Zstar=2.802658 resolved 7.5 s
[15/48] T=50 eV rho=65 g/cc Zbar=3.421701 Zstar=2.830396 resolved 7.6 s
[16/48] T=50 eV rho=70 g/cc Zbar=3.553468 Zstar=2.855966 resolved 8.2 s
[17/48] T=100 eV rho=1 g/cc Zbar=3.639929 Zstar=3.370990 resolved 69.1 s
[18/48] T=100 eV rho=3 g/cc Zbar=3.463554 Zstar=3.112195 resolved 25.5 s
[19/48] T=100 eV rho=6 g/cc Zbar=3.272257 Zstar=2.979115 resolved 17.1 s
[20/48] T=100 eV rho=10 g/cc Zbar=3.113323 Zstar=2.902815 resolved 12.8 s
[21/48] T=100 eV rho=15 g/cc Zbar=2.983562 Zstar=2.862399 resolved 11.4 s
[22/48] T=100 eV rho=20 g/cc Zbar=2.897873 Zstar=2.847469 resolved 10.1 s
[23/48] T=100 eV rho=25 g/cc Zbar=2.842111 Zstar=2.844952 resolved 9.9 s
[24/48] T=100 eV rho=30 g/cc Zbar=2.808916 Zstar=2.849205 resolved 9.6 s
[25/48] T=100 eV rho=35 g/cc Zbar=2.793742 Zstar=2.857388 resolved 9.3 s
[26/48] T=100 eV rho=40 g/cc Zbar=2.794529 Zstar=2.867909 resolved 8.9 s
[27/48] T=100 eV rho=45 g/cc Zbar=2.816337 Zstar=2.879817 resolved 8.8 s
[28/48] T=100 eV rho=50 g/cc Zbar=2.871971 Zstar=2.892497 resolved 8.5 s
[29/48] T=100 eV rho=55 g/cc Zbar=2.968219 Zstar=2.905580 resolved 9.3 s
[30/48] T=100 eV rho=60 g/cc Zbar=3.097090 Zstar=2.918809 resolved 8.7 s
[31/48] T=100 eV rho=65 g/cc Zbar=3.241862 Zstar=2.932012 resolved 8.3 s
[32/48] T=100 eV rho=70 g/cc Zbar=3.386276 Zstar=2.945076 resolved 9.0 s
[33/48] T=150 eV rho=1 g/cc Zbar=3.732699 Zstar=3.616604 resolved 54.6 s
[34/48] T=150 eV rho=3 g/cc Zbar=3.734954 Zstar=3.441902 resolved 34.2 s
[35/48] T=150 eV rho=6 g/cc Zbar=3.591192 Zstar=3.319597 resolved 24.5 s
[36/48] T=150 eV rho=10 g/cc Zbar=3.455964 Zstar=3.233148 resolved 17.7 s
[37/48] T=150 eV rho=15 g/cc Zbar=3.333723 Zstar=3.172562 resolved 15.4 s
[38/48] T=150 eV rho=20 g/cc Zbar=3.243999 Zstar=3.136608 resolved 12.6 s
[39/48] T=150 eV rho=25 g/cc Zbar=3.177586 Zstar=3.113996 resolved 11.8 s
[40/48] T=150 eV rho=30 g/cc Zbar=3.129483 Zstar=3.099566 resolved 11.2 s
[41/48] T=150 eV rho=35 g/cc Zbar=3.096499 Zstar=3.090544 resolved 10.4 s
[42/48] T=150 eV rho=40 g/cc Zbar=3.076381 Zstar=3.085258 resolved 9.8 s
[43/48] T=150 eV rho=45 g/cc Zbar=3.067486 Zstar=3.082658 resolved 9.5 s
[44/48] T=150 eV rho=50 g/cc Zbar=3.069509 Zstar=3.082035 resolved 6.9 s
[45/48] T=150 eV rho=55 g/cc Zbar=3.087031 Zstar=3.082886 resolved 5.8 s
[46/48] T=150 eV rho=60 g/cc Zbar=3.130020 Zstar=3.084858 resolved 5.9 s
[47/48] T=150 eV rho=65 g/cc Zbar=3.204567 Zstar=3.087689 resolved 6.3 s
[48/48] T=150 eV rho=70 g/cc Zbar=3.305774 Zstar=3.091180 resolved 5.7 s
Saved 48/48 converged full SCFs (candidate): benchmarks/outputs/doppner_2023_be_ionization/recomputed/Be_ionization.npz
Finished: 48/48 full SCFs; 640.6 s
QM ion-sphere, Dirac exchange.
Full SCF: 48/48; threshold: resolved=48; failed=0
T [eV] rho [g/cc] Zbar AA,WS Zstar mu [Ha] threshold time [s]
50 1 2.964777 2.555507 -4.588931 resolved 23.32
50 3 2.723907 2.373393 -2.613666 resolved 17.47
50 6 2.563412 2.328884 -1.236708 resolved 12.36
50 10 2.452816 2.348229 -0.092580 resolved 11.20
50 15 2.384223 2.402004 0.942865 resolved 9.93
50 20 2.358773 2.461425 1.775348 resolved 8.78
50 25 2.361489 2.518648 2.494270 resolved 8.58
50 30 2.388167 2.571559 3.139782 resolved 7.98
50 35 2.450928 2.619794 3.733472 resolved 8.29
50 40 2.564377 2.663580 4.288327 resolved 8.13
50 45 2.723122 2.703326 4.812779 resolved 7.49
50 50 2.906421 2.739482 5.312640 resolved 7.44
50 55 3.092768 2.772465 5.792078 resolved 7.53
50 60 3.267179 2.802658 6.254206 resolved 7.52
50 65 3.421701 2.830396 6.701409 resolved 7.60
50 70 3.553468 2.855966 7.135558 resolved 8.18
100 1 3.639929 3.370990 -12.036350 resolved 69.08
100 3 3.463554 3.112195 -8.206947 resolved 25.49
100 6 3.272257 2.979115 -5.698006 resolved 17.08
100 10 3.113323 2.902815 -3.756857 resolved 12.80
100 15 2.983562 2.862399 -2.120504 resolved 11.38
100 20 2.897873 2.847469 -0.884025 resolved 10.14
100 25 2.842111 2.844952 0.132687 resolved 9.89
100 30 2.808916 2.849205 1.009780 resolved 9.60
100 35 2.793742 2.857388 1.790066 resolved 9.26
100 40 2.794529 2.867909 2.499128 resolved 8.94
100 45 2.816337 2.879817 3.153493 resolved 8.76
100 50 2.871971 2.892497 3.764428 resolved 8.53
100 55 2.968219 2.905580 4.340039 resolved 9.28
100 60 3.097090 2.918809 4.886285 resolved 8.66
100 65 3.241862 2.932012 5.407710 resolved 8.34
100 70 3.386276 2.945076 5.907847 resolved 8.97
150 1 3.732699 3.616604 -21.049697 resolved 54.59
150 3 3.734954 3.441902 -15.187947 resolved 34.23
150 6 3.591192 3.319597 -11.454144 resolved 24.53
150 10 3.455964 3.233148 -8.638597 resolved 17.71
150 15 3.333723 3.172562 -6.329835 resolved 15.43
150 20 3.243999 3.136608 -4.630653 resolved 12.62
150 25 3.177586 3.113996 -3.265030 resolved 11.79
150 30 3.129483 3.099566 -2.110366 resolved 11.17
150 35 3.096499 3.090544 -1.101369 resolved 10.40
150 40 3.076381 3.085258 -0.199150 resolved 9.84
150 45 3.067486 3.082658 0.621396 resolved 9.53
150 50 3.069509 3.082035 1.377402 resolved 6.92
150 55 3.087031 3.082886 2.081063 resolved 5.84
150 60 3.130020 3.084858 2.741413 resolved 5.88
150 65 3.204567 3.087689 3.365306 resolved 6.30
150 70 3.305774 3.091180 3.958068 resolved 5.70
from __future__ import annotations
import csv
import hashlib
import json
import os
from pathlib import Path
import subprocess
import sys
import matplotlib.pyplot as plt
import numpy as np
from otter.plotting import PALETTES, grid_figsize, save_figure, set_style
RECOMPUTE_WITH_OTTER = True
USE_CANDIDATES = True
if os.environ.get("OTTER_USE_CANDIDATE_BE_IONIZATION", "0") == "1":
USE_CANDIDATES = True
PACKAGE = "doppner_2023_be_ionization"
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()
def load_reference() -> list[dict]:
"""Verify and parse ragged two-column series, never treating blanks as zero."""
directory = ROOT / "benchmarks/reference_data/doppner_2023_Be_ionization_fig3a"
manifest = json.loads((directory / "manifest.json").read_text())
path = directory / manifest["file"]
if hashlib.sha256(path.read_bytes()).hexdigest() != manifest["sha256"]:
raise ValueError("Be reference checksum mismatch")
with path.open(newline="") as stream:
rows = list(csv.reader(stream))
curves = []
for series in manifest["series"]:
i, j = series["columns"]
if rows[0][i] != series["header"] or [rows[1][i], rows[1][j]] != ["rho", "Zbar"]:
raise ValueError("Be reference column mapping changed")
values = []
for row in rows[2:]:
if bool(row[i]) != bool(row[j]):
raise ValueError("Incomplete reference coordinate pair")
if row[i]:
values.append((float(row[i]), float(row[j])))
array = np.asarray(values)
if (len(array) != series["points"] or not np.all(np.isfinite(array))
or not np.all(np.diff(array[:, 0]) > 0)):
raise ValueError("Invalid Be reference series")
curves.append({**series, "rho": array[:, 0], "z": array[:, 1]})
return curves
def load_scan() -> dict[str, np.ndarray]:
directory = ROOT / "benchmarks" / ("outputs" if USE_CANDIDATES else "baselines") / PACKAGE
if USE_CANDIDATES:
directory /= "recomputed"
manifest = json.loads((directory / "manifest.json").read_text())
record = manifest["states"][0]
path = directory / record["baseline_file"]
if hashlib.sha256(path.read_bytes()).hexdigest() != record["baseline_sha256"]:
raise ValueError("Be Otter checksum mismatch")
with np.load(path, allow_pickle=False) as archive:
data = {k: archive[k] for k in archive.files}
audit = manifest["scientific_audit"]
if (data["zbar"].size != audit["full_scf_converged_states"]
or not np.all(data["stage2_converged"])):
raise ValueError("Be scan disagrees with its convergence audit")
return data
def plot_comparison(data: dict[str, np.ndarray], references: list[dict]):
"""Use the Bethkenhagen gallery's Otter/reference visual convention."""
set_style("thesis", palette="bing")
colors = PALETTES["bing"]
reference_styles = {
"DFT-MD": (colors[2], ":", "o"),
"Stewart-Pyatt": (colors[3], "-.", None),
"OPAL": (colors[4], "--", None),
}
fig, axes = plt.subplots(1, 3, figsize=grid_figsize(1, 3), sharey=True)
for ax, temperature in zip(axes, (50, 100, 150), strict=True):
mask = data["te_ev"] == temperature
rho = np.asarray([1, 3, 6, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70])
indices = np.searchsorted(rho, data["rho_g_cc"][mask])
resolved = data["threshold_state_status"][mask] == "resolved"
for field, label, color in (
("zbar", r"Otter $\bar Z_{\rm AA,WS}$", "#131313"),
("zstar", r"Otter $Z^*$", colors[1]),
):
y = data[field][mask]
plotted = np.full(rho.size, np.nan)
plotted[indices] = np.where(resolved, y, np.nan)
ax.plot(rho, plotted, color=color, ls="-", lw=2.3,
alpha=0.45, marker="o", ms=4.0, label=label)
if np.any(~resolved):
ax.scatter(rho[indices[~resolved]], y[~resolved], color=color, marker="x", zorder=5)
for reference in references:
t_ref = reference["temperature_ev"]
if t_ref != temperature and not (temperature == 150 and t_ref == 160):
continue
method = reference["method"]
label = method + (" (160 eV)" if t_ref == 160 else "")
color, line_style, marker = reference_styles[method]
ax.plot(reference["rho"], reference["z"],
color=color, ls=line_style, lw=1.75,
marker=marker, ms=5.0 if marker else None,
markerfacecolor="white" if marker else None,
markeredgewidth=1.1 if marker else None,
alpha=0.75, label=label)
ax.set(title=rf"Be: $T_e={temperature}$ eV", xlabel=r"$\rho$ (g cm$^{-3}$)",
xlim=(0, 72))
ax.legend(ncol=1 if temperature == 150 else 2, fontsize="small",
loc="best", columnspacing=0.8, handlelength=1.8)
axes[0].set_ylabel("Ionization (electrons per Be)")
fig.tight_layout()
return fig
def main() -> None:
if RECOMPUTE_WITH_OTTER:
subprocess.run([sys.executable, str(ROOT / "benchmarks/runners/regenerate_doppner_2023_be_ionization.py")], check=True)
data = load_scan()
references = load_reference()
fig = plot_comparison(data, references)
save_figure(fig, ROOT / "benchmarks/outputs" / PACKAGE / "figures/be_ionization")
print("QM ion-sphere, Dirac exchange.")
audit = json.loads(data["metadata_json"].item())["convergence"]
counts = ", ".join(f"{key}={value}" for key, value in sorted(audit["threshold_counts"].items()))
print(f"Full SCF: {audit['full_scf_converged_states']}/{audit['requested_states']}; "
f"threshold: {counts}; failed={audit['stage2_nonconverged_states']}")
if np.any(data["threshold_state_status"] != "resolved"):
print("Crosses mark non-resolved threshold states; failed points are omitted.")
print(f"{'T [eV]':>8} {'rho [g/cc]':>12} {'Zbar AA,WS':>12} {'Zstar':>10} {'mu [Ha]':>12} {'threshold':>12} {'time [s]':>10}")
for i in range(data["zbar"].size):
print(f"{data['te_ev'][i]:8g} {data['rho_g_cc'][i]:12g} {data['zbar'][i]:12.6f} "
f"{data['zstar'][i]:10.6f} {data['mu_ha'][i]:12.6f} "
f"{data['threshold_state_status'][i]:>12} {data['elapsed_s'][i]:10.2f}")
for failure in audit.get("failures", []):
print(f"FAILED: T={failure['te_ev']} eV, rho={failure['rho_g_cc']} g/cc: {failure['error']}")
if __name__ == "__main__":
main()
Reproduction
The runner called above computes all 48 states, one AA at a time. Results are saved in
benchmarks/outputs/doppner_2023_be_ionization/recomputed.