anchor.diff:       call apply_pbc(solvent, solvent_cells% edges)
        call apply_pbc(colloids, solvent_cells% edges)
        call solvent% sort(solvent_cells)
  +     loop_i_last_sort = N_MD_steps*(i+1)
        call neigh% update_list(colloids, solvent, sigma_cut + skin, solvent_cells)
        solvent% pos_old = solvent% pos
        colloids% pos_old = colloids% pos
  
positive.code:   call neigh% update_list(colloids, solvent, sigma_cut+skin, solvent_cells)
  
    n_extra_sorting = 0
  
    e1 = compute_force(colloids, solvent, neigh, solvent_cells% edges, solvent_colloid_lj)
    e2 = compute_force_n2(colloids, solvent_cells% edges, colloid_lj)
  
negative.code:   double precision :: tau, dt, T, T_init, T_final
    integer :: N_MD_steps, N_loop
    integer :: colloid_sampling
    integer :: n_extra_sorting, loop_i_last_sort
    integer :: n_threads
    logical :: do_hydro, do_thermostat
  
  