pgen_select_orb_guga_mol Subroutine

    subroutine pgen_select_orb_guga_mol(csf_i, occ_orbs, orb_b, orb_a, cpt, &
                                        cum_sum, orb_res, range_flag)
        ! routine to recalculate the pgen contribution if orbital (a) and (b)
        ! could have been picked in the opposite order
        ! additional GUGA-restrictions on orbitals are again dealt with
        ! optional input paramters
        type(CSF_Info_t), intent(in) :: csf_i
        integer, intent(in) :: occ_orbs(2), orb_b, orb_a
        real(dp), intent(out) :: cpt, cum_sum
        integer, intent(in), optional :: orb_res
        logical, intent(in), optional :: range_flag
        character(*), parameter :: this_routine = "pgen_select_orb_guga_mol"

        integer :: cc_a, label_index, nOrbs, i, orb
        real(dp) :: tmp
        logical :: tSingle

        ! first get the orbitals from the correct symmetry
        ! again have to settle on a certain spin -> if beta spin was occupied
        ! originally for a -> take alpha spin component
        ! otherwise choose beta
        ! UPDATE: change to only pick spatial orbitals -> only check if orbital
        ! was singly occupied
        if (csf_i%Occ_int(orb_b) == 1) then
            ! then i have to exclude orb_a in the recalculation of p(a|b) prob
            tSingle = .true.
        else
            tSingle = .false.
        end if

        ! need all allowed orbitals independent of actual "spin" of
        ! target orbita
        cc_a = ClassCountInd(2 * orb_a)

        label_index = SymLabelCounts2(1, cc_a)

        nOrbs = OrbClassCount(cc_a)

        cum_sum = 0.0_dp

        if (present(range_flag)) then
            ! if range flag is present, orb_b is already off-limits so no
            ! need to take tSingle into account
            ASSERT(present(orb_res))
            ASSERT(orb_res /= 0)

            if (orb_res < 0) then
                ! only orbitals below restriction
                do i = 1, nOrbs

                    orb = sym_label_list_spat(label_index + i - 1)

                    if (csf_i%stepvector(orb) /= 3 .and. orb < -orb_res) then

                        tmp = get_guga_integral_contrib(occ_orbs, orb_b, orb)

                        cum_sum = cum_sum + tmp
                        if (orb == orb_a) cpt = tmp
                    end if
                end do

            else
                ! only orbitals above restriction
                do i = 1, nOrbs
                    orb = sym_label_list_spat(label_index + i - 1)
                    if (csf_i%stepvector(orb) /= 3 .and. orb > orb_res) then
                        tmp = get_guga_integral_contrib(occ_orbs, orb_b, orb)

                        cum_sum = cum_sum + tmp
                        if (orb == orb_a) cpt = tmp
                    end if
                end do
            end if
        else
            if (present(orb_res)) then
                ASSERT(orb_res /= 0)
                ASSERT(orb_res > 0)
                ! the orbital associated with orb_res is off-limits

                ! also consider if orb b is singly occupied
                if (tSingle) then
                    do i = 1, nOrbs
                        orb = sym_label_list_spat(label_index + i - 1)

                        if (csf_i%stepvector(orb) /= 3 .and. orb /= orb_res &
                            .and. orb /= orb_b) then

                            tmp = get_guga_integral_contrib(occ_orbs, orb_b, orb)

                            cum_sum = cum_sum + tmp
                            if (orb == orb_a) cpt = tmp
                        end if
                    end do
                else
                    ! orb b is also allowed!
                    do i = 1, nOrbs
                        orb = sym_label_list_spat(label_index + i - 1)

                        if (csf_i%stepvector(orb) /= 3 .and. orb /= orb_res) then

                            tmp = get_guga_integral_contrib(occ_orbs, orb_b, orb)

                            cum_sum = cum_sum + tmp
                            if (orb == orb_a) cpt = tmp
                        end if
                    end do
                end if
            else
                ! no restrictions execpt single occupancy maybe
                if (tSingle) then
                    do i = 1, nOrbs
                        orb = sym_label_list_spat(label_index + i - 1)

                        if (csf_i%stepvector(orb) /= 3 .and. orb /= orb_b) then

                            tmp = get_guga_integral_contrib(occ_orbs, orb_b, orb)

                            cum_sum = cum_sum + tmp
                            if (orb == orb_a) cpt = tmp
                        end if
                    end do
                else
                    do i = 1, nOrbs
                        orb = sym_label_list_spat(label_index + i - 1)

                        if (csf_i%stepvector(orb) /= 3) then
                            tmp = get_guga_integral_contrib(occ_orbs, orb_b, orb)

                            cum_sum = cum_sum + tmp

                            if (orb == orb_a) cpt = tmp
                        end if
                    end do
                end if
            end if
        end if

    end subroutine pgen_select_orb_guga_mol