generate_intersections_from_atom_lists Subroutine

   subroutine generate_intersections_from_atom_lists(atom_lists, n_atoms_list, n_sets, max_k_level, &
                                                   intersections, intersection_sets, intersection_levels, n_intersections)
      !! Generate k-way intersections from arbitrary atom lists (not tied to sys_geom)
      !! max_k_level limits the maximum intersection order to prevent combinatorial explosion
      use pic_logger, only: logger => global_logger
      use pic_io, only: to_char
      integer, intent(in) :: atom_lists(:, :)  !! (max_atoms, n_sets)
      integer, intent(in) :: n_atoms_list(:)   !! Number of atoms in each set
      integer, intent(in) :: n_sets            !! Number of sets (polymers)
      integer, intent(in) :: max_k_level       !! Maximum intersection level to compute
      integer, allocatable, intent(out) :: intersections(:, :)
      integer, allocatable, intent(out) :: intersection_sets(:, :)
      integer, allocatable, intent(out) :: intersection_levels(:)
      integer, intent(out) :: n_intersections

      integer :: max_intersections, max_atoms
      integer, allocatable :: temp_intersections(:, :), temp_sets(:, :), temp_levels(:)
      integer :: intersection_count, k, idx, actual_max_k
      integer, allocatable :: combination(:)

      if (n_sets < 2) then
         n_intersections = 0
         return
      end if

      ! Limit k-way intersections to min(max_k_level, n_sets)
      actual_max_k = min(max_k_level, n_sets)

      max_intersections = 2**n_sets - n_sets - 1
      max_atoms = maxval(n_atoms_list)

      allocate (temp_intersections(max_atoms, max_intersections))
      allocate (temp_sets(n_sets, max_intersections))
      allocate (temp_levels(max_intersections))
      temp_intersections = 0
      temp_sets = 0
      intersection_count = 0

      call logger%info("Generating k-way intersections (k=2 to "//to_char(actual_max_k)//")")

      ! Loop over intersection levels k from 2 to actual_max_k
      do k = 2, actual_max_k
         allocate (combination(k))
         call generate_k_way_intersections_from_lists(atom_lists, n_atoms_list, n_sets, k, &
                                                      combination, max_atoms, &
                                                      temp_intersections, temp_sets, temp_levels, intersection_count)
         deallocate (combination)
      end do

      n_intersections = intersection_count

      ! Allocate output arrays
      if (n_intersections > 0) then
         allocate (intersections(max_atoms, n_intersections))
         allocate (intersection_sets(n_sets, n_intersections))
         allocate (intersection_levels(n_intersections))
         intersections = temp_intersections(1:max_atoms, 1:n_intersections)
         intersection_sets = temp_sets(1:n_sets, 1:n_intersections)
         intersection_levels = temp_levels(1:n_intersections)

         call logger%info("Generated "//to_char(n_intersections)//" total intersections:")
         do k = 2, actual_max_k
            idx = count(intersection_levels == k)
            if (idx > 0) then
               call logger%info("  "//to_char(idx)//" intersections at level "//to_char(k))
            end if
         end do
      else
         allocate (intersections(1, 0))
         allocate (intersection_sets(1, 0))
         allocate (intersection_levels(0))
         call logger%info("No intersections found")
      end if

      deallocate (temp_intersections, temp_sets, temp_levels)
   end subroutine generate_intersections_from_atom_lists