diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp
index a4b9c6d5..a75b7e41 100644
--- a/src/CabanaPD_ForceModels.hpp
+++ b/src/CabanaPD_ForceModels.hpp
@@ -17,11 +17,11 @@
 namespace CabanaPD
 {
 
-template <typename... MemorySpace>
+template <typename MemorySpace>
 struct BaseForceModel;
 
 template <>
-struct BaseForceModel<>
+struct BaseForceModel<SingleSpecies>
 {
     using species_type = SingleSpecies;
     double delta;
@@ -38,14 +38,13 @@ struct BaseForceModel<>
     void thermalStretch( double&, const int, const int ) const {}
 };
 
-template <typename... MemorySpace>
+template <typename MemorySpace>
 struct BaseForceModel
 {
     using species_type = MultiSpecies;
 
     // Only allow one memory space.
-    using memory_space =
-        typename std::tuple_element<0, std::tuple<MemorySpace...>>::type;
+    using memory_space = MemorySpace;
     using view_type_1d = Kokkos::View<double*, memory_space>;
     view_type_1d delta;
     double max_delta;
@@ -86,8 +85,11 @@ struct BaseForceModel
     void thermalStretch( double&, const int, const int ) const {}
 };
 
+template <typename TemperatureType, typename SpeciesType>
+struct BaseTemperatureModel;
+
 template <typename TemperatureType>
-struct BaseTemperatureModel
+struct BaseTemperatureModel<TemperatureType, SingleSpecies>
 {
     using species_type = SingleSpecies;
     using memory_space = typename TemperatureType::memory_space;
@@ -115,6 +117,53 @@ struct BaseTemperatureModel
     }
 };
 
+template <typename TemperatureType, typename ParticleType>
+struct BaseTemperatureModel
+{
+    using species_type = MultiSpecies;
+    using memory_space = typename TemperatureType::memory_space;
+    using view_type_1d = Kokkos::View<double*, memory_space>;
+    view_type_1d alpha;
+    view_type_1d temp0;
+
+    // Particle fields
+    TemperatureType temperature;
+    ParticleType type;
+
+    template <typename ArrayType>
+    BaseTemperatureModel( const TemperatureType& _temp, const ArrayType _alpha,
+                          const ArrayType _temp0, const ParticleType& _type )
+        : alpha( view_type_1d( "delta", _alpha.size() ) )
+        , temp0( view_type_1d( "delta", _temp0.size() ) )
+        , temperature( _temp )
+        , type( _type )
+    {
+        auto init_func = KOKKOS_CLASS_LAMBDA( const int i )
+        {
+            alpha( i ) = _alpha[i];
+            temp0( i ) = _temp0[i];
+        };
+        using exec_space = typename memory_space::execution_space;
+        Kokkos::RangePolicy<exec_space> policy( 0, alpha.size() );
+        Kokkos::parallel_for( "CabanaPD::Model::Init", policy, init_func );
+    }
+
+    void update( const TemperatureType& _temp, const ParticleType& _type )
+    {
+        temperature = _temp;
+        type = _type;
+    }
+
+    // Update stretch with temperature effects.
+    KOKKOS_INLINE_FUNCTION
+    void thermalStretch( double& s, const int i, const int j ) const
+    {
+        double temp_avg =
+            0.5 * ( temperature( i ) + temperature( j ) ) - temp0( type( i ) );
+        s -= alpha( type( i ) ) * temp_avg;
+    }
+};
+
 template <typename ModelType, typename DamageType,
           typename ThermalType = TemperatureIndependent, typename... DataTypes>
 struct ForceModel;
diff --git a/src/force/CabanaPD_ForceModels_LPS.hpp b/src/force/CabanaPD_ForceModels_LPS.hpp
index 14bd9be8..b42eb30c 100644
--- a/src/force/CabanaPD_ForceModels_LPS.hpp
+++ b/src/force/CabanaPD_ForceModels_LPS.hpp
@@ -19,9 +19,9 @@ namespace CabanaPD
 {
 
 template <>
-struct ForceModel<LPS, Elastic> : public BaseForceModel<>
+struct ForceModel<LPS, Elastic> : public BaseForceModel<SingleSpecies>
 {
-    using base_type = BaseForceModel<>;
+    using base_type = BaseForceModel<SingleSpecies>;
     using species_type = SingleSpecies;
     using base_model = LPS;
     using fracture_type = Elastic;
diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp
index 96e7f5cc..c9bf438e 100644
--- a/src/force/CabanaPD_ForceModels_PMB.hpp
+++ b/src/force/CabanaPD_ForceModels_PMB.hpp
@@ -20,10 +20,10 @@ namespace CabanaPD
 
 template <>
 struct ForceModel<PMB, Elastic, TemperatureIndependent>
-    : public BaseForceModel<>
+    : public BaseForceModel<SingleSpecies>
 {
-    using base_type = BaseForceModel<>;
-    using species_type = typename base_type::species_type;
+    using base_type = BaseForceModel<SingleSpecies>;
+    using species_type = SingleSpecies;
     using base_model = PMB;
     using fracture_type = Elastic;
     using thermal_type = TemperatureIndependent;
@@ -326,12 +326,13 @@ struct ForceModel<LinearPMB, Fracture, TemperatureIndependent, ModelParams...>
 template <typename TemperatureType>
 struct ForceModel<PMB, Elastic, TemperatureDependent, TemperatureType>
     : public ForceModel<PMB, Elastic, TemperatureIndependent>,
-      BaseTemperatureModel<TemperatureType>
+      BaseTemperatureModel<TemperatureType, SingleSpecies>
 {
     using memory_space = typename TemperatureType::memory_space;
     using base_type = ForceModel<PMB, Elastic, TemperatureIndependent>;
     using species_type = typename base_type::species_type;
-    using base_temperature_type = BaseTemperatureModel<TemperatureType>;
+    using base_temperature_type =
+        BaseTemperatureModel<TemperatureType, SingleSpecies>;
     using base_model = PMB;
     using fracture_type = Elastic;
     using thermal_type = TemperatureDependent;
@@ -347,7 +348,6 @@ struct ForceModel<PMB, Elastic, TemperatureDependent, TemperatureType>
     // Explicitly use the temperature-dependent stretch.
     using base_temperature_type::thermalStretch;
 
-    // ForceModel(){};
     ForceModel( const double _delta, const double _K,
                 const TemperatureType _temp, const double _alpha,
                 const double _temp0 = 0.0 )
@@ -357,6 +357,44 @@ struct ForceModel<PMB, Elastic, TemperatureDependent, TemperatureType>
     }
 };
 
+template <typename TemperatureType, typename ParticleType>
+struct ForceModel<PMB, Elastic, TemperatureDependent, TemperatureType,
+                  ParticleType>
+    : public ForceModel<PMB, Elastic, TemperatureIndependent, ParticleType>,
+      BaseTemperatureModel<TemperatureType, ParticleType>
+{
+    using memory_space = typename TemperatureType::memory_space;
+    using base_type =
+        ForceModel<PMB, Elastic, TemperatureIndependent, ParticleType>;
+    using species_type = typename base_type::species_type;
+    using base_temperature_type =
+        BaseTemperatureModel<TemperatureType, ParticleType>;
+    using base_model = PMB;
+    using fracture_type = Elastic;
+    using thermal_type = TemperatureDependent;
+
+    using base_type::c;
+    using base_type::delta;
+    using base_type::K;
+    using base_type::type;
+
+    // Thermal parameters
+    using base_temperature_type::alpha;
+    using base_temperature_type::temp0;
+
+    // Explicitly use the temperature-dependent stretch.
+    using base_temperature_type::thermalStretch;
+
+    template <typename ArrayType>
+    ForceModel( const ArrayType _delta, const ArrayType _K,
+                const TemperatureType& _temp, const ArrayType _alpha,
+                const ArrayType _temp0, ParticleType& _type )
+        : base_type( _delta, _K, _type )
+        , base_temperature_type( _temp, _alpha, _temp0 )
+    {
+    }
+};
+
 template <typename ParticleType>
 auto createForceModel( PMB, Elastic, TemperatureDependent,
                        ParticleType particles, const double delta,
@@ -371,12 +409,13 @@ auto createForceModel( PMB, Elastic, TemperatureDependent,
 template <typename TemperatureType>
 struct ForceModel<PMB, Fracture, TemperatureDependent, TemperatureType>
     : public ForceModel<PMB, Fracture, TemperatureIndependent>,
-      BaseTemperatureModel<TemperatureType>
+      BaseTemperatureModel<TemperatureType, SingleSpecies>
 {
     using memory_space = typename TemperatureType::memory_space;
     using base_type = ForceModel<PMB, Fracture, TemperatureIndependent>;
     using species_type = typename base_type::species_type;
-    using base_temperature_type = BaseTemperatureModel<TemperatureType>;
+    using base_temperature_type =
+        BaseTemperatureModel<TemperatureType, SingleSpecies>;
     using base_model = typename base_type::base_model;
     using fracture_type = typename base_type::fracture_type;
     using thermal_type = TemperatureDependent;
@@ -397,7 +436,6 @@ struct ForceModel<PMB, Fracture, TemperatureDependent, TemperatureType>
     // Explicitly use the temperature-dependent stretch.
     using base_temperature_type::thermalStretch;
 
-    // ForceModel(){};
     ForceModel( const double _delta, const double _K, const double _G0,
                 const TemperatureType _temp, const double _alpha,
                 const double _temp0 = 0.0 )
@@ -417,6 +455,60 @@ struct ForceModel<PMB, Fracture, TemperatureDependent, TemperatureType>
     }
 };
 
+template <typename TemperatureType, typename ParticleType>
+struct ForceModel<PMB, Fracture, TemperatureDependent, TemperatureType,
+                  ParticleType>
+    : public ForceModel<PMB, Fracture, TemperatureIndependent, ParticleType>,
+      BaseTemperatureModel<TemperatureType, ParticleType>
+{
+    using memory_space = typename TemperatureType::memory_space;
+    using base_type =
+        ForceModel<PMB, Fracture, TemperatureIndependent, ParticleType>;
+    using species_type = typename base_type::species_type;
+    using base_temperature_type =
+        BaseTemperatureModel<TemperatureType, ParticleType>;
+    using base_model = typename base_type::base_model;
+    using fracture_type = typename base_type::fracture_type;
+    using thermal_type = TemperatureDependent;
+
+    using base_type::c;
+    using base_type::delta;
+    using base_type::K;
+    using base_type::type;
+
+    // Does not use the base bond_break_coeff.
+    using base_type::G0;
+    using base_type::s0;
+
+    // Thermal parameters
+    using base_temperature_type::alpha;
+    using base_temperature_type::temp0;
+    using base_temperature_type::temperature;
+
+    // Explicitly use the temperature-dependent stretch.
+    using base_temperature_type::thermalStretch;
+
+    template <typename ArrayType>
+    ForceModel( const ArrayType _delta, const ArrayType _K, const ArrayType _G0,
+                const TemperatureType _temp, const ArrayType _alpha,
+                const ArrayType _temp0, ParticleType& _type )
+        : base_type( _delta, _K, _G0, _type )
+        , base_temperature_type( _temp, _alpha, _temp0 )
+    {
+    }
+
+    KOKKOS_INLINE_FUNCTION
+    bool criticalStretch( const int i, const int j, const double r,
+                          const double xi ) const
+    {
+        double temp_avg =
+            0.5 * ( temperature( i ) + temperature( j ) ) - temp0( type( i ) );
+        double bond_break_coeff = ( 1.0 + s0 + alpha( type( i ) ) * temp_avg ) *
+                                  ( 1.0 + s0 + alpha( type( i ) ) * temp_avg );
+        return r * r >= bond_break_coeff * xi * xi;
+    }
+};
+
 template <typename ParticleType>
 auto createForceModel( PMB, Fracture, ParticleType particles,
                        const double delta, const double K, const double G0,