# Build your own nuclear DFT -- step 5: the Skyrme force becomes fields.
#
# Here the fixed Woods-Saxon well retires. In Skyrme-Hartree-Fock the potential
# is BUILT FROM THE DENSITIES: the energy density H(rho, tau, J) is an ordinary
# function of the local densities, and differentiating it hands each nucleon
# three fields,
#
# U(r) = dH/drho central potential
# B(r) = dH/dtau = hbar^2/2m*(r), a position-dependent effective mass
# W(r) = dH/dJ spin-orbit form factor
#
# This step computes all three from step 4's densities (SIII force, N = Z, no
# Coulomb yet) and shows the two signatures worth staring at:
# - U(r) comes out ~ -90 MeV deep with a Woods-Saxon-ish shape NOBODY put in;
# the flat-bottomed well EMERGES from the force + the density (deeper than
# the self-consistent one, because step 4 densities are too compact).
# - m*/m = hbar^2/2m / B drops to ~0.75 inside: nucleons in nuclear matter
# act lighter. That number controls level spacings.
#
# The derivative bookkeeping: the t1/t2 terms need the laplacian of rho, and
# W needs its gradient -- hence grad1 applied twice for lap = f'' + (2/r) f'.
pi = 3.141592653589793
hbar2_2m = 20.7355
h = 0.1
nr = 140
amass = 16.0
# SIII
t0 = -1128.75
t1 = 395.0
t2 = -95.0
t3 = 14000.0
x0 = 0.45
w0 = 120.0
r = zeros(nr)
for i to nr { r[i] = (i + 1.0) * h }
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
}
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
}
# ---- solve the fixed well once, exactly as in step 4, to get densities -------
v0 = -51.0
rws = 1.27 * pow(amass, 1.0 / 3.0)
aws = 0.67
u_pot = zeros(nr)
wso = zeros(nr)
for i to nr {
f = 1.0 / (1.0 + exp((r[i] - rws) / aws))
u_pot[i] = v0 * f
wso[i] = 22.0 * (0.0 - f * (1.0 - f) / aws)
}
diag = zeros(nr)
off = zeros(nr)
u_new = zeros(nr)
u_cur = zeros(nr)
tc = zeros(nr)
td = zeros(nr)
ev = zeros(1024)
eps_k = zeros(2)
def build_h(l, sdl) {
for i to nr {
diag[i] = 2.0 * hbar2_2m / (h * h) + hbar2_2m * l * (l + 1.0) / (r[i] * r[i]) + u_pot[i] + sdl * wso[i] / r[i]
if i < nr - 1.0 { off[i] = 0.0 - hbar2_2m / (h * h) }
}
return 0.0
}
def lowest_k(n, kk) {
sigma = 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 < sigma { sigma = g }
}
sigma = sigma - 1.0
for s to kk {
for i to n { u_new[i] = sin((s + 1.0) * pi * (i + 1.0) / (n + 1.0)) }
lam = 0.0
for sweep to 300 {
for i to n { u_cur[i] = u_new[i] }
tc[0] = off[0] / (diag[0] - sigma)
td[0] = u_cur[0] / (diag[0] - sigma)
for i to n {
if i > 0 {
m = diag[i] - sigma - off[i - 1.0] * tc[i - 1.0]
if i < n - 1.0 { tc[i] = off[i] / m }
td[i] = (u_cur[i] - off[i - 1.0] * td[i - 1.0]) / m
}
}
u_new[n - 1.0] = td[n - 1.0]
k = n - 2.0
while k >= 0.0 {
u_new[k] = td[k] - tc[k] * u_new[k + 1.0]
k = k - 1.0
}
for t to s {
dot = 0.0
for i to n { dot += ev[t * 512.0 + i] * u_new[i] }
for i to n { u_new[i] = u_new[i] - dot * ev[t * 512.0 + 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 }
newlam = 0.0
for i to n {
newlam += diag[i] * u_new[i] * u_new[i]
if i < n - 1.0 { newlam += 2.0 * off[i] * u_new[i] * u_new[i + 1.0] }
}
dl = abs(newlam - lam)
lam = newlam
if sweep > 2.0 {
if dl < 1.0e-12 { break }
}
}
eps_k[s] = lam
for i to n { ev[s * 512.0 + i] = u_new[i] }
}
return 0.0
}
rho = zeros(nr)
tau = zeros(nr)
jso = zeros(nr)
R_tmp = zeros(nr)
dR = zeros(nr)
def add_orbital(deg, l, sdl) {
for i to nr { R_tmp[i] = ev[i] / sqrt(h) / r[i] }
grad1(R_tmp, dR, nr, h)
w = deg / (4.0 * pi)
for i to nr {
rho[i] = rho[i] + w * R_tmp[i] * R_tmp[i]
tau[i] = tau[i] + w * (dR[i] * dR[i] + l * (l + 1.0) / (r[i] * r[i]) * R_tmp[i] * R_tmp[i])
jso[i] = jso[i] + w * sdl / r[i] * R_tmp[i] * R_tmp[i]
}
return 0.0
}
build_h(0.0, 0.0)
lowest_k(nr, 1.0)
add_orbital(2.0, 0.0, 0.0)
build_h(1.0, 1.0)
lowest_k(nr, 1.0)
add_orbital(4.0, 1.0, 1.0)
build_h(1.0, -2.0)
lowest_k(nr, 1.0)
add_orbital(2.0, 1.0, -2.0)
for i to nr {
rho[i] = 2.0 * rho[i]
tau[i] = 2.0 * tau[i]
jso[i] = 2.0 * jso[i]
}
# ---- NEW: the Skyrme fields from (rho, tau, J), N = Z symmetric --------------
# with rho_n = rho_p = rho/2, tau_q = tau/2, J_q = J/2 and no Coulomb.
drho = zeros(nr)
d1 = zeros(nr)
d2 = zeros(nr)
lap = zeros(nr)
grad1(rho, drho, nr, h)
grad1(rho, d1, nr, h)
grad1(d1, d2, nr, h)
for i to nr { lap[i] = d2[i] + 2.0 / r[i] * d1[i] }
b_field = zeros(nr)
u_field = zeros(nr)
w_field = zeros(nr)
mstar = zeros(nr)
djso = zeros(nr)
grad1(jso, djso, nr, h)
for i to nr {
# effective mass field: dH/dtau
b_field[i] = hbar2_2m + 0.125 * (2.0 * t1 + 2.0 * t2) * rho[i] + 0.125 * (t2 - t1) * 0.5 * rho[i]
mstar[i] = hbar2_2m / b_field[i]
# central field: dH/drho (t0 + t1 + t2 + genuine three-body terms)
u_field[i] = t0 * ((1.0 + x0 / 2.0) * rho[i] - (x0 + 0.5) * 0.5 * rho[i])
u_field[i] = u_field[i] + 0.25 * t1 * ((tau[i] - 1.5 * lap[i]) - 0.5 * (0.5 * tau[i] - 0.75 * lap[i]))
u_field[i] = u_field[i] + 0.25 * t2 * ((tau[i] + 0.5 * lap[i]) + 0.5 * (0.5 * tau[i] + 0.25 * lap[i]))
u_field[i] = u_field[i] + 0.25 * t3 * 0.5 * rho[i] * (rho[i] + 0.5 * rho[i])
# spin-orbit central piece: s = J + J_q = 1.5 J
u_field[i] = u_field[i] - 0.5 * w0 * (1.5 * jso[i] / r[i] + 0.75 * djso[i])
# spin-orbit form factor
w_field[i] = 0.5 * w0 * (drho[i] + 0.5 * drho[i]) + 0.125 * (t1 - t2) * 0.5 * jso[i]
}
u_depth = u_field[0]
show u_depth # ~ -90 MeV: the well the force digs, nobody drew it
mstar_center = mstar[0]
show mstar_center # ~ 0.75: SIII's effective mass, matching the literature
show u_field # Woods-Saxon-shaped -- EMERGED from the density
show mstar # dips inside the nucleus, 1 in the vacuum
show w_field # surface-peaked, like the hand-drawn one in step 3