anchor.diff:      integer :: m
       double precision :: x(3)
   
  +    !$omp parallel do private(x, dist_to_C_sq, dist_to_N_sq)
       do m = 1, solvent% Nmax
          if (solvent% flag(m) == 1) then
             x = rel_pos(colloids% pos(:,1), solvent% pos(:,m), solvent_cells% edges)
  
positive.code: 
      far = (L(1)*0.45)**2
  
  
      do n = 1,solvent% Nmax 
         if (solvent% pos(1,n) > bufferlength) then
            if (solvent% species(n) == 2) then
               x = rel_pos(colloids% pos(:,1), solvent% pos(:,n), solvent_cells% edges)
               dist_to_C_sq = dot_product(x, x)
               x= rel_pos(colloids% pos(:,2), solvent% pos(:,n), solvent_cells% edges)
               dist_to_N_sq = dot_product(x, x)
               if ((dist_to_C_sq > far) .and. (dist_to_N_sq > far)) then
                  solvent% species(n) = 1
               end if
            end if
         end if
      end do
  
negative.code: contains
  
    subroutine flag_particles
      double precision :: dist_to_C_sq
      integer :: r, s
      double precision :: x(3)
  
      do s = 1,neigh% n(1)
         r = neigh%list(s,1)
         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 (threefry_double(state(1)) <= prob) then
                  solvent% flag(r) = 1
               end if
            end if
         end if
      end do
  
    end subroutine flag_particles
    
    
  