create_cum_list_rs_hubbard_transcorr_double Subroutine

    subroutine create_cum_list_rs_hubbard_transcorr_double(nI, ilutI, src, orb_a, &
                                                           cum_arr, cum_sum, tgt, p_orb)
        ! routine to create the cum-list to pick the second orbital given
        ! the electrons (src) and the first orbital (orb_a)
        ! if second orbital (tgt) is given it calculates the probability (p_orb)
        ! to have picked this orbital.
        integer, intent(in) :: nI(nel), src(2), orb_a
        integer(n_int), intent(in) :: ilutI(0:NIfTot)
        real(dp), intent(out) :: cum_arr(nBasis / 2), cum_sum
        integer, intent(in), optional :: tgt
        real(dp), intent(out), optional :: p_orb
#ifdef DEBUG_
        character(*), parameter :: this_routine = "create_cum_list_rs_hubbard_transcorr_double"
#endif
        integer :: ex(2, 2), spin, b, nJ(nel), orb_b
        integer, allocatable :: ex2(:, :)
        real(dp) :: elem

        ASSERT(.not. same_spin(src(1), src(2)))
        ASSERT(IsNotOcc(ilutI, orb_a))
        ASSERT(IsOcc(ilutI, src(1)))
        ASSERT(IsOcc(ilutI, src(2)))

        cum_arr = 0.0_dp
        cum_sum = 0.0_dp

        ! make a spin factor for the orbital conversion
        ! 1...alpha
        ! 2...beta
        spin = get_spin(orb_a)

        ex(1, :) = src
        ex(2, 1) = orb_a

        if (present(tgt)) then
            ASSERT(present(p_orb))
            ASSERT(.not. same_spin(orb_a, tgt))

            p_orb = 0.0_dp

            do b = 0, nBasis / 2 - 1
                elem = 0.0_dp

                ! add the spin to get correct anti-parallel spin-orbtial to (a)
                ! if (a) is alpha, spin = 1 -> so add
                orb_b = 2 * b + spin

                if (IsNotOcc(ilutI, orb_b)) then
                    ! with an occupancy everything is fine.. since a == b
                    ! is not possible due to opposite spin

                    ex(2, 2) = orb_b
                    call swap_excitations(nI, ex, nJ, ex2)
                    ! elem = abs(get_double_helem_rs_hub_transcorr(ex, .false.))
                    elem = abs(get_double_helem_rs_hub_transcorr(ex2, .false.))

                end if
                cum_sum = cum_sum + elem

                if (tgt == orb_b) then
                    p_orb = elem
                end if
            end do
            if (cum_sum < EPS) then
                p_orb = 0.0_dp
            else
                p_orb = p_orb / cum_sum
            end if
        else
            do b = 0, nBasis / 2 - 1

                elem = 0.0_dp
                orb_b = 2 * b + spin

                if (IsNotOcc(ilutI, orb_b)) then

                    ex(2, 2) = orb_b
                    call swap_excitations(nI, ex, nJ, ex2)
                    ! elem = abs(get_double_helem_rs_hub_transcorr(ex, .false.))
                    elem = abs(get_double_helem_rs_hub_transcorr(ex2, .false.))

                end if

                cum_sum = cum_sum + elem
                cum_arr(b + 1) = cum_sum
            end do
        end if

    end subroutine create_cum_list_rs_hubbard_transcorr_double