recursive subroutine generate_k_way_intersections_for_level(sys_geom, monomers, n_monomers, k, &
combination, max_atoms, &
temp_intersections, temp_sets, temp_levels, intersection_count)
!! Helper to generate all k-way intersections at a specific level k
use mqc_physical_fragment, only: system_geometry_t
type(system_geometry_t), intent(in) :: sys_geom
integer, intent(in) :: monomers(:), n_monomers, k, max_atoms
integer, intent(inout) :: combination(:)
integer, intent(inout) :: temp_intersections(:, :), temp_sets(:, :), temp_levels(:)
integer, intent(inout) :: intersection_count
integer, allocatable :: current_intersection(:), temp_intersection(:)
integer :: current_n_intersect, temp_n_intersect
integer :: i, mono_idx, frag_size
logical :: has_intersection
! Generate combinations using an iterative approach
call next_combination_init(combination, k)
do
! Compute intersection of all fragments in this combination
has_intersection = .false.
! Start with the first fragment in the combination
mono_idx = monomers(combination(1))
frag_size = sys_geom%fragment_sizes(mono_idx)
allocate (current_intersection(frag_size))
current_intersection = sys_geom%fragment_atoms(1:frag_size, mono_idx)
current_n_intersect = frag_size
! Intersect with each subsequent fragment
do i = 2, k
mono_idx = monomers(combination(i))
frag_size = sys_geom%fragment_sizes(mono_idx)
has_intersection = find_fragment_intersection( &
current_intersection, current_n_intersect, &
sys_geom%fragment_atoms(1:frag_size, mono_idx), frag_size, &
temp_intersection, temp_n_intersect)
if (.not. has_intersection) then
! Intersection is empty, break early
deallocate (current_intersection)
if (allocated(temp_intersection)) deallocate (temp_intersection)
exit
end if
! Replace current_intersection with the new intersection
deallocate (current_intersection)
allocate (current_intersection(temp_n_intersect))
current_intersection = temp_intersection
current_n_intersect = temp_n_intersect
deallocate (temp_intersection)
end do
! If we have a non-empty intersection, store it
if (has_intersection .and. allocated(current_intersection)) then
intersection_count = intersection_count + 1
temp_intersections(1:current_n_intersect, intersection_count) = current_intersection
temp_sets(1:k, intersection_count) = monomers(combination)
temp_levels(intersection_count) = k
deallocate (current_intersection)
end if
! Get next combination
if (.not. next_combination(combination, k, n_monomers)) exit
end do
end subroutine generate_k_way_intersections_for_level