create_ab_list_par_hubbard Subroutine

    subroutine create_ab_list_par_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, 2)
        real(dp), intent(out) :: cum_arr(nbasis / 2)
        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_par_hubbard"
#endif
        integer :: a, b, ex(2, 2), spin, orb_a
        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

        ! this routine only checks for parallel spins, depending on src
        ASSERT(same_spin(src(1), src(2)))

        ! make a spin factor for the orbital conversion
        ! 0...alpha
        ! 1...beta
        spin = get_spin(src(1)) - 1

        ! and only loop over the correct spin
        if (present(tgt)) then
            ASSERT(present(cpt))

            cpt = 0.0_dp

            ! if target does not have the same spin, do we abort or return?
            if (.not. same_spin(src(1), tgt)) return

            do a = 1, nbasis / 2
                elem = 0.0_dp

                orb_a = 2 * a - spin

                if (IsNotOcc(ilutI, orb_a)) then
                    ! modify get_orb_from_kpoints to give spin-parallel
                    ! it already does i think!
                    b = get_orb_from_kpoints(src(1), src(2), orb_a)

                    if (b /= orb_a .and. IsNotOcc(ilutI, b)) then

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

                        call swap_excitations(nI, ex, nJ, ex2)
                        elem = abs(get_offdiag_helement_k_sp_hub(nJ, ex2, .false.))

                    end if
                end if
                cum_sum = cum_sum + elem
                if (tgt == orb_a) 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 a = 1, nbasis / 2
                orb_a = 2 * a - spin

                elem = 0.0_dp
                b = 0

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

                    if (b /= orb_a .and. IsNotOcc(ilutI, b)) then
                        ex(2, :) = [orb_a, b]
                        call swap_excitations(nI, ex, nJ, ex2)
                        elem = abs(get_offdiag_helement_k_sp_hub(nJ, ex2, .false.))
                    end if
                end if
                cum_sum = cum_sum + elem
                cum_arr(a) = cum_sum
                orb_list(a, :) = [orb_a, b]
            end do
        end if

    end subroutine create_ab_list_par_hubbard