run_calculation Subroutine

   subroutine run_calculation(resources, config, sys_geom, bonds, result_out, all_ranks_write_json)
      !! Main calculation dispatcher - routes to fragmented or unfragmented calculation
      !!
      !! Determines calculation type based on nlevel and dispatches to appropriate
      !! calculation routine with proper MPI setup and validation.
      !! If result_out is present, returns result instead of writing JSON (for dynamics/optimization)
      type(resources_t), intent(in) :: resources  !! Resources container (MPI comms, etc.)
      type(driver_config_t), intent(in) :: config  !! Driver configuration
      type(system_geometry_t), intent(in) :: sys_geom  !! System geometry and fragment info
      type(bond_t), intent(in), optional :: bonds(:)  !! Bond connectivity information
      type(calculation_result_t), intent(out), optional :: result_out  !! Optional result output
      logical, intent(in), optional :: all_ranks_write_json  !! If true, all ranks write JSON (for multi-molecule)

      ! Local variables
      integer :: max_level   !! Maximum fragment level (nlevel from config)
      integer :: i  !! Loop counter
      type(json_output_data_t) :: json_data  !! Cached output data for centralized JSON writing
      logical :: should_write_json  !! Whether this rank should write JSON

      ! Set max_level from config
      max_level = config%nlevel

      ! Log method-specific settings (rank 0 only)
      if (resources%mpi_comms%world_comm%rank() == 0) then
         call config%method_config%log_settings()
      end if

      if (resources%mpi_comms%world_comm%rank() == 0 .and. max_level > 0) then
         call logger%info("============================================")
         call logger%info("Loaded geometry:")
         call logger%info("  Total monomers: "//to_char(sys_geom%n_monomers))
         call logger%info("  Atoms per monomer: "//to_char(sys_geom%atoms_per_monomer))
         call logger%info("  Fragment level: "//to_char(max_level))
         call logger%info("  Total atoms: "//to_char(sys_geom%total_atoms))
         call logger%info("============================================")
      end if

      ! Warn if overlapping fragments flag is set but nlevel=0
      if (config%allow_overlapping_fragments .and. max_level == 0) then
         if (resources%mpi_comms%world_comm%rank() == 0) then
            call logger%warning("allow_overlapping_fragments is set to true, but nlevel=0")
            call logger%warning("Running unfragmented calculation - overlapping fragments flag will be ignored")
         end if
      end if

      ! GMBE (overlapping fragments) with inclusion-exclusion principle
      ! GMBE(1): Base fragments are monomers
      ! GMBE(N): Base fragments are N-mers (e.g., dimers for N=2)
      ! Algorithm: Generate primaries, use DFS to enumerate overlapping cliques,
      ! accumulate PIE coefficients per unique atom set, evaluate each once

      if (max_level == 0) then
         call omp_set_num_threads(1)
         if (present(result_out)) then
            ! For dynamics/optimization: return result directly, no JSON output
            call run_unfragmented_calculation(resources%mpi_comms%world_comm, sys_geom, config, result_out)
         else
            ! Normal mode: collect json_data for centralized output
            call run_unfragmented_calculation(resources%mpi_comms%world_comm, sys_geom, config, json_data=json_data)
         end if
      else
         if (present(result_out)) then
            ! For fragmented calculations with result_out (future use)
            call run_fragmented_calculation(resources, config, sys_geom, bonds)
         else
            ! Normal mode: collect json_data for centralized output
            call run_fragmented_calculation(resources, config, sys_geom, bonds, json_data)
         end if
      end if

      ! Centralized JSON output (rank 0 only by default, or all ranks if all_ranks_write_json is set)
      if (.not. present(result_out)) then
         ! Check if JSON output should be skipped
         if (config%skip_json_output) then
            if (resources%mpi_comms%world_comm%rank() == 0) then
               call logger%info("Skipping JSON output (skip_json_output = true)")
            end if
         else
            ! Determine if this rank should write JSON
            should_write_json = (resources%mpi_comms%world_comm%rank() == 0)
            if (present(all_ranks_write_json)) then
               if (all_ranks_write_json) should_write_json = .true.
            end if

            if (should_write_json) then
               if (json_data%output_mode /= OUTPUT_MODE_NONE) then
                  call write_json_output(json_data)
                  call json_data%destroy()
               end if
            end if
         end if
      end if

   end subroutine run_calculation