gmbe_enumerate_pie_terms Subroutine

   subroutine gmbe_enumerate_pie_terms(sys_geom, primaries, n_primaries, polymer_level, max_k_level, &
                                       pie_atom_sets, pie_coefficients, n_pie_terms, error, initial_max_terms)
      !! Enumerate all unique intersections via DFS and accumulate PIE coefficients
      !! This implements the GMBE(N) algorithm with inclusion-exclusion principle
      !!
      !! Algorithm:
      !! 1. For each primary i, start DFS with clique=[i]
      !! 2. Recursively grow cliques by adding overlapping primaries
      !! 3. For each clique of size k, compute intersection and add PIE coefficient:
      !!    coefficient = (+1) if k odd, (-1) if k even
      !! 4. Accumulate coefficients for each unique atom set
      use mqc_physical_fragment, only: system_geometry_t

      type(system_geometry_t), intent(in) :: sys_geom
      integer, intent(in) :: primaries(:, :)   !! Primary polymers (n_primaries, polymer_level)
      integer, intent(in) :: n_primaries       !! Number of primary polymers
      integer, intent(in) :: polymer_level     !! Level of primaries (1=monomers, 2=dimers, etc.)
      integer, intent(in) :: max_k_level       !! Maximum clique size (intersection depth limit)
      integer, allocatable, intent(out) :: pie_atom_sets(:, :)  !! Unique atom sets (max_atoms, n_terms)
      integer, allocatable, intent(out) :: pie_coefficients(:)  !! PIE coefficient for each term
      integer(int64), intent(out) :: n_pie_terms      !! Number of unique PIE terms
      type(error_t), intent(out) :: error             !! Error status
      integer(int64), intent(in), optional :: initial_max_terms  !! Initial PIE storage capacity

      ! Temporary storage for PIE terms (allocate generously)
      integer(int64), parameter :: INITIAL_MAX_PIE_TERMS = 100000_int64  ! Adjust if needed
      integer(int64) :: max_terms
      integer(default_int) :: max_terms_i
      integer, allocatable :: temp_atom_sets(:, :)
      integer, allocatable :: temp_coefficients(:)
      integer, allocatable :: primary_atoms(:, :)    !! Precomputed atom lists for each primary
      integer, allocatable :: primary_n_atoms(:)     !! Atom counts for each primary
      integer, allocatable :: clique(:)              !! Current clique being built
      integer, allocatable :: current_atoms(:)       !! Current intersection atoms
      integer, allocatable :: candidates(:)          !! Candidate primaries to add
      integer :: i, j, max_atoms, n_candidates

      call logger%info("Enumerating GMBE PIE terms via DFS...")

      ! Find maximum atoms in any fragment
      max_atoms = sys_geom%total_atoms

      ! Allocate temporary storage
      if (present(initial_max_terms)) then
         max_terms = initial_max_terms
      else
         max_terms = INITIAL_MAX_PIE_TERMS
      end if
      if (max_terms < 1_int64) then
         call error%set(ERROR_VALIDATION, "Initial PIE term capacity must be positive")
         return
      end if
      if (max_terms > int(huge(0_default_int), int64)) then
         call error%set(ERROR_VALIDATION, "Initial PIE term capacity exceeds default integer limit")
         return
      end if
      max_terms_i = int(max_terms, default_int)

      allocate (temp_atom_sets(max_atoms, max_terms_i))
      allocate (temp_coefficients(max_terms_i))
      temp_atom_sets = -1
      temp_coefficients = 0
      n_pie_terms = 0_int64

      ! Precompute atom lists for all primaries
      allocate (primary_atoms(max_atoms, n_primaries))
      allocate (primary_n_atoms(n_primaries))
      primary_atoms = -1

      do i = 1, n_primaries
         block
            integer, allocatable :: atom_list(:)
            integer :: n_atoms
            call compute_polymer_atoms(sys_geom, primaries(i, :), polymer_level, atom_list, n_atoms)
            primary_n_atoms(i) = n_atoms
            primary_atoms(1:n_atoms, i) = atom_list(1:n_atoms)
            deallocate (atom_list)
         end block
      end do

      ! Allocate work arrays
      allocate (clique(max_k_level))
      allocate (current_atoms(max_atoms))
      allocate (candidates(n_primaries))

      ! Start DFS from each primary
      do i = 1, n_primaries
         ! Initial clique: just primary i
         clique(1) = i
         current_atoms(1:primary_n_atoms(i)) = primary_atoms(1:primary_n_atoms(i), i)

         ! Candidates: all primaries after i (to avoid duplicates)
         n_candidates = n_primaries - i
         if (n_candidates > 0) then
            candidates(1:n_candidates) = [(i + j, j=1, n_candidates)]
         end if

         ! DFS from this primary
         call dfs_pie_accumulate(primary_atoms, primary_n_atoms, n_primaries, max_atoms, &
                                 clique, 1, current_atoms(1:primary_n_atoms(i)), primary_n_atoms(i), &
                                 candidates, n_candidates, max_k_level, &
                                 temp_atom_sets, temp_coefficients, n_pie_terms, max_terms, error)
         if (error%has_error()) return
      end do

      ! Copy to output arrays
      if (n_pie_terms > 0_int64) then
         if (n_pie_terms > int(huge(0_default_int), int64)) then
            call error%set(ERROR_VALIDATION, "n_pie_terms exceeds default integer limit")
            return
         end if
         allocate (pie_atom_sets(max_atoms, int(n_pie_terms, default_int)))
         allocate (pie_coefficients(int(n_pie_terms, default_int)))
         pie_atom_sets = temp_atom_sets(:, 1:n_pie_terms)
         pie_coefficients = temp_coefficients(1:n_pie_terms)
      end if

      call logger%info("Generated "//to_char(n_pie_terms)//" unique PIE terms")

      ! Cleanup
      deallocate (temp_atom_sets, temp_coefficients, primary_atoms, primary_n_atoms)
      deallocate (clique, current_atoms, candidates)

   end subroutine gmbe_enumerate_pie_terms