anchor.diff:    skin = 0.8
     n_extra_sorting = 0
   
  -  call neigh% update_list(colloids, solvent, sigma_cut + skin, solvent_cells)
  +  call neigh% make_stencil(solvent_cells, sigma_cut+skin)
  +  call neigh% update_list(colloids, solvent, sigma_cut+skin, solvent_cells)
   
     e1 = compute_force(colloids, solvent, neigh, solvent_cells% edges, solvent_colloid_lj)
     e2 = compute_force_n2(colloids, solvent_cells% edges, colloid_lj)
  
positive.code:         e1 = compute_force(colloids, solvent, neigh, solvent_cells% edges, solvent_colloid_lj)
          e2 = compute_force_n2(colloids, solvent_cells% edges, colloid_lj)
  
          !$omp parallel do private(k)
          do k = 1, solvent% Nmax
             solvent% vel(:,k) = solvent% vel(:,k) + dt * ( solvent% force(:,k) + solvent% force_old(:,k) ) / 2
          end do
          colloids% vel = colloids% vel + dt * ( colloids% force + colloids% force_old ) / (2 * mass)
  
          write(15,*) colloids% pos + colloids% image * spread(solvent_cells% edges, dim=2, ncopies=colloids% Nmax)
  
       end do md_loop
  
       call solvent% sort(solvent_cells)
       call neigh% update_list(colloids, solvent, sigma_cut + skin, solvent_cells)
       solvent% pos_old = solvent% pos
  
negative.code: 
       call simple_mpcd_step(solvent, solvent_cells, state)
  
       if (modulo(i,10)==0) &
       write(*,'(6f16.3)') e1, e2, mass*sum(colloids% vel**2)/2, sum(solvent% vel**2)/2, &
           e1+e2+mass*sum(colloids% vel**2)/2+sum(solvent% vel**2)/2, &
           compute_temperature(solvent, solvent_cells)
  