gen_all_doubles_rs_hub_hop_transcorr Subroutine

    subroutine gen_all_doubles_rs_hub_hop_transcorr(nI, n_excits, det_list)
        integer, intent(in) :: nI(nel)
        integer, intent(out) :: n_excits
        integer(n_int), intent(out), allocatable :: det_list(:, :)
        integer(n_int) :: ilut(0:NIfTot), ilutJ(0:NIfTot)
        integer :: n_bound, i, src(2), j, ex(2, 2), pos, a, b
        integer(n_int), allocatable :: temp_list(:, :)
        real(dp) :: elem

        call EncodeBitDet(nI, ilut)

        n_excits = 1

        n_bound = nel * (nel - 1) * (nbasis - nel) * (nbasis - nel - 1)

        allocate(temp_list(0:NIfTot, n_bound))
        temp_list = 0_n_int

        ! we just have opposite spin doubles!
        do i = 1, nel - 1
            do j = i + 1, nel

                src = nI([i, j])

                if (same_spin(src(1), src(2))) cycle

                ex(1, :) = src
                ! and now loop over all possible empty spin opposite
                ! orbitals in a unique way..
                do a = 1, nbasis - 1
                    if (IsOcc(ilut, a)) cycle
                    do b = a + 1, nbasis
                        if (IsOcc(ilut, b) .or. same_spin(a, b)) cycle

                        ex(2, :) = [a, b]

                        elem = abs(get_helement_lattice(nI, 2, ex, .false.))
                        if (elem > EPS) then

                            ilutJ = make_ilutJ(ilut, ex, 2)
                            ! actually a search is not really necessary.. since
                            ! all the single excitations are unique.. but
                            ! just to be sure
                            pos = binary_search_ilut(temp_list(:, 1:n_excits), ilutJ, nifd + 1)

                            if (pos < 0) then

                                temp_list(:, n_excits) = ilutJ
                                call sort(temp_list(:, 1:n_excits), ilut_lt, ilut_gt)
                                n_excits = n_excits + 1
                                ! damn.. i have to sort everytime i guess..
                            end if
                        end if
                    end do
                end do
            end do
        end do

        n_excits = n_excits - 1
        allocate(det_list(0:NIfTot, n_excits), source=temp_list(:, 1:n_excits))

        call sort(det_list, ilut_lt, ilut_gt)

    end subroutine gen_all_doubles_rs_hub_hop_transcorr