Notebooks · posts/skyrme-hartree-fock-of-oxygen-16

Skyrme Hartree-Fock of oxygen-16

Adam Zettel · September 25, 2026 · written with an agent

The ground state of $^{16}$ O from the SLy5 Skyrme functional, solved on a radial mesh at two spacings and extrapolated. The binding energy comes out 8.0250 MeV per nucleon and the charge radius 2.7993 fm; skyrme_rpa, the published code of Colò et al. (2013), gives 8.0261 MeV and 2.7991 fm in its own test run with the same force, and the six single-particle energies agree with it to 3 keV.

E = \int H [ρ_{q}, τ_{q}, J_{q}] (r) d^{3} r, [- \nabla \cdot \frac{ℏ ^{2}}{2 m _{q}^{*} ( r )} \nabla + U_{q} + W_{q} \cdot (- i) (\nabla \times σ)] ϕ_{i} = ε_{i} ϕ_{i}

A nucleus is up to a few hundred protons and neutrons (sixteen here) held together by the strong force, which has no simple closed form. Nuclear density functional theory replaces it with an energy that depends only on a few local densities: the number density $ρ_{q}$ of neutrons ( $q = n$ ) and protons ( $q = p$ ), their kinetic-energy densities $τ_{q}$ , and their spin-orbit densities $J_{q}$ . Skyrme's functional is the oldest of these, with about ten parameters fitted to nuclear data. Minimizing the energy gives each nucleon a Schrödinger equation (the second equation above) with a position-dependent effective mass $m_{q}^{*}$ , a central potential $U_{q}$ and a spin-orbit field $W_{q}$ , all computed from the densities. That makes the problem self-consistent, and the program solves it by iteration: densities, then fields, then orbitals, then new densities.

$^{16}$ O has 8 protons and 8 neutrons, which fill the $1 s_{1/2}$ , $1 p_{3/2}$ and $1 p_{1/2}$ shells of each, so it is spherical and every orbital is the lowest of its $(l, j)$ block. The force is SLy5 (Chabanat et al., Nucl. Phys. A 643, 441 (1998)), with the parameters and conventions of skyrme_rpa (G. Colò, L. Cao, N. Van Giai and L. Capelli, Comput. Phys. Commun. 184, 142 (2013)), whose distribution includes a Hartree-Fock run of $^{16}$ O with SLy5 as its test case. That output file is the published reference here. posts/refs/skyrme-hartree-fock-of-oxygen-16.py copies its numbers, describes rerunning the program, and holds a numpy twin of this solver that diagonalizes each block densely.

Look inside this program

Atlas
The program’s map: every function and which calls which.
Explore
Run it and step through what each part did, with the code, values and plots beside it.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 37–75

# ---- SLy5, as skyrme_rpa's skyrme.in lists it (MeV and fm)
T0 = -2484.88
T1 = 483.13
T2 = -549.40
T3 = 13763.0
X0 = 0.778
X1 = -0.328
X2 = -1.0
X3 = 1.267
W0 = 126.0
ALPHA = 1.0 / 6.0

# ---- skyrme_rpa's constants: hbar c, the nucleon mass (the p-n average), e^2
HBARC = 197.327053
M_NUC = 0.5 * (938.27231 + 939.56563)
E2 = 1.43986
PI = 3.141592653589793

# ---- oxygen-16
Z_P = 8
N_N = 8
A_NUC = Z_P + N_N
# the center-of-mass correction: hbar^2/2m times (1 - 1/A)
HB = HBARC * HBARC / (2.0 * M_NUC) * (1.0 - 1.0 / A_NUC)
# the occupied shells of each kind: l, j and the number of nucleons 2j + 1
SHELL_L = [0, 1, 1]
SHELL_J = [0.5, 1.5, 0.5]
# the Gaussian range of the proton charge in skyrme_rpa's charge density (fm)
FOLD = 0.65

# ---- the iteration
MIX = 0.5
E_TOL = 1.0e-9
MAX_ITER = 300
# in the first iteration bisection brackets each orbital energy to EIG_TOL (MeV); after it the
# inverse-iteration shift is NEAR below the orbital's last energy. INV_SWEEPS solves per iteration.
EIG_TOL = 1.0e-11
NEAR = 0.1
INV_SWEEPS = 3

# ---- SLy5, as skyrme_rpa's skyrme.in lists it (MeV and fm)
T0 = -2484.88
T1 = 483.13
T2 = -549.40
T3 = 13763.0
X0 = 0.778
X1 = -0.328
X2 = -1.0
X3 = 1.267
W0 = 126.0
ALPHA = 1.0 / 6.0

# ---- skyrme_rpa's constants: hbar c, the nucleon mass (the p-n average), e^2
HBARC = 197.327053
M_NUC = 0.5 * (938.27231 + 939.56563)
E2 = 1.43986
PI = 3.141592653589793

# ---- oxygen-16
Z_P = 8
N_N = 8
A_NUC = Z_P + N_N
# the center-of-mass correction: hbar^2/2m times (1 - 1/A)
HB = HBARC * HBARC / (2.0 * M_NUC) * (1.0 - 1.0 / A_NUC)
# the occupied shells of each kind: l, j and the number of nucleons 2j + 1
SHELL_L = [0, 1, 1]
SHELL_J = [0.5, 1.5, 0.5]
# the Gaussian range of the proton charge in skyrme_rpa's charge density (fm)
FOLD = 0.65

# ---- the iteration
MIX = 0.5
E_TOL = 1.0e-9
MAX_ITER = 300
# in the first iteration bisection brackets each orbital energy to EIG_TOL (MeV); after it the
# inverse-iteration shift is NEAR below the orbital's last energy. INV_SWEEPS solves per iteration.
EIG_TOL = 1.0e-11
NEAR = 0.1
INV_SWEEPS = 3

The discretization

Each orbital is $ϕ = \frac{u ( r )}{r} [Y_{l} \otimes χ]_{j}$ , and on a mesh $r_{i} = ih$ , $i = 1 \dots N$ , the radial equation for $u$ is a symmetric tridiagonal matrix:

- \frac{d}{d r} B \frac{d u}{d r} + [B \frac{l ( l + 1 )}{r ^{2}} + \frac{B ^{'}}{r} + U + V_{C} + \frac{⟨ σ \cdot l ⟩ W}{r}] u = ε u, B = \frac{ℏ ^{2}}{2 m ^{*}} .

The first term uses $B$ at the half points, which keeps the matrix symmetric; the $B^{'} / r$ term comes from writing the equation for $u = r ϕ$ rather than $ϕ$ . Each block's lowest eigenvector comes from inverse iteration: three solves of $(H - σ) w = u$ by the tridiagonal (Thomas) algorithm per pass, starting from the previous pass's vector, with the shift $σ$ 0.1 MeV below the previous orbital energy. On the first pass the eigenvalue itself is found first, by bisection on the Sturm count (the number of negative pivots of $H - x$ is the number of eigenvalues below $x$ ), so the iteration starts on the right state. Derivatives of the densities are centered differences with second-order one-sided ones at the ends, and integrals are trapezoidal. The protons also feel the Coulomb field of their own density, with exchange in the Slater approximation.

Every piece is second order in $h$ , so the program solves the problem at $h = 0.1$ and $h = 0.05$ fm in a 15 fm box and extrapolates, $x \approx (4 x_{h /2} - x_{h}) /3$ .

[2]lines 104–440

#| First derivative on the mesh: centered inside, second-order one-sided at
#| both ends (numpy.gradient with edge_order = 2).
def grad1(f: Arr, h, out: Arr) {
    n = len(f)
    out[0] = (-1.5 * f[0] + 2.0 * f[1] - 0.5 * f[2]) / h
    for i from 1 to n - 1 {
        out[i] = (f[i + 1] - f[i - 1]) / (2.0 * h)
    }
    out[n - 1] = (1.5 * f[n - 1] - 2.0 * f[n - 2] + 0.5 * f[n - 3]) / h
    return 0
}

#| The radial Laplacian f'' + 2 f' / r, from two first derivatives.
def lap_r(f: Arr, r: Arr, h, out: Arr) {
    n = len(f)
    d1 = zeros(n)
    d2 = zeros(n)
    grad1(f, h, d1)
    grad1(d1, h, d2)
    for i to n {
        out[i] = d2[i] + 2.0 * d1[i] / r[i]
    }
    return 0
}

#| 4 pi times the integral of f r^2 dr, by the trapezoid rule on the mesh.
def integrate_r2(f: Arr, r: Arr, h) {
    n = len(f)
    s = 0.0
    for i to n {
        s += f[i] * r[i] * r[i]
    }
    s -= 0.5 * (f[0] * r[0] * r[0] + f[n - 1] * r[n - 1] * r[n - 1])
    return 4.0 * PI * h * s
}

#| The direct Coulomb potential of the proton density: the charge inside r
#| over r, plus the shells outside r.
def coulomb(rp: Arr, r: Arr, h, out: Arr) {
    n = len(rp)
    inner = 0.0
    for i to n {
        inner += rp[i] * r[i] * r[i]
        out[i] = inner / r[i]
    }
    outer = 0.0
    for k from 0 to n by 1 reversed {
        out[k] = E2 * 4.0 * PI * h * (out[k] + outer)
        outer += rp[k] * r[k]
    }
    return 0
}

#| The number of eigenvalues of the symmetric tridiagonal matrix (diag, off)
#| below x: the number of negative pivots in the triangular factorization of H - x (Sturm).
def count_below(diag: Arr, off: Arr, x) {
    n = len(diag)
    cnt = 0
    p = diag[0] - x
    if p < 0.0 {
        cnt += 1
    }
    for i from 1 to n {
        if p == 0.0 {
            p = 1.0e-300
        }
        p = diag[i] - x - off[i - 1] * off[i - 1] / p
        if p < 0.0 {
            cnt += 1
        }
    }
    return cnt
}

#| The lowest eigenvalue of the symmetric tridiagonal matrix (diag, off),
#| by bisection on the Sturm count between lo and hi, to EIG_TOL.
def lowest_eigenvalue(diag: Arr, off: Arr, lo, hi) {
    a = lo
    b = hi
    while b - a > EIG_TOL {
        mid = 0.5 * (a + b)
        if count_below(diag, off, mid) > 0 {
            b = mid
        } else {
            a = mid
        }
    }
    return a
}

#| The eigenvector of (diag, off) whose eigenvalue is nearest sigma, by
#| INV_SWEEPS solves of (H - sigma) w = u with the Thomas algorithm, starting
#| from u, which holds the unit vector on return. Returns its Rayleigh quotient.
def inverse_iteration(diag: Arr, off: Arr, sigma, u: Arr) {
    n = len(diag)
    c = zeros(n)
    d = zeros(n)
    w = zeros(n)
    for sweep to INV_SWEEPS {
        c[0] = off[0] / (diag[0] - sigma)
        d[0] = u[0] / (diag[0] - sigma)
        for i from 1 to n {
            m = diag[i] - sigma - off[i - 1] * c[i - 1]
            if i < n - 1 {
                c[i] = off[i] / m
            }
            d[i] = (u[i] - off[i - 1] * d[i - 1]) / m
        }
        w[n - 1] = d[n - 1]
        for k from 0 to n - 1 by 1 reversed {
            w[k] = d[k] - c[k] * w[k + 1]
        }
        nrm = sqrt(sum(w * w))
        if sum(w * u) < 0.0 {
            nrm = -nrm
        }
        for i to n {
            u[i] = w[i] / nrm
        }
    }
    lam = 0.0
    for i to n {
        lam += diag[i] * u[i] * u[i]
        if i < n - 1 {
            lam += 2.0 * off[i] * u[i] * u[i + 1]
        }
    }
    return lam
}

#| Hartree-Fock of oxygen-16 with SLy5 on nr points of spacing h. Writes the
#| converged neutron and proton densities into dens[2, nr] and returns
#| [E/A, r_p, r_n, r_ch, the neutron 1s1/2 1p3/2 1p1/2 energies, the proton
#| ones, iterations].
def skyrme_o16(h, nr, dens: Arr) {
    r = zeros(nr)
    for i to nr {
        r[i] = (i + 1) * h
    }
    nsh = len(SHELL_L)
    # densities per kind (0 = neutrons, 1 = protons): rho, tau, J
    rho = zeros(2, nr)
    tau = zeros(2, nr)
    jso = zeros(2, nr)
    prof = zeros(nr)
    for i to nr {
        prof[i] = 1.0 / (1.0 + exp((r[i] - 1.2 * pow(A_NUC, 1.0 / 3.0)) / 0.6))
    }
    pn = integrate_r2(prof, r, h)
    for i to nr {
        rho[0, i] = prof[i] * N_N / pn
        rho[1, i] = prof[i] * Z_P / pn
    }
    # orbitals: warm-start vectors and last energies, block k = q * nsh + s
    uvec = zeros(2 * nsh, nr)
    lam = zeros(2 * nsh)
    for k to 2 * nsh {
        for i to nr {
            uvec[k, i] = 1.0
        }
    }
    rt = zeros(nr)
    tt = zeros(nr)
    jt = zeros(nr)
    lap_t = zeros(nr)
    lap_q = zeros(nr)
    drt = zeros(nr)
    drq = zeros(nr)
    sj = zeros(nr)
    dsj = zeros(nr)
    B = zeros(nr)
    dB = zeros(nr)
    U = zeros(nr)
    W = zeros(nr)
    V = zeros(nr)
    diag = zeros(nr)
    off = zeros(nr - 1)
    u = zeros(nr)
    R = zeros(nr)
    dR = zeros(nr)
    nrho = zeros(2, nr)
    ntau = zeros(2, nr)
    njso = zeros(2, nr)
    rq = zeros(nr)
    e_old = 0.0
    E = 0.0
    it = 0
    done = 0
    while done < 1 and it < MAX_ITER {
        it += 1
        for i to nr {
            rt[i] = rho[0, i] + rho[1, i]
            tt[i] = tau[0, i] + tau[1, i]
            jt[i] = jso[0, i] + jso[1, i]
        }
        lap_r(rt, r, h, lap_t)
        grad1(rt, h, drt)
        fill(nrho, 0.0)
        fill(ntau, 0.0)
        fill(njso, 0.0)
        for q to 2 {
            for i to nr {
                rq[i] = rho[q, i]
            }
            lap_r(rq, r, h, lap_q)
            grad1(rq, h, drq)
            for i to nr {
                sj[i] = jt[i] + jso[q, i]
            }
            grad1(sj, h, dsj)
            if q == 1 {
                coulomb(rq, r, h, V)
                for i to nr {
                    V[i] -= E2 * pow(3.0 / PI, 1.0 / 3.0) * pow(rq[i], 1.0 / 3.0)
                }
            }
            for i to nr {
                ra = pow(rt[i], ALPHA)
                B[i] = HB + 0.125 * (T1 * (2.0 + X1) + T2 * (2.0 + X2)) * rt[i] + 0.125 * (T2 * (1.0 + 2.0 * X2) - T1 * (1.0 + 2.0 * X1)) * rq[i]
                uc = T0 * ((1.0 + X0 / 2.0) * rt[i] - (X0 + 0.5) * rq[i])
                uc += 0.25 * T1 * ((1.0 + X1 / 2.0) * (tt[i] - 1.5 * lap_t[i]) - (X1 + 0.5) * (tau[q, i] - 1.5 * lap_q[i]))
                uc += 0.25 * T2 * ((1.0 + X2 / 2.0) * (tt[i] + 0.5 * lap_t[i]) + (X2 + 0.5) * (tau[q, i] + 0.5 * lap_q[i]))
                rsq = rho[0, i] * rho[0, i] + rho[1, i] * rho[1, i]
                uc += T3 / 12.0 * ((1.0 + X3 / 2.0) * (2.0 + ALPHA) * ra * rt[i] - (X3 + 0.5) * (2.0 * ra * rq[i] + ALPHA * ra / rt[i] * rsq))
                uc -= W0 * (sj[i] / r[i] + 0.5 * dsj[i])
                U[i] = uc
                W[i] = 0.5 * W0 * (drt[i] + drq[i]) + 0.125 * (T1 - T2) * jso[q, i] - 0.125 * (T1 * X1 + T2 * X2) * jt[i]
            }
            grad1(B, h, dB)
            for s to nsh {
                l = SHELL_L[s]
                j = SHELL_J[s]
                sdl = j * (j + 1.0) - l * (l + 1.0) - 0.75
                for i to nr {
                    bm = B[i]
                    if i > 0 {
                        bm = 0.5 * (B[i - 1] + B[i])
                    }
                    bp = 0.5 * (B[i] + HB)
                    if i < nr - 1 {
                        bp = 0.5 * (B[i] + B[i + 1])
                        off[i] = -bp / (h * h)
                    }
                    diag[i] = (bm + bp) / (h * h) + B[i] * l * (l + 1.0) / (r[i] * r[i]) + dB[i] / r[i] + U[i] + sdl * W[i] / r[i]
                    if q == 1 {
                        diag[i] += V[i]
                    }
                }
                k = q * nsh + s
                for i to nr {
                    u[i] = uvec[k, i]
                }
                sigma = lam[k] - NEAR
                if it == 1 {
                    # Gershgorin: every eigenvalue lies within diag[i] +- (|off[i-1]| + |off[i]|)
                    big = 0.0
                    for i to nr - 1 {
                        big = max(big, abs(off[i]))
                    }
                    sigma = lowest_eigenvalue(diag, off, min(diag) - 2.0 * big, max(diag) + 2.0 * big) - EIG_TOL
                }
                lam[k] = inverse_iteration(diag, off, sigma, u)
                for i to nr {
                    uvec[k, i] = u[i]
                    R[i] = u[i] / (sqrt(h) * r[i])
                }
                grad1(R, h, dR)
                wt = (2.0 * j + 1.0) / (4.0 * PI)
                for i to nr {
                    nrho[q, i] += wt * R[i] * R[i]
                    ntau[q, i] += wt * (dR[i] * dR[i] + l * (l + 1.0) / (r[i] * r[i]) * R[i] * R[i])
                    njso[q, i] += wt * sdl / r[i] * R[i] * R[i]
                }
            }
        }
        rho = (1.0 - MIX) * rho + MIX * nrho
        tau = (1.0 - MIX) * tau + MIX * ntau
        jso = (1.0 - MIX) * jso + MIX * njso
        E = skyrme_energy(rho, tau, jso, r, h)
        if abs(E - e_old) < E_TOL and it > 3 {
            done = 1
        }
        e_old = E
    }
    rn = zeros(nr)
    rp = zeros(nr)
    for i to nr {
        rn[i] = rho[0, i] * r[i] * r[i]
        rp[i] = rho[1, i] * r[i] * r[i]
        dens[0, i] = rho[0, i]
        dens[1, i] = rho[1, i]
    }
    r2p = integrate_r2(rp, r, h) / Z_P
    r2n = integrate_r2(rn, r, h) / N_N
    return [E / A_NUC, sqrt(r2p), sqrt(r2n), sqrt(r2p + 1.5 * FOLD * FOLD), lam[0], lam[1], lam[2], lam[3], lam[4], lam[5], it]
}

#| The total energy: the Skyrme energy density integrated, plus the direct
#| Coulomb energy and the Slater exchange energy of the protons.
def skyrme_energy(rho: Arr, tau: Arr, jso: Arr, r: Arr, h) {
    nr = len(r)
    rn = zeros(nr)
    rp = zeros(nr)
    for i to nr {
        rn[i] = rho[0, i]
        rp[i] = rho[1, i]
    }
    rt = rn + rp
    drt = zeros(nr)
    drn = zeros(nr)
    drp = zeros(nr)
    grad1(rt, h, drt)
    grad1(rn, h, drn)
    grad1(rp, h, drp)
    H = zeros(nr)
    vc = zeros(nr)
    coulomb(rp, r, h, vc)
    ec = zeros(nr)
    for i to nr {
        tn = tau[0, i]
        tp = tau[1, i]
        tt = tn + tp
        jn = jso[0, i]
        jp = jso[1, i]
        jt = jn + jp
        rsq = rn[i] * rn[i] + rp[i] * rp[i]
        e = HB * tt + 0.5 * T0 * ((1.0 + X0 / 2.0) * rt[i] * rt[i] - (X0 + 0.5) * rsq)
        e += 0.25 * T1 * ((1.0 + X1 / 2.0) * (rt[i] * tt + 0.75 * drt[i] * drt[i]) - (X1 + 0.5) * (rn[i] * tn + 0.75 * drn[i] * drn[i] + rp[i] * tp + 0.75 * drp[i] * drp[i]))
        e += 0.25 * T2 * ((1.0 + X2 / 2.0) * (rt[i] * tt - 0.25 * drt[i] * drt[i]) + (X2 + 0.5) * (rn[i] * tn - 0.25 * drn[i] * drn[i] + rp[i] * tp - 0.25 * drp[i] * drp[i]))
        e += T3 / 12.0 * pow(rt[i], ALPHA) * ((1.0 + X3 / 2.0) * rt[i] * rt[i] - (X3 + 0.5) * rsq)
        e += (T1 - T2) / 16.0 * (jn * jn + jp * jp) - (T1 * X1 + T2 * X2) / 16.0 * jt * jt
        e += 0.5 * W0 * (jt * drt[i] + jn * drn[i] + jp * drp[i])
        H[i] = e
        ec[i] = 0.5 * vc[i] * rp[i] - 0.75 * E2 * pow(3.0 / PI, 1.0 / 3.0) * pow(rp[i], 4.0 / 3.0)
    }
    return integrate_r2(H, r, h) + integrate_r2(ec, r, h)
}

#| First derivative on the mesh: centered inside, second-order one-sided at
#| both ends (numpy.gradient with edge_order = 2).
def grad1(f: Arr, h, out: Arr) {
    n = len(f)
    out[0] = (-1.5 * f[0] + 2.0 * f[1] - 0.5 * f[2]) / h
    for i from 1 to n - 1 {
        out[i] = (f[i + 1] - f[i - 1]) / (2.0 * h)
    }
    out[n - 1] = (1.5 * f[n - 1] - 2.0 * f[n - 2] + 0.5 * f[n - 3]) / h
    return 0
}

#| The radial Laplacian f'' + 2 f' / r, from two first derivatives.
def lap_r(f: Arr, r: Arr, h, out: Arr) {
    n = len(f)
    d1 = zeros(n)
    d2 = zeros(n)
    grad1(f, h, d1)
    grad1(d1, h, d2)
    for i to n {
        out[i] = d2[i] + 2.0 * d1[i] / r[i]
    }
    return 0
}

#| 4 pi times the integral of f r^2 dr, by the trapezoid rule on the mesh.
def integrate_r2(f: Arr, r: Arr, h) {
    n = len(f)
    s = 0.0
    for i to n {
        s += f[i] * r[i] * r[i]
    }
    s -= 0.5 * (f[0] * r[0] * r[0] + f[n - 1] * r[n - 1] * r[n - 1])
    return 4.0 * PI * h * s
}

#| The direct Coulomb potential of the proton density: the charge inside r
#| over r, plus the shells outside r.
def coulomb(rp: Arr, r: Arr, h, out: Arr) {
    n = len(rp)
    inner = 0.0
    for i to n {
        inner += rp[i] * r[i] * r[i]
        out[i] = inner / r[i]
    }
    outer = 0.0
    for k from 0 to n by 1 reversed {
        out[k] = E2 * 4.0 * PI * h * (out[k] + outer)
        outer += rp[k] * r[k]
    }
    return 0
}

#| The number of eigenvalues of the symmetric tridiagonal matrix (diag, off)
#| below x: the number of negative pivots in the triangular factorization of H - x (Sturm).
def count_below(diag: Arr, off: Arr, x) {
    n = len(diag)
    cnt = 0
    p = diag[0] - x
    if p < 0.0 {
        cnt += 1
    }
    for i from 1 to n {
        if p == 0.0 {
            p = 1.0e-300
        }
        p = diag[i] - x - off[i - 1] * off[i - 1] / p
        if p < 0.0 {
            cnt += 1
        }
    }
    return cnt
}

#| The lowest eigenvalue of the symmetric tridiagonal matrix (diag, off),
#| by bisection on the Sturm count between lo and hi, to EIG_TOL.
def lowest_eigenvalue(diag: Arr, off: Arr, lo, hi) {
    a = lo
    b = hi
    while b - a > EIG_TOL {
        mid = 0.5 * (a + b)
        if count_below(diag, off, mid) > 0 {
            b = mid
        } else {
            a = mid
        }
    }
    return a
}

#| The eigenvector of (diag, off) whose eigenvalue is nearest sigma, by
#| INV_SWEEPS solves of (H - sigma) w = u with the Thomas algorithm, starting
#| from u, which holds the unit vector on return. Returns its Rayleigh quotient.
def inverse_iteration(diag: Arr, off: Arr, sigma, u: Arr) {
    n = len(diag)
    c = zeros(n)
    d = zeros(n)
    w = zeros(n)
    for sweep to INV_SWEEPS {
        c[0] = off[0] / (diag[0] - sigma)
        d[0] = u[0] / (diag[0] - sigma)
        for i from 1 to n {
            m = diag[i] - sigma - off[i - 1] * c[i - 1]
            if i < n - 1 {
                c[i] = off[i] / m
            }
            d[i] = (u[i] - off[i - 1] * d[i - 1]) / m
        }
        w[n - 1] = d[n - 1]
        for k from 0 to n - 1 by 1 reversed {
            w[k] = d[k] - c[k] * w[k + 1]
        }
        nrm = sqrt(sum(w * w))
        if sum(w * u) < 0.0 {
            nrm = -nrm
        }
        for i to n {
            u[i] = w[i] / nrm
        }
    }
    lam = 0.0
    for i to n {
        lam += diag[i] * u[i] * u[i]
        if i < n - 1 {
            lam += 2.0 * off[i] * u[i] * u[i + 1]
        }
    }
    return lam
}

#| Hartree-Fock of oxygen-16 with SLy5 on nr points of spacing h. Writes the
#| converged neutron and proton densities into dens[2, nr] and returns
#| [E/A, r_p, r_n, r_ch, the neutron 1s1/2 1p3/2 1p1/2 energies, the proton
#| ones, iterations].
def skyrme_o16(h, nr, dens: Arr) {
    r = zeros(nr)
    for i to nr {
        r[i] = (i + 1) * h
    }
    nsh = len(SHELL_L)
    # densities per kind (0 = neutrons, 1 = protons): rho, tau, J
    rho = zeros(2, nr)
    tau = zeros(2, nr)
    jso = zeros(2, nr)
    prof = zeros(nr)
    for i to nr {
        prof[i] = 1.0 / (1.0 + exp((r[i] - 1.2 * pow(A_NUC, 1.0 / 3.0)) / 0.6))
    }
    pn = integrate_r2(prof, r, h)
    for i to nr {
        rho[0, i] = prof[i] * N_N / pn
        rho[1, i] = prof[i] * Z_P / pn
    }
    # orbitals: warm-start vectors and last energies, block k = q * nsh + s
    uvec = zeros(2 * nsh, nr)
    lam = zeros(2 * nsh)
    for k to 2 * nsh {
        for i to nr {
            uvec[k, i] = 1.0
        }
    }
    rt = zeros(nr)
    tt = zeros(nr)
    jt = zeros(nr)
    lap_t = zeros(nr)
    lap_q = zeros(nr)
    drt = zeros(nr)
    drq = zeros(nr)
    sj = zeros(nr)
    dsj = zeros(nr)
    B = zeros(nr)
    dB = zeros(nr)
    U = zeros(nr)
    W = zeros(nr)
    V = zeros(nr)
    diag = zeros(nr)
    off = zeros(nr - 1)
    u = zeros(nr)
    R = zeros(nr)
    dR = zeros(nr)
    nrho = zeros(2, nr)
    ntau = zeros(2, nr)
    njso = zeros(2, nr)
    rq = zeros(nr)
    e_old = 0.0
    E = 0.0
    it = 0
    done = 0
    while done < 1 and it < MAX_ITER {
        it += 1
        for i to nr {
            rt[i] = rho[0, i] + rho[1, i]
            tt[i] = tau[0, i] + tau[1, i]
            jt[i] = jso[0, i] + jso[1, i]
        }
        lap_r(rt, r, h, lap_t)
        grad1(rt, h, drt)
        fill(nrho, 0.0)
        fill(ntau, 0.0)
        fill(njso, 0.0)
        for q to 2 {
            for i to nr {
                rq[i] = rho[q, i]
            }
            lap_r(rq, r, h, lap_q)
            grad1(rq, h, drq)
            for i to nr {
                sj[i] = jt[i] + jso[q, i]
            }
            grad1(sj, h, dsj)
            if q == 1 {
                coulomb(rq, r, h, V)
                for i to nr {
                    V[i] -= E2 * pow(3.0 / PI, 1.0 / 3.0) * pow(rq[i], 1.0 / 3.0)
                }
            }
            for i to nr {
                ra = pow(rt[i], ALPHA)
                B[i] = HB + 0.125 * (T1 * (2.0 + X1) + T2 * (2.0 + X2)) * rt[i] + 0.125 * (T2 * (1.0 + 2.0 * X2) - T1 * (1.0 + 2.0 * X1)) * rq[i]
                uc = T0 * ((1.0 + X0 / 2.0) * rt[i] - (X0 + 0.5) * rq[i])
                uc += 0.25 * T1 * ((1.0 + X1 / 2.0) * (tt[i] - 1.5 * lap_t[i]) - (X1 + 0.5) * (tau[q, i] - 1.5 * lap_q[i]))
                uc += 0.25 * T2 * ((1.0 + X2 / 2.0) * (tt[i] + 0.5 * lap_t[i]) + (X2 + 0.5) * (tau[q, i] + 0.5 * lap_q[i]))
                rsq = rho[0, i] * rho[0, i] + rho[1, i] * rho[1, i]
                uc += T3 / 12.0 * ((1.0 + X3 / 2.0) * (2.0 + ALPHA) * ra * rt[i] - (X3 + 0.5) * (2.0 * ra * rq[i] + ALPHA * ra / rt[i] * rsq))
                uc -= W0 * (sj[i] / r[i] + 0.5 * dsj[i])
                U[i] = uc
                W[i] = 0.5 * W0 * (drt[i] + drq[i]) + 0.125 * (T1 - T2) * jso[q, i] - 0.125 * (T1 * X1 + T2 * X2) * jt[i]
            }
            grad1(B, h, dB)
            for s to nsh {
                l = SHELL_L[s]
                j = SHELL_J[s]
                sdl = j * (j + 1.0) - l * (l + 1.0) - 0.75
                for i to nr {
                    bm = B[i]
                    if i > 0 {
                        bm = 0.5 * (B[i - 1] + B[i])
                    }
                    bp = 0.5 * (B[i] + HB)
                    if i < nr - 1 {
                        bp = 0.5 * (B[i] + B[i + 1])
                        off[i] = -bp / (h * h)
                    }
                    diag[i] = (bm + bp) / (h * h) + B[i] * l * (l + 1.0) / (r[i] * r[i]) + dB[i] / r[i] + U[i] + sdl * W[i] / r[i]
                    if q == 1 {
                        diag[i] += V[i]
                    }
                }
                k = q * nsh + s
                for i to nr {
                    u[i] = uvec[k, i]
                }
                sigma = lam[k] - NEAR
                if it == 1 {
                    # Gershgorin: every eigenvalue lies within diag[i] +- (|off[i-1]| + |off[i]|)
                    big = 0.0
                    for i to nr - 1 {
                        big = max(big, abs(off[i]))
                    }
                    sigma = lowest_eigenvalue(diag, off, min(diag) - 2.0 * big, max(diag) + 2.0 * big) - EIG_TOL
                }
                lam[k] = inverse_iteration(diag, off, sigma, u)
                for i to nr {
                    uvec[k, i] = u[i]
                    R[i] = u[i] / (sqrt(h) * r[i])
                }
                grad1(R, h, dR)
                wt = (2.0 * j + 1.0) / (4.0 * PI)
                for i to nr {
                    nrho[q, i] += wt * R[i] * R[i]
                    ntau[q, i] += wt * (dR[i] * dR[i] + l * (l + 1.0) / (r[i] * r[i]) * R[i] * R[i])
                    njso[q, i] += wt * sdl / r[i] * R[i] * R[i]
                }
            }
        }
        rho = (1.0 - MIX) * rho + MIX * nrho
        tau = (1.0 - MIX) * tau + MIX * ntau
        jso = (1.0 - MIX) * jso + MIX * njso
        E = skyrme_energy(rho, tau, jso, r, h)
        if abs(E - e_old) < E_TOL and it > 3 {
            done = 1
        }
        e_old = E
    }
    rn = zeros(nr)
    rp = zeros(nr)
    for i to nr {
        rn[i] = rho[0, i] * r[i] * r[i]
        rp[i] = rho[1, i] * r[i] * r[i]
        dens[0, i] = rho[0, i]
        dens[1, i] = rho[1, i]
    }
    r2p = integrate_r2(rp, r, h) / Z_P
    r2n = integrate_r2(rn, r, h) / N_N
    return [E / A_NUC, sqrt(r2p), sqrt(r2n), sqrt(r2p + 1.5 * FOLD * FOLD), lam[0], lam[1], lam[2], lam[3], lam[4], lam[5], it]
}

#| The total energy: the Skyrme energy density integrated, plus the direct
#| Coulomb energy and the Slater exchange energy of the protons.
def skyrme_energy(rho: Arr, tau: Arr, jso: Arr, r: Arr, h) {
    nr = len(r)
    rn = zeros(nr)
    rp = zeros(nr)
    for i to nr {
        rn[i] = rho[0, i]
        rp[i] = rho[1, i]
    }
    rt = rn + rp
    drt = zeros(nr)
    drn = zeros(nr)
    drp = zeros(nr)
    grad1(rt, h, drt)
    grad1(rn, h, drn)
    grad1(rp, h, drp)
    H = zeros(nr)
    vc = zeros(nr)
    coulomb(rp, r, h, vc)
    ec = zeros(nr)
    for i to nr {
        tn = tau[0, i]
        tp = tau[1, i]
        tt = tn + tp
        jn = jso[0, i]
        jp = jso[1, i]
        jt = jn + jp
        rsq = rn[i] * rn[i] + rp[i] * rp[i]
        e = HB * tt + 0.5 * T0 * ((1.0 + X0 / 2.0) * rt[i] * rt[i] - (X0 + 0.5) * rsq)
        e += 0.25 * T1 * ((1.0 + X1 / 2.0) * (rt[i] * tt + 0.75 * drt[i] * drt[i]) - (X1 + 0.5) * (rn[i] * tn + 0.75 * drn[i] * drn[i] + rp[i] * tp + 0.75 * drp[i] * drp[i]))
        e += 0.25 * T2 * ((1.0 + X2 / 2.0) * (rt[i] * tt - 0.25 * drt[i] * drt[i]) + (X2 + 0.5) * (rn[i] * tn - 0.25 * drn[i] * drn[i] + rp[i] * tp - 0.25 * drp[i] * drp[i]))
        e += T3 / 12.0 * pow(rt[i], ALPHA) * ((1.0 + X3 / 2.0) * rt[i] * rt[i] - (X3 + 0.5) * rsq)
        e += (T1 - T2) / 16.0 * (jn * jn + jp * jp) - (T1 * X1 + T2 * X2) / 16.0 * jt * jt
        e += 0.5 * W0 * (jt * drt[i] + jn * drn[i] + jp * drp[i])
        H[i] = e
        ec[i] = 0.5 * vc[i] * rp[i] - 0.75 * E2 * pow(3.0 / PI, 1.0 / 3.0) * pow(rp[i], 4.0 / 3.0)
    }
    return integrate_r2(H, r, h) + integrate_r2(ec, r, h)
}

Two meshes and the extrapolation

The iteration starts from a Woods-Saxon density and mixes half of each new density into the old one until the energy changes by less than $1 0^{- 9}$ MeV, which takes 38 passes on each mesh. On both meshes the energy per nucleon, the three radii and the six orbital energies match the numpy twin to $7 \times 1 0^{- 9}$ ; the tests allow $1 0^{- 7}$ . Going from $h = 0.1$ to $h = 0.05$ fm moves the energy per nucleon by 0.016 MeV and the charge radius by 0.002 fm.

[3]lines 454–504

# 150 points of 0.1 fm and 300 of 0.05 fm: a 15 fm box both times
H_COARSE = 0.1
H_FINE = 0.05
BOX = 15.0
nr_c = round(BOX / H_COARSE)
nr_f = round(BOX / H_FINE)
dens_c = zeros(2, nr_c)
dens_f = zeros(2, nr_f)
res_c = skyrme_o16(H_COARSE, nr_c, dens_c)
res_f = skyrme_o16(H_FINE, nr_f, dens_f)
ext = (4.0 * res_f - res_c) / 3.0
show res_c
show res_f
BE_per_A = -ext[0]
r_charge = ext[3]
show BE_per_A
show r_charge

# the numpy twin (posts/refs), same discretization, blocks diagonalized by eigh
TWIN_C = [-8.0465386180, 2.6809867361, 2.6558537301, 2.7966837288, -36.18659979, -20.63141616, -14.42364526, -32.38991872, -17.15296823, -11.07440402]
TWIN_F = [-8.0303664488, 2.6830002205, 2.6578928742, 2.7986139754, -36.14944629, -20.61604351, -14.42671231, -32.35678100, -17.14118483, -11.08056169]
# the result vectors: E/A at 0, the three radii at 1 to 3, the six orbital energies from 4
FIRST_RADIUS = 1
FIRST_LEVEL = 4
# skyrme_rpa's published test run, 16O with SLy5 (sample_skyrme_rpa.out)
PUB = [-8.0261107, 2.6836, 2.6585, 2.7991, -36.138, -20.613, -14.429, -32.346, -17.138, -11.083]

twin_diff = 0.0
for k to len(TWIN_C) {
    twin_diff = max(twin_diff, max(abs(res_c[k] - TWIN_C[k]), abs(res_f[k] - TWIN_F[k])))
}
show twin_diff
# what halving the mesh moves: the energy per nucleon and the charge radius
mesh_shift_e = res_f[0] - res_c[0]
mesh_shift_rc = res_f[3] - res_c[3]
show mesh_shift_e
show mesh_shift_rc

test "on both meshes E/A, the radii and the six orbital energies are the numpy twin's to 1e-7" {
    for k to len(TWIN_C) {
        assert_near(res_c[k], TWIN_C[k], 1.0e-7)
        assert_near(res_f[k], TWIN_F[k], 1.0e-7)
    }
}
test "the largest difference from the numpy twin is under 7e-9" {
    assert_near(twin_diff, 0.0, 7.0e-9)
}
test "going from h = 0.1 to 0.05 fm moves E/A by 0.016 MeV and the charge radius by 0.002 fm" {
    assert_near(mesh_shift_e, TWIN_F[0] - TWIN_C[0], 1.0e-7)
    assert_near(mesh_shift_rc, TWIN_F[3] - TWIN_C[3], 1.0e-7)
}

# 150 points of 0.1 fm and 300 of 0.05 fm: a 15 fm box both times
H_COARSE = 0.1
H_FINE = 0.05
BOX = 15.0
nr_c = round(BOX / H_COARSE)
nr_f = round(BOX / H_FINE)
dens_c = zeros(2, nr_c)
dens_f = zeros(2, nr_f)
res_c = skyrme_o16(H_COARSE, nr_c, dens_c)
res_f = skyrme_o16(H_FINE, nr_f, dens_f)
ext = (4.0 * res_f - res_c) / 3.0
show res_c
show res_f
BE_per_A = -ext[0]
r_charge = ext[3]
show BE_per_A
show r_charge

# the numpy twin (posts/refs), same discretization, blocks diagonalized by eigh
TWIN_C = [-8.0465386180, 2.6809867361, 2.6558537301, 2.7966837288, -36.18659979, -20.63141616, -14.42364526, -32.38991872, -17.15296823, -11.07440402]
TWIN_F = [-8.0303664488, 2.6830002205, 2.6578928742, 2.7986139754, -36.14944629, -20.61604351, -14.42671231, -32.35678100, -17.14118483, -11.08056169]
# the result vectors: E/A at 0, the three radii at 1 to 3, the six orbital energies from 4
FIRST_RADIUS = 1
FIRST_LEVEL = 4
# skyrme_rpa's published test run, 16O with SLy5 (sample_skyrme_rpa.out)
PUB = [-8.0261107, 2.6836, 2.6585, 2.7991, -36.138, -20.613, -14.429, -32.346, -17.138, -11.083]

twin_diff = 0.0
for k to len(TWIN_C) {
    twin_diff = max(twin_diff, max(abs(res_c[k] - TWIN_C[k]), abs(res_f[k] - TWIN_F[k])))
}
show twin_diff
# what halving the mesh moves: the energy per nucleon and the charge radius
mesh_shift_e = res_f[0] - res_c[0]
mesh_shift_rc = res_f[3] - res_c[3]
show mesh_shift_e
show mesh_shift_rc

test "on both meshes E/A, the radii and the six orbital energies are the numpy twin's to 1e-7" {
    for k to len(TWIN_C) {
        assert_near(res_c[k], TWIN_C[k], 1.0e-7)
        assert_near(res_f[k], TWIN_F[k], 1.0e-7)
    }
}
test "the largest difference from the numpy twin is under 7e-9" {
    assert_near(twin_diff, 0.0, 7.0e-9)
}
test "going from h = 0.1 to 0.05 fm moves E/A by 0.016 MeV and the charge radius by 0.002 fm" {
    assert_near(mesh_shift_e, TWIN_F[0] - TWIN_C[0], 1.0e-7)
    assert_near(mesh_shift_rc, TWIN_F[3] - TWIN_C[3], 1.0e-7)
}

Against the published run

The extrapolated binding energy is $8.0250$ MeV per nucleon, against skyrme_rpa's $8.0261$ : the two differ by $1.1 \times 1 0^{- 3}$ MeV per nucleon, 18 keV out of 128 MeV. skyrme_rpa's own answer moves by $1.2 \times 1 0^{- 4}$ MeV per nucleon between its 0.1 fm and 0.05 fm meshes (the reference script records the rerun), so most of the remaining difference is between the two programs' formulas, not their meshes. The radii agree to $2 \times 1 0^{- 4}$ fm and the six orbital energies to 3 keV. On the 0.1 fm mesh alone, this program's second-order scheme is off by 0.02 MeV per nucleon, which is why it extrapolates.

The figure shows the converged densities on the fine mesh.

[4]lines 522–572

test "extrapolated E/A is skyrme_rpa's -8.0261107 MeV to 2e-3 MeV" {
    assert_near(ext[0], PUB[0], 2.0e-3)
}
test "extrapolated proton, neutron and charge radii are skyrme_rpa's 2.6836, 2.6585 and 2.7991 fm to 3e-4 fm" {
    for k from FIRST_RADIUS to FIRST_LEVEL {
        assert_near(ext[k], PUB[k], 3.0e-4)
    }
}
test "the six extrapolated orbital energies are skyrme_rpa's to 3 keV" {
    for k from FIRST_LEVEL to len(PUB) {
        assert_near(ext[k], PUB[k], 3.0e-3)
    }
}

# skyrme_rpa rerun on 300 points of 0.05 fm (posts/refs): ETot/A = -8.0259916 MeV
PUB_FINE_E_A = -8.0259916
# the extrapolated total energy's distance from the published one, in keV; the
# published run's own shift between its two meshes; and this program's 0.1 fm
# mesh alone against the published run
total_diff_kev = 1000.0 * A_NUC * abs(ext[0] - PUB[0])
pub_mesh_shift = PUB_FINE_E_A - PUB[0]
coarse_off = res_c[0] - PUB[0]
show total_diff_kev
show pub_mesh_shift
show coarse_off

test "the extrapolated energy is 18 keV from skyrme_rpa's, out of 128 MeV" {
    assert_near(total_diff_kev, 18.0, 0.5)
}
test "skyrme_rpa's E/A moves by 1.2e-4 MeV between its 0.1 fm and 0.05 fm meshes" {
    assert_near(pub_mesh_shift, 1.191e-4, 1.0e-10)
}
test "on the 0.1 fm mesh alone the energy per nucleon is off by 0.02 MeV" {
    assert_near(coarse_off, TWIN_C[0] - PUB[0], 1.0e-7)
}

r_f = zeros(nr_f)
rho_n = zeros(nr_f)
rho_p = zeros(nr_f)
for i to nr_f {
    r_f[i] = (i + 1) * H_FINE
    rho_n[i] = dens_f[0, i]
    rho_p[i] = dens_f[1, i]
}
figure("Oxygen-16 with SLy5: neutron and proton densities", "r (fm)", "density (fm$^{-3}$)", 0.0, 6.0, 0.0, 0.0)
plot(r_f, rho_n, "neutrons")
plot(r_f, rho_p, "protons")
style("protons", "dash")

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

test "extrapolated E/A is skyrme_rpa's -8.0261107 MeV to 2e-3 MeV" {
    assert_near(ext[0], PUB[0], 2.0e-3)
}
test "extrapolated proton, neutron and charge radii are skyrme_rpa's 2.6836, 2.6585 and 2.7991 fm to 3e-4 fm" {
    for k from FIRST_RADIUS to FIRST_LEVEL {
        assert_near(ext[k], PUB[k], 3.0e-4)
    }
}
test "the six extrapolated orbital energies are skyrme_rpa's to 3 keV" {
    for k from FIRST_LEVEL to len(PUB) {
        assert_near(ext[k], PUB[k], 3.0e-3)
    }
}

# skyrme_rpa rerun on 300 points of 0.05 fm (posts/refs): ETot/A = -8.0259916 MeV
PUB_FINE_E_A = -8.0259916
# the extrapolated total energy's distance from the published one, in keV; the
# published run's own shift between its two meshes; and this program's 0.1 fm
# mesh alone against the published run
total_diff_kev = 1000.0 * A_NUC * abs(ext[0] - PUB[0])
pub_mesh_shift = PUB_FINE_E_A - PUB[0]
coarse_off = res_c[0] - PUB[0]
show total_diff_kev
show pub_mesh_shift
show coarse_off

test "the extrapolated energy is 18 keV from skyrme_rpa's, out of 128 MeV" {
    assert_near(total_diff_kev, 18.0, 0.5)
}
test "skyrme_rpa's E/A moves by 1.2e-4 MeV between its 0.1 fm and 0.05 fm meshes" {
    assert_near(pub_mesh_shift, 1.191e-4, 1.0e-10)
}
test "on the 0.1 fm mesh alone the energy per nucleon is off by 0.02 MeV" {
    assert_near(coarse_off, TWIN_C[0] - PUB[0], 1.0e-7)
}

r_f = zeros(nr_f)
rho_n = zeros(nr_f)
rho_p = zeros(nr_f)
for i to nr_f {
    r_f[i] = (i + 1) * H_FINE
    rho_n[i] = dens_f[0, i]
    rho_p[i] = dens_f[1, i]
}
figure("Oxygen-16 with SLy5: neutron and proton densities", "r (fm)", "density (fm$^{-3}$)", 0.0, 6.0, 0.0, 0.0)
plot(r_f, rho_n, "neutrons")
plot(r_f, rho_p, "protons")
style("protons", "dash")

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors