Note
Go to the end to download the full example code.
Single-species Starrett–Saumon ion-structure benchmarks
This benchmark compares Otter IS-QOZ/HNC pair-distribution functions with Starrett TF numerical extractions from Starrett and Saumon [2014] (doi: 10.1016/j.hedp.2013.12.001). Published data are shown only as open circles. Otter KS-DFT (QM) calculations are requested for every thermodynamic state; the tungsten panels additionally request Thomas–Fermi (TF), matching the source comparison. Curves are shown only when they pass the electronic, HNC, and transform-closure gates.
Every figure is saved as both a documentation PNG and a vector PDF suitable
for slides. The publication extractions are attributed reference data with
license status NOASSERTION; see the provenance and reuse notice.
Reproduction
From the root of the complete Otter checkout, using Poetry, run:
poetry run python benchmarks/examples/plot_starrett_single_species_2013_2014.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] Fe | rho=22.5 g/cm^3 | Te=10 eV | Ti=10 eV | QM -> QOZ/HNC
[AA] Fe (Z=26) | Rws=1.879546 a0 | full + external
[AA/full SCF: stage 1]
iter | d_n | d_v
0 | nan | 5.546e-02
1 | 2.391e-02 | 4.921e-02
2 | 2.266e-02 | 4.332e-02
3 | 1.938e-02 | 3.782e-02
4 | 1.623e-02 | 3.275e-02
[AA/full SCF: stage 2]
iter | d_n | d_v
0 | nan | 2.814e-02
1 | 1.089e-02 | 2.402e-02
2 | 8.709e-03 | 2.043e-02
3 | 6.794e-03 | 1.759e-02
4 | 5.150e-03 | 1.593e-02
5 | 3.848e-03 | 1.592e-02
6 | 2.940e-03 | 1.837e-02
7 | 2.691e-03 | 2.481e-02
8 | 3.587e-03 | 1.638e-02
9 | 2.268e-02 | 9.056e-03
10 | 2.025e-02 | 1.257e-02
11 | 1.179e-02 | 4.228e-03
12 | 3.160e-03 | 2.885e-03
13 | 4.962e-03 | 1.482e-03
14 | 7.618e-04 | 1.664e-03
15 | 2.243e-03 | 4.738e-04
16 | 4.337e-04 | 8.620e-05
17 | 1.798e-04 | 1.042e-04
18 | 1.405e-04 | 2.096e-05
19 | 3.082e-05 | 6.476e-06
20 | 1.132e-05 | 8.191e-06
21 | 1.355e-05 | 6.216e-06
22 | 9.513e-06 | 3.230e-06
23 | 5.206e-06 | 3.560e-06
24 | 5.137e-06 | 1.554e-06
[AA/external SCF]
iter | d_n | d_v
0 | 1.083e-07 | 4.940e+12
1 | 1.008e-01 | 1.102e+00
2 | 1.034e-01 | 7.780e-01
3 | 3.172e-02 | 3.349e-01
4 | 8.881e-03 | 1.125e-01
5 | 3.420e-03 | 1.226e-02
6 | 4.990e-04 | 8.863e-03
7 | 3.095e-04 | 4.717e-03
8 | 1.968e-04 | 1.115e-03
9 | 2.455e-05 | 1.217e-03
10 | 1.957e-05 | 6.066e-04
11 | 1.014e-05 | 1.014e-04
12 | 3.797e-06 | 5.506e-05
13 | 9.599e-07 | 5.040e-05
14 | 7.885e-07 | 7.325e-06
15 | 2.710e-07 | 5.723e-07
16 | 6.513e-09 | 8.070e-07
17 | 1.362e-08 | 6.221e-08
[AA] converged=True | mu=1.41943393 Ha | Zbar=8.31155805 | time=29.712 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -252.38612 | 1.00000
2s | -27.89842 | 1.00000
2p | -23.89132 | 1.00000
3s | -1.88263 | 0.99987
3p | -0.74051 | 0.99721
[ion] HNC iterations=12 QOZ=0.830 s HNC=0.021 s
[done] total time=30.576 s
Otter 0.4.1
[run] Fe | rho=22.5 g/cm^3 | Te=10 eV | Ti=10 eV | TF -> QOZ/HNC
[accepted] fe_10_qm
[AA/TF] Fe (Z=26) Te=10 eV rho=22.5 g/cc Rws=1.879546 Bohr run_mode=full+ext
[AA/TF full SCF]
iter | d_n | d_v
1 | inf | 6.846e-01
10 | 6.927e-02 | 9.930e-03
20 | 2.170e-04 | 1.628e-04
30 | 1.478e-04 | 9.578e-05
[AA/TF external SCF]
iter | d_n | d_v
1 | inf | 2.266e+01
10 | 1.986e-02 | 5.335e-01
20 | 2.047e-06 | 1.755e-04
[AA/TF] converged=True mu=1.62153017 Ha Zbar=5.85298260 time=0.288 s
[ion] HNC iterations=16 QOZ=0.799 s HNC=0.027 s
[done] total time=1.121 s
Otter 0.4.1
[run] H | rho=80 g/cm^3 | Te=5 eV | Ti=5 eV | QM -> QOZ/HNC
[AA] H (Z=1) | Rws=0.323030 a0 | full + external
[AA/full SCF: stage 1]
iter | d_n | d_v
[accepted] fe_10_tf
0 | nan | 2.708e-01
1 | 7.969e-03 | 2.695e-01
2 | 1.170e-03 | 1.394e-01
3 | 2.619e-02 | 6.567e-01
4 | 5.063e-02 | 1.720e-01
[AA/full SCF: stage 2]
iter | d_n | d_v
0 | nan | 2.229e-01
1 | 2.105e-02 | 2.255e-01
2 | 1.421e-02 | 2.161e-01
3 | 8.264e-03 | 1.892e-01
4 | 3.338e-03 | 1.462e-01
5 | 1.021e-03 | 1.030e-01
6 | 9.324e-04 | 5.329e-02
7 | 1.024e-01 | 1.460e-01
8 | 1.864e-02 | 1.037e-01
9 | 2.029e-02 | 1.963e-02
10 | 3.979e-03 | 3.992e-03
11 | 4.670e-03 | 4.608e-03
12 | 5.912e-03 | 3.108e-03
13 | 4.392e-03 | 2.018e-04
14 | 2.616e-04 | 3.461e-04
15 | 3.761e-04 | 4.549e-05
16 | 3.137e-05 | 2.522e-05
17 | 2.624e-05 | 1.786e-06
18 | 6.472e-07 | 6.254e-07
19 | 3.949e-07 | 3.465e-06
20 | 1.561e-06 | 3.694e-08
[AA/external SCF]
iter | d_n | d_v
0 | 1.990e-07 | 1.131e+12
1 | 7.045e-03 | 1.163e+00
2 | 6.962e-03 | 1.102e+00
3 | 5.327e-03 | 3.894e-01
4 | 1.804e-03 | 7.646e-02
5 | 4.036e-04 | 9.911e-03
6 | 6.364e-05 | 8.716e-03
7 | 3.667e-05 | 8.016e-04
8 | 3.084e-06 | 1.229e-03
9 | 4.654e-06 | 9.186e-05
10 | 2.410e-07 | 5.830e-05
11 | 3.678e-07 | 2.468e-05
12 | 7.102e-08 | 7.889e-06
13 | 3.994e-08 | 5.832e-06
14 | 2.417e-08 | 8.071e-07
15 | 5.338e-09 | 4.063e-07
[AA] converged=True | mu=16.63032936 Ha | Zbar=1.00000000 | time=21.448 s
Bound levels
none
[ion] HNC iterations=18 QOZ=0.805 s HNC=0.031 s
[done] total time=22.294 s
[accepted] h_5_qm
Otter 0.4.1
[run] H | rho=80 g/cm^3 | Te=5 eV | Ti=5 eV | TF -> QOZ/HNC
[AA/TF] H (Z=1) Te=5 eV rho=80 g/cc Rws=0.323030 Bohr run_mode=full+ext
[AA/TF full SCF]
iter | d_n | d_v
1 | inf | 5.841e-01
10 | 1.085e-04 | 2.867e-04
20 | 3.345e-07 | 1.249e-04
[AA/TF external SCF]
iter | d_n | d_v
1 | inf | 9.849e+00
10 | 2.154e-06 | 2.764e-04
[AA/TF] converged=True mu=16.65155489 Ha Zbar=0.91661650 time=0.217 s
[ion] HNC iterations=17 QOZ=0.797 s HNC=0.028 s
[done] total time=1.049 s
[accepted] h_5_tf
Otter 0.4.1
[run] H | rho=80 g/cm^3 | Te=172 eV | Ti=172 eV | QM -> QOZ/HNC
[AA] H (Z=1) | Rws=0.323030 a0 | full + external
[AA/full SCF: stage 1]
iter | d_n | d_v
0 | nan | 1.513e-01
1 | 1.940e-01 | 1.866e-01
2 | 3.873e-02 | 8.251e-02
3 | 3.229e-02 | 5.375e-02
4 | 3.958e-02 | 5.259e-02
[AA/full SCF: stage 2]
iter | d_n | d_v
0 | nan | 1.365e-02
1 | 6.870e-03 | 7.363e-03
2 | 5.252e-04 | 1.349e-02
3 | 1.478e-02 | 8.006e-03
4 | 9.554e-03 | 2.666e-03
5 | 3.153e-03 | 1.385e-03
6 | 1.307e-03 | 4.221e-04
7 | 3.294e-04 | 3.981e-05
8 | 2.452e-05 | 5.355e-05
9 | 2.622e-05 | 1.105e-05
10 | 5.027e-06 | 1.741e-06
11 | 1.180e-06 | 5.568e-06
12 | 3.573e-06 | 1.014e-06
[AA/external SCF]
iter | d_n | d_v
0 | 1.221e-06 | 1.149e+12
1 | 6.022e-03 | 1.149e+00
2 | 5.854e-03 | 1.145e+00
3 | 4.823e-03 | 3.535e-01
4 | 1.510e-03 | 7.172e-02
5 | 3.269e-04 | 8.836e-03
6 | 4.912e-05 | 8.159e-03
7 | 2.997e-05 | 6.331e-04
8 | 2.104e-06 | 1.017e-03
9 | 3.307e-06 | 5.877e-05
10 | 1.856e-07 | 3.873e-05
11 | 1.563e-07 | 2.271e-05
12 | 4.849e-08 | 2.163e-06
13 | 8.436e-09 | 2.483e-06
14 | 1.145e-08 | 1.072e-06
[AA] converged=True | mu=14.48604573 Ha | Zbar=1.00000000 | time=8.004 s
Bound levels
none
[ion] HNC iterations=4 QOZ=0.807 s HNC=0.008 s
[done] total time=8.829 s
[accepted] h_172_qm
Otter 0.4.1
[run] H | rho=80 g/cm^3 | Te=172 eV | Ti=172 eV | TF -> QOZ/HNC
[AA/TF] H (Z=1) Te=172 eV rho=80 g/cc Rws=0.323030 Bohr run_mode=full+ext
[AA/TF full SCF]
iter | d_n | d_v
1 | inf | 5.545e-01
10 | 1.120e-04 | 2.597e-04
[AA/TF external SCF]
iter | d_n | d_v
1 | inf | 8.560e+00
10 | 1.052e-06 | 2.516e-04
[AA/TF] converged=True mu=14.41270308 Ha Zbar=0.92608629 time=0.151 s
[ion] HNC iterations=4 QOZ=0.799 s HNC=0.008 s
[done] total time=0.965 s
Otter 0.4.1
[run] C | rho=12.64 g/cm^3 | Te=64.64 eV | Ti=64.64 eV | QM -> QOZ/HNC
[AA] C (Z=6) | Rws=1.364781 a0 | full + external
[AA/full SCF: stage 1]
iter | d_n | d_v
[accepted] h_172_tf
0 | nan | 6.198e-02
1 | 6.714e-02 | 6.108e-02
2 | 7.372e-02 | 6.044e-02
3 | 7.566e-02 | 5.827e-02
4 | 1.555e-02 | 5.657e-02
[AA/full SCF: stage 2]
iter | d_n | d_v
0 | nan | 2.347e-02
1 | 3.337e-02 | 1.718e-02
2 | 2.253e-02 | 3.271e-02
3 | 8.808e-02 | 1.606e-02
4 | 4.018e-02 | 1.562e-02
5 | 4.424e-02 | 6.854e-04
6 | 1.107e-03 | 1.375e-03
7 | 1.729e-03 | 1.790e-03
8 | 4.475e-03 | 6.282e-05
9 | 2.643e-04 | 2.216e-04
10 | 7.650e-04 | 3.537e-06
11 | 5.564e-06 | 6.180e-06
12 | 2.524e-05 | 2.605e-06
13 | 4.250e-06 | 1.643e-06
14 | 4.449e-06 | 6.220e-07
15 | 2.076e-06 | 1.283e-06
[AA/external SCF]
iter | d_n | d_v
0 | 2.163e-06 | 3.426e+12
1 | 4.041e-02 | 9.888e-01
2 | 3.474e-02 | 1.025e+00
3 | 2.051e-02 | 2.201e-01
4 | 5.150e-03 | 4.601e-02
5 | 1.242e-03 | 3.380e-02
6 | 8.024e-04 | 8.688e-03
7 | 1.477e-04 | 2.961e-03
8 | 6.734e-05 | 9.735e-04
9 | 1.818e-05 | 5.779e-05
10 | 7.633e-07 | 1.509e-04
11 | 1.770e-06 | 1.535e-05
12 | 2.805e-07 | 9.012e-06
13 | 2.227e-07 | 4.094e-06
14 | 6.748e-08 | 1.798e-06
15 | 3.299e-08 | 5.095e-07
[AA] converged=True | mu=-0.32802883 Ha | Zbar=4.07907607 | time=11.681 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -7.91966 | 0.96068
[ion] HNC iterations=6 QOZ=0.799 s HNC=0.011 s
[done] total time=12.501 s
[accepted] c_64p64_qm
Otter 0.4.1
[run] C | rho=12.64 g/cm^3 | Te=64.64 eV | Ti=64.64 eV | TF -> QOZ/HNC
[AA/TF] C (Z=6) Te=64.64 eV rho=12.64 g/cc Rws=1.364781 Bohr run_mode=full+ext
[AA/TF full SCF]
iter | d_n | d_v
1 | inf | 7.122e-01
10 | 3.345e-03 | 8.808e-02
20 | 1.441e-05 | 8.365e-05
[AA/TF external SCF]
iter | d_n | d_v
1 | inf | 2.077e+01
10 | 2.064e-04 | 1.356e-02
[AA/TF] converged=True mu=-0.19562887 Ha Zbar=3.54769601 time=0.208 s
[ion] HNC iterations=6 QOZ=0.785 s HNC=0.011 s
[done] total time=1.011 s
[accepted] c_64p64_tf
Otter 0.4.1
[run] W | rho=40 g/cm^3 | Te=10 eV | Ti=10 eV | QM -> QOZ/HNC
[AA] W (Z=74) | Rws=2.308053 a0 | full + external
[AA/full SCF: stage 1]
iter | d_n | d_v
0 | nan | 6.062e-02
1 | 1.306e-02 | 5.399e-02
2 | 1.650e-02 | 4.785e-02
3 | 1.399e-02 | 4.206e-02
4 | 1.166e-02 | 3.667e-02
[AA/full SCF: stage 2]
iter | d_n | d_v
0 | nan | 3.172e-02
1 | 7.813e-03 | 2.724e-02
2 | 6.285e-03 | 2.324e-02
3 | 4.996e-03 | 1.970e-02
4 | 3.924e-03 | 1.662e-02
5 | 3.047e-03 | 1.396e-02
6 | 2.340e-03 | 1.171e-02
7 | 1.778e-03 | 9.847e-03
8 | 1.341e-03 | 8.376e-03
9 | 1.009e-03 | 7.322e-03
10 | 7.596e-04 | 6.742e-03
11 | 5.754e-04 | 6.776e-03
12 | 4.492e-04 | 7.864e-03
13 | 4.059e-04 | 9.774e-03
14 | 4.691e-04 | 6.420e-03
15 | 6.670e-04 | 3.991e-03
16 | 2.194e-04 | 2.780e-03
17 | 9.865e-04 | 7.737e-04
18 | 2.016e-04 | 1.974e-03
19 | 1.145e-03 | 2.854e-04
20 | 1.127e-04 | 5.327e-04
21 | 3.975e-04 | 4.750e-04
22 | 1.512e-04 | 7.551e-05
23 | 3.684e-05 | 6.570e-05
24 | 4.196e-05 | 2.949e-05
25 | 2.063e-05 | 5.353e-06
26 | 1.976e-06 | 6.745e-06
27 | 4.566e-06 | 1.325e-06
28 | 1.132e-06 | 3.341e-07
[AA/external SCF]
iter | d_n | d_v
0 | 9.753e-08 | 6.143e+12
1 | 1.187e-01 | 1.117e+00
2 | 1.295e-01 | 7.608e-01
3 | 3.756e-02 | 4.065e-01
4 | 1.061e-02 | 1.882e-01
5 | 5.832e-03 | 7.996e-02
6 | 3.000e-03 | 3.246e-02
7 | 1.391e-03 | 4.656e-03
8 | 1.793e-04 | 2.633e-03
9 | 5.149e-05 | 1.682e-03
10 | 3.495e-05 | 6.723e-04
11 | 1.422e-05 | 3.440e-04
12 | 5.870e-06 | 2.658e-04
13 | 5.164e-06 | 1.178e-04
14 | 4.912e-06 | 5.992e-05
15 | 1.496e-06 | 4.180e-05
16 | 8.180e-07 | 1.715e-06
17 | 5.738e-08 | 5.299e-05
18 | 6.293e-07 | 3.093e-05
19 | 4.753e-07 | 2.531e-05
20 | 5.994e-07 | 2.055e-06
21 | 4.465e-08 | 4.987e-07
22 | 8.854e-09 | 2.282e-07
[AA] converged=True | mu=1.03907760 Ha | Zbar=9.38617324 | time=73.629 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -2336.88203 | 1.00000
2s | -383.28517 | 1.00000
2p | -367.75534 | 1.00000
3s | -86.58219 | 1.00000
3p | -79.28425 | 1.00000
3d | -65.50102 | 1.00000
4s | -16.33083 | 1.00000
4p | -13.27103 | 1.00000
4d | -7.65390 | 1.00000
4f | -0.32564 | 0.97619
5s | -1.18266 | 0.99764
5p | -0.33966 | 0.97706
[ion] HNC iterations=13 QOZ=0.796 s HNC=0.022 s
[done] total time=74.460 s
[accepted] w_10_qm
Otter 0.4.1
[run] W | rho=40 g/cm^3 | Te=10 eV | Ti=10 eV | TF -> QOZ/HNC
[AA/TF] W (Z=74) Te=10 eV rho=40 g/cc Rws=2.308053 Bohr run_mode=full+ext
[AA/TF full SCF]
iter | d_n | d_v
1 | inf | 8.021e-01
10 | 8.401e-01 | 7.057e-02
20 | 8.634e-01 | 4.549e-01
30 | 1.796e-02 | 1.588e-03
40 | 6.004e-04 | 2.760e-04
50 | 1.482e-04 | 1.330e-04
60 | 2.284e-04 | 1.401e-04
70 | 1.084e-04 | 5.399e-04
80 | 1.452e-04 | 7.872e-04
90 | 3.136e-04 | 3.278e-04
100 | 2.039e-04 | 8.855e-04
110 | 4.677e-05 | 3.559e-04
120 | 1.202e-04 | 1.710e-04
130 | 3.194e-05 | 2.624e-04
140 | 1.064e-04 | 5.220e-04
150 | 7.208e-05 | 6.981e-04
160 | 2.046e-04 | 3.674e-04
170 | 3.494e-04 | 4.393e-04
180 | 1.830e-04 | 1.114e-03
190 | 1.536e-04 | 3.571e-04
200 | 2.313e-04 | 1.318e-04
210 | 1.710e-04 | 1.827e-04
220 | 1.306e-04 | 1.088e-03
230 | 2.262e-05 | 3.454e-04
240 | 1.892e-05 | 4.700e-04
250 | 2.351e-04 | 6.464e-04
260 | 2.194e-04 | 2.778e-04
270 | 2.608e-05 | 2.937e-04
280 | 1.606e-04 | 3.678e-04
290 | 2.643e-04 | 4.869e-04
300 | 1.437e-04 | 7.767e-04
310 | 7.794e-05 | 4.604e-04
320 | 1.464e-04 | 8.771e-04
330 | 3.399e-04 | 2.446e-04
340 | 2.313e-04 | 7.703e-05
350 | 1.618e-04 | 9.423e-05
360 | 9.661e-05 | 3.517e-04
370 | 1.498e-04 | 8.050e-04
380 | 3.012e-04 | 4.823e-04
390 | 8.485e-05 | 1.238e-03
400 | 1.733e-04 | 2.551e-04
[AA/TF external SCF]
iter | d_n | d_v
1 | inf | 1.305e+01
10 | 9.262e-03 | 9.099e-01
20 | 3.084e-05 | 1.014e-04
[AA/TF] converged=True mu=1.25150578 Ha Zbar=7.67562725 time=1.087 s
[ion] HNC iterations=10 QOZ=0.796 s HNC=0.484 s
[done] total time=2.373 s
[accepted] w_10_tf
Otter 0.4.1
[run] W | rho=40 g/cm^3 | Te=60 eV | Ti=60 eV | QM -> QOZ/HNC
[AA] W (Z=74) | Rws=2.308053 a0 | full + external
[AA/full SCF: stage 1]
iter | d_n | d_v
0 | nan | 6.062e-02
1 | 1.629e-02 | 5.410e-02
2 | 1.497e-02 | 4.790e-02
3 | 1.323e-02 | 4.207e-02
4 | 1.126e-02 | 3.665e-02
[AA/full SCF: stage 2]
iter | d_n | d_v
0 | nan | 3.167e-02
1 | 7.677e-03 | 2.715e-02
2 | 6.198e-03 | 2.309e-02
3 | 4.934e-03 | 1.952e-02
4 | 3.869e-03 | 1.647e-02
5 | 2.988e-03 | 1.402e-02
6 | 2.291e-03 | 1.224e-02
7 | 1.775e-03 | 1.112e-02
8 | 1.430e-03 | 1.058e-02
9 | 1.242e-03 | 1.043e-02
10 | 1.166e-03 | 1.015e-02
11 | 1.097e-03 | 8.470e-03
12 | 5.185e-04 | 3.869e-03
13 | 2.487e-04 | 2.015e-03
14 | 4.104e-04 | 2.038e-03
15 | 7.524e-04 | 1.563e-03
16 | 7.188e-04 | 1.908e-03
17 | 7.191e-04 | 5.152e-04
18 | 1.706e-04 | 3.435e-04
19 | 1.467e-04 | 4.634e-05
20 | 7.086e-06 | 3.579e-05
21 | 2.084e-05 | 1.548e-05
22 | 8.199e-06 | 2.245e-06
23 | 2.083e-06 | 3.491e-06
24 | 3.362e-06 | 1.280e-06
[AA/external SCF]
iter | d_n | d_v
0 | 2.320e-06 | 5.434e+12
1 | 4.427e-02 | 1.190e+00
2 | 4.767e-02 | 1.356e+00
3 | 3.464e-02 | 3.246e-01
4 | 8.039e-03 | 9.038e-02
5 | 2.414e-03 | 2.202e-02
6 | 7.801e-04 | 6.052e-03
7 | 1.661e-04 | 1.921e-03
8 | 7.569e-05 | 8.129e-04
9 | 1.776e-05 | 4.295e-04
10 | 8.326e-06 | 7.709e-05
11 | 1.614e-06 | 2.017e-05
12 | 5.883e-07 | 3.677e-05
13 | 7.111e-07 | 8.242e-06
14 | 2.360e-07 | 2.653e-06
15 | 5.062e-08 | 3.992e-07
16 | 7.946e-09 | 4.433e-07
[AA] converged=True | mu=-1.12786972 Ha | Zbar=14.30080681 | time=76.891 s
Bound levels
level | E [Ha] | f_FD
------+------------+----------
1s | -2339.79054 | 1.00000
2s | -386.27184 | 1.00000
2p | -370.72064 | 1.00000
3s | -89.69117 | 1.00000
3p | -82.39004 | 1.00000
3d | -68.59507 | 1.00000
4s | -18.99627 | 0.99970
4p | -15.91637 | 0.99878
4d | -10.25720 | 0.98433
4f | -2.82358 | 0.68331
5s | -2.29714 | 0.62955
5p | -1.27666 | 0.51686
[ion] HNC iterations=11 QOZ=0.788 s HNC=0.019 s
[done] total time=77.715 s
Otter 0.4.1
[run] W | rho=40 g/cm^3 | Te=60 eV | Ti=60 eV | TF -> QOZ/HNC
[accepted] w_60_qm
[AA/TF] W (Z=74) Te=60 eV rho=40 g/cc Rws=2.308053 Bohr run_mode=full+ext
[AA/TF full SCF]
iter | d_n | d_v
1 | inf | 8.017e-01
10 | 4.325e-01 | 1.556e-01
20 | 8.864e-04 | 1.873e-04
30 | 5.570e-05 | 3.920e-05
[AA/TF external SCF]
iter | d_n | d_v
1 | inf | 1.080e+01
10 | 3.428e-03 | 5.974e-02
[AA/TF] converged=True mu=-0.99677085 Ha Zbar=11.31855997 time=0.276 s
[ion] HNC iterations=5 QOZ=0.796 s HNC=0.201 s
[done] total time=1.280 s
[accepted] w_60_tf
fe_10_qm Zbar=8.28279690 HNC=3.419e-05 closure=1.999e-03
fe_10_tf Zbar=8.80781669 HNC=7.197e-05 closure=3.089e-04
h_5_qm Zbar=1.00000000 HNC=7.958e-06 closure=5.123e-04
h_5_tf Zbar=0.97263161 HNC=7.466e-05 closure=9.122e-04
h_172_qm Zbar=1.00000000 HNC=7.771e-06 closure=1.121e-05
h_172_tf Zbar=0.97422595 HNC=5.635e-06 closure=7.893e-06
c_64p64_qm Zbar=4.07902606 HNC=5.860e-06 closure=8.458e-05
c_64p64_tf Zbar=4.16575793 HNC=7.906e-06 closure=1.424e-04
w_10_qm Zbar=9.33008924 HNC=6.624e-05 closure=3.553e-04
w_10_tf Zbar=12.73137517 HNC=4.357e-09 closure=8.342e-07
w_60_qm Zbar=14.29310823 HNC=9.574e-05 closure=7.486e-04
w_60_tf Zbar=15.28524162 HNC=1.990e-10 closure=5.361e-11
from __future__ import annotations
from concurrent.futures import ProcessPoolExecutor, as_completed
import hashlib
import json
import os
from pathlib import Path
import subprocess
import time
from typing import Any
import matplotlib.pyplot as plt
import numpy as np
from otter import PlasmaWorkflowConfig, solve_plasma_workflow
from otter.numerics.constants import (
ATOMIC_MASS_UNIT_TO_G,
BOHR_TO_ANGSTROM,
BOHR_TO_CM,
)
from otter.plotting import grid_figsize, save_figure, set_style
# =============================================================================
# User input
# =============================================================================
USE_PRECOMPUTED_DATA = False
# Serial states, each with one AA worker, also on memory-constrained hosts.
MAX_STATE_WORKERS = 1
_QOZ_POINTS_ENV = os.environ.get("OTTER_STARRETT_SINGLE_QOZ_POINTS")
QOZ_N_POINTS_OVERRIDE = (
None if _QOZ_POINTS_ENV is None else int(_QOZ_POINTS_ENV)
)
HNC_CLOSURE_TOL = 2.5e-3
R_RETAIN_MAX_BOHR = 20.0
K_RETAIN_MAX_BOHR_INV = 20.0
# An empty tuple means all state/model combinations. Users may list stable
# IDs from CALCULATIONS below for a shorter live calculation.
STATE_IDS_TO_COMPUTE: tuple[str, ...] = ()
# Maintenance override; normal users only need to edit USE_PRECOMPUTED_DATA.
if os.environ.get("OTTER_RECOMPUTE_STARRETT_SINGLE", "0") == "1":
USE_PRECOMPUTED_DATA = False
if selected_ids := os.environ.get("OTTER_STARRETT_SINGLE_STATE_IDS"):
STATE_IDS_TO_COMPUTE = tuple(
item.strip() for item in selected_ids.split(",") if item.strip()
)
# =============================================================================
BENCHMARK_ID = "starrett_single_species_2013_2014"
SCHEMA = "otter_starrett_single_species_state_v1"
PHYSICAL_STATES: tuple[dict[str, Any], ...] = (
{
"panel_id": "fe_10",
"element": "Fe",
"atomic_mass": 55.845,
"rho_g_cc": 22.5,
"te_ev": 10.0,
"ti_ev": 10.0,
"title": r"Fe, 22.5 g cm$^{-3}$, 10 eV",
"x_max": 5.0,
},
{
"panel_id": "h_5",
"element": "H",
"atomic_mass": 1.008,
"rho_g_cc": 80.0,
"te_ev": 5.0,
"ti_ev": 5.0,
"title": r"H, 80 g cm$^{-3}$, 5 eV",
"x_max": 6.2,
},
{
"panel_id": "h_172",
"element": "H",
"atomic_mass": 1.008,
"rho_g_cc": 80.0,
"te_ev": 172.0,
"ti_ev": 172.0,
"title": r"H, 80 g cm$^{-3}$, 172 eV",
"x_max": 6.2,
},
{
"panel_id": "c_64p64",
"element": "C",
"atomic_mass": 12.011,
"rho_g_cc": 12.64,
"te_ev": 64.64,
"ti_ev": 64.64,
"title": r"C, 12.64 g cm$^{-3}$, 64.64 eV",
"x_max": 4.0,
},
{
"panel_id": "w_10",
"element": "W",
"atomic_mass": 183.84,
"rho_g_cc": 40.0,
"te_ev": 10.0,
"ti_ev": 10.0,
"title": r"W, 40 g cm$^{-3}$, 10 eV",
"x_max": 6.2,
},
{
"panel_id": "w_60",
"element": "W",
"atomic_mass": 183.84,
"rho_g_cc": 40.0,
"te_ev": 60.0,
"ti_ev": 60.0,
"title": r"W, 40 g cm$^{-3}$, 60 eV",
"x_max": 6.2,
},
)
CALCULATIONS: tuple[dict[str, Any], ...] = tuple(
{
**state,
"model": model,
"state_id": f"{state['panel_id']}_{model}",
}
for state in PHYSICAL_STATES
for model in ("qm", "tf")
)
REFERENCE_SERIES: dict[str, tuple[dict[str, str], ...]] = {
"fe_10": (
{
"file": "gii_Fe_22.5gcc_10.0ev_starrett.csv",
"label": "Starrett TF",
"coordinate": "r_over_rws",
},
),
"h_5": (
{
"file": "gii_H_80gcc_5.0ev_starrett.csv",
"label": "Starrett TF",
"coordinate": "bohr",
},
),
"h_172": (
{
"file": "gii_H_80gcc_172.0ev_starrett.csv",
"label": "Starrett TF",
"coordinate": "bohr",
},
),
"c_64p64": (
{
"file": "gii_C_12.64gcc_64.64ev_starrett.csv",
"label": "Starrett TF",
"coordinate": "r_over_rws",
},
),
"w_10": (
{
"file": "gii_W_40gcc_10.0ev_starrett_TF.csv",
"label": "Starrett TF",
"coordinate": "bohr",
},
),
"w_60": (
{
"file": "gii_W_40gcc_60.0ev_starrett.csv",
"label": "Starrett TF",
"coordinate": "bohr",
},
),
}
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" / BENCHMARK_ID
REFERENCE_DIR = ROOT / "benchmarks" / "reference_data" / BENCHMARK_ID
OUTPUT_DIR = ROOT / "benchmarks" / "outputs" / BENCHMARK_ID
CANDIDATE_DIR = OUTPUT_DIR / "gallery_recomputed"
FIGURE_DIR = OUTPUT_DIR / "figures"
CONTROLLER_PATH = (
ROOT
/ "benchmarks"
/ "examples"
/ "plot_starrett_single_species_2013_2014.py"
)
_REFERENCE_RECORDS_CACHE: dict[str, dict[str, Any]] | None = None
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 only portable, pickle-free arrays."""
with np.load(path, allow_pickle=False) as archive:
result = {key: np.asarray(archive[key]) for key in archive.files}
if any(value.dtype.hasobject for value in result.values()):
raise TypeError(f"Object arrays are forbidden in {path}.")
return result
def reference_records_by_file() -> dict[str, dict[str, Any]]:
"""Verify the attributed reference manifest and return file records."""
global _REFERENCE_RECORDS_CACHE
if _REFERENCE_RECORDS_CACHE is not None:
return _REFERENCE_RECORDS_CACHE
manifest_path = REFERENCE_DIR / "manifest.json"
manifest = json.loads(manifest_path.read_text(encoding="utf-8"))
if manifest.get("schema_version") != "otter_reference_manifest_v1":
raise ValueError(f"Unexpected reference schema in {manifest_path}.")
if manifest.get("reference_id") != BENCHMARK_ID:
raise ValueError(f"Unexpected reference identifier in {manifest_path}.")
if manifest.get("redistribution_status") != (
"published_by_maintainer_with_attribution"
):
raise ValueError("Reference-data release status was not reviewed.")
records = {
str(record["path"]): dict(record) for record in manifest["files"]
}
expected = {
str(record["file"])
for panel_records in REFERENCE_SERIES.values()
for record in panel_records
}
if not expected.issubset(records):
raise RuntimeError(
"Reference manifest lacks a curve listed in REFERENCE_SERIES."
)
for filename, record in records.items():
path = REFERENCE_DIR / filename
if not path.is_file():
raise FileNotFoundError(path)
if sha256_file(path) != str(record["sha256"]):
raise RuntimeError(f"Reference checksum mismatch for {path}.")
_REFERENCE_RECORDS_CACHE = records
return records
def validated_two_column_curve(path: Path) -> np.ndarray:
"""Load one finite two-column curve and return it in radius order.
The inherited W 60 eV ``60.0ev`` extraction contains one adjacent
coordinate inversion. The source CSV remains byte-for-byte unchanged;
plotting applies a stable coordinate sort after rejecting duplicate
radii.
"""
raw = np.asarray(np.loadtxt(path, delimiter=","), dtype=float)
if raw.ndim != 2 or raw.shape[1] != 2 or raw.shape[0] < 2:
raise ValueError(f"{path} must contain an N x 2 curve with N >= 2.")
if not np.all(np.isfinite(raw)):
raise ValueError(f"{path} contains non-finite values.")
if np.unique(raw[:, 0]).size != raw.shape[0]:
raise ValueError(f"{path} contains duplicate radius coordinates.")
order = np.argsort(raw[:, 0], kind="stable")
ordered = raw[order]
if not np.all(np.diff(ordered[:, 0]) > 0.0):
raise ValueError(f"{path} radius coordinates cannot be ordered.")
return ordered
def ion_sphere_radius_bohr(state: dict[str, Any]) -> float:
"""Return R_WS from mass density and the standard atomic mass."""
ion_density_cm3 = float(state["rho_g_cc"]) / (
float(state["atomic_mass"]) * ATOMIC_MASS_UNIT_TO_G
)
radius_cm = (3.0 / (4.0 * np.pi * ion_density_cm3)) ** (1.0 / 3.0)
return float(radius_cm / BOHR_TO_CM)
def load_reference(panel_id: str) -> list[dict[str, Any]]:
"""Load raw publication coordinates and map only at plot time."""
manifest_records = reference_records_by_file()
state = next(
item for item in PHYSICAL_STATES if item["panel_id"] == panel_id
)
rws = ion_sphere_radius_bohr(state)
loaded: list[dict[str, Any]] = []
for record in REFERENCE_SERIES[panel_id]:
filename = str(record["file"])
manifest_record = manifest_records[filename]
raw = validated_two_column_curve(REFERENCE_DIR / filename)
coordinate = str(record["coordinate"])
if coordinate == "bohr":
if manifest_record.get("column_1_unit") != "Bohr":
raise ValueError(f"Unit mismatch for {filename}.")
x = raw[:, 0] / rws
elif coordinate == "r_over_rws":
if (
manifest_record.get("column_1_unit")
!= "r_over_R_WS_dimensionless"
):
raise ValueError(f"Unit mismatch for {filename}.")
x = raw[:, 0]
else:
raise ValueError(f"Unknown coordinate convention: {coordinate}")
loaded.append(
{
**record,
"r_over_rws": np.asarray(x, dtype=float),
"gii": np.asarray(raw[:, 1], dtype=float),
}
)
return loaded
def validate_baseline_payload(
payload: dict[str, np.ndarray],
state: dict[str, Any],
record: dict[str, Any],
path: Path,
) -> None:
"""Validate scientific metadata and arrays in one compact baseline."""
required = {
"schema_version",
"state_id",
"panel_id",
"element",
"electronic_model",
"rho_g_cc",
"te_ev",
"ti_ev",
"r_ws_bohr",
"r_bohr",
"gii_r",
"k_bohr_inv",
"sii_k",
"zbar_partition",
"hnc_output_residual",
"closure_transform_max_abs",
"threshold_state_status",
}
missing = required - set(payload)
if missing:
raise ValueError(f"{path} is missing fields: {sorted(missing)}")
expected_strings = {
"schema_version": SCHEMA,
"state_id": str(state["state_id"]),
"panel_id": str(state["panel_id"]),
"element": str(state["element"]),
"electronic_model": str(state["model"]),
}
for key, expected in expected_strings.items():
value = np.asarray(payload[key])
if value.shape != () or str(value.item()) != expected:
raise ValueError(f"{path}: invalid {key}.")
expected_scalars = {
"rho_g_cc": float(state["rho_g_cc"]),
"te_ev": float(state["te_ev"]),
"ti_ev": float(state["ti_ev"]),
}
for key, expected in expected_scalars.items():
value = np.asarray(payload[key])
if value.shape != () or not np.isfinite(float(value)):
raise ValueError(f"{path}: invalid {key}.")
if not np.isclose(float(value), expected, rtol=0.0, atol=1.0e-12):
raise ValueError(f"{path}: metadata mismatch for {key}.")
r_ws = float(np.asarray(payload["r_ws_bohr"]))
expected_r_ws = ion_sphere_radius_bohr(state)
if not np.isfinite(r_ws) or not np.isclose(
r_ws, expected_r_ws, rtol=1.0e-8, atol=1.0e-10
):
raise ValueError(f"{path}: inconsistent ion-sphere radius.")
r = np.asarray(payload["r_bohr"], dtype=float)
g = np.asarray(payload["gii_r"], dtype=float)
k = np.asarray(payload["k_bohr_inv"], dtype=float)
s = np.asarray(payload["sii_k"], dtype=float)
for name, grid, values in (
("real-space", r, g),
("reciprocal-space", k, s),
):
if grid.ndim != 1 or values.ndim != 1 or grid.size < 2:
raise ValueError(f"{path}: invalid {name} array dimensions.")
if grid.shape != values.shape:
raise ValueError(f"{path}: mismatched {name} arrays.")
if not np.all(np.isfinite(grid)) or not np.all(np.isfinite(values)):
raise ValueError(f"{path}: non-finite {name} data.")
if not np.all(np.diff(grid) > 0.0) or grid[0] <= 0.0:
raise ValueError(f"{path}: {name} grid is not positive/ordered.")
if r[-1] > R_RETAIN_MAX_BOHR + 1.0e-10:
raise ValueError(f"{path}: real-space data exceed the retention limit.")
if k[-1] > K_RETAIN_MAX_BOHR_INV + 1.0e-10:
raise ValueError(
f"{path}: reciprocal-space data exceed the retention limit."
)
residual = float(np.asarray(payload["hnc_output_residual"]))
closure = float(np.asarray(payload["closure_transform_max_abs"]))
zbar = float(np.asarray(payload["zbar_partition"]))
threshold = str(np.asarray(payload["threshold_state_status"]).item()).lower()
if not np.isfinite(zbar) or zbar <= 0.0:
raise ValueError(f"{path}: invalid pseudoatom-partition Zbar.")
configured_hnc_tol = float(workflow_config(state).hnc_tol)
if not np.isfinite(residual) or residual > configured_hnc_tol:
raise ValueError(f"{path}: HNC residual does not pass the strict gate.")
if not np.isfinite(closure) or closure > HNC_CLOSURE_TOL:
raise ValueError(f"{path}: transform closure does not pass the gate.")
if threshold == "unresolved":
raise ValueError(f"{path}: unresolved threshold representation.")
manifest_scalars = {
"zbar_partition": zbar,
"hnc_output_residual": residual,
"closure_transform_max_abs": closure,
}
for key, payload_value in manifest_scalars.items():
if not np.isclose(
float(record[key]), payload_value, rtol=1.0e-12, atol=1.0e-14
):
raise ValueError(f"{path}: manifest mismatch for {key}.")
def load_precomputed_states() -> dict[str, dict[str, np.ndarray]]:
"""Verify the manifest and load all strictly accepted states."""
manifest = json.loads(
(BASELINE_DIR / "manifest.json").read_text(encoding="utf-8")
)
if manifest.get("benchmark_id") != BENCHMARK_ID:
raise ValueError("Unexpected baseline manifest.")
controller_hash = manifest.get("producer", {}).get(
"current_controller_sha256"
)
# This opt-in reader displays a historical run. Its producer hash need
# not equal the current plotting controller; array hashes and physical
# state/convergence invariants are still checked below.
if not isinstance(controller_hash, str) or len(controller_hash) != 64:
raise RuntimeError("Missing recorded baseline controller fingerprint.")
expected = {str(item["state_id"]) for item in CALCULATIONS}
definitions = {str(item["state_id"]): item for item in CALCULATIONS}
records = {
str(record["state_id"]): record for record in manifest["states"]
}
if set(records) != expected:
raise RuntimeError("Baseline manifest state coverage is incomplete.")
terminal = {
"accepted",
"strict_aa_rejected",
"strict_hnc_rejected",
"not_calculated",
}
bad = {
key: record.get("status")
for key, record in records.items()
if record.get("status") not in terminal
}
if bad:
raise RuntimeError(f"Unreviewed baseline records: {bad}")
loaded: dict[str, dict[str, np.ndarray]] = {}
for state_id, record in records.items():
if record["status"] != "accepted":
if record.get("baseline_file") is not None:
raise RuntimeError(
f"Rejected state {state_id} names a baseline file."
)
continue
path = BASELINE_DIR / str(record["baseline_file"])
if sha256_file(path) != str(record["baseline_sha256"]):
raise RuntimeError(f"Checksum mismatch for {path}.")
payload = load_npz(path)
validate_baseline_payload(
payload, definitions[state_id], record, path
)
loaded[state_id] = payload
return loaded
def workflow_config(state: dict[str, Any]) -> PlasmaWorkflowConfig:
"""Build one complete public Otter IS-QOZ/HNC calculation."""
model = str(state["model"])
model_override = {} if model == "qm" else {"electronic_model": model}
qoz_override = (
{}
if QOZ_N_POINTS_OVERRIDE is None
else {"qoz_linear_n_points": QOZ_N_POINTS_OVERRIDE}
)
return PlasmaWorkflowConfig(
elements=[str(state["element"])],
temperature_ev=float(state["te_ev"]),
ion_temperature_ev=float(state["ti_ev"]),
rho_g_cc=float(state["rho_g_cc"]),
**model_override,
**qoz_override,
hnc_closure_transform_tol=float(HNC_CLOSURE_TOL),
)
def pack_result(
workflow: dict[str, Any],
state: dict[str, Any],
elapsed_s: float,
) -> dict[str, np.ndarray]:
"""Apply strict gates and retain compact numerical output."""
electronic = dict(workflow["electronic"]["result"])
ion = dict(workflow["ion"])
if electronic.get("stage2_converged") is not True:
raise RuntimeError("full average-atom stage 2 did not converge")
if dict(electronic.get("ext_status", {})).get("converged") is not True:
raise RuntimeError("external fixed-mu average atom did not converge")
threshold = str(
electronic.get("threshold_state_status", "")
).lower()
if threshold == "unresolved":
raise RuntimeError("threshold-state representation is unresolved")
if ion.get("hnc_converged") is not True:
raise RuntimeError("HNC did not reach a physical fixed point")
residual = float(ion["hnc_output_residual"])
closure = float(ion["closure_transform_max_abs"])
configured_hnc_tol = float(workflow["configuration"]["hnc_tol"])
if residual > configured_hnc_tol:
raise RuntimeError(
f"HNC residual {residual:.6e} exceeds {configured_hnc_tol:.6e}"
)
if closure > HNC_CLOSURE_TOL:
raise RuntimeError(
f"transform closure {closure:.6e} exceeds "
f"{HNC_CLOSURE_TOL:.6e}"
)
r = np.asarray(ion["r"], dtype=float)
k = np.asarray(ion["k"], dtype=float)
r_mask = r <= R_RETAIN_MAX_BOHR
k_mask = k <= K_RETAIN_MAX_BOHR_INV
rws = float(electronic["r_ws"])
return {
"schema_version": np.asarray(SCHEMA),
"storage_profile": np.asarray("benchmark_analysis"),
"state_id": np.asarray(str(state["state_id"])),
"panel_id": np.asarray(str(state["panel_id"])),
"element": np.asarray(str(state["element"])),
"electronic_model": np.asarray(str(state["model"])),
"rho_g_cc": np.asarray(float(state["rho_g_cc"])),
"te_ev": np.asarray(float(state["te_ev"])),
"ti_ev": np.asarray(float(state["ti_ev"])),
"producer_elapsed_s": np.asarray(float(elapsed_s)),
"r_ws_bohr": np.asarray(rws),
"r_bohr": r[r_mask],
"gii_r": np.asarray(ion["gij_r"][0, 0], dtype=float)[r_mask],
"k_bohr_inv": k[k_mask],
"sii_k": np.asarray(ion["sij_k"][0, 0], dtype=float)[k_mask],
"zbar_partition": np.asarray(np.asarray(ion["zbar_partition"]).item()),
"hnc_output_residual": np.asarray(residual),
"closure_transform_max_abs": np.asarray(closure),
"threshold_state_status": np.asarray(threshold),
}
def solve_state(state: dict[str, Any]) -> dict[str, np.ndarray]:
"""Calculate one state/model in a worker process."""
started = time.perf_counter()
workflow = solve_plasma_workflow(workflow_config(state))
return pack_result(workflow, state, time.perf_counter() - started)
def git_revision() -> dict[str, Any]:
"""Record the local revision without claiming dirty work is committed."""
try:
head = subprocess.run(
["git", "rev-parse", "HEAD"],
cwd=ROOT,
check=True,
capture_output=True,
text=True,
).stdout.strip()
dirty = bool(
subprocess.run(
["git", "status", "--porcelain"],
cwd=ROOT,
check=True,
capture_output=True,
text=True,
).stdout.strip()
)
except (OSError, subprocess.CalledProcessError):
return {"git_commit": None, "worktree_dirty": None}
return {"git_commit": head, "worktree_dirty": dirty}
def save_candidates(
solved: dict[str, dict[str, np.ndarray]],
failures: dict[str, str],
attempted: set[str],
) -> None:
"""Save accepted candidates and explicit non-promoted failures."""
CANDIDATE_DIR.mkdir(parents=True, exist_ok=True)
records: list[dict[str, Any]] = []
for state in CALCULATIONS:
state_id = str(state["state_id"])
if state_id not in attempted:
records.append(
{
"state_id": state_id,
"status": "not_attempted_by_state_filter",
"candidate_file": None,
"candidate_sha256": None,
}
)
elif state_id in failures:
records.append(
{
"state_id": state_id,
"status": "strict_calculation_rejected",
"candidate_file": None,
"candidate_sha256": None,
"reason": failures[state_id],
}
)
else:
filename = f"{state_id}.npz"
path = CANDIDATE_DIR / filename
np.savez_compressed(path, **solved[state_id])
records.append(
{
"state_id": state_id,
"status": "candidate_passed_strict_gates",
"candidate_file": filename,
"candidate_sha256": sha256_file(path),
}
)
controller = CONTROLLER_PATH
manifest = {
"schema_version": "otter_candidate_manifest_v1",
"benchmark_id": BENCHMARK_ID,
"producer": {
"script_relative_path": str(controller.relative_to(ROOT)),
"script_sha256": sha256_file(controller),
**git_revision(),
},
"states": records,
}
(CANDIDATE_DIR / "manifest.json").write_text(
json.dumps(manifest, indent=2) + "\n",
encoding="utf-8",
)
def solve_selected_states() -> dict[str, dict[str, np.ndarray]]:
"""Run selected state/models with bounded process-level parallelism."""
selected = set(STATE_IDS_TO_COMPUTE)
jobs = [
state
for state in CALCULATIONS
if not selected or str(state["state_id"]) in selected
]
unknown = selected - {str(item["state_id"]) for item in CALCULATIONS}
if unknown:
raise ValueError(f"Unknown STATE_IDS_TO_COMPUTE: {sorted(unknown)}")
solved: dict[str, dict[str, np.ndarray]] = {}
failures: dict[str, str] = {}
with ProcessPoolExecutor(max_workers=MAX_STATE_WORKERS) as pool:
future_map = {pool.submit(solve_state, state): state for state in jobs}
for future in as_completed(future_map):
state = future_map[future]
state_id = str(state["state_id"])
try:
solved[state_id] = future.result()
print(f"[accepted] {state_id}")
except Exception as exc: # scientific rejection is recorded
failures[state_id] = f"{type(exc).__name__}: {exc}"
print(f"[rejected] {state_id}: {failures[state_id]}")
save_candidates(
solved,
failures,
{str(state["state_id"]) for state in jobs},
)
if failures:
raise RuntimeError(f"Incomplete Starrett scan: {failures}; successful candidates retained")
return solved
def main() -> None:
states = (
load_precomputed_states()
if USE_PRECOMPUTED_DATA
else solve_selected_states()
)
print(
"Using "
+ (
"checksummed, precomputed Otter states."
if USE_PRECOMPUTED_DATA
else "new states calculated directly by this gallery script."
)
)
for definition in CALCULATIONS:
state_id = str(definition["state_id"])
if state_id not in states:
continue
payload = states[state_id]
print(
f"{state_id:18s} "
f"Zbar={float(payload['zbar_partition']):.8f} "
f"HNC={float(payload['hnc_output_residual']):.3e} "
f"closure={float(payload['closure_transform_max_abs']):.3e}"
)
# %%
# All panels use :math:`r/R_{\rm WS}`. Explicitly Bohr-valued reference
# coordinates are converted at plot time without modifying the stored CSV
# files.
set_style("thesis", palette="bing")
fig, axes = plt.subplots(
2,
3,
figsize=grid_figsize(2, 3),
squeeze=False,
)
axes_flat = axes.ravel()
reference_colors = ("#333333", "#666666")
model_styles = {
"qm": {
"label": "Otter QM",
"color": "#0072B2",
"ls": "-",
},
"tf": {"label": "Otter TF", "color": "#D55E00", "ls": "--"},
}
for axis, physical in zip(axes_flat, PHYSICAL_STATES, strict=True):
panel_id = str(physical["panel_id"])
for index, reference in enumerate(load_reference(panel_id)):
axis.plot(
reference["r_over_rws"],
reference["gii"],
ls="none",
marker="o",
ms=5.4,
mfc="none",
mew=1.15,
color=reference_colors[index % len(reference_colors)],
label=str(reference["label"]),
)
for model in ("qm", "tf"):
state_id = f"{panel_id}_{model}"
if state_id not in states:
continue
payload = states[state_id]
style = model_styles[model]
axis.plot(
np.asarray(payload["r_bohr"], dtype=float)
/ float(payload["r_ws_bohr"]),
np.asarray(payload["gii_r"], dtype=float),
color=str(style["color"]),
ls=str(style["ls"]),
lw=2.0,
label=str(style["label"]),
)
missing_models = [
model
for model in ("qm", "tf")
if f"{panel_id}_{model}" not in states
]
if missing_models:
axis.text(
0.03,
0.05,
"No accepted Otter "
+ "/".join(model.upper() for model in missing_models)
+ " result",
transform=axis.transAxes,
fontsize="small",
color="0.35",
)
axis.axhline(1.0, color="#777777", ls=":", lw=0.8)
axis.set_xlim(-0.5, float(physical["x_max"]))
axis.set_ylim(-0.04, 1.86)
axis.set_title(str(physical["title"]))
axis.set_xlabel(r"$r/R_{\rm WS}$")
axis.set_ylabel(r"$g_{ii}(r)$")
axis.legend(fontsize="small", loc="best")
fig.suptitle(
"Single-species ion structure: Otter and Starrett–Saumon references",
y=0.985,
)
fig.text(
0.5,
0.006,
"Reference data: Starrett and Saumon (2014), "
"doi:10.1016/j.hedp.2013.12.001.",
ha="center",
va="bottom",
fontsize=7.5,
)
fig.tight_layout(rect=(0.0, 0.035, 1.0, 0.965), pad=0.55)
saved = save_figure(
fig,
FIGURE_DIR / "starrett_single_species_2013_2014",
formats=("png", "pdf"),
)
print(
"[saved] "
+ ", ".join(str(path.relative_to(ROOT)) for path in saved.values())
)
if __name__ == "__main__":
main()