anchor.diff:      double precision :: dist_to_N_sq
       integer :: m
       double precision :: x(3)
  -
  -    !$omp parallel do private(x, dist_to_C_sq, dist_to_N_sq)
  +    integer :: thread_id
  +
  +    !$omp parallel private(thread_id)
  +    thread_id = omp_get_thread_num() + 1
  +    !$omp 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:                (dist_to_N_sq > solvent_colloid_lj%cut_sq(1,2)) &
                 ) &
                 then
               solvent% species(m) = 2
               solvent% flag(m) = 0
            end if
         end if
      end do
  
    end subroutine change_species
    
    subroutine refuel
  
negative.code:           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
               solvent% flag(r) = 1
            end if
         end if
      end do
  