compute_gmbe Subroutine

   subroutine compute_gmbe(monomers, n_monomers, monomer_results, &
                           n_intersections, intersection_results, &
                           intersection_sets, intersection_levels, &
                           total_energy, &
                           sys_geom, total_gradient, total_hessian, bonds, world_comm)
      !! Compute generalized many-body expansion (GMBE) energy with optional gradient and/or hessian
      !!
      !! This is the core routine that handles all GMBE computations using
      !! the inclusion-exclusion principle for overlapping fragments.
      !! Optional arguments control what derivatives are computed:
      !!   - If sys_geom and total_gradient are present: compute gradient
      !!   - If total_hessian is also present: compute hessian
      use mqc_result_types, only: calculation_result_t
      use mqc_physical_fragment, only: build_fragment_from_indices, build_fragment_from_atom_list, &
                                       redistribute_cap_gradients, redistribute_cap_hessian
      use mqc_gmbe_utils, only: find_fragment_intersection
      use mqc_config_parser, only: bond_t
      use mqc_error, only: error_t

      ! Required arguments
      integer, intent(in) :: monomers(:)
      integer, intent(in) :: n_monomers
      type(calculation_result_t), intent(in) :: monomer_results(:)
      integer, intent(in) :: n_intersections
      type(calculation_result_t), intent(in), optional :: intersection_results(:)
      integer, intent(in), optional :: intersection_sets(:, :)
      integer, intent(in), optional :: intersection_levels(:)
      real(dp), intent(out) :: total_energy

      ! Optional arguments for gradient computation
      type(system_geometry_t), intent(in), optional :: sys_geom
      real(dp), intent(out), optional :: total_gradient(:, :)

      ! Optional arguments for hessian computation
      real(dp), intent(out), optional :: total_hessian(:, :)

      ! Optional bond information for hydrogen cap handling
      type(bond_t), intent(in), optional :: bonds(:)

      ! Optional MPI communicator for abort
      type(comm_t), intent(in), optional :: world_comm

      ! Local variables
      integer :: i, k, max_level, current_log_level, hess_dim
      real(dp) :: monomer_energy, sign_factor
      real(dp), allocatable :: level_energies(:)
      integer, allocatable :: level_counts(:)
      type(physical_fragment_t) :: fragment
      type(error_t) :: error
      integer, allocatable :: monomer_idx(:)
      logical :: compute_grad, compute_hess, has_intersections

      ! Determine what to compute based on optional arguments
      compute_grad = present(sys_geom) .and. present(total_gradient)
      compute_hess = compute_grad .and. present(total_hessian)
      has_intersections = n_intersections > 0 .and. present(intersection_results) .and. &
                          present(intersection_sets) .and. present(intersection_levels)

      ! Initialize outputs
      if (compute_grad) then
         total_gradient = 0.0_dp
      end if
      if (compute_hess) then
         total_hessian = 0.0_dp
         hess_dim = 3*sys_geom%total_atoms
      end if

      ! Sum monomer energies (and gradients/hessians if requested)
      monomer_energy = 0.0_dp
      do i = 1, n_monomers
         monomer_energy = monomer_energy + monomer_results(i)%energy%total()

         ! Accumulate monomer derivatives if requested
         if (compute_grad .and. monomer_results(i)%has_gradient) then
            allocate (monomer_idx(1))
            monomer_idx(1) = monomers(i)
            call build_fragment_from_indices(sys_geom, monomer_idx, fragment, error, bonds)
            if (error%has_error()) then
               call logger%error(error%get_full_trace())
               if (present(world_comm)) then
                  call abort_comm(world_comm, 1)
               else
                  error stop "Failed to build monomer fragment in GMBE"
               end if
            end if
            call redistribute_cap_gradients(fragment, monomer_results(i)%gradient, total_gradient)
            if (compute_hess .and. monomer_results(i)%has_hessian) then
               call redistribute_cap_hessian(fragment, monomer_results(i)%hessian, total_hessian)
            end if
            call fragment%destroy()
            deallocate (monomer_idx)
         end if
      end do

      ! Start with monomer contribution
      total_energy = monomer_energy

      ! Handle intersections with inclusion-exclusion
      if (has_intersections) then
         max_level = maxval(intersection_levels)

         allocate (level_energies(max_level))
         allocate (level_counts(max_level))
         level_energies = 0.0_dp
         level_counts = 0

         ! Process intersection energies and derivatives
         do i = 1, n_intersections
            k = intersection_levels(i)
            sign_factor = real((-1)**(k + 1), dp)
            level_energies(k) = level_energies(k) + intersection_results(i)%energy%total()
            level_counts(k) = level_counts(k) + 1

            ! Handle intersection derivatives if requested
            if (compute_grad .and. (intersection_results(i)%has_gradient .or. &
                                    (compute_hess .and. intersection_results(i)%has_hessian))) then
               call process_intersection_derivatives(i, k, sign_factor, intersection_results, &
                                                     intersection_sets, sys_geom, total_gradient, total_hessian, bonds, &
                                                     compute_grad, compute_hess, hess_dim, world_comm)
            end if
         end do

         ! Apply inclusion-exclusion to energy
         do k = 2, max_level
            if (level_counts(k) > 0) then
               sign_factor = real((-1)**(k + 1), dp)
               total_energy = total_energy + sign_factor*level_energies(k)
            end if
         end do

         ! Print energy breakdown
         call print_gmbe_energy_breakdown(monomer_energy, n_monomers, level_energies, level_counts, &
                                          max_level, total_energy, .true.)

         ! Print debug info for intersections
         call print_gmbe_intersection_debug(n_intersections, n_monomers, intersection_sets, &
                                            intersection_levels, intersection_results)

         deallocate (level_energies, level_counts)
      else
         ! No intersections - just print monomer sum
         call print_gmbe_energy_breakdown(monomer_energy, n_monomers, level_energies, level_counts, &
                                          0, total_energy, .false.)
      end if

      ! Print gradient/hessian info
      if (compute_grad) then
         call logger%configuration(level=current_log_level)
         call print_gmbe_gradient_info(total_gradient, sys_geom, current_log_level)
      end if

      if (compute_hess) then
         call logger%info("GMBE Hessian computation completed")
         call logger%info("  Total Hessian Frobenius norm: "//to_char(sqrt(sum(total_hessian**2))))

         ! Compute and print full vibrational analysis
         block
            real(dp), allocatable :: frequencies(:), reduced_masses(:), force_constants(:)
            real(dp), allocatable :: cart_disp(:, :), fc_mdyne(:)

            call logger%info("  Computing vibrational analysis (projecting trans/rot modes)...")
            call compute_vibrational_analysis(total_hessian, sys_geom%element_numbers, frequencies, &
                                              reduced_masses, force_constants, cart_disp, &
                                              coordinates=sys_geom%coordinates, &
                                              project_trans_rot=.true., &
                                              force_constants_mdyne=fc_mdyne)

            if (allocated(frequencies)) then
               call print_vibrational_analysis(frequencies, reduced_masses, force_constants, &
                                               cart_disp, sys_geom%element_numbers, &
                                               force_constants_mdyne=fc_mdyne, &
                                               coordinates=sys_geom%coordinates, &
                                               electronic_energy=total_energy)
               deallocate (frequencies, reduced_masses, force_constants, cart_disp, fc_mdyne)
            end if
         end block
      end if

   end subroutine compute_gmbe