diff --git a/src/include/gridfire/engine/types/jacobian.h b/src/include/gridfire/engine/types/jacobian.h
new file mode 100644
index 00000000..e69de29b
diff --git a/src/include/gridfire/reaction/reaction.h b/src/include/gridfire/reaction/reaction.h
index 4b054333..997b8396 100644
--- a/src/include/gridfire/reaction/reaction.h
+++ b/src/include/gridfire/reaction/reaction.h
@@ -217,6 +217,10 @@ namespace gridfire::reaction {
          */
         [[nodiscard]] virtual std::unordered_set<fourdst::atomic::Species> product_species() const = 0;
 
+        [[nodiscard]] virtual size_t countReactantOccurrences(const fourdst::atomic::Species& species) const = 0;
+
+        [[nodiscard]] virtual size_t countProductOccurrences(const fourdst::atomic::Species& species) const = 0;
+
         /**
          * @brief Number of unique species involved in the reaction.
          * @return Count of distinct species across reactants and products.
@@ -525,8 +529,10 @@ namespace gridfire::reaction {
             if (!m_reactantsVec) {
                 m_reactantsVec.emplace(std::vector<fourdst::atomic::Species>());
                 m_reactantsVec->reserve(m_reactants.size());
-                for (const auto& reactant : m_reactants) {
-                    m_reactantsVec->push_back(reactant);
+                for (const auto& [reactant, count] : m_reactants) {
+                    for (size_t i = 0; i < count; ++i) {
+                        m_reactantsVec->push_back(reactant);
+                    }
                 }
             }
             return m_reactantsVec.value();
@@ -540,8 +546,10 @@ namespace gridfire::reaction {
             if (!m_productsVec) {
                 m_productsVec.emplace(std::vector<fourdst::atomic::Species>());
                 m_productsVec->reserve(m_products.size());
-                for (const auto& product : m_products) {
-                    m_productsVec->push_back(product);
+                for (const auto& [product, count] : m_products) {
+                    for (size_t i = 0; i < count; ++i) {
+                        m_productsVec->push_back(product);
+                    }
                 }
             }
             return m_productsVec.value();
@@ -587,14 +595,24 @@ namespace gridfire::reaction {
             return os << "(ReaclibReaction:" << r.m_id << ")";
         }
 
+        size_t countReactantOccurrences(const fourdst::atomic::Species& species) const override {
+            if (!m_reactants.contains(species)) return 0;
+            return m_reactants.at(species);
+        }
+
+        size_t countProductOccurrences(const fourdst::atomic::Species& species) const override {
+            if (!m_products.contains(species)) return 0;
+            return m_products.at(species);
+        }
+
     protected:
         quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
         std::string m_id; ///< Unique identifier for the reaction (e.g., "h1+h1=>h2+e+nu").
         std::string m_peName; ///< Name of the reaction in (projectile, ejectile) notation (e.g. "p(p,g)d").
         int m_chapter;    ///< Chapter number from the REACLIB database, defining the reaction structure.
         double m_qValue = 0.0; ///< Q-value of the reaction in MeV.
-        std::set<fourdst::atomic::Species> m_reactants; ///< Reactants of the reaction.
-        std::set<fourdst::atomic::Species> m_products; ///< Products of the reaction.
+        std::unordered_map<fourdst::atomic::Species, size_t> m_reactants; ///< Reactants of the reaction.
+        std::unordered_map<fourdst::atomic::Species, size_t> m_products; ///< Products of the reaction.
 
         mutable std::optional<std::vector<fourdst::atomic::Species>> m_reactantsVec;
         mutable std::optional<std::vector<fourdst::atomic::Species>> m_productsVec;
@@ -836,6 +854,8 @@ namespace gridfire::reaction {
          */
         [[nodiscard]] size_t size() const { return m_reactions.size(); }
 
+        [[nodiscard]] bool empty() const {return m_reactions.empty(); }
+
         /**
          * @brief Removes all reactions from the set.
          */
diff --git a/src/include/gridfire/reaction/weak/weak.h b/src/include/gridfire/reaction/weak/weak.h
index 00659fe1..0428a1e9 100644
--- a/src/include/gridfire/reaction/weak/weak.h
+++ b/src/include/gridfire/reaction/weak/weak.h
@@ -400,6 +400,11 @@ namespace gridfire::rates::weak {
          */
         double get_log_neutrino_loss_from_payload(const WeakRatePayload& payload) const;
 
+    public:
+        [[nodiscard]] size_t countReactantOccurrences(const fourdst::atomic::Species &species) const override;
+
+        [[nodiscard]] size_t countProductOccurrences(const fourdst::atomic::Species &species) const override;
+
     private:
         /**
          * @brief CppAD atomic that wraps weak-rate interpolation for AD evaluation.
diff --git a/src/lib/engine/types/jacobian.cpp b/src/lib/engine/types/jacobian.cpp
new file mode 100644
index 00000000..e69de29b
diff --git a/src/lib/reaction/reaction.cpp b/src/lib/reaction/reaction.cpp
index 3a227243..6bdd8af7 100644
--- a/src/lib/reaction/reaction.cpp
+++ b/src/lib/reaction/reaction.cpp
@@ -40,11 +40,16 @@ namespace gridfire::reaction {
     m_peName(peName),
     m_chapter(chapter),
     m_qValue(qValue),
-    m_reactants(reactants.begin(), reactants.end()),
-    m_products(products.begin(), products.end()),
     m_sourceLabel(label),
     m_rateCoefficients(sets),
-    m_reverse(reverse) {}
+    m_reverse(reverse) {
+        for (const auto& reactant : reactants) {
+            m_reactants[reactant]++;
+        }
+        for (const auto& product : products) {
+            m_products[product]++;
+        }
+    }
 
     double ReaclibReaction::calculate_rate(
         const double T9,
@@ -107,11 +112,11 @@ namespace gridfire::reaction {
 
 
     bool ReaclibReaction::contains_reactant(const Species& species) const {
-        return std::ranges::any_of(m_reactants, [&](const Species& r) { return r == species; });
+        return m_reactants.contains(species);
     }
 
     bool ReaclibReaction::contains_product(const Species& species) const {
-        return std::ranges::any_of(m_products, [&](const Species& p) { return p == species; });
+        return m_products.contains(species);
     }
 
     std::unordered_set<Species> ReaclibReaction::all_species() const {
@@ -123,7 +128,7 @@ namespace gridfire::reaction {
 
     std::unordered_set<Species> ReaclibReaction::reactant_species() const {
         std::unordered_set<Species> reactantsSet;
-        for (const auto& reactant : m_reactants) {
+        for (const auto& reactant : m_reactants | std::views::keys) {
             reactantsSet.insert(reactant);
         }
         return reactantsSet;
@@ -131,7 +136,7 @@ namespace gridfire::reaction {
 
     std::unordered_set<Species> ReaclibReaction::product_species() const {
         std::unordered_set<Species> productsSet;
-        for (const auto& product : m_products) {
+        for (const auto& product : m_products | std::views::keys) {
             productsSet.insert(product);
         }
         return productsSet;
@@ -139,14 +144,14 @@ namespace gridfire::reaction {
 
     int ReaclibReaction::stoichiometry(const Species& species) const {
         int s = 0;
-        for (const auto& reactant : m_reactants) {
+        for (const auto& [reactant, count] : m_reactants) {
             if (reactant == species) {
-                s--;
+                s -= count;
             }
         }
-        for (const auto& product : m_products) {
+        for (const auto& [product, count] : m_products) {
             if (product == species) {
-                s++;
+                s += count;
             }
         }
         return s;
@@ -158,11 +163,8 @@ namespace gridfire::reaction {
 
     std::unordered_map<Species, int> ReaclibReaction::stoichiometry() const {
         std::unordered_map<Species, int> stoichiometryMap;
-        for (const auto& reactant : m_reactants) {
-            stoichiometryMap[reactant]--;
-        }
-        for (const auto& product : m_products) {
-            stoichiometryMap[product]++;
+        for (const auto& sp : all_species()) {
+            stoichiometryMap[sp] = stoichiometry(sp);
         }
         return stoichiometryMap;
     }
@@ -171,10 +173,10 @@ namespace gridfire::reaction {
         double reactantMass = 0.0;
         double productMass = 0.0;
         constexpr double AMU2MeV = 931.494893; // Conversion factor from atomic mass unit to MeV
-        for (const auto& reactant : m_reactants) {
+        for (const auto& reactant : m_reactants | std::views::keys) {
             reactantMass += reactant.mass();
         }
-        for (const auto& product : m_products) {
+        for (const auto& product : m_products | std::views::keys) {
             productMass += product.mass();
         }
         return (reactantMass - productMass) * AMU2MeV;
diff --git a/src/lib/reaction/weak/weak.cpp b/src/lib/reaction/weak/weak.cpp
index 93e75c08..34d57eae 100644
--- a/src/lib/reaction/weak/weak.cpp
+++ b/src/lib/reaction/weak/weak.cpp
@@ -490,6 +490,26 @@ namespace gridfire::rates::weak {
         return logNeutrinoLoss;
     }
 
+    size_t WeakReaction::countReactantOccurrences(const fourdst::atomic::Species &species) const {
+        size_t count = 0;
+        for (const auto& reactant : m_reactants) {
+            if (reactant == species) {
+                count++;
+            }
+        }
+        return count;
+    }
+
+    size_t WeakReaction::countProductOccurrences(const fourdst::atomic::Species &species) const {
+        size_t count = 0;
+        for (const auto& product : m_products) {
+            if (product == species) {
+                count++;
+            }
+        }
+        return count;
+    }
+
     // Note that the input vector tx is of size 2: [T9, log10(rho*Ye)]
     bool WeakReaction::AtomicWeakRate::forward (
         const size_t p,