anchor.diff:    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
   
  
positive.code:         colloids% pos = colloids% pos + dt * colloids% vel + dt**2 * colloids% force / (2 * mass)
  
          call rattle_dimer_pos(colloids, d, dt)
  
          do k = 1, colloids% Nmax
             jump = floor(colloids% pos(:,k) / solvent_cells% edges)
  
negative.code:   call solvent% sort(solvent_cells)
  
    call neigh% init(colloids% Nmax, int(300*sigma**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
  
  