anchor.diff:  
     call solvent% sort(solvent_cells)
   
  -  call neigh% init(colloids% Nmax, int(300*sigma**3))
  +  call neigh% init(colloids% Nmax, int(300*max(sigma_C,sigma_N)**3))
   
     skin = 1.5
     n_extra_sorting = 0
   
  -  call neigh% make_stencil(solvent_cells, sigma_cut+skin)
  -
  -  call neigh% update_list(colloids, solvent, sigma_cut+skin, solvent_cells)
  -
  -  tau = 1.d0
  -  N_MD_steps = 100
  -  dt = tau / N_MD_steps
  +  call neigh% make_stencil(solvent_cells, max_cut+skin)
  +
  +  call neigh% update_list(colloids, solvent, max_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:      if (solvent% flag(m) == 1) then
         dist_to_C_sq = dot_product(colloids% pos(:,1) - solvent% pos(:,m),colloids% pos(:,1) - solvent% pos(:,m))
         dist_to_N_sq = dot_product(colloids% pos(:,2) - solvent% pos(:,m),colloids% pos(:,2) - solvent% pos(:,m))
         if ((dist_to_C_sq > (sigma_cut)**2) .and. (dist_to_N_sq > (sigma_cut)**2)) then
           solvent% species(m) = 2
           solvent% flag(m) = 0
         end if
  
negative.code:   write(*,*) '    i           |    e co so     |   e co co     |   kin co      |   kin so      |   total       |   temp        |'
    write(*,*) ''
  
    do i = 1, N_loop
       md: do j = 1, N_MD_steps
          call md_pos(solvent, dt)
  
          ! Extra copy for rattle
          colloids% pos_rattle = colloids% pos
          do k=1, colloids% Nmax
             colloids% pos(:,k) = colloids% pos(:,k) + dt * colloids% vel(:,k) + &
                  dt**2 * colloids% force(:,k) / (2 * mass(colloids%species(k)))
          end do
  
          call rattle_dimer_pos(colloids, d, dt)
  
  