# Vautherin-Brink Skyrme-Hartree-Fock for 16O, in Knot.
#
# A port of skyrme_hf.py (browsable at /code/skyrme-hf-vb): spherical nuclear
# DFT with the SIII force -- genuine three-body term, J^2 tensor terms on,
# radial mesh, self-consistency by linear density mixing. Every algorithmic
# choice mirrors the Python (np.gradient edge_order=2 derivatives, the same
# half-point Hermitized effective-mass stencil, the same inner/outer Coulomb
# partition), so the converged numbers agree to ~1e-3 MeV:
#
# E_total = -113.353 MeV BE/A = 7.085 MeV (rmax=14, h=0.1, no c.m. corr.)
# levels: 1s1/2 -33.59 1p3/2 -18.60 1p1/2 -12.91 MeV
#
# Where the Python diagonalizes each (l,j) block densely with numpy.eigh, this
# port exploits that the radial Hamiltonian is TRIDIAGONAL and that every 16O
# orbital is the lowest state of its block: a Gershgorin shift makes H - sigma*I
# positive definite, and inverse iteration (Thomas solves, O(N) each) converges
# to exactly that lowest eigenpair. No dense matrices anywhere.
pi = 3.141592653589793
hbar2_2m = 20.7355 # hbar^2/2m [MeV fm^2]
e2 = 1.439964 # e^2 [MeV fm]
# SIII (Beiner 1975): t0,t1,t2,t3, x0, W0; x1=x2=0, genuine 3-body, J^2 on
t0 = -1128.75
t1 = 395.0
t2 = -95.0
t3 = 14000.0
x0 = 0.45
w0 = 120.0
# grid: r_i = (i+1) h, excludes the origin (u(0)=0 boundary handles r=0)
h = 0.1
nr = 140 # rmax = 14 fm
mix = 0.4
zz = 8.0 # protons
nn_neut = 8.0 # neutrons
amass = 16.0
r = zeros(nr)
for i to nr {
r[i] = (i + 1.0) * h
}
# ---- np.gradient with edge_order=2, exactly (grad1: deriv is a Knot builtin) ---------------------------------
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
}
# trapezoid of f r^2 on the mesh, times 4 pi (the radial volume integral)
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
}
# plain trapezoid of g on the mesh (g already carries its r factors)
def trapz(g, n, hh) {
s = 0.5 * g[0] + 0.5 * g[n - 1.0]
for i to n {
if i > 0 {
if i < n - 1.0 {
s += g[i]
}
}
}
return s * hh
}
# ---- lowest eigenpair of the tridiagonal block ------------------------------
# diag[], off[] (off[i] couples i,i+1). Gershgorin lower bound as shift, then
# inverse iteration via the Thomas algorithm; Rayleigh quotient at the end.
# u_out is normalized so sum u^2 = 1 (matches an eigh unit eigenvector).
tdl_c = zeros(512)
tdl_d = zeros(512)
u_work = zeros(512)
def lowest_eigenpair(diag, off, n, u_out) {
# shift strictly below the spectrum
sigma = diag[0]
for i to n {
g = diag[i]
reach = 0.0
if i > 0 { reach += abs(off[i - 1.0]) }
if i < n - 1.0 { reach += abs(off[i]) }
g = g - reach
if g < sigma { sigma = g }
}
sigma = sigma - 1.0
for i to n { u_out[i] = 1.0 }
lam = 0.0
for sweep to 200 {
# Thomas solve (A - sigma I) w = u -> w overwrites u_work
tdl_c[0] = off[0] / (diag[0] - sigma)
tdl_d[0] = u_out[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_out[i] - off[i - 1.0] * tdl_d[i - 1.0]) / m
}
}
u_work[n - 1.0] = tdl_d[n - 1.0]
k = n - 2.0
while k >= 0.0 {
u_work[k] = tdl_d[k] - tdl_c[k] * u_work[k + 1.0]
k = k - 1.0
}
# normalize
nrm = 0.0
for i to n { nrm += u_work[i] * u_work[i] }
nrm = sqrt(nrm)
# convergence: vector change (up to sign)
dot = 0.0
for i to n { dot += u_work[i] / nrm * u_out[i] }
sgn = 1.0
if dot < 0.0 { sgn = -1.0 }
change = 0.0
for i to n {
d = sgn * u_work[i] / nrm - u_out[i]
change += d * d
u_out[i] = sgn * u_work[i] / nrm
}
if change < 1.0e-24 { break }
}
# Rayleigh quotient lam = u^T H u
lam = 0.0
for i to n {
lam += diag[i] * u_out[i] * u_out[i]
if i < n - 1.0 { lam += 2.0 * off[i] * u_out[i] * u_out[i + 1.0] }
}
return lam
}
# ---- one (l, j) block: build H, get lowest state ----------------------------
# B, U, W, Vc are the fields; writes R = u/(r sqrt(h)) into R_out, returns eps.
bh_diag = zeros(512)
bh_off = zeros(512)
bh_u = zeros(512)
def solve_block(l, sdl, B, U, W, Vc, R_out) {
# half-point effective mass: Bm_i = (B_{i-1}+B_i)/2 with the same padding
# as the Python (B[-1] pad = B[0]; B[N] pad = hbar^2/2m vacuum)
for i to nr {
bm = 0.5 * (B[i] + B[i]) # placeholder, fixed below
if i > 0 { bm = 0.5 * (B[i - 1.0] + B[i]) }
bp = hbar2_2m
if i < nr - 1.0 { bp = 0.5 * (B[i] + B[i + 1.0]) }
if i == nr - 1.0 { bp = 0.5 * (B[i] + hbar2_2m) }
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) }
}
eps = lowest_eigenpair(bh_diag, bh_off, nr, bh_u)
for i to nr {
R_out[i] = bh_u[i] / sqrt(h) / r[i]
}
return eps
}
# ---- densities from the three occupied orbitals of one isospin --------------
# orbitals: 1s1/2 (deg 2, sdl 0), 1p3/2 (deg 4, sdl 1), 1p1/2 (deg 2, sdl -2)
dR_work = zeros(512)
def add_orbital(R, deg, l, sdl, rho, tau, jso) {
w = deg / (4.0 * pi)
grad1(R, dR_work, nr, h)
for i to nr {
rho[i] = rho[i] + w * R[i] * R[i]
tau[i] = tau[i] + w * (dR_work[i] * dR_work[i] + l * (l + 1.0) / (r[i] * r[i]) * R[i] * R[i])
jso[i] = jso[i] + w * sdl / r[i] * R[i] * R[i]
}
return 0.0
}
# ---- direct Coulomb: clean inner/outer shell partition -----------------------
def coulomb_direct(rho_p, out) {
cum = 0.0
for i to nr {
cum += rho_p[i] * r[i] * r[i]
out[i] = cum / r[i] # inner sum / r
}
rev = 0.0
k = nr - 1.0
while k >= 0.0 {
out[k] = e2 * 4.0 * pi * h * (out[k] + rev)
rev += rho_p[k] * r[k] # shells strictly beyond r_k, next round
k = k - 1.0
}
return 0.0
}
# ---- state ------------------------------------------------------------------
rho_n = zeros(nr)
rho_p = zeros(nr)
tau_n = zeros(nr)
tau_p = zeros(nr)
j_n = zeros(nr)
j_p = zeros(nr)
# Woods-Saxon initial guess, normalized to N and Z
r0ws = 1.2 * pow(amass, 1.0 / 3.0)
prof_norm = 0.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
}
# work arrays
rho = zeros(nr)
tau = zeros(nr)
jtot = zeros(nr)
drho = zeros(nr)
drho_n = zeros(nr)
drho_p = zeros(nr)
lap = zeros(nr)
lap_n = zeros(nr)
lap_p = zeros(nr)
d1 = zeros(nr)
d2 = zeros(nr)
b_n = zeros(nr)
b_p = zeros(nr)
u_n = zeros(nr)
u_p = zeros(nr)
w_n = zeros(nr)
w_p = zeros(nr)
vcoul = zeros(nr)
vzero = zeros(nr)
sso = zeros(nr)
dsso = zeros(nr)
R_s = zeros(nr)
R_p3 = zeros(nr)
R_p1 = zeros(nr)
nrho_n = zeros(nr)
nrho_p = zeros(nr)
ntau_n = zeros(nr)
ntau_p = zeros(nr)
nj_n = zeros(nr)
nj_p = zeros(nr)
edens = zeros(nr)
e_old = 0.0
e_total = 0.0
eps_s = 0.0
eps_p3 = 0.0
eps_p1 = 0.0
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
}
# ---- SCF --------------------------------------------------------------------
for iter to 300 {
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) # total density last: stays in lap[]
# fields (x1 = x2 = 0 throughout: SIII)
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])
}
# spin-orbit central term: -(W0/2)[ (J+Jq)/r + 1/2 d/dr(J+Jq) ]
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]) }
# spin-orbit form factor W_q = (W0/2)(drho + drho_q) + (t1-t2)/8 J_q
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 partition + Slater exchange, protons only
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)
}
# solve the three blocks per isospin (all n=1 -> lowest state of each)
eps_s = solve_block(0.0, 0.0, b_n, u_n, w_n, vzero, R_s)
eps_p3 = solve_block(1.0, 1.0, b_n, u_n, w_n, vzero, R_p3)
eps_p1 = solve_block(1.0, -2.0, b_n, u_n, w_n, vzero, R_p1)
for i to nr {
nrho_n[i] = 0.0
ntau_n[i] = 0.0
nj_n[i] = 0.0
}
add_orbital(R_s, 2.0, 0.0, 0.0, nrho_n, ntau_n, nj_n)
add_orbital(R_p3, 4.0, 1.0, 1.0, nrho_n, ntau_n, nj_n)
add_orbital(R_p1, 2.0, 1.0, -2.0, nrho_n, ntau_n, nj_n)
ep_s = solve_block(0.0, 0.0, b_p, u_p, w_p, vcoul, R_s)
ep_p3 = solve_block(1.0, 1.0, b_p, u_p, w_p, vcoul, R_p3)
ep_p1 = solve_block(1.0, -2.0, b_p, u_p, w_p, vcoul, R_p1)
for i to nr {
nrho_p[i] = 0.0
ntau_p[i] = 0.0
nj_p[i] = 0.0
}
add_orbital(R_s, 2.0, 0.0, 0.0, nrho_p, ntau_p, nj_p)
add_orbital(R_p3, 4.0, 1.0, 1.0, nrho_p, ntau_p, nj_p)
add_orbital(R_p1, 2.0, 1.0, -2.0, nrho_p, ntau_p, nj_p)
# linear mixing
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]
}
# total energy from the functional (no c.m. correction)
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-7 {
if iter > 3.0 { break }
}
}
# ---- results -----------------------------------------------------------------
show e_total # -113.353 MeV (Python: -113.352915)
be_per_a = 0.0 - e_total / amass
show be_per_a # 7.085 MeV/A
# neutron single-particle levels (MeV)
show eps_s # 1s1/2 -33.59
show eps_p3 # 1p3/2 -18.60
show eps_p1 # 1p1/2 -12.91
# rms proton radius
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 # 2.692 fm
# densities: plotted as figures by the IDE
show rho_n
show rho_p