NuclearRMF.jl

@@ -1,68 +0,0 @@
-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,3 +3,4 @@ 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

@@ -1,62 +0,0 @@
-"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

@@ -1,17 +0,0 @@
-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,5 +1,23 @@
-include("common.jl")
+using DifferentialEquations
-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))
@@ -10,10 +28,7 @@ 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)))
@@ -26,8 +41,7 @@ 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=15, oscillation_control_parameter=0.3)
+function solveMesonFields!(s::system, iterations=10)
-    (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)
 
@@ -41,32 +55,10 @@ function solveMesonFields!(s::system, iterations=15, oscillation_control_paramet
         @. 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
-        # W0 and B0 keep a fraction of their previous solutions to suppress oscillations
+        B0 .= (solveKG(m_ρ, src_B0, s) .+ B0) ./ 2 # to supress oscillation
-        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,129 +1,88 @@
-using LinearAlgebra, DifferentialEquations, Interpolations
+using DifferentialEquations, Roots
-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,
-    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()),
+    Φ0, W0, B0, A0 are the mean-field potentials (couplings included) in MeV as functions of r in fm,
-    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
+function dirac!(du, u, (κ, p, E, Φ0, W0, B0, A0), r)
+    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
-    @inbounds du[1] = -(κ/(r + r_reg)) * g + (E + f1(r)) * f / ħc
+    du[1] = -(κ/(r + r_reg)) * g + (common1 + common2) * f / ħc
-    @inbounds du[2] =  (κ/(r + r_reg)) * f - (E + f2(r)) * g / ħc
+    du[2] =  (κ/(r + r_reg)) * f - (common1 - common2) * 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::Bool, E, s::system; normalize=true, algo=Vern9())
+function solveNucleonWf(κ, p, E, s::system; shooting=true, normalize=true)
-    (f1, f2) = optimized_dirac_potentials(p, s)
+    (Φ0, W0, B0, A0) = fields_interp(s)
 
-    # partitioning
+    if shooting
-    mid_idx = s.divs ÷ 2
+        @assert s.divs % 2 == 0 "divs must be an even number when shooting=true"
-    r_mid = rs(s)[mid_idx]
+        prob = ODEProblem(dirac!, [0, 1], (s.r_max, s.r_max / 2))
-    left_r = rs(s)[1:mid_idx]
+        sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr(s))
-    right_r = rs(s)[mid_idx:end]
+        wf_right = reverse(hcat(sol.u...); dims=2)
-
+        next_r_max = s.r_max / 2 # for the next step
-    # 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
+        next_r_max = s.r_max
-        wf_right .*= proj
+    end
+    
+    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
-        g2_int = simpsons_integrate(wf[1, :] .^ 2, Δr(s))
+        norm = sum(wf .* wf) * Δr(s) # integration by Reimann sum
-        f2_int = simpsons_integrate(wf[2, :] .^ 2, Δr(s))
+        wf = wf ./ sqrt(norm)
-        wf ./= sqrt(g2_int + f2_int)
     end
 
     return wf
 end
 
-"Returns a function that solves the Dirac equation in two partitions and returns 
+"Returns a function that solves the Dirac equation and returns g(r=r_max) where
-    the determinant of [g_left(r) g_right(r); f_left(r) f_right(r)],
+    r_max is the outer boundary in fm,
-    where is r is in fm."
+    the other parameters are the same from dirac!(...)."
-function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_max/2, algo=Vern9())
+function boundaryValueFunc(κ, p, s::system; dtype=Float64, algo=Tsit5())
-    (f1, f2) = optimized_dirac_potentials(p, s)
+    prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, s.r_max))
-    prob_left  = ODEProblem(dirac!, init_BC(), (0, r))
+    (Φ0, W0, B0, A0) = fields_interp(s)
-    function func(E)
+    func(E) = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=[s.r_max], save_idxs=[1])[1, 1]
-        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
+"Find all bound energies between E_min (=0) and E_max (=mass) 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)
+function findEs(κ, p, s::system, E_min=0, E_max=(p ? M_p : M_n), tol_digits=5)
-    func = determinantFunc(κ, p, s)
+    func = boundaryValueFunc(κ, p, s)
-    Es = find_all_zeros(func, E_min, E_max; partitions=20, tol=1/10^tol_digits)
+    Es = find_zeros(func, (E_min, E_max); xatol=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)
+function findAllOrbitals(p, s::system, E_min=0, E_max=(p ? M_p : M_n))
     κs = Int[]
     Es = Float64[]
     # start from κ=-1 and go both up and down
@@ -138,8 +97,8 @@ 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"
+"For a given list of κ values with corresponding energies, attempt to fill Z_or_N lowest lying orbitals and return occupancy numbers"
-function fillNucleons(Z_or_N::Int, κs, Es)::spectrum
+function fillNucleons(Z_or_N::Int, κs, Es)
     sort_i = sortperm(Es)
 
     occ = zeros(Int, length(κs))
@@ -151,8 +110,8 @@ function fillNucleons(Z_or_N::Int, κs, Es)::spectrum
         Z_or_N -= occ[i]
     end
 
-    @assert Z_or_N == 0 "All orbitals could not be filled"
+    Z_or_N == 0 || @warn "All orbitals could not be filled"
-    return spectrum(κs, Es, occ)
+    return occ
 end
 
 "Total angular momentum j for a given κ value"
@@ -164,14 +123,11 @@ 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(p::Bool, s::system)::Tuple{Vector{Float64}, Vector{Float64}}
+function calculateNucleonDensity(κs, Es, occs, p, s::system)
-    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; normalize=true)
+        wf = solveNucleonWf(κ, p, E, s; shooting=true, normalize=true)
         wf2 = wf .* wf
         ρr2 += (occ / (4 * pi)) * wf2 # 2j+1 factor is accounted in the occupancy number
     end
@@ -186,25 +142,10 @@ function calculateNucleonDensity(p::Bool, s::system)::Tuple{Vector{Float64}, Vec
     return (ρ_s, ρ_v)
 end
 
-"For a nucleon species, solve the Dirac equation and save the spectrum & densities in-place where
+"Solve the Dirac equation and calculate scalar and vector densities of a nucleon species where
     the other parameters are defined above"
-function solveNucleonDensity!(p::Bool, s::system, E_min=850.0, E_max=938.0)
+function solveNucleonDensity(p, s::system, E_min=800, E_max=939)
     κs, Es = findAllOrbitals(p, s, E_min, E_max)
-    spec = fillNucleons(Z_or_N(s, p), κs, Es)
+    occs = fillNucleons(Z_or_N(s, p), κs, Es)
-    if p
+    return calculateNucleonDensity(κs, Es, occs, p, s)
-        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,52 +1,13 @@
-"Stores a set of coupling constants and meson masses that goes into the Lagrangian"
+using Interpolations, PolyLog, Plots
-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}
@@ -58,10 +19,10 @@ mutable struct system
     ρ_vn::Vector{Float64}
 
     "Initialize an unsolved system"
-    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]...)
+    system(Z, N, r_max, divs) = new(Z, N, r_max, divs, [zeros(1 + divs) for _ in 1:8]...)
 
-    "Dummy struct to define the mesh (for testing)"
+    "Dummy struct to define the mesh"
-    system(r_max, divs) = system(0, 0, parameters(), r_max, divs)
+    system(r_max, divs) = system(0, 0, r_max, divs)
 end
 
 "Get mass number of nucleus"
@@ -78,3 +39,66 @@ 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, Plots
+using DelimitedFiles, Interpolations, Plots
-include("../nucleons.jl")
+include("../system.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -13,16 +13,19 @@ As = test_data[:, 5]
 p = false
 r_max = maximum(xs)
 divs = length(xs) - 1
-s = system(8, 8, parameters(), r_max, divs)
+s = system(8, 8, r_max, divs)
 
 s.Φ0 = Ss
 s.W0 = Vs
 s.B0 = Rs
 s.A0 = As
 
-κs, Es = findAllOrbitals(p, s)
+E_min = 860
+E_max = 939
+
+κs, Es = findAllOrbitals(p, s, E_min, E_max)
 Ebinds = (p ? M_p : M_n) .- Es
-spec = fillNucleons(Z_or_N(s, p), κs, Es)
+occ = fillNucleons(Z_or_N(s, p), κs, Es)
 
 # format: proton, kappa, filling, gnodes, fnodes, Ebind
 bench_data, _ = readdlm("test/LinearSpectrum.csv", ','; header=true)
@@ -33,7 +36,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 spec.occ])
+annotate!(κs .+ 0.3, Ebinds .+ 1, [(o, 9) for o in 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,14 +1,5 @@
-include("../NuclearRMF.jl")
+include("../system.jl")
 
-Z = 82
+s = system(82, 126, 20.0, 400)
-N = 126
 
-# Parameter values calibrated by FSUGarnet for Pb-208
+solve_system(s; monitor_plot=true)
-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

@@ -1 +0,0 @@
-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, Plots
+using DelimitedFiles, Interpolations, Plots
-include("../mesons.jl")
+include("../system.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -17,12 +17,9 @@ 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(N_p, N_n, param, r_max, divs)
+s = system(r_max, divs)
 
 s.ρ_sn = test_data[:, 2]
 s.ρ_vn = test_data[:, 3]

test/Pb208_nucleon_basic.jl

@@ -1,5 +1,5 @@
 using DelimitedFiles, Plots
-include("../nucleons.jl")
+include("../system.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.0
+E_min = 850
-E_max = 938.0
+E_max = 939
 
 boundEs = findEs(κ, p, s, E_min, E_max)
 
 binding_Es = (p ? M_p : M_n) .- boundEs
 println("binding energies = $binding_Es")
 
-func = determinantFunc(κ, p, s)
+func = boundaryValueFunc(κ, p, s)
-Es = collect(E_min:0.1:E_max)
+Es = collect(E_min:0.5:E_max)
-dets = [func(E)^2 for E in Es]
+boundaryVals = [func(E)^2 for E in Es]
 
-plot(Es, dets, yscale=:log10, label="determinant^2")
+plot(Es, boundaryVals, yscale=:log10, label="g(r_max)^2")
 vline!(boundEs, label="bound E")
 xlabel!("E (MeV)")

test/Pb208_nucleon_dens.jl

@@ -1,5 +1,5 @@
 using DelimitedFiles, Plots
-include("../nucleons.jl")
+include("../system.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, parameters(), r_max, divs)
+s = system(N_p, N_n, r_max, divs)
 
 s.Φ0 = Ss
 s.W0 = Vs
 s.B0 = Rs
 s.A0 = As
 
-solveNucleonDensity!(true, s)
+(ρ_sp, ρ_vp) = solveNucleonDensity(true, s)
 
-p_sp = plot(rs(s), s.ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
+p_sp = plot(rs(s), ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
-p_vp = plot(rs(s), s.ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
+p_vp = plot(rs(s), ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
 
-solveNucleonDensity!(false, s)
+(ρ_sn, ρ_vn) = solveNucleonDensity(false, s)
 
-p_sn = plot(rs(s), s.ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
+p_sn = plot(rs(s), ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
-p_vn = plot(rs(s), s.ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")
+p_vn = plot(rs(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, Plots
+using DelimitedFiles, Interpolations, Plots
-include("../nucleons.jl")
+include("../system.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -15,16 +15,19 @@ N_p = 82
 N_n = 126
 r_max = maximum(xs)
 divs = length(xs) - 1
-s = system(N_p, N_n, parameters(), r_max, divs)
+s = system(N_p, N_n, r_max, divs)
 
 s.Φ0 = Ss
 s.W0 = Vs
 s.B0 = Rs
 s.A0 = As
 
-κs, Es = findAllOrbitals(p, s)
+E_min = 800
+E_max = 939
+
+κs, Es = findAllOrbitals(p, s, E_min, E_max)
 Ebinds = (p ? M_p : M_n) .- Es
-spec = fillNucleons(Z_or_N(s, p), κs, Es)
+occ = fillNucleons(Z_or_N(s, p), κs, Es)
 
 # format: proton, kappa, filling, gnodes, fnodes, Ebind
 bench_data, _ = readdlm("test/Pb208Spectrum.csv", ','; header=true)
@@ -35,7 +38,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 spec.occ])
+annotate!(κs .+ 0.3, Ebinds .+ 1, [(o, 9) for o in 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, Plots
+using DelimitedFiles, Interpolations, Plots
-include("../nucleons.jl")
+include("../system.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -21,7 +21,10 @@ s.W0 = Vs
 s.B0 = Rs
 s.A0 = As
 
-groundE = findEs(κ, p, s) |> minimum
+E_min = 800
+E_max = 939
+
+groundE = findEs(κ, p, s, E_min, E_max) |> minimum
 println("ground state E = $groundE")
 
 wf = solveNucleonWf(κ, p, groundE, s)