Note
Go to the end to download the full example code.
Al: elastic Rayleigh weight versus density
This self-contained example evaluates the static ionic Rayleigh weight
\[W(k)=\left|q(k)+f(k)\right|^2 S_{ii}(k),\]
for aluminium at \(T_e=T_i=10\,\mathrm{eV}\) and three mass densities.
Here q(k) is Otter’s charge-closed screening-cloud transform and
f(k) is the ionic form factor (the transform of \(n_{\rm ion}\),
including the pressure-ionization and radial-cutoff partition, not simply
the full bound-state density). Both are exported by
otter.io.state.build_state_arrays(),
so the expression uses exactly the quantities supplied to the QOZ/HNC
workflow. The convention is the static form-factor factor used in the
ion-structure construction of Starrett and Saumon (2014),
Starrett and Saumon [2014].
Reproduction
From the root of the complete Otter checkout, using Poetry, run:
poetry run python docs/examples/plot_al_rayleigh_weight.py
Downloads are optional: .ipynb launches this repository script; .zip
contains both formats. See Reproducing examples and benchmarks for setup.
The script calculates the states from their input parameters and then plots the results. No bundled Otter NPZ is required. Numerical outputs are written locally; literature reference tables remain inputs to the comparison.
Recorded results
The figures and output below are from the recorded validation run; running the source recalculates them with the installed Otter version.
Recorded terminal output
[compute] Al, rho=2.7 g/cc, Te=Ti=10 eV
Otter 0.4.1
[run] Al | rho=2.7 g/cm^3 | Te=10 eV | Ti=10 eV | QM -> QOZ/HNC
[AA] Al (Z=13) | Rws=2.990107 a0 | full + external
[AA/full SCF: stage 1]
iter | d_n | d_v
0 | nan | 5.450e-02
1 | 2.761e-02 | 4.831e-02
2 | 2.464e-02 | 4.251e-02
3 | 2.174e-02 | 3.713e-02
4 | 1.875e-02 | 3.226e-02
[AA/full SCF: stage 2]
iter | d_n | d_v
0 | nan | 2.796e-02
1 | 1.311e-02 | 2.431e-02
2 | 1.071e-02 | 2.136e-02
3 | 8.686e-03 | 1.911e-02
4 | 7.062e-03 | 1.757e-02
5 | 5.840e-03 | 1.677e-02
6 | 5.109e-03 | 5.097e-02
7 | 5.542e-03 | 1.612e-02
8 | 1.951e-03 | 5.850e-03
9 | 1.327e-03 | 3.203e-02
10 | 1.943e-02 | 6.165e-04
11 | 2.158e-04 | 1.352e-02
12 | 4.772e-03 | 2.834e-03
13 | 8.687e-04 | 1.965e-03
14 | 1.817e-03 | 2.282e-03
15 | 6.269e-04 | 7.042e-03
16 | 2.492e-03 | 1.146e-03
17 | 4.615e-04 | 5.551e-04
18 | 5.509e-04 | 4.591e-04
19 | 3.113e-04 | 1.640e-04
20 | 1.161e-04 | 1.792e-04
21 | 1.257e-04 | 3.195e-05
22 | 1.956e-05 | 2.689e-05
23 | 1.069e-05 | 4.722e-06
24 | 1.901e-06 | 3.522e-06
25 | 5.912e-06 | 1.392e-06
26 | 2.038e-06 | 7.264e-07
[AA/external SCF]
iter | d_n | d_v
0 | 1.972e-07 | 2.337e+12
1 | 7.396e-02 | 1.080e+00
2 | 7.210e-02 | 8.208e-01
3 | 2.566e-02 | 3.140e-01
4 | 7.203e-03 | 1.398e-01
5 | 2.764e-03 | 5.046e-02
6 | 1.199e-03 | 9.587e-03
7 | 2.265e-04 | 7.774e-03
8 | 1.342e-04 | 1.402e-03
9 | 3.052e-05 | 3.275e-04
10 | 5.860e-06 | 6.840e-04
11 | 8.270e-06 | 2.500e-05
12 | 8.372e-07 | 3.316e-05
13 | 6.986e-07 | 2.646e-05
14 | 2.850e-07 | 3.368e-06
15 | 7.867e-08 | 3.482e-06
16 | 5.003e-08 | 7.467e-07
[AA] converged=True | mu=-0.02168056 Ha | Zbar=3.02673844 | time=58.071 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -54.58452 | 1.00000
2s | -3.38869 | 0.99990
2p | -2.01976 | 0.99567
[ion] HNC iterations=10 QOZ=1.931 s HNC=0.030 s
[done] total time=60.049 s
mu=-0.02168056 Ha, Zbar=3.026684, HNC=2.313e-07
[compute] Al, rho=8.1 g/cc, Te=Ti=10 eV
Otter 0.4.1
[run] Al | rho=8.1 g/cm^3 | Te=10 eV | Ti=10 eV | QM -> QOZ/HNC
[AA] Al (Z=13) | Rws=2.073224 a0 | full + external
[AA/full SCF: stage 1]
iter | d_n | d_v
0 | nan | 5.240e-02
1 | 2.756e-02 | 4.637e-02
2 | 2.453e-02 | 4.083e-02
3 | 2.156e-02 | 3.587e-02
4 | 1.852e-02 | 3.162e-02
[AA/full SCF: stage 2]
iter | d_n | d_v
0 | nan | 2.819e-02
1 | 1.292e-02 | 2.572e-02
2 | 1.062e-02 | 2.447e-02
3 | 8.811e-03 | 2.497e-02
4 | 7.948e-03 | 2.865e-02
5 | 8.671e-03 | 3.152e-02
6 | 6.299e-03 | 9.271e-03
7 | 7.210e-03 | 1.988e-02
8 | 3.137e-02 | 2.541e-02
9 | 4.144e-02 | 6.313e-03
10 | 1.230e-02 | 7.173e-03
11 | 3.418e-03 | 9.972e-03
12 | 1.750e-02 | 1.764e-03
13 | 3.516e-03 | 2.529e-03
14 | 3.710e-03 | 4.845e-04
15 | 5.016e-04 | 4.804e-04
16 | 1.066e-03 | 1.295e-04
17 | 2.083e-04 | 7.334e-05
18 | 7.784e-05 | 2.805e-05
19 | 2.125e-05 | 8.400e-06
20 | 1.039e-05 | 2.920e-06
21 | 5.987e-06 | 6.529e-06
22 | 1.777e-05 | 1.271e-06
23 | 2.266e-06 | 1.686e-06
24 | 3.648e-06 | 4.156e-07
25 | 5.985e-07 | 7.073e-08
[AA/external SCF]
iter | d_n | d_v
0 | 1.054e-07 | 3.439e+12
1 | 9.445e-02 | 1.102e+00
2 | 9.526e-02 | 7.481e-01
3 | 2.708e-02 | 2.715e-01
4 | 7.723e-03 | 1.315e-01
5 | 3.286e-03 | 2.658e-02
6 | 7.578e-04 | 6.706e-03
7 | 2.680e-04 | 4.091e-03
8 | 9.485e-05 | 1.184e-03
9 | 3.577e-05 | 9.267e-04
10 | 1.784e-05 | 2.072e-04
11 | 5.290e-06 | 1.298e-04
12 | 1.759e-06 | 3.905e-05
13 | 6.600e-07 | 1.653e-05
14 | 4.628e-07 | 4.255e-06
15 | 7.187e-08 | 8.561e-07
16 | 8.710e-09 | 7.546e-07
[AA] converged=True | mu=0.83162061 Ha | Zbar=3.03997079 | time=41.404 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -53.89385 | 1.00000
2s | -2.75781 | 0.99994
2p | -1.38861 | 0.99763
[ion] HNC iterations=11 QOZ=1.894 s HNC=0.030 s
[done] total time=43.339 s
mu=0.83162061 Ha, Zbar=3.036860, HNC=4.861e-07
[compute] Al, rho=15 g/cc, Te=Ti=10 eV
Otter 0.4.1
[run] Al | rho=15 g/cm^3 | Te=10 eV | Ti=10 eV | QM -> QOZ/HNC
[AA] Al (Z=13) | Rws=1.688279 a0 | full + external
[AA/full SCF: stage 1]
iter | d_n | d_v
0 | nan | 5.130e-02
1 | 2.733e-02 | 4.558e-02
2 | 2.425e-02 | 4.052e-02
3 | 2.130e-02 | 3.624e-02
4 | 1.833e-02 | 3.289e-02
[AA/full SCF: stage 2]
iter | d_n | d_v
0 | nan | 3.078e-02
1 | 1.298e-02 | 3.059e-02
2 | 1.102e-02 | 3.429e-02
3 | 1.209e-02 | 4.811e-02
4 | 1.826e-02 | 2.676e-02
5 | 2.445e-02 | 1.413e-02
6 | 5.230e-03 | 1.354e-02
7 | 2.792e-02 | 5.447e-03
8 | 5.693e-03 | 1.124e-02
9 | 3.316e-02 | 1.020e-03
10 | 2.648e-03 | 1.654e-03
11 | 1.446e-03 | 7.034e-04
12 | 5.545e-04 | 3.153e-04
13 | 8.143e-04 | 2.769e-04
14 | 9.875e-04 | 1.797e-04
15 | 4.388e-04 | 8.007e-06
16 | 1.487e-05 | 1.928e-05
17 | 3.530e-05 | 7.812e-06
18 | 2.699e-05 | 8.464e-07
19 | 2.792e-06 | 5.686e-06
20 | 1.376e-05 | 9.684e-07
21 | 2.672e-06 | 2.363e-07
22 | 5.614e-07 | 1.573e-07
23 | 1.972e-07 | 9.358e-08
[AA/external SCF]
iter | d_n | d_v
0 | 1.327e-07 | 4.370e+12
1 | 9.147e-02 | 1.093e+00
2 | 9.087e-02 | 7.858e-01
3 | 2.869e-02 | 3.023e-01
4 | 8.162e-03 | 1.051e-01
5 | 3.183e-03 | 2.187e-02
6 | 7.680e-04 | 1.411e-02
7 | 3.014e-04 | 8.387e-03
8 | 2.104e-04 | 3.195e-03
9 | 7.794e-05 | 1.063e-03
10 | 1.840e-05 | 1.362e-04
11 | 3.782e-06 | 4.244e-05
12 | 8.600e-07 | 9.623e-05
13 | 1.803e-06 | 7.063e-06
14 | 2.224e-07 | 1.899e-06
15 | 4.053e-08 | 9.100e-07
16 | 2.580e-08 | 1.091e-06
[AA] converged=True | mu=1.64464132 Ha | Zbar=3.16017905 | time=29.905 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -53.25834 | 1.00000
2s | -2.22313 | 0.99997
2p | -0.85585 | 0.99889
[ion] HNC iterations=11 QOZ=1.853 s HNC=0.032 s
[done] total time=31.800 s
mu=1.64464132 Ha, Zbar=3.143753, HNC=6.955e-07
Saved newly calculated states: benchmarks/outputs/al_rayleigh_weight_10ev/recomputed/Al_Te10eV_rayleigh.npz
rho=2.7 g/cc: mu=-0.02168056 Ha, max W=7.982504e+01
rho=8.1 g/cc: mu=0.83162061 Ha, max W=6.476243e+01
rho=15 g/cc: mu=1.64464132 Ha, max W=5.562608e+01
from __future__ import annotations
import hashlib
import json
import os
from pathlib import Path
import time
from typing import Any
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import numpy as np
from otter import PlasmaWorkflowConfig, solve_plasma_workflow
from otter.io.state import StateExportOptions, build_state_arrays
from otter.plotting import grid_figsize, save_figure, style_context
# =============================================================================
# User input
# =============================================================================
USE_PRECOMPUTED_DATA = False
RECOMPUTE_WITH_OTTER = True
if os.environ.get("OTTER_RECOMPUTE_AL_RAYLEIGH", "0") == "1":
USE_PRECOMPUTED_DATA = False
RECOMPUTE_WITH_OTTER = True
if USE_PRECOMPUTED_DATA == RECOMPUTE_WITH_OTTER:
raise ValueError("Select exactly one data mode: precomputed or recompute.")
ELEMENT = "Al"
RHO_VALUES_G_CC = (2.7, 8.1, 15.0)
TE_EV = 10.0
TI_EV = 10.0
HNC_TOL = 1.0e-6
HNC_CLOSURE_TOL = 1.0e-3
# =============================================================================
SCHEMA = "otter_al_rayleigh_weight_10ev_v2"
def repository_root() -> Path:
"""Locate source and reference inputs, independently of numerical outputs."""
candidates = [Path.cwd().resolve(), *Path.cwd().resolve().parents]
source_file = globals().get("__file__")
if source_file is not None:
candidates.extend(Path(source_file).resolve().parents)
for candidate in candidates:
if (candidate / "pyproject.toml").is_file() and (candidate / "src/otter").is_dir():
return candidate
raise FileNotFoundError("Run from an Otter source checkout with its dependencies installed.")
ROOT = repository_root()
BASELINE_DIR = ROOT / "benchmarks" / "baselines" / "al_rayleigh_weight_10ev"
OUTPUT_DIR = ROOT / "benchmarks" / "outputs" / "al_rayleigh_weight_10ev"
FIGURE_DIR = OUTPUT_DIR / "figures"
RECOMPUTED_PATH = OUTPUT_DIR / "recomputed" / "Al_Te10eV_rayleigh.npz"
def sha256_file(path: Path) -> str:
digest = hashlib.sha256()
with path.open("rb") as handle:
for block in iter(lambda: handle.read(1024 * 1024), b""):
digest.update(block)
return digest.hexdigest()
def workflow_config(rho_g_cc: float) -> PlasmaWorkflowConfig:
"""Build the public Otter QM -> QOZ/HNC configuration for one density."""
return PlasmaWorkflowConfig(
elements=[ELEMENT],
temperature_ev=TE_EV,
ion_temperature_ev=TI_EV,
rho_g_cc=float(rho_g_cc),
hnc_tol=HNC_TOL,
hnc_closure_transform_tol=HNC_CLOSURE_TOL,
)
def rayleigh_weight(state: dict[str, np.ndarray]) -> np.ndarray:
"""Reconstruct the exactly derived elastic weight from q, f, and Sii."""
weight = np.abs(state["q_k"] + state["f_k"]) ** 2 * state["sii_k"]
if not np.all(np.isfinite(weight)) or np.any(weight < -1.0e-10):
raise RuntimeError("Rayleigh weight is not finite and non-negative.")
return np.asarray(weight, dtype=float)
def pack_workflow(
workflow: dict[str, Any], rho_g_cc: float, elapsed_s: float
) -> dict[str, np.ndarray]:
"""Extract q, f, S and W from one completed Otter workflow."""
exported = build_state_arrays(
workflow,
options=StateExportOptions(
profile="ion_structure",
),
)
k = np.asarray(exported["k_bohr_inv"], dtype=float)
q = np.asarray(exported["q_k"], dtype=float)[0]
f = np.asarray(exported["f_k"], dtype=float)[0]
sii = np.asarray(exported["sij_k"], dtype=float)[0, 0]
if not (k.shape == q.shape == f.shape == sii.shape):
raise ValueError("q, f, Sii and k do not share a reciprocal grid.")
electronic = dict(workflow["electronic"]["result"])
ion = dict(workflow["ion"])
return {
"schema_version": np.asarray(SCHEMA),
"storage_profile": np.asarray("gallery_analysis"),
"element": np.asarray(ELEMENT),
"rho_g_cc": np.asarray(float(rho_g_cc)),
"te_ev": np.asarray(TE_EV),
"ti_ev": np.asarray(TI_EV),
"k_bohr_inv": k,
"q_k": q,
"f_k": f,
"sii_k": sii,
"mu_ha": np.asarray(float(electronic["mu"])),
"zbar": np.asarray(np.asarray(ion["zbar_qoz"]).item()),
"elapsed_s": np.asarray(float(elapsed_s)),
}
def calculate_states() -> list[dict[str, np.ndarray]]:
"""Run the three independent Otter calculations and write a portable NPZ."""
states: list[dict[str, np.ndarray]] = []
started = time.perf_counter()
for rho in RHO_VALUES_G_CC:
print(f"[compute] Al, rho={rho:g} g/cc, Te=Ti={TE_EV:g} eV")
state_started = time.perf_counter()
workflow = solve_plasma_workflow(workflow_config(rho))
state = pack_workflow(workflow, rho, time.perf_counter() - state_started)
states.append(state)
print(
f" mu={float(state['mu_ha']):.8f} Ha, "
f"Zbar={float(state['zbar']):.6f}, "
f"HNC={float(workflow['ion']['hnc_best_residual']):.3e}"
)
payload: dict[str, np.ndarray] = {
"schema_version": np.asarray(SCHEMA),
"storage_profile": np.asarray("gallery_analysis"),
"state_count": np.asarray(len(states)),
"rho_values_g_cc": np.asarray(RHO_VALUES_G_CC),
"producer_elapsed_s": np.asarray(time.perf_counter() - started),
}
for index, state in enumerate(states):
for key, value in state.items():
if key == "schema_version":
continue
payload[f"state_{index}_{key}"] = np.asarray(value)
RECOMPUTED_PATH.parent.mkdir(parents=True, exist_ok=True)
np.savez_compressed(RECOMPUTED_PATH, **payload)
print(f"Saved newly calculated states: {RECOMPUTED_PATH}")
return states
def load_states() -> list[dict[str, np.ndarray]]:
"""Verify the checksummed archive and unpack its three states."""
manifest = json.loads((BASELINE_DIR / "manifest.json").read_text(encoding="utf-8"))
record = manifest["state"]
archive_path = BASELINE_DIR / str(record["data_file"])
if sha256_file(archive_path) != str(record["data_sha256"]):
raise RuntimeError(f"Checksum mismatch for {archive_path}.")
with np.load(archive_path, allow_pickle=False) as archive:
if str(archive["schema_version"].item()) != SCHEMA:
raise ValueError("Unsupported Al Rayleigh-weight archive schema.")
count = int(archive["state_count"])
states = []
for index in range(count):
prefix = f"state_{index}_"
state = {
key[len(prefix) :]: np.asarray(archive[key])
for key in archive.files
if key.startswith(prefix)
}
states.append(state)
return states
def main() -> None:
states = calculate_states() if RECOMPUTE_WITH_OTTER else load_states()
print(
"Using "
+ (
"newly calculated Otter states."
if RECOMPUTE_WITH_OTTER
else "the checksummed current-Otter states."
)
)
for state in states:
print(
f"rho={float(state['rho_g_cc']):g} g/cc: "
f"mu={float(state['mu_ha']):.8f} Ha, "
f"max W={float(np.max(rayleigh_weight(state))):.6e}"
)
with style_context("thesis", palette="bing"):
fig, axes = plt.subplots(
2,
2,
figsize=grid_figsize(2, 2),
squeeze=False,
)
colours = ("#2E5EAA", "#E76F51", "#23BB62")
for state, colour in zip(states, colours, strict=True):
rho = float(state["rho_g_cc"])
label = rf"$\rho={rho:g}\ \mathrm{{g\,cm^{{-3}}}}$"
k = np.asarray(state["k_bohr_inv"], dtype=float)
axes[0, 0].plot(k, state["q_k"], color=colour, label=label + r", $q$")
axes[0, 0].plot(k, state["f_k"], color=colour, ls="--", label=label + r", $f$")
axes[0, 1].plot(k, state["sii_k"], color=colour, label=label)
axes[1, 0].plot(k, state["q_k"] + state["f_k"], color=colour, label=label)
axes[1, 1].plot(k, rayleigh_weight(state), color=colour, label=label)
axes[0, 0].set(title=r"Form-factor components", ylabel=r"$q(k),\ f(k)$")
axes[0, 1].set(title=r"Ionic structure", ylabel=r"$S_{ii}(k)$")
axes[1, 0].set(title=r"Total form factor", ylabel=r"$q(k)+f(k)$")
axes[1, 1].set(title=r"Rayleigh weight", ylabel=r"$W(k)$")
for axis in axes.flat:
axis.set_xlabel(r"$k$ [Bohr$^{-1}$]")
axis.set_xlim(-0.05, 12.0)
density_handles = [
Line2D(
[0.0],
[0.0],
color=colour,
label=rf"${float(state['rho_g_cc']):g}\ \mathrm{{g\,cm^{{-3}}}}$",
)
for state, colour in zip(states, colours, strict=True)
]
component_handles = [
Line2D([0.0], [0.0], color="0.25", ls="-", label=r"$q(k)$"),
Line2D([0.0], [0.0], color="0.25", ls="--", label=r"$f(k)$"),
]
fig.legend(
handles=density_handles + component_handles,
loc="upper center",
bbox_to_anchor=(0.5, 0.945),
ncol=5,
frameon=False,
)
fig.suptitle(
r"Al: static Rayleigh weight, $T_e=T_i=10$ eV",
y=0.99,
)
fig.tight_layout(rect=(0.0, 0.0, 1.0, 0.89))
# Keep the figure open so Sphinx-Gallery embeds the same multi-density
# result on the ``Al: elastic Rayleigh weight versus density`` HTML page.
# ``save_figure`` still writes the slide-ready PNG/PDF copies; the open
# figure is only additionally used by the documentation renderer.
save_figure(fig, FIGURE_DIR / "al_rayleigh_weight_10ev", close=False)
if __name__ == "__main__":
main()