diff --git a/examples/mechanics/fragmenting_cylinder.cpp b/examples/mechanics/fragmenting_cylinder.cpp
index 22b7ce6f..afd4d463 100644
--- a/examples/mechanics/fragmenting_cylinder.cpp
+++ b/examples/mechanics/fragmenting_cylinder.cpp
@@ -139,7 +139,7 @@ void fragmentingCylinderExample( const std::string filename )
     {
         double r_c = inputs["contact_horizon_factor"];
         r_c *= dx[0];
-        CabanaPD::NormalRepulsionModel contact_model( delta, r_c, K );
+        CabanaPD::NormalRepulsionModel contact_model( delta, r_c, r_c, K );
 
         auto cabana_pd = CabanaPD::createSolver<memory_space>(
             inputs, particles, force_model, contact_model );
diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp
index ecd48708..35e928f8 100644
--- a/src/CabanaPD_Force.hpp
+++ b/src/CabanaPD_Force.hpp
@@ -179,6 +179,12 @@ class Force<MemorySpace, BaseForceModel>
     {
     }
 
+    // FIXME: should it be possible to update this list?
+    template <class ParticleType>
+    void update( const ParticleType&, const double, const bool = false )
+    {
+    }
+
     unsigned getMaxLocalNeighbors()
     {
         auto neigh = _neigh_list;
@@ -234,6 +240,7 @@ class Force<MemorySpace, BaseForceModel>
 
     auto time() { return _timer.time(); };
     auto timeEnergy() { return _energy_timer.time(); };
+    auto timeNeighbor() { return 0.0; };
 };
 
 template <class MemorySpace>
@@ -276,6 +283,7 @@ class BaseFracture
 ******************************************************************************/
 template <class ForceType, class ParticleType, class ParallelType>
 void computeForce( ForceType& force, ParticleType& particles,
+                   const double max_displacement,
                    const ParallelType& neigh_op_tag, const bool reset = true )
 {
     auto x = particles.sliceReferencePosition();
@@ -294,9 +302,11 @@ void computeForce( ForceType& force, ParticleType& particles,
 
     // 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_op_tag );
+        force.computeForceFull( f_a, x, u, particles, max_displacement,
+                                neigh_op_tag );
     else
-        force.computeForceFull( f, x, u, particles, neigh_op_tag );
+        force.computeForceFull( f, x, u, particles, max_displacement,
+                                neigh_op_tag );
     Kokkos::fence();
 }
 
diff --git a/src/CabanaPD_Integrate.hpp b/src/CabanaPD_Integrate.hpp
index a0190878..c8ffbf13 100644
--- a/src/CabanaPD_Integrate.hpp
+++ b/src/CabanaPD_Integrate.hpp
@@ -85,7 +85,7 @@ class Integrator
     ~Integrator() {}
 
     template <class ParticlesType>
-    void initialHalfStep( ParticlesType& p )
+    double initialHalfStep( ParticlesType& p )
     {
         _timer.start();
 
@@ -93,10 +93,11 @@ class Integrator
         auto v = p.sliceVelocity();
         auto f = p.sliceForce();
         auto rho = p.sliceDensity();
+        auto u_neigh = p.sliceDisplacementNeighborBuild();
 
         auto dt = _dt;
         auto half_dt = _half_dt;
-        auto init_func = KOKKOS_LAMBDA( const int i )
+        auto init_func = KOKKOS_LAMBDA( const int i, double& max_u )
         {
             const double half_dt_m = half_dt / rho( i );
             v( i, 0 ) += half_dt_m * f( i, 0 );
@@ -105,13 +106,26 @@ class Integrator
             u( i, 0 ) += dt * v( i, 0 );
             u( i, 1 ) += dt * v( i, 1 );
             u( i, 2 ) += dt * v( i, 2 );
+
+            // FIXME: only used for contact: get displacement since last
+            // neighbor update.
+            u_neigh( i, 0 ) += u( i, 0 );
+            u_neigh( i, 1 ) += u( i, 1 );
+            u_neigh( i, 2 ) += u( i, 2 );
+            auto u_mag = Kokkos::hypot( u_neigh( i, 0 ), u_neigh( i, 1 ),
+                                        u_neigh( i, 2 ) );
+            if ( u_mag > max_u )
+                max_u = u_mag;
         };
         Kokkos::RangePolicy<exec_space> policy( p.frozenOffset(),
                                                 p.localOffset() );
-        Kokkos::parallel_for( "CabanaPD::Integrator::Initial", policy,
-                              init_func );
+        double max_displacement;
+        Kokkos::parallel_reduce( "CabanaPD::Integrator::Initial", policy,
+                                 init_func,
+                                 Kokkos::Max<double>( max_displacement ) );
 
         _timer.stop();
+        return max_displacement;
     }
 
     template <class ParticlesType>
diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp
index 08d69a7b..2c91c904 100644
--- a/src/CabanaPD_Particles.hpp
+++ b/src/CabanaPD_Particles.hpp
@@ -108,7 +108,7 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
     // no-fail.
     using int_type = Cabana::MemberTypes<int>;
     // v, rho, type.
-    using other_types = Cabana::MemberTypes<double[dim], double, int>;
+    using other_types = Cabana::MemberTypes<double[dim], double>;
 
     // FIXME: add vector length.
     // FIXME: enable variable aosoa.
@@ -116,6 +116,7 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
     using aosoa_y_type = Cabana::AoSoA<vector_type, memory_space, 1>;
     using aosoa_vol_type = Cabana::AoSoA<scalar_type, memory_space, 1>;
     using aosoa_nofail_type = Cabana::AoSoA<int_type, memory_space, 1>;
+    using aosoa_u_neigh_type = Cabana::AoSoA<vector_type, memory_space, 1>;
     using aosoa_other_type = Cabana::AoSoA<other_types, memory_space>;
     // Using grid here for the particle init.
     using plist_x_type =
@@ -284,11 +285,11 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
         auto x = sliceReferencePosition();
         auto v = sliceVelocity();
         auto f = sliceForce();
-        auto type = sliceType();
+        // auto type = sliceType();
         auto rho = sliceDensity();
         auto u = sliceDisplacement();
         auto vol = sliceVolume();
-        auto nofail = sliceNoFail();
+        auto u_neigh = sliceDisplacementNeighborBuild();
 
         // Initialize particles.
         auto create_functor =
@@ -307,6 +308,7 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
                 Cabana::get( particle, CabanaPD::Field::ReferencePosition(),
                              d ) = px[d];
                 u( pid, d ) = 0.0;
+                u_neigh( pid, d ) = 0;
                 v( pid, d ) = 0.0;
                 f( pid, d ) = 0.0;
             }
@@ -314,8 +316,7 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
             vol( pid ) = pv;
 
             // FIXME: hardcoded.
-            type( pid ) = 0;
-            nofail( pid ) = 0;
+            // type( pid ) = 0;
             rho( pid ) = 1.0;
 
             return create;
@@ -353,10 +354,10 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
         auto p_vol = sliceVolume();
         auto v = sliceVelocity();
         auto f = sliceForce();
-        auto type = sliceType();
+        // auto type = sliceType();
         auto rho = sliceDensity();
         auto u = sliceDisplacement();
-        auto nofail = sliceNoFail();
+        auto u_neigh = sliceDisplacementNeighborBuild();
 
         static_assert(
             Cabana::is_accessible_from<
@@ -374,11 +375,11 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
                 {
                     p_x( pid, d ) = x( pid_offset, d );
                     u( pid, d ) = 0.0;
+                    u_neigh( pid, d ) = 0;
                     v( pid, d ) = 0.0;
                     f( pid, d ) = 0.0;
                 }
-                type( pid ) = 0;
-                nofail( pid ) = 0;
+                // type( pid ) = 0;
                 rho( pid ) = 1.0;
             } );
 
@@ -447,6 +448,14 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
     {
         return Cabana::slice<0>( _aosoa_u, "displacements" );
     }
+    auto sliceDisplacementNeighborBuild()
+    {
+        return Cabana::slice<0>( _aosoa_u_neigh, "displacement_since_rebuild" );
+    }
+    auto sliceDisplacementNeighborBuild() const
+    {
+        return Cabana::slice<0>( _aosoa_u_neigh, "displacement_since_rebuild" );
+    }
     auto sliceForce() { return _plist_f.slice( CabanaPD::Field::Force() ); }
     auto sliceForceAtomic()
     {
@@ -474,8 +483,9 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
     {
         return Cabana::slice<1>( _aosoa_other, "density" );
     }
-    auto sliceType() { return Cabana::slice<2>( _aosoa_other, "type" ); }
-    auto sliceType() const { return Cabana::slice<2>( _aosoa_other, "type" ); }
+    // auto sliceType() { return Cabana::slice<2>( _aosoa_other, "type" ); }
+    // auto sliceType() const { return Cabana::slice<2>( _aosoa_other, "type" );
+    // }
 
     auto sliceNoFail()
     {
@@ -524,10 +534,45 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
         _plist_f.aosoa().resize( localOffset() );
         _aosoa_other.resize( localOffset() );
         _aosoa_nofail.resize( referenceOffset() );
+        _aosoa_u_neigh.resize( localOffset() );
+
         size = _plist_x.size();
         _timer.stop();
     };
 
+    void shrink()
+    {
+        _timer.start();
+        _plist_x.aosoa().shrinkToFit();
+        _aosoa_u.shrinkToFit();
+        _aosoa_y.shrinkToFit();
+        _aosoa_vol.shrinkToFit();
+        _plist_f.aosoa().shrinkToFit();
+        _aosoa_other.shrinkToFit();
+        _aosoa_u_neigh.shrinkToFit();
+        _timer.stop();
+    };
+
+    template <typename KeepType>
+    void remove( const int num_keep, const KeepType& keep )
+    {
+        Cabana::remove( execution_space(), num_keep, keep, _plist_x.aosoa(),
+                        numFrozen() );
+        Cabana::remove( execution_space(), num_keep, keep, _plist_f.aosoa(),
+                        numFrozen() );
+        Cabana::remove( execution_space(), num_keep, keep, _aosoa_u,
+                        numFrozen() );
+        Cabana::remove( execution_space(), num_keep, keep, _aosoa_vol,
+                        numFrozen() );
+        Cabana::remove( execution_space(), num_keep, keep, _aosoa_y,
+                        numFrozen() );
+        Cabana::remove( execution_space(), num_keep, keep, _aosoa_other,
+                        numFrozen() );
+        Cabana::remove( execution_space(), num_keep, keep, _aosoa_u_neigh,
+                        numFrozen() );
+        resize( frozen_offset + num_keep, 0 );
+    }
+
     auto getPosition( const bool use_reference )
     {
         if ( use_reference )
@@ -577,6 +622,7 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
 
   protected:
     aosoa_u_type _aosoa_u;
+    aosoa_u_neigh_type _aosoa_u_neigh;
     aosoa_y_type _aosoa_y;
     aosoa_vol_type _aosoa_vol;
     aosoa_nofail_type _aosoa_nofail;
diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp
index 8092f966..1c9aaaf7 100644
--- a/src/CabanaPD_Solver.hpp
+++ b/src/CabanaPD_Solver.hpp
@@ -217,7 +217,7 @@ class Solver
             comm->gatherWeightedVolume();
         }
         // Compute initial internal forces and energy.
-        updateForce();
+        updateForce( 0.0 );
         computeEnergy( *force, *particles, neigh_iter_tag() );
 
         if ( initial_output )
@@ -266,6 +266,8 @@ class Solver
         init( boundary_condition, initial_output );
     }
 
+    void updateNeighbors() { force->update( *particles, 0.0, true ); }
+
     template <typename BoundaryType>
     void run( BoundaryType boundary_condition )
     {
@@ -277,7 +279,7 @@ class Solver
             _step_timer.start();
 
             // Integrate - velocity Verlet first half.
-            integrator->initialHalfStep( *particles );
+            auto max_displacement = integrator->initialHalfStep( *particles );
 
             // Update ghost particles.
             comm->gatherDisplacement();
@@ -300,10 +302,11 @@ class Solver
                 comm->gatherTemperature();
 
             // Compute internal forces.
-            updateForce();
+            updateForce( max_displacement );
 
             if constexpr ( is_contact<contact_model_type>::value )
-                computeForce( *contact, *particles, neigh_iter_tag{}, false );
+                computeForce( *contact, *particles, max_displacement,
+                              neigh_iter_tag{}, false );
 
             // Add force boundary condition.
             if ( boundary_condition.forceUpdate() )
@@ -329,16 +332,17 @@ class Solver
             _step_timer.start();
 
             // Integrate - velocity Verlet first half.
-            integrator->initialHalfStep( *particles );
+            auto max_displacement = integrator->initialHalfStep( *particles );
 
             // Update ghost particles.
             comm->gatherDisplacement();
 
             // Compute internal forces.
-            updateForce();
+            updateForce( max_displacement );
 
             if constexpr ( is_contact<contact_model_type>::value )
-                computeForce( *contact, *particles, neigh_iter_tag{}, false );
+                computeForce( *contact, *particles, max_displacement,
+                              neigh_iter_tag{}, false );
 
             if constexpr ( is_temperature_dependent<
                                typename force_model_type::thermal_type>::value )
@@ -356,7 +360,7 @@ class Solver
 
     // Compute and communicate fields needed for force computation and update
     // forces.
-    void updateForce()
+    void updateForce( const double max_displacement )
     {
         // Compute and communicate weighted volume for LPS (does nothing for
         // PMB). Only computed once without fracture.
@@ -371,7 +375,7 @@ class Solver
         comm->gatherDilatation();
 
         // Compute internal forces.
-        computeForce( *force, *particles, neigh_iter_tag{} );
+        computeForce( *force, *particles, max_displacement, neigh_iter_tag{} );
     }
 
     void output( const int step )
@@ -419,6 +423,7 @@ class Solver
             double integrate_time = integrator->time();
             double force_time = force->time();
             double energy_time = force->timeEnergy();
+            double neigh_time = force->timeNeighbor();
             double output_time = particles->timeOutput();
             _total_time += step_time;
             auto rate = static_cast<double>( particles->numGlobal() *
@@ -428,7 +433,7 @@ class Solver
                  " ", std::scientific, std::setprecision( 2 ), step * dt, " ",
                  W, " ", std::fixed, _total_time, " ", force_time, " ",
                  comm_time, " ", integrate_time, " ", energy_time, " ",
-                 output_time, " ", std::scientific, rate );
+                 neigh_time, " ", output_time, " ", std::scientific, rate );
             out.close();
         }
     }
diff --git a/src/force/CabanaPD_ContactModels.hpp b/src/force/CabanaPD_ContactModels.hpp
index aff456f1..e35eda81 100644
--- a/src/force/CabanaPD_ContactModels.hpp
+++ b/src/force/CabanaPD_ContactModels.hpp
@@ -25,18 +25,23 @@ namespace CabanaPD
 ******************************************************************************/
 struct ContactModel
 {
+    // PD neighbor search radius.
     double delta;
-    double Rc;
+    // Contact neighbor search radius.
+    double radius;
+    // Extend neighbor search radius to reuse lists.
+    double radius_extend;
 
     // PD horizon
     // Contact radius
-    ContactModel( const double _delta, const double _Rc )
+    ContactModel( const double _delta, const double _radius,
+                  const double _radius_extend )
         : delta( _delta )
-        , Rc( _Rc ){};
+        , radius( _radius )
+        , radius_extend( _radius_extend ){};
 };
 
 /* Normal repulsion */
-
 struct NormalRepulsionModel : public ContactModel
 {
     // FIXME: This is for use as the primary force model.
@@ -45,13 +50,15 @@ struct NormalRepulsionModel : public ContactModel
     using thermal_type = TemperatureIndependent;
 
     using ContactModel::delta;
-    using ContactModel::Rc;
+    using ContactModel::radius;
+    using ContactModel::radius_extend;
 
     double c;
     double K;
 
-    NormalRepulsionModel( const double delta, const double Rc, const double _K )
-        : ContactModel( delta, Rc )
+    NormalRepulsionModel( const double delta, const double radius,
+                          const double radius_extend, const double _K )
+        : ContactModel( delta, radius, radius_extend )
         , K( _K )
     {
         K = _K;
@@ -63,7 +70,7 @@ struct NormalRepulsionModel : public ContactModel
     auto forceCoeff( const double r, const double vol ) const
     {
         // Contact "stretch"
-        const double sc = ( r - Rc ) / delta;
+        const double sc = ( r - radius ) / delta;
         // Normal repulsion uses a 15 factor compared to the PMB force
         return 15.0 * c * sc * vol;
     }
diff --git a/src/force/CabanaPD_Force_Contact.hpp b/src/force/CabanaPD_Force_Contact.hpp
index 9f299f5f..0af627bf 100644
--- a/src/force/CabanaPD_Force_Contact.hpp
+++ b/src/force/CabanaPD_Force_Contact.hpp
@@ -37,36 +37,87 @@ KOKKOS_INLINE_FUNCTION void getRelativeNormalVelocityComponents(
     vn /= r;
 };
 
+// Contact forces base class.
+template <class MemorySpace>
+class BaseForceContact : public Force<MemorySpace, BaseForceModel>
+{
+  public:
+    using base_type = Force<MemorySpace, BaseForceModel>;
+    using neighbor_list_type = typename base_type::neighbor_list_type;
+
+    template <class ParticleType, class ModelType>
+    BaseForceContact( const bool half_neigh, const ParticleType& particles,
+                      const ModelType model )
+        : base_type( half_neigh, model.radius + model.radius_extend,
+                     particles.sliceCurrentPosition(), particles.frozenOffset(),
+                     particles.localOffset(), particles.ghost_mesh_lo,
+                     particles.ghost_mesh_hi )
+        , radius( model.radius )
+        , radius_extend( model.radius_extend )
+    {
+        for ( int d = 0; d < particles.dim; d++ )
+        {
+            mesh_min[d] = particles.ghost_mesh_lo[d];
+            mesh_max[d] = particles.ghost_mesh_hi[d];
+        }
+    }
+
+    // Only rebuild neighbor list as needed.
+    template <class ParticleType>
+    void update( const ParticleType& particles, const double max_displacement,
+                 const bool require_update = false )
+    {
+        if ( max_displacement > radius_extend || require_update )
+        {
+            _neigh_timer.start();
+            const auto y = particles.sliceCurrentPosition();
+            _neigh_list.build( y, particles.frozenOffset(),
+                               particles.localOffset(), radius + radius_extend,
+                               1.0, mesh_min, mesh_max );
+            // Reset neighbor update displacement.
+            const auto u_neigh = particles.sliceDisplacementNeighborBuild();
+            Cabana::deep_copy( u_neigh, 0.0 );
+            _neigh_timer.stop();
+        }
+    }
+
+    auto timeNeighbor() { return _neigh_timer.time(); };
+
+  protected:
+    double radius;
+    double radius_extend;
+    Timer _neigh_timer;
+
+    using base_type::_half_neigh;
+    using base_type::_neigh_list;
+    double mesh_max[3];
+    double mesh_min[3];
+};
+
 /******************************************************************************
   Normal repulsion forces
 ******************************************************************************/
 template <class MemorySpace>
 class Force<MemorySpace, NormalRepulsionModel>
-    : public Force<MemorySpace, BaseForceModel>
+    : public BaseForceContact<MemorySpace>
 {
   public:
-    using base_type = Force<MemorySpace, BaseForceModel>;
+    using base_type = BaseForceContact<MemorySpace>;
     using neighbor_list_type = typename base_type::neighbor_list_type;
 
     template <class ParticleType>
     Force( const bool half_neigh, const ParticleType& particles,
            const NormalRepulsionModel model )
-        : base_type( half_neigh, model.Rc, particles.sliceCurrentPosition(),
-                     particles.frozenOffset(), particles.localOffset(),
-                     particles.ghost_mesh_lo, particles.ghost_mesh_hi )
+        : base_type( half_neigh, particles, model )
         , _model( model )
     {
-        for ( int d = 0; d < particles.dim; d++ )
-        {
-            mesh_min[d] = particles.ghost_mesh_lo[d];
-            mesh_max[d] = particles.ghost_mesh_hi[d];
-        }
     }
 
     template <class ForceType, class PosType, class ParticleType,
               class ParallelType>
     void computeForceFull( ForceType& fc, const PosType& x, const PosType& u,
                            const ParticleType& particles,
+                           const double max_displacement,
                            ParallelType& neigh_op_tag )
     {
         auto model = _model;
@@ -75,10 +126,7 @@ class Force<MemorySpace, NormalRepulsionModel>
         const int n_frozen = particles.frozenOffset();
         const int n_local = particles.localOffset();
 
-        _neigh_timer.start();
-        _neigh_list.build( y, n_frozen, n_local, model.Rc, 1.0, mesh_min,
-                           mesh_max );
-        _neigh_timer.stop();
+        base_type::update( particles, max_displacement );
 
         auto contact_full = KOKKOS_LAMBDA( const int i, const int j )
         {
@@ -126,7 +174,6 @@ class Force<MemorySpace, NormalRepulsionModel>
     using base_type::_half_neigh;
     using base_type::_neigh_list;
     using base_type::_timer;
-    Timer _neigh_timer;
 
     double mesh_max[3];
     double mesh_min[3];
diff --git a/src/force/CabanaPD_Force_HertzianContact.hpp b/src/force/CabanaPD_Force_HertzianContact.hpp
index 19712967..43b8e614 100644
--- a/src/force/CabanaPD_Force_HertzianContact.hpp
+++ b/src/force/CabanaPD_Force_HertzianContact.hpp
@@ -26,32 +26,25 @@ namespace CabanaPD
   Normal repulsion forces
 ******************************************************************************/
 template <class MemorySpace>
-class Force<MemorySpace, HertzianModel>
-    : public Force<MemorySpace, BaseForceModel>
+class Force<MemorySpace, HertzianModel> : public BaseForceContact<MemorySpace>
 {
   public:
-    using base_type = Force<MemorySpace, BaseForceModel>;
+    using base_type = BaseForceContact<MemorySpace>;
     using neighbor_list_type = typename base_type::neighbor_list_type;
 
     template <class ParticleType>
     Force( const bool half_neigh, const ParticleType& particles,
            const HertzianModel model )
-        : base_type( half_neigh, model.Rc, particles.sliceCurrentPosition(),
-                     particles.frozenOffset(), particles.localOffset(),
-                     particles.ghost_mesh_lo, particles.ghost_mesh_hi )
+        : base_type( half_neigh, particles, model )
         , _model( model )
     {
-        for ( int d = 0; d < particles.dim; d++ )
-        {
-            mesh_min[d] = particles.ghost_mesh_lo[d];
-            mesh_max[d] = particles.ghost_mesh_hi[d];
-        }
     }
 
     template <class ForceType, class PosType, class ParticleType,
               class ParallelType>
     void computeForceFull( ForceType& fc, const PosType& x, const PosType& u,
                            const ParticleType& particles,
+                           const double max_displacement,
                            ParallelType& neigh_op_tag )
     {
         const int n_frozen = particles.frozenOffset();
@@ -60,13 +53,9 @@ class Force<MemorySpace, HertzianModel>
         auto model = _model;
         const auto vol = particles.sliceVolume();
         const auto rho = particles.sliceDensity();
-        const auto y = particles.sliceCurrentPosition();
         const auto vel = particles.sliceVelocity();
 
-        _neigh_timer.start();
-        _neigh_list.build( y, n_frozen, n_local, model.Rc, 1.0, mesh_min,
-                           mesh_max );
-        _neigh_timer.stop();
+        base_type::update( particles, max_displacement );
 
         auto contact_full = KOKKOS_LAMBDA( const int i, const int j )
         {
@@ -111,10 +100,6 @@ class Force<MemorySpace, HertzianModel>
     using base_type::_half_neigh;
     using base_type::_neigh_list;
     using base_type::_timer;
-    Timer _neigh_timer;
-
-    double mesh_max[3];
-    double mesh_min[3];
 };
 
 } // namespace CabanaPD
diff --git a/src/force/CabanaPD_Force_LPS.hpp b/src/force/CabanaPD_Force_LPS.hpp
index 0e9183a8..e0de601a 100644
--- a/src/force/CabanaPD_Force_LPS.hpp
+++ b/src/force/CabanaPD_Force_LPS.hpp
@@ -160,7 +160,7 @@ class Force<MemorySpace, ForceModel<LPS, Elastic, NoFracture>>
     template <class ForceType, class PosType, class ParticleType,
               class ParallelType>
     void computeForceFull( ForceType& f, const PosType& x, const PosType& u,
-                           const ParticleType& particles,
+                           const ParticleType& particles, const double,
                            ParallelType& neigh_op_tag )
     {
         _timer.start();
@@ -379,7 +379,8 @@ class Force<MemorySpace, ForceModel<LPS, Elastic, Fracture>>
     template <class ForceType, class PosType, class ParticleType,
               class ParallelType>
     void computeForceFull( ForceType& f, const PosType& x, const PosType& u,
-                           const ParticleType& particles, ParallelType )
+                           const ParticleType& particles, const double,
+                           ParallelType )
     {
         _timer.start();
 
@@ -533,7 +534,7 @@ class Force<MemorySpace, ForceModel<LinearLPS, Elastic, NoFracture>>
     template <class ForceType, class PosType, class ParticleType,
               class ParallelType>
     void computeForceFull( ForceType& f, const PosType& x, const PosType& u,
-                           const ParticleType& particles,
+                           const ParticleType& particles, const double,
                            ParallelType& neigh_op_tag )
     {
         _timer.start();
diff --git a/src/force/CabanaPD_Force_PMB.hpp b/src/force/CabanaPD_Force_PMB.hpp
index 4687a25d..a79fb515 100644
--- a/src/force/CabanaPD_Force_PMB.hpp
+++ b/src/force/CabanaPD_Force_PMB.hpp
@@ -101,7 +101,7 @@ class Force<MemorySpace, ForceModel<PMB, Elastic, NoFracture, ModelParams...>>
     template <class ForceType, class PosType, class ParticleType,
               class ParallelType>
     void computeForceFull( ForceType& f, const PosType& x, const PosType& u,
-                           const ParticleType& particles,
+                           const ParticleType& particles, const double,
                            ParallelType& neigh_op_tag )
     {
         _timer.start();
@@ -223,7 +223,8 @@ class Force<MemorySpace, ForceModel<PMB, Elastic, Fracture, ModelParams...>>
     template <class ForceType, class PosType, class ParticleType,
               class ParallelType>
     void computeForceFull( ForceType& f, const PosType& x, const PosType& u,
-                           const ParticleType& particles, ParallelType& )
+                           const ParticleType& particles, const double,
+                           ParallelType& )
     {
         _timer.start();
 
@@ -371,7 +372,8 @@ class Force<MemorySpace,
     template <class ForceType, class PosType, class ParticleType,
               class ParallelType>
     void computeForceFull( ForceType& f, const PosType& x, const PosType& u,
-                           ParticleType& particles, ParallelType& neigh_op_tag )
+                           ParticleType& particles, const double,
+                           ParallelType& neigh_op_tag )
     {
         _timer.start();
 
diff --git a/src/force/CabanaPD_HertzianContact.hpp b/src/force/CabanaPD_HertzianContact.hpp
index 17a4442f..11a02a7c 100644
--- a/src/force/CabanaPD_HertzianContact.hpp
+++ b/src/force/CabanaPD_HertzianContact.hpp
@@ -20,31 +20,29 @@
 
 namespace CabanaPD
 {
-struct HertzianModel : public ContactModel
+struct HertzianModel
 {
     // FIXME: This is for use as the primary force model.
     using base_model = PMB;
     using fracture_type = NoFracture;
     using thermal_type = TemperatureIndependent;
 
-    using ContactModel::Rc; // Contact horizon (should be > 2*radius)
-
-    double nu;     // Poisson's ratio
-    double radius; // Actual radius
-    double Rs;     // Equivalent radius
-    double Es;     // Equivalent Young's modulus
-    double e;      // Coefficient of restitution
-    double beta;   // Damping coefficient
-
+    double nu;            // Poisson's ratio
+    double radius;        // Actual radius
+    double radius_extend; // Search radius extended (skin distance)
+    double Rs;            // Equivalent radius
+    double Es;            // Equivalent Young's modulus
+    double e;             // Coefficient of restitution
+    double beta;          // Damping coefficient
     double coeff_h_n;
     double coeff_h_d;
 
-    HertzianModel( const double _Rc, const double _radius, const double _nu,
+    HertzianModel( const double _radius, const double _extend, const double _nu,
                    const double _E, const double _e )
-        : ContactModel( 1.0, _Rc )
+        : nu( _nu )
+        , radius( _radius )
+        , radius_extend( _extend )
     {
-        nu = _nu;
-        radius = _radius;
         Rs = 0.5 * radius;
         Es = _E / ( 2.0 * Kokkos::pow( 1.0 - nu, 2.0 ) );
         e = _e;
diff --git a/unit_test/tstForce.hpp b/unit_test/tstForce.hpp
index 9fef4f1b..2e9d9b71 100644
--- a/unit_test/tstForce.hpp
+++ b/unit_test/tstForce.hpp
@@ -657,7 +657,7 @@ template <class ForceType, class ParticleType>
 double computeEnergyAndForce( ForceType force, ParticleType& particles,
                               const int )
 {
-    computeForce( force, particles, Cabana::SerialOpTag() );
+    computeForce( force, particles, 0.0, Cabana::SerialOpTag() );
     double Phi = computeEnergy( force, particles, Cabana::SerialOpTag() );
     return Phi;
 }
diff --git a/unit_test/tstHertz.hpp b/unit_test/tstHertz.hpp
index 0853295d..e55b872b 100644
--- a/unit_test/tstHertz.hpp
+++ b/unit_test/tstHertz.hpp
@@ -52,8 +52,6 @@ void testHertzianContact( const std::string filename )
     // ====================================================
     double rho0 = inputs["density"];
     double vol = inputs["volume"];
-    double delta = inputs["horizon"];
-    delta += 1e-10;
     double radius = inputs["radius"];
     double nu = inputs["poisson_ratio"];
     double E = inputs["elastic_modulus"];
@@ -88,7 +86,8 @@ void testHertzianContact( const std::string filename )
     //            Force model
     // ====================================================
     using model_type = CabanaPD::HertzianModel;
-    model_type contact_model( delta, radius, nu, E, e );
+    // No search radius extension.
+    model_type contact_model( radius, 0.0, nu, E, e );
 
     // ====================================================
     //                 Particle generation