Model parameters are no longer hard-coded

NuclearRMF.jl

@@ -3,13 +3,27 @@ include("nucleons.jl")
 include("mesons.jl")
 
 "Total binding energy of the system"
-total_E(s::system) = -(nucleon_E(s) + meson_E(s))
+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, monitor_print=true, monitor_plot=false)
+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
@@ -18,18 +32,19 @@ function solve_system!(s::system; reinitialize_densities=true, monitor_print=tru
         @. s.ρ_vn = s.N * dens_guess
     end
 
-    if monitor_plot
+    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
-        @time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 15)
+        enable_time = print_time
-        @time "Proton densities" solveNucleonDensity!(true, s)
+        @conditional_time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 15)
-        @time "Neutron densities" solveNucleonDensity!(false, s)
+        @conditional_time "Proton densities" solveNucleonDensity!(true, s)
+        @conditional_time "Neutron densities" solveNucleonDensity!(false, s)
 
-        if monitor_plot
+        if live_plots
             for s in p.series_list
                 s.plotattributes[:linecolor] = :gray
             end
@@ -38,10 +53,16 @@ function solve_system!(s::system; reinitialize_densities=true, monitor_print=tru
         end
 
         E_per_A = total_E(s) / A(s)
-        monitor_print && println("Total binding E per nucleon = $E_per_A")
+        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)

common.jl

@@ -1,2 +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,21 +1,6 @@
 include("common.jl")
 include("system.jl")
 
-# 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
-
 "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))
 
@@ -25,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)))
@@ -39,6 +27,7 @@ end
     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)
+    (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)
 
@@ -68,14 +57,16 @@ end
 
 "Calculate the total energy associated with meson fields"
 function meson_E(s::system)
-    int = 0.0
+    (κ_ss, λ, ζ, Λv) = (s.param.κ_ss, s.param.λ, s.param.ζ, s.param.Λv)
-    for (r, Φ0, W0, B0, A0, ρ_sp, ρ_vp, ρ_sn, ρ_vn) in zip(rs(s), s.Φ0, s.W0, s.B0, s.A0, s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)
+
+    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
-        int += (E_σ + E_ω + E_ρ + E_γ + E_ωρ) * r^2
+        E_σ + E_ω + E_ρ + E_γ + E_ωρ
     end
-    return 4π * int * Δr(s) * ħc
+    
+    return simpsons_integrate(E_densities .* rs(s).^2, Δr(s); coefficient=4π * ħc)
 end

nucleons.jl

@@ -89,7 +89,9 @@ function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Vern9())
     wf = hcat(wf_left[:, 1:(end - 1)], wf_right)
 
     if normalize
-        wf ./= norm(wf) * sqrt(Δr(s)) # integration by Reimann sum
+        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

system.jl

@@ -1,3 +1,30 @@
+"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}
@@ -13,6 +40,7 @@ mutable struct system
     Z::Int
     N::Int
 
+    param::parameters
     r_max::Float64
     divs::Int
 
@@ -30,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, unfilled_spectrum(), unfilled_spectrum(), [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"

test/Linear_nucleon_spectrum.jl

@@ -13,7 +13,7 @@ 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

test/Pb208.jl

@@ -2,9 +2,13 @@ include("../NuclearRMF.jl")
 
 Z = 82
 N = 126
+
+# 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, r_max, divs)
+s = system(Z, N, param, r_max, divs)
 
-solve_system!(s; monitor_plot=true)
+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

@@ -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

@@ -14,7 +14,7 @@ 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

test/Pb208_nucleon_spectrum.jl

@@ -15,7 +15,7 @@ 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