nucleons.jl

@@ -138,99 +138,29 @@ function findAllOrbitals(p::Bool, s::system, E_min=850.0, E_max=938.0)
     return (κs, Es)
 end
 
+"For a given list of κ values with corresponding energies, attempt to fill Z_or_N lowest lying orbitals and return the spectrum"
+function fillNucleons(Z_or_N::Int, κs, Es)::spectrum
+    sort_i = sortperm(Es)
+
+    occ = zeros(Int, length(κs))
+
+    for i in sort_i
+        if Z_or_N ≤ 0; break; end;
+        max_occ = 2 * j_κ(κs[i]) + 1
+        occ[i] = min(max_occ, Z_or_N)
+        Z_or_N -= occ[i]
+    end
+
+    @assert Z_or_N == 0 "All orbitals could not be filled"
+    return spectrum(κs, Es, occ)
+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"
@@ -241,7 +171,6 @@ function calculateNucleonDensity(p::Bool, s::system)::Tuple{Vector{Float64}, Vec
     ρ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
@@ -261,7 +190,7 @@ end
     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))
+    spec = fillNucleons(Z_or_N(s, p), κs, Es)
     if p
         s.p_spectrum = spec
     else
@@ -277,12 +206,5 @@ function solveNucleonDensity!(p::Bool, s::system, E_min=850.0, E_max=938.0)
     end
 end
 
-"Total energy of filled nucleons in the system, including the BCS pairing energy
+"Total energy of filled nucleons in the system"
-    E_pair = −Δ²/G for each species [Vretenar et al. 2005, Eq. 34 with constant-G pp-force]."
+nucleon_E(s::system) = sum(s.p_spectrum.occ .* (s.p_spectrum.E .- M_p)) + sum(s.n_spectrum.occ .* (s.n_spectrum.E .- M_n))
-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

system.jl

@@ -25,19 +25,15 @@ struct parameters
     end
 end
 
-"Tabulates a nucleon spectrum (protons or neutrons) containing κ, single-particle energies,
+"Tabulates a nucleon spectrum (protons or neutrons) containing κ and occupancy"
-    occupation numbers (2Ω·v² per κ, possibly fractional under BCS), the pairing gap Δ,
-    and the chemical potential λ."
 struct spectrum
     κ::Vector{Int}
     E::Vector{Float64}
-    occ::Vector{Float64}
+    occ::Vector{Int}
-    Δ::Float64
-    λ::Float64
 end
 
 "Initializes an unfilled spectrum"
-unfilled_spectrum() = spectrum(Int[], Float64[], Float64[], 0.0, 0.0)
+unfilled_spectrum() = spectrum(Int[], Float64[], Int[])
 
 "Defines a nuclear system containing relevant parameters and meson/nucleon densities"
 mutable struct system