anchor.diff:    write(*,*) '    i   |    e co so     |   e co co     |   kin co      |   kin so      |   total       |   temp        |'
     write(*,*) ''
   
  -  do i = 1, 500
  +  do i = 1, 1000
        so_max = 0
        co_max = 0
        md: do j = 1, N_MD_steps
  
positive.code: 
    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:   colloids% vel = 0
    
    call random_number(solvent% vel(:, :))
    solvent% vel = (solvent% vel - 0.5d0)*sqrt(6.d0*2.d0/3.d0)
    solvent% vel = solvent% vel - spread(sum(solvent% vel, dim=2)/solvent% Nmax, 2, solvent% Nmax)
    solvent% force = 0
    solvent% species = 1
  