create_bc_list_hubbard Subroutine

    subroutine create_bc_list_hubbard(nI, ilutI, src, orb_a, orb_list, cum_arr, &
                                      cum_sum, tgt, cpt)
        integer, intent(in) :: nI(nel), src(3), orb_a
        integer(n_int), intent(in) :: ilutI(0:niftot)
        integer, intent(out) :: orb_list(nbasis / 2, 2)
        real(dp), intent(out) :: cum_arr(nbasis / 2), cum_sum
        integer, intent(in), optional :: tgt
        real(dp), intent(out), optional :: cpt
#ifdef DEBUG_
        character(*), parameter :: this_routine = "create_bc_list_hubbard"
#endif
        integer :: b, c, ex(2, 3), spin, orb_b
        real(dp) :: elem
        integer :: nJ(nel)
        integer, allocatable :: ex2(:, :)

        orb_list = -1
        cum_arr = 0.0_dp
        cum_sum = 0.0_dp
        ex(1, :) = src
        ex(2, 1) = orb_a

        ! we want to do a restriction! to make it easier to recalculate the
        ! pgens and stuff!
        ! the restriction is, that the first picked orbital (a) is of the
        ! minority spin of the picked electrons.
        ! so if we have to alpha and 1 beta electron picked, orbital (a) is
        ! a beta orbital and the last 2 orbitals will be alpha and v.v.

        ! decide that the  first orbital orb_a, which is already picked is the
        ! minority spin electron, so now pick two electrons of the opposite
        ! spin
        if (is_beta(orb_a)) then
            ! then we want alpha orbitals
            spin = 0
            ! and also be sure that we did the right thing until now,
            ! otherwise it breaks
            ASSERT(sum(get_spin_pn(src)) == 1)
        else
            ! otherwise we want beta now
            spin = 1
            ASSERT(sum(get_spin_pn(src)) == -1)
        end if

        if (present(tgt)) then
            ASSERT(present(cpt))

            cpt = 0.0_dp

            ! does the spin of tgt fit?
            ! it must be opposite of orb_a!
            if (same_spin(orb_a, tgt)) return

            !TODO: we only need to consider one spin-type!!
            do b = 1, nbasis / 2
                elem = 0.0_dp

                ! convert to the appropriate spin-orbitals
                orb_b = 2 * b - spin

                if (IsNotOcc(ilutI, orb_b)) then
                    ! get the appropriate momentum conserverd orbital c
                    c = get_orb_from_kpoints_three(src, orb_a, orb_b)

                    if (c /= orb_b .and. IsNotOcc(ilutI, c)) then

                        ex(2, 2:3) = [orb_b, c]

                        ! actually i messed up with the non-hermiticity
                        ! i should actually switch the order of the
                        ! determinants in matrix element calculation
                        call swap_excitations(nI, ex, nJ, ex2)
                        elem = abs(get_3_body_helement_ks_hub(ex2, .false.))

                    end if
                end if
                cum_sum = cum_sum + elem
                if (tgt == orb_b) then
                    cpt = elem
                end if
            end do

            if (cum_sum < EPS) then
                cpt = 0.0_dp
            else
                cpt = cpt / cum_sum
            end if
        else
            do b = 1, nbasis / 2
                orb_b = 2 * b - spin

                elem = 0.0_dp
                c = 0

                if (IsNotOcc(ilutI, orb_b)) then
                    c = get_orb_from_kpoints_three(src, orb_a, orb_b)

                    if (c /= orb_b .and. IsNotOcc(ilutI, c)) then

                        ex(2, 2:3) = [orb_b, c]
                        call swap_excitations(nI, ex, nJ, ex2)
                        elem = abs(get_3_body_helement_ks_hub(ex2, .false.))
                    end if
                end if
                cum_sum = cum_sum + elem
                cum_arr(b) = cum_sum
                orb_list(b, :) = [orb_b, c]
            end do
        end if

    end subroutine create_bc_list_hubbard