calculate_monomer_distance Function

   pure function calculate_monomer_distance(sys_geom, monomer_indices) result(min_distance)
      !! Calculate minimal atomic distance between monomers in a fragment
      !! For single monomer (size 1), returns 0.0
      !! For multi-monomer fragments, returns minimal distance between atoms in different monomers
      !! Result is in Angstrom
      type(system_geometry_t), intent(in) :: sys_geom
      integer, intent(in) :: monomer_indices(:)
      real(dp) :: min_distance

      integer :: n_monomers, i, j, iatom, jatom
      integer :: mon_i, mon_j
      integer :: atom_start_i, atom_end_i, atom_start_j, atom_end_j
      real(dp) :: dist, dx, dy, dz
      logical :: is_variable_size

      n_monomers = size(monomer_indices)

      ! Monomers have distance 0
      if (n_monomers == 1) then
         min_distance = 0.0_dp
         return
      end if

      ! Check if we have variable-sized fragments
      is_variable_size = allocated(sys_geom%fragment_sizes)

      ! Initialize with huge value
      min_distance = huge(1.0_dp)

      ! Loop over all pairs of monomers
      do i = 1, n_monomers - 1
         mon_i = monomer_indices(i)
         do j = i + 1, n_monomers
            mon_j = monomer_indices(j)

            if (is_variable_size) then
               ! Variable-sized fragments
               do iatom = 1, sys_geom%fragment_sizes(mon_i)
                  atom_start_i = sys_geom%fragment_atoms(iatom, mon_i) + 1  ! Convert to 1-indexed
                  do jatom = 1, sys_geom%fragment_sizes(mon_j)
                     atom_start_j = sys_geom%fragment_atoms(jatom, mon_j) + 1  ! Convert to 1-indexed

                     ! Calculate distance (coordinates in Bohr)
                     dx = sys_geom%coordinates(1, atom_start_i) - sys_geom%coordinates(1, atom_start_j)
                     dy = sys_geom%coordinates(2, atom_start_i) - sys_geom%coordinates(2, atom_start_j)
                     dz = sys_geom%coordinates(3, atom_start_i) - sys_geom%coordinates(3, atom_start_j)
                     dist = sqrt(dx*dx + dy*dy + dz*dz)

                     if (dist < min_distance) min_distance = dist
                  end do
               end do
            else
               ! Fixed-sized monomers
               atom_start_i = (mon_i - 1)*sys_geom%atoms_per_monomer + 1
               atom_end_i = mon_i*sys_geom%atoms_per_monomer

               atom_start_j = (mon_j - 1)*sys_geom%atoms_per_monomer + 1
               atom_end_j = mon_j*sys_geom%atoms_per_monomer

               ! Loop over all atom pairs
               do iatom = atom_start_i, atom_end_i
                  do jatom = atom_start_j, atom_end_j
                     ! Calculate distance (coordinates in Bohr)
                     dx = sys_geom%coordinates(1, iatom) - sys_geom%coordinates(1, jatom)
                     dy = sys_geom%coordinates(2, iatom) - sys_geom%coordinates(2, jatom)
                     dz = sys_geom%coordinates(3, iatom) - sys_geom%coordinates(3, jatom)
                     dist = sqrt(dx*dx + dy*dy + dz*dz)

                     if (dist < min_distance) min_distance = dist
                  end do
               end do
            end if
         end do
      end do

      ! Convert from Bohr to Angstrom
      min_distance = to_angstrom(min_distance)

   end function calculate_monomer_distance