diff --git a/networks/He-C-Fe-group/He-C-Fe-group.py b/networks/He-C-Fe-group/He-C-Fe-group.py
index e4429b1d1a..dfa374d31f 100644
--- a/networks/He-C-Fe-group/He-C-Fe-group.py
+++ b/networks/He-C-Fe-group/He-C-Fe-group.py
@@ -95,15 +95,10 @@
             d = pyna.DerivedRate(rate=fr, compute_Q=False, use_pf=True)
             all_lib.add_rate(d)
 
-# combine all three libraries into a single network
-
-net = pyna.AmrexAstroCxxNetwork(libraries=[all_lib],
-                                symmetric_screening=False)
-
 # we will have duplicate rates -- we want to remove any ReacLib rates
 # that we have tabular rates for
 
-dupes = net.find_duplicate_links()
+dupes = all_lib.find_duplicate_links()
 
 rates_to_remove = []
 for d in dupes:
@@ -111,7 +106,14 @@
         if isinstance(r, ReacLibRate):
             rates_to_remove.append(r)
 
-net.remove_rates(rates_to_remove)
+for r in rates_to_remove:
+    all_lib.remove_rate(r)
+
+# combine all three libraries into a single network
+
+net = pyna.AmrexAstroCxxNetwork(libraries=[all_lib],
+                                symmetric_screening=False)
+
 
 # now we approximate some (alpha, p)(p, gamma) links
 
diff --git a/networks/He-C-Fe-group/actual_rhs.H b/networks/He-C-Fe-group/actual_rhs.H
index e8eaf5d688..dfa273d450 100644
--- a/networks/He-C-Fe-group/actual_rhs.H
+++ b/networks/He-C-Fe-group/actual_rhs.H
@@ -1269,13 +1269,15 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
 
     [[maybe_unused]] Real rate, drate_dt, edot_nu, edot_gamma;
 
+    rate_eval.enuc_weak = 0.0;
+
     tabular_evaluate(j_Co55_Fe55_meta, j_Co55_Fe55_rhoy, j_Co55_Fe55_temp, j_Co55_Fe55_data,
                      rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
     rate_eval.screened_rates(k_Co55_to_Fe55) = rate;
     if constexpr (std::is_same<T, rate_derivs_t>::value) {
         rate_eval.dscreened_rates_dT(k_Co55_to_Fe55) = drate_dt;
     }
-    rate_eval.add_energy_rate(k_Co55_to_Fe55) = edot_nu + edot_gamma;
+    rate_eval.enuc_weak += C::Legacy::n_A * Y(Co55) * (edot_nu + edot_gamma);
 
     tabular_evaluate(j_Co56_Fe56_meta, j_Co56_Fe56_rhoy, j_Co56_Fe56_temp, j_Co56_Fe56_data,
                      rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
@@ -1283,7 +1285,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
     if constexpr (std::is_same<T, rate_derivs_t>::value) {
         rate_eval.dscreened_rates_dT(k_Co56_to_Fe56) = drate_dt;
     }
-    rate_eval.add_energy_rate(k_Co56_to_Fe56) = edot_nu + edot_gamma;
+    rate_eval.enuc_weak += C::Legacy::n_A * Y(Co56) * (edot_nu + edot_gamma);
 
     tabular_evaluate(j_Co56_Ni56_meta, j_Co56_Ni56_rhoy, j_Co56_Ni56_temp, j_Co56_Ni56_data,
                      rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
@@ -1291,7 +1293,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
     if constexpr (std::is_same<T, rate_derivs_t>::value) {
         rate_eval.dscreened_rates_dT(k_Co56_to_Ni56) = drate_dt;
     }
-    rate_eval.add_energy_rate(k_Co56_to_Ni56) = edot_nu + edot_gamma;
+    rate_eval.enuc_weak += C::Legacy::n_A * Y(Co56) * (edot_nu + edot_gamma);
 
     tabular_evaluate(j_Co57_Ni57_meta, j_Co57_Ni57_rhoy, j_Co57_Ni57_temp, j_Co57_Ni57_data,
                      rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
@@ -1299,7 +1301,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
     if constexpr (std::is_same<T, rate_derivs_t>::value) {
         rate_eval.dscreened_rates_dT(k_Co57_to_Ni57) = drate_dt;
     }
-    rate_eval.add_energy_rate(k_Co57_to_Ni57) = edot_nu + edot_gamma;
+    rate_eval.enuc_weak += C::Legacy::n_A * Y(Co57) * (edot_nu + edot_gamma);
 
     tabular_evaluate(j_Fe55_Co55_meta, j_Fe55_Co55_rhoy, j_Fe55_Co55_temp, j_Fe55_Co55_data,
                      rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
@@ -1307,7 +1309,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
     if constexpr (std::is_same<T, rate_derivs_t>::value) {
         rate_eval.dscreened_rates_dT(k_Fe55_to_Co55) = drate_dt;
     }
-    rate_eval.add_energy_rate(k_Fe55_to_Co55) = edot_nu + edot_gamma;
+    rate_eval.enuc_weak += C::Legacy::n_A * Y(Fe55) * (edot_nu + edot_gamma);
 
     tabular_evaluate(j_Fe55_Mn55_meta, j_Fe55_Mn55_rhoy, j_Fe55_Mn55_temp, j_Fe55_Mn55_data,
                      rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
@@ -1315,7 +1317,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
     if constexpr (std::is_same<T, rate_derivs_t>::value) {
         rate_eval.dscreened_rates_dT(k_Fe55_to_Mn55) = drate_dt;
     }
-    rate_eval.add_energy_rate(k_Fe55_to_Mn55) = edot_nu + edot_gamma;
+    rate_eval.enuc_weak += C::Legacy::n_A * Y(Fe55) * (edot_nu + edot_gamma);
 
     tabular_evaluate(j_Fe56_Co56_meta, j_Fe56_Co56_rhoy, j_Fe56_Co56_temp, j_Fe56_Co56_data,
                      rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
@@ -1323,7 +1325,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
     if constexpr (std::is_same<T, rate_derivs_t>::value) {
         rate_eval.dscreened_rates_dT(k_Fe56_to_Co56) = drate_dt;
     }
-    rate_eval.add_energy_rate(k_Fe56_to_Co56) = edot_nu + edot_gamma;
+    rate_eval.enuc_weak += C::Legacy::n_A * Y(Fe56) * (edot_nu + edot_gamma);
 
     tabular_evaluate(j_Mn55_Fe55_meta, j_Mn55_Fe55_rhoy, j_Mn55_Fe55_temp, j_Mn55_Fe55_data,
                      rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
@@ -1331,7 +1333,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
     if constexpr (std::is_same<T, rate_derivs_t>::value) {
         rate_eval.dscreened_rates_dT(k_Mn55_to_Fe55) = drate_dt;
     }
-    rate_eval.add_energy_rate(k_Mn55_to_Fe55) = edot_nu + edot_gamma;
+    rate_eval.enuc_weak += C::Legacy::n_A * Y(Mn55) * (edot_nu + edot_gamma);
 
     tabular_evaluate(j_Ni56_Co56_meta, j_Ni56_Co56_rhoy, j_Ni56_Co56_temp, j_Ni56_Co56_data,
                      rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
@@ -1339,7 +1341,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
     if constexpr (std::is_same<T, rate_derivs_t>::value) {
         rate_eval.dscreened_rates_dT(k_Ni56_to_Co56) = drate_dt;
     }
-    rate_eval.add_energy_rate(k_Ni56_to_Co56) = edot_nu + edot_gamma;
+    rate_eval.enuc_weak += C::Legacy::n_A * Y(Ni56) * (edot_nu + edot_gamma);
 
     tabular_evaluate(j_Ni57_Co57_meta, j_Ni57_Co57_rhoy, j_Ni57_Co57_temp, j_Ni57_Co57_data,
                      rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
@@ -1347,7 +1349,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
     if constexpr (std::is_same<T, rate_derivs_t>::value) {
         rate_eval.dscreened_rates_dT(k_Ni57_to_Co57) = drate_dt;
     }
-    rate_eval.add_energy_rate(k_Ni57_to_Co57) = edot_nu + edot_gamma;
+    rate_eval.enuc_weak += C::Legacy::n_A * Y(Ni57) * (edot_nu + edot_gamma);
 
     tabular_evaluate(j_n_p_meta, j_n_p_rhoy, j_n_p_temp, j_n_p_data,
                      rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
@@ -1355,7 +1357,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
     if constexpr (std::is_same<T, rate_derivs_t>::value) {
         rate_eval.dscreened_rates_dT(k_n_to_p) = drate_dt;
     }
-    rate_eval.add_energy_rate(k_n_to_p) = edot_nu + edot_gamma;
+    rate_eval.enuc_weak += C::Legacy::n_A * Y(N) * (edot_nu + edot_gamma);
 
     tabular_evaluate(j_p_n_meta, j_p_n_rhoy, j_p_n_temp, j_p_n_data,
                      rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
@@ -1363,7 +1365,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
     if constexpr (std::is_same<T, rate_derivs_t>::value) {
         rate_eval.dscreened_rates_dT(k_p_to_n) = drate_dt;
     }
-    rate_eval.add_energy_rate(k_p_to_n) = edot_nu + edot_gamma;
+    rate_eval.enuc_weak += C::Legacy::n_A * Y(H1) * (edot_nu + edot_gamma);
 
 
 }
@@ -1724,6 +1726,7 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     rate_t rate_eval;
 
     constexpr int do_T_derivatives = 0;
+
     evaluate_rates<do_T_derivatives, rate_t>(state, rate_eval);
 
     rhs_nuc(state, ydot, Y, rate_eval.screened_rates);
@@ -1733,19 +1736,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     Real enuc;
     ener_gener_rate(ydot, enuc);
 
-    // include reaction neutrino losses (non-thermal) and gamma heating
-    enuc += C::Legacy::n_A * Y(Co55) * rate_eval.add_energy_rate(k_Co55_to_Fe55);
-    enuc += C::Legacy::n_A * Y(Co56) * rate_eval.add_energy_rate(k_Co56_to_Fe56);
-    enuc += C::Legacy::n_A * Y(Co56) * rate_eval.add_energy_rate(k_Co56_to_Ni56);
-    enuc += C::Legacy::n_A * Y(Co57) * rate_eval.add_energy_rate(k_Co57_to_Ni57);
-    enuc += C::Legacy::n_A * Y(Fe55) * rate_eval.add_energy_rate(k_Fe55_to_Co55);
-    enuc += C::Legacy::n_A * Y(Fe55) * rate_eval.add_energy_rate(k_Fe55_to_Mn55);
-    enuc += C::Legacy::n_A * Y(Fe56) * rate_eval.add_energy_rate(k_Fe56_to_Co56);
-    enuc += C::Legacy::n_A * Y(Mn55) * rate_eval.add_energy_rate(k_Mn55_to_Fe55);
-    enuc += C::Legacy::n_A * Y(Ni56) * rate_eval.add_energy_rate(k_Ni56_to_Co56);
-    enuc += C::Legacy::n_A * Y(Ni57) * rate_eval.add_energy_rate(k_Ni57_to_Co57);
-    enuc += C::Legacy::n_A * Y(N) * rate_eval.add_energy_rate(k_n_to_p);
-    enuc += C::Legacy::n_A * Y(H1) * rate_eval.add_energy_rate(k_p_to_n);
+    // include any weak rate neutrino losses
+    enuc += rate_eval.enuc_weak;
 
     // Get the thermal neutrino losses
 
@@ -2770,6 +2762,7 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     rate_derivs_t rate_eval;
 
     constexpr int do_T_derivatives = 1;
+
     evaluate_rates<do_T_derivatives, rate_derivs_t>(state, rate_eval);
 
     // Species Jacobian elements with respect to other species
diff --git a/networks/He-C-Fe-group/reaclib_rates.H b/networks/He-C-Fe-group/reaclib_rates.H
index 0c8d705f8d..9e5dd4ea56 100644
--- a/networks/He-C-Fe-group/reaclib_rates.H
+++ b/networks/He-C-Fe-group/reaclib_rates.H
@@ -13,13 +13,13 @@ using namespace Species;
 
 struct rate_t {
     Array1D<Real, 1, NumRates>  screened_rates;
-    Array1D<Real, NrateReaclib+1, NrateReaclib+NrateTabular> add_energy_rate;
+    Real enuc_weak;
 };
 
 struct rate_derivs_t {
     Array1D<Real, 1, NumRates>  screened_rates;
     Array1D<Real, 1, NumRates>  dscreened_rates_dT;
-    Array1D<Real, NrateReaclib+1, NrateReaclib+NrateTabular> add_energy_rate;
+    Real enuc_weak;
 };
 
 
diff --git a/networks/He-C-Fe-group/table_rates.H b/networks/He-C-Fe-group/table_rates.H
index 5b38af6e84..05743c4b44 100644
--- a/networks/He-C-Fe-group/table_rates.H
+++ b/networks/He-C-Fe-group/table_rates.H
@@ -126,6 +126,11 @@ void init_tab_info(const table_t& tf, const std::string& file, R& log_rhoy_table
     std::ifstream table;
     table.open(file);
 
+    if (!table.is_open()) {
+        // the table was not present or we could not open it; abort
+        amrex::Error("table could not be opened");
+    }
+
     std::string line;
 
     // read and skip over the header
@@ -139,6 +144,10 @@ void init_tab_info(const table_t& tf, const std::string& file, R& log_rhoy_table
     for (int j = 1; j <= tf.nrhoy; ++j) {
         for (int i = 1; i <= tf.ntemp; ++i) {
             std::getline(table, line);
+            if (line.empty()) {
+                amrex::Error("Error reading table data");
+            }
+
             std::istringstream sdata(line);
 
             sdata >> log_rhoy_table(j) >> log_temp_table(i);
@@ -305,8 +314,8 @@ evaluate_dr_dtemp(const table_t& table_meta, const R& log_rhoy_table, const T& l
         // Finally, we perform a 1d-linear interpolation between dlogr_dlogt_ip1
         // and dlogr_dlogt_i to compute dlogr_dlogt
 
-        Real log_rhoy_lo = log_temp_table(irhoy_lo);
-        Real log_rhoy_hi = log_temp_table(irhoy_hi);
+        Real log_rhoy_lo = log_rhoy_table(irhoy_lo);
+        Real log_rhoy_hi = log_rhoy_table(irhoy_hi);
 
         Real log_temp_lo = log_temp_table(jtemp_lo);
         Real log_temp_hi = log_temp_table(jtemp_hi);
diff --git a/unit_test/react_util.H b/unit_test/react_util.H
index 0c970fd08f..88a72b9112 100644
--- a/unit_test/react_util.H
+++ b/unit_test/react_util.H
@@ -138,8 +138,8 @@ void get_xn(const int k, const init_t cd, Real *xn_zone, int uniform_composition
 
         for (int n = 0; n < NumSpec; n++) {
             if (n == cd.is1 ||
-                (n == cd.is2 && cd.nprim >= 2) ||
-                (n == cd.is3 && cd.nprim >= 3)) {
+                (cd.nprim >= 2 && n == cd.is2) ||
+                (cd.nprim >= 3 && n == cd.is3)) {
                 continue;
             }
 
@@ -286,4 +286,3 @@ Real get_xn(const int index) {
 }
 
 #endif
-