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).
USE_PRECOMPUTED_DATA = TrueVerify and load reviewed NPZ files generated by this same current-Otter workflow for all five displayed temperatures.
USE_PRECOMPUTED_DATA = FalseCalculate every average-atom and QOZ/HNC state in this same file, save NPZ files and rejection diagnostics under
benchmarks/outputs/argha_roy_carbon_sii/gallery_recomputed, and include the results only when all five states pass the reliability gates.
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].
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.plotting import grid_figsize, save_figure, set_style
# =============================================================================
# User input
# =============================================================================
USE_PRECOMPUTED_DATA = True
if os.environ.get("OTTER_RECOMPUTE_ARGHA_CARBON", "0") == "1":
USE_PRECOMPUTED_DATA = False
# Four state workers x four continuum workers uses at most about 16 workers.
# All five displayed temperatures are attempted; no more than four run
# simultaneously.
# MAX_STATE_WORKERS = 1 for a serial, memory-conservative calculation.
MAX_STATE_WORKERS = 4
CONTINUUM_WORKERS_PER_STATE = 4
QOZ_N_POINTS = 4096
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
LFC_MODEL = "chabrier1990"
HNC_TOL = 1.0e-4
HNC_CLOSURE_TOL = 2.5e-3
K_RETAIN_MAX_BOHR_INV = 20.0
# =============================================================================
BOHR_TO_ANGSTROM = 0.529177210903
SCHEMA = "otter_argha_roy_carbon_sii_state_v1"
def repository_root() -> Path:
"""Locate the Otter checkout when run directly or by Sphinx-Gallery."""
candidates = [Path.cwd().resolve(), *Path.cwd().resolve().parents]
source_file = globals().get("__file__")
if source_file is not None:
source = Path(str(source_file)).resolve()
candidates.extend([source.parent, *source.parents])
for candidate in candidates:
if (
candidate
/ "benchmarks"
/ "reference_data"
/ "argha_roy_carbon_sii"
/ "manifest.json"
).is_file():
return candidate
raise FileNotFoundError("Cannot locate the Otter checkout.")
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"]),
aa_overrides={
"cont_n_jobs": int(CONTINUUM_WORKERS_PER_STATE),
"cont_shards": int(2 * CONTINUUM_WORKERS_PER_STATE),
# Exterior-match a near-zero-energy pole only when the common
# physical SCF boundary is already asymptotic. No artificial
# extended bound-only box is introduced.
"bound_zero_tail_refine": True,
"bound_zero_tail_max_binding_ha": 1.0e-2,
"bound_zero_tail_scan_points": 64,
"bound_zero_tail_edge_rel_tol": 0.1,
},
hnc_tol=float(HNC_TOL),
hnc_closure_transform_tol=float(HNC_CLOSURE_TOL),
hnc_max_iter=1000,
)
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.")
if float(ion["hnc_output_residual"]) > 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)
k = np.asarray(ion["k"], dtype=float)
mask = k <= K_RETAIN_MAX_BOHR_INV
signature = {
"state": state,
"electronic_model": "qm",
"structure_model": "IS",
"aa_n_points": 4096,
"bound_occ_mode": "fd",
"bound_rmax_mult": None,
"bound_zero_tail_refine": True,
"bound_zero_tail_max_binding_ha": 1.0e-2,
"bound_zero_tail_scan_points": 64,
"bound_zero_tail_edge_rel_tol": 0.1,
"b3_tail_model": "full",
"qoz_n_points": QOZ_N_POINTS,
"qoz_zbar_mode": "pseudoatom_partition",
"qoz_renormalize_nscr_to_zbar": True,
"chi0_model": "lindhard_fd",
"lfc_model": LFC_MODEL,
"hnc_tol": HNC_TOL,
"hnc_closure_tol": HNC_CLOSURE_TOL,
}
payload = {
"schema_version": np.asarray(SCHEMA),
"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(float(ion["zbar_partition"])),
"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["sii_k"], dtype=float)[mask],
"vii_k_ha_bohr3": np.asarray(ion["vii_k"], dtype=float)[mask],
"n_scr_k_electrons": np.asarray(ion["n_scr_k"], 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)
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}"
)
Using checksummed, precomputed strict-Otter audit results.
T [eV] curve RMSE MAE max|delta| HNC residual
20.0 current Otter 3.3529e-02 1.9374e-02 1.5597e-01 7.0657e-05
30.0 current Otter 2.3529e-02 1.6936e-02 7.1308e-02 1.7040e-05
40.0 current Otter 1.6295e-02 1.1288e-02 6.0229e-02 4.0242e-06
50.0 current Otter 2.5523e-02 1.5854e-02 1.0621e-01 5.6363e-06
100.0 current Otter 3.4620e-02 2.5773e-02 9.2713e-02 1.3926e-06
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())
)
if "agg" not in plt.get_backend().lower():
plt.show()
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