289 lines
11 KiB
Julia
289 lines
11 KiB
Julia
using LinearAlgebra, DifferentialEquations, Interpolations
|
||
include("bisection.jl")
|
||
include("common.jl")
|
||
include("system.jl")
|
||
|
||
const M_n = 939.0 # MeV/c2
|
||
const M_p = 939.0 # MeV/c2
|
||
|
||
"The spherical Dirac equation that returns du=[dg, df] in-place where
|
||
u=[g, f] are the reduced radial components evaluated at r,
|
||
κ is the generalized angular momentum,
|
||
p is true for proton and false for neutron,
|
||
E in the energy in MeV,
|
||
f1(r) = M-Φ0(r)-W0(r)-(p-0.5)*2B0(r)-p*A0(r) is a function of r in MeV (see optimized_dirac_potentials()),
|
||
f2(r) = -M+Φ0(r)-W0(r)-(p-0.5)*2B0(r)-p*A0(r) is a function of r in MeV (see optimized_dirac_potentials()),
|
||
r is the radius in fm.
|
||
Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
|
||
function dirac!(du::Vector{Float64}, u::Vector{Float64}, (κ, E, f1, f2), r::Float64) # TODO: Static typing
|
||
(g, f) = u
|
||
@inbounds du[1] = -(κ/(r + r_reg)) * g + (E + f1(r)) * f / ħc
|
||
@inbounds du[2] = (κ/(r + r_reg)) * f - (E + f2(r)) * g / ħc
|
||
end
|
||
|
||
"Get the potentials f1 and f2 that goes into the Dirac equation, defined as
|
||
f1(r) = M-Φ0(r)-W0(r)-(p-0.5)*2B0(r)-p*A0(r),
|
||
f2(r) = -M+Φ0(r)-W0(r)-(p-0.5)*2B0(r)-p*A0(r)."
|
||
function optimized_dirac_potentials(p::Bool, s::system)
|
||
M = p ? M_p : M_n
|
||
|
||
f1s = zero_array(s)
|
||
f2s = zero_array(s)
|
||
|
||
@. f1s = M - s.Φ0 - s.W0 - (p - 0.5) * 2 * s.B0 - p * s.A0
|
||
@. f2s = -M + s.Φ0 - s.W0 - (p - 0.5) * 2 * s.B0 - p * s.A0
|
||
|
||
f1 = linear_interpolation(rs(s), f1s)
|
||
f2 = linear_interpolation(rs(s), f2s)
|
||
|
||
return (f1, f2)
|
||
end
|
||
|
||
"Approximate boundary condition for u(r)=[g(r), f(r)] at r -> ∞ where
|
||
κ is the generalized angular momentum,
|
||
p is true for proton and false for neutron,
|
||
E is the energy in MeV,
|
||
r is the radius in fm."
|
||
function asymp_BC(κ::Int, p::Bool, E::Float64, r::Float64)
|
||
M = p ? M_p : M_n
|
||
g = 1.0
|
||
f = ħc / (E + M) * (-√(M^2 - E^2) + κ/r) * g
|
||
return [g, f]
|
||
end
|
||
|
||
"Initial boundary condition for u(r)=[g(r), f(r)] at r=0"
|
||
init_BC() = [1.0, 1.0] # TODO: Why not [0.0, 1.0]?
|
||
|
||
"Solve the Dirac equation and return the wave function u(r)=[g(r), f(r)] where
|
||
matching_point marks the partitioning position (between 0.0 and 1.0) for shooting method,
|
||
the solution would be returned as a 2×(1+divs) matrix."
|
||
function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Vern9(), matching_point=0.15)
|
||
(f1, f2) = optimized_dirac_potentials(p, s)
|
||
|
||
# partitioning
|
||
mid_idx = s.divs * matching_point |> round |> Int
|
||
r_mid = rs(s)[mid_idx]
|
||
left_r = rs(s)[1:mid_idx]
|
||
right_r = rs(s)[mid_idx:end]
|
||
|
||
# left partition
|
||
prob = ODEProblem(dirac!, init_BC(), (0, r_mid))
|
||
sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=left_r)
|
||
wf_left = hcat(sol.u...)
|
||
|
||
# right partition
|
||
prob = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_max), (s.r_max, r_mid))
|
||
sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=right_r)
|
||
wf_right = reverse(hcat(sol.u...); dims=2)
|
||
|
||
# join two segments
|
||
u1 = wf_left[:, end]
|
||
u2 = wf_right[:, 1]
|
||
if norm(u2) < 1e-10
|
||
@warn "Right partition too small to rescale, setting to zero"
|
||
wf_right .= 0.0
|
||
else
|
||
proj = only(u1' * u2) / norm(u2)^2
|
||
wf_right .*= proj
|
||
end
|
||
wf = hcat(wf_left[:, 1:(end - 1)], wf_right)
|
||
|
||
if normalize
|
||
g2_int = simpsons_integrate(wf[1, :] .^ 2, Δr(s))
|
||
f2_int = simpsons_integrate(wf[2, :] .^ 2, Δr(s))
|
||
wf ./= sqrt(g2_int + f2_int)
|
||
end
|
||
|
||
return wf
|
||
end
|
||
|
||
"Returns a function that solves the Dirac equation in two partitions and returns
|
||
the determinant of [g_left(r) g_right(r); f_left(r) f_right(r)],
|
||
where is r is in fm."
|
||
function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_max/2, algo=Vern9())
|
||
(f1, f2) = optimized_dirac_potentials(p, s)
|
||
prob_left = ODEProblem(dirac!, init_BC(), (0, r))
|
||
function func(E)
|
||
prob_right = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_max), (s.r_max, r))
|
||
u_left = solve(prob_left, algo, p=(κ, E, f1, f2), saveat=[r])
|
||
u_right = solve(prob_right, algo, p=(κ, E, f1, f2), saveat=[r])
|
||
return u_left[1, 1] * u_right[2, 1] - u_right[1, 1] * u_left[2, 1]
|
||
end
|
||
return func
|
||
end
|
||
|
||
"Find all bound energies between E_min (=850.0) and E_max (=938.0) where
|
||
tol_digits determines the precision for root finding and the threshold for identifying duplicates,
|
||
the other parameters are the same from dirac!(...)."
|
||
function findEs(κ, p::Bool, s::system, E_min=850.0, E_max=938.0; tol_digits=8)
|
||
func = determinantFunc(κ, p, s)
|
||
Es = find_all_zeros(func, E_min, E_max; partitions=200, tol=1/10^tol_digits)
|
||
return unique(E -> round(E; digits=tol_digits), Es)
|
||
end
|
||
|
||
"Find all orbitals and return two lists containing κ values and corresponding energies for a single species where
|
||
the other parameters are defined above"
|
||
function findAllOrbitals(p::Bool, s::system, E_min=850.0, E_max=938.0)
|
||
κs = Int[]
|
||
Es = Float64[]
|
||
# start from κ=-1 and go both up and down
|
||
for direction in [-1, 1]
|
||
for κ in direction * (1:100) # cutoff is 100
|
||
new_Es = findEs(κ, p, s, E_min, E_max)
|
||
if isempty(new_Es); break; end
|
||
append!(Es, new_Es)
|
||
append!(κs, fill(κ, length(new_Es)))
|
||
end
|
||
end
|
||
return (κs, Es)
|
||
end
|
||
|
||
"Total angular momentum j for a given κ value"
|
||
j_κ(κ::Int) = abs(κ) - 1/2
|
||
|
||
"Orbital angular momentum l for a given κ value"
|
||
l_κ(κ::Int) = abs(κ) - (κ < 0) # since true = 1 and false = 0
|
||
|
||
"Pair degeneracy Ω = (2j+1)/2 for a given κ (equals |κ|)"
|
||
Ω_κ(κ::Int) = abs(κ)
|
||
|
||
"Constant-G monopole pairing strength: G = 20/A MeV [Bohr-Mottelson Vol. 1]"
|
||
pairing_G(s::system)::Float64 = 20.0 / A(s)
|
||
|
||
"Energy window (MeV) around the Fermi level within which BCS pairing is active;
|
||
orbitals below (above) the window are treated as fully filled (empty) core states.
|
||
This cut-off regularizes the otherwise ill-defined constant-G gap equation."
|
||
const E_pair_window = 5.0
|
||
|
||
"Sharp (non-pairing) filling of the Z_or_N lowest orbitals — kept for diagnostics."
|
||
function fillNucleons(Z_or_N::Int, κs, Es)::spectrum
|
||
occ = zeros(Float64, length(κs))
|
||
rem = Z_or_N
|
||
for i in sortperm(Es)
|
||
rem ≤ 0 && break
|
||
take = min(2 * Ω_κ(κs[i]), rem)
|
||
occ[i] = take
|
||
rem -= take
|
||
end
|
||
@assert rem == 0 "All orbitals could not be filled"
|
||
return spectrum(κs, Es, occ, 0.0, 0.0)
|
||
end
|
||
|
||
"Solve BCS equations for a constant-G monopole pairing force. Returns (λ, Δ, occ) where
|
||
occ[i] = 2 Ω_i v²_i sums to N_target, and the canonical-basis occupations are
|
||
v²_i = ½(1 - (ε_i − λ)/√((ε_i − λ)² + Δ²)) [Vretenar et al. 2005, Eq. 50].
|
||
Orbitals outside ±E_window of the Fermi level are frozen as fully filled / empty."
|
||
function solveBCS(κs::Vector{Int}, Es::Vector{Float64}, G::Float64, N_target::Int;
|
||
E_window::Float64=E_pair_window, tol::Float64=1e-10, max_iter::Int=200)
|
||
isempty(κs) && return (0.0, 0.0, Float64[])
|
||
|
||
Ω = Float64.(Ω_κ.(κs))
|
||
occ = zeros(Float64, length(κs))
|
||
|
||
# Sharp-filling estimate of the Fermi level
|
||
idx = sortperm(Es)
|
||
cum = cumsum(2 .* Ω[idx])
|
||
i_HOMO = findfirst(≥(N_target), cum)
|
||
i_HOMO === nothing && error("Not enough orbitals for N_target=$N_target particles")
|
||
λ_F = Es[idx[i_HOMO]]
|
||
|
||
# Core (frozen filled), active (paired), and above-window (empty)
|
||
core = findall(e -> e < λ_F - E_window, Es)
|
||
active = findall(e -> abs(e - λ_F) ≤ E_window, Es)
|
||
occ[core] .= 2 .* Ω[core]
|
||
N_active = N_target - sum(occ[core])
|
||
max_active = 2 * sum(Ω[active])
|
||
|
||
# Closed / trivial case: no pairing to solve in the active window
|
||
if N_active ≤ 0 || N_active ≥ max_active
|
||
rem = N_active
|
||
for i in sort(active, by = i -> Es[i])
|
||
take = clamp(rem, 0.0, 2 * Ω[i])
|
||
occ[i] = take
|
||
rem -= take
|
||
end
|
||
return (λ_F, 0.0, occ)
|
||
end
|
||
|
||
Es_a, Ω_a = Es[active], Ω[active]
|
||
v²(λ, Δ) = @. 0.5 * (1 - (Es_a - λ) / sqrt((Es_a - λ)^2 + Δ^2))
|
||
N_of(λ, Δ) = sum(2 .* Ω_a .* v²(λ, Δ))
|
||
gap_f(λ, Δ) = G/2 * sum(@. Ω_a / sqrt((Es_a - λ)^2 + Δ^2)) - 1
|
||
|
||
λ, Δ = λ_F, 1.0
|
||
λ_lo, λ_hi = minimum(Es_a) - 10.0, maximum(Es_a) + 10.0
|
||
for _ in 1:max_iter
|
||
λ_new = bisection(λ -> N_of(λ, Δ) - N_active, λ_lo, λ_hi; tol=tol)
|
||
Δ_new = gap_f(λ_new, 0.0) < 0 ? 0.0 :
|
||
bisection(Δ -> gap_f(λ_new, Δ), 0.0, 50.0; tol=tol)
|
||
converged = abs(λ_new - λ) < tol && abs(Δ_new - Δ) < tol
|
||
λ, Δ = λ_new, Δ_new
|
||
converged && break
|
||
end
|
||
|
||
occ[active] .= 2 .* Ω_a .* v²(λ, Δ)
|
||
return (λ, Δ, occ)
|
||
end
|
||
|
||
"Fill orbitals via BCS with constant-G monopole pairing and return the spectrum."
|
||
function fillNucleonsBCS(N_target::Int, κs, Es, G::Float64)::spectrum
|
||
(λ, Δ, occ) = solveBCS(κs, Es, G, N_target)
|
||
return spectrum(κs, Es, occ, Δ, λ)
|
||
end
|
||
|
||
"Calculate scalar and vector densities of a nucleon species on [0,r_max] divided into (divs+1) points and returns them as vectors (ρ_s, ρ_v) where
|
||
the arrays κs, Es, occs tabulate the energies and occupation numbers corresponding to each κ,
|
||
the other parameters are defined above"
|
||
function calculateNucleonDensity(p::Bool, s::system)::Tuple{Vector{Float64}, Vector{Float64}}
|
||
spectrum = p ? s.p_spectrum : s.n_spectrum
|
||
(κs, Es, occs) = (spectrum.κ, spectrum.E, spectrum.occ)
|
||
|
||
ρr2 = zeros(2, s.divs + 1) # ρ×r² values
|
||
|
||
for (κ, E, occ) in zip(κs, Es, occs)
|
||
occ < 1e-10 && continue
|
||
wf = solveNucleonWf(κ, p, E, s; normalize=true)
|
||
wf2 = wf .* wf
|
||
ρr2 += (occ / (4 * pi)) * wf2 # 2j+1 factor is accounted in the occupancy number
|
||
end
|
||
|
||
r2s = rs(s).^2
|
||
ρ = ρr2 ./ transpose(r2s)
|
||
ρ[:, 1] .= ρ[:, 2] # dirty fix for NaN at r=0
|
||
|
||
ρ_s = ρ[1, :] - ρ[2, :] # g^2 - f^2
|
||
ρ_v = ρ[1, :] + ρ[2, :] # g^2 + f^2
|
||
|
||
return (ρ_s, ρ_v)
|
||
end
|
||
|
||
"For a nucleon species, solve the Dirac equation and save the spectrum & densities in-place where
|
||
the other parameters are defined above"
|
||
function solveNucleonDensity!(p::Bool, s::system, E_min=850.0, E_max=938.0)
|
||
κs, Es = findAllOrbitals(p, s, E_min, E_max)
|
||
spec = fillNucleonsBCS(Z_or_N(s, p), κs, Es, pairing_G(s))
|
||
if p
|
||
s.p_spectrum = spec
|
||
else
|
||
s.n_spectrum = spec
|
||
end
|
||
(ρ_s, ρ_v) = calculateNucleonDensity(p, s)
|
||
if p
|
||
s.ρ_sp = ρ_s
|
||
s.ρ_vp = ρ_v
|
||
else
|
||
s.ρ_sn = ρ_s
|
||
s.ρ_vn = ρ_v
|
||
end
|
||
end
|
||
|
||
"Total energy of filled nucleons in the system, including the BCS pairing energy
|
||
E_pair = −Δ²/G for each species [Vretenar et al. 2005, Eq. 34 with constant-G pp-force]."
|
||
function nucleon_E(s::system)
|
||
G = pairing_G(s)
|
||
E_sp = sum(s.p_spectrum.occ .* (s.p_spectrum.E .- M_p)) +
|
||
sum(s.n_spectrum.occ .* (s.n_spectrum.E .- M_n))
|
||
E_pair = -(s.p_spectrum.Δ^2 + s.n_spectrum.Δ^2) / G
|
||
return E_sp + E_pair
|
||
end
|