Note
Go to the end to download the full example code.
Carbon structure factors: Otter PA-HNC and DFT-MD
This benchmark calculates carbon static ion structure factors with Otter’s average-atom -> pseudoatom -> QOZ/HNC workflow and compares them with unpublished DFT-MD data provided by Dr. Argha Roy (private communication).
The DFT-MD files contain wave number in inverse ångström, Sii(k), and a
reported uncertainty. The uncertainty is not called one standard deviation
because that interpretation has not been confirmed. Redistribution
permission is not recorded, so the numerical files carry attribution and
license status NOASSERTION.
The electronic and ionic construction follows Starrett and Saumon [2014], Starrett and Saumon [2013]; the default finite-temperature jellium local-field correction follows Chabrier [1990].
Reproduction
From the root of the complete Otter checkout, using Poetry, run:
poetry run python benchmarks/examples/plot_argha_roy_carbon_sii.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 | rho=3.51538 g/cm^3 | Te=20 eV | Ti=20 eV | QM -> QOZ/HNC
[AA] C (Z=6) | Rws=2.090843 a0 | full + external
[AA/full SCF: stage 1]
iter | d_n | d_v
0 | nan | 4.878e-02
1 | 3.217e-02 | 4.456e-02
2 | 2.855e-02 | 4.168e-02
3 | 2.632e-02 | 4.026e-02
4 | 2.614e-02 | 4.024e-02
[AA/full SCF: stage 2]
iter | d_n | d_v
0 | nan | 4.135e-02
1 | 2.683e-02 | 5.961e-02
2 | 1.983e-02 | 1.337e-02
3 | 2.139e-02 | 3.481e-02
4 | 7.547e-02 | 1.094e-02
5 | 1.853e-02 | 2.318e-02
6 | 7.629e-02 | 7.954e-04
7 | 1.535e-03 | 5.618e-03
8 | 1.031e-02 | 4.970e-04
9 | 9.475e-04 | 1.477e-03
10 | 4.216e-04 | 2.476e-04
11 | 4.182e-04 | 2.802e-04
12 | 9.177e-04 | 1.190e-04
13 | 2.266e-04 | 1.540e-04
14 | 4.520e-04 | 8.620e-06
15 | 9.632e-06 | 8.918e-07
16 | 1.051e-06 | 2.823e-06
17 | 4.050e-06 | 7.415e-07
[AA] unresolved threshold state; retrying with zero-tail matching through l=0; A3 extends to 15 R_geometry
[AA] C (Z=6) | Rws=2.090843 a0 | full + external
[AA/full SCF: stage 1]
iter | d_n | d_v
[AA/full SCF: stage 2]
iter | d_n | d_v
0 | nan | 4.878e-02
1 | 3.217e-02 | 4.456e-02
2 | 2.855e-02 | 4.169e-02
3 | 2.632e-02 | 4.027e-02
4 | 2.914e-02 | 1.147e-01
5 | 3.339e-02 | 2.567e-02
6 | 2.949e-02 | 5.130e-02
7 | 2.645e-02 | 9.476e-02
8 | 1.484e-01 | 1.851e-02
9 | 5.151e-02 | 1.399e-02
10 | 2.267e-02 | 9.566e-03
11 | 2.018e-02 | 1.038e-02
12 | 1.997e-02 | 7.788e-03
13 | 1.129e-02 | 8.950e-04
14 | 1.070e-03 | 3.097e-04
15 | 8.565e-04 | 3.078e-04
16 | 7.697e-04 | 8.602e-05
17 | 1.281e-04 | 1.398e-04
18 | 1.561e-04 | 3.192e-05
19 | 4.530e-05 | 4.969e-06
20 | 1.777e-05 | 2.231e-06
21 | 3.566e-06 | 2.199e-06
22 | 2.256e-06 | 3.880e-07
23 | 6.249e-07 | 4.888e-07
24 | 1.267e-06 | 2.538e-08
25 | 3.747e-08 | 6.983e-08
26 | 1.017e-07 | 1.020e-08
27 | 2.438e-08 | 1.131e-08
[AA/external SCF]
iter | d_n | d_v
0 | 1.119e-06 | 2.458e+12
1 | 5.562e-02 | 1.041e+00
2 | 4.993e-02 | 8.909e-01
3 | 2.186e-02 | 2.894e-01
4 | 6.315e-03 | 8.785e-02
5 | 1.628e-03 | 2.433e-02
6 | 6.548e-04 | 3.519e-03
7 | 9.559e-05 | 4.177e-03
8 | 6.621e-05 | 2.644e-04
9 | 4.545e-06 | 2.752e-04
10 | 4.048e-06 | 2.637e-04
11 | 2.777e-06 | 3.045e-05
12 | 6.453e-07 | 6.872e-06
13 | 9.457e-08 | 1.239e-05
14 | 1.208e-07 | 4.021e-07
15 | 5.659e-09 | 2.169e-07
16 | 1.802e-09 | 4.620e-08
[AA] converged=True | mu=0.03486878 Ha | Zbar=3.99545673 | time=107.674 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -9.13874 | 1.00000
2s | -0.00497 | 0.51355
[ion] HNC iterations=8 QOZ=0.809 s HNC=0.015 s
[done] total time=118.445 s
Otter 0.4.1
[run] C | rho=3.51538 g/cm^3 | Te=30 eV | Ti=30 eV | QM -> QOZ/HNC
[AA] C (Z=6) | Rws=2.090843 a0 | full + external
[AA/full SCF: stage 1]
iter | d_n | d_v
[computed] c_rho3p51538_te020_ti020
0 | nan | 5.005e-02
1 | 3.477e-02 | 4.652e-02
2 | 3.309e-02 | 4.416e-02
3 | 3.316e-02 | 4.286e-02
4 | 3.459e-02 | 4.228e-02
[AA/full SCF: stage 2]
iter | d_n | d_v
0 | nan | 4.162e-02
1 | 1.833e-02 | 6.124e-02
2 | 2.199e-02 | 1.829e-02
3 | 1.426e-02 | 3.375e-02
4 | 7.421e-02 | 7.113e-03
5 | 9.737e-03 | 8.398e-03
6 | 6.567e-03 | 7.702e-03
7 | 1.598e-02 | 1.293e-03
8 | 2.227e-03 | 1.960e-03
9 | 4.167e-03 | 4.348e-04
10 | 5.630e-04 | 2.027e-04
11 | 3.793e-04 | 4.884e-05
12 | 1.004e-04 | 1.257e-05
13 | 3.239e-05 | 6.542e-07
14 | 5.791e-07 | 2.102e-06
15 | 4.417e-06 | 3.824e-07
16 | 6.922e-07 | 1.432e-06
[AA/external SCF]
iter | d_n | d_v
0 | 2.439e-06 | 2.771e+12
1 | 4.877e-02 | 1.026e+00
2 | 4.443e-02 | 9.626e-01
3 | 2.227e-02 | 2.675e-01
4 | 5.650e-03 | 6.710e-02
5 | 1.021e-03 | 1.888e-02
6 | 3.921e-04 | 3.669e-03
7 | 9.668e-05 | 2.703e-03
8 | 4.792e-05 | 3.606e-04
9 | 6.809e-06 | 1.010e-04
10 | 1.483e-06 | 2.264e-04
11 | 2.654e-06 | 1.507e-05
12 | 2.772e-07 | 2.011e-05
13 | 4.378e-07 | 1.861e-06
14 | 1.769e-08 | 4.018e-06
15 | 3.446e-08 | 7.952e-08
[AA] converged=True | mu=-0.60327426 Ha | Zbar=3.99455421 | time=14.083 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -9.19061 | 0.99959
2s | -0.00679 | 0.36795
[ion] HNC iterations=6 QOZ=0.805 s HNC=0.011 s
[done] total time=14.911 s
[computed] c_rho3p51538_te030_ti030
Otter 0.4.1
[run] C | rho=3.51538 g/cm^3 | Te=40 eV | Ti=40 eV | QM -> QOZ/HNC
[AA] C (Z=6) | Rws=2.090843 a0 | full + external
[AA/full SCF: stage 1]
iter | d_n | d_v
0 | nan | 5.081e-02
1 | 3.878e-02 | 4.737e-02
2 | 3.797e-02 | 4.484e-02
3 | 3.787e-02 | 4.301e-02
4 | 3.771e-02 | 4.146e-02
[AA/full SCF: stage 2]
iter | d_n | d_v
0 | nan | 3.939e-02
1 | 1.685e-02 | 6.368e-02
2 | 2.222e-02 | 1.976e-02
3 | 1.205e-02 | 3.063e-02
4 | 6.062e-02 | 1.064e-03
5 | 9.670e-04 | 1.032e-02
6 | 2.365e-02 | 4.830e-03
7 | 7.122e-03 | 4.645e-04
8 | 3.394e-04 | 1.379e-03
9 | 2.646e-03 | 1.374e-04
10 | 1.990e-04 | 1.666e-04
11 | 3.746e-04 | 3.561e-05
12 | 6.749e-05 | 2.351e-06
13 | 7.096e-06 | 1.137e-06
14 | 2.458e-06 | 3.094e-06
15 | 5.180e-06 | 8.519e-07
[AA/external SCF]
iter | d_n | d_v
0 | 2.635e-06 | 2.984e+12
1 | 4.233e-02 | 1.004e+00
2 | 3.785e-02 | 1.025e+00
3 | 2.142e-02 | 2.377e-01
4 | 5.096e-03 | 3.500e-02
5 | 6.727e-04 | 1.432e-02
6 | 3.815e-04 | 4.129e-03
7 | 7.226e-05 | 1.697e-03
8 | 3.993e-05 | 1.830e-04
9 | 3.973e-06 | 8.866e-05
10 | 1.363e-06 | 1.435e-04
11 | 1.767e-06 | 1.039e-05
12 | 1.741e-07 | 1.192e-05
13 | 2.278e-07 | 1.390e-06
14 | 1.386e-08 | 1.436e-06
15 | 1.109e-08 | 8.132e-08
[AA] converged=True | mu=-1.40473509 Ha | Zbar=3.99999941 | time=15.430 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -9.25613 | 0.99523
2s | -0.00984 | 0.27910
[ion] HNC iterations=6 QOZ=0.815 s HNC=0.011 s
[done] total time=16.267 s
[computed] c_rho3p51538_te040_ti040
Otter 0.4.1
[run] C | rho=3.51538 g/cm^3 | Te=50 eV | Ti=50 eV | QM -> QOZ/HNC
[AA] C (Z=6) | Rws=2.090843 a0 | full + external
[AA/full SCF: stage 1]
iter | d_n | d_v
0 | nan | 5.105e-02
1 | 4.244e-02 | 4.746e-02
2 | 4.105e-02 | 4.458e-02
3 | 3.965e-02 | 4.219e-02
4 | 3.764e-02 | 3.990e-02
[AA/full SCF: stage 2]
iter | d_n | d_v
0 | nan | 3.714e-02
1 | 2.207e-02 | 6.648e-02
2 | 2.165e-02 | 1.969e-02
3 | 1.185e-02 | 3.486e-02
4 | 6.576e-02 | 3.109e-03
5 | 3.640e-03 | 8.055e-03
6 | 1.412e-02 | 5.187e-03
7 | 5.941e-03 | 4.454e-04
8 | 5.360e-04 | 1.218e-03
9 | 1.852e-03 | 1.475e-04
10 | 2.470e-04 | 1.424e-04
11 | 3.772e-04 | 1.809e-05
12 | 4.421e-05 | 3.338e-06
13 | 1.396e-05 | 9.572e-07
14 | 3.468e-06 | 1.862e-06
15 | 4.648e-06 | 1.982e-08
16 | 4.779e-08 | 8.359e-08
[AA/external SCF]
iter | d_n | d_v
0 | 2.777e-06 | 3.154e+12
1 | 3.755e-02 | 9.831e-01
2 | 3.281e-02 | 1.066e+00
3 | 2.024e-02 | 2.270e-01
4 | 4.590e-03 | 2.509e-02
5 | 3.815e-04 | 1.341e-02
6 | 3.204e-04 | 3.339e-03
7 | 5.150e-05 | 9.816e-04
8 | 2.653e-05 | 5.034e-05
9 | 1.071e-06 | 5.855e-05
10 | 8.103e-07 | 1.151e-04
11 | 1.330e-06 | 9.606e-07
12 | 2.973e-08 | 1.424e-06
13 | 3.278e-08 | 1.393e-06
[AA] converged=True | mu=-2.32322145 Ha | Zbar=4.02710216 | time=17.036 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -9.36366 | 0.97878
2s | -0.01460 | 0.22159
[ion] HNC iterations=6 QOZ=0.809 s HNC=0.011 s
[done] total time=17.869 s
[computed] c_rho3p51538_te050_ti050
Otter 0.4.1
[run] C | rho=3.51538 g/cm^3 | Te=100 eV | Ti=100 eV | QM -> QOZ/HNC
[AA] C (Z=6) | Rws=2.090843 a0 | full + external
[AA/full SCF: stage 1]
iter | d_n | d_v
0 | nan | 4.907e-02
1 | 6.060e-02 | 4.464e-02
2 | 5.501e-02 | 4.066e-02
3 | 4.911e-02 | 3.703e-02
4 | 4.314e-02 | 3.356e-02
[AA/full SCF: stage 2]
iter | d_n | d_v
0 | nan | 3.035e-02
1 | 3.117e-02 | 2.635e-02
2 | 7.190e-02 | 6.082e-02
3 | 1.771e-02 | 1.554e-02
4 | 1.183e-02 | 2.499e-02
5 | 4.704e-02 | 2.845e-03
6 | 5.316e-02 | 5.594e-03
7 | 8.074e-03 | 8.642e-03
8 | 1.573e-02 | 2.680e-03
9 | 4.565e-03 | 7.596e-04
10 | 1.769e-03 | 3.145e-04
11 | 5.190e-04 | 6.791e-05
12 | 1.107e-04 | 2.775e-05
13 | 2.473e-05 | 1.473e-06
14 | 4.046e-06 | 1.216e-06
15 | 5.210e-06 | 1.208e-06
16 | 8.284e-06 | 7.908e-08
[AA/external SCF]
iter | d_n | d_v
0 | 3.013e-06 | 3.956e+12
1 | 2.622e-02 | 9.094e-01
2 | 2.125e-02 | 1.193e+00
3 | 1.690e-02 | 1.741e-01
4 | 3.214e-03 | 2.370e-02
5 | 4.811e-04 | 1.682e-02
6 | 2.294e-04 | 1.754e-03
7 | 3.061e-05 | 5.632e-04
8 | 1.071e-05 | 5.374e-05
9 | 1.034e-06 | 1.441e-05
10 | 1.717e-07 | 2.945e-05
11 | 3.291e-07 | 7.647e-06
12 | 9.548e-08 | 3.699e-06
13 | 4.006e-08 | 1.090e-06
[AA] converged=True | mu=-7.81908300 Ha | Zbar=4.49451417 | time=26.844 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -10.72561 | 0.68803
2s | -0.10767 | 0.10925
[ion] HNC iterations=6 QOZ=0.805 s HNC=0.011 s
[done] total time=27.671 s
[computed] c_rho3p51538_te100_ti100
T [eV] curve RMSE MAE max|delta| HNC residual
20.0 current Otter 3.1917e-02 1.8810e-02 1.4566e-01 1.9123e-06
30.0 current Otter 2.3228e-02 1.6865e-02 6.7512e-02 1.0549e-05
40.0 current Otter 1.6283e-02 1.1263e-02 6.0607e-02 3.3617e-06
50.0 current Otter 2.5546e-02 1.5872e-02 1.0632e-01 6.1170e-06
100.0 current Otter 3.4663e-02 2.5785e-02 9.2833e-02 1.5315e-06
Saved slide/web figures: benchmarks/outputs/argha_roy_carbon_sii/figures/argha_roy_carbon_sii.png, benchmarks/outputs/argha_roy_carbon_sii/figures/argha_roy_carbon_sii.pdf
from __future__ import annotations
from concurrent.futures import ProcessPoolExecutor, as_completed
import hashlib
import json
import os
from pathlib import Path
import time
from typing import Any
import matplotlib.pyplot as plt
import numpy as np
from otter import PlasmaWorkflowConfig, solve_plasma_workflow
from otter.numerics.constants import BOHR_TO_ANGSTROM
from otter.plotting import grid_figsize, save_figure, set_style
# =============================================================================
# User input
# =============================================================================
USE_PRECOMPUTED_DATA = False
if os.environ.get("OTTER_RECOMPUTE_ARGHA_CARBON", "0") == "1":
USE_PRECOMPUTED_DATA = False
# Independent states run serially; each AA keeps its one-worker default.
# All five displayed temperatures are attempted, one at a time.
MAX_STATE_WORKERS = 1
PLOT_TEMPERATURES_EV = (20.0, 30.0, 40.0, 50.0, 100.0)
REFERENCE_TEMPERATURES_EV = PLOT_TEMPERATURES_EV
STACK_OFFSET = 0.3
ACCEPTED_BASELINE_TEMPERATURES_EV = PLOT_TEMPERATURES_EV
RHO_G_CC = 3.51538
HNC_CLOSURE_TOL = 2.5e-3
K_RETAIN_MAX_BOHR_INV = 20.0
# =============================================================================
SCHEMA = "otter_argha_roy_carbon_sii_state_v1"
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" / "argha_roy_carbon_sii"
REFERENCE_DIR = (
ROOT / "benchmarks" / "reference_data" / "argha_roy_carbon_sii"
)
OUTPUT_DIR = (
ROOT
/ "benchmarks"
/ "outputs"
/ "argha_roy_carbon_sii"
/ "gallery_recomputed"
)
FIGURE_DIR = (
ROOT / "benchmarks" / "outputs" / "argha_roy_carbon_sii" / "figures"
)
def state_id(temperature_ev: float) -> str:
"""Return the stable file identifier for one temperature."""
return f"c_rho3p51538_te{int(temperature_ev):03d}_ti{int(temperature_ev):03d}"
STATES = tuple(
{
"state_id": state_id(temperature),
"element": "C",
"rho_g_cc": RHO_G_CC,
"te_ev": temperature,
"ti_ev": temperature,
}
for temperature in REFERENCE_TEMPERATURES_EV
)
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 load_npz(path: Path) -> dict[str, np.ndarray]:
"""Load a portable archive without permitting pickle/object arrays."""
with np.load(path, allow_pickle=False) as archive:
state = {key: np.asarray(archive[key]) for key in archive.files}
if any(value.dtype.hasobject for value in state.values()):
raise TypeError(f"Object arrays are forbidden in {path}.")
return state
def load_precomputed_states() -> dict[str, dict[str, np.ndarray]]:
"""Verify the separate, modern strict-Otter audit states."""
manifest_path = PRECOMPUTED_DIR / "manifest.json"
manifest = json.loads(manifest_path.read_text(encoding="utf-8"))
if manifest.get("benchmark_id") != "argha_roy_carbon_sii":
raise ValueError("Unexpected Argha-Roy benchmark manifest.")
loaded: dict[str, dict[str, np.ndarray]] = {}
for record in manifest["states"]:
path = PRECOMPUTED_DIR / str(record["baseline_file"])
if sha256_file(path) != str(record["baseline_sha256"]):
raise RuntimeError(f"Checksum mismatch for {path}.")
state = load_npz(path)
archive_id = str(state["state_id"].item())
if archive_id != str(record["state_id"]):
raise ValueError(f"State identifier mismatch in {path}.")
loaded[archive_id] = state
expected = {
state_id(temperature)
for temperature in ACCEPTED_BASELINE_TEMPERATURES_EV
}
if set(loaded) != expected:
raise RuntimeError(
"The accepted manifest does not contain the complete temperature "
f"grid: expected {sorted(expected)}, found {sorted(loaded)}."
)
return loaded
def workflow_config(state: dict[str, Any]) -> PlasmaWorkflowConfig:
"""Construct the full public Otter workflow for one state."""
return PlasmaWorkflowConfig(
elements=[str(state["element"])],
temperature_ev=float(state["te_ev"]),
ion_temperature_ev=float(state["ti_ev"]),
rho_g_cc=float(state["rho_g_cc"]),
hnc_closure_transform_tol=float(HNC_CLOSURE_TOL),
)
def strict_result(workflow: dict[str, Any]) -> tuple[dict[str, Any], dict[str, Any]]:
"""Reject unconverged electronic or HNC best-effort output."""
electronic = dict(workflow["electronic"]["result"])
ion = dict(workflow["ion"])
if electronic.get("stage2_converged") is not True:
raise RuntimeError("Full average-atom stage 2 did not converge.")
if dict(electronic.get("ext_status", {})).get("converged") is not True:
raise RuntimeError("External fixed-mu average atom did not converge.")
if str(electronic.get("threshold_state_status", "")).lower() == "unresolved":
raise RuntimeError("The threshold-state representation is unresolved.")
if ion.get("hnc_converged") is not True:
raise RuntimeError("HNC did not reach a physical fixed point.")
resolved = dict(workflow["configuration"])
if float(ion["hnc_output_residual"]) > float(resolved["hnc_tol"]):
raise RuntimeError("HNC residual exceeds the configured tolerance.")
if float(ion["closure_transform_max_abs"]) > HNC_CLOSURE_TOL:
raise RuntimeError("The g/S transform-closure audit failed.")
return electronic, ion
def solve_state(state: dict[str, Any]) -> tuple[str, dict[str, np.ndarray]]:
"""Calculate and pack one independent thermodynamic state."""
identifier = str(state["state_id"])
started = time.perf_counter()
try:
workflow = solve_plasma_workflow(workflow_config(state))
except Exception as exc:
raise RuntimeError(f"{identifier}: {exc}") from exc
elapsed_s = time.perf_counter() - started
electronic, ion = strict_result(workflow)
resolved = dict(workflow["configuration"])
k = np.asarray(ion["k"], dtype=float)
mask = k <= K_RETAIN_MAX_BOHR_INV
signature = {
"state": state,
"structure_model": "IS",
"resolved_configuration": resolved,
}
payload = {
"schema_version": np.asarray(SCHEMA),
"storage_profile": np.asarray("benchmark_analysis"),
"state_id": np.asarray(str(state["state_id"])),
"element": np.asarray(str(state["element"])),
"rho_g_cc": np.asarray(float(state["rho_g_cc"])),
"te_ev": np.asarray(float(state["te_ev"])),
"ti_ev": np.asarray(float(state["ti_ev"])),
"producer_signature_json": np.asarray(
json.dumps(signature, sort_keys=True, separators=(",", ":"))
),
"producer_elapsed_s": np.asarray(float(elapsed_s)),
"mu_ha": np.asarray(float(electronic["mu"])),
"zbar_aa": np.asarray(float(electronic["zbar"])),
"zbar_partition": np.asarray(np.asarray(ion["zbar_partition"]).item()),
"threshold_state_status": np.asarray(
str(electronic.get("threshold_state_status", "none"))
),
"threshold_state_representation": np.asarray(
str(electronic.get("threshold_state_representation", "none"))
),
"shallowest_bound_energy_ha": np.asarray(
float(electronic.get("shallowest_bound_energy_ha", np.nan))
),
"bound_zero_tail_finite_wall_energy_ha": np.asarray(
float(electronic.get("bound_zero_tail_finite_wall_energy_ha", np.nan))
),
"bound_zero_tail_matched_energy_ha": np.asarray(
float(electronic.get("bound_zero_tail_matched_energy_ha", np.nan))
),
"bound_zero_tail_exterior_probability": np.asarray(
float(electronic.get("bound_zero_tail_exterior_probability", np.nan))
),
"k_bohr_inv": k[mask],
"sii_k": np.asarray(ion["sij_k"][0, 0], dtype=float)[mask],
"hnc_output_residual": np.asarray(float(ion["hnc_output_residual"])),
"closure_transform_max_abs": np.asarray(
float(ion["closure_transform_max_abs"])
),
}
return identifier, payload
def solve_all_states() -> dict[str, dict[str, np.ndarray]]:
"""Attempt the full temperature grid with a bounded process pool.
Strictly rejected states are recorded as diagnostics. The benchmark is
accepted only when all five displayed temperatures pass every gate.
"""
loaded: dict[str, dict[str, np.ndarray]] = {}
failures: list[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]
try:
identifier, payload = future.result()
except Exception as exc:
message = f"{state['state_id']}: {exc}"
failures.append(message)
print(f"[rejected] {message}")
continue
loaded[identifier] = payload
print(f"[computed] {identifier}")
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
for identifier, payload in loaded.items():
path = OUTPUT_DIR / f"{identifier}.npz"
np.savez_compressed(path, **payload)
print(f"[saved] {path}")
rejection_path = OUTPUT_DIR / "rejections.json"
rejection_path.write_text(
json.dumps(
{
"schema_version": "otter_strict_rejections_v1",
"failures": failures,
},
indent=2,
)
+ "\n",
encoding="utf-8",
)
print(f"[saved] {rejection_path}")
required = {state_id(temperature) for temperature in PLOT_TEMPERATURES_EV}
missing_required = sorted(required.difference(loaded))
if missing_required:
detail = "\n- ".join(failures) if failures else "no error captured"
raise RuntimeError(
"A reviewed comparison state failed the strict live calculation: "
f"{missing_required}.\n- {detail}"
)
return loaded
def load_reference(temperature_ev: float) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Load one author-provided three-column DFT-MD file."""
path = REFERENCE_DIR / f"avg_S_ii_{int(temperature_ev)}eV.dat"
values = np.asarray(np.genfromtxt(path, comments="#"), dtype=float)
if values.ndim != 2 or values.shape[1] != 3:
raise ValueError(f"Expected three columns in {path}.")
mask = np.all(np.isfinite(values), axis=1)
return values[mask, 0], values[mask, 1], values[mask, 2]
def inverse_bohr_sii(
state: dict[str, np.ndarray],
) -> tuple[np.ndarray, np.ndarray]:
"""Convert one inverse-Bohr structure-factor curve to inverse ångström."""
k_angstrom_inv = (
np.asarray(state["k_bohr_inv"], dtype=float) / BOHR_TO_ANGSTROM
)
return k_angstrom_inv, np.asarray(state["sii_k"], dtype=float)
def main() -> None:
current_otter_states = (
load_precomputed_states()
if USE_PRECOMPUTED_DATA
else solve_all_states()
)
print(
"Using "
+ (
"checksummed, precomputed strict-Otter audit results."
if USE_PRECOMPUTED_DATA
else "new strict-Otter audit results calculated by this gallery script."
)
)
print(
f"{'T [eV]':>8s} {'curve':>15s} {'RMSE':>12s} {'MAE':>12s} "
f"{'max|delta|':>12s} {'HNC residual':>14s}"
)
for temperature in PLOT_TEMPERATURES_EV:
identifier = state_id(temperature)
state = current_otter_states[identifier]
k_ref, sii_ref, _ = load_reference(temperature)
k_otter, sii_otter = inverse_bohr_sii(state)
mask = (k_ref >= k_otter[0]) & (k_ref <= k_otter[-1])
delta = np.interp(k_ref[mask], k_otter, sii_otter) - sii_ref[mask]
print(
f"{temperature:8.1f} {'current Otter':>15s} "
f"{np.sqrt(np.mean(delta**2)):12.4e} "
f"{np.mean(np.abs(delta)):12.4e} {np.max(np.abs(delta)):12.4e} "
f"{float(state['hnc_output_residual']):14.4e}"
)
# %%
# Static ion structure factors
# ----------------------------
#
# Black curves are generated by the current Otter PA-HNC workflow. Orange
# curves and shaded bands are the DFT-MD values and reported uncertainty
# provided by Dr. Argha Roy. Curves are vertically offset by 0.3.
set_style("thesis", palette="bing")
fig, axis = plt.subplots(
figsize=grid_figsize(1, 1, cell_width=5.2, cell_height=5.8)
)
colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
otter_color = colors[0]
dft_md_color = colors[2]
for index in range(len(PLOT_TEMPERATURES_EV)):
axis.axhline(
1.0 + index * STACK_OFFSET,
color="0.3",
ls=":",
lw=0.7,
zorder=1,
)
for index, temperature in enumerate(PLOT_TEMPERATURES_EV):
offset = index * STACK_OFFSET
k_ref, sii_ref, uncertainty = load_reference(temperature)
state = current_otter_states[state_id(temperature)]
k_otter, sii_otter = inverse_bohr_sii(state)
axis.plot(
k_otter,
sii_otter + offset,
color=otter_color,
lw=2.2,
label="Otter PA-HNC" if index == 0 else None,
zorder=4,
)
axis.fill_between(
k_ref,
sii_ref + offset - uncertainty,
sii_ref + offset + uncertainty,
color=dft_md_color,
alpha=0.60,
linewidth=0.0,
zorder=3,
)
axis.plot(
k_ref,
sii_ref + offset,
color=dft_md_color,
lw=1.5,
label="DFT-MD" if index == 0 else None,
zorder=5,
)
axis.text(
7.5,
offset + 0.85,
rf"$T_e=T_i={temperature:g}\ \mathrm{{eV}}$",
ha="right",
va="center",
fontsize=9,
)
axis.set_xlabel(r"$k\ (\mathrm{\AA}^{-1})$")
axis.set_ylabel(r"$S_{ii}(k)$")
axis.set_xlim(0.1, 12.0)
axis.set_ylim(0.45, 2.35)
axis.set_title(r"C: $\rho = 3.515\ \mathrm{g\,cc^{-1}}$", pad=4.0)
axis.legend(loc="lower right")
fig.text(
0.5,
0.008,
"DFT-MD data: Dr. Argha Roy (private communication).",
ha="center",
va="bottom",
fontsize=7.5,
)
fig.tight_layout(rect=(0.0, 0.035, 1.0, 1.0), pad=0.45)
saved_paths = save_figure(
fig,
FIGURE_DIR / "argha_roy_carbon_sii",
formats=("png", "pdf"),
)
print(
"Saved slide/web figures: "
+ ", ".join(str(path.relative_to(ROOT)) for path in saved_paths.values())
)
plt.show()
if __name__ == "__main__":
main()