generate_intersections Subroutine

   subroutine generate_intersections(sys_geom, monomers, polymers, n_monomers, max_intersection_level, &
                                     intersections, intersection_sets, intersection_levels, n_intersections)
      !! Generate all k-way intersections for k=2 to min(max_intersection_level, n_monomers)
      !!
      !! For a system with overlapping fragments, this computes k-way intersections
      !! following the inclusion-exclusion principle for GMBE.
      !! The max_intersection_level parameter controls the maximum depth to avoid combinatorial explosion.
      !!
      !! Algorithm:
      !! - For each k from 2 to min(max_intersection_level, n_monomers):
      !!   - Generate all C(n_monomers, k) combinations
      !!   - For each combination, compute intersection of all k fragments
      !!   - Store non-empty intersections with their level k
      use mqc_physical_fragment, only: system_geometry_t
      type(system_geometry_t), intent(in) :: sys_geom
      integer, intent(in) :: monomers(:)           !! Monomer indices
      integer, intent(inout) :: polymers(:, :)     !! Output: monomers stored here
      integer, intent(in) :: n_monomers            !! Number of monomers
      integer, intent(in) :: max_intersection_level  !! Maximum k-way intersection depth
      integer, allocatable, intent(out) :: intersections(:, :)  !! Intersection atom lists
      integer, allocatable, intent(out) :: intersection_sets(:, :)  !! Which k-tuple created each intersection
      integer, allocatable, intent(out) :: intersection_levels(:)  !! Level (k) of each intersection
      integer, intent(out) :: n_intersections      !! Number of intersections found

      ! Temporaries for storing intersections
      integer, allocatable :: temp_intersections(:, :)
      integer, allocatable :: temp_sets(:, :)
      integer, allocatable :: temp_levels(:)
      integer, allocatable :: temp_intersection(:)
      integer, allocatable :: current_intersection(:)
      integer :: temp_n_intersect, current_n_intersect
      logical :: has_intersection
      integer :: k, intersection_count, max_atoms, max_intersections, max_k_level
      integer :: i, idx
      integer, allocatable :: combination(:)

      ! Store monomers in polymers array
      polymers(1:n_monomers, 1) = monomers(1:n_monomers)

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

      ! Count maximum possible intersections: sum of C(n,k) for k=2 to n
      ! For small n, this is 2^n - n - 1
      max_intersections = 2**n_monomers - n_monomers - 1

      ! Find maximum atoms in any fragment for allocation
      max_atoms = maxval(sys_geom%fragment_sizes(1:n_monomers))

      ! Allocate temporary arrays
      allocate (temp_intersections(max_atoms, max_intersections))
      allocate (temp_sets(n_monomers, max_intersections))
      allocate (temp_levels(max_intersections))
      temp_intersections = 0
      temp_sets = 0
      intersection_count = 0

      ! Determine actual maximum intersection level to use
      max_k_level = min(max_intersection_level, n_monomers)

      if (max_k_level < n_monomers) then
         call logger%info("Generating k-way intersections up to k="//to_char(max_k_level)// &
                          " (limited by max_intersection_level)")
      else
         call logger%info("Generating all k-way intersections for GMBE (inclusion-exclusion principle)")
      end if

      ! Loop over intersection levels k from 2 to max_k_level
      do k = 2, max_k_level
         ! Generate all C(n_monomers, k) combinations
         allocate (combination(k))
         call generate_k_way_intersections_for_level(sys_geom, monomers, n_monomers, 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_monomers, n_intersections))
         allocate (intersection_levels(n_intersections))
         intersections = temp_intersections(1:max_atoms, 1:n_intersections)
         intersection_sets = temp_sets(1:n_monomers, 1:n_intersections)
         intersection_levels = temp_levels(1:n_intersections)

         call logger%info("Generated "//to_char(n_intersections)//" total intersections:")
         do k = 2, max_k_level
            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
         call logger%info("No intersections found (fragments are non-overlapping)")
      end if

      deallocate (temp_intersections, temp_sets, temp_levels)

   end subroutine generate_intersections