anchor.diff:          end if
        end do
   
  -     ! Rebuild the neighbor list after possible unbindings
  -     call neigh% make_stencil(solvent_cells, max_cut+skin)
  +     call solvent_cells%random_shift(state(1))
        call apply_pbc(solvent, solvent_cells% edges)
        call apply_pbc(colloids, solvent_cells% edges)
        call solvent% sort(solvent_cells)
  -     call neigh% update_list(colloids, solvent, max_cut+skin, solvent_cells, solvent_colloid_lj)
  +     call neigh% update_list(colloids, solvent, enzyme_capture_radius+skin, solvent_cells, solvent_colloid_lj)
        call varia%tic()
        !$omp parallel do
        do k = 1, solvent%Nmax
  
positive.code:         end do
          call varia%tac()
  
       end do md_loop
  
       current_time = (i+1-equilibration_loops)*tau
  
       call solvent_cells%random_shift(state(1))
       call apply_pbc(solvent, solvent_cells% edges)
       call apply_pbc(colloids, solvent_cells% edges)
       call solvent% sort(solvent_cells)
       call neigh% update_list(colloids, solvent, max_cut+skin, solvent_cells, solvent_colloid_lj)
       call varia%tic()
       !$omp parallel do
       do k = 1, solvent%Nmax
          solvent% pos_old(:,k) = solvent% pos(:,k)
       end do
       colloids% pos_old = colloids% pos
       call varia%tac()
  
       !! TODO: add thermostat and hydro option
       call simple_mpcd_step(solvent, solvent_cells, state)
  
       if (bulk_rmpcd) then
          call bulk_reaction(solvent, solvent_cells, 1, 2, bulk_rate(1), tau, state)
          call bulk_reaction(solvent, solvent_cells, 2, 1, bulk_rate(2), tau, state)
       end if
  
       call neigh% make_stencil(solvent_cells, enzyme_capture_radius)
       call neigh% update_list(colloids, solvent, enzyme_capture_radius, solvent_cells)
       if (sampling) then
          n_plus = 0
          n_minus = 0
          call reset_enzyme_region_bit
          call select_substrate
       end if
  
       ! Check if unbinding should occur
       do enzyme_i = 1, N_enzymes
          if (current_time > next_reaction_time(enzyme_i)) then
             if (threefry_double(state(1)) < rate_release_s / total_rate) then
                ! release substrate
                call unbind_molecule(enzyme_i, 1)
  
negative.code:   skin = 1
    n_extra_sorting = 0
  
    call neigh% make_stencil(solvent_cells, enzyme_capture_radius+skin)
  
    call neigh% update_list(colloids, solvent, enzyme_capture_radius+skin, solvent_cells, solvent_colloid_lj)
  
    e1 = compute_force(colloids, solvent, neigh, solvent_cells% edges, solvent_colloid_lj)
    e2 = compute_force_n2(colloids, solvent_cells% edges, colloid_lj)
  