calc_orbital_pgen_contrib_start_def Subroutine

    subroutine calc_orbital_pgen_contrib_start_def(csf_i, occ_orbs, orb_a, orb_pgen)
        ! write a combined function for both r2l and l2r since its only
        ! one difference -> only one if condition to adjust for both!
        type(CSF_Info_t), intent(in) :: csf_i
        integer, intent(in) :: occ_orbs(2), orb_a
        real(dp), intent(out) :: orb_pgen

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

        ! electron indices
        i = gtID(occ_orbs(1))
        j = gtID(occ_orbs(2))

        ! damn... i need the probability of the elec-pair picking too or?

        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 mix up i and j here below.. what do i assume picked already
        ! here? if its really p(b|a) i should switch up i and j..
        ! hm.. the inputted j is the holeInd i which is fixed.. i gets looped
        ! over on the outside and is assumed picked first or 2nd here??
        ! taking i and j here is wrong! i is the open orbital, but j
        ! is the already picked electron! it has to be orb_a here or?
        call pgen_select_orb_guga_mol(csf_i, occ_orbs, i, orb_a, cpt_ba, ba_sum, i, .true.)

        ! change to the fullstart into fullstop: loop until orbital a
        if (tGen_guga_weighted) then
            do orb = i + 1, orb_a - 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 a also speciallly
        cpt_b = get_guga_integral_contrib(occ_orbs, orb_a, -1)

        cum_sum = cum_sum + cpt_b

        ! and get p(a|b)
        ! here the only difference between r2l and l2r fullstarts come into
        ! play!
        if (orb_a > j) then
            ! then orb_j is off-limits
            call pgen_select_orb_guga_mol(csf_i, occ_orbs, orb_a, i, cpt_ab, ab_sum, j)
        else
            ! in this case there is no restriction guga-wise..
            call pgen_select_orb_guga_mol(csf_i, occ_orbs, orb_a, i, cpt_ab, ab_sum)
        end if

        ! and deal with rest:
        if (tGen_guga_weighted) then
            do orb = orb_a + 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
            orb_pgen = 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

            ! and add them up to the final orbital pgen
            orb_pgen = cpt_a + cpt_b
        end if

    end subroutine calc_orbital_pgen_contrib_start_def