serial_gmbe_pie_processor Subroutine

   subroutine serial_gmbe_pie_processor(pie_atom_sets, pie_coefficients, n_pie_terms, &
                                        sys_geom, method_config, calc_type, json_data)
      !! Serial GMBE processor using PIE coefficients
      !! Evaluates each unique atom set once and sums with PIE coefficients
      !! Supports energy-only, energy+gradient, and energy+gradient+Hessian calculations
      !! If json_data is present, populates it for centralized JSON output
      !! Bond connectivity is accessed via sys_geom%bonds
      use mqc_calc_types, only: CALC_TYPE_GRADIENT, CALC_TYPE_HESSIAN, CALC_TYPE_ENERGY, calc_type_to_string
      use mqc_physical_fragment, only: redistribute_cap_gradients, redistribute_cap_hessian, &
                                       redistribute_cap_dipole_derivatives
      use mqc_error, only: error_t
      use pic_logger, only: info_level
      integer, intent(in) :: pie_atom_sets(:, :)  !! Unique atom sets (max_atoms, n_pie_terms)
      integer, intent(in) :: pie_coefficients(:)  !! PIE coefficient for each term
      integer(int64), intent(in) :: n_pie_terms
      type(system_geometry_t), intent(in) :: sys_geom
      type(method_config_t), intent(in) :: method_config  !! Method configuration
      integer(int32), intent(in) :: calc_type
      type(json_output_data_t), intent(out), optional :: json_data  !! JSON output data

      type(physical_fragment_t) :: phys_frag
      type(calculation_result_t), allocatable :: results(:)
      type(error_t) :: error
      integer :: n_atoms, max_atoms, iatom, current_log_level, hess_dim
      integer(int64) :: term_idx
      integer, allocatable :: atom_list(:)
      real(dp) :: total_energy, term_energy
      real(dp), allocatable :: pie_energies(:)  !! Store individual energies for JSON output
      real(dp), allocatable :: total_gradient(:, :)  !! Total gradient (3, total_atoms)
      real(dp), allocatable :: term_gradient(:, :)  !! Temporary gradient for each term
      real(dp), allocatable :: total_hessian(:, :)  !! Total Hessian (3*total_atoms, 3*total_atoms)
      real(dp), allocatable :: term_hessian(:, :)  !! Temporary Hessian for each term
      real(dp), allocatable :: total_dipole_derivs(:, :)  !! Total dipole derivatives (3, 3*total_atoms)
      real(dp), allocatable :: term_dipole_derivs(:, :)   !! Temporary dipole derivatives for each term
      real(dp), allocatable :: ir_intensities(:)          !! IR intensities in km/mol
      integer :: coeff

      if (int(size(pie_atom_sets, 2), int64) < n_pie_terms .or. &
          int(size(pie_coefficients), int64) < n_pie_terms) then
         call logger%error("PIE term arrays are smaller than n_pie_terms")
         error stop "Invalid PIE term array sizes"
      end if

      call logger%info("Processing "//to_char(n_pie_terms)//" unique PIE terms...")
      call logger%info("  Calculation type: "//calc_type_to_string(calc_type))

      total_energy = 0.0_dp
      max_atoms = size(pie_atom_sets, 1)
      allocate (pie_energies(n_pie_terms))
      allocate (results(n_pie_terms))

      ! Allocate gradient and Hessian arrays if needed
      if (calc_type == CALC_TYPE_GRADIENT .or. calc_type == CALC_TYPE_HESSIAN) then
         allocate (total_gradient(3, sys_geom%total_atoms))
         allocate (term_gradient(3, sys_geom%total_atoms))
         total_gradient = 0.0_dp
      end if

      if (calc_type == CALC_TYPE_HESSIAN) then
         hess_dim = 3*sys_geom%total_atoms
         allocate (total_hessian(hess_dim, hess_dim))
         allocate (term_hessian(hess_dim, hess_dim))
         total_hessian = 0.0_dp
         ! Allocate dipole derivative arrays for IR intensities
         allocate (total_dipole_derivs(3, hess_dim))
         allocate (term_dipole_derivs(3, hess_dim))
         total_dipole_derivs = 0.0_dp
      end if

      do term_idx = 1_int64, n_pie_terms
         coeff = pie_coefficients(term_idx)

         ! Skip terms with zero coefficient (shouldn't happen, but safety check)
         if (coeff == 0) then
            pie_energies(term_idx) = 0.0_dp  ! Mark as skipped
            cycle
         end if

         ! Extract atom list for this term
         n_atoms = 0
         do while (n_atoms < max_atoms .and. pie_atom_sets(n_atoms + 1, term_idx) >= 0)
            n_atoms = n_atoms + 1
         end do

         if (n_atoms == 0) then
            pie_energies(term_idx) = 0.0_dp  ! Mark as skipped
            cycle
         end if

         allocate (atom_list(n_atoms))
         atom_list = pie_atom_sets(1:n_atoms, term_idx)

         ! Build fragment from atom list
         call build_fragment_from_atom_list(sys_geom, atom_list, n_atoms, phys_frag, error, sys_geom%bonds)
         if (error%has_error()) then
            call logger%error(error%get_full_trace())
            error stop "Failed to build intersection fragment"
         end if

         ! Compute energy (and gradient if requested)
         call do_fragment_work(term_idx, results(term_idx), method_config, phys_frag, calc_type)

         ! Check for calculation errors
         if (results(term_idx)%has_error) then
            call logger%error("PIE term "//to_char(term_idx)//" calculation failed: "// &
                              results(term_idx)%error%get_message())
            error stop "PIE term calculation failed in serial processing"
         end if

         term_energy = results(term_idx)%energy%total()

         ! Store energy for JSON output
         pie_energies(term_idx) = term_energy

         ! Accumulate with PIE coefficient
         total_energy = total_energy + real(coeff, dp)*term_energy

         ! Accumulate gradient if present
         if ((calc_type == CALC_TYPE_GRADIENT .or. calc_type == CALC_TYPE_HESSIAN) .and. &
             results(term_idx)%has_gradient) then
            ! Map fragment gradient to system coordinates with proper cap handling
            term_gradient = 0.0_dp
            call redistribute_cap_gradients(phys_frag, results(term_idx)%gradient, term_gradient)

            ! Accumulate with PIE coefficient
            total_gradient = total_gradient + real(coeff, dp)*term_gradient
         end if

         ! Accumulate Hessian if present
         if (calc_type == CALC_TYPE_HESSIAN .and. results(term_idx)%has_hessian) then
            ! Map fragment Hessian to system coordinates with proper cap handling
            term_hessian = 0.0_dp
            call redistribute_cap_hessian(phys_frag, results(term_idx)%hessian, term_hessian)

            ! Accumulate with PIE coefficient
            total_hessian = total_hessian + real(coeff, dp)*term_hessian

            ! Accumulate dipole derivatives if present (for IR intensities)
            if (results(term_idx)%has_dipole_derivatives) then
               term_dipole_derivs = 0.0_dp
               call redistribute_cap_dipole_derivatives(phys_frag, results(term_idx)%dipole_derivatives, &
                                                        term_dipole_derivs)
               total_dipole_derivs = total_dipole_derivs + real(coeff, dp)*term_dipole_derivs
            end if
         end if

         call logger%verbose("PIE term "//to_char(term_idx)//"/"//to_char(n_pie_terms)// &
                             ": "//to_char(n_atoms)//" atoms, coeff="//to_char(coeff)// &
                             ", E="//to_char(term_energy))

         deallocate (atom_list)
         call phys_frag%destroy()
      end do

      call logger%info(" ")
      call logger%info("GMBE PIE calculation completed successfully")
      call logger%info("Final GMBE energy: "//to_char(total_energy)//" Hartree")

      ! Print gradient info if computed
      if (calc_type == CALC_TYPE_GRADIENT .or. calc_type == CALC_TYPE_HESSIAN) then
         call logger%info("GMBE PIE gradient computation completed")
         call logger%info("  Total gradient norm: "//to_char(sqrt(sum(total_gradient**2))))

         ! Print detailed gradient if info level and small system
         call logger%configuration(level=current_log_level)
         if (current_log_level >= info_level .and. sys_geom%total_atoms < 100) then
            call logger%info(" ")
            call logger%info("Total GMBE PIE Gradient (Hartree/Bohr):")
            do iatom = 1, sys_geom%total_atoms
               block
                  character(len=256) :: grad_line
                  write (grad_line, '(a,i5,a,3f20.12)') "  Atom ", iatom, ": ", &
                     total_gradient(1, iatom), total_gradient(2, iatom), total_gradient(3, iatom)
                  call logger%info(trim(grad_line))
               end block
            end do
            call logger%info(" ")
         end if
      end if

      ! Print Hessian info if computed
      if (calc_type == CALC_TYPE_HESSIAN) then
         call logger%info("GMBE PIE Hessian computation completed")
         call logger%info("  Total Hessian Frobenius norm: "//to_char(sqrt(sum(total_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
            type(mbe_result_t) :: gmbe_result
            integer :: n_at, n_modes

            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, &
                                              dipole_derivatives=total_dipole_derivs, &
                                              ir_intensities=ir_intensities)

            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
               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=total_energy)

               ! Build temporary mbe_result for JSON output
               gmbe_result%total_energy = total_energy
               gmbe_result%has_energy = .true.
               gmbe_result%has_hessian = .true.
               if (allocated(total_gradient)) then
                  gmbe_result%has_gradient = .true.
                  allocate (gmbe_result%gradient, source=total_gradient)
               end if
               allocate (gmbe_result%hessian, source=total_hessian)

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

                  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

                  if (allocated(total_gradient)) then
                     allocate (json_data%gradient, source=total_gradient)
                     json_data%has_gradient = .true.
                  end if

                  allocate (json_data%hessian, source=total_hessian)
                  json_data%has_hessian = .true.
               end if

               if (allocated(ir_intensities)) deallocate (ir_intensities)
               call gmbe_result%destroy()
               deallocate (frequencies, reduced_masses, force_constants, cart_disp, fc_mdyne)
            end if
         end block
      end if

      call logger%info(" ")

      ! Populate json_data for non-Hessian case if present
      if (present(json_data) .and. calc_type /= CALC_TYPE_HESSIAN) then
         json_data%output_mode = OUTPUT_MODE_GMBE_PIE
         json_data%total_energy = total_energy
         json_data%has_energy = .true.
         json_data%n_pie_terms = n_pie_terms

         ! Copy PIE data
         allocate (json_data%pie_atom_sets, source=pie_atom_sets(:, 1:n_pie_terms))
         allocate (json_data%pie_coefficients(n_pie_terms))
         json_data%pie_coefficients = pie_coefficients(1:n_pie_terms)
         allocate (json_data%pie_energies(n_pie_terms))
         json_data%pie_energies = pie_energies

         if (allocated(total_gradient)) then
            allocate (json_data%gradient, source=total_gradient)
            json_data%has_gradient = .true.
         end if
         if (allocated(total_hessian)) then
            allocate (json_data%hessian, source=total_hessian)
            json_data%has_hessian = .true.
         end if
      end if

      deallocate (pie_energies, results)
      if (allocated(total_gradient)) deallocate (total_gradient)
      if (allocated(term_gradient)) deallocate (term_gradient)
      if (allocated(total_hessian)) deallocate (total_hessian)
      if (allocated(term_hessian)) deallocate (term_hessian)
      if (allocated(total_dipole_derivs)) deallocate (total_dipole_derivs)
      if (allocated(term_dipole_derivs)) deallocate (term_dipole_derivs)

   end subroutine serial_gmbe_pie_processor