anchor.diff:  
     subroutine flag_particles
     double precision :: dist_to_C_sq
  -  real :: prob 
     real :: rndnumbers(solvent% Nmax)
     integer :: r
  -  prob = 1.0
  -  
  -  call random_number(prob)
  +  double precision :: x(3)
  +  
  +  call random_number(rndnumbers)
     
     do  r = 1,solvent% Nmax
        if (solvent% species(r) == 1) then
  -       dist_to_C_sq = dot_product(colloids% pos(:,1) - solvent% pos(:,r),colloids% pos(:,1) - solvent% pos(:,r))
  -       if (dist_to_C_sq < 3**2) then
  +       x = rel_pos(colloids% pos(:,1),solvent% pos(:,r),solvent_cells% edges) 
  +       dist_to_C_sq = dot_product(x, x)
  +       if (dist_to_C_sq < solvent_colloid_lj%cut_sq(1,1)) then
            if (rndnumbers(r) <= prob) then
              solvent% flag(r) = 1 
            end if
  
positive.code:   call init_genrand(mt, int(clock, c_long))
  
    call h5open_f(error)
  
    L = PTread_ivec(config, 'L', 3)
    rho = PTread_i(config, 'rho')
  
negative.code:     double precision :: x(3)
      integer :: n
  
      far = (L(1)/2)**2
  
  
      do n = 1,solvent% Nmax
         if (solvent% species(n) == 2) then
  