gen_double_back_spawn_hubbard Subroutine

    subroutine gen_double_back_spawn_hubbard(nI, ilutI, nJ, ilutJ, tParity, ExcitMat, &
                                             pgen, part_type)
        integer, intent(in) :: nI(nel)
        integer(n_int), intent(in) :: ilutI(0:niftot)
        integer, intent(out) :: nJ(nel), ExcitMat(2, maxExcit)
        integer(n_int), intent(out) :: ilutJ(0:niftot)
        logical, intent(out) :: tParity
        real(dp), intent(out) :: pgen
        integer, optional :: part_type

        integer :: elec_i, elec_j, iSpn, orb_b, src(2), elecs(2), &
                   loc, temp_part_type, orb_a
        real(dp) :: pAIJ, mult, pgen_elec
        logical :: t_temp_back_spawn

        ! damn.. remember if i initialize stuff above it implicitly assumes
        ! the (save) attribut!
        t_temp_back_spawn = .false.

        if (present(part_type)) then
            temp_part_type = part_type
        else
            temp_part_type = 1
        end if

        ! i know here that back-spawn is on and this is a non-inititator:
        ! so for the beginning implementation we pick the two electrons
        ! randomly
        ! i should make this to a function below: maybe with the hubbard
        ! flag as additional input to provide us with two independent
        ! electrons in any order possible!
        ! maybe for now.. i still want to retain the original
        ! back-spawn functionality
        if (t_back_spawn) then
            call pick_virtual_electrons_double_hubbard(nI, temp_part_type, elecs, src, ispn, &
                                                       pgen_elec)
            ! check if enogh electrons are in the virtual
            if (elecs(1) == 0) then
                nJ(1) = 0
                pgen = 0.0_dp
                return
            end if

            elec_i = elecs(1)
            elec_j = elecs(2)

            loc = -1
        else

            do
                elec_i = 1 + int(genrand_real2_dsfmt() * nel)
                do
                    elec_j = 1 + int(genrand_real2_dsfmt() * nel)
                    if (elec_j /= elec_i) exit
                end do
                ispn = get_ispn([nI(elec_i), nI(elec_j)])

                ! i need opposite spin-excitations in the hubbard model
                if (iSpn == 2) exit
            end do

            ! do i need 2* here?
            pgen_elec = 1.0_dp / real(nOccBeta * nOccAlpha, dp)

            src = nI([elec_i, elec_j])

            loc = check_electron_location(src, 2, temp_part_type)
        end if

        ! i need to call nI(elecs)
        ! this test should be fine since, if back_spawn is active loc should
        ! always be 0 so we are not risking of picking occupied orbitals
        ! below..

        ! wait a minute.. thats incorrect or?
        ! if we have both in the occupied manifold we want to restrict
        ! our orbital choice..
        ! i can decide based on the occ_virt_level or?
        if ((loc == 2) .or. (loc == 1 .and. occ_virt_level /= -1) .or. &
            (loc == 0 .and. occ_virt_level >= 1)) then
            t_temp_back_spawn = .true.
            ! i think i need a new one for the ueg also.. since there is not
            ! such a spin restriction as in the hubbard model
            ! i think just using the "normal" occupied picker should be fine
            ! how does this compile even??
            ! no.. write a new one without the spin-restriction
            call pick_occupied_orbital_hubbard(ilutI, temp_part_type, pAIJ, orb_a)
            ! i should take k-space symmetry into accound in picking
            ! orb a

            ! i guess we can have no possible excitations..
            if (orb_a == 0) then
                nJ(1) = 0
                pgen = 0.0_dp
                return
            end if
        else
            ! otherwise pick freely..
            ! in the hubbard case we always have iSpn = 2
            do
                orb_a = int(nBasis * genrand_real2_dsfmt()) + 1

                ! i guess it is more efficient to check if the momentum is
                ! correctly conserved here..

                if (IsNotOcc(ilutI, orb_a)) exit
            end do

            paij = 1.0_dp / real(nBasis - nel, dp)

        end if

        ! i guess this is not necessary, but anyway..
        IF ((orb_a < 0) .or. (orb_a > nBasis)) THEN
            CALL Stop_All("CreateDoubExcitLattice", "Incorrect basis function generated")
        end if

!         ! Is kb allowed by the size of the space?
!         ! Currently only applies when NMAXX etc. are set by the CELL keyword
!         tAllowedExcit = is_allowed_ueg_k_vector(src(1), src(2), orb_a)
!
!         IF(.not.tAllowedExcit) THEN
!             nJ(1)=0
!             RETURN
!         end if

        orb_b = get_orb_from_kpoints(src(1), src(2), orb_a)

        IF (orb_b <= 0 .or. orb_a == orb_b) THEN
            nJ(1) = 0
            pgen = 0.0_dp
            RETURN
        end if

        ! Is b occupied?
        if (IsOcc(ilutI, orb_b)) then
            nj(1) = 0
            pgen = 0.0_dp
            return
        end if

        ! Find the new determinant
        call make_double(nI, nJ, elec_i, elec_j, orb_a, &
                         orb_b, ExcitMat, tParity)

        ilutJ = make_ilutJ(ilutI, ExcitMat, 2)

        ! we knoe it is back-spawn + hubbard remember! so paij is already above
        !calculate generation probabilities
        ! note, p(b|ij)=p(a|ij) for this system
        ! if b is also in the occupied manifold in the case of back_spawn_flex
        mult = 2.0_dp
        if (t_temp_back_spawn .and. (.not. is_in_ref(orb_b, temp_part_type))) then
            mult = 1.0_dp
        end if

        pgen = mult * pgen_elec * pAIJ

    end subroutine gen_double_back_spawn_hubbard