compute_mbe Subroutine

   subroutine compute_mbe(polymers, fragment_count, max_level, results, &
                          mbe_result, sys_geom, world_comm, json_data)
      !! Compute many-body expansion (MBE) energy with optional gradient, hessian, and dipole
      !!
      !! This is the core routine that handles all MBE computations.
      !! The caller pre-allocates desired components in mbe_result before calling:
      !!   - mbe_result%gradient allocated: compute gradient (requires sys_geom)
      !!   - mbe_result%hessian allocated: compute hessian (requires sys_geom)
      !!   - mbe_result%dipole allocated: compute total dipole moment
      !! If json_data is present, populates it for centralized JSON output
      !! Bond connectivity is accessed via sys_geom%bonds
      use mqc_result_types, only: calculation_result_t, mbe_result_t

      ! Required arguments
      integer(int64), intent(in) :: fragment_count
      integer, intent(in) :: polymers(:, :), max_level
      type(calculation_result_t), intent(in) :: results(:)
      type(mbe_result_t), intent(inout) :: mbe_result  !! Pre-allocated by caller

      ! Optional arguments
      type(system_geometry_t), intent(in), optional :: sys_geom  !! Required for gradient/hessian
      type(comm_t), intent(in), optional :: world_comm           !! MPI communicator for abort
      type(json_output_data_t), intent(out), optional :: json_data  !! JSON output data

      ! Local variables
      integer(int64) :: i
      integer :: fragment_size, nlevel, current_log_level, hess_dim
      real(dp), allocatable :: sum_by_level(:), delta_energies(:), energies(:)
      real(dp), allocatable :: delta_gradients(:, :, :), delta_hessians(:, :, :)
      real(dp), allocatable :: delta_dipoles(:, :)  !! (3, fragment_count)
      real(dp), allocatable :: delta_dipole_derivs(:, :, :)  !! (3, 3*total_atoms, fragment_count)
      real(dp), allocatable :: ir_intensities(:)  !! IR intensities in km/mol
      real(dp) :: delta_E
      logical :: do_detailed_print, compute_grad, compute_hess, compute_dipole, compute_dipole_derivs
      type(fragment_lookup_t) :: lookup

      ! Determine what to compute based on allocated components in mbe_result
      compute_grad = allocated(mbe_result%gradient)
      compute_hess = allocated(mbe_result%hessian)
      compute_dipole = allocated(mbe_result%dipole)
      compute_dipole_derivs = .false.  ! Will be set true if all fragments have dipole_derivatives

      ! Validate inputs for gradient computation
      if (compute_grad) then
         do i = 1_int64, fragment_count
            if (.not. results(i)%has_gradient) then
               call logger%error("Fragment "//to_char(i)//" does not have gradient!")
               if (present(world_comm)) then
                  call abort_comm(world_comm, 1)
               else
                  error stop "Missing gradient in compute_mbe_internal"
               end if
            end if
         end do
      end if

      ! Validate inputs for hessian computation
      if (compute_hess) then
         do i = 1_int64, fragment_count
            if (.not. results(i)%has_hessian) then
               call logger%error("Fragment "//to_char(i)//" does not have Hessian!")
               if (present(world_comm)) then
                  call abort_comm(world_comm, 1)
               else
                  error stop "Missing Hessian in compute_mbe_internal"
               end if
            end if
         end do
         hess_dim = 3*sys_geom%total_atoms
      end if

      ! Validate inputs for dipole computation
      if (compute_dipole) then
         do i = 1_int64, fragment_count
            if (.not. results(i)%has_dipole) then
               call logger%warning("Fragment "//to_char(i)//" does not have dipole - skipping dipole MBE")
               compute_dipole = .false.
               exit
            end if
         end do
      end if

      ! Check if dipole derivatives are available (for IR intensities)
      if (compute_hess) then
         compute_dipole_derivs = .true.
         do i = 1_int64, fragment_count
            if (.not. results(i)%has_dipole_derivatives) then
               compute_dipole_derivs = .false.
               exit
            end if
         end do
      end if

      ! Get logger configuration
      call logger%configuration(level=current_log_level)
      do_detailed_print = (current_log_level >= verbose_level)

      ! Allocate energy arrays (always needed)
      allocate (sum_by_level(max_level))
      allocate (delta_energies(fragment_count))
      allocate (energies(fragment_count))
      sum_by_level = 0.0_dp
      delta_energies = 0.0_dp

      ! Extract total energies from results
      do i = 1_int64, fragment_count
         energies(i) = results(i)%energy%total()
      end do

      ! Allocate gradient delta arrays if needed
      if (compute_grad) then
         allocate (delta_gradients(3, sys_geom%total_atoms, fragment_count))
         delta_gradients = 0.0_dp
         mbe_result%gradient = 0.0_dp
      end if

      ! Allocate hessian delta arrays if needed
      if (compute_hess) then
         allocate (delta_hessians(hess_dim, hess_dim, fragment_count))
         delta_hessians = 0.0_dp
         mbe_result%hessian = 0.0_dp
      end if

      ! Allocate dipole delta arrays if needed
      if (compute_dipole) then
         allocate (delta_dipoles(3, fragment_count))
         delta_dipoles = 0.0_dp
         mbe_result%dipole = 0.0_dp
      end if

      ! Allocate dipole derivative delta arrays if needed (for IR intensities)
      if (compute_dipole_derivs) then
         allocate (delta_dipole_derivs(3, hess_dim, fragment_count))
         delta_dipole_derivs = 0.0_dp
         allocate (mbe_result%dipole_derivatives(3, hess_dim))
         mbe_result%dipole_derivatives = 0.0_dp
      end if

      ! Build hash table for fast fragment lookups
      block
         use mqc_error, only: error_t
         type(error_t) :: lookup_error
         call build_mbe_lookup_table(polymers, fragment_count, max_level, lookup, lookup_error)
         if (lookup_error%has_error()) then
            call logger%error("Failed to build lookup table: "//lookup_error%get_message())
            if (present(world_comm)) then
               call abort_comm(world_comm, 1)
            else
               error stop "Failed to build lookup table"
            end if
         end if
      end block

      ! Bottom-up computation: process fragments by size (1-body, then 2-body, etc.)
      ! This makes the algorithm independent of input fragment order
      ! We process by n-mer level to ensure all subsets are computed before they're needed
      do nlevel = 1, max_level
         do i = 1_int64, fragment_count
            fragment_size = count(polymers(i, :) > 0)

            ! Only process fragments of the current nlevel
            if (fragment_size /= nlevel) cycle

            if (fragment_size == 1) then
               ! 1-body: delta = value (no subsets to subtract)
               delta_energies(i) = energies(i)
               sum_by_level(1) = sum_by_level(1) + delta_energies(i)

               if (compute_grad) then
                  call map_fragment_to_system_gradient(results(i)%gradient, &
                                             polymers(i, 1:fragment_size), sys_geom, delta_gradients(:, :, i), world_comm)
               end if

               if (compute_hess) then
                  call map_fragment_to_system_hessian(results(i)%hessian, &
                                                      polymers(i, 1:fragment_size), sys_geom, delta_hessians(:, :, i))
               end if

               if (compute_dipole) then
                  ! For 1-body, delta dipole is just the fragment dipole
                  delta_dipoles(:, i) = results(i)%dipole
               end if

               if (compute_dipole_derivs) then
                  ! For 1-body, delta dipole derivatives are just the fragment values mapped to system
                  call map_fragment_to_system_dipole_derivatives(results(i)%dipole_derivatives, &
                                                     polymers(i, 1:fragment_size), sys_geom, delta_dipole_derivs(:, :, i))
               end if

            else if (fragment_size >= 2 .and. fragment_size <= max_level) then
               ! n-body: delta = value - sum(all subset deltas)
               delta_E = compute_mbe_delta(i, polymers(i, 1:fragment_size), lookup, &
                                           energies, delta_energies, fragment_size, world_comm)
               delta_energies(i) = delta_E
               sum_by_level(fragment_size) = sum_by_level(fragment_size) + delta_E

               if (compute_grad) then
                  call compute_mbe_gradient(i, polymers(i, 1:fragment_size), lookup, &
                                            results, delta_gradients, fragment_size, sys_geom, world_comm)
               end if

               if (compute_hess) then
                  call compute_mbe_hessian(i, polymers(i, 1:fragment_size), lookup, &
                                           results, delta_hessians, fragment_size, sys_geom)
               end if

               if (compute_dipole) then
                  call compute_mbe_dipole(i, polymers(i, 1:fragment_size), lookup, &
                                          results, delta_dipoles, fragment_size, world_comm)
               end if

               if (compute_dipole_derivs) then
                  call compute_mbe_dipole_derivatives(i, polymers(i, 1:fragment_size), lookup, &
                                                      results, delta_dipole_derivs, fragment_size, sys_geom)
               end if
            end if
         end do
      end do

      ! Clean up lookup table
      call lookup%destroy()

      ! Compute totals and set status flags
      mbe_result%total_energy = sum(sum_by_level)
      mbe_result%has_energy = .true.

      if (compute_grad) then
         do i = 1_int64, fragment_count
            fragment_size = count(polymers(i, :) > 0)
            if (fragment_size <= max_level) then
               mbe_result%gradient = mbe_result%gradient + delta_gradients(:, :, i)
            end if
         end do
         mbe_result%has_gradient = .true.
      end if

      if (compute_hess) then
         do i = 1_int64, fragment_count
            fragment_size = count(polymers(i, :) > 0)
            if (fragment_size <= max_level) then
               mbe_result%hessian = mbe_result%hessian + delta_hessians(:, :, i)
            end if
         end do
         mbe_result%has_hessian = .true.
      end if

      if (compute_dipole) then
         do i = 1_int64, fragment_count
            fragment_size = count(polymers(i, :) > 0)
            if (fragment_size <= max_level) then
               mbe_result%dipole = mbe_result%dipole + delta_dipoles(:, i)
            end if
         end do
         mbe_result%has_dipole = .true.
      end if

      if (compute_dipole_derivs) then
         do i = 1_int64, fragment_count
            fragment_size = count(polymers(i, :) > 0)
            if (fragment_size <= max_level) then
               mbe_result%dipole_derivatives = mbe_result%dipole_derivatives + delta_dipole_derivs(:, :, i)
            end if
         end do
         mbe_result%has_dipole_derivatives = .true.
      end if

      ! Print energy breakdown (always)
      call print_mbe_energy_breakdown(sum_by_level, max_level, mbe_result%total_energy)

      ! Print gradient info if computed
      if (compute_grad) then
         call print_mbe_gradient_info(mbe_result%gradient, sys_geom, current_log_level)
      end if

      ! Print hessian info if computed
      if (compute_hess) then
         call logger%info("MBE Hessian computation completed")
         call logger%info("  Total Hessian Frobenius norm: "//to_char(sqrt(sum(mbe_result%hessian**2))))

         ! Compute and print full vibrational analysis with thermochemistry
         block
            real(dp), allocatable :: frequencies(:), reduced_masses(:), force_constants(:)
            real(dp), allocatable :: cart_disp(:, :), fc_mdyne(:)
            type(thermochemistry_result_t) :: thermo_result
            integer :: n_at, n_modes

            call logger%info("  Computing vibrational analysis (projecting trans/rot modes)...")
            if (compute_dipole_derivs) then
               call compute_vibrational_analysis(mbe_result%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, &
                                                 dipole_derivatives=mbe_result%dipole_derivatives, &
                                                 ir_intensities=ir_intensities)
            else
               call compute_vibrational_analysis(mbe_result%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)
            end if

            if (allocated(frequencies)) then
               ! Compute thermochemistry
               n_at = size(sys_geom%element_numbers)
               n_modes = size(frequencies)
               call compute_thermochemistry(sys_geom%coordinates, sys_geom%element_numbers, &
                                            frequencies, n_at, n_modes, thermo_result)

               ! Print vibrational analysis to log
               if (allocated(ir_intensities)) then
                  call print_vibrational_analysis(frequencies, reduced_masses, force_constants, &
                                                  cart_disp, sys_geom%element_numbers, &
                                                  force_constants_mdyne=fc_mdyne, &
                                                  ir_intensities=ir_intensities, &
                                                  coordinates=sys_geom%coordinates, &
                                                  electronic_energy=mbe_result%total_energy)
               else
                  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=mbe_result%total_energy)
               end if

               ! Populate json_data for vibrational output if present
               if (present(json_data)) then
                  json_data%output_mode = OUTPUT_MODE_MBE
                  json_data%total_energy = mbe_result%total_energy
                  json_data%has_energy = mbe_result%has_energy
                  json_data%has_vibrational = .true.

                  ! Copy vibrational data
                  allocate (json_data%frequencies(n_modes))
                  allocate (json_data%reduced_masses(n_modes))
                  allocate (json_data%force_constants(n_modes))
                  json_data%frequencies = frequencies
                  json_data%reduced_masses = reduced_masses
                  json_data%force_constants = fc_mdyne
                  json_data%thermo = thermo_result

                  if (allocated(ir_intensities)) then
                     allocate (json_data%ir_intensities(n_modes))
                     json_data%ir_intensities = ir_intensities
                     json_data%has_ir_intensities = .true.
                  end if

                  ! Copy dipole if available
                  if (mbe_result%has_dipole) then
                     allocate (json_data%dipole(3))
                     json_data%dipole = mbe_result%dipole
                     json_data%has_dipole = .true.
                  end if

                  ! Copy gradient if available
                  if (mbe_result%has_gradient) then
                     allocate (json_data%gradient, source=mbe_result%gradient)
                     json_data%has_gradient = .true.
                  end if

                  ! Copy hessian if available
                  if (mbe_result%has_hessian) then
                     allocate (json_data%hessian, source=mbe_result%hessian)
                     json_data%has_hessian = .true.
                  end if
               end if

               if (allocated(ir_intensities)) deallocate (ir_intensities)
               deallocate (frequencies, reduced_masses, force_constants, cart_disp, fc_mdyne)
            end if
         end block
      end if

      ! Print dipole info if computed
      if (compute_dipole) then
         block
            character(len=256) :: dipole_line
            real(dp) :: dipole_magnitude
            dipole_magnitude = norm2(mbe_result%dipole)*2.541746_dp  ! Convert e*Bohr to Debye
            call logger%info("MBE Dipole moment:")
            write (dipole_line, '(a,3f15.8)') "  Dipole (e*Bohr): ", mbe_result%dipole
            call logger%info(trim(dipole_line))
            write (dipole_line, '(a,f15.8)') "  Dipole magnitude (Debye): ", dipole_magnitude
            call logger%info(trim(dipole_line))
         end block
      end if

      ! Print detailed breakdown if requested
      if (do_detailed_print) then
         call print_detailed_breakdown(polymers, fragment_count, max_level, energies, delta_energies)
      end if

      ! Populate json_data for non-Hessian case if present
      ! (Hessian case already handled above in the vibrational block)
      if (present(json_data) .and. .not. compute_hess) then
         json_data%output_mode = OUTPUT_MODE_MBE
         json_data%total_energy = mbe_result%total_energy
         json_data%has_energy = mbe_result%has_energy
         json_data%max_level = max_level
         json_data%fragment_count = fragment_count

         ! Copy fragment breakdown data
         allocate (json_data%polymers(fragment_count, max_level))
         json_data%polymers = polymers(1:fragment_count, 1:max_level)

         allocate (json_data%fragment_energies(fragment_count))
         json_data%fragment_energies = energies

         allocate (json_data%delta_energies(fragment_count))
         json_data%delta_energies = delta_energies

         allocate (json_data%sum_by_level(max_level))
         json_data%sum_by_level = sum_by_level

         ! Copy fragment distances if available
         allocate (json_data%fragment_distances(fragment_count))
         do i = 1_int64, fragment_count
            json_data%fragment_distances(i) = results(i)%distance
         end do

         ! Copy dipole if available
         if (mbe_result%has_dipole) then
            allocate (json_data%dipole(3))
            json_data%dipole = mbe_result%dipole
            json_data%has_dipole = .true.
         end if

         ! Copy gradient if available
         if (mbe_result%has_gradient) then
            allocate (json_data%gradient, source=mbe_result%gradient)
            json_data%has_gradient = .true.
         end if
      end if

      ! Cleanup
      deallocate (sum_by_level, delta_energies, energies)
      if (allocated(delta_gradients)) deallocate (delta_gradients)
      if (allocated(delta_hessians)) deallocate (delta_hessians)
      if (allocated(delta_dipoles)) deallocate (delta_dipoles)
      if (allocated(delta_dipole_derivs)) deallocate (delta_dipole_derivs)

   end subroutine compute_mbe