anchor.diff:       end do md
   
        call solvent% sort(solvent_cells)
  -     call neigh% update_list(colloids, solvent, 1.5d0, 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)
   
  
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:   double precision :: e1, e2
    double precision :: tau, dt
    integer :: N_MD_steps
    integer :: n_extra_sorting
  
    type(mt19937ar_t), target :: mt
  
    double precision :: tmp_x(3)
    double precision :: skin
  
    integer :: i, L(3), seed_size, clock
    integer :: jump(3)
  