Note
Go to the end to download the full example code.
Carbon: local-field-correction sensitivity
This complete example keeps the carbon KS average atom and screening cloud
fixed, then changes only the static electron local-field correction (LFC) in
the response, effective ion potential, and HNC calculation. The state is
rho=5 g/cc at T=2 and 100 eV.
The pseudoatom/QOZ construction follows Starrett and Saumon [2014]. The
models are RPA, Hubbard [Hubbard, 1958], Utsumi–Ichimaru
[Utsumi and Ichimaru, 1982], Chabrier-1990 [Chabrier, 1990], and the
Gregori interpolation [Geldart and Vosko, 1966, Gregori et al., 2007]. Otter’s
Gregori/Geldart–Vosko implementation is an attributed NumPy port from
JaXRTS [Lütgert et al., 2026]; the exact revision and BSD-3-Clause notice
are recorded in THIRD_PARTY_NOTICES.md.
Further numerical definitions and interpretation are given in the LFC model-study notes.
Reproduction
From the root of the complete Otter checkout, using Poetry, run:
poetry run python docs/examples/plot_carbon_lfc_sensitivity.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
Otter 0.4.1
[run] C | rho=5 g/cm^3 | Te=2 eV | Ti=-- | QM electronic structure
[AA] C (Z=6) | Rws=1.859183 a0 | full + external
[AA/full SCF: stage 1]
iter | d_n | d_v
0 | nan | 4.740e-02
1 | 1.650e-02 | 4.204e-02
2 | 2.183e-02 | 3.734e-02
3 | 2.251e-02 | 3.327e-02
4 | 2.068e-02 | 2.982e-02
[AA/full SCF: stage 2]
iter | d_n | d_v
0 | nan | 2.700e-02
1 | 1.566e-02 | 2.482e-02
2 | 1.344e-02 | 2.336e-02
3 | 1.157e-02 | 2.366e-02
4 | 1.015e-02 | 3.053e-02
5 | 1.104e-02 | 3.678e-02
6 | 3.169e-02 | 2.515e-02
7 | 1.281e-02 | 4.140e-02
8 | 4.821e-02 | 3.687e-02
9 | 2.067e-01 | 6.818e-02
10 | 1.636e-01 | 1.196e-02
11 | 1.397e-02 | 7.174e-03
12 | 1.921e-02 | 2.933e-03
13 | 1.722e-02 | 3.221e-03
14 | 1.225e-02 | 8.055e-04
15 | 3.578e-03 | 7.028e-04
16 | 1.950e-03 | 1.350e-04
17 | 2.227e-04 | 6.900e-04
18 | 2.698e-03 | 5.339e-05
19 | 8.840e-05 | 2.119e-04
20 | 8.192e-04 | 2.672e-05
21 | 1.384e-04 | 1.378e-05
22 | 6.093e-05 | 9.861e-06
23 | 2.251e-05 | 5.431e-06
24 | 2.969e-05 | 5.166e-07
25 | 1.868e-06 | 9.417e-07
26 | 2.524e-06 | 1.553e-06
27 | 9.172e-06 | 1.144e-07
[AA/external SCF]
iter | d_n | d_v
0 | 2.425e-07 | 1.646e+12
1 | 7.958e-02 | 1.079e+00
2 | 7.570e-02 | 7.903e-01
3 | 2.367e-02 | 3.181e-01
4 | 8.589e-03 | 1.335e-01
5 | 3.084e-03 | 2.328e-02
6 | 8.076e-04 | 9.383e-03
7 | 2.277e-04 | 3.242e-03
8 | 1.073e-04 | 7.808e-04
9 | 2.151e-05 | 6.749e-04
10 | 8.312e-06 | 5.615e-05
11 | 1.066e-06 | 1.123e-05
12 | 2.787e-07 | 1.986e-05
13 | 2.596e-07 | 5.371e-06
14 | 1.034e-07 | 2.229e-06
15 | 2.844e-08 | 1.029e-06
16 | 1.315e-08 | 1.441e-07
[AA] converged=True | mu=0.73168831 Ha | Zbar=4.00001133 | time=20.032 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -8.90073 | 1.00000
[done] total time=20.033 s
Otter 0.4.1
[run] C | rho=5 g/cm^3 | Te=100 eV | Ti=-- | QM electronic structure
[AA] C (Z=6) | Rws=1.859183 a0 | full + external
[AA/full SCF: stage 1]
iter | d_n | d_v
[candidate] benchmarks/outputs/carbon_lfc_sensitivity/gallery_recomputed/C_rho5gcc_T2eV_lfc.npz
0 | nan | 5.150e-02
1 | 6.308e-02 | 4.744e-02
2 | 5.809e-02 | 4.369e-02
3 | 5.238e-02 | 4.007e-02
4 | 4.598e-02 | 3.632e-02
[AA/full SCF: stage 2]
iter | d_n | d_v
0 | nan | 3.233e-02
1 | 6.946e-02 | 6.753e-02
2 | 2.443e-02 | 1.741e-02
3 | 1.578e-02 | 2.941e-02
4 | 6.491e-02 | 3.204e-03
5 | 1.689e-03 | 1.163e-02
6 | 3.109e-02 | 1.931e-03
7 | 1.419e-03 | 8.089e-04
8 | 2.007e-03 | 3.437e-04
9 | 6.199e-04 | 3.707e-05
10 | 5.061e-05 | 6.452e-05
11 | 2.832e-04 | 4.750e-06
12 | 1.887e-05 | 3.937e-06
13 | 3.145e-05 | 2.736e-07
14 | 1.235e-06 | 2.175e-07
15 | 1.118e-06 | 1.686e-07
16 | 5.670e-07 | 1.431e-07
[AA] unresolved threshold state; retrying with zero-tail matching through l=0; A3 extends to 15 R_geometry
[AA] C (Z=6) | Rws=1.859183 a0 | full + external
[AA/full SCF: stage 1]
iter | d_n | d_v
[AA/full SCF: stage 2]
iter | d_n | d_v
0 | nan | 5.165e-02
1 | 6.311e-02 | 4.760e-02
2 | 5.812e-02 | 4.387e-02
3 | 5.242e-02 | 4.026e-02
4 | 7.941e-02 | 9.767e-02
5 | 3.342e-02 | 1.686e-02
6 | 3.301e-02 | 4.424e-02
7 | 8.548e-02 | 2.012e-02
8 | 2.524e-02 | 1.589e-02
9 | 2.234e-02 | 6.443e-03
10 | 8.170e-03 | 3.416e-03
11 | 6.113e-03 | 3.004e-03
12 | 3.570e-03 | 1.824e-04
13 | 5.979e-04 | 4.761e-05
14 | 2.111e-04 | 4.630e-05
15 | 1.724e-04 | 5.439e-06
16 | 1.724e-05 | 1.661e-06
17 | 9.257e-06 | 1.403e-06
18 | 8.879e-06 | 2.121e-06
19 | 1.524e-05 | 1.604e-07
20 | 8.815e-07 | 1.185e-07
21 | 2.349e-07 | 4.542e-08
22 | 1.868e-07 | 3.671e-08
[AA/external SCF]
iter | d_n | d_v
0 | 2.913e-06 | 3.882e+12
1 | 2.613e-02 | 8.870e-01
2 | 1.889e-02 | 1.093e+00
3 | 1.523e-02 | 1.898e-01
4 | 3.328e-03 | 2.311e-02
5 | 5.278e-04 | 1.990e-02
6 | 2.881e-04 | 2.572e-03
7 | 3.957e-05 | 6.897e-04
8 | 1.543e-05 | 5.987e-05
9 | 1.186e-06 | 1.900e-05
10 | 2.110e-07 | 4.044e-05
11 | 4.396e-07 | 7.532e-06
12 | 8.582e-08 | 7.355e-06
13 | 7.781e-08 | 2.925e-07
14 | 4.791e-09 | 1.375e-06
[AA] converged=True | mu=-6.53601783 Ha | Zbar=4.54732225 | time=168.976 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -10.11695 | 0.72600
2s | -0.00213 | 0.14455
[done] total time=183.570 s
[candidate] benchmarks/outputs/carbon_lfc_sensitivity/gallery_recomputed/C_rho5gcc_T100eV_lfc.npz
T[eV] model max|dg| max|dS| max|dV|[Ha Bohr^3]
2.0 none 0.40140 0.15411 19.55532
2.0 hubbard 0.06154 0.06145 17.72958
2.0 utsumiichimaru 0.00118 0.00054 0.08059
2.0 chabrier1990 0.00000 0.00000 0.00000
2.0 gregori2007 0.00197 0.00188 0.36033
100.0 none 0.00430 0.00537 24.62183
100.0 hubbard 0.00146 0.00400 19.14946
100.0 utsumiichimaru 0.00025 0.00033 1.56503
100.0 chabrier1990 0.00000 0.00000 0.00000
100.0 gregori2007 0.00044 0.00031 1.34053
Low-k audit at T=2 eV, k=0.05633257 Bohr^-1
model chi_ee dchi/chi_ref[%] G/k^2[Bohr^2] V_LFC V_ii
none -2.5214329e-04 -3.0909701e-02 0.0000000e+00 0.00000 93.71404
hubbard -2.5229197e-04 2.8040400e-02 1.8599682e-01 37.28490 130.99895
utsumiichimaru -2.5222157e-04 1.2742161e-04 9.7954327e-02 19.63591 113.34996
chabrier1990 -2.5222125e-04 0.0000000e+00 9.7552303e-02 19.55532 113.26937
gregori2007 -2.5222269e-04 5.6972064e-04 9.9349800e-02 19.91565 113.62969
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,
continue_plasma_workflow_from_electronic_result,
solve_plasma_workflow,
)
from otter.ionic import decompose_effective_pair_potential_k
from otter.plotting import grid_figsize, save_figure, style_context
# =============================================================================
# User input
# =============================================================================
RECOMPUTE_WITH_OTTER = True
if os.environ.get("OTTER_RECOMPUTE_CARBON_LFC", "0") == "1":
RECOMPUTE_WITH_OTTER = True
# Use the ionic comparison as the gallery-card preview.
sphinx_gallery_thumbnail_number = 3
RHO_G_CC = 5.0
TEMPERATURES_EV = (2.0, 100.0)
LFC_MODELS = (
"none",
"hubbard",
"utsumiichimaru",
"chabrier1990",
"gregori2007",
)
REFERENCE_LFC = "chabrier1990"
MAX_STATE_WORKERS = 1
HNC_TOL = 1.0e-5
# Audit limit from the default max(10*hnc_tol, 1e-6) transform check.
HNC_CLOSURE_TOL = 1.0e-4
R_RETAIN_MAX_BOHR = 20.0
K_RETAIN_MAX_BOHR_INV = 20.0
LOW_K_PLOT_MAX_BOHR_INV = 1.2
# =============================================================================
MODEL_LABELS = {
"none": "RPA",
"hubbard": "Hubbard",
"utsumiichimaru": "Utsumi–Ichimaru",
"chabrier1990": "Chabrier-1990",
"gregori2007": "Gregori interpolation",
}
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" / "carbon_lfc_sensitivity"
OUTPUT_DIR = ROOT / "benchmarks" / "outputs" / "carbon_lfc_sensitivity"
CANDIDATE_DIR = OUTPUT_DIR / "gallery_recomputed"
FIGURE_DIR = OUTPUT_DIR / "figures"
def sha256_file(path: Path) -> str:
"""Return a streaming SHA-256 digest."""
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 enabling object deserialization."""
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_reviewed_states() -> list[dict[str, np.ndarray]]:
"""Verify and load the accepted current-Otter calculations."""
manifest_path = BASELINE_DIR / "manifest.json"
manifest = json.loads(manifest_path.read_text(encoding="utf-8"))
if manifest.get("benchmark_id") != "carbon_lfc_sensitivity":
raise ValueError("Unexpected carbon LFC manifest.")
states: list[dict[str, np.ndarray]] = []
for record in manifest["states"]:
path = BASELINE_DIR / str(record["data_file"])
if sha256_file(path) != str(record["data_sha256"]):
raise RuntimeError(f"Checksum mismatch for {path}.")
state = load_npz(path)
if tuple(str(value) for value in state["model_labels"]) != LFC_MODELS:
raise ValueError(f"Unexpected LFC model order in {path}.")
states.append(state)
return states
def workflow_config(
temperature_ev: float,
*,
ion_temperature_ev: float | None,
lfc_model: str,
) -> PlasmaWorkflowConfig:
"""Construct one public Otter workflow configuration."""
return PlasmaWorkflowConfig(
elements=["C"],
temperature_ev=float(temperature_ev),
ion_temperature_ev=(
None if ion_temperature_ev is None else float(ion_temperature_ev)
),
rho_g_cc=RHO_G_CC,
qoz_response_lfc_model=str(lfc_model),
hnc_tol=HNC_TOL,
)
def _validate_electronic(electronic: dict[str, Any]) -> None:
"""Refuse to propagate an unreliable average atom into QOZ/HNC."""
if electronic.get("stage2_converged") is not True:
raise RuntimeError("The full carbon average atom did not converge.")
if dict(electronic.get("ext_status", {})).get("converged") is not True:
raise RuntimeError("The external carbon average atom did not converge.")
status = str(electronic.get("threshold_state_status", "")).lower()
if status == "unresolved":
raise RuntimeError("The carbon threshold-state representation is unresolved.")
def solve_temperature(temperature_ev: float) -> dict[str, np.ndarray]:
"""Solve one KS state, then reuse it for every LFC/QOZ/HNC branch."""
electronic_started = time.perf_counter()
shared = solve_plasma_workflow(
workflow_config(
temperature_ev,
ion_temperature_ev=None,
lfc_model=REFERENCE_LFC,
)
)
electronic_elapsed_s = time.perf_counter() - electronic_started
electronic_kind = str(shared["electronic"]["kind"])
electronic = dict(shared["electronic"]["result"])
_validate_electronic(electronic)
ion_results: list[dict[str, Any]] = []
ion_elapsed_s: list[float] = []
for model in LFC_MODELS:
ion_started = time.perf_counter()
workflow = continue_plasma_workflow_from_electronic_result(
workflow_config(
temperature_ev,
ion_temperature_ev=temperature_ev,
lfc_model=model,
),
electronic_kind=electronic_kind,
electronic_result=electronic,
)
ion_elapsed_s.append(time.perf_counter() - ion_started)
ion = dict(workflow["ion"])
if float(ion["hnc_output_residual"]) > HNC_TOL:
raise RuntimeError(f"{model}: HNC residual exceeds tolerance.")
if float(ion["closure_transform_max_abs"]) > HNC_CLOSURE_TOL:
raise RuntimeError(f"{model}: g/S closure exceeds tolerance.")
ion_results.append(ion)
electronic_r = np.asarray(electronic["r"], dtype=float)
electronic_mask = electronic_r <= R_RETAIN_MAX_BOHR
source_r = np.asarray(ion_results[0]["r"], dtype=float)
source_k = np.asarray(ion_results[0]["k"], dtype=float)
r = source_r[source_r <= R_RETAIN_MAX_BOHR]
k = source_k[source_k <= K_RETAIN_MAX_BOHR_INV]
def on_r(ion: dict[str, Any], key: str) -> np.ndarray:
return np.interp(
r,
np.asarray(ion["r"], dtype=float),
np.asarray(ion[key], dtype=float).reshape(-1),
)
def on_k(ion: dict[str, Any], key: str) -> np.ndarray:
return np.interp(
k,
np.asarray(ion["k"], dtype=float),
np.asarray(ion[key], dtype=float).reshape(-1),
)
return {
"schema_version": np.asarray("otter_gallery_carbon_lfc_v1"),
"rho_g_cc": np.asarray(RHO_G_CC),
"temperature_ev": np.asarray(float(temperature_ev)),
"model_labels": np.asarray(LFC_MODELS),
"reference_model": np.asarray(REFERENCE_LFC),
"electronic_elapsed_s": np.asarray(electronic_elapsed_s),
"ion_elapsed_s": np.asarray(ion_elapsed_s),
"electronic_r_bohr": electronic_r[electronic_mask],
"n_bound_bohr3": np.asarray(electronic["n_bound"])[electronic_mask],
"n_ion_bohr3": np.asarray(electronic["n_ion"])[electronic_mask],
"n_scr_bohr3": np.asarray(electronic["n_scr"])[electronic_mask],
"r_bohr": r,
"k_bohr_inv": k,
"gii_r": np.asarray([on_r(ion, "gii_r") for ion in ion_results]),
"sii_k": np.asarray([on_k(ion, "sii_k") for ion in ion_results]),
"vii_k_ha_bohr3": np.asarray(
[on_k(ion, "vii_k") for ion in ion_results]
),
"chi_ee_k": np.asarray(
[on_k(ion, "chi_ee_k") for ion in ion_results]
),
"gee_k": np.asarray([on_k(ion, "gee_k") for ion in ion_results]),
"n_scr_k": on_k(ion_results[0], "n_scr_k"),
"chi0_k": on_k(ion_results[0], "chi0_k"),
"zbar_qoz": np.asarray(
[np.asarray(ion["zbar_qoz"]).item() for ion in ion_results]
),
"hnc_output_residual": np.asarray(
[float(ion["hnc_output_residual"]) for ion in ion_results]
),
}
def low_k_decomposition(
state: dict[str, np.ndarray],
) -> dict[str, np.ndarray]:
"""Decompose the finite low-k potential for every LFC branch."""
k = np.asarray(state["k_bohr_inv"], dtype=float)
q = np.asarray(state["n_scr_k"], dtype=float)
chi0 = np.asarray(state["chi0_k"], dtype=float)
gee = np.asarray(state["gee_k"], dtype=float)
zbar = np.asarray(state["zbar_qoz"], dtype=float)
state_models = tuple(str(value) for value in state["model_labels"])
if zbar.ndim == 0:
zbar = np.full(len(state_models), float(zbar))
if zbar.shape != (len(state_models),):
raise ValueError("zbar_qoz must contain one value per LFC branch.")
charge_terms: list[np.ndarray] = []
lfc_terms: list[np.ndarray] = []
chi0_terms: list[np.ndarray] = []
for model_index in range(len(state_models)):
charge, lfc, chi0_term = decompose_effective_pair_potential_k(
n_scr_k=q[np.newaxis, :],
zbar=np.asarray([zbar[model_index]], dtype=float),
k=k,
chi0_k=chi0,
gee_k=gee[model_index],
)
charge_terms.append(charge[0, 0])
lfc_terms.append(lfc[0, 0])
chi0_terms.append(chi0_term[0, 0])
charge_array = np.asarray(charge_terms)
lfc_array = np.asarray(lfc_terms)
chi0_array = np.asarray(chi0_terms)
reconstructed = charge_array + lfc_array + chi0_array
stored = np.asarray(state["vii_k_ha_bohr3"], dtype=float)
if not np.allclose(reconstructed, stored, rtol=2.0e-11, atol=2.0e-9):
raise RuntimeError("Stable low-k decomposition does not reproduce V_ii(k).")
return {
"g_over_k2": gee / k[np.newaxis, :] ** 2,
"charge": charge_array,
"lfc": lfc_array,
"chi0": chi0_array,
"shared_non_lfc": charge_array + chi0_array,
}
def solve_all_states() -> list[dict[str, np.ndarray]]:
"""Calculate independent temperatures with bounded process parallelism."""
CANDIDATE_DIR.mkdir(parents=True, exist_ok=True)
solved: dict[float, dict[str, np.ndarray]] = {}
with ProcessPoolExecutor(max_workers=MAX_STATE_WORKERS) as pool:
futures = {
pool.submit(solve_temperature, float(temperature)): float(temperature)
for temperature in TEMPERATURES_EV
}
for future in as_completed(futures):
temperature = futures[future]
state = future.result()
solved[temperature] = state
token = f"{temperature:g}".replace(".", "p")
path = CANDIDATE_DIR / f"C_rho5gcc_T{token}eV_lfc.npz"
np.savez_compressed(path, **state)
print(f"[candidate] {path}")
return [solved[float(value)] for value in TEMPERATURES_EV]
def main() -> None:
states = (
solve_all_states()
if RECOMPUTE_WITH_OTTER
else load_reviewed_states()
)
model_labels = LFC_MODELS
reference_index = model_labels.index(REFERENCE_LFC)
print(
"Using "
+ (
"new results calculated directly by this example."
if RECOMPUTE_WITH_OTTER
else "checksummed, reviewed Otter results."
)
)
print(
f"{'T[eV]':>7s} {'model':>18s} {'max|dg|':>11s} "
f"{'max|dS|':>11s} {'max|dV|[Ha Bohr^3]':>20s}"
)
for state in states:
r = np.asarray(state["r_bohr"], dtype=float)
k = np.asarray(state["k_bohr_inv"], dtype=float)
for model_index, model in enumerate(model_labels):
dg = np.asarray(state["gii_r"])[model_index] - np.asarray(
state["gii_r"]
)[reference_index]
ds = np.asarray(state["sii_k"])[model_index] - np.asarray(
state["sii_k"]
)[reference_index]
dv = np.asarray(state["vii_k_ha_bohr3"])[model_index] - np.asarray(
state["vii_k_ha_bohr3"]
)[reference_index]
print(
f"{float(state['temperature_ev']):7.1f} "
f"{model:>18s} "
f"{np.max(np.abs(dg[r <= 20.0])):11.5f} "
f"{np.max(np.abs(ds[k <= 6.0])):11.5f} "
f"{np.max(np.abs(dv[k <= 6.0])):20.5f}"
)
cold_state = min(states, key=lambda item: float(item["temperature_ev"]))
cold_low_k = low_k_decomposition(cold_state)
cold_k = np.asarray(cold_state["k_bohr_inv"], dtype=float)
cold_chi = np.asarray(cold_state["chi_ee_k"], dtype=float)
cold_vii = np.asarray(cold_state["vii_k_ha_bohr3"], dtype=float)
first_mode = 0
print(
"\nLow-k audit at "
f"T={float(cold_state['temperature_ev']):g} eV, "
f"k={cold_k[first_mode]:.8f} Bohr^-1"
)
print(
f"{'model':>18s} {'chi_ee':>14s} {'dchi/chi_ref[%]':>18s} "
f"{'G/k^2[Bohr^2]':>18s} {'V_LFC':>12s} {'V_ii':>12s}"
)
for model_index, model in enumerate(model_labels):
relative_chi = 100.0 * (
cold_chi[model_index, first_mode]
/ cold_chi[reference_index, first_mode]
- 1.0
)
print(
f"{model:>18s} "
f"{cold_chi[model_index, first_mode]:14.7e} "
f"{relative_chi:18.7e} "
f"{cold_low_k['g_over_k2'][model_index, first_mode]:18.7e} "
f"{cold_low_k['lfc'][model_index, first_mode]:12.5f} "
f"{cold_vii[model_index, first_mode]:12.5f}"
)
# %%
# All LFC branches share the same electronic density at each temperature.
# The shell-weighted profiles make the radial electron numbers legible.
with style_context("thesis", palette="bing"):
fig_density, density_axes = plt.subplots(
1,
len(states),
figsize=grid_figsize(1, len(states)),
squeeze=False,
)
for axis, state in zip(density_axes[0], states, strict=True):
electronic_r = np.asarray(state["electronic_r_bohr"], dtype=float)
shell_weight = 4.0 * np.pi * electronic_r**2
mask = electronic_r <= 12.0
for key, label in (
("n_bound_bohr3", r"$n_{\rm bound}$"),
("n_ion_bohr3", r"$n_{\rm ion}$"),
("n_scr_bohr3", r"$n_{\rm scr}$"),
):
axis.plot(
electronic_r[mask],
(shell_weight * np.asarray(state[key], dtype=float))[mask],
label=label,
)
axis.set(
xlim=(-0.5, 12.0),
xlabel=r"$r$ [Bohr]",
ylabel=r"$4\pi r^2n(r)$ [Bohr$^{-1}$]",
title=rf"$T={float(state['temperature_ev']):g}$ eV",
)
density_axes[0, 0].legend(frameon=False)
fig_density.suptitle(
rf"C, $\rho={RHO_G_CC:g}$ g cm$^{{-3}}$: shared KS densities",
y=0.985,
)
fig_density.tight_layout(rect=(0.0, 0.0, 1.0, 0.94))
save_figure(
fig_density,
FIGURE_DIR / "carbon_lfc_electronic_structure",
close=False,
)
# %%
# The response chain exposes the effect of each LFC on
# :math:`G_{ee}`, :math:`\chi_{ee}`, :math:`G_{ee}/k^2`, and
# :math:`V_{ii}`. The finite small-:math:`k` limit of
# :math:`G_{ee}/k^2` is the quantity that enters the LFC contribution to the
# effective potential.
with style_context("thesis", palette="bing"):
colors = dict(
zip(model_labels, plt.rcParams["axes.prop_cycle"].by_key()["color"])
)
fig_response, response_axes = plt.subplots(
len(states),
4,
figsize=grid_figsize(len(states), 4),
squeeze=False,
)
for row, state in enumerate(states):
k = np.asarray(state["k_bohr_inv"], dtype=float)
decomposition = low_k_decomposition(state)
response_axes[row, 1].plot(
k,
np.asarray(state["chi0_k"], dtype=float),
color="0.25",
ls=":",
label=r"shared $\chi^0_{ee}$",
)
for model_index, model in enumerate(model_labels):
label = MODEL_LABELS[model]
color = colors[model]
response_axes[row, 0].plot(
k,
np.asarray(state["gee_k"], dtype=float)[model_index],
color=color,
label=label,
)
response_axes[row, 1].plot(
k,
np.asarray(state["chi_ee_k"], dtype=float)[model_index],
color=color,
label=label,
)
response_axes[row, 2].plot(
k,
decomposition["g_over_k2"][model_index],
color=color,
label=label,
)
response_axes[row, 3].plot(
k,
np.asarray(state["vii_k_ha_bohr3"], dtype=float)[model_index],
color=color,
label=label,
)
temperature = float(state["temperature_ev"])
response_axes[row, 0].set_ylabel(
rf"{temperature:g} eV" + "\n" + r"$G_{ee}(k)$"
)
response_axes[row, 1].set_ylabel(
r"$\chi_{ee}$ [Bohr$^{-3}$ Ha$^{-1}$]"
)
response_axes[row, 2].set_ylabel(
r"$G_{ee}(k)/k^2$ [Bohr$^2$]"
)
response_axes[row, 3].set_ylabel(
r"$V_{ii}$ [Ha Bohr$^3$]"
)
for axis in response_axes[row]:
axis.set_xlim(0.0, 6.0)
for axis, title in zip(
response_axes[0],
(
r"$G_{ee}(k)$",
r"$\chi_{ee}(k)$",
r"$G_{ee}(k)/k^2$",
r"$V_{ii}(k)$",
),
strict=True,
):
axis.set_title(title)
for axis in response_axes[-1]:
axis.set_xlabel(r"$k$ [Bohr$^{-1}$]")
handles, labels = response_axes[0, 0].get_legend_handles_labels()
fig_response.legend(
handles,
labels,
loc="upper center",
bbox_to_anchor=(0.5, 0.945),
ncol=5,
frameon=False,
)
fig_response.suptitle(
rf"C, $\rho={RHO_G_CC:g}$ g cm$^{{-3}}$: LFC response chain",
y=0.985,
)
fig_response.tight_layout(rect=(0.0, 0.0, 1.0, 0.91))
save_figure(
fig_response,
FIGURE_DIR / "carbon_lfc_response",
close=False,
)
# %%
# Why can nearly identical :math:`\chi_{ee}` produce different
# :math:`V_{ii}`?
#
# The Starrett--Saumon reduction contains the *inverse* response,
#
# .. math::
#
# V_{ii}(k)=\frac{4\pi\bar Z^2}{k^2}
# +\frac{q(k)^2}{\chi_{ee}(k)}.
#
# Every Coulomb LFC considered here has
# :math:`G(k)=a k^2+O(k^4)`. Consequently all models share the leading
# :math:`\chi_{ee}(k)=-k^2/(4\pi)+O(k^4)`, so the response curves overlap.
# The LFC-dependent part of the potential nevertheless has the finite limit
#
# .. math::
#
# V_{\rm LFC}(k)=\frac{4\pi}{k^2}G(k)q(k)^2
# \longrightarrow 4\pi a\bar Z^2.
#
# The panels below show the cold-state response, the LFC contribution, and
# the resulting total potential. The remaining terms are identical for all
# branches because :math:`q(k)`, :math:`\bar Z`, and :math:`\chi_0(k)` are
# held fixed.
with style_context("thesis", palette="bing"):
colors = dict(
zip(model_labels, plt.rcParams["axes.prop_cycle"].by_key()["color"])
)
fig_low_k, low_k_axes = plt.subplots(
2,
2,
figsize=grid_figsize(2, 2),
squeeze=False,
)
low_k_mask = cold_k <= LOW_K_PLOT_MAX_BOHR_INV
shared_non_lfc = np.asarray(cold_low_k["shared_non_lfc"], dtype=float)
if not np.allclose(
shared_non_lfc,
shared_non_lfc[reference_index],
rtol=1.0e-12,
atol=1.0e-10,
):
raise RuntimeError("Expected a shared low-k charge/chi0 term.")
for model_index, model in enumerate(model_labels):
label = MODEL_LABELS[model]
color = colors[model]
plot_options = {
"color": color,
"label": label,
"marker": "o",
"markevery": [0],
"ms": 4.0,
}
low_k_axes[0, 0].plot(
cold_k[low_k_mask],
cold_chi[model_index, low_k_mask],
**plot_options,
)
low_k_axes[0, 1].plot(
cold_k[low_k_mask],
100.0
* (
cold_chi[model_index, low_k_mask]
/ cold_chi[reference_index, low_k_mask]
- 1.0
),
**plot_options,
)
low_k_axes[1, 0].plot(
cold_k[low_k_mask],
cold_low_k["lfc"][model_index, low_k_mask],
**plot_options,
)
low_k_axes[1, 1].plot(
cold_k[low_k_mask],
cold_vii[model_index, low_k_mask],
**plot_options,
)
low_k_axes[0, 0].set(
ylabel=r"$\chi_{ee}$ [Bohr$^{-3}$ Ha$^{-1}$]",
title=r"Interacting response $\chi_{ee}(k)$",
)
low_k_axes[0, 1].set(
ylabel=r"relative $\chi_{ee}$ difference [%]",
title="Relative difference from Chabrier-1990",
)
low_k_axes[1, 0].set(
xlabel=r"$k$ [Bohr$^{-1}$]",
ylabel=r"$V_{\rm LFC}$ [Ha Bohr$^3$]",
title=r"LFC term in $V_{ii}(k)$",
)
low_k_axes[1, 1].set(
xlabel=r"$k$ [Bohr$^{-1}$]",
ylabel=r"$V_{ii}$ [Ha Bohr$^3$]",
title=r"Total $V_{ii}(k)$",
)
for axis in low_k_axes.flat:
axis.set_xlim(0.0, LOW_K_PLOT_MAX_BOHR_INV)
low_k_axes[0, 1].axhline(0.0, color="0.55", lw=0.8, ls=":")
handles, labels = low_k_axes[1, 1].get_legend_handles_labels()
fig_low_k.legend(
handles,
labels,
loc="upper center",
bbox_to_anchor=(0.5, 0.94),
ncol=3,
frameon=False,
)
fig_low_k.suptitle(
rf"C, $\rho={RHO_G_CC:g}$ g cm$^{{-3}}$, "
rf"$T={float(cold_state['temperature_ev']):g}$ eV: "
r"small-$k$ LFC decomposition",
y=0.99,
)
fig_low_k.tight_layout(rect=(0.0, 0.0, 1.0, 0.85))
save_figure(
fig_low_k,
FIGURE_DIR / "carbon_lfc_low_k_audit",
close=False,
)
# %%
# The same branches now feed HNC. Comparing 2 and 100 eV displays both a
# strongly LFC-sensitive regime and a regime in which the curves nearly
# coincide, without treating any approximation as exact reference data.
with style_context("thesis", palette="bing"):
colors = dict(
zip(model_labels, plt.rcParams["axes.prop_cycle"].by_key()["color"])
)
fig_ionic, ionic_axes = plt.subplots(
len(states),
2,
figsize=grid_figsize(len(states), 2),
squeeze=False,
)
for row, state in enumerate(states):
r = np.asarray(state["r_bohr"], dtype=float)
k = np.asarray(state["k_bohr_inv"], dtype=float)
for model_index, model in enumerate(model_labels):
label = MODEL_LABELS[model]
color = colors[model]
ionic_axes[row, 0].plot(
r,
np.asarray(state["gii_r"], dtype=float)[model_index],
color=color,
label=label,
)
ionic_axes[row, 1].plot(
k,
np.asarray(state["sii_k"], dtype=float)[model_index],
color=color,
label=label,
)
temperature = float(state["temperature_ev"])
ionic_axes[row, 0].set_ylabel(
rf"{temperature:g} eV" + "\n" + r"$g_{ii}(r)$"
)
ionic_axes[row, 0].set_xlim(-0.5, 12.0)
ionic_axes[row, 1].set_xlim(0.0, 6.0)
ionic_axes[row, 1].axhline(1.0, color="0.55", lw=0.8, ls=":")
ionic_axes[0, 0].set_title(r"$g_{ii}(r)$")
ionic_axes[0, 1].set_title(r"$S_{ii}(k)$")
ionic_axes[-1, 0].set_xlabel(r"$r$ [Bohr]")
ionic_axes[-1, 1].set_xlabel(r"$k$ [Bohr$^{-1}$]")
handles, labels = ionic_axes[0, 0].get_legend_handles_labels()
fig_ionic.legend(
handles,
labels,
loc="upper center",
bbox_to_anchor=(0.5, 0.945),
ncol=5,
frameon=False,
)
fig_ionic.suptitle(
rf"C, $\rho={RHO_G_CC:g}$ g cm$^{{-3}}$: ionic LFC sensitivity",
y=0.985,
)
fig_ionic.tight_layout(rect=(0.0, 0.0, 1.0, 0.91))
save_figure(
fig_ionic,
FIGURE_DIR / "carbon_lfc_ionic_structure",
close=False,
)
plt.show()
if __name__ == "__main__":
main()