anchor.diff:  
     call neigh% update_list(colloids, solvent, 1.5d0*sigma, solvent_cells)
   
  -  tau = 0.1d0 !was 0.1d0 but we changed it to Kaprals' value 
  +  tau = 0.1d0 !was 0.1d0 but we changed it to Kaprals' value(1.0) => doesn't work
     N_MD_steps = 100
     dt = tau / N_MD_steps
   
  
positive.code:   call colloids% init(2,2) !there will be 2 species of colloids
    allocate(pos_rattle(3, colloids% Nmax))
    
    open(15,file ='dimerdata_chem01.txt')
    
    write(*, *) colloids% pos
    colloids% species(1) = 1
  
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
  