serial_fragment_processor Module Subroutine

   module subroutine serial_fragment_processor(total_fragments, polymers, max_level, &
                                               sys_geom, method_config, calc_type, json_data)
      !! Process all fragments serially in single-rank mode
      !! This is used when running with only 1 MPI rank
      !! Bond connectivity is accessed via sys_geom%bonds
      use mqc_error, only: error_t
      use mqc_result_types, only: mbe_result_t
      use mqc_json_output_types, only: json_output_data_t
      integer(int64), intent(in) :: total_fragments
      integer, intent(in) :: polymers(:, :), max_level
      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

      integer(int64) :: frag_idx
      integer :: fragment_size, current_log_level, iatom
      integer, allocatable :: fragment_indices(:)
      type(calculation_result_t), allocatable :: results(:)
      type(mbe_result_t) :: mbe_result
      type(physical_fragment_t) :: phys_frag
      type(timer_type) :: coord_timer
      integer(int32) :: calc_type_local
      type(error_t) :: error

      calc_type_local = calc_type

      call logger%info("Processing "//to_char(total_fragments)//" fragments serially...")
      call logger%info("  Calculation type: "//calc_type_to_string(calc_type_local))

      allocate (results(total_fragments))

      call omp_set_num_threads(1)
      call coord_timer%start()
      do frag_idx = 1_int64, total_fragments
         fragment_size = count(polymers(frag_idx, :) > 0)
         allocate (fragment_indices(fragment_size))
         fragment_indices = polymers(frag_idx, 1:fragment_size)

         call build_fragment_from_indices(sys_geom, fragment_indices, phys_frag, error, sys_geom%bonds)
         if (error%has_error()) then
            call logger%error(error%get_full_trace())
            error stop "Failed to build fragment in serial processing"
         end if

         call do_fragment_work(frag_idx, results(frag_idx), method_config, phys_frag, calc_type=calc_type_local)

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

         ! Debug output for gradients
         if (calc_type_local == CALC_TYPE_GRADIENT .and. results(frag_idx)%has_gradient) then
            call logger%configuration(level=current_log_level)
            if (current_log_level >= verbose_level) then
               block
                  character(len=512) :: debug_line
                  integer :: iatom_local
                  write (debug_line, '(a,i0,a,*(i0,1x))') "Fragment ", frag_idx, " monomers: ", fragment_indices
                  call logger%verbose(trim(debug_line))
                  write (debug_line, '(a,f25.15)') "  Energy: ", results(frag_idx)%energy%total()
                  call logger%verbose(trim(debug_line))
                  write (debug_line, '(a,f25.15)') "  Gradient norm: ", sqrt(sum(results(frag_idx)%gradient**2))
                  call logger%verbose(trim(debug_line))
                  if (size(results(frag_idx)%gradient, 2) <= 20) then
                     call logger%verbose("  Fragment gradient:")
                     do iatom_local = 1, size(results(frag_idx)%gradient, 2)
                        write (debug_line, '(a,i3,a,3f20.12)') "    Atom ", iatom_local, ": ", &
                           results(frag_idx)%gradient(1, iatom_local), &
                           results(frag_idx)%gradient(2, iatom_local), &
                           results(frag_idx)%gradient(3, iatom_local)
                        call logger%verbose(trim(debug_line))
                     end do
                  end if
               end block
            end if
         end if

         call phys_frag%destroy()
         deallocate (fragment_indices)

         if (mod(frag_idx, max(1_int64, total_fragments/10)) == 0 .or. frag_idx == total_fragments) then
            call logger%info("  Processed "//to_char(frag_idx)//"/"//to_char(total_fragments)// &
                             " fragments ["//to_char(coord_timer%get_elapsed_time())//" s]")
         end if
      end do
      call coord_timer%stop()
      call logger%info("Time to evaluate all fragments "//to_char(coord_timer%get_elapsed_time())//" s")
      call omp_set_num_threads(omp_get_max_threads())

      call logger%info("All fragments processed")

      call logger%info(" ")
      call logger%info("Computing Many-Body Expansion (MBE)...")
      call coord_timer%start()

      ! Allocate mbe_result components based on calc_type
      call mbe_result%allocate_dipole()  ! Always compute dipole
      if (calc_type_local == CALC_TYPE_HESSIAN) then
         call mbe_result%allocate_gradient(sys_geom%total_atoms)
         call mbe_result%allocate_hessian(sys_geom%total_atoms)
      else if (calc_type_local == CALC_TYPE_GRADIENT) then
         call mbe_result%allocate_gradient(sys_geom%total_atoms)
      end if

      call compute_mbe(polymers, total_fragments, max_level, results, mbe_result, sys_geom, json_data=json_data)
      call mbe_result%destroy()

      call coord_timer%stop()
      call logger%info("Time to compute MBE "//to_char(coord_timer%get_elapsed_time())//" s")

      deallocate (results)

   end subroutine serial_fragment_processor