frozen_diagonal_entry Subroutine

    pure subroutine frozen_diagonal_entry(indices, prefactor)
        integer(int64), intent(in) :: indices(num_inds)
        real(dp), intent(out) :: prefactor

        integer :: ct, directs
        logical :: t_quad
        ! Only one info is needed here, that is the parity and number of contributing spin
        ! configs.

        ! If there are five or more repeated indices, the LMat entry does not contribute
        ! (it would require three same-spin electrons
        t_quad = .false.
        do ct = 1, 3
            if(count(indices(ct) == indices) > 4) then
                prefactor = 0.0_dp
                return
            end if
            if(count(indices(ct) == indices) > 3) t_quad = .true.
        end do

        ! Both relate to the number of different direct excits, so count these
        directs = 0
        do ct = 1, step
            if(is_repeated_pair(indices, ct)) directs = directs + 1
        end do

        ! If there is a quadruple index, the spin of it is fixed (both have to be occupied), then,
        ! the prefactor is always +/- 2, depending on if there is more than one direct pair
        if(t_quad) then
            ! The only options with a quadruple index are
            ! a) one direct excitation => -2
            ! b) three direct exctiations => +2
            prefactor = merge(-2.0_dp, 2.0_dp, directs == 1)
        else
            ! In the other case, there are three relevant cases:
            ! a) All excitations are direct => prefactor 8 = 2**3 (all spins are free)
            ! b) One excitation is direct => prefactor -4 (from the spin freedom of the direct orb + overall spin)
            ! c) No excitation is direct => prefactor 4
            !    the factor of 4 comes from two sources: a factor of 2 from the overall spin
            !    and another factor of two from permutational symmetry: both indirect terms
            !    are identical in this case, so this entry has to be counted as both
            if(directs == 0) then
                prefactor = 4.0_dp
            elseif(directs == 1) then
                prefactor = -4.0_dp
            else ! directs == 3
                prefactor = 8.0_dp
            endif
        endif

    end subroutine frozen_diagonal_entry