# Vautherin-Brink Skyrme-Hartree-Fock across the doubly-magic chain, in Knot.
#
# The generalization of skyrme_vb_o16.knot to heavier nuclei. Same physics
# (SIII, genuine three-body, J^2 terms, half-point Hermitized effective mass,
# np.gradient-matched derivatives, inner/outer Coulomb) -- what changes is the
# machinery: a data-driven shell table (every doubly-magic filling is a prefix
# of one master list), and an eigensolver that can return the k lowest states
# of a tridiagonal block, because heavy nuclei occupy n = 2, 3 radial states.
#
# Deflated inverse iteration: state s is kept orthogonal to states < s after
# every Thomas solve, and each state's shift is warm-started from its own
# eigenvalue at the previous SCF iteration -- so each solve takes a handful
# of O(N) sweeps even where Pb's least-bound levels sit ~1 MeV apart.
#
# Validated against skyrme_hf.py at identical grid/mixing (tol 1e-9, no c.m.):
# nucleus E_knot E_python BE/A
# 16O -113.353 -113.3529 7.0846
# 40Ca -326.345 -326.3451 8.1586
# 48Ca -399.029 -399.0288 8.3131
# 90Zr -763.846 -763.8456 8.4872
# 208Pb -1616.547 -1616.5471 7.7719
#
# Pick the nucleus below. 16O/40Ca/48Ca run in the browser interpreter
# (~5-30 s); 90Zr and 208Pb are better taken native: Download .c, then
# gcc -O2 -o vb program.c -lm && ./vb (fractions of a second)
nucleus = 2.0 # 1 = 16O 2 = 40Ca 3 = 48Ca 4 = 90Zr 5 = 208Pb
pi = 3.141592653589793
hbar2_2m = 20.7355
e2 = 1.439964
# SIII (x1 = x2 = 0, genuine 3-body, J^2 terms on)
t0 = -1128.75
t1 = 395.0
t2 = -95.0
t3 = 14000.0
x0 = 0.45
w0 = 120.0
h = 0.1
nrmax = 512
# ---- per-nucleus configuration ----------------------------------------------
zz = 8.0
nn_neut = 8.0
amass = 16.0
nr = 140
mix = 0.4
if nucleus == 2.0 { zz = 20.0 nn_neut = 20.0 amass = 40.0 nr = 160 mix = 0.3 }
if nucleus == 3.0 { zz = 20.0 nn_neut = 28.0 amass = 48.0 nr = 160 mix = 0.3 }
if nucleus == 4.0 { zz = 40.0 nn_neut = 50.0 amass = 90.0 nr = 200 mix = 0.25 }
if nucleus == 5.0 { zz = 82.0 nn_neut = 126.0 amass = 208.0 nr = 220 mix = 0.25 }
r = zeros(nrmax)
for i to nr { r[i] = (i + 1.0) * h }
# ---- master shell list (spectroscopic filling order; every magic number is a
# prefix): (n, l, 2j) with degeneracy 2j+1 ------------------------------------
m_nq = zeros(22)
m_l = zeros(22)
m_j2 = zeros(22)
m_cnt = zeros(1)
def mput(nq, l, j2) {
i = m_cnt[0]
m_nq[i] = nq
m_l[i] = l
m_j2[i] = j2
m_cnt[0] = i + 1.0
return 0.0
}
mput(1.0, 0.0, 1.0) # 1s1/2 -> 2
mput(1.0, 1.0, 3.0) mput(1.0, 1.0, 1.0) # 1p -> 8
mput(1.0, 2.0, 5.0) mput(2.0, 0.0, 1.0) mput(1.0, 2.0, 3.0) # -> 20
mput(1.0, 3.0, 7.0) # 1f7/2 -> 28
mput(2.0, 1.0, 3.0) mput(1.0, 3.0, 5.0) mput(2.0, 1.0, 1.0) # -> 40
mput(1.0, 4.0, 9.0) # 1g9/2 -> 50
mput(1.0, 4.0, 7.0) mput(2.0, 2.0, 5.0) mput(2.0, 2.0, 3.0)
mput(3.0, 0.0, 1.0) mput(1.0, 5.0, 11.0) # -> 82
mput(1.0, 5.0, 9.0) mput(2.0, 3.0, 7.0) mput(2.0, 3.0, 5.0)
mput(3.0, 1.0, 3.0) mput(3.0, 1.0, 1.0) mput(1.0, 6.0, 13.0) # -> 126
# ---- species orbital tables + (l, j) block tables ----------------------------
# per species: orbital list (prefix of master to the right particle count),
# blocks = unique (l, j2) with kmax = highest radial index n needed.
MAXB = 24
p_ol = zeros(MAXB) p_oj2 = zeros(MAXB) p_onq = zeros(MAXB) p_ocnt = zeros(1)
n_ol = zeros(MAXB) n_oj2 = zeros(MAXB) n_onq = zeros(MAXB) n_ocnt = zeros(1)
p_bl = zeros(MAXB) p_bj2 = zeros(MAXB) p_bk = zeros(MAXB) p_bcnt = zeros(1)
n_bl = zeros(MAXB) n_bj2 = zeros(MAXB) n_bk = zeros(MAXB) n_bcnt = zeros(1)
p_oblk = zeros(MAXB)
n_oblk = zeros(MAXB)
def build_species(target, s_ol, s_oj2, s_onq, s_ocnt, s_bl, s_bj2, s_bk, s_bcnt, s_oblk) {
acc = 0.0
i = 0.0
while acc < target {
s_ol[i] = m_l[i]
s_oj2[i] = m_j2[i]
s_onq[i] = m_nq[i]
acc += m_j2[i] + 1.0
i += 1.0
}
s_ocnt[0] = i
s_bcnt[0] = 0.0
for o to s_ocnt[0] {
found = -1.0
for b to s_bcnt[0] {
if s_bl[b] == s_ol[o] {
if s_bj2[b] == s_oj2[o] { found = b }
}
}
if found < 0.0 {
b = s_bcnt[0]
s_bl[b] = s_ol[o]
s_bj2[b] = s_oj2[o]
s_bk[b] = s_onq[o]
s_bcnt[0] = b + 1.0
found = b
}
if s_onq[o] > s_bk[found] { s_bk[found] = s_onq[o] }
s_oblk[o] = found
}
return 0.0
}
build_species(zz, p_ol, p_oj2, p_onq, p_ocnt, p_bl, p_bj2, p_bk, p_bcnt, p_oblk)
build_species(nn_neut, n_ol, n_oj2, n_onq, n_ocnt, n_bl, n_bj2, n_bk, n_bcnt, n_oblk)
# ---- numerics helpers (identical to the 16O port) ----------------------------
def grad1(f, out, n, hh) {
out[0] = (-1.5 * f[0] + 2.0 * f[1] - 0.5 * f[2]) / hh
for i to n {
if i > 0 {
if i < n - 1.0 { out[i] = (f[i + 1] - f[i - 1]) / (2.0 * hh) }
}
}
out[n - 1.0] = (1.5 * f[n - 1.0] - 2.0 * f[n - 2.0] + 0.5 * f[n - 3.0]) / hh
return 0.0
}
def integrate_r2(f, n, hh) {
s = 0.5 * f[0] * r[0] * r[0] + 0.5 * f[n - 1.0] * r[n - 1.0] * r[n - 1.0]
for i to n {
if i > 0 {
if i < n - 1.0 { s += f[i] * r[i] * r[i] }
}
}
return 4.0 * pi * s * hh
}
# ---- k lowest states of a tridiagonal block ----------------------------------
# Deflated inverse iteration. ev[s*nrmax + i] holds state s; eps_k[s] its
# eigenvalue. shift_in[s] < 0 warm-starts state s near its previous eigenvalue;
# pass +1e9 to fall back to the Gershgorin bound (first SCF iteration).
tdl_c = zeros(512)
tdl_d = zeros(512)
u_cur = zeros(512)
u_new = zeros(512)
ev = zeros(1536) # 3 * 512: up to k = 3 radial states per block
eps_k = zeros(3)
def thomas_solve(diag, off, n, sigma) {
tdl_c[0] = off[0] / (diag[0] - sigma)
tdl_d[0] = u_cur[0] / (diag[0] - sigma)
for i to n {
if i > 0 {
m = diag[i] - sigma - off[i - 1.0] * tdl_c[i - 1.0]
if i < n - 1.0 { tdl_c[i] = off[i] / m }
tdl_d[i] = (u_cur[i] - off[i - 1.0] * tdl_d[i - 1.0]) / m
}
}
u_new[n - 1.0] = tdl_d[n - 1.0]
k = n - 2.0
while k >= 0.0 {
u_new[k] = tdl_d[k] - tdl_c[k] * u_new[k + 1.0]
k = k - 1.0
}
return 0.0
}
def project_normalize(n, s) {
# remove components along states < s, then normalize u_new
for t to s {
dot = 0.0
for i to n { dot += ev[t * nrmax + i] * u_new[i] }
for i to n { u_new[i] = u_new[i] - dot * ev[t * nrmax + i] }
}
nrm = 0.0
for i to n { nrm += u_new[i] * u_new[i] }
nrm = sqrt(nrm)
for i to n { u_new[i] = u_new[i] / nrm }
return 0.0
}
def rayleigh(diag, off, n) {
lam = 0.0
for i to n {
lam += diag[i] * u_new[i] * u_new[i]
if i < n - 1.0 { lam += 2.0 * off[i] * u_new[i] * u_new[i + 1.0] }
}
return lam
}
def lowest_k(diag, off, n, kk, shift_in) {
# Gershgorin lower bound (the kinetic diagonal cancels against the
# off-diagonals, so this lands near the potential floor -- a tight shift)
gersh = diag[0]
for i to n {
g = diag[i]
if i > 0 { g = g - abs(off[i - 1.0]) }
if i < n - 1.0 { g = g - abs(off[i]) }
if g < gersh { gersh = g }
}
gersh = gersh - 1.0
for s to kk {
sigma = shift_in[s] - 0.25
if shift_in[s] > 1.0e8 { sigma = gersh }
# sin profile with s nodes: overlaps the target state well
for i to n { u_new[i] = sin((s + 1.0) * pi * (i + 1.0) / (n + 1.0)) }
project_normalize(n, s)
# The shift stays FIXED for the whole solve. With sigma below the
# spectrum (Gershgorin) deflated inverse iteration provably converges
# to the lowest remaining state; with a warm-start sigma it sits ~0.25
# MeV from its own state while the nearest same-block state is an
# oscillator spacing (~7 MeV) away. A moving Rayleigh shift, tried
# first, could land between Pb's crowded levels mid-flight and capture
# the wrong state -- which the warm start then locks in for good.
lam = 0.0
for sweep to 400 {
for i to n { u_cur[i] = u_new[i] }
thomas_solve(diag, off, n, sigma)
project_normalize(n, s)
newlam = rayleigh(diag, off, n)
dl = abs(newlam - lam)
lam = newlam
if sweep > 2.0 {
if dl < 1.0e-11 { break }
}
}
eps_k[s] = lam
for i to n { ev[s * nrmax + i] = u_new[i] }
}
# bubble-sort the k states ascending (deflation can misorder warm starts)
for a to kk {
for b to kk {
if b > 0.0 {
if eps_k[b - 1.0] > eps_k[b] {
tmpe = eps_k[b - 1.0]
eps_k[b - 1.0] = eps_k[b]
eps_k[b] = tmpe
for i to n {
tmpv = ev[(b - 1.0) * nrmax + i]
ev[(b - 1.0) * nrmax + i] = ev[b * nrmax + i]
ev[b * nrmax + i] = tmpv
}
}
}
}
}
return 0.0
}
# ---- block Hamiltonian --------------------------------------------------------
bh_diag = zeros(512)
bh_off = zeros(512)
def build_block_h(l, sdl, B, U, W, Vc) {
for i to nr {
bm = B[0]
if i > 0 { bm = 0.5 * (B[i - 1.0] + B[i]) }
bp = 0.5 * (B[i] + hbar2_2m)
if i < nr - 1.0 { bp = 0.5 * (B[i] + B[i + 1.0]) }
bh_diag[i] = (bm + bp) / (h * h) + B[i] * l * (l + 1.0) / (r[i] * r[i]) + U[i] + Vc[i] + sdl * W[i] / r[i]
if i < nr - 1.0 { bh_off[i] = -bp / (h * h) }
}
return 0.0
}
# ---- densities ------------------------------------------------------------------
dR_work = zeros(512)
R_tmp = zeros(512)
def add_orbital_R(deg, l, sdl, rho_a, tau_a, jso_a) {
w = deg / (4.0 * pi)
grad1(R_tmp, dR_work, nr, h)
for i to nr {
rho_a[i] = rho_a[i] + w * R_tmp[i] * R_tmp[i]
tau_a[i] = tau_a[i] + w * (dR_work[i] * dR_work[i] + l * (l + 1.0) / (r[i] * r[i]) * R_tmp[i] * R_tmp[i])
jso_a[i] = jso_a[i] + w * sdl / r[i] * R_tmp[i] * R_tmp[i]
}
return 0.0
}
def coulomb_direct(rp, out) {
cum = 0.0
for i to nr {
cum += rp[i] * r[i] * r[i]
out[i] = cum / r[i]
}
rev = 0.0
k = nr - 1.0
while k >= 0.0 {
out[k] = e2 * 4.0 * pi * h * (out[k] + rev)
rev += rp[k] * r[k]
k = k - 1.0
}
return 0.0
}
# ---- state ----------------------------------------------------------------------
rho_n = zeros(nrmax)
rho_p = zeros(nrmax)
tau_n = zeros(nrmax)
tau_p = zeros(nrmax)
j_n = zeros(nrmax)
j_p = zeros(nrmax)
r0ws = 1.2 * pow(amass, 1.0 / 3.0)
for i to nr { rho_n[i] = 1.0 / (1.0 + exp((r[i] - r0ws) / 0.6)) }
prof_norm = integrate_r2(rho_n, nr, h)
for i to nr {
p = rho_n[i]
rho_n[i] = p * nn_neut / prof_norm
rho_p[i] = p * zz / prof_norm
}
rho = zeros(nrmax)
tau = zeros(nrmax)
jtot = zeros(nrmax)
drho = zeros(nrmax)
drho_n = zeros(nrmax)
drho_p = zeros(nrmax)
lap = zeros(nrmax)
lap_n = zeros(nrmax)
lap_p = zeros(nrmax)
d1 = zeros(nrmax)
d2 = zeros(nrmax)
b_n = zeros(nrmax)
b_p = zeros(nrmax)
u_n = zeros(nrmax)
u_p = zeros(nrmax)
w_n = zeros(nrmax)
w_p = zeros(nrmax)
vcoul = zeros(nrmax)
vzero = zeros(nrmax)
sso = zeros(nrmax)
dsso = zeros(nrmax)
nrho_n = zeros(nrmax)
nrho_p = zeros(nrmax)
ntau_n = zeros(nrmax)
ntau_p = zeros(nrmax)
nj_n = zeros(nrmax)
nj_p = zeros(nrmax)
edens = zeros(nrmax)
# warm-start shift memory: per species, per block, per state (MAXB * 3)
p_shift = zeros(72)
n_shift = zeros(72)
for i to 72 { p_shift[i] = 1.0e9 n_shift[i] = 1.0e9 }
# level output: eigenvalue per orbital
lev_n = zeros(MAXB)
lev_p = zeros(MAXB)
def laplacian(f) {
grad1(f, d1, nr, h)
grad1(d1, d2, nr, h)
for i to nr { lap[i] = d2[i] + 2.0 / r[i] * d1[i] }
return 0.0
}
# solve all blocks of one species, rebuild its densities, record levels
def species_step(s_bl, s_bj2, s_bk, s_bcnt, s_ol, s_oj2, s_onq, s_ocnt, s_oblk,
B, U, W, Vc, shifts, rho_a, tau_a, jso_a, lev) {
for i to nr {
rho_a[i] = 0.0
tau_a[i] = 0.0
jso_a[i] = 0.0
}
for b to s_bcnt[0] {
l = s_bl[b]
jj = s_bj2[b] / 2.0
sdl = jj * (jj + 1.0) - l * (l + 1.0) - 0.75
build_block_h(l, sdl, B, U, W, Vc)
kk = s_bk[b]
for s to kk { eps_k[s] = 0.0 }
# load warm shifts for this block
sh3 = zeros(3)
for s to kk { sh3[s] = shifts[b * 3.0 + s] }
lowest_k(bh_diag, bh_off, nr, kk, sh3)
for s to kk { shifts[b * 3.0 + s] = eps_k[s] }
# distribute states to this block's orbitals
for o to s_ocnt[0] {
if s_oblk[o] == b {
st = s_onq[o] - 1.0
for i to nr { R_tmp[i] = ev[st * nrmax + i] / sqrt(h) / r[i] }
add_orbital_R(s_oj2[o] + 1.0, l, sdl, rho_a, tau_a, jso_a)
lev[o] = eps_k[st]
}
}
}
return 0.0
}
# ---- SCF --------------------------------------------------------------------
e_old = 0.0
e_total = 0.0
for iter to 400 {
for i to nr {
rho[i] = rho_n[i] + rho_p[i]
tau[i] = tau_n[i] + tau_p[i]
jtot[i] = j_n[i] + j_p[i]
}
grad1(rho, drho, nr, h)
grad1(rho_n, drho_n, nr, h)
grad1(rho_p, drho_p, nr, h)
laplacian(rho_n)
for i to nr { lap_n[i] = lap[i] }
laplacian(rho_p)
for i to nr { lap_p[i] = lap[i] }
laplacian(rho)
for i to nr {
b_n[i] = hbar2_2m + 0.125 * (2.0 * t1 + 2.0 * t2) * rho[i] + 0.125 * (t2 - t1) * rho_n[i]
b_p[i] = hbar2_2m + 0.125 * (2.0 * t1 + 2.0 * t2) * rho[i] + 0.125 * (t2 - t1) * rho_p[i]
u_n[i] = t0 * ((1.0 + x0 / 2.0) * rho[i] - (x0 + 0.5) * rho_n[i])
u_n[i] = u_n[i] + 0.25 * t1 * ((tau[i] - 1.5 * lap[i]) - 0.5 * (tau_n[i] - 1.5 * lap_n[i]))
u_n[i] = u_n[i] + 0.25 * t2 * ((tau[i] + 0.5 * lap[i]) + 0.5 * (tau_n[i] + 0.5 * lap_n[i]))
u_n[i] = u_n[i] + 0.25 * t3 * rho_p[i] * (2.0 * rho_n[i] + rho_p[i])
u_p[i] = t0 * ((1.0 + x0 / 2.0) * rho[i] - (x0 + 0.5) * rho_p[i])
u_p[i] = u_p[i] + 0.25 * t1 * ((tau[i] - 1.5 * lap[i]) - 0.5 * (tau_p[i] - 1.5 * lap_p[i]))
u_p[i] = u_p[i] + 0.25 * t2 * ((tau[i] + 0.5 * lap[i]) + 0.5 * (tau_p[i] + 0.5 * lap_p[i]))
u_p[i] = u_p[i] + 0.25 * t3 * rho_n[i] * (2.0 * rho_p[i] + rho_n[i])
}
for i to nr { sso[i] = jtot[i] + j_n[i] }
grad1(sso, dsso, nr, h)
for i to nr { u_n[i] = u_n[i] - 0.5 * w0 * (sso[i] / r[i] + 0.5 * dsso[i]) }
for i to nr { sso[i] = jtot[i] + j_p[i] }
grad1(sso, dsso, nr, h)
for i to nr { u_p[i] = u_p[i] - 0.5 * w0 * (sso[i] / r[i] + 0.5 * dsso[i]) }
for i to nr {
w_n[i] = 0.5 * w0 * (drho[i] + drho_n[i]) + 0.125 * (t1 - t2) * j_n[i]
w_p[i] = 0.5 * w0 * (drho[i] + drho_p[i]) + 0.125 * (t1 - t2) * j_p[i]
}
coulomb_direct(rho_p, vcoul)
for i to nr { vcoul[i] = vcoul[i] - e2 * pow(3.0 / pi, 1.0 / 3.0) * pow(rho_p[i], 1.0 / 3.0) }
species_step(n_bl, n_bj2, n_bk, n_bcnt, n_ol, n_oj2, n_onq, n_ocnt, n_oblk,
b_n, u_n, w_n, vzero, n_shift, nrho_n, ntau_n, nj_n, lev_n)
species_step(p_bl, p_bj2, p_bk, p_bcnt, p_ol, p_oj2, p_onq, p_ocnt, p_oblk,
b_p, u_p, w_p, vcoul, p_shift, nrho_p, ntau_p, nj_p, lev_p)
for i to nr {
rho_n[i] = (1.0 - mix) * rho_n[i] + mix * nrho_n[i]
rho_p[i] = (1.0 - mix) * rho_p[i] + mix * nrho_p[i]
tau_n[i] = (1.0 - mix) * tau_n[i] + mix * ntau_n[i]
tau_p[i] = (1.0 - mix) * tau_p[i] + mix * ntau_p[i]
j_n[i] = (1.0 - mix) * j_n[i] + mix * nj_n[i]
j_p[i] = (1.0 - mix) * j_p[i] + mix * nj_p[i]
}
for i to nr {
rho[i] = rho_n[i] + rho_p[i]
tau[i] = tau_n[i] + tau_p[i]
}
grad1(rho, drho, nr, h)
grad1(rho_n, drho_n, nr, h)
grad1(rho_p, drho_p, nr, h)
for i to nr {
hd = hbar2_2m * tau[i]
hd += 0.5 * t0 * ((1.0 + x0 / 2.0) * rho[i] * rho[i] - (x0 + 0.5) * (rho_n[i] * rho_n[i] + rho_p[i] * rho_p[i]))
hd += 0.25 * t1 * ((rho[i] * tau[i] + 0.75 * drho[i] * drho[i])
- 0.5 * (rho_n[i] * tau_n[i] + 0.75 * drho_n[i] * drho_n[i]
+ rho_p[i] * tau_p[i] + 0.75 * drho_p[i] * drho_p[i]))
hd += 0.25 * t2 * ((rho[i] * tau[i] - 0.25 * drho[i] * drho[i])
+ 0.5 * (rho_n[i] * tau_n[i] - 0.25 * drho_n[i] * drho_n[i]
+ rho_p[i] * tau_p[i] - 0.25 * drho_p[i] * drho_p[i]))
hd += 0.25 * t3 * rho_n[i] * rho_p[i] * rho[i]
hd += (1.0 / 16.0) * (t1 - t2) * (j_n[i] * j_n[i] + j_p[i] * j_p[i])
hd += 0.5 * w0 * ((j_n[i] + j_p[i]) * drho[i] + j_n[i] * drho_n[i] + j_p[i] * drho_p[i])
edens[i] = hd
}
e_nuc = integrate_r2(edens, nr, h)
coulomb_direct(rho_p, vcoul)
for i to nr { edens[i] = 0.5 * vcoul[i] * rho_p[i] }
e_cdir = integrate_r2(edens, nr, h)
for i to nr { edens[i] = pow(rho_p[i], 4.0 / 3.0) }
e_cex = -0.75 * e2 * pow(3.0 / pi, 1.0 / 3.0) * integrate_r2(edens, nr, h)
e_total = e_nuc + e_cdir + e_cex
de = e_total - e_old
e_old = e_total
if abs(de) < 1.0e-9 {
if iter > 3.0 { break }
}
}
# ---- results ------------------------------------------------------------------
show e_total
be_per_a = 0.0 - e_total / amass
show be_per_a
for i to nr { edens[i] = rho_p[i] * r[i] * r[i] }
rms_p = sqrt(integrate_r2(edens, nr, h) / zz)
show rms_p
# single-particle levels (MeV), in shell-filling order -- plotted as figures
show lev_n
show lev_p
# converged densities
show rho_n
show rho_p