anchor.diff:  
     solvent% pos_old = solvent% pos
     colloids% pos_old = colloids% pos
  -  do i = 1, 100
  -     so_max = 0
  -     co_max = 0
  +  do i = 1, 2000
        md_loop: do j = 1, N_MD_steps
  -        !$omp parallel do private(k, tmp_x)
  -        do k = 1, solvent% Nmax
  -           solvent% pos(:,k) = solvent% pos(:,k) + dt * solvent% vel(:,k) + dt**2 * solvent% force(:,k) / 2
  -        end do
  +        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/2) .or. (so_max >= skin/2) ) then
  +        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)
  +           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
  -           write(*,*) 'extra sorting'
           end if
   
           call switch(solvent% force, solvent% force_old)
  
positive.code:   skin = 0.8
    n_extra_sorting = 0
  
    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)
  
negative.code:         e1 = compute_force(colloids, solvent, neigh, solvent_cells% edges, solvent_colloid_lj)
          e2 = compute_force_n2(colloids, solvent_cells% edges, colloid_lj)
  
          call md_vel(solvent, dt)
          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 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
  