anchor.diff:  
     call solvent% init(N,2) !there will be 2 species of solvent particles
   
  -  call colloids% init(2,2, mass=[mass, mass]) !there will be 2 species of colloids
  +  call colloids% init(2,2, mass) !there will be 2 species of colloids
   
     call PTkill(config)
     
  
positive.code: 
  
       write(15,*) colloids% pos + colloids% image * spread(solvent_cells% edges, dim=2, ncopies=colloids% Nmax), &
                   colloids% vel, e1+e2+mass*sum(colloids% vel**2)/2+sum(solvent% vel**2)/2
       
       solvent_cells% origin(1) = genrand_real1(mt) - 1
       solvent_cells% origin(2) = genrand_real1(mt) - 1
  
negative.code: 
       call simple_mpcd_step(solvent, solvent_cells, mt)
  
       kin_co = (colloids% mass(1)*sum(colloids% vel(:,1)**2)+ colloids% mass(2)*sum(colloids% vel(:,2)**2))/2
       write(*,'(1i16,6f16.3,1e16.8)') i, e1, e2, &
            kin_co, sum(solvent% vel**2)/2, &
            e1+e2+kin_co+sum(solvent% vel**2)/2, &
  