diff --git a/doc/_static/zotero.bib b/doc/_static/zotero.bib
index 837492d408..36c84fa9a9 100644
--- a/doc/_static/zotero.bib
+++ b/doc/_static/zotero.bib
@@ -3760,3 +3760,19 @@ @article{yamashita-muonicspin-2022
     url = {https://www.nature.com/articles/s41598-022-09487-0},
 	pages = {6393}
 }
+
+@article{faifman-mucfformation-1996,
+	title = {Quadrupole corrections to matrix elements of transitions in resonant reactions of muonic molecule formation},
+    author = {Faifman, M. P. and Strizh, T. A. and Armour, E. A. G. and Harston, M. R.},
+    year = 1996,
+    month = dec,
+	volume = {101-102},
+	rights = {http://www.springer.com/tdm},
+	issn = {0304-3834, 1572-9540},
+	url = {http://link.springer.com/10.1007/BF02227621},
+	doi = {10.1007/BF02227621},
+	pages = {179--189},
+	number = {1},
+	journaltitle = {Hyperfine Interactions},
+	shortjournal = {Hyperfine Interact},
+}
diff --git a/doc/implementation/mucf-physics.rst b/doc/implementation/mucf-physics.rst
index ff014d715d..819c263f85 100644
--- a/doc/implementation/mucf-physics.rst
+++ b/doc/implementation/mucf-physics.rst
@@ -127,7 +127,7 @@ formation, molecule formation, and fusion.
 .. doxygenclass:: celeritas::DTMixMucfModel
 
 Most of the data is material-dependent, being calculated and cached during model
-construction. All of the cached quantities are calculated and added to
+construction. All the cached quantities are calculated and added to
 host/device data via :cpp:class:`celeritas::detail::MucfMaterialInserter`.
 
 .. doxygenclass:: celeritas::detail::MucfMaterialInserter
diff --git a/src/celeritas/CMakeLists.txt b/src/celeritas/CMakeLists.txt
index 6188974647..e32b537b1c 100644
--- a/src/celeritas/CMakeLists.txt
+++ b/src/celeritas/CMakeLists.txt
@@ -107,6 +107,8 @@ list(APPEND SOURCES
   mat/MaterialParams.cc
   mat/MaterialParamsOutput.cc
   mat/detail/Utils.cc
+  mucf/model/detail/EquilibrateDensitiesCalculator.cc
+  mucf/model/detail/InterpolatorHelper.cc
   mucf/model/detail/MucfMaterialInserter.cc
   mucf/process/MucfProcess.cc
   neutron/model/CascadeOptions.cc
diff --git a/src/celeritas/inp/MucfPhysics.cc b/src/celeritas/inp/MucfPhysics.cc
index 6a51bf5f5e..d201688905 100644
--- a/src/celeritas/inp/MucfPhysics.cc
+++ b/src/celeritas/inp/MucfPhysics.cc
@@ -16,7 +16,8 @@ namespace
 /*!
  * Muon energy CDF data for muon-catalyzed fusion.
  *
- * Data is extracted from https://doi.org/10.1103/PhysRevC.107.034607 .
+ * Data is extracted from \citet{kamimura-mucf-2023,
+ * https://doi.org/10.1103/PhysRevC.107.034607}.
  */
 Grid mucf_muon_energy_cdf()
 {
@@ -77,14 +78,480 @@ Grid mucf_muon_energy_cdf()
     return cdf;
 }
 
+//---------------------------------------------------------------------------//
+// DD cycle rate data
+//---------------------------------------------------------------------------//
+/*!
+ * Muon-catalyzed fusion cycle rate data for \em dd fusion with F = 1/2 state.
+ */
+MucfCycleRate dd_1_over_2_cycle_data()
+{
+    MucfCycleRate result;
+    result.type = CycleTableType::deuterium_deuterium;
+    result.spin_state = units::HalfSpinInt{1};  // F = 1/2
+
+    result.rate.interpolation.type = InterpolationType::cubic_spline;
+    // Temperature [K]
+    result.rate.x = {
+        0.0,
+        36.43103676856282,
+        49.84753686046977,
+        62.30013430590361,
+        72.8249264583911,
+        94.80254726485413,
+        118.65929788178391,
+        134.8637502794038,
+        174.9275151614521,
+        212.16105692405534,
+        260.90506935936025,
+        325.03565799325077,
+        380.6089687138337,
+        434.29674495520254,
+        486.0836155123639,
+        564.7393523983244,
+        668.3489596576321,
+        764.2819309563899,
+        860.2110594538995,
+        935.0316806939586,
+        998.3392693938174,
+    };
+    // Mean cycle rate [1/s]
+    result.rate.y = {
+        0.0,
+        0.005464463375221662e6,
+        0.020292872692477815e6,
+        0.04715981762065269e6,
+        0.098103833770665e6,
+        0.2721222062034343e6,
+        0.5362926959236205e6,
+        0.7584940318094278e6,
+        1.2388537948623028e6,
+        1.6110396634730852e6,
+        1.9710284809214524e6,
+        2.2225696043387786e6,
+        2.3059026708109145e6,
+        2.3141137763784285e6,
+        2.283275296360035e6,
+        2.19792475923249e6,
+        2.05207892345204e6,
+        1.921378528091838e6,
+        1.7996962265613146e6,
+        1.7144033314524956e6,
+        1.6473195500592115e6,
+    };
+
+    CELER_ENSURE(result);
+    return result;
+}
+//---------------------------------------------------------------------------//
+/*!
+ * Muon-catalyzed fusion cycle rate data for \em dd fusion with F = 3/2 state.
+ */
+MucfCycleRate dd_3_over_2_cycle_data()
+{
+    MucfCycleRate result;
+    result.type = CycleTableType::deuterium_deuterium;
+    result.spin_state = units::HalfSpinInt{3};  // F = 3/2
+
+    result.rate.interpolation.type = InterpolationType::cubic_spline;
+    // Temperature [K]
+    result.rate.x = {
+        0.0,
+        5.4317995646105715,
+        8.023769006655044,
+        15.563345056011997,
+        21.23403876489806,
+        24.19235525928977,
+        30.970416194589205,
+        36.77176514592074,
+        40.6299375993124,
+        48.3052926261127,
+        54.996250066448454,
+        65.48389514020138,
+        83.57516295078466,
+        104.55813870078731,
+        127.50862868996978,
+        141.88390769335217,
+        157.2358986807023,
+        176.44862398894298,
+        208.18055529758965,
+        248.59273415949275,
+        289.00491302139596,
+        347.6774429488119,
+        437.0821353934164,
+        532.1831402218961,
+        626.2766906564209,
+        748.1620201862349,
+        868.0888020131156,
+        992.7877025236317,
+    };
+    // Mean cycle rate [1/s]
+    result.rate.y = {
+        0.0,
+        3.002924371750031,
+        3.670220083632165,
+        3.976719989803767,
+        4.169019528475477,
+        3.9765902952616323,
+        3.820175795347777,
+        3.705860143810432,
+        3.6516939388136374,
+        3.6395545296699483,
+        3.6875501565621276,
+        3.8256690797323616,
+        4.119986345246231,
+        4.378188003927343,
+        4.519125621813277,
+        4.533939620625851,
+        4.506654771061913,
+        4.419191653948062,
+        4.202281855381299,
+        3.8650011112100278,
+        3.5277203670387562,
+        3.0879607598755188,
+        2.527498765500099,
+        2.0992156842652023,
+        1.785176202044242,
+        1.500779127953113,
+        1.312603875721246,
+        1.1754591026674808,
+    };
+
+    CELER_ENSURE(result);
+    return result;
+}
+
+//---------------------------------------------------------------------------//
+// DT cycle rate data
+//---------------------------------------------------------------------------//
+/*!
+ * Muon-catalyzed fusion cycle rate data for \em dt fusion with F = 0 state.
+ */
+MucfCycleRate dd_0_cycle_data()
+{
+    MucfCycleRate result;
+    result.type = CycleTableType::deuterium_deuterium;
+    result.spin_state = units::HalfSpinInt{0};  // F = 0
+
+    result.rate.interpolation.type = InterpolationType::cubic_spline;
+    // Temperature [K]
+    result.rate.x = {
+        9.719,
+        61.224,
+        113.677,
+        149.956,
+        203.175,
+        301.385,
+        407.915,
+        491.843,
+        609.109,
+        731.772,
+        859.935,
+        1041.085,
+        1205.668,
+        1370.401,
+        1487.748,
+    };
+    // Mean cycle rate [1/s]
+    result.rate.y = {
+        0.073e8,
+        1.614e8,
+        2.294e8,
+        2.404e8,
+        2.386e8,
+        2.195e8,
+        2.077e8,
+        2.143e8,
+        2.447e8,
+        2.934e8,
+        3.514e8,
+        4.286e8,
+        4.847e8,
+        5.271e8,
+        5.502e8,
+    };
+
+    CELER_ENSURE(result);
+    return result;
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Muon-catalyzed fusion cycle rate data for \em dt fusion with F = 0 state.
+ */
+MucfCycleRate dt_0_cycle_data()
+{
+    MucfCycleRate result;
+    result.type = CycleTableType::deuterium_tritium;
+    result.spin_state = units::HalfSpinInt{0};  // F = 0
+
+    result.rate.interpolation.type = InterpolationType::cubic_spline;
+    // Temperature [K]
+    result.rate.x = {
+        2.800,
+        40.599,
+        90.979,
+        155.357,
+        194.476,
+        233.474,
+        277.910,
+        319.385,
+        392.354,
+        462.421,
+        554.546,
+        669.030,
+        831.184,
+        991.254,
+        1175.567,
+        1349.074,
+        1500.362,
+    };
+    // Mean cycle rate [1/s]
+    result.rate.y = {
+        0.000e8,
+        0.001e8,
+        0.020e8,
+        0.039e8,
+        0.113e8,
+        0.296e8,
+        0.627e8,
+        1.104e8,
+        2.224e8,
+        3.435e8,
+        4.958e8,
+        6.518e8,
+        8.015e8,
+        8.860e8,
+        9.303e8,
+        9.388e8,
+        9.307e8,
+    };
+
+    CELER_ENSURE(result);
+    return result;
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Muon-catalyzed fusion cycle rate data for \em dt fusion with F = 0 state.
+ */
+MucfCycleRate hd_0_cycle_data()
+{
+    MucfCycleRate result;
+    result.type = CycleTableType::protium_deuterium;
+    result.spin_state = units::HalfSpinInt{0};  // F = 0
+
+    result.rate.interpolation.type = InterpolationType::cubic_spline;
+    // Temperature [K]
+    result.rate.x = {
+        1.400,   65.789,   141.358,  205.626,  237.603,  273.559,  308.014,
+        344.996, 409.061,  474.303,  555.983,  628.326,  714.709,  813.915,
+        926.103, 1058.264, 1183.899, 1262.237, 1368.676, 1500.476,
+    };
+    // Mean cycle rate [1/s]
+    result.rate.y = {
+        0.000e8,  0.010e8,  0.039e8,  0.159e8,  0.361e8,  0.765e8,  1.260e8,
+        2.003e8,  3.581e8,  5.361e8,  7.470e8,  9.158e8,  10.810e8, 12.260e8,
+        13.361e8, 14.124e8, 14.456e8, 14.512e8, 14.476e8, 14.295e8,
+    };
+
+    CELER_ENSURE(result);
+    return result;
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Muon-catalyzed fusion cycle rate data for \em dt fusion with F = 1 state.
+ */
+MucfCycleRate dd_1_cycle_data()
+{
+    MucfCycleRate result;
+    result.type = CycleTableType::deuterium_deuterium;
+    result.spin_state = units::HalfSpinInt{2};  // F = 1
+
+    result.rate.interpolation.type = InterpolationType::cubic_spline;
+    // Temperature [K]
+    result.rate.x = {
+        0.0,
+        40.97238,
+        86.18093,
+        154.00251,
+        216.16937,
+        309.43132,
+        391.40911,
+        484.73130,
+        592.21748,
+        733.66109,
+        900.53147,
+        1070.18645,
+        1251.11005,
+        1383.94917,
+        1501.23102,
+    };
+    // Mean cycle rate [1/s]
+    result.rate.y = {
+        0.0e8,
+        0.02590e8,
+        0.03788e8,
+        0.11774e8,
+        0.17032e8,
+        0.33172e8,
+        0.61735e8,
+        1.20516e8,
+        2.05375e8,
+        3.27240e8,
+        4.47630e8,
+        5.39123e8,
+        6.07188e8,
+        6.38302e8,
+        6.57099e8,
+    };
+
+    CELER_ENSURE(result);
+    return result;
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Muon-catalyzed fusion cycle rate data for \em dt fusion with F = 1 state.
+ */
+MucfCycleRate dt_1_cycle_data()
+{
+    MucfCycleRate result;
+    result.type = CycleTableType::deuterium_tritium;
+    result.spin_state = units::HalfSpinInt{2};  // F = 1
+
+    result.rate.interpolation.type = InterpolationType::cubic_spline;
+    // Temperature [K]
+    result.rate.x = {
+        0.0,
+        42.3405,
+        90.3261,
+        149.6118,
+        193.4004,
+        257.0266,
+        312.1948,
+        377.2730,
+        445.1485,
+        527.1292,
+        630.2549,
+        731.9188,
+        825.0831,
+        935.1662,
+        1079.0999,
+        1268.1724,
+        1502.3866,
+    };
+    // Mean cycle rate [1/s]
+    result.rate.y = {
+        0.0,
+        0.0254e8,
+        0.0505e8,
+        0.2126e8,
+        0.7603e8,
+        2.4341e8,
+        4.2458e8,
+        6.4969e8,
+        8.3905e8,
+        10.1734e8,
+        11.6117e8,
+        12.3352e8,
+        12.6055e8,
+        12.6414e8,
+        12.3869e8,
+        11.8140e8,
+        10.9640e8,
+    };
+
+    CELER_ENSURE(result);
+    return result;
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Muon-catalyzed fusion cycle rate data for \em dt fusion with F = 1 state.
+ */
+MucfCycleRate hd_1_cycle_data()
+{
+    MucfCycleRate result;
+    result.type = CycleTableType::protium_deuterium;
+    result.spin_state = units::HalfSpinInt{2};  // F = 1
+
+    result.rate.interpolation.type = InterpolationType::cubic_spline;
+    // Temperature [K]
+    result.rate.x = {
+        0.0,        57.92486,   108.78618,  151.19847, 182.34991, 197.94409,
+        227.75862,  274.69662,  320.26078,  380.06082, 442.66685, 515.12296,
+        595.97369,  690.86987,  795.54230,  883.16703, 976.39400, 1089.35418,
+        1261.56836, 1405.53035, 1505.73924,
+    };
+    // Mean cycle rate [1/s]
+    result.rate.y = {
+        0.0,        0.02523e8,  0.05075e8,  0.26916e8,  0.77688e8,  1.16141e8,
+        2.20563e8,  4.45963e8,  6.98879e8,  10.28768e8, 13.44890e8, 16.33464e8,
+        18.64233e8, 20.37175e8, 21.30299e8, 21.56089e8, 21.47469e8, 21.07136e8,
+        20.10173e8, 19.14696e8, 18.48278e8,
+    };
+
+    CELER_ENSURE(result);
+    return result;
+}
+
+//---------------------------------------------------------------------------//
+// TT cycle rate data
+//---------------------------------------------------------------------------//
+/*!
+ * Muon-catalyzed fusion cycle rate data for \em tt fusion with F = 1/2 state.
+ *
+ * The tt fusion cycle rate data is currently a constant value at 2.8e6 1/s.
+ */
+MucfCycleRate tt_1_over_2_cycle_data()
+{
+    MucfCycleRate result;
+    result.type = CycleTableType::tritium_tritium;
+    result.spin_state = units::HalfSpinInt{1};  // F = 1/2
+
+    result.rate.interpolation.type = InterpolationType::linear;
+    // Temperature [K]
+    result.rate.x = {
+        0.0,
+        std::numeric_limits<double>::max(),
+    };
+    // Mean cycle rate [1/s]
+    result.rate.y = {2.8e6, 2.8e6};
+
+    CELER_ENSURE(result);
+    return result;
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Cycle rate data for muon-catalyzed fusion.
+ *
+ * Data is extracted from \citet{faifman-mucfformation-1996,
+ * https://doi.org/10.1007/BF02227621}.
+ *
+ * \todo Use native units.
  */
 std::vector<MucfCycleRate> mucf_cycle_rates()
 {
-    //! \todo Implement
-    return {};
+    std::vector<MucfCycleRate> result;
+    // DD fusion
+    result.push_back(dd_1_over_2_cycle_data());  // F = 1/2
+    result.push_back(dd_3_over_2_cycle_data());  // F = 3/2
+    // DT fusion
+    // F = 0
+    result.push_back(hd_0_cycle_data());
+    result.push_back(dd_0_cycle_data());
+    result.push_back(dt_0_cycle_data());
+    // F = 1
+    result.push_back(hd_1_cycle_data());
+    result.push_back(dd_1_cycle_data());
+    result.push_back(dt_1_cycle_data());
+    // TT fusion
+    result.push_back(tt_1_over_2_cycle_data());  // F = 1/2
+
+    return result;
 }
 
 //---------------------------------------------------------------------------//
@@ -143,20 +610,6 @@ MucfPhysics MucfPhysics::from_default()
     result.atom_transfer = mucf_atom_transfer_rates();
     result.atom_spin_flip = mucf_atom_spin_flip_rates();
 
-    // Temporary test dummy data to verify correct import
-    {
-        MucfCycleRate dt_cycle;
-        dt_cycle.molecule = MucfMuonicMolecule::deuterium_tritium;
-
-        dt_cycle.spin_label = "F=0";
-        dt_cycle.rate = Grid::from_constant(2.0);
-        result.cycle_rates.push_back(dt_cycle);
-
-        dt_cycle.spin_label = "F=1";
-        dt_cycle.rate = Grid::from_constant(3.0);
-        result.cycle_rates.push_back(dt_cycle);
-    }
-
     return result;
 }
 
diff --git a/src/celeritas/inp/MucfPhysics.hh b/src/celeritas/inp/MucfPhysics.hh
index 401e58b525..a1925af47e 100644
--- a/src/celeritas/inp/MucfPhysics.hh
+++ b/src/celeritas/inp/MucfPhysics.hh
@@ -11,6 +11,7 @@
 
 #include "corecel/inp/Grid.hh"
 #include "corecel/math/Quantity.hh"
+#include "celeritas/Quantities.hh"
 #include "celeritas/UnitTypes.hh"
 #include "celeritas/mucf/Types.hh"
 
@@ -21,8 +22,6 @@ namespace inp
 //---------------------------------------------------------------------------//
 /*!
  * Muon-catalyzed fusion scalars.
- *
- * Default values are the same used by Acceleron.
  */
 struct MucfScalars
 {
@@ -47,6 +46,20 @@ struct MucfScalars
     static MucfScalars from_default();
 };
 
+//---------------------------------------------------------------------------//
+/*!
+ * Safely access cycle rate tables.
+ */
+enum class CycleTableType
+{
+    protium_deuterium,
+    protium_tritium,
+    deuterium_deuterium,
+    deuterium_tritium,
+    tritium_tritium,
+    size_
+};
+
 //---------------------------------------------------------------------------//
 /*!
  * Muon-catalyzed fusion mean cycle rate data.
@@ -65,15 +78,14 @@ struct MucfScalars
  */
 struct MucfCycleRate
 {
-    MucfMuonicMolecule molecule;
+    CycleTableType type;  //!< Cycle table type
     Grid rate;  //!< x = temperature [K], y = mean cycle rate [1/time]
-    std::string spin_label;
+    units::HalfSpinInt spin_state;  // Spin state F, in unit of hbar/2
 
     //! True if data is assigned
     explicit operator bool() const
     {
-        return molecule < MucfMuonicMolecule::size_ && rate
-               && !spin_label.empty();
+        return type < CycleTableType::size_ && rate;
     }
 };
 
diff --git a/src/celeritas/mucf/Types.hh b/src/celeritas/mucf/Types.hh
index 1fd0d46d3f..6580ea7ce5 100644
--- a/src/celeritas/mucf/Types.hh
+++ b/src/celeritas/mucf/Types.hh
@@ -68,6 +68,7 @@ enum class MucfIsoprotologueMolecule
     protium_protium,
     protium_deuterium,
     protium_tritium,
+    deuterium_deuterium,
     deuterium_tritium,
     tritium_tritium,
     size_
diff --git a/src/celeritas/mucf/data/DTMixMucfData.hh b/src/celeritas/mucf/data/DTMixMucfData.hh
index 3586060459..ff6aeb7c66 100644
--- a/src/celeritas/mucf/data/DTMixMucfData.hh
+++ b/src/celeritas/mucf/data/DTMixMucfData.hh
@@ -50,8 +50,8 @@ struct MucfParticleIds
     CELER_FUNCTION explicit operator bool() const
     {
         return mu_minus && proton && triton && neutron && alpha && he3
-               && muonic_hydrogen && muonic_alpha && muonic_triton
-               && muonic_he3;
+               && muonic_hydrogen && muonic_deuteron && muonic_triton
+               && muonic_alpha && muonic_he3;
     }
 };
 
@@ -89,8 +89,9 @@ struct MucfParticleMasses
                && triton > zero_quantity() && neutron > zero_quantity()
                && alpha > zero_quantity() && he3 > zero_quantity()
                && muonic_hydrogen > zero_quantity()
-               && muonic_alpha > zero_quantity()
+               && muonic_deuteron > zero_quantity()
                && muonic_triton > zero_quantity()
+               && muonic_alpha > zero_quantity()
                && muonic_he3 > zero_quantity();
     }
 };
@@ -134,14 +135,9 @@ struct DTMixMucfData
     //! Check whether the data are assigned
     explicit CELER_FUNCTION operator bool() const
     {
-        return true;
-#if 0
-        // Re-enable once full data assignment is implemented
         return particle_ids && particle_masses && muon_energy_cdf
-               && !mucfmatid_to_matid.empty() && !cycle_times.empty()
-               && (mucfmatid_to_matid.size() == isotopic_fractions.size())
-               && (mucfmatid_to_matid.size() == cycle_times.size());
-#endif
+               && !reals.empty() && !mucfmatid_to_matid.empty()
+               && !isotopic_fractions.empty() && !cycle_times.empty();
     }
 
     //! Assign from another set of data
@@ -153,8 +149,8 @@ struct DTMixMucfData
         //! \todo Finish implementation
         this->particle_ids = other.particle_ids;
         this->particle_masses = other.particle_masses;
-        this->reals = other.reals;
         this->muon_energy_cdf = other.muon_energy_cdf;
+        this->reals = other.reals;
         this->mucfmatid_to_matid = other.mucfmatid_to_matid;
         this->isotopic_fractions = other.isotopic_fractions;
         this->cycle_times = other.cycle_times;
diff --git a/src/celeritas/mucf/interactor/DTMucfInteractor.hh b/src/celeritas/mucf/interactor/DTMucfInteractor.hh
index ac6c5dec9f..39932258a8 100644
--- a/src/celeritas/mucf/interactor/DTMucfInteractor.hh
+++ b/src/celeritas/mucf/interactor/DTMucfInteractor.hh
@@ -59,7 +59,6 @@ class DTMucfInteractor
     Channel channel_{Channel::size_};
     // Allocate space for secondary particles
     StackAllocator<Secondary>& allocate_;
-
     // Number of secondaries per channel
     EnumArray<Channel, size_type> num_secondaries_{
         3,  // alpha_muon_neutron
diff --git a/src/celeritas/mucf/model/DTMixMucfModel.cc b/src/celeritas/mucf/model/DTMixMucfModel.cc
index b8e57efe60..d32706e1a9 100644
--- a/src/celeritas/mucf/model/DTMixMucfModel.cc
+++ b/src/celeritas/mucf/model/DTMixMucfModel.cc
@@ -50,19 +50,16 @@ from_params(ParticleParams const& particles)
     }
 
     MP_ADD(mu_minus);
-    MP_ADD(neutron);
     MP_ADD(proton);
+    MP_ADD(neutron);
+    MP_ADD(triton);
     MP_ADD(alpha);
     MP_ADD(he3);
+    MP_ADD(muonic_hydrogen);
     MP_ADD(muonic_deuteron);
     MP_ADD(muonic_triton);
     MP_ADD(muonic_alpha);
-
-    //! \todo Decide whether to implement these PDGs in PDGNumber.hh
-#if 0
-    MP_ADD(muonic_hydrogen);
     MP_ADD(muonic_he3);
-#endif
 
     CELER_VALIDATE(missing.empty(),
                    << "missing particles required for muon-catalyzed fusion: "
@@ -121,7 +118,7 @@ DTMixMucfModel::DTMixMucfModel(ActionId id,
     host_data.muon_energy_cdf = build_grid_record(inp_data.muon_energy_cdf);
 
     // Calculate and cache quantities for all materials with dt mixtures
-    detail::MucfMaterialInserter insert(&host_data);
+    detail::MucfMaterialInserter insert(&host_data, inp_data);
     for (auto const& matid : range(materials.num_materials()))
     {
         auto const& mat_view = materials.get(PhysMatId{matid});
diff --git a/src/celeritas/mucf/model/detail/EquilibrateDensitiesCalculator.cc b/src/celeritas/mucf/model/detail/EquilibrateDensitiesCalculator.cc
new file mode 100644
index 0000000000..0eae18f1b6
--- /dev/null
+++ b/src/celeritas/mucf/model/detail/EquilibrateDensitiesCalculator.cc
@@ -0,0 +1,246 @@
+//------------------------------- -*- C++ -*- -------------------------------//
+// Copyright Celeritas contributors: see top-level COPYRIGHT file for details
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file celeritas/mucf/model/detail/EquilibrateDensitiesCalculator.cc
+//---------------------------------------------------------------------------//
+#include "EquilibrateDensitiesCalculator.hh"
+
+#include <cmath>
+
+#include "corecel/cont/Range.hh"
+#include "corecel/io/Logger.hh"
+#include "corecel/math/Algorithms.hh"
+
+namespace celeritas
+{
+namespace detail
+{
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with material information.
+ */
+EquilibrateDensitiesCalculator::EquilibrateDensitiesCalculator(
+    LhdArray const& lhd_densities, real_type const temperature)
+    : lhd_densities_(lhd_densities), temperature_(temperature)
+{
+    CELER_EXPECT(temperature > 0);
+    CELER_EXPECT(lhd_densities[MucfIsotope::protium]
+                     + lhd_densities[MucfIsotope::deuterium]
+                     + lhd_densities[MucfIsotope::tritium]
+                 > 0);
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Return equilibrated isoprotologue values.
+ */
+EquilibrateDensitiesCalculator::EquilibriumArray
+EquilibrateDensitiesCalculator::operator()()
+{
+    using Iso = MucfIsotope;
+    using IsoProt = MucfIsoprotologueMolecule;
+
+    real_type const total_density = lhd_densities_[Iso::protium]
+                                    + lhd_densities_[Iso::deuterium]
+                                    + lhd_densities_[Iso::tritium];
+    CELER_ASSERT(total_density > 0);
+    real_type const inv_tot_density = real_type{1} / total_density;
+
+    // Cache equilibrium constants for this temperature for the while loop
+    real_type const k_hd = this->calc_hd_equilibrium_constant();
+    real_type const k_ht = this->calc_ht_equilibrium_constant();
+    real_type const k_dt = this->calc_dt_equilibrium_constant();
+
+    // Initialize result and set homonuclear molecules values
+    EquilibriumArray result;
+    result[IsoProt::protium_protium] = lhd_densities_[Iso::protium]
+                                       * inv_tot_density;
+    result[IsoProt::deuterium_deuterium] = lhd_densities_[Iso::deuterium]
+                                           * inv_tot_density;
+    result[IsoProt::tritium_tritium] = lhd_densities_[Iso::tritium]
+                                       * inv_tot_density;
+
+    EquilibriumArray previous_equilib_dens;
+    real_type iter_diff{0};
+    size_type iter{0};
+    do
+    {
+        iter++;
+        iter_diff = 0;
+        previous_equilib_dens = result;
+
+        // Equilibrate DT
+        this->equilibrate_pair(IsoProt::deuterium_deuterium,
+                               IsoProt::tritium_tritium,
+                               IsoProt::deuterium_tritium,
+                               k_dt,
+                               result);
+        // Equilibrate HT
+        this->equilibrate_pair(IsoProt::protium_protium,
+                               IsoProt::tritium_tritium,
+                               IsoProt::protium_tritium,
+                               k_ht,
+                               result);
+        // Equilibrate HD
+        this->equilibrate_pair(IsoProt::protium_protium,
+                               IsoProt::deuterium_deuterium,
+                               IsoProt::protium_deuterium,
+                               k_hd,
+                               result);
+
+        for (auto i : range(MucfIsoprotologueMolecule::size_))
+        {
+            // Calculate difference between current and previous densities
+            real_type diff = std::abs(result[i] - previous_equilib_dens[i]);
+            if (diff > iter_diff)
+            {
+                // Select maximum difference for convergence check
+                iter_diff = diff;
+            }
+        }
+    } while ((iter_diff > this->convergence_err())
+             && (iter < this->max_iterations()));
+
+    if (iter == this->max_iterations())
+    {
+        CELER_LOG(warning) << "Equilibration did not converge after "
+                           << max_iterations()
+                           << " iterations. Current error is " << iter_diff;
+    }
+
+    for (auto& val : result)
+    {
+        val *= total_density;
+    }
+
+    return result;
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Calculate equilibrium constant for the
+ * \f$ H_2 + D_2 \rightleftharpoons 2HD \f$ reaction.
+ */
+real_type EquilibrateDensitiesCalculator::calc_hd_equilibrium_constant()
+{
+    real_type result;
+
+    if (temperature_ < 30)
+    {
+        result = 6.785 * exp(-654.3 / (r_gas_constant_ * temperature_));
+    }
+    else
+    {
+        real_type const c_hd = r_gas_constant_ * 30 * (log(4) - log(0.49));
+        result = (4.0 * exp(-c_hd / (r_gas_constant_ * temperature_)));
+    }
+
+    return result;
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Calculate equilibrium constant for the
+ * \f$ H_2 + T_2 \rightleftharpoons 2HT \f$ reaction.
+ */
+real_type EquilibrateDensitiesCalculator::calc_ht_equilibrium_constant()
+{
+    real_type result;
+
+    if (temperature_ < 30)
+    {
+        result = 10.22 * exp(-1423 / (r_gas_constant_ * temperature_));
+    }
+    else
+    {
+        real_type const c_ht = r_gas_constant_ * 30 * (log(4) - log(0.034));
+        result = (4.0 * exp(-c_ht / (r_gas_constant_ * temperature_)));
+    }
+
+    return result;
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Calculate equilibrium constant for the
+ * \f$ D_2 + T_2 \rightleftharpoons 2DT \f$ reaction.
+ */
+real_type EquilibrateDensitiesCalculator::calc_dt_equilibrium_constant()
+{
+    real_type result;
+
+    if (temperature_ < 15)
+    {
+        result = 5.924 * exp(-168.3 / (r_gas_constant_ * temperature_));
+    }
+    else if (temperature_ < 30)
+    {
+        result = 2.995 * exp(-89.96 / (r_gas_constant_ * temperature_));
+    }
+    else if (temperature_ < 100)
+    {
+        real_type const c_dt = r_gas_constant_ * 30 * (log(4) - log(2.09));
+        result = 4.0 * exp(-c_dt / (r_gas_constant_ * temperature_));
+    }
+    else
+    {
+        real_type const c_dt = r_gas_constant_ * 100 * (log(4) - log(3.29));
+        result = 4.0 * exp(-c_dt / (r_gas_constant_ * temperature_));
+    }
+
+    return result;
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Equilibrate a pair of isotopes and return the new density.
+ *
+ * Since there are 3 isotopes (H, D, and T), and 6 molecular combinations, the
+ * equilibrium cannot be solved at once and has to be done iteratively for each
+ * pair until a convergence criterion is met.
+ *
+ * Therefore, this function takes 2 isotope combinations (e.g. DD, TT, and DT),
+ * the equilibrium constant for this temperature, and calculates how much of
+ * the homonuclear molecules (e.g. DD and TT) convert to the heteronuclear
+ * molecule (e.g. DT).
+ *
+ * The new densities are written into the input array.
+ */
+void EquilibrateDensitiesCalculator::equilibrate_pair(
+    MucfIsoprotologueMolecule molecule_aa,
+    MucfIsoprotologueMolecule molecule_bb,
+    MucfIsoprotologueMolecule molecule_ab,
+    real_type eq_constant_ab,
+    EquilibriumArray& input)
+{
+    CELER_EXPECT(molecule_aa < MucfIsoprotologueMolecule::size_);
+    CELER_EXPECT(molecule_bb < MucfIsoprotologueMolecule::size_);
+    CELER_EXPECT(molecule_ab < MucfIsoprotologueMolecule::size_);
+    CELER_EXPECT(eq_constant_ab > 0);
+
+    auto const& dens_aa = input[molecule_aa];
+    auto const& dens_bb = input[molecule_bb];
+    auto const& dens_ab = input[molecule_ab];
+
+    // AA + AB / 2
+    real_type const mix_a = dens_aa + dens_ab * real_type{0.5};
+    // BB + AB / 2
+    real_type const mix_b = dens_bb + dens_ab * real_type{0.5};
+
+    real_type sigma
+        = ((mix_a + mix_b)
+           - std::sqrt(ipow<2>(mix_a - mix_b)
+                       + 16 * mix_a * mix_b
+                             / (eq_constant_ab - this->convergence_err())))
+          / (2 * (1 - 4 / (eq_constant_ab - this->convergence_err())));
+
+    // Write new density into the equilibrium array
+    input[molecule_aa] = mix_a - sigma;
+    input[molecule_ab] = 2 * sigma;
+    input[molecule_bb] = mix_b - sigma;
+}
+
+//---------------------------------------------------------------------------//
+}  // namespace detail
+}  // namespace celeritas
diff --git a/src/celeritas/mucf/model/detail/EquilibrateDensitiesCalculator.hh b/src/celeritas/mucf/model/detail/EquilibrateDensitiesCalculator.hh
new file mode 100644
index 0000000000..99f21aacb7
--- /dev/null
+++ b/src/celeritas/mucf/model/detail/EquilibrateDensitiesCalculator.hh
@@ -0,0 +1,83 @@
+//------------------------------- -*- C++ -*- -------------------------------//
+// Copyright Celeritas contributors: see top-level COPYRIGHT file for details
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file celeritas/mucf/model/detail/EquilibrateDensitiesCalculator.hh
+//---------------------------------------------------------------------------//
+#pragma once
+
+#include "corecel/cont/EnumArray.hh"
+#include "celeritas/Constants.hh"
+#include "celeritas/mucf/Types.hh"
+
+namespace celeritas
+{
+namespace detail
+{
+//---------------------------------------------------------------------------//
+/*!
+ * Calculate dt mixture densities after reaching thermodynamic
+ * equilibrium based isotopic fraction, density, and material temperature.
+ *
+ * Based on the theory from https://www.osti.gov/biblio/6205719.
+ *
+ * The equilibrated densities are needed to correctly calculate the cycle time
+ * of dd, dt, and tt fusion cycles.
+ *
+ * \sa MucfMaterialInserter
+ */
+class EquilibrateDensitiesCalculator
+{
+  public:
+    //!@{
+    //! \name Type aliases
+    using LhdArray = EnumArray<MucfIsotope, real_type>;
+    using EquilibriumArray = EnumArray<MucfIsoprotologueMolecule, real_type>;
+    //!@}
+
+    // Construct with material information
+    EquilibrateDensitiesCalculator(LhdArray const& lhd_densities,
+                                   real_type const temperature);
+
+    // Calculate equilibrated isoprotologue densities
+    EquilibriumArray operator()();
+
+  private:
+    //// DATA ////
+
+    LhdArray lhd_densities_;  //!< Number densities in units of LHD
+    real_type temperature_;  //!< Material temperature [K]
+    real_type const r_gas_constant_ = constants::k_boltzmann.value()
+                                      * constants::na_avogadro.value();
+
+    //// HELPER FUNCTIONS ////
+
+    // {
+    // Convergence limit parameters
+    // Acceptance error between current and previous equilibrium iteration
+    static real_type constexpr convergence_err() { return 1e-6; }
+
+    // Maximum number of iterations to reach convergence
+    static size_type constexpr max_iterations() { return 1000; }
+    // }
+
+    // Calculate equilibrium constant: \f$ H_2 + D_2 \rightleftharpoons 2HD \f$
+    real_type calc_hd_equilibrium_constant();
+
+    // Calculate equilibrium constant: \f$ H_2 + T_2 \rightleftharpoons 2HT \f$
+    real_type calc_ht_equilibrium_constant();
+
+    // Calculate equilibrium constant: \f$ D_2 + T_2 \rightleftharpoons 2DT \f$
+    real_type calc_dt_equilibrium_constant();
+
+    // Equilibrate a pair of isotopes and write the new density in the array
+    void equilibrate_pair(MucfIsoprotologueMolecule molecule_aa,
+                          MucfIsoprotologueMolecule molecule_bb,
+                          MucfIsoprotologueMolecule molecule_ab,
+                          real_type eq_constant_ab,
+                          EquilibriumArray& input);
+};
+
+//---------------------------------------------------------------------------//
+}  // namespace detail
+}  // namespace celeritas
diff --git a/src/celeritas/mucf/model/detail/InterpolatorHelper.cc b/src/celeritas/mucf/model/detail/InterpolatorHelper.cc
new file mode 100644
index 0000000000..3d9beb6ead
--- /dev/null
+++ b/src/celeritas/mucf/model/detail/InterpolatorHelper.cc
@@ -0,0 +1,39 @@
+//------------------------------- -*- C++ -*- -------------------------------//
+// Copyright Celeritas contributors: see top-level COPYRIGHT file for details
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file celeritas/mucf/model/detail/InterpolatorHelper.cc
+//---------------------------------------------------------------------------//
+#include "InterpolatorHelper.hh"
+
+#include "corecel/Assert.hh"
+
+namespace celeritas
+{
+namespace detail
+{
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with grid input data.
+ */
+InterpolatorHelper::InterpolatorHelper(inp::Grid const& input)
+    : grid_record_(NonuniformGridBuilder(&reals_)(input))
+    , reals_ref_(reals_)
+    , interpolate_(grid_record_, reals_ref_)
+{
+    CELER_EXPECT(input);
+    CELER_ENSURE(grid_record_);
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Interpolate data at given point.
+ */
+real_type InterpolatorHelper::operator()(real_type value) const
+{
+    return interpolate_(value);
+}
+
+//---------------------------------------------------------------------------//
+}  // namespace detail
+}  // namespace celeritas
diff --git a/src/celeritas/mucf/model/detail/InterpolatorHelper.hh b/src/celeritas/mucf/model/detail/InterpolatorHelper.hh
new file mode 100644
index 0000000000..7f4e496049
--- /dev/null
+++ b/src/celeritas/mucf/model/detail/InterpolatorHelper.hh
@@ -0,0 +1,50 @@
+//------------------------------- -*- C++ -*- -------------------------------//
+// Copyright Celeritas contributors: see top-level COPYRIGHT file for details
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file celeritas/mucf/model/detail/InterpolatorHelper.hh
+//---------------------------------------------------------------------------//
+#pragma once
+
+#include "corecel/Macros.hh"
+#include "corecel/grid/NonuniformGridData.hh"
+#include "corecel/inp/Grid.hh"
+#include "celeritas/grid/NonuniformGridBuilder.hh"
+#include "celeritas/grid/NonuniformGridCalculator.hh"
+#include "celeritas/inp/MucfPhysics.hh"
+
+namespace celeritas
+{
+namespace detail
+{
+//---------------------------------------------------------------------------//
+/*!
+ * Host-only interpolator wrapper class for muCF input data.
+ *
+ * \sa MucfMaterialInserter
+ */
+class InterpolatorHelper
+{
+  public:
+    // Construct with grid input data
+    explicit InterpolatorHelper(inp::Grid const& input);
+    // Prevent copy and move
+    CELER_DELETE_COPY_MOVE(InterpolatorHelper);
+
+    // Interpolate data at given point
+    real_type operator()(real_type value) const;
+
+  private:
+    using Items = Collection<real_type, Ownership::value, MemSpace::host>;
+    using ItemsRef
+        = Collection<real_type, Ownership::const_reference, MemSpace::host>;
+
+    Items reals_;
+    NonuniformGridRecord grid_record_;
+    ItemsRef reals_ref_;
+    NonuniformGridCalculator interpolate_;
+};
+
+//---------------------------------------------------------------------------//
+}  // namespace detail
+}  // namespace celeritas
diff --git a/src/celeritas/mucf/model/detail/MucfMaterialInserter.cc b/src/celeritas/mucf/model/detail/MucfMaterialInserter.cc
index f859297ab9..282ee5ddac 100644
--- a/src/celeritas/mucf/model/detail/MucfMaterialInserter.cc
+++ b/src/celeritas/mucf/model/detail/MucfMaterialInserter.cc
@@ -16,23 +16,42 @@ namespace detail
 /*!
  * Construct with \c DTMixMucfModel model data.
  */
-MucfMaterialInserter::MucfMaterialInserter(HostVal<DTMixMucfData>* host_data)
+MucfMaterialInserter::MucfMaterialInserter(HostVal<DTMixMucfData>* host_data,
+                                           inp::MucfPhysics const& data)
     : mucfmatid_to_matid_(&host_data->mucfmatid_to_matid)
+    , isotopic_fractions_(&host_data->isotopic_fractions)
     , cycle_times_(&host_data->cycle_times)
+    , data_(data)
 {
-    CELER_EXPECT(host_data);
+    CELER_EXPECT(data_);
+
+    // Initialize interpolators for cycle time tables
+    for (auto const& cycle_data : data_.cycle_rates)
+    {
+        // Use emplace to avoid copy/move of InterpolatorHelper objects
+        interpolators_.emplace(
+            std::pair{cycle_data.type, cycle_data.spin_state}, cycle_data.rate);
+    }
 }
 
 //---------------------------------------------------------------------------//
 /*!
  * Insert material information if applicable.
- *
- * Calculates and caches material-dependent properties needed by the
- * \c DTMixMucfModel . If the material does not contain deuterium and/or
- * tritium the operator will return false.
  */
 bool MucfMaterialInserter::operator()(MaterialView const& material)
 {
+    using LhdArray = EquilibrateDensitiesCalculator::LhdArray;
+
+    MaterialFractionsArray isotopic_fractions;
+    CycleTimesArray cycle_times;
+    LhdArray lhd_densities{};
+
+    auto from_mass_number = [&](AtomicMassNumber mass) -> MucfIsotope {
+        auto it = mass_isotope_map_.find(mass);
+        return (it != mass_isotope_map_.end()) ? it->second
+                                               : MucfIsotope::size_;
+    };
+
     for (auto elcompid : range(material.num_elements()))
     {
         auto const& element_view
@@ -43,199 +62,181 @@ bool MucfMaterialInserter::operator()(MaterialView const& material)
             continue;
         }
 
-        // Found hydrogen; Check for d and t isotopes
-        IsotopeChecker has_isotope{false, false};
+        // Found hydrogen; Check isotopes
+        auto const elem_rel_abundance = material.elements()[elcompid].fraction;
         for (auto el_comp : range(element_view.num_isotopes()))
         {
             auto iso_view
                 = element_view.isotope_record(IsotopeComponentId{el_comp});
-            auto mass = iso_view.atomic_mass_number();
-            if (mass == AtomicMassNumber{1})
-            {
-                // Skip protium
-                continue;
-            }
-
-            auto const atom = this->from_mass_number(mass);
-            if (atom < MucfMuonicAtom::size_)
-            {
-                // Mark d or t isotope as found
-                has_isotope[atom] = true;
-            }
-        }
+            auto const atom = from_mass_number(iso_view.atomic_mass_number());
+            CELER_ASSERT(atom < MucfIsotope::size_);
 
-        if (!has_isotope[MucfMuonicAtom::deuterium]
-            && !has_isotope[MucfMuonicAtom::tritium])
-        {
-            // No deuterium or tritium found; skip material
-            return false;
-        }
+            auto const iso_frac = element_view.isotopes()[el_comp].fraction;
 
-        // Temporary data needed to calculate model data, such as cycle times
-        lhd_densities_ = this->calc_lhd_densities(element_view);
-        equilibrium_densities_ = this->calc_equilibrium_densities(element_view);
+            // Cache density for this hydrogen isotope
+            isotopic_fractions[atom] = iso_frac;
+            lhd_densities[atom]
+                = iso_frac
+                  * (elem_rel_abundance * material.number_density()
+                     / data_.scalars.liquid_hydrogen_density.value());
+        }
+    }
 
-        // Calculate and insert muCF material data into model data
-        mucfmatid_to_matid_.push_back(material.material_id());
-        cycle_times_.push_back(
-            this->calc_cycle_times(element_view, has_isotope));
-        //! \todo Store mean atom spin flip and transfer times
+    if (!lhd_densities[MucfIsotope::deuterium]
+        && !lhd_densities[MucfIsotope::tritium])
+    {
+        // No deuterium or tritium densities; skip material
+        return false;
     }
-    return true;
-}
 
-//---------------------------------------------------------------------------//
-/*!
- * Return \c MucfMuonicAtom from a given atomic mass number.
- */
-MucfMuonicAtom MucfMaterialInserter::from_mass_number(AtomicMassNumber mass)
-{
-    if (mass == AtomicMassNumber{2})
+    // Found d and/or t, calculate and insert data into collection
+
+    auto equilibrium_densities = EquilibrateDensitiesCalculator(
+        lhd_densities, material.temperature())();
+
+    if (lhd_densities[MucfIsotope::deuterium])
     {
-        return MucfMuonicAtom::deuterium;
+        cycle_times[MucfMuonicMolecule::deuterium_deuterium]
+            = this->calc_dd_cycle(equilibrium_densities,
+                                  material.temperature());
     }
-    if (mass == AtomicMassNumber{3})
+    if (lhd_densities[MucfIsotope::tritium])
     {
-        return MucfMuonicAtom::tritium;
+        cycle_times[MucfMuonicMolecule::tritium_tritium] = this->calc_tt_cycle(
+            equilibrium_densities, material.temperature());
+    }
+    if (lhd_densities[MucfIsotope::deuterium]
+        && lhd_densities[MucfIsotope::tritium])
+    {
+        cycle_times[MucfMuonicMolecule::deuterium_tritium]
+            = this->calc_dt_cycle(equilibrium_densities,
+                                  material.temperature());
     }
-    return MucfMuonicAtom::size_;
-}
 
-//---------------------------------------------------------------------------//
-/*!
- * Convert dt mixture densities to units of liquid hydrogen density.
- *
- * Used during cycle time calculations.
- */
-MucfMaterialInserter::LhdArray
-MucfMaterialInserter::calc_lhd_densities(ElementView const&)
-{
-    LhdArray result;
+    // Add muCF material to the model's host/device data
+    mucfmatid_to_matid_.push_back(material.material_id());
+    isotopic_fractions_.push_back(std::move(isotopic_fractions));
+    cycle_times_.push_back(std::move(cycle_times));
 
-    //! \todo Implement
+    //! \todo Store mean atom spin flip and transfer times
 
-    return result;
+    return true;
 }
 
 //---------------------------------------------------------------------------//
 /*!
- * Calculate dt mixture densities after reaching thermodynamical
- * equilibrium.
+ * Calculate dd muonic molecules cycle times.
  *
- * Used during cycle time calculations.
+ * F = 1/2 and F = 3/2 are the reactive spin states for dd fusion.
  */
-MucfMaterialInserter::EquilibriumArray
-MucfMaterialInserter::calc_equilibrium_densities(ElementView const&)
+MucfMaterialInserter::MoleculeCycles
+MucfMaterialInserter::calc_dd_cycle(EquilibriumArray const& eq_dens,
+                                    real_type const temperature)
 {
-    EquilibriumArray result;
+    using IsoProt = MucfIsoprotologueMolecule;
+    using CTT = inp::CycleTableType;
+    using units::HalfSpinInt;
 
-    //! \todo Implement
+    auto const& dd_dens = eq_dens[IsoProt::deuterium_deuterium];
 
-    return result;
-}
+    auto const& dd_1_over_2_interpolate
+        = this->interpolator(CTT::deuterium_deuterium, HalfSpinInt{1});
+    auto const& dd_3_over_2_interpolate
+        = this->interpolator(CTT::deuterium_deuterium, HalfSpinInt{3});
 
-//---------------------------------------------------------------------------//
-/*!
- * Calculate fusion mean cycle times.
- *
- * This is designed to work with the user's material definition being either:
- * - Single element, multiple isotopes (H element, with H, d, and t isotopes);
- * or
- * - Multiple elements, single isotope each (separate H, d, and t elements).
- */
-MucfMaterialInserter::CycleTimesArray
-MucfMaterialInserter::calc_cycle_times(ElementView const& element,
-                                       IsotopeChecker const& has_isotope)
-{
-    CycleTimesArray result;
-    for (auto el_comp : range(element.num_isotopes()))
-    {
-        auto iso_view = element.isotope_record(IsotopeComponentId{el_comp});
+    MoleculeCycles result;
+    result[0] = real_type{1}
+                / (dd_dens * dd_1_over_2_interpolate(temperature));  // F = 1/2
+    result[1] = real_type{1}
+                / (dd_dens * dd_3_over_2_interpolate(temperature));  // F = 3/2
 
-        // Select possible muonic atom based on the isotope/element mass number
-        auto atom = this->from_mass_number(iso_view.atomic_mass_number());
-        switch (atom)
-        {
-            // Calculate cycle times for dd molecules
-            case MucfMuonicAtom::deuterium: {
-                result[MucfMuonicMolecule::deuterium_deuterium]
-                    = this->calc_dd_cycle(element);
-                if (has_isotope[MucfMuonicAtom::tritium])
-                {
-                    // Calculate cycle times for dt molecules
-                    result[MucfMuonicMolecule::deuterium_tritium]
-                        = this->calc_dt_cycle(element);
-                }
-                break;
-            }
-            // Calculate cycle times for tt molecules
-            case MucfMuonicAtom::tritium: {
-                result[MucfMuonicMolecule::tritium_tritium]
-                    = this->calc_tt_cycle(element);
-                break;
-            }
-            default:
-                CELER_ASSERT_UNREACHABLE();
-        }
-    }
+    CELER_ENSURE(result[0] >= 0 && result[1] >= 0);
     return result;
 }
 
 //---------------------------------------------------------------------------//
 /*!
- * Calculate dd muonic molecules cycle times from material properties and grid
- * data.
+ * Calculate dt muonic molecules cycle times.
  *
- * Cycle times for dd molecules come from F = 0 and F = 1 spin states.
+ * F = 0 and F = 1 are the reactive spin states for dt fusion.
  */
 MucfMaterialInserter::MoleculeCycles
-MucfMaterialInserter::calc_dd_cycle(ElementView const&)
+MucfMaterialInserter::calc_dt_cycle(EquilibriumArray const& eq_dens,
+                                    real_type const temperature)
 {
+    CELER_EXPECT(temperature > 0);
+
+    using IsoProt = MucfIsoprotologueMolecule;
+    using CTT = inp::CycleTableType;
+    using units::HalfSpinInt;
+
+    auto const& dd_dens = eq_dens[IsoProt::deuterium_deuterium];
+    auto const& dt_dens = eq_dens[IsoProt::deuterium_tritium];
+    auto const& hd_dens = eq_dens[IsoProt::protium_deuterium];
+
+    // F = 0 interpolators
+    auto const& hd0_interpolate
+        = this->interpolator(CTT::protium_deuterium, HalfSpinInt{0});
+    auto const& dd0_interpolate
+        = this->interpolator(CTT::deuterium_deuterium, HalfSpinInt{0});
+    auto const& dt0_interpolate
+        = this->interpolator(CTT::deuterium_tritium, HalfSpinInt{0});
+    // F = 1 interpolators
+    auto const& hd1_interpolate
+        = this->interpolator(CTT::protium_deuterium, HalfSpinInt{2});
+    auto const& dd1_interpolate
+        = this->interpolator(CTT::deuterium_deuterium, HalfSpinInt{2});
+    auto const& dt1_interpolate
+        = this->interpolator(CTT::deuterium_tritium, HalfSpinInt{2});
+
+    // Interpolate over rates, store final cycle time (1/rate)
     MoleculeCycles result;
+    result[0] = real_type{1}
+                / (hd_dens * hd0_interpolate(temperature)
+                   + dd_dens * dd0_interpolate(temperature)
+                   + dt_dens * dt0_interpolate(temperature));  // F = 0
+    result[1] = real_type{1}
+                / (hd_dens * hd1_interpolate(temperature)
+                   + dd_dens * dd1_interpolate(temperature)
+                   + dt_dens * dt1_interpolate(temperature));  // F = 1
 
-    //! \todo Implement
-
-    // Reactive states are F = 0 and F = 1
     CELER_ENSURE(result[0] >= 0 && result[1] >= 0);
     return result;
 }
 
 //---------------------------------------------------------------------------//
 /*!
- * Calculate dt muonic molecules cycle times from material properties and grid
- * data.
+ * Calculate tt muonic molecules cycle times.
  *
- * Cycle times for dt molecules come from F = 1/2 and F = 3/2 spin states.
+ * F = 1/2 is the only reactive spin state for tt fusion.
  */
 MucfMaterialInserter::MoleculeCycles
-MucfMaterialInserter::calc_dt_cycle(ElementView const&)
+MucfMaterialInserter::calc_tt_cycle(EquilibriumArray const& eq_dens,
+                                    real_type const temperature)
 {
-    MoleculeCycles result;
+    using IsoProt = MucfIsoprotologueMolecule;
+    using CTT = inp::CycleTableType;
+    using units::HalfSpinInt;
 
-    //! \todo Implement
+    auto const& tt_dens = eq_dens[IsoProt::tritium_tritium];
+    auto const& tt_interpolate
+        = this->interpolator(CTT::tritium_tritium, HalfSpinInt{1});
 
-    // Reactive states are F = 1/2 and F = 3/2
-    CELER_ENSURE(result[0] >= 0 && result[1] >= 0);
+    MoleculeCycles result;
+    result[0] = real_type{1} / (tt_dens * tt_interpolate(temperature));
+
+    CELER_ENSURE(result[0] >= 0 && result[1] == 0);
     return result;
 }
 
 //---------------------------------------------------------------------------//
-/*!
- * Calculate tt muonic molecules cycle times from material properties and grid
- * data.
- *
- * Cycle times for tt molecules come only from the F = 1/2 spin state.
- */
-MucfMaterialInserter::MoleculeCycles
-MucfMaterialInserter::calc_tt_cycle(ElementView const&)
+InterpolatorHelper const&
+MucfMaterialInserter::interpolator(inp::CycleTableType type,
+                                   units::HalfSpinInt spin) const
 {
-    MoleculeCycles result;
-
-    //! \todo Implement
-
-    // Only F = 1/2 is reactive
-    CELER_ENSURE(result[0] >= 0 && result[1] == 0);
-    return result;
+    auto it = interpolators_.find({type, spin});
+    CELER_ASSERT(it != interpolators_.end());
+    return it->second;
 }
 
 //---------------------------------------------------------------------------//
diff --git a/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh b/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh
index ad7c307b9b..77a44496a3 100644
--- a/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh
+++ b/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh
@@ -6,11 +6,17 @@
 //---------------------------------------------------------------------------//
 #pragma once
 
+#include <map>
+
 #include "corecel/data/CollectionBuilder.hh"
+#include "celeritas/inp/MucfPhysics.hh"
 #include "celeritas/mat/MaterialView.hh"
 #include "celeritas/mucf/Types.hh"
 #include "celeritas/mucf/data/DTMixMucfData.hh"
 
+#include "EquilibrateDensitiesCalculator.hh"
+#include "InterpolatorHelper.hh"
+
 namespace celeritas
 {
 namespace detail
@@ -18,58 +24,70 @@ namespace detail
 //---------------------------------------------------------------------------//
 /*!
  * Helper class to calculate and insert muCF material-dependent data into
- * \c DTMixMucfData .
+ * \c DTMixMucfData . If the material does not contain deuterium and/or
+ * tritium the operator will return false.
+ *
+ * This is designed to work with the user's material definition being either:
+ * - Single element, multiple isotopes (H element, with H, d, and t isotopes);
+ * or
+ * - Multiple elements, single isotope each (separate H, d, and t elements).
+ *
+ * The \c inp:: data has cycle \em rate (\f$\lambda\f$) tables, while the
+ * host/device cached data is the cycle \em time \f$\tau = 1/\lambda\f$.
  */
 class MucfMaterialInserter
 {
   public:
     // Construct with muCF host data
-    explicit MucfMaterialInserter(HostVal<DTMixMucfData>* host_data);
+    explicit MucfMaterialInserter(HostVal<DTMixMucfData>* host_data,
+                                  inp::MucfPhysics const& data);
 
     // Insert material if it is a valid d-t mixture
     bool operator()(MaterialView const& material);
 
   private:
-    //// DATA ////
-
     using MoleculeCycles = Array<real_type, 2>;
     using CycleTimesArray = EnumArray<MucfMuonicMolecule, MoleculeCycles>;
-    using LhdArray = EnumArray<MucfMuonicAtom, real_type>;
-    using EquilibriumArray = EnumArray<MucfIsoprotologueMolecule, real_type>;
+    using EquilibriumArray = EquilibrateDensitiesCalculator::EquilibriumArray;
+    using MaterialFractionsArray = EnumArray<MucfIsotope, real_type>;
     using AtomicMassNumber = AtomicNumber;
-    using MucfIsotope = MucfMuonicAtom;
-    using IsotopeChecker = EnumArray<MucfIsotope, bool>;
+    using InterpolatorsMap
+        = std::map<std::pair<inp::CycleTableType, units::HalfSpinInt>,
+                   InterpolatorHelper>;
+
+    //// DATA ////
 
     // DTMixMucfModel host data references populated by operator()
     CollectionBuilder<PhysMatId, MemSpace::host, MuCfMatId> mucfmatid_to_matid_;
+    CollectionBuilder<MaterialFractionsArray, MemSpace::host, MuCfMatId>
+        isotopic_fractions_;
     CollectionBuilder<CycleTimesArray, MemSpace::host, MuCfMatId> cycle_times_;
-    // Temporary quantities needed for calculating the model data
-    LhdArray lhd_densities_;
-    EquilibriumArray equilibrium_densities_;
+    // Const data
+    std::map<AtomicMassNumber, MucfIsotope> const mass_isotope_map_{
+        {AtomicMassNumber{1}, MucfIsotope::protium},
+        {AtomicMassNumber{2}, MucfIsotope::deuterium},
+        {AtomicMassNumber{3}, MucfIsotope::tritium},
+    };
+    inp::MucfPhysics const& data_;
+    InterpolatorsMap interpolators_;
 
     //// HELPER FUNCTIONS ////
 
-    // Return muonic atom from given atomic mass number
-    MucfMuonicAtom from_mass_number(AtomicMassNumber mass);
-
-    // Calculate dt mixture densities in units of liquid hydrogen density
-    LhdArray calc_lhd_densities(ElementView const&);
-
-    // Calculate thermal equilibrium densities
-    EquilibriumArray calc_equilibrium_densities(ElementView const&);
-
-    // Calculate mean fusion cycle times for all reactive muonic molecules
-    CycleTimesArray calc_cycle_times(ElementView const& element,
-                                     IsotopeChecker const& has_isotope);
-
     // Calculate mean fusion cycle times for dd muonic molecules
-    Array<real_type, 2> calc_dd_cycle(ElementView const&);
+    Array<real_type, 2> calc_dd_cycle(EquilibriumArray const& eq_dens,
+                                      real_type const temperature);
 
     // Calculate mean fusion cycle times for dt muonic molecules
-    Array<real_type, 2> calc_dt_cycle(ElementView const&);
+    Array<real_type, 2> calc_dt_cycle(EquilibriumArray const& eq_dens,
+                                      real_type const temperature);
 
     // Calculate mean fusion cycle times for tt muonic molecules
-    Array<real_type, 2> calc_tt_cycle(ElementView const&);
+    Array<real_type, 2> calc_tt_cycle(EquilibriumArray const& eq_dens,
+                                      real_type const temperature);
+
+    // Get interpolator for given cycle type and spin
+    InterpolatorHelper const&
+    interpolator(inp::CycleTableType type, units::HalfSpinInt spin) const;
 };
 
 //---------------------------------------------------------------------------//
diff --git a/test/celeritas/CMakeLists.txt b/test/celeritas/CMakeLists.txt
index 5a21748cbb..2a9b6c4fac 100644
--- a/test/celeritas/CMakeLists.txt
+++ b/test/celeritas/CMakeLists.txt
@@ -153,6 +153,7 @@ celeritas_add_test(decay/MuDecay.test.cc ${_needs_double})
 
 celeritas_add_test(mucf/DDMucfInteractor.test.cc)
 celeritas_add_test(mucf/DTMucfInteractor.test.cc)
+celeritas_add_test(mucf/DTMixMucfModel.test.cc)
 celeritas_add_test(mucf/MuonicAtomSelector.test.cc)
 
 #-----------------------------------------------------------------------------#
diff --git a/test/celeritas/ext/GeantImporter.test.cc b/test/celeritas/ext/GeantImporter.test.cc
index a42578192a..2d5128d98b 100644
--- a/test/celeritas/ext/GeantImporter.test.cc
+++ b/test/celeritas/ext/GeantImporter.test.cc
@@ -2152,22 +2152,7 @@ TEST_F(MucfBox, static_data)
     EXPECT_SOFT_EQ(0.55157437567861023, average(mucf.muon_energy_cdf.x));
     EXPECT_SOFT_EQ(11.250286274435437, average(mucf.muon_energy_cdf.y));
 
-    // Dummy data
-    auto const& cycle_f0 = mucf.cycle_rates[0];
-    static double const expected_cycle_rate_f0_y[] = {2, 2};
-    EXPECT_TRUE(cycle_f0);
-    EXPECT_EQ(cycle_f0.molecule, MucfMuonicMolecule::deuterium_tritium);
-    EXPECT_EQ("F=0", cycle_f0.spin_label);
-    EXPECT_EQ(2, cycle_f0.rate.x.size());
-    EXPECT_VEC_EQ(expected_cycle_rate_f0_y, cycle_f0.rate.y);
-
-    auto const& cycle_f1 = mucf.cycle_rates[1];
-    static double const expected_cycle_rate_f1_y[] = {3, 3};
-    EXPECT_TRUE(cycle_f1);
-    EXPECT_EQ(cycle_f1.molecule, MucfMuonicMolecule::deuterium_tritium);
-    EXPECT_EQ("F=1", cycle_f1.spin_label);
-    EXPECT_EQ(2, cycle_f1.rate.x.size());
-    EXPECT_VEC_EQ(expected_cycle_rate_f1_y, cycle_f1.rate.y);
+    //! \todo Add real cycle rate data test
 
     EXPECT_TRUE(mucf.atom_transfer.empty());
     EXPECT_TRUE(mucf.atom_spin_flip.empty());
diff --git a/test/celeritas/mucf/DDMucfInteractor.test.cc b/test/celeritas/mucf/DDMucfInteractor.test.cc
index 283534290e..522e656da0 100644
--- a/test/celeritas/mucf/DDMucfInteractor.test.cc
+++ b/test/celeritas/mucf/DDMucfInteractor.test.cc
@@ -9,7 +9,6 @@
 #include "corecel/cont/Range.hh"
 #include "corecel/math/ArrayUtils.hh"
 #include "celeritas/Quantities.hh"
-#include "celeritas/grid/NonuniformGridBuilder.hh"
 #include "celeritas/inp/MucfPhysics.hh"
 
 #include "MucfInteractorHostTestBase.hh"
@@ -33,39 +32,8 @@ class DDMucfInteractorTest : public MucfInteractorHostTestBase
   protected:
     void SetUp() override
     {
-        auto const& params = *this->particle_params();
-        this->set_material("hdt_fuel");
-
-        HostVal<DTMixMucfData> host_data;
-
-        // Set up particle IDs
-        host_data.particle_ids.mu_minus = params.find(pdg::mu_minus());
-        host_data.particle_ids.neutron = params.find(pdg::neutron());
-        host_data.particle_ids.proton = params.find(pdg::proton());
-        host_data.particle_ids.he3 = params.find(pdg::he3());
-        host_data.particle_ids.triton = params.find(pdg::triton());
-        host_data.particle_ids.muonic_he3 = params.find(pdg::muonic_he3());
-
-        // Set up particle masses
-        host_data.particle_masses.mu_minus
-            = params.get(host_data.particle_ids.mu_minus).mass();
-        host_data.particle_masses.neutron
-            = params.get(host_data.particle_ids.neutron).mass();
-        host_data.particle_masses.proton
-            = params.get(host_data.particle_ids.proton).mass();
-        host_data.particle_masses.he3
-            = params.get(host_data.particle_ids.he3).mass();
-        host_data.particle_masses.triton
-            = params.get(host_data.particle_ids.triton).mass();
-        host_data.particle_masses.muonic_he3
-            = params.get(host_data.particle_ids.muonic_he3).mass();
-
-        // Set up muon energy CDF
-        auto const inp_data = inp::MucfPhysics::from_default();
-        NonuniformGridBuilder build_grid_record{&host_data.reals};
-        host_data.muon_energy_cdf = build_grid_record(inp_data.muon_energy_cdf);
-
         // Construct collection
+        auto host_data = this->make_host_data();
         data_ = ParamsDataStore<DTMixMucfData>{std::move(host_data)};
 
         // At-rest muon primary
diff --git a/test/celeritas/mucf/DTMixMucfModel.test.cc b/test/celeritas/mucf/DTMixMucfModel.test.cc
new file mode 100644
index 0000000000..37c3886dfc
--- /dev/null
+++ b/test/celeritas/mucf/DTMixMucfModel.test.cc
@@ -0,0 +1,115 @@
+//------------------------------- -*- C++ -*- -------------------------------//
+// Copyright Celeritas contributors: see top-level COPYRIGHT file for details
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file celeritas/mucf/DTMixMucfModel.test.cc
+//---------------------------------------------------------------------------//
+
+#include "celeritas/mucf/model/DTMixMucfModel.hh"
+
+#include "MucfInteractorHostTestBase.hh"
+#include "celeritas_test.hh"
+
+namespace celeritas
+{
+namespace test
+{
+//---------------------------------------------------------------------------//
+
+class DTMixMucfModelTest : public MucfInteractorHostTestBase
+{
+  protected:
+    void SetUp() override
+    {
+        particles_ = this->particle_params();
+        this->set_material("hdt_fuel");
+        materials_ = this->material_params();
+    }
+
+  protected:
+    std::shared_ptr<ParticleParams const> particles_;
+    std::shared_ptr<MaterialParams const> materials_;
+};
+
+//---------------------------------------------------------------------------//
+TEST_F(DTMixMucfModelTest, data)
+{
+    using Molecule = MucfMuonicMolecule;
+
+    auto model = DTMixMucfModel(ActionId{0}, *particles_, *materials_);
+    auto const& data = model.host_ref();
+    auto const& pids = data.particle_ids;
+    auto const& masses = data.particle_masses;
+
+#define EXPECT_PDG_EQ(MEMBER) \
+    EXPECT_EQ(pdg::MEMBER().get(), particles_->id_to_pdg(pids.MEMBER).get())
+
+    EXPECT_PDG_EQ(mu_minus);
+    EXPECT_PDG_EQ(neutron);
+    EXPECT_PDG_EQ(proton);
+    EXPECT_PDG_EQ(alpha);
+    EXPECT_PDG_EQ(he3);
+    EXPECT_PDG_EQ(muonic_hydrogen);
+    EXPECT_PDG_EQ(muonic_deuteron);
+    EXPECT_PDG_EQ(muonic_triton);
+    EXPECT_PDG_EQ(muonic_alpha);
+    EXPECT_PDG_EQ(muonic_he3);
+
+#undef EXPECT_PDG_EQ
+
+#define EXPECT_MASS_EQ(MEMBER) \
+    EXPECT_EQ(masses.MEMBER, particles_->get(pids.MEMBER).mass())
+
+    EXPECT_MASS_EQ(mu_minus);
+    EXPECT_MASS_EQ(neutron);
+    EXPECT_MASS_EQ(proton);
+    EXPECT_MASS_EQ(alpha);
+    EXPECT_MASS_EQ(he3);
+    EXPECT_MASS_EQ(muonic_hydrogen);
+    EXPECT_MASS_EQ(muonic_deuteron);
+    EXPECT_MASS_EQ(muonic_triton);
+    EXPECT_MASS_EQ(muonic_alpha);
+    EXPECT_MASS_EQ(muonic_he3);
+
+#undef EXPECT_MASS_EQ
+
+    EXPECT_EQ(21, data.muon_energy_cdf.grid.size());
+
+    // Check sotopic fractions
+    // Single material with 50/50 d and t fractions
+    EXPECT_SOFT_EQ(
+        0.5, data.isotopic_fractions[MuCfMatId{0}][MucfIsotope::deuterium]);
+    EXPECT_SOFT_EQ(
+        0.5, data.isotopic_fractions[MuCfMatId{0}][MucfIsotope::tritium]);
+
+    // Cycle times are in seconds
+    // DD (reactivity of F = 3/2 is almost negligible, with huge cycle times)
+    static double const expected_dd_1_over_2_cycle_time{1.8312922823566493e-06};
+    static double const expected_dd_3_over_2_cycle_time{1.1439517165483279};
+    // DT
+    static double const expected_dt_0_cycle_time{1.0182824459351898e-08};
+    static double const expected_dt_1_cycle_time{5.098478246172425e-09};
+    // TT
+    static double const expected_tt_1_over_2_cycle_time{1.4056833511329384e-06};
+
+    auto const& cycles = data.cycle_times;
+
+    // DD cycle times
+    EXPECT_SOFT_EQ(expected_dd_1_over_2_cycle_time,
+                   cycles[MuCfMatId{0}][Molecule::deuterium_deuterium][0]);
+    EXPECT_SOFT_EQ(expected_dd_3_over_2_cycle_time,
+                   cycles[MuCfMatId{0}][Molecule::deuterium_deuterium][1]);
+    // DT cycle times
+    EXPECT_SOFT_EQ(expected_dt_0_cycle_time,
+                   cycles[MuCfMatId{0}][Molecule::deuterium_tritium][0]);
+    EXPECT_SOFT_EQ(expected_dt_1_cycle_time,
+                   cycles[MuCfMatId{0}][Molecule::deuterium_tritium][1]);
+    // TT cycle times
+    EXPECT_SOFT_EQ(expected_tt_1_over_2_cycle_time,
+                   cycles[MuCfMatId{0}][Molecule::tritium_tritium][0]);
+    EXPECT_SOFT_EQ(0, cycles[MuCfMatId{0}][Molecule::tritium_tritium][1]);
+}
+
+//---------------------------------------------------------------------------//
+}  // namespace test
+}  // namespace celeritas
diff --git a/test/celeritas/mucf/DTMucfInteractor.test.cc b/test/celeritas/mucf/DTMucfInteractor.test.cc
index 1d2e47fd71..994b00cadd 100644
--- a/test/celeritas/mucf/DTMucfInteractor.test.cc
+++ b/test/celeritas/mucf/DTMucfInteractor.test.cc
@@ -33,33 +33,8 @@ class DTMucfInteractorTest : public MucfInteractorHostTestBase
   protected:
     void SetUp() override
     {
-        auto const& params = *this->particle_params();
-        this->set_material("hdt_fuel");
-
-        HostVal<DTMixMucfData> host_data;
-
-        // Set up particle IDs
-        host_data.particle_ids.mu_minus = params.find(pdg::mu_minus());
-        host_data.particle_ids.neutron = params.find(pdg::neutron());
-        host_data.particle_ids.alpha = params.find(pdg::alpha());
-        host_data.particle_ids.muonic_alpha = params.find(pdg::muonic_alpha());
-
-        // Set up particle masses
-        host_data.particle_masses.mu_minus
-            = params.get(host_data.particle_ids.mu_minus).mass();
-        host_data.particle_masses.neutron
-            = params.get(host_data.particle_ids.neutron).mass();
-        host_data.particle_masses.alpha
-            = params.get(host_data.particle_ids.alpha).mass();
-        host_data.particle_masses.muonic_alpha
-            = params.get(host_data.particle_ids.muonic_alpha).mass();
-
-        // Set up muon energy CDF
-        auto const inp_data = inp::MucfPhysics::from_default();
-        NonuniformGridBuilder build_grid_record{&host_data.reals};
-        host_data.muon_energy_cdf = build_grid_record(inp_data.muon_energy_cdf);
-
         // Construct collection
+        auto host_data = this->make_host_data();
         data_ = ParamsDataStore<DTMixMucfData>{std::move(host_data)};
 
         // At-rest muon primary
diff --git a/test/celeritas/mucf/MucfInteractorHostTestBase.cc b/test/celeritas/mucf/MucfInteractorHostTestBase.cc
index 1889f4b3be..43cde61501 100644
--- a/test/celeritas/mucf/MucfInteractorHostTestBase.cc
+++ b/test/celeritas/mucf/MucfInteractorHostTestBase.cc
@@ -7,6 +7,7 @@
 #include "MucfInteractorHostTestBase.hh"
 
 #include "celeritas/Units.hh"
+#include "celeritas/grid/NonuniformGridBuilder.hh"
 #include "celeritas/inp/MucfPhysics.hh"
 
 namespace celeritas
@@ -69,27 +70,27 @@ MucfInteractorHostBase::MucfInteractorHostBase()
          protium_mass,
          ElementaryCharge{1},
          stable_decay_constant},
-        {"neutron",
-         pdg::neutron(),
-         neutron_mass,
-         zero_quantity(),
-         stable_decay_constant},
-        {"deuterium",
-         pdg::deuteron(),
-         deuterium_mass,
-         ElementaryCharge{1},
-         stable_decay_constant},
         {"tritium",
          pdg::triton(),
          tritium_mass,
          ElementaryCharge{1},
          native_value_from(tritium_decay_constant)},
+        {"neutron",
+         pdg::neutron(),
+         neutron_mass,
+         zero_quantity(),
+         stable_decay_constant},
         {"alpha",
          pdg::alpha(),
          alpha_mass,
          ElementaryCharge{2},
          stable_decay_constant},
         {"he3", pdg::he3(), he3_mass, ElementaryCharge{2}, stable_decay_constant},
+        {"deuterium",
+         pdg::deuteron(),
+         deuterium_mass,
+         ElementaryCharge{1},
+         stable_decay_constant},
 
         // Muonic atoms
         {"muonic_hydrogen",
@@ -110,12 +111,12 @@ MucfInteractorHostBase::MucfInteractorHostBase()
         {"muonic_alpha",
          pdg::muonic_alpha(),
          alpha_mass + muon_mass,
-         ElementaryCharge{1},
+         ElementaryCharge{2},
          native_value_from(muon_decay_constant)},
         {"muonic_he3",
          pdg::muonic_he3(),
          he3_mass + muon_mass,
-         ElementaryCharge{1},
+         ElementaryCharge{2},
          native_value_from(muon_decay_constant)},
     };
     this->set_particle_params(std::move(par_inp));
@@ -179,6 +180,119 @@ MucfInteractorHostBase::MucfInteractorHostBase()
     this->set_material_params(std::move(mat_inp));
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Return a populated \c DTMixMucfData host data.
+ */
+HostVal<DTMixMucfData> MucfInteractorHostBase::make_host_data()
+{
+    using AtomicMassNumber = AtomicNumber;
+    using MaterialFractionsArray = EnumArray<MucfIsotope, real_type>;
+    using MoleculeCycles = Array<real_type, 2>;
+    using CycleTimesArray = EnumArray<MucfMuonicMolecule, MoleculeCycles>;
+
+    auto const& particles = *this->particle_params();
+    this->set_material("hdt_fuel");
+
+    HostVal<DTMixMucfData> host_data;
+
+    // Set up particle IDs
+    host_data.particle_ids.mu_minus = particles.find(pdg::mu_minus());
+
+    host_data.particle_ids.proton = particles.find(pdg::proton());
+    host_data.particle_ids.triton = particles.find(pdg::triton());
+    host_data.particle_ids.neutron = particles.find(pdg::neutron());
+    host_data.particle_ids.alpha = particles.find(pdg::alpha());
+    host_data.particle_ids.he3 = particles.find(pdg::he3());
+
+    host_data.particle_ids.muonic_hydrogen
+        = particles.find(pdg::muonic_hydrogen());
+    host_data.particle_ids.muonic_deuteron
+        = particles.find(pdg::muonic_deuteron());
+    host_data.particle_ids.muonic_triton = particles.find(pdg::muonic_triton());
+    host_data.particle_ids.muonic_alpha = particles.find(pdg::muonic_alpha());
+    host_data.particle_ids.muonic_he3 = particles.find(pdg::muonic_he3());
+
+    // Set up particle masses
+    host_data.particle_masses.mu_minus
+        = particles.get(host_data.particle_ids.mu_minus).mass();
+
+    host_data.particle_masses.proton
+        = particles.get(host_data.particle_ids.proton).mass();
+    host_data.particle_masses.triton
+        = particles.get(host_data.particle_ids.triton).mass();
+    host_data.particle_masses.neutron
+        = particles.get(host_data.particle_ids.neutron).mass();
+    host_data.particle_masses.alpha
+        = particles.get(host_data.particle_ids.alpha).mass();
+    host_data.particle_masses.he3
+        = particles.get(host_data.particle_ids.he3).mass();
+
+    host_data.particle_masses.muonic_hydrogen
+        = particles.get(host_data.particle_ids.muonic_hydrogen).mass();
+    host_data.particle_masses.muonic_deuteron
+        = particles.get(host_data.particle_ids.muonic_deuteron).mass();
+    host_data.particle_masses.muonic_triton
+        = particles.get(host_data.particle_ids.muonic_triton).mass();
+    host_data.particle_masses.muonic_alpha
+        = particles.get(host_data.particle_ids.muonic_alpha).mass();
+    host_data.particle_masses.muonic_he3
+        = particles.get(host_data.particle_ids.muonic_he3).mass();
+
+    // Set up muon energy CDF
+    auto const inp_data = inp::MucfPhysics::from_default();
+    NonuniformGridBuilder build_grid_record{&host_data.reals};
+    host_data.muon_energy_cdf = build_grid_record(inp_data.muon_energy_cdf);
+
+    auto const& material = *this->material_params();
+    auto const& el_view = material.get(ElementId{0});  // Only one element
+
+    MaterialFractionsArray iso_fractions_array;
+    for (auto const& frac : el_view.isotopes())
+    {
+        auto const& iso_view = material.get(frac.isotope);
+        if (iso_view.atomic_number() != AtomicNumber{1})
+        {
+            // Skip non-hydrogen (if added later to test)
+            continue;
+        }
+
+        // Set up isotopic fractions for D and T
+
+        if (iso_view.atomic_mass_number() == AtomicMassNumber{2})
+        {
+            iso_fractions_array[MucfIsotope::deuterium] = frac.fraction;
+        }
+        if (iso_view.atomic_mass_number() == AtomicMassNumber{3})
+        {
+            iso_fractions_array[MucfIsotope::tritium] = frac.fraction;
+        }
+    }
+
+    // Set up fractions
+    CollectionBuilder<MaterialFractionsArray, MemSpace::host, MuCfMatId>
+        host_iso_frac(&host_data.isotopic_fractions);
+    host_iso_frac.push_back(std::move(iso_fractions_array));
+
+    // Set up mucf material id to physics material id mapping
+    CollectionBuilder<PhysMatId, MemSpace::host, MuCfMatId> host_matid(
+        &host_data.mucfmatid_to_matid);
+    auto const& mat_view = material.get(PhysMatId{0});  // Only one material
+    host_matid.push_back(mat_view.material_id());
+
+    // Set up cycle times (numbers from DTMixMucfModel test)
+    CycleTimesArray ct_array;
+    ct_array[MucfMuonicMolecule::deuterium_deuterium] = {1.83e-6, 1.14};
+    ct_array[MucfMuonicMolecule::deuterium_tritium] = {1.018e-8, 5.098e-9};
+    ct_array[MucfMuonicMolecule::deuterium_tritium] = {1.40e-6, 0};
+
+    CollectionBuilder<CycleTimesArray, MemSpace::host, MuCfMatId> host_ct(
+        &host_data.cycle_times);
+    host_ct.push_back(std::move(ct_array));
+
+    return host_data;
+}
+
 //---------------------------------------------------------------------------//
 }  // namespace test
 }  // namespace celeritas
diff --git a/test/celeritas/mucf/MucfInteractorHostTestBase.hh b/test/celeritas/mucf/MucfInteractorHostTestBase.hh
index 5160828030..c89af305c3 100644
--- a/test/celeritas/mucf/MucfInteractorHostTestBase.hh
+++ b/test/celeritas/mucf/MucfInteractorHostTestBase.hh
@@ -6,6 +6,7 @@
 //---------------------------------------------------------------------------//
 #pragma once
 
+#include "celeritas/mucf/data/DTMixMucfData.hh"
 #include "celeritas/phys/InteractorHostTestBase.hh"
 
 namespace celeritas
@@ -27,6 +28,9 @@ class MucfInteractorHostBase : public InteractorHostBase
     MucfInteractorHostBase();
     ~MucfInteractorHostBase() = default;
     //!@}
+
+    // Construct MuCF data from test values for interactors
+    HostVal<DTMixMucfData> make_host_data();
 };
 
 class MucfInteractorHostTestBase : public MucfInteractorHostBase, public Test