"""Single-particle level diagrams: Skyrme-HF (SIII) vs experimental single-particle
energies from Schwierz, Wiedenhover & Volya, arXiv:0709.3525 (Table II),
themselves extracted from NNDC/ENSDF transfer- and pickup-reaction data.
"""
import numpy as np
import matplotlib.pyplot as plt
import skyrme_hf as m
from skyrme_hf import Orbital as O
plt.rcParams.update({
"figure.dpi": 140, "savefig.dpi": 140, "font.size": 9.5,
"axes.spines.top": False, "axes.spines.right": False,
"font.family": "DejaVu Sans",
})
CALC_C, EXP_C = "#1a1a1a", "#c0392b"
# ---- experimental single-particle energies (MeV), SWV Table II ---------------
# 16O 1p3/2 holes added from the first 3/2- excitation in 15N/15O (ENSDF).
EXP = {
"16O": {
"p": {"1p3/2": -18.45, "1p1/2": -12.13, "1d5/2": -0.60, "2s1/2": -0.11},
"n": {"1p3/2": -21.84, "1p1/2": -15.66, "1d5/2": -4.14, "2s1/2": -3.27},
},
"40Ca": {
"p": {"1d5/2": -15.07, "2s1/2": -10.92, "1d3/2": -8.33,
"1f7/2": -1.09, "2p3/2": 0.69, "2p1/2": 2.38, "1f5/2": 4.96},
"n": {"1d5/2": -22.39, "2s1/2": -18.19, "1d3/2": -15.64,
"1f7/2": -8.36, "2p3/2": -5.84, "2p1/2": -4.20, "1f5/2": -1.56},
},
"208Pb": {
"p": {"1g7/2": -11.49, "2d5/2": -9.70, "1h11/2": -9.36, "2d3/2": -8.36,
"3s1/2": -8.01, "1h9/2": -3.80, "2f7/2": -2.90, "1i13/2": -2.10,
"2f5/2": -0.97, "3p3/2": -0.68, "3p1/2": -0.16},
"n": {"1h9/2": -11.40, "2f7/2": -9.81, "1i13/2": -9.24, "3p3/2": -8.26,
"2f5/2": -7.94, "3p1/2": -7.37, "2g9/2": -3.94, "1i11/2": -3.16,
"1j15/2": -2.51, "3d5/2": -2.37, "4s1/2": -1.90, "2g7/2": -1.44,
"3d3/2": -1.40},
},
}
# orbital (n,l,2j) lookup for every label we need
def parse(label):
spec = "spdfghijklm"
n = int(label[0]); l = spec.index(label[1]); j2 = int(label.split("/")[0][2:])
return O(n, l, j2)
NUCF = {"16O": m.closed_shell_16O, "40Ca": m.closed_shell_40Ca,
"208Pb": m.closed_shell_208Pb}
# where the Fermi level sits (last occupied), for shading
FERMI = {"16O": {"p": -12.13, "n": -15.66}, "40Ca": {"p": -8.33, "n": -15.64},
"208Pb": {"p": -8.01, "n": -7.37}}
def calc_levels(name):
res = m.run_hf(NUCF[name](), force=m.SIII) if False else \
m.run_hf(m.SIII, NUCF[name](), grid=m.Grid(rmax=22, h=0.1),
mix=0.25, maxiter=400, com_correction=True, verbose=False)
out = {}
for iso in ("p", "n"):
wanted = [parse(lb) for lb in EXP[name][iso]]
lp, ln = m.sp_levels(res, wanted, wanted)
levs = lp if iso == "p" else ln
out[iso] = {o.label(): e for o, e in levs}
return out, res
def _declutter(items, min_gap):
"""items: list of (energy, label). Return label-y positions nudged apart."""
items = sorted(items) # by energy ascending
ys = [e for e, _ in items]
# push upward pass
for i in range(1, len(ys)):
if ys[i] - ys[i-1] < min_gap:
ys[i] = ys[i-1] + min_gap
return {lab: (e, y) for (e, lab), y in zip(items, ys)}
def panel(ax, name, iso, calc, title):
exp = EXP[name][iso]
ef = FERMI[name][iso]
lo = min(min(exp.values()), min(calc.values())) - 1
hi = max(max(exp.values()), max(calc.values())) + 1
ax.axhspan(lo, ef, color="#4a90d9", alpha=0.05)
ax.axhline(ef, color="#4a90d9", lw=0.8, ls=":", alpha=0.7)
for lab, e in exp.items():
ax.hlines(e, 0.55, 0.95, color=EXP_C, lw=1.6)
if lab in calc:
ax.plot([0.45, 0.55], [calc[lab], e], color="#ccc", lw=0.6, zorder=0)
gap = (hi - lo) * 0.033
lpos = _declutter([(e, lab) for lab, e in calc.items()], gap)
for lab, e in calc.items():
ax.hlines(e, 0.05, 0.45, color=CALC_C, lw=1.6)
ytext = lpos[lab][1]
ax.text(0.02, ytext, lab, ha="right", va="center", fontsize=6.6, color=CALC_C)
if abs(ytext - e) > 1e-6:
ax.plot([0.02, 0.05], [ytext, e], color="#bbb", lw=0.4)
ax.set_xlim(-0.4, 1.02); ax.set_xticks([0.25, 0.75])
ax.set_xticklabels(["SIII", "exp"], fontsize=8)
ax.tick_params(axis="x", length=0)
ax.set_ylim(lo, hi)
ax.set_title(title, fontsize=9.5)
ax.spines["bottom"].set_visible(False)
shared = [lab for lab in calc if lab in exp]
rms = np.sqrt(np.mean([(calc[lab] - exp[lab])**2 for lab in shared]))
ax.text(0.5, 1.005, f"rms {rms:.2f} MeV", transform=ax.transAxes,
ha="center", va="bottom", fontsize=7.5, color="#666")
return rms
# =============================== 16O + 40Ca figure
fig, axes = plt.subplots(1, 4, figsize=(12, 5.2))
data = {}
for name in ("16O", "40Ca"):
data[name] = calc_levels(name)
rmslist = []
rmslist.append(panel(axes[0], "16O", "p", data["16O"][0]["p"], r"$^{16}$O protons"))
rmslist.append(panel(axes[1], "16O", "n", data["16O"][0]["n"], r"$^{16}$O neutrons"))
rmslist.append(panel(axes[2], "40Ca", "p", data["40Ca"][0]["p"], r"$^{40}$Ca protons"))
rmslist.append(panel(axes[3], "40Ca", "n", data["40Ca"][0]["n"], r"$^{40}$Ca neutrons"))
axes[0].set_ylabel("single-particle energy (MeV)")
fig.suptitle("Single-particle spectra: Skyrme-HF (SIII) vs experiment "
"[Schwierz-Wiedenhover-Volya, arXiv:0709.3525]",
y=1.00, fontsize=11, fontweight="bold")
fig.text(0.5, -0.02, "shaded band = occupied (below Fermi level); grey lines connect "
"matched orbitals", ha="center", fontsize=8, color="#666")
fig.tight_layout()
fig.savefig("./fig5_levels_16O_40Ca.png", bbox_inches="tight")
plt.close(fig)
# =============================== 208Pb figure
fig, axes = plt.subplots(1, 2, figsize=(8.5, 7.2))
d208, res208 = calc_levels("208Pb")
r1 = panel(axes[0], "208Pb", "p", d208["p"], r"$^{208}$Pb protons")
r2 = panel(axes[1], "208Pb", "n", d208["n"], r"$^{208}$Pb neutrons")
axes[0].set_ylabel("single-particle energy (MeV)")
fig.suptitle("$^{208}$Pb single-particle spectrum: SIII vs experiment "
"[SWV, arXiv:0709.3525]", y=0.98, fontsize=11, fontweight="bold")
fig.tight_layout()
fig.savefig("./fig6_levels_208Pb.png", bbox_inches="tight")
plt.close(fig)
print("16O/40Ca rms (p,n,p,n):", [f"{x:.2f}" for x in rmslist])
print(f"208Pb rms proton {r1:.2f} neutron {r2:.2f} MeV")
print("wrote fig5_levels_16O_40Ca.png, fig6_levels_208Pb.png")