anchor.diff:  
     end do
   
  +  write(*,*) 'n extra sorting', n_extra_sorting
   
     call h5close_f(error)
   
  
positive.code:   write(*,*) '    e co so     |   e co co     |   kin co      |   kin so      |   total       |   temp        |'
    write(*,*) ''
  
    do i = 1, 100
       so_max = 0
       co_max = 0
       md: do j = 1, N_MD_steps
          solvent% pos_old = solvent% pos + dt * solvent% vel + dt**2 * solvent% force / 2
          !$omp parallel do private(k)
          do k = 1, solvent% Nmax
             solvent% pos(:,k) = modulo(solvent% pos_old(:,k), solvent_cells% edges)
          end do
          colloids% pos_old = colloids% pos + dt * colloids% vel + dt**2 * colloids% force / (2 * mass)
          do k = 1, colloids% Nmax
             jump = floor(colloids% pos_old(:,k) / solvent_cells% edges)
             colloids% image(:,k) = colloids% image(:,k) + jump
             colloids% pos(:,k) = colloids% pos_old(:,k) - jump*solvent_cells% edges
          end do
          so_max = max(maxval(sqrt(sum((solvent% pos - solvent% pos_old)**2, dim=1))), so_max)
          co_max = max(maxval(sqrt(sum((colloids% pos - colloids% pos_old)**2, dim=1))), co_max)
  
          call switch(solvent% force, solvent% force_old)
          call switch(colloids% force, colloids% force_old)
  
negative.code:      end do md
  
       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
  
       call simple_mpcd_step(solvent, solvent_cells, mt)
  
  