Note
Go to the end to download the full example code.
Ion-structure literature library
This benchmark compares Otter ion-structure results with literature curves for
aluminium, beryllium, and carbon. USE_PRECOMPUTED_DATA = True verifies
and loads checksummed Otter NPZ files. With False, this file
constructs otter.PlasmaWorkflowConfig, evaluates every average atom
and ion structure, saves new NPZ files below
benchmarks/outputs/ion_structure_library/gallery_recomputed, and plots
those results.
In the four-panel \(S_{ii}(k)\) figure, panel 1 uses Gill et al., Fig. 3 [Gill et al., 2015]; panels 2 and 3 use Clérouin et al., Fig. 1 [Clérouin et al., 2015]; and panel 4 uses Wünsch et al., Fig. 2 [Wünsch et al., 2009]. For Wünsch Be, ordinary HNC, Rosenfeld–Ashcroft VMHNC, and same-potential LAMMPS MD reuse one IS-QOZ pair potential. The MD error band is twice the standard error across independent RDF blocks or saved-frame reciprocal-shell averages. The real-space comparison uses Fig. 1(c). The carbon PA-HNC data were provided by C. E. Starrett (private communication; unpublished). Reference coordinates remain separate data files because they are digitized publication data; all unit conversions are explicit below.
The average-atom/pseudoatom construction follows
Starrett and Saumon [2014], Starrett and Saumon [2013], and the default
finite-temperature jellium LFC follows Chabrier [1990]. See
the provenance, coordinate units, and data notice for the scientific interpretation.
Both benchmark figures are exported as matching PNG and vector PDF files
under benchmarks/outputs/ion_structure_library/figures.
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 sys
import time
from typing import Any
import matplotlib.pyplot as plt
import numpy as np
from scipy.constants import physical_constants
from otter import (
PlasmaWorkflowConfig,
continue_plasma_workflow_from_electronic_result,
solve_plasma_workflow,
)
from otter.plotting import (
MODEL_STYLES,
grid_figsize,
save_figure,
style_context,
)
# =============================================================================
# User input
# =============================================================================
USE_PRECOMPUTED_DATA = True
if os.environ.get("OTTER_RECOMPUTE_ION_STRUCTURE_LIBRARY", "0") == "1":
USE_PRECOMPUTED_DATA = False
# Three independent state groups, each with six continuum workers. The two
# Al 8.1-g/cc states share one electronic calculation because only the ion
# temperature differs.
MAX_STATE_WORKERS = 3
CONTINUUM_WORKERS_PER_STATE = 6
HNC_TOL = 1.0e-4
HNC_CLOSURE_TOL = 2.5e-3
R_RETAIN_MAX_BOHR = 20.0
K_RETAIN_MAX_BOHR_INV = 20.0
VMHNC_ETA_TOL = 1.0e-6
RUN_WUNSCH_SAME_POTENTIAL_MD = True
WUNSCH_MD_MPI_PROCESSES = 10
WUNSCH_MD_ATOMS = 2048
WUNSCH_MD_K_MAX_ANGSTROM_INV = 10.2
# The fundamental cubic-cell shell has only three independent half-space
# vectors and is visibly direction-starved. Keep it in the NPZ audit data,
# but require at least the next shell's four modes when drawing MD S(k).
MD_MIN_HALF_SPACE_MODES_PER_BIN = 4
# =============================================================================
BOHR_TO_ANGSTROM = 0.529177210903
STATE_GROUPS: dict[str, tuple[dict[str, Any], ...]] = {
"al_gill": (
{
"state_id": "al_gill_rho2p7_te5_ti5",
"element": "Al",
"rho_g_cc": 2.7,
"te_ev": 5.0,
"ti_ev": 5.0,
},
),
"al_clerouin": (
{
"state_id": "al_clerouin_rho8p1_te10_ti10",
"element": "Al",
"rho_g_cc": 8.1,
"te_ev": 10.0,
"ti_ev": 10.0,
},
{
"state_id": "al_clerouin_rho8p1_te10_ti2",
"element": "Al",
"rho_g_cc": 8.1,
"te_ev": 10.0,
"ti_ev": 2.0,
},
),
"al_clerouin_tf": (
{
"state_id": "al_clerouin_rho8p1_te10_ti10_tf",
"element": "Al",
"rho_g_cc": 8.1,
"te_ev": 10.0,
"ti_ev": 10.0,
"electronic_model": "tf",
},
{
"state_id": "al_clerouin_rho8p1_te10_ti2_tf",
"element": "Al",
"rho_g_cc": 8.1,
"te_ev": 10.0,
"ti_ev": 2.0,
"electronic_model": "tf",
},
),
"be_wunsch": (
{
"state_id": "be_wunsch_rho5p544_te13_ti13",
"element": "Be",
"rho_g_cc": 5.544,
"te_ev": 13.0,
"ti_ev": 13.0,
},
),
"c_starrett_hot": (
{
"state_id": "c_starrett_rho20_te50_ti50",
"element": "C",
"rho_g_cc": 20.0,
"te_ev": 50.0,
"ti_ev": 50.0,
},
),
}
REFERENCE_SERIES: dict[str, tuple[dict[str, str], ...]] = {
"al_gill_rho2p7_te5_ti5": (
{
"observable": "sii",
"label": "Gill KS-PAMD",
"file": "gill_et_al_2015/Sii_Al_T5ev_rho2.7_KS-PAMD_Gill.csv",
"x_unit": "angstrom^-1",
},
{
"observable": "sii",
"label": "Gill TF-PAMD",
"file": "gill_et_al_2015/Sii_Al_T5ev_rho2.7_TF-PAMD_Gill.csv",
"x_unit": "angstrom^-1",
},
{
"observable": "sii",
"label": "Gill TF-DFT-MD",
"file": "gill_et_al_2015/Sii_Al_T5ev_rho2.7_QMD_Gill.csv",
"x_unit": "angstrom^-1",
},
),
"al_clerouin_rho8p1_te10_ti10": (
{
"observable": "sii",
"label": "Clérouin OFMD",
"file": (
"clerouin_et_al_2015/"
"Jean_2015_Al_rho8.1_Te10.0_Ti10.csv"
),
"x_unit": "angstrom^-1",
},
),
"al_clerouin_rho8p1_te10_ti2": (
{
"observable": "sii",
"label": "Clérouin OFMD",
"file": (
"clerouin_et_al_2015/"
"Jean_2015_Al_rho8.1_Te10.0_Ti2.0_OFMD.csv"
),
"x_unit": "angstrom^-1",
},
{
"observable": "sii",
"label": "HNC-Y-SRR",
"file": (
"clerouin_et_al_2015/"
"Jean_2015_Al_rho8.1_Te10.0_Ti2.0_SRR.csv"
),
"x_unit": "angstrom^-1",
},
),
"be_wunsch_rho5p544_te13_ti13": (
{
"observable": "sii",
"label": "Wünsch DFT-MD",
"file": (
"wunsch_et_al_2009/"
"Sii_DFTMD_Be_3rho0_T13_Z2_wunsch2009.csv"
),
"x_unit": "angstrom^-1",
},
{
"observable": "sii",
"label": "HNC-Y-SRR",
"file": (
"wunsch_et_al_2009/"
"Sii_HNC-YSRR_Be_3rho0_T13_Z2_wunsch2009.csv"
),
"x_unit": "angstrom^-1",
},
{
"observable": "sii",
"label": "HNC-KK",
"file": (
"wunsch_et_al_2009/"
"Sii_HNCKK_Be_3rho0_T13_Z2_wunsch2009.csv"
),
"x_unit": "angstrom^-1",
},
{
"observable": "sii",
"label": "HNC-Y",
"file": (
"wunsch_et_al_2009/"
"Sii_HNCY_Be_3rho0_T13_Z2_wunsch2009.csv"
),
"x_unit": "angstrom^-1",
},
{
"observable": "gii",
"label": "Wünsch DFT-MD",
"file": (
"wunsch_et_al_2009/"
"gii_DFTMD_Be_3rho0_T13_Z2_wunsch2009.csv"
),
"x_unit": "angstrom",
},
{
"observable": "gii",
"label": "HNC-Y-SRR",
"file": (
"wunsch_et_al_2009/"
"gii_HNC-YSRR_Be_3rho0_T13_Z2_wunsch2009.csv"
),
"x_unit": "angstrom",
},
{
"observable": "gii",
"label": "HNC-KK",
"file": (
"wunsch_et_al_2009/"
"gii_HNC-KK_Be_3rho0_T13_Z2_wunsch2009.csv"
),
"x_unit": "angstrom",
},
{
"observable": "gii",
"label": "HNC-Y",
"file": (
"wunsch_et_al_2009/"
"gii_HNC-Y_Be_3rho0_T13_Z2_wunsch2009.csv"
),
"x_unit": "angstrom",
},
),
"c_starrett_rho20_te50_ti50": (
{
"observable": "gii",
"label": "Starrett PA-HNC",
"file": (
"starrett_saumon_2013/"
"gii_C_20gcc_50.0ev_starrett.csv"
),
"x_unit": "bohr",
},
),
}
STATE_TITLES = {
"al_gill_rho2p7_te5_ti5": (
r"Al: $\rho=2.7$ g cm$^{-3}$, $T_e=T_i=5$ eV"
),
"al_clerouin_rho8p1_te10_ti10": (
r"Al: $\rho=8.1$ g cm$^{-3}$, $T_e=T_i=10$ eV"
),
"al_clerouin_rho8p1_te10_ti2": (
r"Al: $\rho=8.1$ g cm$^{-3}$, $T_e=10$, $T_i=2$ eV"
),
"be_wunsch_rho5p544_te13_ti13": (
r"Be: $\rho=5.544$ g cm$^{-3}$, $T_e=T_i=13$ eV"
),
"c_starrett_rho20_te50_ti50": (
r"C: $\rho=20$ g cm$^{-3}$, $T_e=T_i=50$ eV"
),
}
OTTER_SERIES = {
state_id: ((state_id, "Otter KS", "-", ""),)
for state_id in REFERENCE_SERIES
}
OTTER_SERIES.update(
{
"al_clerouin_rho8p1_te10_ti10": (
("al_clerouin_rho8p1_te10_ti10", "Otter KS", "-", ""),
("al_clerouin_rho8p1_te10_ti10_tf", "Otter TF", "--", ""),
),
"al_clerouin_rho8p1_te10_ti2": (
("al_clerouin_rho8p1_te10_ti2", "Otter KS", "-", ""),
("al_clerouin_rho8p1_te10_ti2_tf", "Otter TF", "--", ""),
),
"be_wunsch_rho5p544_te13_ti13": (
("be_wunsch_rho5p544_te13_ti13", "Otter-HNC", "-", ""),
(
"be_wunsch_rho5p544_te13_ti13",
"Otter-VMHNC",
"--",
"vmhnc_",
),
("be_wunsch_rho5p544_te13_ti13", "Otter-MD", "-.", "md_"),
),
}
)
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"
/ "baselines"
/ "ion_structure_library"
/ "manifest.json"
).is_file():
return candidate
raise FileNotFoundError("Cannot locate the Otter checkout.")
ROOT = repository_root()
if str(ROOT) not in sys.path:
sys.path.insert(0, str(ROOT))
PRECOMPUTED_DIR = (
ROOT / "benchmarks" / "baselines" / "ion_structure_library"
)
REFERENCE_DIR = (
ROOT / "benchmarks" / "reference_data" / "ion_structure_library"
)
OUTPUT_DIR = (
ROOT
/ "benchmarks"
/ "outputs"
/ "ion_structure_library"
/ "gallery_recomputed"
)
FIGURE_DIR = (
ROOT / "benchmarks" / "outputs" / "ion_structure_library" / "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 load_npz(path: Path) -> dict[str, np.ndarray]:
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 every accepted result checksum before plotting."""
manifest = json.loads(
(PRECOMPUTED_DIR / "manifest.json").read_text(encoding="utf-8")
)
if manifest.get("benchmark_id") != "ion_structure_library":
raise ValueError("Unexpected ion-structure 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)
state_id = str(state["state_id"].item())
if state_id != str(record["state_id"]):
raise ValueError(f"State identifier mismatch in {path}.")
loaded[state_id] = state
return loaded
def workflow_config(
state: dict[str, Any],
*,
ion_temperature_ev: float | None,
bridge_model: str = "none",
) -> PlasmaWorkflowConfig:
"""Build the complete public Otter workflow for one thermodynamic state."""
return PlasmaWorkflowConfig(
elements=[str(state["element"])],
temperature_ev=float(state["te_ev"]),
ion_temperature_ev=(
None if ion_temperature_ev is None else float(ion_temperature_ev)
),
rho_g_cc=float(state["rho_g_cc"]),
electronic_model=str(state.get("electronic_model", "qm")),
aa_overrides={
"cont_n_jobs": int(CONTINUUM_WORKERS_PER_STATE),
"cont_shards": int(2 * CONTINUUM_WORKERS_PER_STATE),
},
hnc_tol=float(HNC_TOL),
hnc_closure_transform_tol=float(HNC_CLOSURE_TOL),
hnc_max_iter=500,
hnc_bridge_model=bridge_model,
vmhnc_eta_tol=VMHNC_ETA_TOL,
)
def strict_check(
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 pack_result(
workflow: dict[str, Any],
state: dict[str, Any],
*,
elapsed_s: float,
) -> dict[str, np.ndarray]:
"""Keep the portable arrays needed by this benchmark and downstream use."""
electronic, ion = strict_check(workflow)
r_e = np.asarray(electronic["r"], dtype=float)
r = np.asarray(ion["r"], dtype=float)
k = np.asarray(ion["k"], dtype=float)
e_mask = r_e <= R_RETAIN_MAX_BOHR
r_mask = r <= R_RETAIN_MAX_BOHR
k_mask = k <= K_RETAIN_MAX_BOHR_INV
return {
"schema_version": np.asarray("otter_gallery_ion_library_v1"),
"state_id": np.asarray(str(state["state_id"])),
"electronic_model": np.asarray(
str(state.get("electronic_model", "qm"))
),
"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_elapsed_s": np.asarray(float(elapsed_s)),
"r_e_bohr": r_e[e_mask],
"n_full_bohr3": np.asarray(electronic["n_full"])[e_mask],
"n_scr_bohr3": np.asarray(electronic["n_scr"])[e_mask],
"r_bohr": r[r_mask],
"gii_r": np.asarray(ion["gii_r"])[r_mask],
"k_bohr_inv": k[k_mask],
"sii_k": np.asarray(ion["sii_k"])[k_mask],
"vii_k_ha_bohr3": np.asarray(ion["vii_k"])[k_mask],
"n_scr_k_electrons": np.asarray(ion["n_scr_k"])[k_mask],
"zbar_partition": np.asarray(float(ion["zbar_partition"])),
"hnc_best_residual": np.asarray(float(ion["hnc_output_residual"])),
"hnc_closure_mismatch": np.asarray(
float(ion["closure_transform_max_abs"])
),
}
def add_wunsch_vmhnc(
payload: dict[str, np.ndarray],
hnc_workflow: dict[str, Any],
vmhnc_workflow: dict[str, Any],
) -> None:
"""Store the second closure after checking that its potential is unchanged."""
_, hnc = strict_check(hnc_workflow)
_, vmhnc = strict_check(vmhnc_workflow)
if str(hnc["hnc_bridge_model"]) != "none":
raise RuntimeError("The ordinary-HNC Wünsch curve used a bridge.")
if str(vmhnc["hnc_bridge_model"]) != "rosenfeld_ashcroft":
raise RuntimeError("The Wünsch VMHNC curve used the wrong bridge.")
if not np.allclose(hnc["vii_r"], vmhnc["vii_r"], rtol=0.0, atol=0.0):
raise RuntimeError("Wünsch HNC and VMHNC did not reuse one potential.")
r = np.asarray(vmhnc["r"], dtype=float)
k = np.asarray(vmhnc["k"], dtype=float)
r_mask, k_mask = r <= R_RETAIN_MAX_BOHR, k <= K_RETAIN_MAX_BOHR_INV
payload.update(
{
"vmhnc_r_bohr": r[r_mask],
"vmhnc_gii_r": np.asarray(vmhnc["gii_r"])[r_mask],
"vmhnc_k_bohr_inv": k[k_mask],
"vmhnc_sii_k": np.asarray(vmhnc["sii_k"])[k_mask],
"vmhnc_best_residual": np.asarray(float(vmhnc["hnc_output_residual"])),
"vmhnc_closure_mismatch": np.asarray(
float(vmhnc["closure_transform_max_abs"])
),
"vmhnc_eta": np.asarray(float(vmhnc["vmhnc_eta"])),
"vmhnc_variational_residual": np.asarray(
float(vmhnc["vmhnc_variational_residual"])
),
}
)
def add_wunsch_md(
payload: dict[str, np.ndarray],
workflow: dict[str, Any],
) -> None:
"""Run the reusable single-/multi-species LAMMPS driver for Wünsch Be."""
from tools.otter_lammps_md import (
MDConfig,
MDSpecies,
PairPotential,
run_otter_lammps_md,
)
ion = dict(workflow["ion"])
mass_u = 9.0121831
mass_ratio = (
physical_constants["atomic mass constant"][0]
/ physical_constants["electron mass"][0]
)
atomic_time_ps = physical_constants["atomic unit of time"][0] * 1.0e12
omega_p = np.sqrt(
4.0
* np.pi
* float(ion["n_i"])
* float(ion["zbar_partition"]) ** 2
/ (mass_u * mass_ratio)
)
state_id = str(payload["state_id"].item())
result = run_otter_lammps_md(
MDConfig(
output_dir=OUTPUT_DIR / "md_work" / state_id,
species=(MDSpecies("Be", mass_u, WUNSCH_MD_ATOMS),),
ion_density_bohr3=float(ion["n_i"]),
ion_temperature_ev=float(payload["ti_ev"]),
timestep_ps=0.005 / omega_p * atomic_time_ps,
thermostat_damp_ps=0.5 / omega_p * atomic_time_ps,
equilibration_steps=10_000,
production_steps=100_000,
rdf_bins=500,
rdf_every=100,
rdf_repeat=50,
trajectory_every=5_000,
table_points=8_192,
k_max_angstrom_inv=WUNSCH_MD_K_MAX_ANGSTROM_INV,
random_seed=20_260_825,
mpi_processes=WUNSCH_MD_MPI_PROCESSES,
structure_factor_workers=WUNSCH_MD_MPI_PROCESSES,
),
[PairPotential("Be", "Be", np.asarray(ion["r"]), np.asarray(ion["vii_r"]))],
)
payload.update(result)
payload["md_gii_r"] = np.asarray(result["md_gij_r"])[:, 0]
payload["md_gii_block_sem"] = np.asarray(result["md_gij_block_sem"])[:, 0]
payload["md_sii_k"] = np.asarray(result["md_snn_k"])
payload["md_sii_frame_sem"] = np.asarray(result["md_snn_frame_sem"])
payload["md_sii_vectors_per_bin"] = np.asarray(
result["md_vectors_per_k_bin"]
)
def solve_group(
group: tuple[dict[str, Any], ...],
) -> dict[str, dict[str, np.ndarray]]:
"""Solve one electronic state and reuse it when only Ti changes."""
first = group[0]
started = time.perf_counter()
electronic_only = solve_plasma_workflow(
workflow_config(first, ion_temperature_ev=None)
)
electronic_elapsed = time.perf_counter() - started
electronic_kind = str(electronic_only["electronic"]["kind"])
electronic = dict(electronic_only["electronic"]["result"])
solved: dict[str, dict[str, np.ndarray]] = {}
for state in group:
ion_started = time.perf_counter()
workflow = continue_plasma_workflow_from_electronic_result(
workflow_config(
state,
ion_temperature_ev=float(state["ti_ev"]),
),
electronic_kind=electronic_kind,
electronic_result=electronic,
)
elapsed = electronic_elapsed + time.perf_counter() - ion_started
payload = pack_result(workflow, state, elapsed_s=elapsed)
if str(state["state_id"]).startswith("be_wunsch"):
vmhnc = continue_plasma_workflow_from_electronic_result(
workflow_config(
state,
ion_temperature_ev=float(state["ti_ev"]),
bridge_model="rosenfeld_ashcroft",
),
electronic_kind=electronic_kind,
electronic_result=electronic,
)
add_wunsch_vmhnc(payload, workflow, vmhnc)
if RUN_WUNSCH_SAME_POTENTIAL_MD:
add_wunsch_md(payload, workflow)
solved[str(state["state_id"])] = payload
return solved
def solve_all_states() -> dict[str, dict[str, np.ndarray]]:
"""Calculate all state groups with a bounded process pool."""
groups = tuple(STATE_GROUPS.values())
loaded: dict[str, dict[str, np.ndarray]] = {}
with ProcessPoolExecutor(max_workers=MAX_STATE_WORKERS) as pool:
futures = {pool.submit(solve_group, group): group for group in groups}
for future in as_completed(futures):
group_result = future.result()
loaded.update(group_result)
print("[computed] " + ", ".join(sorted(group_result)))
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
for state_id, state in loaded.items():
path = OUTPUT_DIR / f"{state_id}.npz"
np.savez_compressed(path, **state)
print(f"[saved] {path}")
return loaded
def load_reference(series: dict[str, str]) -> tuple[np.ndarray, np.ndarray]:
"""Load one two-column digitization without changing its coordinate unit."""
path = (REFERENCE_DIR / series["file"]).resolve()
if not path.is_relative_to(REFERENCE_DIR.resolve()):
raise ValueError("Reference path escapes its package directory.")
values = np.asarray(
np.genfromtxt(path, delimiter=",", comments="#"),
dtype=float,
)
if values.ndim != 2 or values.shape[1] < 2:
raise ValueError(f"Expected two reference columns in {path}.")
mask = np.isfinite(values[:, 0]) & np.isfinite(values[:, 1])
return values[mask, 0], values[mask, 1]
def otter_curve(
state: dict[str, np.ndarray],
observable: str,
x_unit: str,
prefix: str = "",
) -> tuple[np.ndarray, np.ndarray]:
"""Convert only the Otter coordinate to the publication's stated unit."""
if observable == "sii":
x = np.asarray(state[f"{prefix}k_bohr_inv"], dtype=float)
y = np.asarray(state[f"{prefix}sii_k"], dtype=float)
if prefix == "md_":
reliable = (
np.asarray(state["md_sii_vectors_per_bin"], dtype=int)
>= MD_MIN_HALF_SPACE_MODES_PER_BIN
)
x, y = x[reliable], y[reliable]
if x_unit == "angstrom^-1":
x = x / BOHR_TO_ANGSTROM
elif x_unit != "bohr^-1":
raise ValueError(f"Unsupported reciprocal unit {x_unit!r}.")
return x, y
if observable == "gii":
x = np.asarray(state[f"{prefix}r_bohr"], dtype=float)
y = np.asarray(state[f"{prefix}gii_r"], dtype=float)
if x_unit == "angstrom":
x = x * BOHR_TO_ANGSTROM
elif x_unit != "bohr":
raise ValueError(f"Unsupported radius unit {x_unit!r}.")
return x, y
raise ValueError(f"Unknown observable {observable!r}.")
def print_metrics(states: dict[str, dict[str, np.ndarray]]) -> None:
print(
f"{'state':42s} {'model':10s} {'obs':>3s} {'reference':19s} "
f"{'RMSE':>10s} {'MAE':>10s} {'max':>10s}"
)
for state_id, series_list in REFERENCE_SERIES.items():
for result_id, model_label, _, prefix in OTTER_SERIES[state_id]:
for series in series_list:
x_ref, y_ref = load_reference(series)
x_otter, y_otter = otter_curve(
states[result_id],
str(series["observable"]),
str(series["x_unit"]),
prefix,
)
mask = (x_ref >= x_otter[0]) & (x_ref <= x_otter[-1])
delta = (
np.interp(x_ref[mask], x_otter, y_otter) - y_ref[mask]
)
print(
f"{state_id:42s} {model_label:10s} "
f"{series['observable']:>3s} "
f"{series['label'][:19]:19s} "
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}"
)
states = (
load_precomputed_states()
if USE_PRECOMPUTED_DATA
else solve_all_states()
)
print(
"Using "
+ (
"checksummed, precomputed Otter results."
if USE_PRECOMPUTED_DATA
else "new results calculated directly by this gallery script."
)
)
print_metrics(states)
def plot_observable(
*,
observable: str,
state_ids: tuple[str, ...],
) -> plt.Figure:
"""Draw one self-contained comparison figure."""
ncols = 2
nrows = (len(state_ids) + ncols - 1) // ncols
fig, axes = plt.subplots(
nrows,
ncols,
figsize=grid_figsize(nrows, ncols),
squeeze=False,
)
marker_cycle = ("o", "s", "^", "x")
colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
for panel, state_id in enumerate(state_ids):
axis = axes.ravel()[panel]
series_list = [
item
for item in REFERENCE_SERIES[state_id]
if item["observable"] == observable
]
display_unit = str(series_list[0]["x_unit"])
for model_index, (result_id, label, line_style, prefix) in enumerate(
OTTER_SERIES[state_id]
):
x_otter, y_otter = otter_curve(
states[result_id],
observable,
display_unit,
prefix,
)
style = dict(MODEL_STYLES["otter"])
style["linestyle"] = line_style
style["alpha"] = 0.82
if model_index == 1:
style["color"] = colors[0]
elif model_index == 2:
style["color"] = colors[1]
axis.plot(x_otter, y_otter, label=label, **style)
sem_key = (
"md_sii_frame_sem"
if prefix == "md_" and observable == "sii"
else "md_gii_block_sem"
)
if prefix == "md_" and sem_key in states[result_id]:
sem = np.asarray(states[result_id][sem_key], dtype=float)
if observable == "sii":
reliable = (
np.asarray(
states[result_id]["md_sii_vectors_per_bin"],
dtype=int,
)
>= MD_MIN_HALF_SPACE_MODES_PER_BIN
)
sem = sem[reliable]
axis.fill_between(
x_otter,
y_otter - 2.0 * sem,
y_otter + 2.0 * sem,
color=style["color"],
alpha=0.14,
linewidth=0.0,
)
reference_x: list[np.ndarray] = []
for index, series in enumerate(series_list):
x_ref, y_ref = load_reference(series)
reference_x.append(x_ref)
marker = marker_cycle[index % len(marker_cycle)]
scatter_options: dict[str, Any] = {
"s": 25,
"marker": marker,
"linewidths": 1.2,
"label": series["label"],
"zorder": 3,
}
color = colors[index % len(colors)]
if marker == "x":
scatter_options["color"] = color
else:
scatter_options["facecolors"] = "none"
scatter_options["edgecolors"] = color
axis.scatter(x_ref, y_ref, **scatter_options)
axis.set_title(STATE_TITLES[state_id], fontsize=10)
axis.set_ylabel(
r"$S_{ii}(k)$" if observable == "sii" else r"$g_{ii}(r)$"
)
if display_unit == "angstrom^-1":
axis.set_xlabel(r"$k$ [$\mathrm{\AA}^{-1}$]")
elif display_unit == "angstrom":
axis.set_xlabel(r"$r$ [$\mathrm{\AA}$]")
else:
axis.set_xlabel(r"$r$ [Bohr]")
all_reference_x = np.concatenate(reference_x)
span = float(np.ptp(all_reference_x))
margin = max(0.03 * span, 1.0e-6)
left = (
-0.5
if observable == "gii"
else max(0.0, float(np.min(all_reference_x)) - margin)
)
axis.set_xlim(left, float(np.max(all_reference_x)) + margin)
axis.axhline(1.0, color="0.55", lw=0.8, ls=":")
axis.legend(fontsize="small")
for panel in range(len(state_ids), axes.size):
axes.ravel()[panel].set_visible(False)
fig.suptitle("Otter ion structure versus curated literature curves", y=0.985)
source_line = (
"Reference data: Gill et al. (2015); Clérouin et al. (2015); "
"Wunsch et al. (2009)."
if observable == "sii"
else "Reference data: Wunsch et al. (2009); "
"C. E. Starrett (private communication)."
)
fig.text(
0.5,
0.006,
source_line,
ha="center",
va="bottom",
fontsize=7.5,
)
fig.tight_layout(rect=(0.0, 0.035, 1.0, 0.965), pad=0.55)
return fig
Using checksummed, precomputed Otter results.
state model obs reference RMSE MAE max
al_gill_rho2p7_te5_ti5 Otter KS sii Gill KS-PAMD 1.1377e-02 8.5291e-03 3.2706e-02
al_gill_rho2p7_te5_ti5 Otter KS sii Gill TF-PAMD 1.6740e-02 1.2652e-02 3.9127e-02
al_gill_rho2p7_te5_ti5 Otter KS sii Gill TF-DFT-MD 4.2915e-02 3.2646e-02 7.9026e-02
al_clerouin_rho8p1_te10_ti10 Otter KS sii Clérouin OFMD 3.8487e-02 2.8896e-02 7.4072e-02
al_clerouin_rho8p1_te10_ti10 Otter TF sii Clérouin OFMD 1.0938e-02 8.1710e-03 3.0392e-02
al_clerouin_rho8p1_te10_ti2 Otter KS sii Clérouin OFMD 2.4625e-02 1.4583e-02 8.6002e-02
al_clerouin_rho8p1_te10_ti2 Otter KS sii HNC-Y-SRR 4.2045e-02 3.2284e-02 1.1667e-01
al_clerouin_rho8p1_te10_ti2 Otter TF sii Clérouin OFMD 1.9277e-02 1.1336e-02 8.0555e-02
al_clerouin_rho8p1_te10_ti2 Otter TF sii HNC-Y-SRR 3.1915e-02 2.4906e-02 7.2462e-02
be_wunsch_rho5p544_te13_ti13 Otter-HNC sii Wünsch DFT-MD 2.5938e-02 2.0125e-02 5.2334e-02
be_wunsch_rho5p544_te13_ti13 Otter-HNC sii HNC-Y-SRR 1.6678e-02 1.3711e-02 3.2287e-02
be_wunsch_rho5p544_te13_ti13 Otter-HNC sii HNC-KK 7.4336e-02 4.6061e-02 1.7745e-01
be_wunsch_rho5p544_te13_ti13 Otter-HNC sii HNC-Y 7.3473e-02 5.4381e-02 1.3616e-01
be_wunsch_rho5p544_te13_ti13 Otter-HNC gii Wünsch DFT-MD 6.3554e-02 4.6868e-02 1.3475e-01
be_wunsch_rho5p544_te13_ti13 Otter-HNC gii HNC-Y-SRR 5.1948e-02 3.4890e-02 1.2638e-01
be_wunsch_rho5p544_te13_ti13 Otter-HNC gii HNC-KK 4.6212e-02 3.3795e-02 8.8908e-02
be_wunsch_rho5p544_te13_ti13 Otter-HNC gii HNC-Y 9.4323e-02 5.9712e-02 2.0424e-01
be_wunsch_rho5p544_te13_ti13 Otter-VMHNC sii Wünsch DFT-MD 2.4059e-02 1.9083e-02 4.6395e-02
be_wunsch_rho5p544_te13_ti13 Otter-VMHNC sii HNC-Y-SRR 1.4503e-02 1.1955e-02 2.8401e-02
be_wunsch_rho5p544_te13_ti13 Otter-VMHNC sii HNC-KK 7.3093e-02 4.5863e-02 1.7529e-01
be_wunsch_rho5p544_te13_ti13 Otter-VMHNC sii HNC-Y 7.6669e-02 5.7311e-02 1.4034e-01
be_wunsch_rho5p544_te13_ti13 Otter-VMHNC gii Wünsch DFT-MD 5.9294e-02 4.4025e-02 1.2050e-01
be_wunsch_rho5p544_te13_ti13 Otter-VMHNC gii HNC-Y-SRR 4.6008e-02 3.1068e-02 1.1127e-01
be_wunsch_rho5p544_te13_ti13 Otter-VMHNC gii HNC-KK 4.2729e-02 3.1024e-02 8.4241e-02
be_wunsch_rho5p544_te13_ti13 Otter-VMHNC gii HNC-Y 1.0098e-01 6.3375e-02 2.2004e-01
be_wunsch_rho5p544_te13_ti13 Otter-MD sii Wünsch DFT-MD 2.5642e-02 2.0517e-02 5.4641e-02
be_wunsch_rho5p544_te13_ti13 Otter-MD sii HNC-Y-SRR 1.2866e-02 1.0447e-02 2.5718e-02
be_wunsch_rho5p544_te13_ti13 Otter-MD sii HNC-KK 6.2927e-02 3.9791e-02 1.6995e-01
be_wunsch_rho5p544_te13_ti13 Otter-MD sii HNC-Y 7.3797e-02 5.3634e-02 1.4091e-01
be_wunsch_rho5p544_te13_ti13 Otter-MD gii Wünsch DFT-MD 6.1362e-02 4.5209e-02 1.2454e-01
be_wunsch_rho5p544_te13_ti13 Otter-MD gii HNC-Y-SRR 4.7961e-02 3.2206e-02 1.1615e-01
be_wunsch_rho5p544_te13_ti13 Otter-MD gii HNC-KK 4.0131e-02 2.9536e-02 8.0302e-02
be_wunsch_rho5p544_te13_ti13 Otter-MD gii HNC-Y 1.0101e-01 6.3456e-02 2.2009e-01
c_starrett_rho20_te50_ti50 Otter KS gii Starrett PA-HNC 6.5747e-03 5.1104e-03 1.6421e-02
Static ion structure factors
Literature wave numbers are stored in inverse ångström. Otter’s native
inverse-Bohr grid is converted explicitly inside otter_curve.
plot_style = ExitStack()
plot_style.enter_context(style_context("thesis", palette="bing"))
fig_sii = plot_observable(
observable="sii",
state_ids=(
"al_gill_rho2p7_te5_ti5",
"al_clerouin_rho8p1_te10_ti10",
"al_clerouin_rho8p1_te10_ti2",
"be_wunsch_rho5p544_te13_ti13",
),
)
Pair distribution functions
The Be coordinates are ångström, while the carbon digitization uses Bohr.
fig_gii = plot_observable(
observable="gii",
state_ids=(
"be_wunsch_rho5p544_te13_ti13",
"c_starrett_rho20_te50_ti50",
),
)
save_figure(
fig_sii,
FIGURE_DIR / "ion_structure_library_sii",
close=False,
)
save_figure(
fig_gii,
FIGURE_DIR / "ion_structure_library_gii",
close=False,
)
plot_style.close()
if "agg" not in plt.get_backend().lower():
plt.show()
