build_fragment_from_indices Subroutine

   subroutine build_fragment_from_indices(sys_geom, monomer_indices, fragment, error, bonds)
      !! Build a fragment on-the-fly from monomer indices with hydrogen capping for broken bonds
      !!
      !! Extracts atoms from specified monomers and adds hydrogen caps where bonds are broken.
      !! Caps are always added at the end of the atom list.
      !! Supports both fixed-size (identical monomers) and variable-sized fragments.
      !!
      !! Example: monomer_indices = [1, 3, 5] extracts waters 1, 3, and 5
      !!          If connectivity shows broken bonds, hydrogens are capped at positions of missing atoms
      type(system_geometry_t), intent(in) :: sys_geom
      integer, intent(in) :: monomer_indices(:)
      type(physical_fragment_t), intent(out) :: fragment
      type(error_t), intent(out) :: error
      type(bond_t), intent(in), optional :: bonds(:)  !! Connectivity information for capping

      integer :: n_monomers_in_frag, atoms_per_monomer, n_atoms_no_caps
      integer :: i, j, mono_idx, atom_start, atom_end, frag_atom_idx
      integer :: atom_i, atom_j, n_caps
      integer, allocatable :: atoms_in_fragment(:)  !! List of all atom indices in this fragment
      integer :: iatom, atom_global_idx
      logical :: use_explicit_fragments

      n_monomers_in_frag = size(monomer_indices)

      ! Determine if we're using explicit fragment definitions or regular monomer-based
      use_explicit_fragments = allocated(sys_geom%fragment_atoms)

      if (use_explicit_fragments) then
         ! Variable-sized fragments: count total atoms from fragment definitions
         n_atoms_no_caps = 0
         do i = 1, n_monomers_in_frag
            mono_idx = monomer_indices(i)
            n_atoms_no_caps = n_atoms_no_caps + sys_geom%fragment_sizes(mono_idx)
         end do

         ! Build list of atom indices (0-indexed) from explicit fragment definitions
         allocate (atoms_in_fragment(n_atoms_no_caps))
         iatom = 0
         do i = 1, n_monomers_in_frag
            mono_idx = monomer_indices(i)
            do j = 1, sys_geom%fragment_sizes(mono_idx)
               iatom = iatom + 1
               atoms_in_fragment(iatom) = sys_geom%fragment_atoms(j, mono_idx)
            end do
         end do
      else
         ! Fixed-size monomers: use atoms_per_monomer
         atoms_per_monomer = sys_geom%atoms_per_monomer
         n_atoms_no_caps = n_monomers_in_frag*atoms_per_monomer

         ! Build list of atom indices in this fragment (0-indexed to match bond indices)
         allocate (atoms_in_fragment(n_atoms_no_caps))
         iatom = 0
         do i = 1, n_monomers_in_frag
            mono_idx = monomer_indices(i)
            atom_start = (mono_idx - 1)*atoms_per_monomer
            do atom_i = 0, atoms_per_monomer - 1
               iatom = iatom + 1
               atoms_in_fragment(iatom) = atom_start + atom_i
            end do
         end do
      end if

      ! Count how many caps we need
      call count_hydrogen_caps(atoms_in_fragment, bonds, n_caps)

      ! Allocate arrays with space for original atoms + caps
      fragment%n_atoms = n_atoms_no_caps + n_caps
      fragment%n_caps = n_caps
      allocate (fragment%element_numbers(fragment%n_atoms))
      allocate (fragment%coordinates(3, fragment%n_atoms))
      if (n_caps > 0) allocate (fragment%cap_replaces_atom(n_caps))
      allocate (fragment%local_to_global(n_atoms_no_caps))  ! Only non-cap atoms

      ! Copy original atoms and build local→global mapping
      frag_atom_idx = 0
      if (use_explicit_fragments) then
         ! Variable-sized: copy atoms based on explicit fragment definitions
         do i = 1, n_monomers_in_frag
            mono_idx = monomer_indices(i)
            do j = 1, sys_geom%fragment_sizes(mono_idx)
               frag_atom_idx = frag_atom_idx + 1
               ! fragment_atoms is 0-indexed, so +1 for Fortran arrays
               atom_global_idx = sys_geom%fragment_atoms(j, mono_idx) + 1
               fragment%element_numbers(frag_atom_idx) = sys_geom%element_numbers(atom_global_idx)
               fragment%coordinates(:, frag_atom_idx) = sys_geom%coordinates(:, atom_global_idx)
               fragment%local_to_global(frag_atom_idx) = atom_global_idx  ! Store 1-indexed global position
            end do
         end do
      else
         ! Fixed-size: use atoms_per_monomer
         do i = 1, n_monomers_in_frag
            mono_idx = monomer_indices(i)
            atom_start = (mono_idx - 1)*atoms_per_monomer + 1
            atom_end = mono_idx*atoms_per_monomer

            ! Copy coordinates and elements
            do atom_i = atom_start, atom_end
               frag_atom_idx = frag_atom_idx + 1
               fragment%element_numbers(frag_atom_idx) = sys_geom%element_numbers(atom_i)
               fragment%coordinates(:, frag_atom_idx) = sys_geom%coordinates(:, atom_i)
               fragment%local_to_global(frag_atom_idx) = atom_i  ! Store 1-indexed global position
            end do
         end do
      end if

      ! Add hydrogen caps at end (if any)
      if (present(bonds) .and. n_caps > 0) then
         call add_hydrogen_caps(atoms_in_fragment, bonds, sys_geom, fragment, n_atoms_no_caps)
      end if

      ! Set electronic structure properties from system geometry
      if (use_explicit_fragments .and. allocated(sys_geom%fragment_charges) .and. &
          allocated(sys_geom%fragment_multiplicities)) then
         ! Explicit fragments: sum charges and multiplicities from constituent fragments
         fragment%charge = 0
         fragment%multiplicity = 1  ! Start with singlet assumption

         do i = 1, n_monomers_in_frag
            mono_idx = monomer_indices(i)
            fragment%charge = fragment%charge + sys_geom%fragment_charges(mono_idx)
         end do

         ! For single fragment, use its specific multiplicity
         if (n_monomers_in_frag == 1) then
            fragment%multiplicity = sys_geom%fragment_multiplicities(monomer_indices(1))
         else
            ! For multi-fragment composites, multiplicity needs careful treatment
            ! For now, default to system multiplicity (this may need refinement)
            fragment%multiplicity = sys_geom%multiplicity
         end if
      else
         ! Fixed-size monomers: use system defaults
         fragment%charge = sys_geom%charge
         fragment%multiplicity = sys_geom%multiplicity
      end if
      call fragment%compute_nelec()

      ! Validate: check for spatially overlapping atoms
      call check_duplicate_atoms(fragment, error)
      if (error%has_error()) then
         call error%add_context("mqc_physical_fragment:build_fragment_from_indices")
         return
      end if

      ! Calculate minimal distance between monomers in this fragment
      fragment%distance = calculate_monomer_distance(sys_geom, monomer_indices)

      deallocate (atoms_in_fragment)

   end subroutine build_fragment_from_indices