anchor.diff:          colloids% pos = colloids% pos + dt * colloids% vel + dt**2 * colloids% force / (2 * mass)
   
           call rattle_dimer_pos(colloids, d, dt)
  +
  +        so_max = solvent% maximum_displacement(solvent_cells% edges)
  +        co_max = colloids% maximum_displacement(solvent_cells% edges)
  +
  +        if ( (co_max >= skin/2) .or. (so_max >= skin/2) ) then
  +           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
   
           do k = 1, colloids% Nmax
              jump = floor(colloids% pos(:,k) / solvent_cells% edges)
  
positive.code:   call solvent% sort(solvent_cells)
  
    call neigh% init(colloids% Nmax, int(300*sigma**3))
    call neigh% make_stencil(solvent_cells, 1.5d0*sigma)
  
    call neigh% update_list(colloids, solvent, 1.5d0*sigma, solvent_cells)
  
    tau = 0.1d0 !was 0.1d0 but we changed it to Kaprals' value 
    N_MD_steps = 100
    dt = tau / N_MD_steps
  
  
negative.code:   double precision :: e1, e2
    double precision :: tau, dt
    double precision :: d
    double precision :: skin, co_max, so_max
    integer :: N_MD_steps
    integer :: n_extra_sorting
  
    type(mt19937ar_t), target :: mt
  
  