Step 1 — The grid and a guess density
Nuclear Physics
Everything lives on a uniform radial mesh with fm that deliberately excludes the origin — the terms coming in step 3 would blow up there, and the boundary condition handles for us. The only other tool: a volume integral. For a spherically symmetric function, , done by the trapezoid rule. That's enough to write down a nucleus-shaped guess — the Woods–Saxon profile, flat interior and soft surface — normalized so it holds exactly 16 nucleons.
The program
# Build your own nuclear DFT -- step 1: the grid and a guess density.
#
# Everything lives on a uniform radial mesh r_i = (i+1)h that EXCLUDES the
# origin: the 1/r terms coming later would blow up at r = 0, and the boundary
# condition u(0) = 0 handles the origin for us. The only other tool we need is
# a volume integral: for a spherically symmetric f(r),
# integral f d^3r = 4 pi integral f(r) r^2 dr (trapezoid rule).
#
# With those two pieces we can already write down a nucleus-shaped guess: the
# Woods-Saxon profile 1/(1 + exp((r-R0)/a)) -- flat interior, soft surface --
# normalized so it holds exactly A = 16 nucleons.
pi = 3.141592653589793
h = 0.1 # grid spacing [fm]
nr = 140 # 14 fm box
amass = 16.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
}
# Woods-Saxon shape: R0 = 1.2 A^(1/3) fm is the nuclear radius, a = 0.6 fm the
# surface thickness. Normalize it to A particles.
r0ws = 1.2 * pow(amass, 1.0 / 3.0)
rho = zeros(nr)
for i to nr { rho[i] = 1.0 / (1.0 + exp((r[i] - r0ws) / 0.6)) }
norm = integrate_r2(rho, nr, h)
for i to nr { rho[i] = rho[i] * amass / norm }
# check: integrating the density must count the nucleons
a_check = integrate_r2(rho, nr, h)
show a_check # 16 -- the normalization worked
# the central density, in nucleons/fm^3 -- compare nuclear saturation ~0.16
rho_central = rho[0]
show rho_central
show rho # the guess density: flat middle, soft edge Run it:
a_check: 16 (the normalization closes),
rho_central: 0.0987 nucleons/fm³ — compare nuclear saturation,
0.16 fm⁻³; the guess is deliberately too diffuse, and the SCF will fix it —
and a figure of the density with the flat middle and soft edge.