Model parameters are no longer hard-coded

NuclearRMF.jl

@@ -0,0 +1,68 @@
+using PolyLog, Plots
+include("nucleons.jl")
+include("mesons.jl")
+
+"Total binding energy of the system"
+function total_E(s::system)
+    E_cm = (3/4) * 41.0 * A(s)^(-1/3) # Center-of-mass correction [Ring and Schuck, Sec. 2.3]
+    return -(nucleon_E(s) + meson_E(s)) + E_cm
+end
+
+"Normalized Woods-Saxon form used for constructing an initial solution"
+Woods_Saxon(r::Float64; R::Float64=7.0, a::Float64=0.5) = -1 / (8π * a^3 * reli3(-exp(R / a)) * (1 + exp((r - R) / a)))
+
+"Conditionally toggle @time if enable_time is true"
+macro conditional_time(label, expr)
+    quote
+        if $(esc(:enable_time))
+            @time $label $(esc(expr))
+        else
+            $(esc(expr))
+        end
+    end
+end
+
+"Run the full program by self-consistent solution of nucleon and meson densities"
+function solve_system!(s::system; reinitialize_densities=true, print_E=true, print_time=false, live_plots=false)
+    if reinitialize_densities
+        dens_guess = Woods_Saxon.(rs(s))
+        @. s.ρ_sp = s.Z * dens_guess
+        @. s.ρ_vp = s.Z * dens_guess
+        @. s.ρ_sn = s.N * dens_guess
+        @. s.ρ_vn = s.N * dens_guess
+    end
+
+    if live_plots
+        p = plot(legends=false, size=(1024, 768), layout=(2, 4), title=["ρₛₚ" "ρᵥₚ" "ρₛₙ" "ρᵥₙ" "Φ₀" "W₀" "B₀" "A₀"])
+    end
+
+    previous_E_per_A = NaN
+
+    while true
+        enable_time = print_time
+        @conditional_time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 15)
+        @conditional_time "Proton densities" solveNucleonDensity!(true, s)
+        @conditional_time "Neutron densities" solveNucleonDensity!(false, s)
+
+        if live_plots
+            for s in p.series_list
+                s.plotattributes[:linecolor] = :gray
+            end
+            plot!(p, rs(s), hcat(s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn, s.Φ0, s.W0, s.B0, s.A0), linecolor=:red)
+            display(p)
+        end
+
+        E_per_A = total_E(s) / A(s)
+        print_E && println("Total binding E per nucleon = $E_per_A")
+
+        # check convergence
+        abs(previous_E_per_A - E_per_A) < 0.0001 && break
+        previous_E_per_A = E_per_A
+    end
+end
+
+"Calculate RMS radius from density"
+rms_radius(p::Bool, s::system) = 4pi * Δr(s) * sum((rs(s) .^ 4) .* (p ? s.ρ_vp : s.ρ_vn)) / (p ? s.Z : s.N) |> sqrt
+
+"Calculate neutron skin thickness"
+R_skin(s::system) = rms_radius(false, s) - rms_radius(true, s)

Project.toml

@@ -3,4 +3,3 @@ DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PolyLog = "85e3b03c-9856-11eb-0374-4dc1f8670e7f"
-Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"

bisection.jl

@@ -0,0 +1,62 @@
+"Bisection method: Finds a zero in the interval [a, b]"
+function bisection(f::Function, a::Float64, b::Float64; tol::Float64 = 1e-8, max_iter::Int = 1000)::Float64
+    fa::Float64 = f(a)
+    fb::Float64 = f(b)
+    
+    # Ensure the endpoints bracket a root.
+    if fa * fb > 0.0
+        error("No sign change detected in the interval: f(a)*f(b) > 0")
+    end
+
+    c::Float64 = a
+    for _ in 1:max_iter
+        c = (a + b) / 2.0
+        fc::Float64 = f(c)
+        
+        # Check for convergence using the function value or the interval width.
+        if abs(fc) < tol || abs(b - a) < tol
+            return c
+        elseif fa * fc < 0.0
+            b = c
+            fb = fc
+        else
+            a = c
+            fa = fc
+        end
+    end
+    return c  # Return the current approximation if max iterations reached.
+end
+
+"Function to find all zeros in [x_start, x_end] by scanning for sign changes."
+function find_all_zeros(f::Function, x_start::Float64, x_end::Float64; partitions::Int = 1000, tol::Float64 = 1e-8)::Vector{Float64}
+    zeros_list::Vector{Float64} = Float64[]
+    Δ::Float64 = (x_end - x_start) / partitions
+    x_prev::Float64 = x_start
+    f_prev::Float64 = f(x_prev)
+    
+    for i in 1:partitions
+        x_curr::Float64 = x_start + i * Δ
+        f_curr::Float64 = f(x_curr)
+        
+        # Check for a sign change or an exact zero.
+        if f_prev * f_curr < 0.0 || f_prev == 0.0
+            try
+                root::Float64 = bisection(f, x_prev, x_curr; tol=tol)
+                # Add the root if it's new (avoid duplicates from neighboring intervals)
+                if isempty(zeros_list) || abs(root - zeros_list[end]) > tol
+                    push!(zeros_list, root)
+                end
+            catch err
+                # Skip this interval if bisection fails (should not occur if a sign change exists)
+            end
+        elseif f_curr == 0.0
+            if isempty(zeros_list) || abs(x_curr - zeros_list[end]) > tol
+                push!(zeros_list, x_curr)
+            end
+        end
+        
+        x_prev = x_curr
+        f_prev = f_curr
+    end
+    return zeros_list
+end

common.jl

@@ -0,0 +1,17 @@
+const ħc = 197.33 # MeVfm
+const r_reg = 1E-8 # fm # regulator for R
+
+"Integrate a uniformly discretized function f using Simpson's rule where h is the step size and coefficient is an optional scaling factor."
+function simpsons_integrate(f::AbstractVector{Float64}, h::Float64; coefficient::Float64 = 1.0)
+    @assert length(f) % 2 == 1 "Number of mesh divisions must be even for Simpson's rule"
+    s = sum(enumerate(f)) do (i, fi)
+        if i == 1 || i == length(f)
+            return fi
+        elseif i % 2 == 0
+            return 4fi
+        else
+            return 2fi
+        end
+    end
+    return (h/3) * coefficient * s
+end

mesons.jl

@@ -1,23 +1,5 @@
-using DifferentialEquations
+include("common.jl")
-
+include("system.jl")
-const ħc = 197.33 # MeVfm
-
-# Values defined in C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001)
-# Values taken from Hartree.f (FSUGarnet)
-const m_s = 496.939473213388 # MeV/c2
-const m_ω = 782.5 # MeV/c2
-const m_ρ = 763.0 # MeV/c2
-const m_γ = 0.000001000 # MeV/c2 # should be 0?
-const g2_s = 110.349189097820 # dimensionless
-const g2_v = 187.694676506801 # dimensionless
-const g2_ρ = 192.927428365698 # dimensionless
-const g2_γ = 0.091701236 # dimensionless # equal to 4πα
-const κ_ss = 3.260178893447 # MeV
-const λ = -0.003551486718 # dimensionless # LambdaSS
-const ζ = 0.023499504053 # dimensionless # LambdaVV
-const Λv = 0.043376933644 # dimensionless # LambdaVR
-
-const r_reg = 1E-8 # fm # regulator for Green's functions
 
 "Green's function for Klein-Gordon equation in natural units"
 greensFunctionKG(m, r, rp) = -1 / (m * (r + r_reg) * (rp + r_reg)) * sinh(m * min(r, rp)) * exp(-m * max(r, rp))
@@ -28,7 +10,10 @@ greensFunctionP(r, rp) = -1 / (max(r, rp) + r_reg)
 "Solve the Klein-Gordon equation (or Poisson's equation when m=0) and return an array in MeV for a source array given in fm⁻³ where
     m is the mass of the meson in MeV/c2."
 function solveKG(m, source, s::system)
-    int_measure = ħc .* Δr(s) .* rs(s) .^ 2
+
+    @assert s.divs % 2 == 0 "Number of mesh divisions must be even for Simpson's rule"
+    simpsons_weights = (Δr(s)/3) .* [1; repeat([2,4], s.divs ÷ 2)[2:end]; 1]
+    int_measure = ħc .* simpsons_weights .* rs(s) .^ 2
 
     greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m / ħc, r, rp)
     greenMat = greensFunction.(rs(s), transpose(rs(s)))
@@ -41,7 +26,8 @@ end
     r is the radius in fm,
     the inital solutions are read from s and the final solutions are saved in-place.
 Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
-function solveMesonFields!(s::system, iterations=10)
+function solveMesonFields!(s::system, iterations=15, oscillation_control_parameter=0.3)
+    (m_s, m_ω, m_ρ, m_γ, g2_s, g2_v, g2_ρ, g2_γ, κ_ss, λ, ζ, Λv) = (s.param.m_s, s.param.m_ω, s.param.m_ρ, s.param.m_γ, s.param.g2_s, s.param.g2_v, s.param.g2_ρ, s.param.g2_γ, s.param.κ_ss, s.param.λ, s.param.ζ, s.param.Λv)
     (Φ0, W0, B0, A0) = (s.Φ0, s.W0, s.B0, s.A0)
     (ρ_sp, ρ_vp, ρ_sn, ρ_vn) = (s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)
 
@@ -55,10 +41,32 @@ function solveMesonFields!(s::system, iterations=10)
         @. src_Φ0 = g2_s * ((κ_ss/ħc)/2 * (Φ0/ħc)^2 + (λ/6) * (Φ0/ħc)^3) - g2_s * (ρ_sp + ρ_sn)
         @. src_W0 = g2_v * ((ζ/6) * (W0/ħc)^3 + 2Λv * (2B0/ħc)^2 * (W0/ħc)) - g2_v * (ρ_vp + ρ_vn)
         @. src_B0 = (2Λv * g2_ρ * (W0/ħc)^2 * (2B0/ħc) - g2_ρ/2 * (ρ_vp - ρ_vn)) / 2
+
         Φ0 .= solveKG(m_s, src_Φ0, s)
-        W0 .= (solveKG(m_ω, src_W0, s) .+ W0) ./ 2 # to supress oscillation
+
-        B0 .= (solveKG(m_ρ, src_B0, s) .+ B0) ./ 2 # to supress oscillation
+        # W0 and B0 keep a fraction of their previous solutions to suppress oscillations
+        W0 .*= (1 - oscillation_control_parameter)
+        B0 .*= (1 - oscillation_control_parameter)
+
+        W0 .+= solveKG(m_ω, src_W0, s) .* oscillation_control_parameter
+        B0 .+= solveKG(m_ρ, src_B0, s) .* oscillation_control_parameter
     end
 
     return (Φ0, W0, B0, A0)
 end
+
+"Calculate the total energy associated with meson fields"
+function meson_E(s::system)
+    (κ_ss, λ, ζ, Λv) = (s.param.κ_ss, s.param.λ, s.param.ζ, s.param.Λv)
+
+    E_densities = map(zip(s.Φ0, s.W0, s.B0, s.A0, s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)) do (Φ0, W0, B0, A0, ρ_sp, ρ_vp, ρ_sn, ρ_vn)
+        E_σ = (1/2) * (Φ0/ħc) * (ρ_sp + ρ_sn) - ((κ_ss/ħc)/12 * (Φ0/ħc)^3 + (λ/24) * (Φ0/ħc)^4)
+        E_ω = -(1/2) * (W0/ħc) * (ρ_vp + ρ_vn) + (ζ/24) * (W0/ħc)^4
+        E_ρ = -(1/4) * (2B0/ħc) * (ρ_vp - ρ_vn)
+        E_γ = -(1/2) * (A0/ħc) * ρ_vp
+        E_ωρ = Λv * (W0/ħc)^2 * (2B0/ħc)^2
+        E_σ + E_ω + E_ρ + E_γ + E_ωρ
+    end
+    
+    return simpsons_integrate(E_densities .* rs(s).^2, Δr(s); coefficient=4π * ħc)
+end

nucleons.jl

@@ -1,88 +1,129 @@
-using DifferentialEquations, Roots
+using LinearAlgebra, DifferentialEquations, Interpolations
+include("bisection.jl")
+include("common.jl")
+include("system.jl")
 
-const ħc = 197.33 # MeVfm
 const M_n = 939.0 # MeV/c2
 const M_p = 939.0 # MeV/c2
 
-const r_reg = 1E-8 # fm # regulator for the centrifugal term
-
 "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,
-    Φ0, W0, B0, A0 are the mean-field potentials (couplings included) in MeV as functions of r in fm,
+    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, u, (κ, p, E, Φ0, W0, B0, A0), r)
+function dirac!(du::Vector{Float64}, u::Vector{Float64}, (κ, E, f1, f2), r::Float64) # TODO: Static typing
-    M = p ? M_p : M_n
-    common1 = E - W0(r) - (p - 0.5) * 2B0(r) - p * A0(r)
-    common2 = M - Φ0(r)
     (g, f) = u
-    du[1] = -(κ/(r + r_reg)) * g + (common1 + common2) * f / ħc
+    @inbounds du[1] = -(κ/(r + r_reg)) * g + (E + f1(r)) * f / ħc
-    du[2] =  (κ/(r + r_reg)) * f - (common1 - common2) * g / ħ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
     divs is the number of mesh divisions so solution would be returned as a 2×(1+divs) matrix,
-    shooting method divides the interval into two partitions at r_max/2, ensuring convergence at both r=0 and r=r_max,
     the other parameters are the same from dirac!(...)."
-function solveNucleonWf(κ, p, E, s::system; shooting=true, normalize=true)
+function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Vern9())
-    (Φ0, W0, B0, A0) = fields_interp(s)
+    (f1, f2) = optimized_dirac_potentials(p, s)
 
-    if shooting
+    # partitioning
-        @assert s.divs % 2 == 0 "divs must be an even number when shooting=true"
+    mid_idx = s.divs ÷ 2
-        prob = ODEProblem(dirac!, [0, 1], (s.r_max, s.r_max / 2))
+    r_mid = rs(s)[mid_idx]
-        sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr(s))
+    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)    
-        next_r_max = s.r_max / 2 # for the next step
+    
+    # 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
-        next_r_max = s.r_max
+        proj = only(u1' * u2) / norm(u2)^2
-    end
+        wf_right .*= proj
-    
-    prob = ODEProblem(dirac!, [0, 1], (0, next_r_max))
-    sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr(s))
-    wf = hcat(sol.u...)
-
-    if shooting # join two segments
-        rescale_factor_g = wf[1, end] / wf_right[1, 1]
-        rescale_factor_f = wf[2, end] / wf_right[2, 1]
-        @assert isfinite(rescale_factor_g) && isfinite(rescale_factor_f) "Cannot rescale the right partition"
-        isapprox(rescale_factor_g, rescale_factor_f; rtol=0.03) || @warn "Discontinuity between two partitions"
-        wf_right_rescaled = wf_right .* ((rescale_factor_g + rescale_factor_f) / 2)
-        wf = hcat(wf[:, 1:(end - 1)], wf_right_rescaled)
     end
+    wf = hcat(wf_left[:, 1:(end - 1)], wf_right)
 
     if normalize
-        norm = sum(wf .* wf) * Δr(s) # integration by Reimann sum
+        g2_int = simpsons_integrate(wf[1, :] .^ 2, Δr(s))
-        wf = wf ./ sqrt(norm)
+        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 and returns g(r=r_max) where
+"Returns a function that solves the Dirac equation in two partitions and returns 
-    r_max is the outer boundary in fm,
+    the determinant of [g_left(r) g_right(r); f_left(r) f_right(r)],
-    the other parameters are the same from dirac!(...)."
+    where is r is in fm."
-function boundaryValueFunc(κ, p, s::system; dtype=Float64, algo=Tsit5())
+function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_max/2, algo=Vern9())
-    prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, s.r_max))
+    (f1, f2) = optimized_dirac_potentials(p, s)
-    (Φ0, W0, B0, A0) = fields_interp(s)
+    prob_left  = ODEProblem(dirac!, init_BC(), (0, r))
-    func(E) = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=[s.r_max], save_idxs=[1])[1, 1]
+    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 (=0) and E_max (=mass) where
+"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, s::system, E_min=0, E_max=(p ? M_p : M_n), tol_digits=5)
+function findEs(κ, p::Bool, s::system, E_min=850.0, E_max=938.0; tol_digits=8)
-    func = boundaryValueFunc(κ, p, s)
+    func = determinantFunc(κ, p, s)
-    Es = find_zeros(func, (E_min, E_max); xatol=1/10^tol_digits)
+    Es = find_all_zeros(func, E_min, E_max; partitions=20, 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, s::system, E_min=0, E_max=(p ? M_p : M_n))
+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
@@ -97,8 +138,8 @@ function findAllOrbitals(p, s::system, E_min=0, E_max=(p ? M_p : M_n))
     return (κs, Es)
 end
 
-"For a given list of κ values with corresponding energies, attempt to fill Z_or_N lowest lying orbitals and return occupancy numbers"
+"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)
+function fillNucleons(Z_or_N::Int, κs, Es)::spectrum
     sort_i = sortperm(Es)
 
     occ = zeros(Int, length(κs))
@@ -110,8 +151,8 @@ function fillNucleons(Z_or_N::Int, κs, Es)
         Z_or_N -= occ[i]
     end
 
-    Z_or_N == 0 || @warn "All orbitals could not be filled"
+    @assert Z_or_N == 0 "All orbitals could not be filled"
-    return occ
+    return spectrum(κs, Es, occ)
 end
 
 "Total angular momentum j for a given κ value"
@@ -123,11 +164,14 @@ l_κ(κ::Int) = abs(κ) - (κ < 0) # since true = 1 and false = 0
 "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(κs, Es, occs, p, s::system)
+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)
-        wf = solveNucleonWf(κ, p, E, s; shooting=true, normalize=true)
+        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
@@ -142,10 +186,25 @@ function calculateNucleonDensity(κs, Es, occs, p, s::system)
     return (ρ_s, ρ_v)
 end
 
-"Solve the Dirac equation and calculate scalar and vector densities of a nucleon species where
+"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, s::system, E_min=800, E_max=939)
+function solveNucleonDensity!(p::Bool, s::system, E_min=850.0, E_max=938.0)
     κs, Es = findAllOrbitals(p, s, E_min, E_max)
-    occs = fillNucleons(Z_or_N(s, p), κs, Es)
+    spec = fillNucleons(Z_or_N(s, p), κs, Es)
-    return calculateNucleonDensity(κs, Es, occs, p, 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"
+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))

system.jl

@@ -1,13 +1,52 @@
-using Interpolations, PolyLog, Plots
+"Stores a set of coupling constants and meson masses that goes into the Lagrangian"
+struct parameters
+    # Values defined in C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001)
+    m_s::Float64 # MeV/c2
+    m_ω::Float64 # MeV/c2
+    m_ρ::Float64 # MeV/c2
+    m_γ::Float64 # MeV/c2
+    g2_s::Float64 # dimensionless
+    g2_v::Float64 # dimensionless
+    g2_ρ::Float64 # dimensionless
+    g2_γ::Float64 # dimensionless
+    κ_ss::Float64 # MeV
+    λ::Float64 # dimensionless, aka LambdaSS
+    ζ::Float64 # dimensionless, aka LambdaVV
+    Λv::Float64 # dimensionless, aka LambdaVR
+
+    "Dummy struct when parameters are not needed (for testing)"
+    parameters() = new(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
+
+    "Initialize parameters from a string with values provided in order of struct definition separated by commas"
+    function parameters(s::String)
+        values = parse.(Float64, strip.(split(s, ',')))
+        @assert length(values) == 12 "String must contain exactly 12 values separated by commas"
+        return new(values...)
+    end
+end
+
+"Tabulates a nucleon spectrum (protons or neutrons) containing κ and occupancy"
+struct spectrum
+    κ::Vector{Int}
+    E::Vector{Float64}
+    occ::Vector{Int}
+end
+
+"Initializes an unfilled spectrum"
+unfilled_spectrum() = spectrum(Int[], Float64[], Int[])
 
 "Defines a nuclear system containing relevant parameters and meson/nucleon densities"
 mutable struct system
     Z::Int
     N::Int
 
+    param::parameters
     r_max::Float64
     divs::Int
 
+    p_spectrum::spectrum
+    n_spectrum::spectrum
+
     Φ0::Vector{Float64}
     W0::Vector{Float64}
     B0::Vector{Float64}
@@ -19,10 +58,10 @@ mutable struct system
     ρ_vn::Vector{Float64}
 
     "Initialize an unsolved system"
-    system(Z, N, r_max, divs) = new(Z, N, r_max, divs, [zeros(1 + divs) for _ in 1:8]...)
+    system(Z, N, parameters, r_max, divs) = new(Z, N, parameters, r_max, divs, unfilled_spectrum(), unfilled_spectrum(), [zeros(1 + divs) for _ in 1:8]...)
 
-    "Dummy struct to define the mesh"
+    "Dummy struct to define the mesh (for testing)"
-    system(r_max, divs) = system(0, 0, r_max, divs)
+    system(r_max, divs) = system(0, 0, parameters(), r_max, divs)
 end
 
 "Get mass number of nucleus"
@@ -39,66 +78,3 @@ Z_or_N(s::system, p::Bool)::Int = p ? s.Z : s.N
 
 "Create an empty array for the size of the mesh"
 zero_array(s::system) = zeros(1 + s.divs)
-
-"Normalized Woods-Saxon form used for constructing an initial solution"
-Woods_Saxon(r::Float64; R::Float64=7.0, a::Float64=0.5) = -1 / (8π * a^3 * reli3(-exp(R / a)) * (1 + exp((r - R) / a)))
-
-"Create linear interpolations for the meson fields"
-function fields_interp(s::system)
-    S_interp = linear_interpolation(rs(s), s.Φ0)
-    V_interp = linear_interpolation(rs(s), s.W0)
-    R_interp = linear_interpolation(rs(s), s.B0)
-    A_interp = linear_interpolation(rs(s), s.A0)
-    return (S_interp, V_interp, R_interp, A_interp)
-end
-
-include("nucleons.jl")
-include("mesons.jl")
-
-"Run the full program by self-consistent solution of nucleon and meson densities"
-function solve_system(s::system; reinitialize_densities=true, monitor_print=true, monitor_plot=false)
-    if reinitialize_densities
-        dens_guess = Woods_Saxon.(rs(s))
-        @. s.ρ_sp = s.Z * dens_guess
-        @. s.ρ_vp = s.Z * dens_guess
-        @. s.ρ_sn = s.N * dens_guess
-        @. s.ρ_vn = s.N * dens_guess
-    end
-
-    if monitor_plot
-        p = plot(legends=false, size=(1024, 768), layout=(2, 4), title=["ρₛₚ" "ρᵥₚ" "ρₛₙ" "ρᵥₙ" "Φ₀" "W₀" "B₀" "A₀"])
-    end
-
-    E_total_previous = NaN
-
-    while true
-        # mesons
-        @time "Meson fields" solveMesonFields!(s, isnan(E_total_previous) ? 50 : 5)
-
-        # protons
-        @time "Proton spectrum" (κs_p, Es_p) = findAllOrbitals(true, s)
-        occs_p = fillNucleons(s.Z, κs_p, Es_p)
-        @time "Proton densities" (s.ρ_sp, s.ρ_vp) = calculateNucleonDensity(κs_p, Es_p, occs_p, true, s)
-
-        # neutrons
-        @time "Neutron spectrum" (κs_n, Es_n) = findAllOrbitals(false, s)
-        occs_n = fillNucleons(s.N, κs_n, Es_n)
-        @time "Neutron densities" (s.ρ_sn, s.ρ_vn) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, s)
-
-        if monitor_plot
-            for s in p.series_list
-                s.plotattributes[:linecolor] = :gray
-            end
-            plot!(p, rs(s), hcat(s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn, s.Φ0, s.W0, s.B0, s.A0), linecolor=:red)
-            display(p)
-        end
-
-        # total binding energy of nucleons
-        E_total = sum(M_p .- Es_p) + sum(M_n .- Es_n)
-        monitor_print && println("Total binding E per nucleon = $(E_total/A(s))")
-
-        # check convergence
-        abs(E_total_previous - E_total) < 0.1 && break
-        E_total_previous = E_total
-    end
-end

test/Linear_nucleon_spectrum.jl

@@ -1,5 +1,5 @@
-using DelimitedFiles, Interpolations, Plots
+using DelimitedFiles, Plots
-include("../system.jl")
+include("../nucleons.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -13,19 +13,16 @@ As = test_data[:, 5]
 p = false
 r_max = maximum(xs)
 divs = length(xs) - 1
-s = system(8, 8, r_max, divs)
+s = system(8, 8, parameters(), r_max, divs)
 
 s.Φ0 = Ss
 s.W0 = Vs
 s.B0 = Rs
 s.A0 = As
 
-E_min = 860
+κs, Es = findAllOrbitals(p, s)
-E_max = 939
-
-κs, Es = findAllOrbitals(p, s, E_min, E_max)
 Ebinds = (p ? M_p : M_n) .- Es
-occ = fillNucleons(Z_or_N(s, p), κs, Es)
+spec = fillNucleons(Z_or_N(s, p), κs, Es)
 
 # format: proton, kappa, filling, gnodes, fnodes, Ebind
 bench_data, _ = readdlm("test/LinearSpectrum.csv", ','; header=true)
@@ -36,7 +33,7 @@ bench_filling = bench_data[mask, 3]
 bench_occ = Int.(bench_filling .* (2 .* j_κ.(Int.(bench_κs)) .+ 1))
 
 scatter(κs, Ebinds, label=(p ? "proton" : "neutron") * " spectrum")
-annotate!(κs .+ 0.3, Ebinds .+ 1, [(o, 9) for o in occ])
+annotate!(κs .+ 0.3, Ebinds .+ 1, [(o, 9) for o in spec.occ])
 
 scatter!(bench_κs, bench_Ebinds, label="benchmark"; markershape=:cross, markercolor="red")
 annotate!(bench_κs .+ 0.3, bench_Ebinds .- 1, [(o, 9, :red) for o in bench_occ])

test/Pb208.jl

@@ -1,5 +1,14 @@
-include("../system.jl")
+include("../NuclearRMF.jl")
 
-s = system(82, 126, 20.0, 400)
+Z = 82
+N = 126
 
-solve_system(s; monitor_plot=true)
+# Parameter values calibrated by FSUGarnet for Pb-208
+param = parameters("496.939473213388, 782.5, 763.0, 0.000001000, 110.349189097820, 187.694676506801, 192.927428365698, 0.091701236, 3.260178893447, -0.003551486718, 0.023499504053, 0.043376933644")
+
+r_max = 20.0
+divs = 400
+
+s = system(Z, N, param, r_max, divs)
+
+solve_system!(s; live_plots=false)

test/Pb208ParamsFSUGarnet.txt

@@ -0,0 +1 @@
+496.939473213388, 782.5, 763.0, 0.000001000, 110.349189097820, 187.694676506801, 192.927428365698, 0.091701236, 3.260178893447, -0.003551486718, 0.023499504053, 0.043376933644

test/Pb208_meson_flds.jl

@@ -1,5 +1,5 @@
-using DelimitedFiles, Interpolations, Plots
+using DelimitedFiles, Plots
-include("../system.jl")
+include("../mesons.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -17,9 +17,12 @@ plot(xs_bench, hcat(Φ0_bench, W0_bench, B0_bench, A0_bench), layout=4, label=["
 test_data = readdlm("test/Pb208DensFSUGarnet.csv")
 xs = test_data[:, 1]
 
+N_p = 82
+N_n = 126
+param = parameters(read("test/Pb208ParamsFSUGarnet.txt", String))
 r_max = maximum(xs)
 divs = length(xs) - 1
-s = system(r_max, divs)
+s = system(N_p, N_n, param, r_max, divs)
 
 s.ρ_sn = test_data[:, 2]
 s.ρ_vn = test_data[:, 3]

test/Pb208_nucleon_dens.jl

@@ -1,5 +1,5 @@
 using DelimitedFiles, Plots
-include("../system.jl")
+include("../nucleons.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -14,22 +14,22 @@ N_p = 82
 N_n = 126
 r_max = maximum(xs)
 divs = length(xs) - 1
-s = system(N_p, N_n, r_max, divs)
+s = system(N_p, N_n, parameters(), r_max, divs)
 
 s.Φ0 = Ss
 s.W0 = Vs
 s.B0 = Rs
 s.A0 = As
 
-(ρ_sp, ρ_vp) = solveNucleonDensity(true, s)
+solveNucleonDensity!(true, s)
 
-p_sp = plot(rs(s), ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
+p_sp = plot(rs(s), s.ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
-p_vp = plot(rs(s), ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
+p_vp = plot(rs(s), s.ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
 
-(ρ_sn, ρ_vn) = solveNucleonDensity(false, s)
+solveNucleonDensity!(false, s)
 
-p_sn = plot(rs(s), ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
+p_sn = plot(rs(s), s.ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
-p_vn = plot(rs(s), ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")
+p_vn = plot(rs(s), s.ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")
 
 # benchmark data generated from Hartree.f
 # format: x  Rhos(n)  Rhov(n)  Rhot(n)  Rhos(p)  Rhov(p)  Rhot(p)

test/Pb208_nucleon_spectrum.jl

@@ -1,5 +1,5 @@
-using DelimitedFiles, Interpolations, Plots
+using DelimitedFiles, Plots
-include("../system.jl")
+include("../nucleons.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -15,19 +15,16 @@ N_p = 82
 N_n = 126
 r_max = maximum(xs)
 divs = length(xs) - 1
-s = system(N_p, N_n, r_max, divs)
+s = system(N_p, N_n, parameters(), r_max, divs)
 
 s.Φ0 = Ss
 s.W0 = Vs
 s.B0 = Rs
 s.A0 = As
 
-E_min = 800
+κs, Es = findAllOrbitals(p, s)
-E_max = 939
-
-κs, Es = findAllOrbitals(p, s, E_min, E_max)
 Ebinds = (p ? M_p : M_n) .- Es
-occ = fillNucleons(Z_or_N(s, p), κs, Es)
+spec = fillNucleons(Z_or_N(s, p), κs, Es)
 
 # format: proton, kappa, filling, gnodes, fnodes, Ebind
 bench_data, _ = readdlm("test/Pb208Spectrum.csv", ','; header=true)
@@ -38,7 +35,7 @@ bench_filling = bench_data[mask, 3]
 bench_occ = Int.(bench_filling .* (2 .* j_κ.(Int.(bench_κs)) .+ 1))
 
 scatter(κs, Ebinds, label=(p ? "proton" : "neutron") * " spectrum")
-annotate!(κs .+ 0.3, Ebinds .+ 1, [(o, 9) for o in occ])
+annotate!(κs .+ 0.3, Ebinds .+ 1, [(o, 9) for o in spec.occ])
 
 scatter!(bench_κs, bench_Ebinds, label="benchmark"; markershape=:cross, markercolor="red")
 annotate!(bench_κs .+ 0.3, bench_Ebinds .- 1, [(o, 9, :red) for o in bench_occ])

test/Pb208_nucleon_wf.jl

@@ -1,5 +1,5 @@
-using DelimitedFiles, Interpolations, Plots
+using DelimitedFiles, Plots
-include("../system.jl")
+include("../nucleons.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -21,10 +21,7 @@ s.W0 = Vs
 s.B0 = Rs
 s.A0 = As
 
-E_min = 800
+groundE = findEs(κ, p, s) |> minimum
-E_max = 939
-
-groundE = findEs(κ, p, s, E_min, E_max) |> minimum
 println("ground state E = $groundE")
 
 wf = solveNucleonWf(κ, p, groundE, s)

test/Pb208_shooting.jl

@@ -1,5 +1,5 @@
 using DelimitedFiles, Plots
-include("../system.jl")
+include("../nucleons.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -21,18 +21,18 @@ s.W0 = Vs
 s.B0 = Rs
 s.A0 = As
 
-E_min = 850
+E_min = 850.0
-E_max = 939
+E_max = 938.0
 
 boundEs = findEs(κ, p, s, E_min, E_max)
 
 binding_Es = (p ? M_p : M_n) .- boundEs
 println("binding energies = $binding_Es")
 
-func = boundaryValueFunc(κ, p, s)
+func = determinantFunc(κ, p, s)
-Es = collect(E_min:0.5:E_max)
+Es = collect(E_min:0.1:E_max)
-boundaryVals = [func(E)^2 for E in Es]
+dets = [func(E)^2 for E in Es]
 
-plot(Es, boundaryVals, yscale=:log10, label="g(r_max)^2")
+plot(Es, dets, yscale=:log10, label="determinant^2")
 vline!(boundEs, label="bound E")
 xlabel!("E (MeV)")