diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp
index 9ea0e040..6dab7f08 100644
--- a/src/CabanaPD_Force.hpp
+++ b/src/CabanaPD_Force.hpp
@@ -126,11 +126,11 @@ getLinearizedDistance( const PosType& x, const PosType& u, const int i,
 }
 
 // Forward declaration.
-template <class ExecutionSpace, class ForceType>
+template <class MemorySpace, class ForceType>
 class Force;
 
-template <class ExecutionSpace>
-class Force<ExecutionSpace, BaseForceModel>
+template <class MemorySpace>
+class Force<MemorySpace, BaseForceModel>
 {
   protected:
     bool _half_neigh;
@@ -139,19 +139,79 @@ class Force<ExecutionSpace, BaseForceModel>
     Timer _energy_timer;
 
   public:
-    Force( const bool half_neigh )
+    using neighbor_list_type =
+        Cabana::VerletList<MemorySpace, Cabana::FullNeighborTag,
+                           Cabana::VerletLayout2D, Cabana::TeamOpTag>;
+    neighbor_list_type _neigh_list;
+
+    // Primary constructor: use positions and construct neighbors.
+    template <class ParticleType>
+    Force( const bool half_neigh, const double delta,
+           const ParticleType& particles, const double tol = 1e-14 )
         : _half_neigh( half_neigh )
+        , _neigh_list( neighbor_list_type( particles.sliceReferencePosition(),
+                                           0, particles.n_local, delta + tol,
+                                           1.0, particles.ghost_mesh_lo,
+                                           particles.ghost_mesh_hi ) )
     {
     }
 
-    template <class ParticleType, class NeighListType, class ParallelType>
-    void computeWeightedVolume( ParticleType&, const NeighListType&,
-                                const ParallelType ) const
+    // Primary constructor: use positions and construct neighbors.
+    template <class NeighborListType>
+    Force( const bool half_neigh, const NeighborListType& neighbors )
+        : _half_neigh( half_neigh )
+        , _neigh_list( neighbors )
+    {
+    }
+
+    unsigned getMaxLocalNeighbors()
     {
+        auto neigh = _neigh_list;
+        unsigned local_max_neighbors;
+        auto neigh_max = KOKKOS_LAMBDA( const int, unsigned& max_n )
+        {
+            max_n =
+                Cabana::NeighborList<neighbor_list_type>::maxNeighbor( neigh );
+        };
+        using exec_space = typename MemorySpace::execution_space;
+        Kokkos::RangePolicy<exec_space> policy( 0, 1 );
+        Kokkos::parallel_reduce( policy, neigh_max, local_max_neighbors );
+        Kokkos::fence();
+        return local_max_neighbors;
     }
-    template <class ParticleType, class NeighListType, class ParallelType>
-    void computeDilatation( ParticleType&, const NeighListType&,
-                            const ParallelType ) const
+
+    void getNeighborStatistics( unsigned& max_neighbors,
+                                unsigned long long& total_neighbors )
+    {
+        auto neigh = _neigh_list;
+        unsigned local_max_neighbors;
+        unsigned long long local_total_neighbors;
+        auto neigh_stats = KOKKOS_LAMBDA( const int, unsigned& max_n,
+                                          unsigned long long& total_n )
+        {
+            max_n =
+                Cabana::NeighborList<neighbor_list_type>::maxNeighbor( neigh );
+            total_n = Cabana::NeighborList<neighbor_list_type>::totalNeighbor(
+                neigh );
+        };
+        using exec_space = typename MemorySpace::execution_space;
+        Kokkos::RangePolicy<exec_space> policy( 0, 1 );
+        Kokkos::parallel_reduce( policy, neigh_stats, local_max_neighbors,
+                                 local_total_neighbors );
+        Kokkos::fence();
+        MPI_Reduce( &local_max_neighbors, &max_neighbors, 1, MPI_UNSIGNED,
+                    MPI_MAX, 0, MPI_COMM_WORLD );
+        MPI_Reduce( &local_total_neighbors, &total_neighbors, 1,
+                    MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD );
+    }
+
+    // Default to no-op for pair models.
+    template <class ParticleType, class ParallelType>
+    void computeWeightedVolume( ParticleType&, const ParallelType ) const
+    {
+    }
+    template <class ParticleType, class ParallelType>
+    void computeDilatation( ParticleType&, const ParallelType ) const
     {
     }
 
@@ -162,10 +222,8 @@ class Force<ExecutionSpace, BaseForceModel>
 /******************************************************************************
   Force free functions.
 ******************************************************************************/
-template <class ForceType, class ParticleType, class NeighListType,
-          class ParallelType>
+template <class ForceType, class ParticleType, class ParallelType>
 void computeForce( ForceType& force, ParticleType& particles,
-                   const NeighListType& neigh_list,
                    const ParallelType& neigh_op_tag )
 {
     auto n_local = particles.n_local;
@@ -179,23 +237,19 @@ void computeForce( ForceType& force, ParticleType& particles,
 
     // if ( half_neigh )
     // Forces must be atomic for half list
-    // computeForce_half( f_a, x, u, neigh_list, n_local,
+    // computeForce_half( f_a, x, u, n_local,
     //                    neigh_op_tag );
 
     // Forces only atomic if using team threading.
     if ( std::is_same<decltype( neigh_op_tag ), Cabana::TeamOpTag>::value )
-        force.computeForceFull( f_a, x, u, particles, neigh_list, n_local,
-                                neigh_op_tag );
+        force.computeForceFull( f_a, x, u, particles, n_local, neigh_op_tag );
     else
-        force.computeForceFull( f, x, u, particles, neigh_list, n_local,
-                                neigh_op_tag );
+        force.computeForceFull( f, x, u, particles, n_local, neigh_op_tag );
     Kokkos::fence();
 }
 
-template <class ForceType, class ParticleType, class NeighListType,
-          class ParallelType>
+template <class ForceType, class ParticleType, class ParallelType>
 double computeEnergy( ForceType& force, ParticleType& particles,
-                      const NeighListType& neigh_list,
                       const ParallelType& neigh_op_tag )
 {
     auto n_local = particles.n_local;
@@ -210,21 +264,20 @@ double computeEnergy( ForceType& force, ParticleType& particles,
 
     double energy;
     // if ( _half_neigh )
-    //    energy = computeEnergy_half( force, x, u, neigh_list,
+    //    energy = computeEnergy_half( force, x, u,
     //                                  n_local, neigh_op_tag );
     // else
-    energy = force.computeEnergyFull( W, x, u, particles, neigh_list, n_local,
-                                      neigh_op_tag );
+    energy =
+        force.computeEnergyFull( W, x, u, particles, n_local, neigh_op_tag );
     Kokkos::fence();
 
     return energy;
 }
 
 // Forces with bond breaking.
-template <class ForceType, class ParticleType, class NeighListType,
-          class NeighborView, class ParallelType>
-void computeForce( ForceType& force, ParticleType& particles,
-                   const NeighListType& neigh_list, NeighborView& mu,
+template <class ForceType, class ParticleType, class NeighborView,
+          class ParallelType>
+void computeForce( ForceType& force, ParticleType& particles, NeighborView& mu,
                    const ParallelType& neigh_op_tag )
 {
     auto n_local = particles.n_local;
@@ -238,25 +291,23 @@ void computeForce( ForceType& force, ParticleType& particles,
 
     // if ( half_neigh )
     // Forces must be atomic for half list
-    // computeForce_half( f_a, x, u, neigh_list, n_local,
+    // computeForce_half( f_a, x, u, n_local,
     //                    neigh_op_tag );
 
     // Forces only atomic if using team threading.
     if ( std::is_same<decltype( neigh_op_tag ), Cabana::TeamOpTag>::value )
-        force.computeForceFull( f_a, x, u, particles, neigh_list, mu, n_local,
+        force.computeForceFull( f_a, x, u, particles, mu, n_local,
                                 neigh_op_tag );
     else
-        force.computeForceFull( f, x, u, particles, neigh_list, mu, n_local,
-                                neigh_op_tag );
+        force.computeForceFull( f, x, u, particles, mu, n_local, neigh_op_tag );
     Kokkos::fence();
 }
 
 // Energy and damage.
-template <class ForceType, class ParticleType, class NeighListType,
-          class NeighborView, class ParallelType>
+template <class ForceType, class ParticleType, class NeighborView,
+          class ParallelType>
 double computeEnergy( ForceType& force, ParticleType& particles,
-                      const NeighListType& neigh_list, NeighborView& mu,
-                      const ParallelType& neigh_op_tag )
+                      NeighborView& mu, const ParallelType& neigh_op_tag )
 {
     auto n_local = particles.n_local;
     auto x = particles.sliceReferencePosition();
@@ -270,11 +321,11 @@ double computeEnergy( ForceType& force, ParticleType& particles,
 
     double energy;
     // if ( _half_neigh )
-    //    energy = computeEnergy_half( force, x, u, neigh_list,
+    //    energy = computeEnergy_half( force, x, u,
     //                                  n_local, neigh_op_tag );
     // else
-    energy = force.computeEnergyFull( W, x, u, phi, particles, neigh_list, mu,
-                                      n_local, neigh_op_tag );
+    energy = force.computeEnergyFull( W, x, u, phi, particles, mu, n_local,
+                                      neigh_op_tag );
     Kokkos::fence();
 
     return energy;
diff --git a/src/CabanaPD_HeatTransfer.hpp b/src/CabanaPD_HeatTransfer.hpp
index 75a7531c..b3818fd0 100644
--- a/src/CabanaPD_HeatTransfer.hpp
+++ b/src/CabanaPD_HeatTransfer.hpp
@@ -19,34 +19,40 @@ namespace CabanaPD
 
 // Peridynamic heat transfer with forward-Euler time integration.
 // Inherits only because this is a similar neighbor-based kernel.
-template <class ExecutionSpace, class ModelType>
-class HeatTransfer : public Force<ExecutionSpace, BaseForceModel>
+template <class MemorySpace, class ModelType>
+class HeatTransfer : public Force<MemorySpace, BaseForceModel>
 {
   protected:
-    using base_type = Force<ExecutionSpace, BaseForceModel>;
+    using base_type = Force<MemorySpace, BaseForceModel>;
     using base_type::_half_neigh;
     using base_type::_timer;
-    using model_type = ModelType;
-    static_assert(
-        std::is_same_v<typename model_type::fracture_type, Elastic> );
 
     Timer _euler_timer = base_type::_energy_timer;
     ModelType _model;
+    // Using the default exec_space.
+    using exec_space = typename MemorySpace::execution_space;
 
   public:
-    HeatTransfer( const bool half_neigh, const ModelType model )
-        : base_type( half_neigh )
+    using base_type::_neigh_list;
+    using model_type = ModelType;
+    static_assert(
+        std::is_same_v<typename model_type::fracture_type, Elastic> );
+
+    // Running with mechanics as well; no reason to rebuild neighbors.
+    template <class NeighborType>
+    HeatTransfer( const bool half_neigh, const NeighborType& neighbors,
+                  const ModelType model )
+        : base_type( half_neigh, neighbors )
         , _model( model )
     {
     }
 
     template <class TemperatureType, class PosType, class ParticleType,
-              class NeighListType, class ParallelType>
+              class ParallelType>
     void
     computeHeatTransferFull( TemperatureType& conduction, const PosType& x,
                              const PosType& u, const ParticleType& particles,
-                             const NeighListType& neigh_list, const int n_local,
-                             ParallelType& neigh_op_tag )
+                             const int n_local, ParallelType& neigh_op_tag )
     {
         _timer.start();
 
@@ -64,9 +70,9 @@ class HeatTransfer : public Force<ExecutionSpace, BaseForceModel>
                 coeff * ( temp( j ) - temp( i ) ) / xi / xi * vol( j );
         };
 
-        Kokkos::RangePolicy<ExecutionSpace> policy( 0, n_local );
+        Kokkos::RangePolicy<exec_space> policy( 0, n_local );
         Cabana::neighbor_parallel_for(
-            policy, temp_func, neigh_list, Cabana::FirstNeighborsTag(),
+            policy, temp_func, _neigh_list, Cabana::FirstNeighborsTag(),
             neigh_op_tag, "CabanaPD::HeatTransfer::computeFull" );
 
         _timer.stop();
@@ -85,7 +91,7 @@ class HeatTransfer : public Force<ExecutionSpace, BaseForceModel>
         {
             temp( i ) += dt / rho( i ) / model.cp * conduction( i );
         };
-        Kokkos::RangePolicy<ExecutionSpace> policy( 0, n_local );
+        Kokkos::RangePolicy<exec_space> policy( 0, n_local );
         Kokkos::parallel_for( "CabanaPD::HeatTransfer::forwardEuler", policy,
                               euler_func );
         _euler_timer.stop();
@@ -93,11 +99,9 @@ class HeatTransfer : public Force<ExecutionSpace, BaseForceModel>
 };
 
 // Heat transfer free function.
-template <class HeatTransferType, class ParticleType, class NeighListType,
-          class ParallelType>
+template <class HeatTransferType, class ParticleType, class ParallelType>
 void computeHeatTransfer( HeatTransferType& heat_transfer,
                           ParticleType& particles,
-                          const NeighListType& neigh_list,
                           const ParallelType& neigh_op_tag, const double dt )
 {
     auto n_local = particles.n_local;
@@ -111,11 +115,11 @@ void computeHeatTransfer( HeatTransferType& heat_transfer,
 
     // Temperature only needs to be atomic if using team threading.
     if ( std::is_same<decltype( neigh_op_tag ), Cabana::TeamOpTag>::value )
-        heat_transfer.computeHeatTransferFull(
-            conduction_a, x, u, particles, neigh_list, n_local, neigh_op_tag );
+        heat_transfer.computeHeatTransferFull( conduction_a, x, u, particles,
+                                               n_local, neigh_op_tag );
     else
-        heat_transfer.computeHeatTransferFull(
-            conduction, x, u, particles, neigh_list, n_local, neigh_op_tag );
+        heat_transfer.computeHeatTransferFull( conduction, x, u, particles,
+                                               n_local, neigh_op_tag );
     Kokkos::fence();
 
     heat_transfer.forwardEuler( particles, dt );
diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp
index 02843e57..52cc36ed 100644
--- a/src/CabanaPD_Particles.hpp
+++ b/src/CabanaPD_Particles.hpp
@@ -134,8 +134,9 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
     std::array<double, dim> local_mesh_ext;
     std::array<double, dim> local_mesh_lo;
     std::array<double, dim> local_mesh_hi;
-    std::array<double, dim> ghost_mesh_lo;
-    std::array<double, dim> ghost_mesh_hi;
+    // FIXME: this is for neighborlist construction.
+    double ghost_mesh_lo[dim];
+    double ghost_mesh_hi[dim];
     std::shared_ptr<
         Cabana::Grid::LocalGrid<Cabana::Grid::UniformMesh<double, dim>>>
         local_grid;
@@ -399,9 +400,9 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
     auto getForce() { return _plist_f; }
     auto getReferencePosition() { return _plist_x; }
 
-    void updateCurrentPosition()
+    void updateCurrentPosition() const
     {
-        _timer.start();
+        //_timer.start();
         // Not using slice function because this is called inside.
         auto y = Cabana::slice<0>( _aosoa_y, "current_positions" );
         auto x = sliceReferencePosition();
@@ -414,7 +415,7 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
         };
         Kokkos::parallel_for( "CabanaPD::CalculateCurrentPositions", policy,
                               sum_x_u );
-        _timer.stop();
+        //_timer.stop();
     }
 
     void resize( int new_local, int new_ghost )
diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp
index 44c7ab91..d7a99321 100644
--- a/src/CabanaPD_Solver.hpp
+++ b/src/CabanaPD_Solver.hpp
@@ -103,17 +103,14 @@ class SolverElastic
     using particle_type = ParticleType;
     using integrator_type = Integrator<exec_space>;
     using force_model_type = ForceModel;
-    using force_type = Force<exec_space, force_model_type>;
+    using force_type = Force<memory_space, force_model_type>;
     using comm_type = Comm<particle_type, typename force_model_type::base_model,
                            typename particle_type::thermal_type>;
-    using neighbor_type =
-        Cabana::VerletList<memory_space, Cabana::FullNeighborTag,
-                           Cabana::VerletLayout2D, Cabana::TeamOpTag>;
     using neigh_iter_tag = Cabana::SerialOpTag;
     using input_type = InputType;
 
     // Optional module types.
-    using heat_transfer_type = HeatTransfer<exec_space, force_model_type>;
+    using heat_transfer_type = HeatTransfer<memory_space, force_model_type>;
 
     SolverElastic( input_type _inputs,
                    std::shared_ptr<particle_type> _particles,
@@ -138,35 +135,24 @@ class SolverElastic
                            typename force_model_type::thermal_type>::value )
             force_model.update( particles->sliceTemperature() );
 
+        _neighbor_timer.start();
+        force = std::make_shared<force_type>( inputs["half_neigh"], *particles,
+                                              force_model );
+        _neighbor_timer.stop();
+
+        _init_timer.start();
+        unsigned max_neighbors;
+        unsigned long long total_neighbors;
+        force->getNeighborStatistics( max_neighbors, total_neighbors );
+
         // Create heat transfer if needed.
         if constexpr ( is_heat_transfer<
                            typename force_model_type::thermal_type>::value )
         {
             thermal_subcycle_steps = inputs["thermal_subcycle_steps"];
             heat_transfer = std::make_shared<heat_transfer_type>(
-                inputs["half_neigh"], force_model );
+                inputs["half_neigh"], force->_neigh_list, force_model );
         }
-        force =
-            std::make_shared<force_type>( inputs["half_neigh"], force_model );
-
-        // Create the neighbor list.
-        _neighbor_timer.start();
-        double mesh_min[3] = { particles->ghost_mesh_lo[0],
-                               particles->ghost_mesh_lo[1],
-                               particles->ghost_mesh_lo[2] };
-        double mesh_max[3] = { particles->ghost_mesh_hi[0],
-                               particles->ghost_mesh_hi[1],
-                               particles->ghost_mesh_hi[2] };
-        auto x = particles->sliceReferencePosition();
-        neighbors = std::make_shared<neighbor_type>( x, 0, particles->n_local,
-                                                     force_model.delta, 1.0,
-                                                     mesh_min, mesh_max );
-        _neighbor_timer.stop();
-
-        _init_timer.start();
-        unsigned max_neighbors;
-        unsigned long long total_neighbors;
-        getNeighborStatistics( max_neighbors, total_neighbors );
 
         print = print_rank();
         if ( print )
@@ -195,40 +181,16 @@ class SolverElastic
         _init_timer.stop();
     }
 
-    void getNeighborStatistics( unsigned& max_neighbors,
-                                unsigned long long& total_neighbors )
-    {
-        auto neigh = *neighbors;
-        unsigned local_max_neighbors;
-        unsigned long long local_total_neighbors;
-        auto neigh_stats = KOKKOS_LAMBDA( const int, unsigned& max_n,
-                                          unsigned long long& total_n )
-        {
-            max_n = Cabana::NeighborList<neighbor_type>::maxNeighbor( neigh );
-            total_n =
-                Cabana::NeighborList<neighbor_type>::totalNeighbor( neigh );
-        };
-        Kokkos::RangePolicy<exec_space> policy( 0, 1 );
-        Kokkos::parallel_reduce( policy, neigh_stats, local_max_neighbors,
-                                 local_total_neighbors );
-        Kokkos::fence();
-        MPI_Reduce( &local_max_neighbors, &max_neighbors, 1, MPI_UNSIGNED,
-                    MPI_MAX, 0, MPI_COMM_WORLD );
-        MPI_Reduce( &local_total_neighbors, &total_neighbors, 1,
-                    MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD );
-    }
-
     void init( const bool initial_output = true )
     {
         // Compute and communicate weighted volume for LPS (does nothing for
         // PMB). Only computed once.
-        force->computeWeightedVolume( *particles, *neighbors,
-                                      neigh_iter_tag{} );
+        force->computeWeightedVolume( *particles, neigh_iter_tag{} );
         comm->gatherWeightedVolume();
 
         // Compute initial internal forces and energy.
         updateForce();
-        computeEnergy( *force, *particles, *neighbors, neigh_iter_tag() );
+        computeEnergy( *force, *particles, neigh_iter_tag() );
 
         if ( initial_output )
             particles->output( 0, 0.0, output_reference );
@@ -278,7 +240,7 @@ class SolverElastic
                                typename force_model_type::thermal_type>::value )
             {
                 if ( step % thermal_subcycle_steps == 0 )
-                    computeHeatTransfer( *heat_transfer, *particles, *neighbors,
+                    computeHeatTransfer( *heat_transfer, *particles,
                                          neigh_iter_tag{}, dt );
             }
 
@@ -344,11 +306,11 @@ class SolverElastic
     void updateForce()
     {
         // Compute and communicate dilatation for LPS (does nothing for PMB).
-        force->computeDilatation( *particles, *neighbors, neigh_iter_tag{} );
+        force->computeDilatation( *particles, neigh_iter_tag{} );
         comm->gatherDilatation();
 
         // Compute internal forces.
-        computeForce( *force, *particles, *neighbors, neigh_iter_tag{} );
+        computeForce( *force, *particles, neigh_iter_tag{} );
     }
 
     void output( const int step )
@@ -356,8 +318,7 @@ class SolverElastic
         // Print output.
         if ( step % output_frequency == 0 )
         {
-            auto W = computeEnergy( *force, *particles, *neighbors,
-                                    neigh_iter_tag() );
+            auto W = computeEnergy( *force, *particles, neigh_iter_tag() );
 
             particles->output( step / output_frequency, step * dt,
                                output_reference );
@@ -459,7 +420,6 @@ class SolverElastic
     std::shared_ptr<comm_type> comm;
     std::shared_ptr<integrator_type> integrator;
     std::shared_ptr<force_type> force;
-    std::shared_ptr<neighbor_type> neighbors;
     // Optional modules.
     std::shared_ptr<heat_transfer_type> heat_transfer;
 
@@ -490,7 +450,6 @@ class SolverFracture
     using particle_type = typename base_type::particle_type;
     using integrator_type = typename base_type::integrator_type;
     using comm_type = typename base_type::comm_type;
-    using neighbor_type = typename base_type::neighbor_type;
     using force_model_type = ForceModel;
     using force_type = typename base_type::force_type;
     using neigh_iter_tag = Cabana::SerialOpTag;
@@ -505,7 +464,7 @@ class SolverFracture
         init_mu();
 
         // Create prenotch.
-        prenotch.create( exec_space{}, mu, *particles, *neighbors );
+        prenotch.create( exec_space{}, mu, *particles, force->_neigh_list );
         _init_time += prenotch.time();
     }
 
@@ -521,8 +480,7 @@ class SolverFracture
     {
         _init_timer.start();
         // Create View to track broken bonds.
-        int max_neighbors =
-            Cabana::NeighborList<neighbor_type>::maxNeighbor( *neighbors );
+        auto max_neighbors = force->getMaxLocalNeighbors();
         mu = NeighborView(
             Kokkos::ViewAllocateWithoutInitializing( "broken_bonds" ),
             particles->n_local, max_neighbors );
@@ -534,7 +492,7 @@ class SolverFracture
     {
         // Compute initial internal forces and energy.
         updateForce();
-        computeEnergy( *force, *particles, *neighbors, mu, neigh_iter_tag() );
+        computeEnergy( *force, *particles, mu, neigh_iter_tag() );
 
         if ( initial_output )
             particles->output( 0, 0.0, output_reference );
@@ -643,15 +601,15 @@ class SolverFracture
     {
         // Compute and communicate weighted volume for LPS (does nothing for
         // PMB).
-        force->computeWeightedVolume( *particles, *neighbors, mu );
+        force->computeWeightedVolume( *particles, mu );
         comm->gatherWeightedVolume();
 
         // Compute and communicate dilatation for LPS (does nothing for PMB).
-        force->computeDilatation( *particles, *neighbors, mu );
+        force->computeDilatation( *particles, mu );
         comm->gatherDilatation();
 
         // Compute internal forces.
-        computeForce( *force, *particles, *neighbors, mu, neigh_iter_tag{} );
+        computeForce( *force, *particles, mu, neigh_iter_tag{} );
     }
 
     void output( const int step )
@@ -659,8 +617,7 @@ class SolverFracture
         // Print output.
         if ( step % output_frequency == 0 )
         {
-            auto W = computeEnergy( *force, *particles, *neighbors, mu,
-                                    neigh_iter_tag() );
+            auto W = computeEnergy( *force, *particles, mu, neigh_iter_tag() );
 
             particles->output( step / output_frequency, step * dt,
                                output_reference );
@@ -683,7 +640,6 @@ class SolverFracture
     using base_type::force;
     using base_type::inputs;
     using base_type::integrator;
-    using base_type::neighbors;
     using base_type::particles;
 
     using NeighborView = typename Kokkos::View<int**, memory_space>;
diff --git a/src/force/CabanaPD_Force_LPS.hpp b/src/force/CabanaPD_Force_LPS.hpp
index 78f25f81..ae28bc76 100644
--- a/src/force/CabanaPD_Force_LPS.hpp
+++ b/src/force/CabanaPD_Force_LPS.hpp
@@ -69,28 +69,34 @@
 
 namespace CabanaPD
 {
-template <class ExecutionSpace>
-class Force<ExecutionSpace, ForceModel<LPS, Elastic>>
+template <class MemorySpace>
+class Force<MemorySpace, ForceModel<LPS, Elastic>>
+    : public Force<MemorySpace, BaseForceModel>
 {
   protected:
-    bool _half_neigh;
+    using base_type = Force<MemorySpace, BaseForceModel>;
+    using base_type::_half_neigh;
     ForceModel<LPS, Elastic> _model;
 
-    Timer _timer;
-    Timer _energy_timer;
+    using base_type::_energy_timer;
+    using base_type::_timer;
 
   public:
-    using exec_space = ExecutionSpace;
-
-    Force( const bool half_neigh, const ForceModel<LPS, Elastic> model )
-        : _half_neigh( half_neigh )
+    // Using the default exec_space.
+    using exec_space = typename MemorySpace::execution_space;
+    using neighbor_list_type = typename base_type::neighbor_list_type;
+    using base_type::_neigh_list;
+
+    template <class ParticleType>
+    Force( const bool half_neigh, ParticleType& particles,
+           const ForceModel<LPS, Elastic> model )
+        : base_type( half_neigh, model.delta, particles )
         , _model( model )
     {
     }
 
-    template <class ParticleType, class NeighListType, class ParallelType>
+    template <class ParticleType, class ParallelType>
     void computeWeightedVolume( ParticleType& particles,
-                                const NeighListType& neigh_list,
                                 const ParallelType neigh_op_tag )
     {
         _timer.start();
@@ -114,15 +120,14 @@ class Force<ExecutionSpace, ForceModel<LPS, Elastic>>
 
         Kokkos::RangePolicy<exec_space> policy( 0, n_local );
         Cabana::neighbor_parallel_for(
-            policy, weighted_volume, neigh_list, Cabana::FirstNeighborsTag(),
+            policy, weighted_volume, _neigh_list, Cabana::FirstNeighborsTag(),
             neigh_op_tag, "CabanaPD::ForceLPS::computeWeightedVolume" );
 
         _timer.stop();
     }
 
-    template <class ParticleType, class NeighListType, class ParallelType>
+    template <class ParticleType, class ParallelType>
     void computeDilatation( ParticleType& particles,
-                            const NeighListType& neigh_list,
                             const ParallelType neigh_op_tag )
     {
         _timer.start();
@@ -148,17 +153,16 @@ class Force<ExecutionSpace, ForceModel<LPS, Elastic>>
 
         Kokkos::RangePolicy<exec_space> policy( 0, n_local );
         Cabana::neighbor_parallel_for(
-            policy, dilatation, neigh_list, Cabana::FirstNeighborsTag(),
+            policy, dilatation, _neigh_list, Cabana::FirstNeighborsTag(),
             neigh_op_tag, "CabanaPD::ForceLPS::computeDilatation" );
 
         _timer.stop();
     }
 
     template <class ForceType, class PosType, class ParticleType,
-              class NeighListType, class ParallelType>
+              class ParallelType>
     void computeForceFull( ForceType& f, const PosType& x, const PosType& u,
-                           const ParticleType& particles,
-                           const NeighListType& neigh_list, const int n_local,
+                           const ParticleType& particles, const int n_local,
                            ParallelType& neigh_op_tag )
     {
         _timer.start();
@@ -194,24 +198,24 @@ class Force<ExecutionSpace, ForceModel<LPS, Elastic>>
 
         Kokkos::RangePolicy<exec_space> policy( 0, n_local );
         Cabana::neighbor_parallel_for(
-            policy, force_full, neigh_list, Cabana::FirstNeighborsTag(),
+            policy, force_full, _neigh_list, Cabana::FirstNeighborsTag(),
             neigh_op_tag, "CabanaPD::ForceLPS::computeFull" );
 
         _timer.stop();
     }
 
     template <class PosType, class WType, class ParticleType,
-              class NeighListType, class ParallelType>
+              class ParallelType>
     double computeEnergyFull( WType& W, const PosType& x, const PosType& u,
-                              const ParticleType& particles,
-                              const NeighListType& neigh_list,
-                              const int n_local, ParallelType& neigh_op_tag )
+                              const ParticleType& particles, const int n_local,
+                              ParallelType& neigh_op_tag )
     {
         _energy_timer.start();
 
         auto theta_coeff = _model.theta_coeff;
         auto s_coeff = _model.s_coeff;
         auto model = _model;
+        auto neigh_list = _neigh_list;
 
         const auto vol = particles.sliceVolume();
         const auto theta = particles.sliceDilatation();
@@ -224,8 +228,8 @@ class Force<ExecutionSpace, ForceModel<LPS, Elastic>>
             getDistance( x, u, i, j, xi, r, s );
 
             double num_neighbors = static_cast<double>(
-                Cabana::NeighborList<NeighListType>::numNeighbor( neigh_list,
-                                                                  i ) );
+                Cabana::NeighborList<neighbor_list_type>::numNeighbor(
+                    neigh_list, i ) );
 
             double w = ( 1.0 / num_neighbors ) * 0.5 * theta_coeff / 3.0 *
                            ( theta( i ) * theta( i ) ) +
@@ -240,7 +244,7 @@ class Force<ExecutionSpace, ForceModel<LPS, Elastic>>
 
         Kokkos::RangePolicy<exec_space> policy( 0, n_local );
         Cabana::neighbor_parallel_reduce(
-            policy, energy_full, neigh_list, Cabana::FirstNeighborsTag(),
+            policy, energy_full, _neigh_list, Cabana::FirstNeighborsTag(),
             neigh_op_tag, strain_energy,
             "CabanaPD::ForceLPS::computeEnergyFull" );
 
@@ -252,12 +256,12 @@ class Force<ExecutionSpace, ForceModel<LPS, Elastic>>
     auto timeEnergy() { return _energy_timer.time(); };
 };
 
-template <class ExecutionSpace>
-class Force<ExecutionSpace, ForceModel<LPS, Fracture>>
-    : public Force<ExecutionSpace, ForceModel<LPS, Elastic>>
+template <class MemorySpace>
+class Force<MemorySpace, ForceModel<LPS, Fracture>>
+    : public Force<MemorySpace, ForceModel<LPS, Elastic>>
 {
   protected:
-    using base_type = Force<ExecutionSpace, ForceModel<LPS, Elastic>>;
+    using base_type = Force<MemorySpace, ForceModel<LPS, Elastic>>;
     using base_type::_half_neigh;
     ForceModel<LPS, Fracture> _model;
 
@@ -265,18 +269,21 @@ class Force<ExecutionSpace, ForceModel<LPS, Fracture>>
     using base_type::_timer;
 
   public:
-    using exec_space = ExecutionSpace;
-
-    Force( const bool half_neigh, const ForceModel<LPS, Fracture> model )
-        : base_type( half_neigh, model )
+    // Using the default exec_space.
+    using exec_space = typename MemorySpace::execution_space;
+    using neighbor_list_type = typename base_type::neighbor_list_type;
+    using base_type::_neigh_list;
+
+    template <class ParticleType>
+    Force( const bool half_neigh, const ParticleType& particles,
+           const ForceModel<LPS, Fracture> model )
+        : base_type( half_neigh, particles, model )
         , _model( model )
     {
     }
 
-    template <class ParticleType, class NeighListType, class MuView>
-    void computeWeightedVolume( ParticleType& particles,
-                                const NeighListType& neigh_list,
-                                const MuView& mu )
+    template <class ParticleType, class MuView>
+    void computeWeightedVolume( ParticleType& particles, const MuView& mu )
     {
         _timer.start();
 
@@ -287,16 +294,17 @@ class Force<ExecutionSpace, ForceModel<LPS, Fracture>>
         auto m = particles.sliceWeightedVolume();
         Cabana::deep_copy( m, 0.0 );
         auto model = _model;
+        auto neigh_list = _neigh_list;
 
         auto weighted_volume = KOKKOS_LAMBDA( const int i )
         {
             std::size_t num_neighbors =
-                Cabana::NeighborList<NeighListType>::numNeighbor( neigh_list,
-                                                                  i );
+                Cabana::NeighborList<neighbor_list_type>::numNeighbor(
+                    neigh_list, i );
             for ( std::size_t n = 0; n < num_neighbors; n++ )
             {
                 std::size_t j =
-                    Cabana::NeighborList<NeighListType>::getNeighbor(
+                    Cabana::NeighborList<neighbor_list_type>::getNeighbor(
                         neigh_list, i, n );
 
                 // Get the reference positions and displacements.
@@ -316,9 +324,8 @@ class Force<ExecutionSpace, ForceModel<LPS, Fracture>>
         _timer.stop();
     }
 
-    template <class ParticleType, class NeighListType, class MuView>
-    void computeDilatation( ParticleType& particles,
-                            const NeighListType& neigh_list, const MuView& mu )
+    template <class ParticleType, class MuView>
+    void computeDilatation( ParticleType& particles, const MuView& mu )
     {
         _timer.start();
 
@@ -329,17 +336,18 @@ class Force<ExecutionSpace, ForceModel<LPS, Fracture>>
         auto m = particles.sliceWeightedVolume();
         auto theta = particles.sliceDilatation();
         auto model = _model;
+        auto neigh_list = _neigh_list;
         Cabana::deep_copy( theta, 0.0 );
 
         auto dilatation = KOKKOS_LAMBDA( const int i )
         {
             std::size_t num_neighbors =
-                Cabana::NeighborList<NeighListType>::numNeighbor( neigh_list,
-                                                                  i );
+                Cabana::NeighborList<neighbor_list_type>::numNeighbor(
+                    neigh_list, i );
             for ( std::size_t n = 0; n < num_neighbors; n++ )
             {
                 std::size_t j =
-                    Cabana::NeighborList<NeighListType>::getNeighbor(
+                    Cabana::NeighborList<neighbor_list_type>::getNeighbor(
                         neigh_list, i, n );
 
                 // Get the bond distance, displacement, and stretch.
@@ -363,11 +371,10 @@ class Force<ExecutionSpace, ForceModel<LPS, Fracture>>
         _timer.stop();
     }
 
-    template <class ForceType, class PosType, class ParticleType,
-              class NeighListType, class MuView, class ParallelType>
+    template <class ForceType, class PosType, class ParticleType, class MuView,
+              class ParallelType>
     void computeForceFull( ForceType& f, const PosType& x, const PosType& u,
-                           const ParticleType& particles,
-                           const NeighListType& neigh_list, MuView& mu,
+                           const ParticleType& particles, MuView& mu,
                            const int n_local, ParallelType& )
     {
         _timer.start();
@@ -376,6 +383,7 @@ class Force<ExecutionSpace, ForceModel<LPS, Fracture>>
         auto theta_coeff = _model.theta_coeff;
         auto s_coeff = _model.s_coeff;
         auto model = _model;
+        auto neigh_list = _neigh_list;
 
         const auto vol = particles.sliceVolume();
         auto theta = particles.sliceDilatation();
@@ -384,8 +392,8 @@ class Force<ExecutionSpace, ForceModel<LPS, Fracture>>
         auto force_full = KOKKOS_LAMBDA( const int i )
         {
             std::size_t num_neighbors =
-                Cabana::NeighborList<NeighListType>::numNeighbor( neigh_list,
-                                                                  i );
+                Cabana::NeighborList<neighbor_list_type>::numNeighbor(
+                    neigh_list, i );
             for ( std::size_t n = 0; n < num_neighbors; n++ )
             {
                 double fx_i = 0.0;
@@ -393,7 +401,7 @@ class Force<ExecutionSpace, ForceModel<LPS, Fracture>>
                 double fz_i = 0.0;
 
                 std::size_t j =
-                    Cabana::NeighborList<NeighListType>::getNeighbor(
+                    Cabana::NeighborList<neighbor_list_type>::getNeighbor(
                         neigh_list, i, n );
 
                 // Get the reference positions and displacements.
@@ -437,17 +445,17 @@ class Force<ExecutionSpace, ForceModel<LPS, Fracture>>
     }
 
     template <class PosType, class WType, class DamageType, class ParticleType,
-              class NeighListType, class MuView, class ParallelType>
+              class MuView, class ParallelType>
     double computeEnergyFull( WType& W, const PosType& x, const PosType& u,
                               DamageType& phi, const ParticleType& particles,
-                              const NeighListType& neigh_list, MuView& mu,
-                              const int n_local, ParallelType& )
+                              MuView& mu, const int n_local, ParallelType& )
     {
         _energy_timer.start();
 
         auto theta_coeff = _model.theta_coeff;
         auto s_coeff = _model.s_coeff;
         auto model = _model;
+        auto neigh_list = _neigh_list;
 
         const auto vol = particles.sliceVolume();
         const auto theta = particles.sliceDilatation();
@@ -456,8 +464,8 @@ class Force<ExecutionSpace, ForceModel<LPS, Fracture>>
         auto energy_full = KOKKOS_LAMBDA( const int i, double& Phi )
         {
             std::size_t num_neighbors =
-                Cabana::NeighborList<NeighListType>::numNeighbor( neigh_list,
-                                                                  i );
+                Cabana::NeighborList<neighbor_list_type>::numNeighbor(
+                    neigh_list, i );
             double num_bonds = 0.0;
             for ( std::size_t n = 0; n < num_neighbors; n++ )
                 num_bonds += static_cast<double>( mu( i, n ) );
@@ -467,7 +475,7 @@ class Force<ExecutionSpace, ForceModel<LPS, Fracture>>
             for ( std::size_t n = 0; n < num_neighbors; n++ )
             {
                 std::size_t j =
-                    Cabana::NeighborList<NeighListType>::getNeighbor(
+                    Cabana::NeighborList<neighbor_list_type>::getNeighbor(
                         neigh_list, i, n );
                 // Get the bond distance, displacement, and stretch.
                 double xi, r, s;
@@ -498,32 +506,35 @@ class Force<ExecutionSpace, ForceModel<LPS, Fracture>>
     }
 };
 
-template <class ExecutionSpace>
-class Force<ExecutionSpace, ForceModel<LinearLPS, Elastic>>
-    : public Force<ExecutionSpace, ForceModel<LPS, Elastic>>
+template <class MemorySpace>
+class Force<MemorySpace, ForceModel<LinearLPS, Elastic>>
+    : public Force<MemorySpace, ForceModel<LPS, Elastic>>
 {
   protected:
-    using base_type = Force<ExecutionSpace, ForceModel<LPS, Elastic>>;
-    using base_type::_half_neigh;
+    using base_type = Force<MemorySpace, ForceModel<LPS, Elastic>>;
     ForceModel<LinearLPS, Elastic> _model;
 
     using base_type::_energy_timer;
     using base_type::_timer;
 
   public:
-    using exec_space = ExecutionSpace;
-
-    Force( const bool half_neigh, const ForceModel<LinearLPS, Elastic> model )
-        : base_type( half_neigh, model )
+    // Using the default exec_space.
+    using exec_space = typename MemorySpace::execution_space;
+    using neighbor_list_type = typename base_type::neighbor_list_type;
+    using base_type::_neigh_list;
+
+    template <class ParticleType>
+    Force( const bool half_neigh, ParticleType& particles,
+           const ForceModel<LinearLPS, Elastic> model )
+        : base_type( half_neigh, particles, model )
         , _model( model )
     {
     }
 
     template <class ForceType, class PosType, class ParticleType,
-              class NeighListType, class ParallelType>
+              class ParallelType>
     void computeForceFull( ForceType& f, const PosType& x, const PosType& u,
-                           const ParticleType& particles,
-                           const NeighListType& neigh_list, const int n_local,
+                           const ParticleType& particles, const int n_local,
                            ParallelType& neigh_op_tag )
     {
         _timer.start();
@@ -565,24 +576,24 @@ class Force<ExecutionSpace, ForceModel<LinearLPS, Elastic>>
 
         Kokkos::RangePolicy<exec_space> policy( 0, n_local );
         Cabana::neighbor_parallel_for(
-            policy, force_full, neigh_list, Cabana::FirstNeighborsTag(),
+            policy, force_full, _neigh_list, Cabana::FirstNeighborsTag(),
             neigh_op_tag, "CabanaPD::ForceLPS::computeFull" );
 
         _timer.stop();
     }
 
     template <class PosType, class WType, class ParticleType,
-              class NeighListType, class ParallelType>
+              class ParallelType>
     double computeEnergyFull( WType& W, const PosType& x, const PosType& u,
-                              const ParticleType& particles,
-                              const NeighListType& neigh_list,
-                              const int n_local, ParallelType& neigh_op_tag )
+                              const ParticleType& particles, const int n_local,
+                              ParallelType& neigh_op_tag )
     {
         _energy_timer.start();
 
         auto theta_coeff = _model.theta_coeff;
         auto s_coeff = _model.s_coeff;
         auto model = _model;
+        auto neigh_list = _neigh_list;
 
         const auto vol = particles.sliceVolume();
         const auto linear_theta = particles.sliceDilatation();
@@ -597,8 +608,8 @@ class Force<ExecutionSpace, ForceModel<LinearLPS, Elastic>>
             getLinearizedDistance( x, u, i, j, xi, linear_s );
 
             double num_neighbors = static_cast<double>(
-                Cabana::NeighborList<NeighListType>::numNeighbor( neigh_list,
-                                                                  i ) );
+                Cabana::NeighborList<neighbor_list_type>::numNeighbor(
+                    neigh_list, i ) );
 
             double w = ( 1.0 / num_neighbors ) * 0.5 * theta_coeff / 3.0 *
                            ( linear_theta( i ) * linear_theta( i ) ) +
@@ -613,7 +624,7 @@ class Force<ExecutionSpace, ForceModel<LinearLPS, Elastic>>
 
         Kokkos::RangePolicy<exec_space> policy( 0, n_local );
         Cabana::neighbor_parallel_reduce(
-            policy, energy_full, neigh_list, Cabana::FirstNeighborsTag(),
+            policy, energy_full, _neigh_list, Cabana::FirstNeighborsTag(),
             neigh_op_tag, strain_energy,
             "CabanaPD::ForceLPS::computeEnergyFull" );
 
diff --git a/src/force/CabanaPD_Force_PMB.hpp b/src/force/CabanaPD_Force_PMB.hpp
index 5b8cd537..07552d59 100644
--- a/src/force/CabanaPD_Force_PMB.hpp
+++ b/src/force/CabanaPD_Force_PMB.hpp
@@ -70,14 +70,17 @@
 
 namespace CabanaPD
 {
-template <class ExecutionSpace, class... ModelParams>
-class Force<ExecutionSpace, ForceModel<PMB, Elastic, ModelParams...>>
-    : public Force<ExecutionSpace, BaseForceModel>
+template <class MemorySpace, class... ModelParams>
+class Force<MemorySpace, ForceModel<PMB, Elastic, ModelParams...>>
+    : public Force<MemorySpace, BaseForceModel>
 {
   public:
-    using exec_space = ExecutionSpace;
+    // Using the default exec_space.
+    using exec_space = typename MemorySpace::execution_space;
     using model_type = ForceModel<PMB, Elastic, ModelParams...>;
-    using base_type = Force<exec_space, BaseForceModel>;
+    using base_type = Force<MemorySpace, BaseForceModel>;
+    using neighbor_list_type = typename base_type::neighbor_list_type;
+    using base_type::_neigh_list;
 
   protected:
     using base_type::_half_neigh;
@@ -87,17 +90,18 @@ class Force<ExecutionSpace, ForceModel<PMB, Elastic, ModelParams...>>
     using base_type::_timer;
 
   public:
-    Force( const bool half_neigh, const model_type model )
-        : base_type( half_neigh )
+    template <class ParticleType>
+    Force( const bool half_neigh, const ParticleType& particles,
+           const model_type model )
+        : base_type( half_neigh, model.delta, particles )
         , _model( model )
     {
     }
 
     template <class ForceType, class PosType, class ParticleType,
-              class NeighListType, class ParallelType>
+              class ParallelType>
     void computeForceFull( ForceType& f, const PosType& x, const PosType& u,
-                           const ParticleType& particles,
-                           const NeighListType& neigh_list, const int n_local,
+                           const ParticleType& particles, const int n_local,
                            ParallelType& neigh_op_tag )
     {
         _timer.start();
@@ -129,18 +133,17 @@ class Force<ExecutionSpace, ForceModel<PMB, Elastic, ModelParams...>>
 
         Kokkos::RangePolicy<exec_space> policy( 0, n_local );
         Cabana::neighbor_parallel_for(
-            policy, force_full, neigh_list, Cabana::FirstNeighborsTag(),
+            policy, force_full, _neigh_list, Cabana::FirstNeighborsTag(),
             neigh_op_tag, "CabanaPD::ForcePMB::computeFull" );
 
         _timer.stop();
     }
 
     template <class PosType, class WType, class ParticleType,
-              class NeighListType, class ParallelType>
+              class ParallelType>
     double computeEnergyFull( WType& W, const PosType& x, const PosType& u,
-                              const ParticleType& particles,
-                              const NeighListType& neigh_list,
-                              const int n_local, ParallelType& neigh_op_tag )
+                              const ParticleType& particles, const int n_local,
+                              ParallelType& neigh_op_tag )
     {
         _energy_timer.start();
 
@@ -166,7 +169,7 @@ class Force<ExecutionSpace, ForceModel<PMB, Elastic, ModelParams...>>
         double strain_energy = 0.0;
         Kokkos::RangePolicy<exec_space> policy( 0, n_local );
         Cabana::neighbor_parallel_reduce(
-            policy, energy_full, neigh_list, Cabana::FirstNeighborsTag(),
+            policy, energy_full, _neigh_list, Cabana::FirstNeighborsTag(),
             neigh_op_tag, strain_energy,
             "CabanaPD::ForcePMB::computeEnergyFull" );
 
@@ -175,17 +178,20 @@ class Force<ExecutionSpace, ForceModel<PMB, Elastic, ModelParams...>>
     }
 };
 
-template <class ExecutionSpace, class... ModelParams>
-class Force<ExecutionSpace, ForceModel<PMB, Fracture, ModelParams...>>
-    : public Force<ExecutionSpace, BaseForceModel>
+template <class MemorySpace, class... ModelParams>
+class Force<MemorySpace, ForceModel<PMB, Fracture, ModelParams...>>
+    : public Force<MemorySpace, BaseForceModel>
 {
   public:
-    using exec_space = ExecutionSpace;
+    // Using the default exec_space.
+    using exec_space = typename MemorySpace::execution_space;
     using model_type = ForceModel<PMB, Fracture, ModelParams...>;
+    using base_type = Force<MemorySpace, BaseForceModel>;
+    using neighbor_list_type = typename base_type::neighbor_list_type;
+    using base_type::_neigh_list;
 
   protected:
     using base_model_type = typename model_type::base_type;
-    using base_type = Force<ExecutionSpace, BaseForceModel>;
     using base_type::_half_neigh;
     model_type _model;
 
@@ -193,30 +199,32 @@ class Force<ExecutionSpace, ForceModel<PMB, Fracture, ModelParams...>>
     using base_type::_timer;
 
   public:
-    Force( const bool half_neigh, const model_type model )
-        : base_type( half_neigh )
+    template <class ParticleType>
+    Force( const bool half_neigh, const ParticleType& particles,
+           const model_type model )
+        : base_type( half_neigh, model.delta, particles )
         , _model( model )
     {
     }
 
-    template <class ForceType, class PosType, class ParticleType,
-              class NeighListType, class MuView, class ParallelType>
+    template <class ForceType, class PosType, class ParticleType, class MuView,
+              class ParallelType>
     void computeForceFull( ForceType& f, const PosType& x, const PosType& u,
-                           const ParticleType& particles,
-                           const NeighListType& neigh_list, MuView& mu,
+                           const ParticleType& particles, MuView& mu,
                            const int n_local, ParallelType& )
     {
         _timer.start();
 
         auto model = _model;
+        auto neigh_list = _neigh_list;
         const auto vol = particles.sliceVolume();
         const auto nofail = particles.sliceNoFail();
 
         auto force_full = KOKKOS_LAMBDA( const int i )
         {
             std::size_t num_neighbors =
-                Cabana::NeighborList<NeighListType>::numNeighbor( neigh_list,
-                                                                  i );
+                Cabana::NeighborList<neighbor_list_type>::numNeighbor(
+                    neigh_list, i );
             for ( std::size_t n = 0; n < num_neighbors; n++ )
             {
                 double fx_i = 0.0;
@@ -224,7 +232,7 @@ class Force<ExecutionSpace, ForceModel<PMB, Fracture, ModelParams...>>
                 double fz_i = 0.0;
 
                 std::size_t j =
-                    Cabana::NeighborList<NeighListType>::getNeighbor(
+                    Cabana::NeighborList<neighbor_list_type>::getNeighbor(
                         neigh_list, i, n );
 
                 // Get the reference positions and displacements.
@@ -264,28 +272,28 @@ class Force<ExecutionSpace, ForceModel<PMB, Fracture, ModelParams...>>
     }
 
     template <class PosType, class WType, class DamageType, class ParticleType,
-              class NeighListType, class MuView, class ParallelType>
+              class MuView, class ParallelType>
     double computeEnergyFull( WType& W, const PosType& x, const PosType& u,
                               DamageType& phi, const ParticleType& particles,
-                              const NeighListType& neigh_list, MuView& mu,
-                              const int n_local, ParallelType& )
+                              MuView& mu, const int n_local, ParallelType& )
     {
         _energy_timer.start();
 
         auto model = _model;
+        auto neigh_list = _neigh_list;
         const auto vol = particles.sliceVolume();
 
         auto energy_full = KOKKOS_LAMBDA( const int i, double& Phi )
         {
             std::size_t num_neighbors =
-                Cabana::NeighborList<NeighListType>::numNeighbor( neigh_list,
-                                                                  i );
+                Cabana::NeighborList<neighbor_list_type>::numNeighbor(
+                    neigh_list, i );
             double phi_i = 0.0;
             double vol_H_i = 0.0;
             for ( std::size_t n = 0; n < num_neighbors; n++ )
             {
                 std::size_t j =
-                    Cabana::NeighborList<NeighListType>::getNeighbor(
+                    Cabana::NeighborList<neighbor_list_type>::getNeighbor(
                         neigh_list, i, n );
                 // Get the bond distance, displacement, and stretch.
                 double xi, r, s;
@@ -315,16 +323,19 @@ class Force<ExecutionSpace, ForceModel<PMB, Fracture, ModelParams...>>
     }
 };
 
-template <class ExecutionSpace, class... ModelParams>
-class Force<ExecutionSpace, ForceModel<LinearPMB, Elastic, ModelParams...>>
-    : public Force<ExecutionSpace, BaseForceModel>
+template <class MemorySpace, class... ModelParams>
+class Force<MemorySpace, ForceModel<LinearPMB, Elastic, ModelParams...>>
+    : public Force<MemorySpace, BaseForceModel>
 {
   public:
-    using exec_space = ExecutionSpace;
+    // Using the default exec_space.
+    using exec_space = typename MemorySpace::execution_space;
     using model_type = ForceModel<LinearPMB, Elastic, TemperatureIndependent>;
+    using base_type = Force<MemorySpace, BaseForceModel>;
+    using neighbor_list_type = typename base_type::neighbor_list_type;
+    using base_type::_neigh_list;
 
   protected:
-    using base_type = Force<ExecutionSpace, BaseForceModel>;
     using base_type::_half_neigh;
     model_type _model;
 
@@ -332,17 +343,18 @@ class Force<ExecutionSpace, ForceModel<LinearPMB, Elastic, ModelParams...>>
     using base_type::_timer;
 
   public:
-    Force( const bool half_neigh, const model_type model )
-        : base_type( half_neigh )
+    template <class ParticleType>
+    Force( const bool half_neigh, const ParticleType& particles,
+           const model_type model )
+        : base_type( half_neigh, model.delta, particles )
         , _model( model )
     {
     }
 
     template <class ForceType, class PosType, class ParticleType,
-              class NeighListType, class ParallelType>
+              class ParallelType>
     void computeForceFull( ForceType& f, const PosType& x, const PosType& u,
-                           ParticleType& particles,
-                           const NeighListType& neigh_list, const int n_local,
+                           ParticleType& particles, const int n_local,
                            ParallelType& neigh_op_tag )
     {
         _timer.start();
@@ -376,18 +388,17 @@ class Force<ExecutionSpace, ForceModel<LinearPMB, Elastic, ModelParams...>>
 
         Kokkos::RangePolicy<exec_space> policy( 0, n_local );
         Cabana::neighbor_parallel_for(
-            policy, force_full, neigh_list, Cabana::FirstNeighborsTag(),
+            policy, force_full, _neigh_list, Cabana::FirstNeighborsTag(),
             neigh_op_tag, "CabanaPD::ForceLinearPMB::computeFull" );
 
         _timer.stop();
     }
 
     template <class PosType, class WType, class ParticleType,
-              class NeighListType, class ParallelType>
+              class ParallelType>
     double computeEnergyFull( WType& W, const PosType& x, const PosType& u,
-                              ParticleType& particles,
-                              const NeighListType& neigh_list,
-                              const int n_local, ParallelType& neigh_op_tag )
+                              ParticleType& particles, const int n_local,
+                              ParallelType& neigh_op_tag )
     {
         _energy_timer.start();
 
@@ -413,7 +424,7 @@ class Force<ExecutionSpace, ForceModel<LinearPMB, Elastic, ModelParams...>>
         double strain_energy = 0.0;
         Kokkos::RangePolicy<exec_space> policy( 0, n_local );
         Cabana::neighbor_parallel_reduce(
-            policy, energy_full, neigh_list, Cabana::FirstNeighborsTag(),
+            policy, energy_full, _neigh_list, Cabana::FirstNeighborsTag(),
             neigh_op_tag, strain_energy,
             "CabanaPD::ForceLinearPMB::computeEnergyFull" );
 
diff --git a/unit_test/tstForce.hpp b/unit_test/tstForce.hpp
index 8e7d2a10..30dc8e1f 100644
--- a/unit_test/tstForce.hpp
+++ b/unit_test/tstForce.hpp
@@ -623,52 +623,44 @@ struct NoDamageTag
 {
 };
 
-template <class ForceType, class ParticleType, class NeighborList>
+template <class ForceType, class ParticleType>
 double computeEnergyAndForce( NoDamageTag, ForceType force,
-                              ParticleType& particles,
-                              const NeighborList& neigh_list, const int )
+                              ParticleType& particles, const int )
 {
-    computeForce( force, particles, neigh_list, Cabana::SerialOpTag() );
-    double Phi =
-        computeEnergy( force, particles, neigh_list, Cabana::SerialOpTag() );
+    computeForce( force, particles, Cabana::SerialOpTag() );
+    double Phi = computeEnergy( force, particles, Cabana::SerialOpTag() );
     return Phi;
 }
-template <class ForceType, class ParticleType, class NeighborList>
-double
-computeEnergyAndForce( DamageTag, ForceType force, ParticleType& particles,
-                       const NeighborList& neigh_list, const int max_neighbors )
+template <class ForceType, class ParticleType>
+double computeEnergyAndForce( DamageTag, ForceType force,
+                              ParticleType& particles, const int max_neighbors )
 {
     Kokkos::View<int**, TEST_MEMSPACE> mu(
         Kokkos::ViewAllocateWithoutInitializing( "broken_bonds" ),
         particles.n_local, max_neighbors );
     Kokkos::deep_copy( mu, 1 );
-    computeForce( force, particles, neigh_list, mu, Cabana::SerialOpTag() );
-    double Phi = computeEnergy( force, particles, neigh_list, mu,
-                                Cabana::SerialOpTag() );
+    computeForce( force, particles, mu, Cabana::SerialOpTag() );
+    double Phi = computeEnergy( force, particles, mu, Cabana::SerialOpTag() );
     return Phi;
 }
 
-template <class ModelType, class ForceType, class ParticleType,
-          class NeighborList>
-void initializeForce( ModelType, ForceType& force, ParticleType& particles,
-                      const NeighborList& neigh_list )
+template <class ModelType, class ForceType, class ParticleType>
+void initializeForce( ModelType, ForceType& force, ParticleType& particles )
 {
-    force.computeWeightedVolume( particles, neigh_list, Cabana::SerialOpTag() );
-    force.computeDilatation( particles, neigh_list, Cabana::SerialOpTag() );
+    force.computeWeightedVolume( particles, Cabana::SerialOpTag() );
+    force.computeDilatation( particles, Cabana::SerialOpTag() );
 }
 
-template <class ForceType, class ParticleType, class NeighborList>
+template <class ForceType, class ParticleType>
 void initializeForce( CabanaPD::ForceModel<CabanaPD::LPS, CabanaPD::Fracture>,
-                      ForceType& force, ParticleType& particles,
-                      const NeighborList& neigh_list )
+                      ForceType& force, ParticleType& particles )
 {
-    int max_neighbors =
-        Cabana::NeighborList<NeighborList>::maxNeighbor( neigh_list );
+    auto max_neighbors = force.getMaxLocalNeighbors();
     Kokkos::View<int**, TEST_MEMSPACE> mu( "broken_bonds", particles.n_local,
                                            max_neighbors );
     Kokkos::deep_copy( mu, 1 );
-    force.computeWeightedVolume( particles, neigh_list, mu );
-    force.computeDilatation( particles, neigh_list, mu );
+    force.computeWeightedVolume( particles, mu );
+    force.computeDilatation( particles, mu );
 }
 
 template <class ParticleType, class AoSoAType>
@@ -698,33 +690,21 @@ void testForce( ModelType model, const DamageType damage_tag, const double dx,
 
     // This needs to exactly match the mesh spacing to compare with the single
     // particle calculation.
-    CabanaPD::Force<TEST_EXECSPACE, ModelType> force( true, model );
-
-    double mesh_min[3] = { particles.ghost_mesh_lo[0],
-                           particles.ghost_mesh_lo[1],
-                           particles.ghost_mesh_lo[2] };
-    double mesh_max[3] = { particles.ghost_mesh_hi[0],
-                           particles.ghost_mesh_hi[1],
-                           particles.ghost_mesh_hi[2] };
-    using verlet_list =
-        Cabana::VerletList<TEST_MEMSPACE, Cabana::FullNeighborTag,
-                           Cabana::VerletLayout2D, Cabana::TeamOpTag>;
-    // Add to delta to make sure neighbors are found.
-    auto x = particles.sliceReferencePosition();
-    verlet_list neigh_list( x, 0, particles.n_local, model.delta + 1e-14, 1.0,
-                            mesh_min, mesh_max );
-    int max_neighbors =
-        Cabana::NeighborList<verlet_list>::maxNeighbor( neigh_list );
+    CabanaPD::Force<TEST_MEMSPACE, ModelType> force( true, particles, model );
 
+    auto x = particles.sliceReferencePosition();
     auto f = particles.sliceForce();
     auto W = particles.sliceStrainEnergy();
     auto vol = particles.sliceVolume();
     //  No communication needed (as in the main solver) since this test is only
     //  intended for one rank.
-    initializeForce( model, force, particles, neigh_list );
+    initializeForce( model, force, particles );
 
-    double Phi = computeEnergyAndForce( damage_tag, force, particles,
-                                        neigh_list, max_neighbors );
+    unsigned int max_neighbors;
+    unsigned long long total_neighbors;
+    force.getNeighborStatistics( max_neighbors, total_neighbors );
+    double Phi =
+        computeEnergyAndForce( damage_tag, force, particles, max_neighbors );
 
     // Make a copy of final results on the host
     std::size_t num_particle = x.size();