calc_orbital_pgen_contrib_end_def Subroutine

    subroutine calc_orbital_pgen_contrib_end_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))

        cum_sum = 0.0_dp
        if (tGen_guga_weighted) then
            do orb = 1, orb_a - 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 (a) in specific way:
        cpt_a = get_guga_integral_contrib(occ_orbs, orb_a, -1)

        cum_sum = cum_sum + cpt_a

        ! also get p(b|a)
        ! depending if its a r2l or l2r full-stop:
        if (i < orb_a) then
            ! its a L2R -> so no restrictions
            call pgen_select_orb_guga_mol(csf_i, occ_orbs, orb_a, j, cpt_ba, ba_sum)
        else
            ! its a R2L so orbital i is off-limits
            call pgen_select_orb_guga_mol(csf_i, occ_orbs, orb_a, j, cpt_ba, ba_sum, i)
        end if

        ! change to the fullstart into fullstop: loop until orbital j for the
        ! fullstop implementattion
        if (tGen_guga_weighted) then
            do orb = orb_a + 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)
        ! only orbitals below j are allowed!
        call pgen_select_orb_guga_mol(csf_i, occ_orbs, j, orb_a, 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
            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_end_def