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: 
       call simple_mpcd_step(solvent, solvent_cells, state)
  
       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)
  
negative.code: 
    solvent% pos_old = solvent% pos
    colloids% pos_old = colloids% pos
    do i = 1, 2000
       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) + dt**2 * colloids% force(:,k) / (2 * mass)
          end do
          so_max = solvent% maximum_displacement()
          co_max = colloids% maximum_displacement()
  
          if ( (co_max >= skin*0.1) .or. (so_max >= skin*0.9) ) then
             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, sigma_cut+skin, solvent_cells)
             solvent% pos_old = solvent% pos
             colloids% pos_old = colloids% pos
             n_extra_sorting = n_extra_sorting + 1
          end if
  
          call switch(solvent% force, solvent% force_old)
  