calc_orbital_pgen_contr_mol Subroutine

    subroutine calc_orbital_pgen_contr_mol(csf_i, occ_orbs, cpt_a, cpt_b)
        ! calculates the cumulatice probability list for different
        ! full-start -> full-stop mixed excitations, used in the recalculation of
        ! contrbuting pgens from different picked orbitals
        type(CSF_Info_t), intent(in) :: csf_i
        integer, intent(in) :: occ_orbs(2)
        real(dp), intent(out) :: cpt_a, cpt_b

        integer :: i, j, orb
        real(dp) :: cum_sum, cpt_ab, cpt_ba, ba_sum, ab_sum

        ! given 2 already picked electrons, this routine creates a list of
        ! p(a)*p(b|a) probabilities to pick the already determined holes
        ! is this so easy??
        ! what do i know? i know there is a possible switch between i and j
        ! or otherwise those indices would not have been taken, since an
        ! excitation was already created
        ! i know electrons i and j are from different orbitals and singly
        ! occupied
        ! have to loop once over all orbitals to create list for p(x)
        ! most def, but that should be easy since there are not so many
        ! restrictions, but to determine p(b|a) i have to take into account
        ! the symmetries and the guga restrictions..
        ! but i need the whole information to get cum_sum correct, for the
        ! given electrons i and j..
        ! maybe split up the loop to different sectors, where i know more info

        ! UPDATE: calculate it more directly! also for UEG models!
        ! just output the nececcary p(a)*p(b|a) and vv. and not a whole
        ! list!

        ! chanve that routine now, to use the general pre-generated cum-list
        ! for the current CSF
        i = gtID(minval(occ_orbs))
        j = gtID(maxval(occ_orbs))

        cum_sum = 0.0_dp

        if (tGen_guga_weighted) then
            do orb = 1, i - 1
                ! calc. the p(a)
                if (csf_i%stepvector(orb) /= 3) then
                    cum_sum = cum_sum + get_guga_integral_contrib(occ_orbs, orb, -1)

                end if

            end do
        end if

        ! deal with orb (i) in specific way:
        cpt_a = get_guga_integral_contrib(occ_orbs, i, -1)

        cum_sum = cum_sum + cpt_a

        ! also get p(b|a)
        ! did i get that the wrong way around??
        call pgen_select_orb_guga_mol(csf_i, occ_orbs, i, j, cpt_ba, ba_sum, i, .true.)

        if (tGen_guga_weighted) then
            do orb = i + 1, j - 1
                if (csf_i%stepvector(orb) /= 3) then
                    cum_sum = cum_sum + get_guga_integral_contrib(occ_orbs, orb, -1)
                end if
            end do
        end if

        ! deal with j also speciallly
        cpt_b = get_guga_integral_contrib(occ_orbs, j, -1)

        cum_sum = cum_sum + cpt_b

        ! and get p(a|b)
        call pgen_select_orb_guga_mol(csf_i, occ_orbs, j, i, cpt_ab, ab_sum, -j, .true.)

        ! and deal with rest:

        if (tGen_guga_weighted) then
            do orb = j + 1, nSpatOrbs
                if (csf_i%stepvector(orb) /= 3) then
                    cum_sum = cum_sum + get_guga_integral_contrib(occ_orbs, orb, -1)
                end if
            end do
        end if

        if (.not. tGen_guga_weighted) then
            cum_sum = csf_i%cum_list(nSpatOrbs)
        end if

        if (near_zero(cum_sum) .or. near_zero(ab_sum) .or. near_zero(ba_sum)) then
            cpt_a = 0.0_dp
            cpt_b = 0.0_dp
        else
            ! and get hopefully correct final values:
            cpt_a = cpt_a / cum_sum * cpt_ba / ba_sum
            cpt_b = cpt_b / cum_sum * cpt_ab / ab_sum
        end if

    end subroutine calc_orbital_pgen_contr_mol