"""Figures for the Skyrme-HF / Vautherin-Brink reproduction."""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import gridspec
import skyrme_hf as m
plt.rcParams.update({
"figure.dpi": 140, "savefig.dpi": 140,
"font.size": 10, "axes.titlesize": 11, "axes.labelsize": 10,
"axes.spines.top": False, "axes.spines.right": False,
"axes.grid": True, "grid.alpha": 0.25, "grid.linewidth": 0.6,
"legend.frameon": False, "font.family": "DejaVu Sans",
})
INK = "#1a1a1a"
COL = {"SIII": "#c0392b", "SI": "#2c7fb8", "SII": "#31a354", "exp": INK}
EXPT = {"16O": (7.976, 2.699), "40Ca": (8.551, 3.478), "48Ca": (8.666, 3.477),
"90Zr": (8.710, 4.269), "208Pb": (7.867, 5.501)}
ORDER = ["16O", "40Ca", "48Ca", "90Zr", "208Pb"]
AOF = {"16O": 16, "40Ca": 40, "48Ca": 48, "90Zr": 90, "208Pb": 208}
NUCF = {"16O": m.closed_shell_16O, "40Ca": m.closed_shell_40Ca,
"48Ca": m.closed_shell_48Ca, "90Zr": m.closed_shell_90Zr,
"208Pb": m.closed_shell_208Pb}
def solve(force, name, h=0.1, rmax=22.0):
g = m.Grid(rmax=rmax, h=h)
res = m.run_hf(force, NUCF[name](), grid=g, mix=0.25, maxiter=400,
com_correction=True, verbose=False)
res["grid"] = g
return res
# ---------------------------------------------------------------- gather data
data = {f.name: {nm: solve(f, nm) for nm in ORDER} for f in (m.SIII, m.SI, m.SII)}
# ============================================================ FIG 1: validation
fig = plt.figure(figsize=(11, 4.4))
gs = gridspec.GridSpec(1, 2, wspace=0.28)
A = [AOF[n] for n in ORDER]
ax = fig.add_subplot(gs[0])
ax.plot(A, [EXPT[n][0] for n in ORDER], "o-", color=INK, lw=1.6, ms=7,
label="experiment", zorder=5)
for fn in ("SIII", "SI", "SII"):
ax.plot(A, [data[fn][n]["BE_per_A"] for n in ORDER], "s--",
color=COL[fn], lw=1.2, ms=5, alpha=0.9, label=fn)
for n in ORDER:
ax.annotate(n, (AOF[n], EXPT[n][0]), textcoords="offset points",
xytext=(0, 9), ha="center", fontsize=8, color="#555")
ax.set_xlabel("mass number $A$"); ax.set_ylabel("binding energy per nucleon $B/A$ (MeV)")
ax.set_title("Binding energy across the doubly-magic chain")
ax.legend(loc="lower right")
ax2 = fig.add_subplot(gs[1])
ax2.plot(A, [EXPT[n][1] for n in ORDER], "o-", color=INK, lw=1.6, ms=7,
label="experiment", zorder=5)
for fn in ("SIII", "SI", "SII"):
ax2.plot(A, [data[fn][n]["rms_ch"] for n in ORDER], "s--",
color=COL[fn], lw=1.2, ms=5, alpha=0.9, label=fn)
ax2.set_xlabel("mass number $A$"); ax2.set_ylabel(r"rms charge radius $r_{\rm ch}$ (fm)")
ax2.set_title("Charge radii")
ax2.legend(loc="lower right")
fig.suptitle("Skyrme-Hartree-Fock reproduction of Vautherin & Brink (1972)",
y=1.02, fontsize=12, fontweight="bold")
fig.savefig("./fig1_validation.png", bbox_inches="tight")
plt.close(fig)
# residual panel: BE/A deviation
fig = plt.figure(figsize=(7.5, 4.2))
ax = fig.add_subplot(111)
ax.axhline(0, color=INK, lw=1)
for fn in ("SIII", "SI", "SII"):
ax.plot(A, [data[fn][n]["BE_per_A"] - EXPT[n][0] for n in ORDER], "o-",
color=COL[fn], lw=1.3, ms=6, label=fn)
for n in ORDER:
ax.annotate(n, (AOF[n], 0), textcoords="offset points", xytext=(0, -14),
ha="center", fontsize=8, color="#555")
ax.set_xlabel("mass number $A$")
ax.set_ylabel(r"$(B/A)_{\rm calc}-(B/A)_{\rm exp}$ (MeV)")
ax.set_title("Binding-energy residual vs experiment")
ax.legend()
ax.set_ylim(-0.55, 0.55)
fig.savefig("./fig2_residual.png", bbox_inches="tight")
plt.close(fig)
# ============================================================ FIG 3: densities
fig = plt.figure(figsize=(11, 4.4))
gs = gridspec.GridSpec(1, 2, wspace=0.24)
for k, nm in enumerate(["16O", "208Pb"]):
ax = fig.add_subplot(gs[k])
res = data["SIII"][nm]; g = res["grid"]; r = g.r
rho_n, rho_p = res["dens"]
ax.plot(r, rho_p, color="#c0392b", lw=1.8, label=r"proton $\rho_p$")
ax.plot(r, rho_n, color="#2c7fb8", lw=1.8, label=r"neutron $\rho_n$")
ax.plot(r, rho_n + rho_p, color=INK, lw=1.2, ls="--", label=r"total $\rho$")
ax.set_xlim(0, 10 if nm == "16O" else 12)
ax.set_xlabel("$r$ (fm)"); ax.set_ylabel(r"density (fm$^{-3}$)")
ax.set_title(f"{nm} density profile (SIII)")
ax.axhline(0.16, color="#888", lw=0.8, ls=":")
ax.text(0.98*ax.get_xlim()[1], 0.163, r"$\rho_0$", ha="right", va="bottom",
fontsize=8, color="#888")
ax.legend(loc="upper right")
fig.suptitle("Self-consistent nucleon densities (note surface diffuseness & interior shell structure)",
y=1.02, fontsize=11, fontweight="bold")
fig.savefig("./fig3_densities.png", bbox_inches="tight")
plt.close(fig)
# ============================================================ FIG 4: 208Pb levels
# experimental single-particle energies (MeV), occupied/hole states near surface
EXP_P = {"3s1/2": -8.01, "2d3/2": -8.36, "1h11/2": -9.36, "2d5/2": -9.70,
"1g7/2": -11.49, "2p1/2": -14.9}
EXP_N = {"3p1/2": -7.37, "2f5/2": -7.94, "3p3/2": -8.27, "1i13/2": -9.00,
"2f7/2": -9.71, "1h9/2": -10.78, "3s1/2": -11.4}
res = data["SIII"]["208Pb"]
lp = {o.label(): e for o, e in res["levels_p"]}
ln = {o.label(): e for o, e in res["levels_n"]}
fig = plt.figure(figsize=(9, 6.2))
ax = fig.add_subplot(111)
def draw(col_calc, col_exp, calc, exp, xc):
# calc levels (left tick), exp (right tick), connect known partners
for lab, e in calc.items():
ax.hlines(e, xc - 0.34, xc - 0.02, color=INK, lw=1.4)
ax.text(xc - 0.37, e, lab, ha="right", va="center", fontsize=7.5, color=INK)
for lab, e in exp.items():
ax.hlines(e, xc + 0.02, xc + 0.34, color="#c0392b", lw=1.4)
if lab in calc:
ax.plot([xc - 0.02, xc + 0.02], [calc[lab], e], color="#bbb", lw=0.6, zorder=0)
draw("k", "r", lp, EXP_P, 0.0)
draw("k", "r", ln, EXP_N, 1.6)
ax.set_xticks([-0.18, 0.18, 1.42, 1.78])
ax.set_xticklabels(["calc", "exp", "calc", "exp"], fontsize=9)
ax.text(0.0, 2.0, "PROTONS", ha="center", fontweight="bold", fontsize=10)
ax.text(1.6, 2.0, "NEUTRONS", ha="center", fontweight="bold", fontsize=10)
ax.set_ylabel("single-particle energy (MeV)")
ax.set_title(r"$^{208}$Pb single-particle spectrum near the Fermi surface (SIII vs experiment)")
ax.set_ylim(-16, 3.0)
ax.set_xlim(-0.7, 2.3)
ax.spines["bottom"].set_visible(False)
ax.tick_params(axis="x", length=0)
ax.grid(axis="x")
fig.savefig("./fig4_pb208_levels.png", bbox_inches="tight")
plt.close(fig)
print("wrote fig1_validation.png, fig2_residual.png, fig3_densities.png, fig4_pb208_levels.png")