Oxygen-16 from a Skyrme energy functional

The ground state

import numpy as np
import vipyr as vp
from vipyr.figure import figure

r = vp.nuclear.skyrme_hf()
print("converged in %d iterations: %s" % (r.iterations, r.converged))
print("E / A = %.4f MeV   E = %.3f MeV" % (r.e_per_a, r.energy))
print(
    "rms radii: protons %.4f fm, neutrons %.4f fm, charge %.4f fm"
    % (r.r_p, r.r_n, r.r_ch)
)
# r.labels names each entry of r.levels: "n 1s1/2", ..., "p 1p1/2"
labels = r.labels.split(",")
levels = [
    [e for l, e in zip(labels, r.levels) if l.startswith(k)]
    for k in ("n ", "p ")
]
names = [l[2:] for l in labels if l.startswith("n ")]
for q, who in enumerate(("neutron", "proton")):
    print(
        "%-8s levels: " % who
        + "   ".join("%s %.3f" % (n, e) for n, e in zip(names, levels[q]))
    )

converged in 25 iterations: True
E / A = -8.0465 MeV   E = -128.745 MeV
rms radii: protons 2.6810 fm, neutrons 2.6559 fm, charge 2.7967 fm
neutron  levels: 1s1/2 -36.187   1p3/2 -20.631   1p1/2 -14.424
proton   levels: 1s1/2 -32.390   1p3/2 -17.153   1p1/2 -11.074

Densities

mesh = np.asarray(r.mesh)
fig = figure("Oxygen-16 densities, SLy5", "$r$ (fm)", "fm$^{-3}$")
fig.plot(mesh[:70], r.rho[0][:70], "neutrons")
fig.plot(mesh[:70], r.rho[1][:70], "protons")
fig.plot(mesh[:70], r.tau[1][:70], "protons: $\\tau$")
fig.style("neutrons", "s1 nomarker")
fig.style("protons", "s2 dash nomarker")
fig.style("protons: $\\tau$", "s3 dot nomarker")
print(
    "densities at the center: neutrons %.4f, protons %.4f fm^-3"
    % (r.rho[0][0], r.rho[1][0])
)
print(
    "neutrons counted on the mesh: %.6f"
    % (4 * np.pi * np.sum(mesh**2 * r.rho[0]) * (mesh[1] - mesh[0]))
)

densities at the center: neutrons 0.0769, protons 0.0753 fm^-3
neutrons counted on the mesh: 8.000000

Single-particle levels

fig = figure(
    "Single-particle levels, SLy5",
    "",
    "e (MeV)",
    window=(0.0, 3.0, -40.0, 0.0),
)
for q, who in enumerate(("neutrons", "protons")):
    x0 = 0.3 + 1.4 * q
    xs, ys = [], []
    for e in levels[q]:
        xs += [x0, x0 + 1.0, float("nan")]
        ys += [e, e, float("nan")]
    fig.plot(xs, ys, who)
    fig.style(who, "s%d thick nomarker" % (q + 1))
    for n, e in zip(names, levels[q]):
        fig.label("%s %s" % (who, n), n, x0 + 0.3, e + 1.2)
fig.xtick("neutrons", 0.8)
fig.xtick("protons", 2.2)

Mesh convergence

spacings = [0.2, 0.15, 0.1, 0.075, 0.05]
conv = [vp.nuclear.skyrme_hf(spacing=h) for h in spacings]
e_h = {h: c.e_per_a for h, c in zip(spacings, conv)}
extrap = e_h[0.05] + (e_h[0.05] - e_h[0.1]) / 3.0
for h, c in zip(spacings, conv):
    print(
        "h = %.3f fm   E/A = %.6f MeV   r_ch = %.4f fm   |E/A - E/A(0)| = "
        "%.1e" % (h, c.e_per_a, c.r_ch, abs(c.e_per_a - extrap))
    )
halving = np.log2((e_h[0.2] - e_h[0.1]) / (e_h[0.1] - e_h[0.05]))
slope = np.polyfit(
    np.log(spacings),
    np.log([abs(c.e_per_a - extrap) for c in conv]),
    1,
)[0]
print("order between halvings 0.2, 0.1, 0.05 fm: %.3f" % halving)
print(
    "fitted order against the limit: %.3f   E/A extrapolated from 0.1 and 0.05 fm: %.4f MeV"
    % (slope, extrap)
)
assert 1.9 < halving < 2.1, "the mesh error is no longer second order: %.3f" % halving
fig = figure(
    "$E/A$ against mesh spacing", "$h$ (fm)", "$|E/A(h) - E/A(0)|$ (MeV)"
)
fig.plot(spacings, [abs(c.e_per_a - extrap) for c in conv], "SLy5")
fig.xscale("log")
fig.yscale("log")

h = 0.200 fm   E/A = -8.110660 MeV   r_ch = 2.7890 fm   |E/A - E/A(0)| = 8.6e-02
h = 0.150 fm   E/A = -8.073391 MeV   r_ch = 2.7935 fm   |E/A - E/A(0)| = 4.8e-02
h = 0.100 fm   E/A = -8.046539 MeV   r_ch = 2.7967 fm   |E/A - E/A(0)| = 2.2e-02
h = 0.075 fm   E/A = -8.037108 MeV   r_ch = 2.7978 fm   |E/A - E/A(0)| = 1.2e-02
h = 0.050 fm   E/A = -8.030366 MeV   r_ch = 2.7986 fm   |E/A - E/A(0)| = 5.4e-03
order between halvings 0.2, 0.1, 0.05 fm: 1.987
fitted order against the limit: 1.996   E/A extrapolated from 0.1 and 0.05 fm: -8.0250 MeV

Against experiment

ba_exp, rch_exp = 7.976, 2.699
print(
    "B/A: SLy5 %.3f MeV (h -> 0), experiment %.3f MeV (%+.1f%%)"
    % (-extrap, ba_exp, 100 * (-extrap / ba_exp - 1))
)
print(
    "r_ch: SLy5 %.3f fm (h = 0.05 fm), experiment %.3f fm (%+.1f%%)"
    % (conv[-1].r_ch, rch_exp, 100 * (conv[-1].r_ch / rch_exp - 1))
)

B/A: SLy5 8.025 MeV (h -> 0), experiment 7.976 MeV (+0.6%)
r_ch: SLy5 2.799 fm (h = 0.05 fm), experiment 2.699 fm (+3.7%)

How it's checked

print("check():", r.check().passed)

check(): True

Source

gallery/science/plot_oxygen16_skyrme.py, run as written when these pages were built. The output and the figures above are from that run.

The Viv program it calls is posts/2026-09-25-skyrme-hartree-fock-of-oxygen-16.viv: its atlas, the defs and the calls between them, or open it in the workbench.