""" 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 :class:`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 :math:`S_{ii}(k)` figure, panel 1 uses Gill *et al.*, Fig. 3 :cite:p:`GillEtAl2015`; panels 2 and 3 use Clérouin *et al.*, Fig. 1 :cite:p:`ClerouinEtAl2015`; and panel 4 uses Wünsch *et al.*, Fig. 2 :cite:p:`WunschEtAl2009`. 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 :cite:t:`StarrettSaumon2013,StarrettSaumon2014`, and the default finite-temperature jellium LFC follows :cite:t:`Chabrier1990`. See :doc:`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 # %% # 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()