anchor.diff:    skin = 1
     n_extra_sorting = 0
   
  -  call neigh% make_stencil(solvent_cells, max_cut+skin)
  -
  -  call neigh% update_list(colloids, solvent, max_cut+skin, solvent_cells, solvent_colloid_lj)
  +  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)
  
positive.code:   do i = 0, N_loop+equilibration_loops
       if (i==equilibration_loops) sampling = .true.
       if (modulo(i,64) == 0) write(*,'(i05)',advance='no') i
       md_loop: do j = 1, N_MD_steps
          call md_pos(solvent, dt)
          do k=1, colloids% Nmax
             colloids% pos(:,k) = colloids% pos(:,k) + dt * colloids% vel(:,k) + &
  
negative.code:      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 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()
  
       call varia%tic()
  
  