diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp
index deb3b688..8acd7f7e 100644
--- a/src/force/CabanaPD_ForceModels_PMB.hpp
+++ b/src/force/CabanaPD_ForceModels_PMB.hpp
@@ -53,6 +53,81 @@ struct ForceModel<PMB, NoFracture, Elastic, TemperatureIndependent>
         K = _K;
         c = 18.0 * K / ( pi * delta * delta * delta * delta );
     }
+
+    KOKKOS_INLINE_FUNCTION
+    bool forceCoeff( const double s, const double vol ) const
+    {
+        return c * s * vol;
+    }
+
+    KOKKOS_INLINE_FUNCTION
+    bool energyCoeff( const double s, const double xi, const double vol ) const
+    {
+        // 0.25 factor is due to 1/2 from outside the integral and 1/2 from
+        // the integrand (pairwise potential).
+        return 0.25 * c * s * s * xi * vol;
+    }
+};
+
+template <>
+struct ForceModel<PMB, NoFracture, Plastic, TemperatureIndependent>
+    : public ForceModel<PMB, NoFracture, Elastic, TemperatureIndependent>
+{
+    using base_type = ForceModel<PMB, NoFracture>;
+    using base_model = PMB;
+    using fracture_type = NoFracture;
+    using plasticity_type = Plastic;
+    using thermal_type = TemperatureIndependent;
+
+    using base_type::c;
+    using base_type::delta;
+    using base_type::K;
+    // FIXME: hardcoded
+    const double s_Y = 0.0014;
+
+    ForceModel(){};
+    ForceModel( const double delta, const double K )
+        : base_type( delta, K )
+    {
+        set_param( delta, K );
+    }
+
+    ForceModel( const ForceModel& model )
+        : base_type( model )
+    {
+        c = model.c;
+        K = model.K;
+    }
+
+    void set_param( const double _delta, const double _K )
+    {
+        delta = _delta;
+        K = _K;
+        c = 18.0 * K / ( pi * delta * delta * delta * delta );
+    }
+
+    KOKKOS_INLINE_FUNCTION
+    bool coeff( const double s, const double vol ) const
+    {
+        if ( s < s_Y )
+            return c * s * vol;
+        else
+            return c * s_Y * vol;
+    }
+
+    KOKKOS_INLINE_FUNCTION
+    bool energyCoeff( const double s, const double xi, const double vol ) const
+    {
+        // 0.25 factor is due to 1/2 from outside the integral and 1/2 from
+        // the integrand (pairwise potential).
+        double stretch_term = 0.0;
+        if ( s < s_Y )
+            stretch_term = s * s;
+        else
+            stretch_term = s_Y * ( 2 * s - s_Y );
+
+        return 0.25 * c * stretch_term * xi * vol;
+    }
 };
 
 template <>
@@ -105,8 +180,7 @@ struct ForceModel<PMB, Fracture, Elastic, TemperatureIndependent>
 
 template <>
 struct ForceModel<PMB, Fracture, Plastic, TemperatureIndependent>
-    : public ForceModel<PMB, Fracture, Elastic, TemperatureIndependent>,
-      ForceModel<PMB, NoFracture, Plastic, TemperatureIndependent>
+    : public ForceModel<PMB, Fracture, Elastic, TemperatureIndependent>
 {
     using base_type = ForceModel<PMB, Fracture>;
     using base_model = typename base_type::base_model;
@@ -120,7 +194,9 @@ struct ForceModel<PMB, Fracture, Plastic, TemperatureIndependent>
     using base_type::G0;
     using base_type::K;
     using base_type::s0;
-    double alpha = 0.25;
+
+    // FIXME: hardcoded
+    const double s_Y = 0.0014;
 
     ForceModel() {}
     ForceModel( const double delta, const double K, const double G0 )
@@ -128,16 +204,28 @@ struct ForceModel<PMB, Fracture, Plastic, TemperatureIndependent>
     {
     }
 
-    using base_type::minLocalStretch;
+    // FIXME: avoiding multiple inheritance..
+    KOKKOS_INLINE_FUNCTION
+    bool coeff( const double s, const double vol ) const
+    {
+        if ( s < s_Y )
+            return c * s * vol;
+        else
+            return c * s_Y * vol;
+    }
 
-    // Use the dynamic s0 from the minLocalStretch function, reused through the
-    // neighbor loop.
     KOKKOS_INLINE_FUNCTION
-    bool criticalStretch( const int, const int, const double r, const double xi,
-                          const double s0_i ) const
+    bool energyCoeff( const double s, const double xi, const double vol ) const
     {
-        double bond_break_coeff = ( 1.0 + s0_i ) * ( 1.0 + s0_i );
-        return r * r >= bond_break_coeff * xi * xi;
+        // 0.25 factor is due to 1/2 from outside the integral and 1/2 from
+        // the integrand (pairwise potential).
+        double stretch_term = 0.0;
+        if ( s < s_Y )
+            stretch_term = s * s;
+        else
+            stretch_term = s_Y * ( 2 * s - s_Y );
+
+        return 0.25 * c * stretch_term * xi * vol;
     }
 };
 
diff --git a/src/force/CabanaPD_Force_PMB.hpp b/src/force/CabanaPD_Force_PMB.hpp
index bd04bd90..02679321 100644
--- a/src/force/CabanaPD_Force_PMB.hpp
+++ b/src/force/CabanaPD_Force_PMB.hpp
@@ -117,7 +117,7 @@ class Force<ExecutionSpace, ForceModel<PMB, NoFracture, ModelParams...>>
 
             model.thermalStretch( s, i, j );
 
-            const double coeff = model.c * s * vol( j );
+            const double coeff = model.forceCoeff( s, vol( j ) );
             fx_i = coeff * rx / r;
             fy_i = coeff * ry / r;
             fz_i = coeff * rz / r;
@@ -158,7 +158,7 @@ class Force<ExecutionSpace, ForceModel<PMB, NoFracture, ModelParams...>>
 
             // 0.25 factor is due to 1/2 from outside the integral and 1/2 from
             // the integrand (pairwise potential).
-            double w = 0.25 * model.c * s * s * xi * vol( j );
+            double w = 0.25 * model.energyCoeff( s, xi, vol( j ) );
             W( i ) += w;
             Phi += w * vol( i );
         };
@@ -243,19 +243,7 @@ class Force<ExecutionSpace, ForceModel<PMB, Fracture, ModelParams...>>
                 // Else if statement is only for performance.
                 else if ( mu( i, n ) > 0 )
                 {
-                    const double s_Y = 0.0014;
-                    // const double coeff = model.c * s * vol( j );
-
-                    double coeff = 0;
-
-                    if ( s < s_Y )
-                    {
-                        coeff = model.c * s * vol( j );
-                    }
-                    else
-                    {
-                        coeff = model.c * s_Y * vol( j );
-                    }
+                    const double coeff = model.forceCoeff( s, vol( j ) );
 
                     double muij = mu( i, n );
                     fx_i = muij * coeff * rx / r;
@@ -306,23 +294,7 @@ class Force<ExecutionSpace, ForceModel<PMB, Fracture, ModelParams...>>
 
                 model.thermalStretch( s, i, j );
 
-                // 0.25 factor is due to 1/2 from outside the integral and 1/2
-                // from the integrand (pairwise potential).
-                // double w = mu( i, n ) * 0.25 * model.c * s * s * xi * vol( j
-                // );
-
-                const double s_Y = 0.0014;
-
-                if ( s < s_Y )
-                {
-                    w = mu( i, n ) * 0.25 * model.c * s * s * xi * vol( j );
-                }
-                else
-                {
-                    w = mu( i, n ) * 0.25 * model.c * s_Y * ( 2 * s - s_Y ) *
-                        xi * vol( j );
-                }
-
+                double w = mu( i, n ) * model.energyCoeff( s, xi, vol( j ) );
                 W( i ) += w;
 
                 phi_i += mu( i, n ) * vol( j );
@@ -391,7 +363,7 @@ class Force<ExecutionSpace, ForceModel<LinearPMB, NoFracture, ModelParams...>>
 
             model.thermalStretch( linear_s, i, j );
 
-            const double coeff = model.c * linear_s * vol( j );
+            const double coeff = model.forceCoeff( linear_s, vol( j ) );
             fx_i = coeff * xi_x / xi;
             fy_i = coeff * xi_y / xi;
             fz_i = coeff * xi_z / xi;