anchor.diff:    call init_genrand(mt, int(clock, c_long))
   
     call h5open_f(error)
  +
  +  prob = PTread_d(config,'probability')
   
     L = PTread_ivec(config, 'L', 3)
     rho = PTread_i(config, 'rho')
  
positive.code:     double precision :: x(3)
      integer :: n
  
      far = 25.d0
  
      do n = 1,solvent% Nmax
         if (solvent% species(n) == 2) then
  
negative.code: 
    subroutine flag_particles
    double precision :: dist_to_C_sq
    real :: rndnumbers(solvent% Nmax)
    integer :: r
    double precision :: x(3)
    
    call random_number(rndnumbers)
    
    do  r = 1,solvent% Nmax
       if (solvent% species(r) == 1) 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
  