.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "benchmarks/gen_benchmarks/plot_ion_structure_library.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_benchmarks_gen_benchmarks_plot_ion_structure_library.py: 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``. .. GENERATED FROM PYTHON SOURCE LINES 33-883 .. code-block:: Python 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 .. rst-class:: sphx-glr-script-out .. code-block:: none 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 .. GENERATED FROM PYTHON SOURCE LINES 884-889 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``. .. GENERATED FROM PYTHON SOURCE LINES 889-903 .. code-block:: Python 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", ), ) .. image-sg:: /benchmarks/gen_benchmarks/images/sphx_glr_plot_ion_structure_library_001.png :alt: Otter ion structure versus curated literature curves, Al: $\rho=2.7$ g cm$^{-3}$, $T_e=T_i=5$ eV, Al: $\rho=8.1$ g cm$^{-3}$, $T_e=T_i=10$ eV, Al: $\rho=8.1$ g cm$^{-3}$, $T_e=10$, $T_i=2$ eV, Be: $\rho=5.544$ g cm$^{-3}$, $T_e=T_i=13$ eV :srcset: /benchmarks/gen_benchmarks/images/sphx_glr_plot_ion_structure_library_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 904-908 Pair distribution functions --------------------------- The Be coordinates are ångström, while the carbon digitization uses Bohr. .. GENERATED FROM PYTHON SOURCE LINES 908-931 .. code-block:: Python 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() .. image-sg:: /benchmarks/gen_benchmarks/images/sphx_glr_plot_ion_structure_library_002.png :alt: Otter ion structure versus curated literature curves, Be: $\rho=5.544$ g cm$^{-3}$, $T_e=T_i=13$ eV, C: $\rho=20$ g cm$^{-3}$, $T_e=T_i=50$ eV :srcset: /benchmarks/gen_benchmarks/images/sphx_glr_plot_ion_structure_library_002.png :class: sphx-glr-single-img .. _sphx_glr_download_benchmarks_gen_benchmarks_plot_ion_structure_library.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_ion_structure_library.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_ion_structure_library.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_ion_structure_library.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_