mimic_sort_fragments Subroutine

subroutine mimic_sort_fragments(subsystem, quantum_fragment, cutoff_distance, sorting_method)

    intrinsic :: norm2

    !> array of subsystems associated with client codes
    type(subsystem_type), intent(inout) :: subsystem
    !> the quantum fragment
    type(quantum_fragment_type), intent(in) :: quantum_fragment
    !> distance cutoff
    real(dp), intent(in) :: cutoff_distance
    !> sorting method (centroid, center_of_mass, center_of_charge, minimum_distance)
    integer, intent(in) :: sorting_method

    integer :: i, j, k
    integer :: num_atoms
    integer, dimension(:), allocatable :: sr_indices
    integer, dimension(:), allocatable :: lr_indices
    integer, dimension(:), allocatable :: temporary_indices
    real(dp) :: distance
    real(dp), dimension(3) :: fragment_center
    real(dp), dimension(3) :: qm_fragment_center

    call timer_start("mimic_sort_fragments")

    if (allocated(subsystem%sr_atoms)) then
        deallocate(subsystem%sr_atoms)
    end if
    if (allocated(subsystem%lr_atoms)) then
        deallocate(subsystem%lr_atoms)
    end if
    if (allocated(subsystem%sr_fragments)) then
        deallocate(subsystem%sr_fragments)
    end if
    if (allocated(subsystem%lr_fragments)) then
        deallocate(subsystem%lr_fragments)
    end if

    allocate(sr_indices(0))
    allocate(lr_indices(0))

    select case (sorting_method)
    ! sort according to distance between centroids
    case (CENTROID)
        qm_fragment_center = quantum_fragment%centroid()
        do i = 1, subsystem%num_fragments
            fragment_center = subsystem%fragments(i)%centroid()
            distance = norm2(fragment_center - qm_fragment_center)
            if (distance <= cutoff_distance) then
                call move_alloc(sr_indices, temporary_indices)
                allocate(sr_indices, source=[temporary_indices, i])
            else
                call move_alloc(lr_indices, temporary_indices)
                allocate(lr_indices, source=[temporary_indices, i])
            end if
        end do
    ! sort according to distance between center-of-mass
    case (CENTER_OF_MASS)
        qm_fragment_center = quantum_fragment%center_of_mass()
        do i = 1, subsystem%num_fragments
            fragment_center = subsystem%fragments(i)%center_of_mass()
            distance = norm2(fragment_center - qm_fragment_center)
            if (distance <= cutoff_distance) then
                call move_alloc(sr_indices, temporary_indices)
                allocate(sr_indices, source=[temporary_indices, i])
            else
                call move_alloc(lr_indices, temporary_indices)
                allocate(lr_indices, source=[temporary_indices, i])
            end if
        end do
    ! sort according to distance between center-of-charge for charged fragments and center-of-mass for other fragments
    case (CENTER_OF_CHARGE)
        qm_fragment_center = quantum_fragment%center_of_charge()
        do i = 1, subsystem%num_fragments
            if (nint(subsystem%fragments(i)%charge) /= 0) then
                fragment_center = subsystem%fragments(i)%center_of_charge()
            else
                fragment_center = subsystem%fragments(i)%center_of_mass()
            end if
            distance = norm2(fragment_center - qm_fragment_center)
            if (distance <= cutoff_distance) then
                call move_alloc(sr_indices, temporary_indices)
                allocate(sr_indices, source=[temporary_indices, i])
            else
                call move_alloc(lr_indices, temporary_indices)
                allocate(lr_indices, source=[temporary_indices, i])
            end if
        end do
    ! sort according to minimum distance
    case (MINIMUM_DISTANCE)
        do i = 1, subsystem%num_fragments
            distance = compute_minimum_distance(subsystem%fragments(i)%atoms, quantum_fragment%nuclei)
            if (distance <= cutoff_distance) then
                call move_alloc(sr_indices, temporary_indices)
                allocate(sr_indices, source=[temporary_indices, i])
            else
                call move_alloc(lr_indices, temporary_indices)
                allocate(lr_indices, source=[temporary_indices, i])
            end if
        end do
    ! sort individual atoms into one short-range and one long-range fragment
    case (ATOM_WISE)
        qm_fragment_center = quantum_fragment%centroid()
        do i = 1, subsystem%num_atoms
            distance = norm2(qm_fragment_center - subsystem%atoms(i)%coordinate)
            if (distance <= cutoff_distance) then
                call move_alloc(sr_indices, temporary_indices)
                allocate(sr_indices, source=[temporary_indices, i])
            else
                call move_alloc(lr_indices, temporary_indices)
                allocate(lr_indices, source=[temporary_indices, i])
            end if
        end do
    case default
        call handle_error(SEVERITY_FATAL, &
                          TYPE_INCORRECT_ARG, &
                          "Unrecognized atom sorting method", &
                          __FILE__, __LINE__)
    end select

    if (sorting_method == ATOM_WISE) then
        subsystem%num_sr_fragments = 1
        allocate(subsystem%sr_fragments(subsystem%num_sr_fragments))
        subsystem%num_sr_atoms = size(sr_indices)
        allocate(subsystem%sr_atoms(size(sr_indices)))
        do i = 1, subsystem%num_sr_atoms
            subsystem%sr_atoms(i) = subsystem%atoms(sr_indices(i))
        end do
        call subsystem%sr_fragments(1)%init(1, subsystem%sr_atoms)
        subsystem%num_lr_fragments = 1
        allocate(subsystem%lr_fragments(subsystem%num_lr_fragments))
        subsystem%num_lr_atoms = size(lr_indices)
        allocate(subsystem%lr_atoms(size(lr_indices)))
        do i = 1, subsystem%num_lr_atoms
            subsystem%lr_atoms(i) = subsystem%atoms(lr_indices(i))
        end do
        call subsystem%lr_fragments(1)%init(1, subsystem%lr_atoms)
    else
        subsystem%num_sr_fragments = size(sr_indices)
        allocate(subsystem%sr_fragments(subsystem%num_sr_fragments))
        num_atoms = 0
        do i = 1, subsystem%num_sr_fragments
            subsystem%sr_fragments(i) = subsystem%fragments(sr_indices(i))
            num_atoms = num_atoms + subsystem%sr_fragments(i)%num_atoms
        end do
        subsystem%num_sr_atoms = num_atoms
        allocate(subsystem%sr_atoms(subsystem%num_sr_atoms))
        k = 0
        do i = 1, subsystem%num_sr_fragments
            j = k + 1
            k = k + subsystem%sr_fragments(i)%num_atoms
            subsystem%sr_atoms(j:k) = subsystem%sr_fragments(i)%atoms
        end do
        subsystem%num_lr_fragments = size(lr_indices)
        allocate(subsystem%lr_fragments(subsystem%num_lr_fragments))
        num_atoms = 0
        do i = 1, subsystem%num_lr_fragments
            subsystem%lr_fragments(i) = subsystem%fragments(lr_indices(i))
            num_atoms = num_atoms + subsystem%lr_fragments(i)%num_atoms
        end do
        subsystem%num_lr_atoms = num_atoms
        allocate(subsystem%lr_atoms(subsystem%num_lr_atoms))
        k = 0
        do i = 1, subsystem%num_lr_fragments
            j = k + 1
            k = k + subsystem%lr_fragments(i)%num_atoms
            subsystem%lr_atoms(j:k) = subsystem%lr_fragments(i)%atoms
        end do
    end if

    call timer_stop

end subroutine mimic_sort_fragments