create_ab_list_hubbard Subroutine

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

        ! do the cum_arr for the k-space hubbard
        ! i think here i might really use flags.. and not just do the
        ! influence over the matrix elements.. since without transcorrelation
        ! i would waste alot of effort if i calculate the matrix elements
        ! here all the time..
        orb_list = -1
        cum_arr = 0.0_dp
        cum_sum = 0.0_dp

        ex(1, :) = src

        ! we are also using this routine for the parallel excitations due to
        ! the transcorrelation factor.. this is nice, since it is easily
        ! reusable, but, loses alot of efficiency, since the we are looping
        ! over all spin-orbital, although we know we only want to loop over a
        ! certain spin!
        ! todo
        if (present(tgt)) then
            ASSERT(present(cpt))

            cpt = 0.0_dp

            ! OPTIMIZATION: Do not loop over nbasis here, but over a pre-computed
            ! lookup table of excitations for src (if possible, is not an option
            ! if the matrix element depends on nI)
            do a = 1, nbasis
                elem = 0.0_dp
                ! if a is empty
                if (IsNotOcc(ilutI, a)) then
                    ! i have to rewrite get_orb, so it gives me the same
                    ! spin if src has the same spin! todo
                    ! to take into account spin-parallel double
                    ! excitations!

                    ! get the excitation
                    b = get_orb_from_kpoints(src(1), src(2), a)
                    ! we have yet to check if b is unoccupied
                    if (b /= a .and. IsNotOcc(ilutI, b)) then
                        ! assert that we hit opposite spins
                        if (.not. t_trans_corr_2body) then
                            ASSERT(.not. same_spin(a, b))
                        end if
                        ! get the matrix element (from storage)
                        if (.not. (t_trans_corr .or. t_trans_corr_2body)) then
                            elem = excit_cache(src(1), src(2), a)
                        else
                            ex(2, :) = [a, b]
                            call swap_excitations(nI, ex, nJ, ex2)
                            elem = abs(get_offdiag_helement_k_sp_hub(nJ, ex2, .false.))
                        end if
                    end if
                end if
                cum_sum = cum_sum + elem

                if (tgt == a) then
                    cpt = elem
                end if
            end do
            if (cum_sum < EPS) then
                cpt = 0.0_dp
            else
                ! todo: maybe i have to multiply by 2 here.. since both
                ! direction possible..
                cpt = cpt / cum_sum
            end if
        else
            do a = 1, nbasis
                elem = 0.0_dp
                b = -1

                if (IsNotOcc(ilutI, a)) then
                    b = get_orb_from_kpoints(src(1), src(2), a)

                    if (b /= a .and. IsNotOcc(ilutI, b)) then
                        ! get the matrix element (from storage)
                        if (.not. (t_trans_corr .or. t_trans_corr_2body)) then
                            elem = excit_cache(src(1), src(2), a)
                        else
                            ex(2, :) = [a, b]
                            call swap_excitations(nI, ex, nJ, ex2)
                            elem = abs(get_offdiag_helement_k_sp_hub(nJ, ex2, .false.))
                        end if
                    end if
                end if
                cum_sum = cum_sum + elem
                cum_arr(a) = cum_sum
                orb_list(a, :) = [a, b]
            end do
        end if

    end subroutine create_ab_list_hubbard