anchor.diff:    colloids% vel = 0
     
     call random_number(solvent% vel(:, :))
  -  solvent% vel = (solvent% vel - 0.5d0)*sqrt(4.d0*3.d0*T)
  +  solvent% vel = (solvent% vel - 0.5d0)*sqrt(12*T)
     solvent% vel = solvent% vel - spread(sum(solvent% vel, dim=2)/solvent% Nmax, 2, solvent% Nmax)
     solvent% force = 0
     solvent% species = 1
  
positive.code:   call h5close_f(error)
    
  contains
  
    subroutine flag_particles
    double precision :: dist_to_C_sq
  
negative.code:      call simple_mpcd_step(solvent, solvent_cells, mt)
  
       kin_co = (mass(1)*sum(colloids% vel(:,1)**2)+mass(2)*sum(colloids% vel(:,2)**2))/2
       call thermo_write
  
    end do
  
  