.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gen_examples/plot_carbon_ionization_levels.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_gen_examples_plot_carbon_ionization_levels.py: Carbon ionization and pressure-ionization levels ================================================ This example scans carbon at :math:`T_e=100` eV with Otter's orbital finite-temperature average atom. The same full-AA solutions provide * :math:`\bar{Z}=Z-Q_{\mathrm{ion}}(R_{\mathrm{WS}})`; * :math:`Z^*=n_e^0/n_i`; and * the carbon 1s, 2s, 2p, 3s, 3p, and 3d energies when localized, relative to Otter's local numerical continuum edge. The right-hand axis in each level panel shows how many electrons from that shell are assigned to the ionic density inside :math:`R_{\mathrm{WS}}`: .. math:: Q^{\rm ion}_{nl}(R_{\rm WS}) = 2(2l+1)f_{\rm FD}(E_{nl})M(E_{nl}) \int_0^{R_{\rm WS}} f_{\rm cut}(r)|P_{nl}(r)|^2\,dr. Summing these shell contributions gives :math:`Q_{\mathrm{ion}}(R_{\mathrm{WS}})`, which is subtracted from the nuclear charge to obtain :math:`\bar Z`. Solid curves with circle markers are level energies; dashed curves with cross markers are the shell contributions. Unlike the total orbital occupation :math:`\mathrm{OCC}_{nl}=2(2l+1)f_{\rm FD}`, this quantity contains the pressure-ionization weight and radial partition that enter :math:`\bar Z=Z-Q_{\rm ion}(R_{\rm WS})`, following :cite:t:`StarrettSaumon2013`. Mean ionization is not unique. Section 4.2 and Eq. (64) of :cite:t:`StarrettEtAl2019` discuss these two definitions: :math:`\bar Z` has an intuitive bound-state count but can jump at pressure ionization, whereas :math:`Z^*` is normally smooth but need not reproduce an intuitive chemical valence. Neither definition changes the self-consistent AA solution. A shallow level is plotted only while it lies below that edge and its threshold classification is resolved; Otter does not infer a precise disappearance density from such points. The bound/continuum construction and ionic-density partition follow :cite:t:`StarrettSaumon2014`; the negative-energy exterior matching used for shallow states follows the boundary-matching construction discussed by :cite:t:`StarrettEtAl2019`. For this scan the displayed edge is :math:`E_{\mathrm{cut}}=V_{\mathrm{eff}}(0.70R_{\mathrm{max}})`. This is the local continuum reference used consistently by the finite-domain orbital partition. It need not be zero because the finite numerical potential has not necessarily reached its asymptotic gauge value at that radius. Energies shown in the plot are :math:`E_{nl}-E_{\mathrm{cut}}`. For context, the ionization figure overlays the model-dependent :math:`Z^{\mathrm{free}}` curves digitized from Fig. 3(a) of :cite:t:`BethkenhagenEtAl2020`. They use different electron partitions and are not equivalent to either Otter :math:`\bar{Z}` or :math:`Z^*`. The default verifies and loads a checksummed 4096-point Otter scan. If the requested density grid contains new points, Otter reuses the accepted states and calculates only the missing densities. New files are staged under ``benchmarks/outputs`` and do not overwrite accepted data. States that fail the SCF or threshold-state checks are recorded as failures rather than plotted as physical results. Both figures are exported as PNG and PDF. .. GENERATED FROM PYTHON SOURCE LINES 65-1501 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /gen_examples/images/sphx_glr_plot_carbon_ionization_levels_001.png :alt: Carbon ionization at $T_e=100$ eV, Mean ionization, Electron chemical potential :srcset: /gen_examples/images/sphx_glr_plot_carbon_ionization_levels_001.png :class: sphx-glr-multi-img * .. image-sg:: /gen_examples/images/sphx_glr_plot_carbon_ionization_levels_002.png :alt: Carbon orbital levels at $T_e=100$ eV, Core level, Outer levels :srcset: /gen_examples/images/sphx_glr_plot_carbon_ionization_levels_002.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none Using checksummed Otter data from benchmarks/baselines/carbon_ionization_levels/C_Te100eV_density_scan.npz. C full-AA density scan at Te=100 eV rho [g/cc] Zbar Zstar mu [Ha] threshold 0.1 5.347648 5.096653 -20.077621 resolved 0.15 5.192316 4.964547 -18.681541 unresolved 0.2 5.069386 4.862945 -17.697859 resolved 0.25 4.969773 4.781682 -16.937347 resolved 0.3 4.997119 4.707238 -16.322653 unresolved 0.35 4.983817 4.658695 -15.791898 resolved 0.4 4.953475 4.608725 -15.338514 unresolved 0.45 4.919845 4.565573 -14.937964 unresolved 0.5 4.881229 4.527497 -14.579291 resolved 0.55 4.842348 4.492967 -14.254931 resolved 0.6 4.805973 4.461667 -13.958641 resolved 0.65 4.771880 4.433164 -13.685836 resolved 0.7 4.739721 4.406902 -13.433137 resolved 0.75 4.706743 4.405650 -13.178272 resolved 0.8 4.677790 4.382826 -12.958007 resolved 0.9 4.623762 4.342043 -12.555188 resolved 1 4.574209 4.305888 -12.194433 resolved 1.05 4.550861 4.289365 -12.027127 resolved 1.1 4.528354 4.273708 -11.867482 resolved 1.15 4.506666 4.258904 -11.714757 resolved 1.2 4.485746 4.244793 -11.568437 resolved 1.25 4.465504 4.231359 -11.427962 resolved 1.3 4.445940 4.218515 -11.292900 resolved 1.4 4.408632 4.194438 -11.037410 resolved 1.5 4.373850 4.172386 -10.799074 resolved 1.6 4.341300 4.152103 -10.575662 resolved 1.7 4.311034 4.133309 -10.365414 resolved 1.8 4.283116 4.116068 -10.166603 resolved 1.9 4.257428 4.099799 -9.978363 resolved 2 4.234006 4.084347 -9.799658 resolved 2.1 4.213707 4.070071 -9.629149 resolved 2.2 4.200870 4.056708 -9.466210 resolved 2.3 4.213833 4.043974 -9.310354 resolved 2.35 4.238909 4.037899 -9.234821 resolved 2.4 4.278830 4.032026 -9.160779 resolved 2.5 4.392770 4.020849 -9.016922 resolved 2.55 4.457192 4.015991 -8.946560 resolved 2.6 4.514166 4.012302 -8.876519 resolved 2.7 4.528657 4.006097 -8.739375 resolved 2.8 4.520863 3.989647 -8.617032 resolved 2.9 4.512592 3.981126 -8.491903 resolved 3 4.505313 3.972490 -8.371292 resolved 3.2 4.494753 3.956192 -8.141238 resolved 3.5 4.495470 3.934285 -7.820383 resolved 3.8 4.515756 3.914984 -7.524331 resolved 4 4.532322 3.903484 -7.338718 resolved 4.1 4.539562 3.897948 -7.249237 resolved 4.2 4.545598 3.892697 -7.161681 resolved 4.3 4.550236 3.887723 -7.075952 resolved 4.4 4.553388 3.882781 -6.992197 resolved 4.5 4.555127 3.878081 -6.910116 marginal 4.6 4.555503 3.873439 -6.829809 marginal 4.7 4.554647 3.869148 -6.750902 resolved 4.8 4.552740 3.864990 -6.673540 resolved 4.9 4.549996 3.860638 -6.597983 resolved 5 4.546859 3.857073 -6.523178 resolved 5.2 4.539173 3.848712 -6.379201 resolved 5.4 4.531304 3.842415 -6.238624 resolved 5.5 4.527508 3.839297 -6.170233 resolved 5.6 4.523781 3.836166 -6.103080 resolved 5.8 4.516520 3.830013 -5.972167 resolved 6 4.509556 3.824412 -5.845098 resolved 8 4.451877 3.783328 -4.749072 resolved 10 4.408982 3.761482 -3.871643 resolved 15 4.334834 3.747765 -2.196440 resolved 20 4.285781 3.761076 -0.924183 resolved 30 4.222816 3.838853 1.062906 resolved 50 4.165752 3.973376 3.936767 resolved 75 4.153607 4.109237 6.668707 resolved 100 4.176308 4.236630 8.991046 resolved 116 4.202605 4.305316 10.338352 resolved 136 4.245488 4.379977 11.915323 resolved 150 4.282169 4.425682 12.962277 resolved 185 4.415951 4.522451 15.424237 resolved 200 4.497026 4.557814 16.424448 resolved 216 4.599831 4.592163 17.461067 resolved 250 4.855347 4.655771 19.575658 resolved 270 5.016159 4.688937 20.773438 resolved 300 5.244634 4.732497 22.511204 resolved 320 5.382566 4.758596 23.636622 resolved 350 5.560044 4.794206 25.281480 resolved 370 5.658188 4.815847 26.351731 resolved 400 5.775171 4.845293 27.920256 resolved 450 5.903251 4.888753 30.451466 resolved | .. 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.electronic import FullExternalConfig, solve_full_only from otter.plotting import PALETTES, grid_figsize, save_figure, style_context # --------------------------------------------------------------------------- # User inputs # --------------------------------------------------------------------------- # Set this one switch in the script, then run the file directly. Incremental # reuse below ensures that only densities absent from the accepted scan/cache # are calculated. RECOMPUTE_WITH_OTTER = False if os.environ.get("OTTER_RECOMPUTE_CARBON_IONIZATION", "0") == "1": RECOMPUTE_WITH_OTTER = True RETRY_NONCONVERGED_POINTS = ( os.environ.get("OTTER_RETRY_NONCONVERGED_CARBON_IONIZATION", "0") == "1" ) ELEMENT = "C" TEMPERATURE_EV = 100.0 # This user-editable grid is deliberately denser where shallow-state branches # change rapidly. Points already present in the accepted scan are reused; # only newly requested densities are calculated. The grid diagnoses # finite-grid level disappearance rather than an exact ionization threshold. DENSITIES_G_CC = np.asarray( ( 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.90, 1.00, 1.05, 1.10, 1.15, 1.20, 1.25, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, 1.90, 2.00, 2.10, 2.20, 2.30, 2.35, 2.40, 2.50, 2.55, 2.60, 2.70, 2.80, 2.90, 3.00, 3.20, 3.50, 3.80, 4.00, 4.10, 4.20, 4.30, 4.40, 4.50, 4.60, 4.70, 4.80, 4.90, 5.00, 5.20, 5.40, 5.50, 5.60, 5.80, 6.00, 8.00, 10.0, 15.0, 20.0, 30.0, 50.0, 75.0, 100.0, 116.0, 136.0, 150.0, 185.0, 200.0, 216.0, 250.0, 270.0, 300.0, 320.0, 350.0, 370.0, 400.0, 450.0, ), dtype=float, ) # Two independent AA states, each with four continuum workers, use eight # explicit workers. This is faster than a sequential scan without the large # memory peak caused by nested state and continuum process pools. MAX_STATE_WORKERS = 2 CONTINUUM_WORKERS_PER_STATE = 2 # Incremental extension is the normal workflow: reuse every requested point # already present in the checksummed accepted scan, then calculate only new # densities. Set this to False only to force an independent full scan. REUSE_ACCEPTED_POINTS_WHEN_RECOMPUTING = ( os.environ.get("OTTER_REUSE_ACCEPTED_CARBON_IONIZATION", "1") == "1" ) AA_N_POINTS = 2**12 BOUND_ENERGY_CUT_MODE = "v_frac" BOUND_ENERGY_CUT_VALUE = 0.70 SCHEMA = "otter_carbon_ionization_levels_v3" LEGACY_BASELINE_SCHEMAS = { "otter_carbon_ionization_levels_v1", "otter_carbon_ionization_levels_v2", } class DensityGridMismatchError(RuntimeError): """The requested density grid differs from the accepted archive.""" HARTREE_TO_EV = 27.211386245988 ORBITAL_LETTERS = ("s", "p", "d", "f", "g", "h") DISPLAYED_SHELLS = ("1s", "2s", "2p", "3s", "3p", "3d") def _repository_root() -> Path: """Locate the source tree from either the gallery source or generated copy.""" 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 / "src" / "otter").is_dir() and ( candidate / "pyproject.toml" ).is_file(): return candidate raise FileNotFoundError("Cannot locate the Otter repository root.") ROOT = _repository_root() BASELINE_DIR = ROOT / "benchmarks" / "baselines" / "carbon_ionization_levels" BASELINE_PATH = BASELINE_DIR / "C_Te100eV_density_scan.npz" BASELINE_MANIFEST = BASELINE_DIR / "manifest.json" OUTPUT_DIR = ROOT / "benchmarks" / "outputs" / "carbon_ionization_levels" CANDIDATE_PATH = OUTPUT_DIR / "C_Te100eV_density_scan.npz" POINT_CACHE_DIR = OUTPUT_DIR / "point_cache" FIGURE_DIR = OUTPUT_DIR / "figures" REFERENCE_DIR = ( ROOT / "benchmarks" / "reference_data" / "bethkenhagen_et_al_2020_carbon_ionization" ) REFERENCE_MANIFEST = REFERENCE_DIR / "manifest.json" def _sha256(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_bethkenhagen_reference() -> dict[str, tuple[np.ndarray, np.ndarray]]: """Load checksummed Fig. 3(a) ``Z_free`` model curves for context.""" manifest = json.loads(REFERENCE_MANIFEST.read_text(encoding="utf-8")) if ( manifest.get("schema_version") != "otter_reference_manifest_v1" or manifest.get("dataset_id") != "bethkenhagen_et_al_2020_carbon_ionization" or str(manifest["source"]["doi"]) != "10.1103/PhysRevResearch.2.023260" ): raise ValueError("Unexpected Bethkenhagen et al. reference manifest.") curves: dict[str, tuple[np.ndarray, np.ndarray]] = {} for record in manifest["files"]: path = REFERENCE_DIR / str(record["file"]) if _sha256(path) != str(record["sha256"]): raise RuntimeError(f"Reference checksum mismatch: {path.name}.") values = np.loadtxt(path, delimiter=",", comments="#") if ( values.ndim != 2 or values.shape[1] != 2 or np.any(~np.isfinite(values)) or np.any(values[:, 0] <= 0.0) or np.any(np.diff(values[:, 0]) <= 0.0) ): raise ValueError(f"Malformed reference curve: {path.name}.") curves[str(record["label"])] = (values[:, 0], values[:, 1]) return curves def _level_label(l_value: int, radial_index: int) -> str: principal_n = int(radial_index) + int(l_value) if 0 <= int(l_value) < len(ORBITAL_LETTERS): return f"{principal_n}{ORBITAL_LETTERS[int(l_value)]}" return f"n={principal_n},l={int(l_value)}" def _configuration(rho_g_cc: float) -> FullExternalConfig: return FullExternalConfig( element=ELEMENT, temperature_ev=float(TEMPERATURE_EV), rho_g_cc=float(rho_g_cc), run_mode="full", # Retain headroom on the steep high-density ionization branch. The # convergence criterion itself is unchanged. stage2_max_iter=180, cont_n_jobs=int(CONTINUUM_WORKERS_PER_STATE), cont_shards=int(2 * CONTINUUM_WORKERS_PER_STATE), # Match a shallow negative-energy orbital at the common outer SCF # boundary, with no separate enlarged bound-only box. This optional # numerical refinement is motivated by the exterior matching in # Starrett et al. (2019), Eqs. (21)-(22), but is not identical to the # ion-sphere-boundary implementation in that work. bound_zero_tail_refine=True, bound_zero_tail_max_binding_ha=1.0e-2, bound_zero_tail_scan_points=64, bound_zero_tail_l_max=1, bound_zero_tail_edge_rel_tol=0.1, ) def _finite_levels( result: dict[str, Any], *, continuum_edge_ha: float, ) -> dict[str, float]: """Return levels below the same continuum edge used by the AA density.""" energies = np.asarray( result.get("bound_energy_ha", np.empty((0, 0))), dtype=float, ) l_values = np.asarray( result.get("bound_l_list", np.arange(energies.shape[0])), dtype=int, ) levels: dict[str, float] = {} for l_index, l_value in enumerate(l_values): for state_index in range(energies.shape[1]): energy_ha = float(energies[l_index, state_index]) if np.isfinite(energy_ha) and energy_ha < continuum_edge_ha: label = _level_label(int(l_value), int(state_index + 1)) levels[label] = ( energy_ha - float(continuum_edge_ha) ) * HARTREE_TO_EV return levels def _finite_level_ion_charges( result: dict[str, Any], ) -> dict[str, float]: r"""Return shell contributions to ``Q_ion(R_WS)`` in electrons. ``bound_q_ion_ws`` is assembled from the exact final orbitals used by the electronic solver and includes ``2(2l+1) f_FD M(E) f_cut(r)``. Keeping this reduction in the producer avoids reconstructing a spatially weighted quantity from energies alone. """ charges = np.asarray( result.get("bound_q_ion_ws", np.empty((0, 0))), dtype=float, ) l_values = np.asarray( result.get("bound_l_list", np.arange(charges.shape[0])), dtype=int, ) if charges.ndim != 2 or l_values.shape != (charges.shape[0],): raise RuntimeError("Malformed shell-resolved Q_ion table.") total = float(np.nansum(charges)) q_ion_ws = float(result.get("q_ion_ws", np.nan)) if ( not np.isfinite(q_ion_ws) or abs(total - q_ion_ws) > 2.0e-8 * max(1.0, abs(q_ion_ws)) ): raise RuntimeError( "Shell-resolved Q_ion does not close the total WS ionic charge: " f"sum={total:.12g}, total={q_ion_ws:.12g}." ) shell_charges: dict[str, float] = {} for l_index, l_value in enumerate(l_values): for state_index in range(charges.shape[1]): charge = float(charges[l_index, state_index]) if not np.isfinite(charge): continue if charge < -1.0e-12: raise RuntimeError("A shell-resolved ionic charge is negative.") label = _level_label(int(l_value), int(state_index + 1)) shell_charges[label] = max(charge, 0.0) return shell_charges def _solve_density(rho_g_cc: float) -> dict[str, Any]: """Compute and reduce one independent full-AA state.""" started = time.perf_counter() result = solve_full_only(_configuration(float(rho_g_cc))) elapsed_s = time.perf_counter() - started stage2_history = list(result.get("history", ())) stage2_error = ( float(stage2_history[-1].get("err", np.nan)) if stage2_history else np.nan ) if not bool(result.get("stage2_converged", False)): bound_charge = np.asarray( [entry.get("charge_bound", np.nan) for entry in stage2_history], dtype=float, ) finite_bound_charge = bound_charge[np.isfinite(bound_charge)] return { "record_type": "stage2_nonconvergence", "rho_g_cc": float(rho_g_cc), "elapsed_s": float(elapsed_s), "stage": "full_aa_stage2", "stage2_converged": False, "stage2_error": stage2_error, "stage2_iters": int( result.get("stage2_iters", len(stage2_history)) ), "bound_charge_min_e": ( float(np.min(finite_bound_charge)) if finite_bound_charge.size else np.nan ), "bound_charge_max_e": ( float(np.max(finite_bound_charge)) if finite_bound_charge.size else np.nan ), "message": ( f"C rho={rho_g_cc:g} g/cc: production full-AA stage 2 " "did not reach a physical fixed point." ), "point_source": "fresh_otter_solve", } threshold_status = str(result.get("threshold_state_status", "none")).lower() if not np.isfinite(stage2_error): raise RuntimeError( f"C rho={rho_g_cc:g} g/cc: final stage-2 error is not finite." ) meta = dict(result.get("meta", {})) continuum_edge_ha = float(meta["bound_energy_cut_ha"]) n_i_bohr3 = float(meta["n_i_bohr3"]) n0_bohr3 = float(meta["n0_final_bohr3"]) zbar = float(result["zbar"]) zstar = n0_bohr3 / n_i_bohr3 if not np.all( np.isfinite( ( continuum_edge_ha, n_i_bohr3, n0_bohr3, zbar, zstar, float(result["mu"]), ) ) ): raise RuntimeError(f"C rho={rho_g_cc:g} g/cc produced non-finite output.") raw_levels = _finite_levels( result, continuum_edge_ha=continuum_edge_ha, ) level_q_ion_ws = _finite_level_ion_charges(result) levels = raw_levels if threshold_status in {"marginal", "unresolved"}: # The flagged state concerns the shallow outer branch. Retain the # deeply bound 1s diagnostic, but never turn an unreliable shallow # eigenvalue into a pressure-ionization datum. levels = { label: energy for label, energy in levels.items() if label == "1s" } return { "rho_g_cc": float(rho_g_cc), "zbar": zbar, "zstar": zstar, "mu_ha": float(result["mu"]), "continuum_edge_ha": continuum_edge_ha, "threshold_status": threshold_status, "threshold_representation": str( result.get("threshold_state_representation", "none") ), "levels_ev": levels, # Keep every shell charge that actually entered Q_ion, including a # shallow branch whose energy is hidden by a marginal/unresolved # threshold flag. This makes the smooth M(E) partition auditable. "level_q_ion_ws": level_q_ion_ws, "elapsed_s": float(elapsed_s), "stage2_converged": True, "stage2_error": stage2_error, "stage2_iters": int(result.get("stage2_iters", len(stage2_history))), "point_source": "fresh_otter_solve", } def _point_cache_path(rho_g_cc: float) -> Path: token = f"{float(rho_g_cc):012.6f}".replace(".", "p") return POINT_CACHE_DIR / f"C_Te100eV_rho{token}.json" def _point_failure_cache_path(rho_g_cc: float) -> Path: return _point_cache_path(rho_g_cc).with_suffix(".failure.json") def _cache_payload(record: dict[str, Any]) -> dict[str, Any]: """Return method metadata shared by converged and failed point caches.""" return { "schema_version": SCHEMA, "temperature_ev": TEMPERATURE_EV, "aa_n_points": AA_N_POINTS, "continuum_workers": CONTINUUM_WORKERS_PER_STATE, "bound_occ_mode": "fd", "bound_energy_cut_mode": BOUND_ENERGY_CUT_MODE, "bound_energy_cut_value": BOUND_ENERGY_CUT_VALUE, "record": record, } def _cache_metadata_matches(payload: dict[str, Any]) -> bool: """Check the immutable scientific controls represented by a point cache.""" return bool( payload.get("schema_version") == SCHEMA and np.isclose( float(payload.get("temperature_ev", np.nan)), TEMPERATURE_EV, ) and int(payload.get("aa_n_points", -1)) == AA_N_POINTS # The first 61 v3 caches predate this explicit field; v3 always used # the production FD state sum, so a missing value is unambiguous. and payload.get("bound_occ_mode", "fd") == "fd" and payload.get("bound_energy_cut_mode") == BOUND_ENERGY_CUT_MODE and np.isclose( float(payload.get("bound_energy_cut_value", np.nan)), BOUND_ENERGY_CUT_VALUE, ) ) def _save_point_cache(row: dict[str, Any]) -> None: """Save one completed AA point so interrupted scans can resume safely.""" POINT_CACHE_DIR.mkdir(parents=True, exist_ok=True) payload = _cache_payload(row) # Preserve the row key used by existing v3 caches. payload["row"] = payload.pop("record") path = _point_cache_path(float(row["rho_g_cc"])) temporary = path.with_suffix(".tmp") temporary.write_text( json.dumps(payload, indent=2, sort_keys=True) + "\n", encoding="utf-8", ) temporary.replace(path) def _load_point_cache(rho_g_cc: float) -> dict[str, Any] | None: path = _point_cache_path(rho_g_cc) if not path.is_file(): return None payload = json.loads(path.read_text(encoding="utf-8")) if not _cache_metadata_matches(payload): return None row = dict(payload["row"]) if not np.isclose(float(row["rho_g_cc"]), float(rho_g_cc)): return None row["levels_ev"] = { str(key): float(value) for key, value in dict(row.get("levels_ev", {})).items() } row["level_q_ion_ws"] = { str(key): float(value) for key, value in dict(row.get("level_q_ion_ws", {})).items() } if not row["level_q_ion_ws"]: return None return row def _save_point_failure(record: dict[str, Any]) -> None: """Checkpoint a scientifically unusable state without storing AA output.""" POINT_CACHE_DIR.mkdir(parents=True, exist_ok=True) payload = _cache_payload(record) path = _point_failure_cache_path(float(record["rho_g_cc"])) temporary = path.with_suffix(".tmp") temporary.write_text( json.dumps(payload, indent=2, sort_keys=True) + "\n", encoding="utf-8", ) temporary.replace(path) def _load_point_failure(rho_g_cc: float) -> dict[str, Any] | None: """Load a same-method failed-point audit unless an explicit retry is set.""" if RETRY_NONCONVERGED_POINTS: return None path = _point_failure_cache_path(rho_g_cc) if not path.is_file(): return None payload = json.loads(path.read_text(encoding="utf-8")) if not _cache_metadata_matches(payload): return None record = dict(payload.get("record", {})) if ( record.get("record_type") != "stage2_nonconvergence" or not np.isclose(float(record.get("rho_g_cc", np.nan)), rho_g_cc) or bool(record.get("stage2_converged", True)) ): return None return record def _accepted_seed_rows(densities: np.ndarray) -> list[dict[str, Any]]: """Reuse accepted points while computing only newly requested states.""" if not REUSE_ACCEPTED_POINTS_WHEN_RECOMPUTING: return [] if not BASELINE_PATH.is_file() or not BASELINE_MANIFEST.is_file(): return [] manifest = json.loads(BASELINE_MANIFEST.read_text(encoding="utf-8")) if ( manifest.get("schema_version") != "otter_example_manifest_v1" or manifest.get("example_id") != "carbon_ionization_levels" or not str(manifest.get("status", "")).startswith("accepted") or dict(manifest.get("state", {})).get("data_sha256") != _sha256(BASELINE_PATH) ): raise RuntimeError("Cannot reuse an unreviewed carbon-ionization state.") with np.load(BASELINE_PATH, allow_pickle=False) as archive: old = {key: np.asarray(archive[key]) for key in archive.files} # Never mix a coarse archive into a production 4096-point regeneration, # and never reconstruct spatial Q_ion data from energies. if ( int(np.asarray(old.get("aa_n_points", -1)).item()) != AA_N_POINTS or "level_q_ion_ws" not in old or "level_q_ion_ws_is_available" not in old ): return [] old = _upgrade_accepted_state(old, manifest) if not bool(np.asarray(old["level_q_ion_ws_available"]).item()): return [] old_rho = np.asarray(old["rho_g_cc"], dtype=float) requested = {float(value) for value in np.asarray(densities, dtype=float)} labels = tuple(np.asarray(old["level_labels"], dtype=str)) rows: list[dict[str, Any]] = [] for index, rho_g_cc in enumerate(old_rho): if float(rho_g_cc) not in requested: continue levels = { label: float(old["level_energy_relative_edge_ev"][index, level_index]) for level_index, label in enumerate(labels) if bool(old["level_is_bound"][index, level_index]) } level_q_ion_ws = { label: float(old["level_q_ion_ws"][index, level_index]) for level_index, label in enumerate(labels) if bool(old["level_q_ion_ws_is_available"][index, level_index]) } rows.append( { "rho_g_cc": float(rho_g_cc), "zbar": float(old["zbar"][index]), "zstar": float(old["zstar"][index]), "mu_ha": float(old["mu_ha"][index]), "continuum_edge_ha": float(old["continuum_edge_ha"][index]), "threshold_status": str(old["threshold_status"][index]), "threshold_representation": str( old["threshold_representation"][index] ), "levels_ev": levels, "level_q_ion_ws": level_q_ion_ws, "elapsed_s": float(old["elapsed_s"][index]), "stage2_converged": True, "stage2_error": float(old["stage2_error"][index]), "stage2_iters": int(old["stage2_iters"][index]), "point_source": "accepted_baseline_seed", } ) return rows def _compute_scan() -> dict[str, np.ndarray]: densities = np.unique(np.asarray(DENSITIES_G_CC, dtype=float)) if densities.size < 2 or np.any(densities <= 0.0): raise ValueError("DENSITIES_G_CC must contain at least two positive values.") rows = _accepted_seed_rows(densities) failures: list[dict[str, Any]] = [] seeded_rho = {float(row["rho_g_cc"]) for row in rows} pending: list[float] = [] for rho in densities: if float(rho) in seeded_rho: print(f"[accepted baseline] C, rho={float(rho):g} g/cc") continue cached = _load_point_cache(float(rho)) if cached is not None: rows.append(cached) print(f"[cached] C, rho={float(rho):g} g/cc") continue failed = _load_point_failure(float(rho)) if failed is not None: failures.append(failed) print(f"[cached nonconverged] C, rho={float(rho):g} g/cc") continue pending.append(float(rho)) worker_count = min(max(int(MAX_STATE_WORKERS), 1), int(densities.size)) if worker_count == 1: for rho in pending: record = _solve_density(float(rho)) if bool(record.get("stage2_converged", False)): _save_point_cache(record) rows.append(record) print(f"[computed] C, rho={rho:g} g/cc") else: _save_point_failure(record) failures.append(record) print(f"[not converged] C, rho={rho:g} g/cc") else: with ProcessPoolExecutor(max_workers=worker_count) as executor: futures = { executor.submit(_solve_density, rho): rho for rho in pending } for future in as_completed(futures): rho = futures[future] record = future.result() if bool(record.get("stage2_converged", False)): _save_point_cache(record) rows.append(record) print(f"[computed] C, rho={rho:g} g/cc") else: _save_point_failure(record) failures.append(record) print(f"[not converged] C, rho={rho:g} g/cc") if len(rows) + len(failures) != densities.size: raise RuntimeError( f"Expected {densities.size} density attempts, obtained " f"{len(rows)} converged and {len(failures)} failed records." ) rows.sort(key=lambda item: float(item["rho_g_cc"])) failures.sort(key=lambda item: float(item["rho_g_cc"])) # Track the K/L branches and the short-lived localized 3s branch. Higher # finite-box Rydberg roots remain internal diagnostics rather than ionic # levels in this gallery. labels = sorted( { label for row in rows for label in ( set(row["levels_ev"]) | set(row["level_q_ion_ws"]) ) if label in DISPLAYED_SHELLS }, key=lambda label: ( int(label[:-1]) if label[-1:] in ORBITAL_LETTERS else 99, ORBITAL_LETTERS.index(label[-1]) if label[-1:] in ORBITAL_LETTERS else 99, ), ) level_energy_ev = np.zeros((len(rows), len(labels)), dtype=float) level_is_bound = np.zeros_like(level_energy_ev, dtype=bool) level_q_ion_ws = np.zeros_like(level_energy_ev, dtype=float) level_q_ion_ws_is_available = np.zeros_like(level_energy_ev, dtype=bool) for density_index, row in enumerate(rows): for level_index, label in enumerate(labels): if label in row["levels_ev"]: level_energy_ev[density_index, level_index] = float( row["levels_ev"][label] ) level_is_bound[density_index, level_index] = True if label in row["level_q_ion_ws"]: level_q_ion_ws[density_index, level_index] = float( row["level_q_ion_ws"][label] ) level_q_ion_ws_is_available[density_index, level_index] = True return { "schema_version": np.asarray(SCHEMA), "element_symbol": np.asarray(ELEMENT), "temperature_ev": np.asarray(TEMPERATURE_EV), "rho_g_cc": np.asarray([row["rho_g_cc"] for row in rows]), "zbar": np.asarray([row["zbar"] for row in rows]), "zstar": np.asarray([row["zstar"] for row in rows]), "mu_ha": np.asarray([row["mu_ha"] for row in rows]), "continuum_edge_ha": np.asarray( [row["continuum_edge_ha"] for row in rows] ), "threshold_status": np.asarray( [row["threshold_status"] for row in rows] ), "threshold_representation": np.asarray( [row["threshold_representation"] for row in rows] ), "level_labels": np.asarray(labels), "level_energy_relative_edge_ev": level_energy_ev, "level_is_bound": level_is_bound, "level_q_ion_ws": level_q_ion_ws, "level_q_ion_ws_is_available": level_q_ion_ws_is_available, "level_q_ion_ws_available": np.asarray(True), "elapsed_s": np.asarray([row["elapsed_s"] for row in rows]), "stage2_converged": np.asarray( [row["stage2_converged"] for row in rows], dtype=bool, ), "stage2_error": np.asarray([row["stage2_error"] for row in rows]), "stage2_iters": np.asarray([row["stage2_iters"] for row in rows]), "point_source": np.asarray([row["point_source"] for row in rows]), "failed_rho_g_cc": np.asarray( [record["rho_g_cc"] for record in failures], dtype=float ), "failed_stage": np.asarray( [record["stage"] for record in failures], dtype=str ), "failed_message": np.asarray( [record["message"] for record in failures], dtype=str ), "failed_stage2_iters": np.asarray( [record["stage2_iters"] for record in failures], dtype=int ), "failed_stage2_error": np.asarray( [record["stage2_error"] for record in failures], dtype=float ), "failed_bound_charge_min_e": np.asarray( [record["bound_charge_min_e"] for record in failures], dtype=float ), "failed_bound_charge_max_e": np.asarray( [record["bound_charge_max_e"] for record in failures], dtype=float ), "failed_point_source": np.asarray( [record["point_source"] for record in failures], dtype=str ), "bound_occ_mode": np.asarray("fd"), "bound_rmax_mult": np.asarray("none"), "bound_zero_tail_refine": np.asarray(True), "bound_energy_cut_mode": np.asarray(BOUND_ENERGY_CUT_MODE), "bound_energy_cut_value": np.asarray(BOUND_ENERGY_CUT_VALUE), "b3_tail_model": np.asarray("full"), "aa_n_points": np.asarray(AA_N_POINTS), } def _upgrade_accepted_state( state: dict[str, np.ndarray], manifest: dict[str, Any], ) -> dict[str, np.ndarray]: """Expose a legacy accepted archive through the current plotting schema. v1 predates per-point SCF errors; v1 and v2 both predate shell-resolved ``Q_ion``. Missing scientific data are marked unavailable and are never reconstructed from energies or OCC values. """ schema = str(state.get("schema_version", np.asarray("")).item()) if schema == SCHEMA: upgraded = dict(state) upgraded.setdefault("failed_rho_g_cc", np.asarray([], dtype=float)) upgraded.setdefault("failed_stage", np.asarray([], dtype=str)) upgraded.setdefault("failed_message", np.asarray([], dtype=str)) upgraded.setdefault("failed_stage2_iters", np.asarray([], dtype=int)) upgraded.setdefault("failed_stage2_error", np.asarray([], dtype=float)) upgraded.setdefault( "failed_bound_charge_min_e", np.asarray([], dtype=float) ) upgraded.setdefault( "failed_bound_charge_max_e", np.asarray([], dtype=float) ) upgraded.setdefault("failed_point_source", np.asarray([], dtype=str)) return upgraded if schema not in LEGACY_BASELINE_SCHEMAS: raise ValueError(f"Unsupported legacy carbon schema {schema!r}.") rho = np.asarray(state.get("rho_g_cc", ()), dtype=float) upgraded = dict(state) if schema == "otter_carbon_ionization_levels_v1": converged_count = int( dict(manifest.get("scientific_audit", {})).get( "full_aa_converged_states", -1, ) ) if rho.ndim != 1 or converged_count != rho.size: raise RuntimeError( "The accepted v1 carbon baseline lacks a complete " "convergence audit." ) upgraded.update({ "stage2_converged": np.ones(rho.shape, dtype=bool), "stage2_error": np.full(rho.shape, np.nan, dtype=float), "stage2_iters": np.full(rho.shape, -1, dtype=int), "point_source": np.full( rho.shape, "accepted_v1_baseline", dtype=" None: required = { "schema_version", "element_symbol", "temperature_ev", "rho_g_cc", "zbar", "zstar", "mu_ha", "continuum_edge_ha", "level_labels", "level_energy_relative_edge_ev", "level_is_bound", "level_q_ion_ws", "level_q_ion_ws_is_available", "level_q_ion_ws_available", "bound_energy_cut_mode", "bound_energy_cut_value", "aa_n_points", "stage2_converged", "stage2_error", "point_source", "failed_rho_g_cc", "failed_stage", "failed_message", "failed_stage2_iters", "failed_stage2_error", "failed_bound_charge_min_e", "failed_bound_charge_max_e", "failed_point_source", } missing = required.difference(state) if missing: raise KeyError(f"Carbon ionization state is missing {sorted(missing)}.") if str(state["schema_version"].item()) != SCHEMA: raise ValueError("Unsupported carbon ionization state schema.") if str(state["element_symbol"].item()) != ELEMENT: raise ValueError("Unexpected element in carbon ionization state.") if not np.isclose(float(state["temperature_ev"]), TEMPERATURE_EV): raise ValueError("Unexpected temperature in carbon ionization state.") if str(state["bound_energy_cut_mode"].item()) != BOUND_ENERGY_CUT_MODE: raise ValueError("Unexpected bound/continuum edge convention.") if not np.isclose( float(state["bound_energy_cut_value"]), BOUND_ENERGY_CUT_VALUE, ): raise ValueError("Unexpected bound/continuum edge parameter.") rho = np.asarray(state["rho_g_cc"], dtype=float) if rho.ndim != 1 or rho.size < 2 or np.any(np.diff(rho) <= 0.0): raise ValueError("Density grid must be a strictly increasing vector.") for key in ("zbar", "zstar", "mu_ha"): values = np.asarray(state[key], dtype=float) if values.shape != rho.shape or not np.all(np.isfinite(values)): raise ValueError(f"{key} must be finite on the density grid.") continuum_edge = np.asarray(state["continuum_edge_ha"], dtype=float) if continuum_edge.shape != rho.shape or not np.all(np.isfinite(continuum_edge)): raise ValueError("continuum_edge_ha must be finite on the density grid.") level_energy = np.asarray( state["level_energy_relative_edge_ev"], dtype=float, ) level_mask = np.asarray(state["level_is_bound"], dtype=bool) labels = np.asarray(state["level_labels"]) if level_energy.shape != level_mask.shape or level_energy.shape != ( rho.size, labels.size, ): raise ValueError("Bound-level arrays are not aligned.") if not np.all(np.isfinite(level_energy)): raise ValueError("Bound-level storage must be finite; mask absent levels.") if np.any(level_energy[level_mask] >= 0.0): raise ValueError("A stored bound level lies above the continuum edge.") level_q_ion_ws = np.asarray(state["level_q_ion_ws"], dtype=float) q_ion_mask = np.asarray( state["level_q_ion_ws_is_available"], dtype=bool, ) q_ion_available = bool(np.asarray(state["level_q_ion_ws_available"]).item()) if level_q_ion_ws.shape != level_energy.shape or q_ion_mask.shape != ( level_energy.shape ): raise ValueError("Shell Q_ion arrays are not aligned with the levels.") if not np.all(np.isfinite(level_q_ion_ws)): raise ValueError("Shell Q_ion storage must be finite; mask absent values.") if np.any(level_q_ion_ws[q_ion_mask] < 0.0): raise ValueError("Shell Q_ion contributions must be non-negative.") if np.any(level_q_ion_ws[~q_ion_mask] != 0.0): raise ValueError("Unavailable shell Q_ion entries must use zero storage.") if q_ion_available != bool(np.any(q_ion_mask)): raise ValueError("Shell Q_ion availability metadata is inconsistent.") n_points = int(np.asarray(state["aa_n_points"]).item()) if q_ion_available and n_points != AA_N_POINTS: raise ValueError( f"Production shell Q_ion data require {AA_N_POINTS} radial points." ) statuses = np.asarray(state.get("threshold_status", ()), dtype=str) if statuses.shape != rho.shape: raise ValueError("threshold_status must align with the density grid.") if not set(statuses).issubset( {"none", "resolved", "marginal", "unresolved"} ): raise ValueError("Unknown threshold-state classification.") converged = np.asarray(state["stage2_converged"], dtype=bool) errors = np.asarray(state["stage2_error"], dtype=float) point_source = np.asarray(state["point_source"], dtype=str) if converged.shape != rho.shape or not np.all(converged): raise ValueError("Every full-AA point must have a reviewed convergence flag.") if errors.shape != rho.shape or point_source.shape != rho.shape: raise ValueError("Convergence audit arrays must align with density.") fresh = point_source == "fresh_otter_solve" if np.any(~np.isfinite(errors[fresh])): raise ValueError("Fresh full-AA points require a finite stage-2 error.") if np.any(fresh) and n_points == AA_N_POINTS and not q_ion_available: raise ValueError("Fresh states require shell-resolved Q_ion diagnostics.") failed_rho = np.asarray(state["failed_rho_g_cc"], dtype=float) failed_size = failed_rho.size if ( failed_rho.ndim != 1 or np.any(~np.isfinite(failed_rho)) or np.any(failed_rho <= 0.0) or np.any(np.diff(failed_rho) <= 0.0) ): raise ValueError("Failed-point densities must be positive and increasing.") failed_fields = ( "failed_stage", "failed_message", "failed_stage2_iters", "failed_stage2_error", "failed_bound_charge_min_e", "failed_bound_charge_max_e", "failed_point_source", ) if any(np.asarray(state[key]).shape != (failed_size,) for key in failed_fields): raise ValueError("Failed-point audit arrays must be aligned.") if failed_size and np.intersect1d(rho, failed_rho).size: raise ValueError("A density cannot be both converged and failed.") if np.any(np.asarray(state["failed_stage2_iters"], dtype=int) <= 0): raise ValueError("Failed-point iteration counts must be positive.") if np.any(np.asarray(state["failed_stage"], dtype=str) == "") or np.any( np.asarray(state["failed_message"], dtype=str) == "" ): raise ValueError("Failed-point stage and message must be recorded.") def _save_candidate(state: dict[str, np.ndarray]) -> Path: _validate_state(state) point_source = np.asarray(state["point_source"], dtype=str) fresh = point_source == "fresh_otter_solve" fresh_errors = np.asarray(state["stage2_error"], dtype=float)[fresh] CANDIDATE_PATH.parent.mkdir(parents=True, exist_ok=True) temporary = CANDIDATE_PATH.with_suffix(".tmp") with temporary.open("wb") as stream: np.savez_compressed(stream, **state) temporary.replace(CANDIDATE_PATH) manifest = { "schema_version": "otter_gallery_manifest_v1", "example_id": "carbon_ionization_levels", "status": "candidate_not_accepted", "state": { "data_file": CANDIDATE_PATH.name, "data_sha256": _sha256(CANDIDATE_PATH), }, "configuration": { "element": ELEMENT, "temperature_ev": TEMPERATURE_EV, "densities_g_cc": np.asarray(DENSITIES_G_CC, dtype=float).tolist(), "reuse_accepted_points_when_recomputing": bool( REUSE_ACCEPTED_POINTS_WHEN_RECOMPUTING ), "bound_occ_mode": "fd", "bound_rmax_mult": None, "bound_zero_tail_refine": True, "bound_energy_cut_mode": BOUND_ENERGY_CUT_MODE, "bound_energy_cut_value": BOUND_ENERGY_CUT_VALUE, "b3_tail_model": "full", "aa_n_points": AA_N_POINTS, "radial_resolution_policy": "production_2**12_no_coarse_seed_reuse", }, "method_references": [ { "citation_key": "StarrettSaumon2013", "doi": "10.1103/PhysRevE.87.013104", "scope": "M(E) pressure-ionization weight and radial cutoff", }, { "citation_key": "StarrettSaumon2014", "doi": "10.1016/j.hedp.2013.12.001", "scope": "average-atom and ionic-density partition", }, ], "shell_charge_diagnostic": { "field": "level_q_ion_ws", "units": "electrons", "definition": ( "2(2l+1) f_FD(E_nl) M(E_nl) integral_0^Rws " "f_cut(r) |P_nl(r)|^2 dr" ), }, "scientific_audit": { "requested_states": int( np.asarray(DENSITIES_G_CC, dtype=float).size ), "stage2_converged_states": int( np.count_nonzero(state["stage2_converged"]) ), "stage2_nonconverged_states": int( np.asarray(state["failed_rho_g_cc"]).size ), "stage2_nonconverged_densities_g_cc": np.asarray( state["failed_rho_g_cc"], dtype=float ).tolist(), "nonconverged_policy": ( "Failed states are audit metadata only and are excluded from " "all electronic and shell-level arrays; diagnostic fd_m " "solutions are never substituted for production fd states." ), "fresh_states": int( np.count_nonzero(fresh) ), "accepted_seed_states": int( np.count_nonzero(point_source == "accepted_baseline_seed") ), "fresh_max_stage2_error": ( float(np.max(fresh_errors)) if fresh_errors.size else None ), }, } manifest_path = CANDIDATE_PATH.with_suffix(".manifest.json") manifest_path.write_text( json.dumps(manifest, indent=2, sort_keys=True) + "\n", encoding="utf-8", ) return CANDIDATE_PATH def _load_precomputed() -> dict[str, np.ndarray]: if not BASELINE_PATH.is_file() or not BASELINE_MANIFEST.is_file(): raise FileNotFoundError( "The checksummed carbon-ionization gallery state is not installed. " "Run with RECOMPUTE_WITH_OTTER=True to calculate it." ) manifest = json.loads(BASELINE_MANIFEST.read_text(encoding="utf-8")) if ( manifest.get("schema_version") != "otter_example_manifest_v1" or manifest.get("example_id") != "carbon_ionization_levels" or not str(manifest.get("status", "")).startswith("accepted") ): raise ValueError("The installed carbon-ionization manifest is not accepted.") state_entry = dict(manifest.get("state", {})) if state_entry.get("data_file") != BASELINE_PATH.name: raise ValueError("Carbon ionization manifest names the wrong data file.") if state_entry.get("data_sha256") != _sha256(BASELINE_PATH): raise RuntimeError("Carbon ionization baseline checksum mismatch.") with np.load(BASELINE_PATH, allow_pickle=False) as archive: state = {key: np.asarray(archive[key]) for key in archive.files} state = _upgrade_accepted_state(state, manifest) requested_rho = np.asarray(DENSITIES_G_CC, dtype=float) stored_rho = np.asarray(state.get("rho_g_cc", ()), dtype=float) if stored_rho.shape != requested_rho.shape or not np.allclose( stored_rho, requested_rho, rtol=0.0, atol=0.0, ): raise DensityGridMismatchError( "The requested density grid differs from the accepted " f"{stored_rho.size}-state grid." ) _validate_state(state) return state def _print_state_table(state: dict[str, np.ndarray]) -> None: print("\nC full-AA density scan at Te=100 eV") print( f"{'rho [g/cc]':>12} {'Zbar':>10} {'Zstar':>10} " f"{'mu [Ha]':>12} {'threshold':>12}" ) statuses = np.asarray(state.get("threshold_status", []), dtype=str) for index, rho in enumerate(np.asarray(state["rho_g_cc"], dtype=float)): status = statuses[index] if statuses.size else "not recorded" print( f"{rho:12.5g} {float(state['zbar'][index]):10.6f} " f"{float(state['zstar'][index]):10.6f} " f"{float(state['mu_ha'][index]):12.6f} {status:>12}" ) failed_rho = np.asarray(state["failed_rho_g_cc"], dtype=float) if failed_rho.size: print("\nRequested states excluded because full-AA did not converge") for index, rho in enumerate(failed_rho): iterations = int(state["failed_stage2_iters"][index]) error = float(state["failed_stage2_error"][index]) error_text = f"{error:.3e}" if np.isfinite(error) else "not recorded" print( f" rho={rho:g} g/cc: stage2_iters={iterations}, " f"last_error={error_text}; {state['failed_message'][index]}" ) def _plot(state: dict[str, np.ndarray]) -> None: rho = np.asarray(state["rho_g_cc"], dtype=float) zbar = np.asarray(state["zbar"], dtype=float) zstar = np.asarray(state["zstar"], dtype=float) mu = np.asarray(state["mu_ha"], dtype=float) all_labels = tuple(np.asarray(state["level_labels"], dtype=str)) displayed_indices = [ index for index, label in enumerate(all_labels) if label in DISPLAYED_SHELLS ] labels = tuple(all_labels[index] for index in displayed_indices) energies = np.asarray( state["level_energy_relative_edge_ev"][:, displayed_indices], dtype=float, ) is_bound = np.asarray( state["level_is_bound"][:, displayed_indices], dtype=bool, ) level_q_ion_ws = np.asarray( state["level_q_ion_ws"][:, displayed_indices], dtype=float, ) q_ion_is_available = np.asarray( state["level_q_ion_ws_is_available"][:, displayed_indices], dtype=bool, ) has_shell_q_ion = bool(np.asarray(state["level_q_ion_ws_available"]).item()) failed_rho = np.asarray(state["failed_rho_g_cc"], dtype=float) colors = PALETTES["bing"] reference_curves = _load_bethkenhagen_reference() with style_context("thesis", palette="bing"): fig_ionization, (ax_z, ax_mu) = plt.subplots( 1, 2, figsize=grid_figsize(1, 2), ) published_styles = { "DFT-MD": (colors[2], ":", "o"), "Purgatorio": (colors[3], "--", None), "OPAL": (colors[4], "--", None), "ATOMIC": (colors[5], ":", None), "BU-EK": (colors[6], ":", None), "BU-SP": (colors[7], "-.", None), "BU-SP + Pauli blocking": ("0.48", "--", None), } for label, (rho_reference, z_reference) in reference_curves.items(): color, line_style, marker = published_styles[label] ax_z.plot( rho_reference, z_reference, color=color, ls=line_style, lw=2.3, marker=marker, ms=3.8 if marker else None, markerfacecolor="white" if marker else None, markeredgewidth=0.9 if marker else None, alpha=0.95, label=rf"{label}", zorder=2, ) ax_z.plot( rho, zbar, color="0.08", lw=2.3, ls="-", marker="o", ms=3.8, label=r"Otter $\bar Z=Z-Q_{\rm ion}(R_{\rm WS})$", zorder=6, alpha=0.4, ) ax_z.plot( rho, zstar, color=colors[1], lw=2.3, ls="-", marker="o", ms=3.8, label=r"Otter $Z^*=n_e^0/n_i$", zorder=5, alpha=0.4, ) ax_z.set( xscale="log", xlabel=r"$\rho$ [g cm$^{-3}$]", ylabel="mean ionization", title="Mean ionization", xlim=(0.09, 500.0), ) ax_z.legend(ncol=2, fontsize=7.0, loc="best") ax_mu.plot( rho, mu, color=colors[2], lw=2.3, ls="-", marker="o", ms=3.8, alpha=0.45, ) ax_mu.axhline(0.0, color="0.35", ls=":", lw=0.9) ax_mu.set( xscale="log", xlabel=r"$\rho$ [g cm$^{-3}$]", ylabel=r"$\mu$ [Ha]", title="Electron chemical potential", xlim=(0.09, 500.0), ) fig_ionization.suptitle( r"Carbon ionization at $T_e=100$ eV", y=0.985, ) fig_ionization.tight_layout(rect=(0.0, 0.0, 1.0, 0.955)) save_figure( fig_ionization, FIGURE_DIR / "carbon_ionization_100ev", ) fig_levels, (ax_core, ax_outer) = plt.subplots( 1, 2, figsize=grid_figsize(1, 2), ) core_indices = [ index for index, label in enumerate(labels) if label == "1s" ] outer_indices = [ index for index, label in enumerate(labels) if label != "1s" ] ax_core_qion = ax_core.twinx() ax_outer_qion = ax_outer.twinx() ax_core_qion.patch.set_visible(False) ax_outer_qion.patch.set_visible(False) def draw_levels( axis: Any, ion_charge_axis: Any, indices: list[int], ) -> None: for index in indices: energy_mask = is_bound[:, index] axis.plot( rho[energy_mask], energies[energy_mask, index], color=colors[index % len(colors)], lw=2.3, ls="-", marker="o", ms=3.8, alpha=0.65, label=labels[index], ) charge_mask = q_ion_is_available[:, index] ion_charge_axis.plot( rho[charge_mask], level_q_ion_ws[charge_mask, index], color=colors[index % len(colors)], lw=1.5, ls=(0, (4, 2)), marker="x", ms=4.2, alpha=0.80, label=rf"$Q^{{\rm ion}}_{{{labels[index]}}}$", ) axis.axhline( 0.0, color="0.25", ls=":", lw=1.0, label="continuum edge" if not indices else None, ) axis.set_xscale("log") axis.set_xlim(0.09, 500.0) axis.set_xlabel(r"$\rho$ [g cm$^{-3}$]") if axis is ax_outer: ion_charge_axis.set_ylim(0.0, 0.55) else: ion_charge_axis.set_ylim(bottom=0.0) axis.set_ylabel("energy [eV]", labelpad=2) ion_charge_axis.set_ylabel( r"$Q^{\rm ion}_{nl}$ [e]", color="0.25", labelpad=2, ) ion_charge_axis.tick_params(axis="y", colors="0.25") if indices: energy_handles, energy_labels = axis.get_legend_handles_labels() charge_handles, charge_labels = ( ion_charge_axis.get_legend_handles_labels() ) axis.legend( energy_handles + charge_handles, energy_labels + charge_labels, ncol=2, fontsize=8.0, loc="best", ) draw_levels(ax_core, ax_core_qion, core_indices) draw_levels(ax_outer, ax_outer_qion, outer_indices) ax_core.set(title="Core level") ax_outer.set(title="Outer levels") if not outer_indices: ax_outer.text( 0.5, 0.5, "No outer bound level in this scan", transform=ax_outer.transAxes, ha="center", va="center", ) if not has_shell_q_ion: ax_core_qion.set_visible(False) ax_outer_qion.set_visible(False) ax_outer.text( 0.98, 0.04, "shell $Q^{\\rm ion}_{nl}$ requires\n4096-point regeneration", transform=ax_outer.transAxes, ha="right", va="bottom", fontsize=8.0, color="0.35", ) fig_levels.suptitle( r"Carbon orbital levels at $T_e=100$ eV", y=0.985, ) fig_levels.tight_layout(rect=(0.0, 0.0, 1.0, 0.925)) save_figure( fig_levels, FIGURE_DIR / "carbon_bound_levels_100ev", ) if "agg" not in plt.get_backend().lower(): plt.show() def _compute_and_stage() -> dict[str, np.ndarray]: """Assemble the requested grid and stage it without changing baseline.""" state = _compute_scan() path = _save_candidate(state) point_source = np.asarray(state["point_source"], dtype=str) accepted_count = int( np.count_nonzero(point_source == "accepted_baseline_seed") ) added_count = int(point_source.size - accepted_count) failed_count = int(np.asarray(state["failed_rho_g_cc"]).size) print( f"Using incrementally assembled Otter data staged at {path}: " f"reused {accepted_count} accepted states; added {added_count} " f"cached or newly calculated states; retained {failed_count} " "nonconverged audit record(s)." ) return state def main() -> None: if RECOMPUTE_WITH_OTTER: state = _compute_and_stage() else: try: state = _load_precomputed() except DensityGridMismatchError as error: print(f"{error} Extending it incrementally.") state = _compute_and_stage() else: print( "Using checksummed Otter data from " f"{BASELINE_PATH.relative_to(ROOT)}." ) _print_state_table(state) _plot(state) if __name__ == "__main__": main() .. _sphx_glr_download_gen_examples_plot_carbon_ionization_levels.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_carbon_ionization_levels.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_carbon_ionization_levels.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_carbon_ionization_levels.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_