diff --git a/NuclearRMF.jl b/NuclearRMF.jl
index cdf314a..ab951fa 100644
--- a/NuclearRMF.jl
+++ b/NuclearRMF.jl
@@ -3,7 +3,7 @@ include("nucleons.jl")
 include("mesons.jl")
 
 "Total binding energy of the system"
-total_E(s::system, Es_p, occs_p, Es_n, occs_n) = sum(occs_p .* (M_p .- Es_p)) + sum(occs_n .* (M_n .- Es_n)) - meson_E(s)
+total_E(s::system) = -(nucleon_E(s) + meson_E(s))
 
 "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)))
@@ -29,14 +29,14 @@ function solve_system(s::system; reinitialize_densities=true, monitor_print=true
         @time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 15)
 
         # 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)
+        @time "Proton spectrum" (s.κ_p, s.E_p) = findAllOrbitals(true, s)
+        s.occ_p = fillNucleons(s.Z, s.κ_p, s.E_p)
+        @time "Proton densities" (s.ρ_sp, s.ρ_vp) = calculateNucleonDensity(s.κ_p, s.E_p, s.occ_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)
+        @time "Neutron spectrum" (s.κ_n, s.E_n) = findAllOrbitals(false, s)
+        s.occ_n = fillNucleons(s.N, s.κ_n, s.E_n)
+        @time "Neutron densities" (s.ρ_sn, s.ρ_vn) = calculateNucleonDensity(s.κ_n, s.E_n, s.occ_n, false, s)
 
         if monitor_plot
             for s in p.series_list
@@ -46,7 +46,7 @@ function solve_system(s::system; reinitialize_densities=true, monitor_print=true
             display(p)
         end
 
-        E_per_A = total_E(s, Es_p, occs_p, Es_n, occs_n) / A(s)
+        E_per_A = total_E(s) / A(s)
         monitor_print && println("Total binding E per nucleon = $E_per_A")
 
         # check convergence
diff --git a/nucleons.jl b/nucleons.jl
index 98f3c0b..095cc2b 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -170,3 +170,6 @@ function solveNucleonDensity(p::Bool, s::system, E_min=850.0, E_max=938.0)
     occs = fillNucleons(Z_or_N(s, p), κs, Es)
     return calculateNucleonDensity(κs, Es, occs, p, s)
 end
+
+"Total energy of filled nucleons in the system"
+nucleon_E(s::system) = sum(s.occ_p .* (s.E_p .- M_p)) + sum(s.occ_n .* (s.E_n .- M_n))
diff --git a/system.jl b/system.jl
index 531e15a..6df5e62 100644
--- a/system.jl
+++ b/system.jl
@@ -8,6 +8,16 @@ mutable struct system
     r_max::Float64
     divs::Int
 
+    # single particle energies and corresponding occupancy numbers for protons
+    κ_p::Vector{Int}
+    E_p::Vector{Float64}
+    occ_p::Vector{Int}
+
+    # single particle energies and corresponding occupancy numbers for neutrons
+    κ_n::Vector{Int}
+    E_n::Vector{Float64}
+    occ_n::Vector{Int}
+
     Φ0::Vector{Float64}
     W0::Vector{Float64}
     B0::Vector{Float64}
@@ -19,7 +29,7 @@ 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, r_max, divs) = new(Z, N, r_max, divs, Int[], Float64[], Int[], Int[], Float64[], Int[], [zeros(1 + divs) for _ in 1:8]...)
 
     "Dummy struct to define the mesh"
     system(r_max, divs) = system(0, 0, r_max, divs)