get_bit_excitmat Subroutine

    pure subroutine get_bit_excitmat(ilutI, iLutJ, ex, IC)

        ! Obatin the excitation matrix between two determinants from their bit
        ! representation without calculating tSign --> a bit quicker.
        !
        ! In:    iLutI, iLutJ - Bit representations of determinants I,J
        ! InOut: IC           - Specify max IC before bailing, and return
        !                       number of orbital I,J differ by
        ! Out:   ex           - Excitation matrix between I,J

        integer(n_int), intent(in) :: iLutI(0:NIfD), iLutJ(0:NIfD)
        integer, intent(inout)  :: IC
        integer, intent(out) :: ex(2, IC)

        integer(n_int) :: ilut(0:NIfD, 2)
        integer :: pos(2), max_ic, i, j, k

        ! Obtain bit representations of I,J containing only unique orbitals
        ilut(:, 1) = ieor(ilutI, ilutJ)
        ilut(:, 2) = iand(ilutJ, ilut(:, 1))
        ilut(:, 1) = iand(ilutI, ilut(:, 1))

        max_ic = IC
        pos = 0
        IC = 0
        do i = 0, NIfD
            do j = 0, bits_n_int - 1
                do k = 1, 2
                    if (pos(k) < max_ic) then
                        if (btest(ilut(i, k), j)) then
                            pos(k) = pos(k) + 1
                            IC = max(IC, pos(k))
                            ex(k, pos(k)) = bits_n_int * i + j + 1
                        end if
                    end if
                end do
                if (pos(1) >= max_ic .and. pos(2) >= max_ic) return
            end do
        end do

    end subroutine get_bit_excitmat