frozen_double_entry Function

    function frozen_double_entry(indices, prefactor) result(idx)
        integer(int64), intent(in) :: indices(num_inds)
        real(dp), intent(out) :: prefactor
        integer(int64) :: idx(num_ex)
        integer :: f_one, f_two, ct, counts, i, rs, marks(3)
        logical :: unfrozen(num_inds)
        integer(int64) :: uf_idx(4)

        idx = 0
        unfrozen = indices > 0
        ! Mark where the unfrozen indices are
        marks = 0

        ! Identify the direct unfrozen pair
        ! Look for the direct pair of unfrozen indices (these shall stay direct)
        do ct = 1, step
            if(unfrozen(ct) .and. unfrozen(ct + step)) then
                uf_idx(1) = indices(ct)
                uf_idx(3) = indices(ct + step)
                ! Mark the extracted indices so they are not taken again
                marks(1) = ct
                marks(2) = ct + step
                exit
            end if
        end do
        ! Now, find the other two (direct or exchange) unfrozen indices
        rs = 2
        do ct = 1, num_inds
            if(unfrozen(ct) .and. .not. any(ct == marks)) then
                uf_idx(rs) = indices(ct)
                marks(3) = ct
                rs = rs + 2
            endif
        end do

        ! Now, get the required UMat indices, including all transposed index pairs
        idx = permute_umat_inds(int(uf_idx(1)), int(uf_idx(2)), int(uf_idx(3)), int(uf_idx(4)))
        ! Check the permutation and possible factor of two due to spin
        ! First, get the two frozen orbitals
        f_one = custom_findloc(unfrozen, .false., back=.false.)
        f_two = custom_findloc(unfrozen, .false., back=.true.)

        ! There are two spin configurations in a close-shell frozen scenario for terms
        ! with a direct frozen orbital (i.e. no permutation required)
        ! All others enter with a prefactor of -1
        if(is_direct(f_one, f_two)) then
            prefactor = 2.0_dp
        else
            prefactor = -1.0_dp
            ! If this is a diagonal term (only three different indices), permutational symmetry
            ! gives an extra factor of 2 for the exchange term. If only two different
            ! indices appear, also the direct term gets this prefactor
            counts = 1
            ! Count the number of different indices
            do ct = 2, num_inds
                if(.not. any(indices(ct) == indices(1:ct - 1))) counts = counts + 1
            end do
            if(counts == 2) then
                prefactor = 2.0_dp * prefactor
            elseif(counts == 3) then
                ! If there are three different indices, only add the extra factor if there is no direct unfrozen pair
                do ct = 1, step
                    if(is_repeated_pair(indices, ct) .and. unfrozen(ct)) return
                end do
                prefactor = 2.0_dp * prefactor
            end if
        end if

    end function frozen_double_entry