.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "benchmarks/gen_benchmarks/plot_argha_roy_carbon_sii.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_argha_roy_carbon_sii.py: 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 = True`` Verify and load reviewed NPZ files generated by this same current-Otter workflow for all five displayed temperatures. ``USE_PRECOMPUTED_DATA = False`` Calculate 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 :cite:t:`StarrettSaumon2013,StarrettSaumon2014`; the default finite-temperature jellium local-field correction follows :cite:t:`Chabrier1990`. .. GENERATED FROM PYTHON SOURCE LINES 29-402 .. code-block:: Python 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}" ) .. rst-class:: sphx-glr-script-out .. code-block:: none 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 .. GENERATED FROM PYTHON SOURCE LINES 403-409 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. .. GENERATED FROM PYTHON SOURCE LINES 409-491 .. code-block:: Python 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() .. image-sg:: /benchmarks/gen_benchmarks/images/sphx_glr_plot_argha_roy_carbon_sii_001.png :alt: C: $\rho = 3.515\ \mathrm{g\,cc^{-1}}$ :srcset: /benchmarks/gen_benchmarks/images/sphx_glr_plot_argha_roy_carbon_sii_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none 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 .. _sphx_glr_download_benchmarks_gen_benchmarks_plot_argha_roy_carbon_sii.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_argha_roy_carbon_sii.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_argha_roy_carbon_sii.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_argha_roy_carbon_sii.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_