pick_b_orb_guga_mol Subroutine

    subroutine pick_b_orb_guga_mol(csf_i, occ_orbs, orb_a, cc_b, int_contrib, &
                                   cum_sum, orb_b, orb_res, range_flag)
        ! restrict the b, if orbital (a) is singly occupied already..
        ! and switch to spatial orbital picking!
        type(CSF_Info_t), intent(in) :: csf_i
        integer, intent(in) :: occ_orbs(2), orb_a, cc_b
        real(dp), intent(out) :: int_contrib, cum_sum
        integer, intent(out) :: orb_b
        integer, intent(in), optional :: orb_res
        logical, intent(in), optional :: range_flag
        character(*), parameter :: this_routine = "pick_b_orb_guga_mol"

        real(dp) :: cum_arr(OrbClassCount(cc_b)), r
        integer :: nOrbs, label_index, i, orb, orb_index
        logical :: tSingle

        ! rewrite this funcitonality new:
        nOrbs = OrbClassCount(cc_b)

        label_index = SymLabelCounts2(1, cc_b)

        cum_sum = 0.0_dp

        ! have to predetermine is already picked orbital is singly occupied
        ! already: if yes, its not allowd to be picked again in here
        if (csf_i%Occ_int(orb_a) == 1) then
            tSingle = .true.
        else
            tSingle = .false.
        end if

        ! make a spatial orbital list of symmetry allowed orbitals!
        ! but that probably should be done globally in the setup phase, since
        ! this does not change
        ! did that in sym_label_list_spat

        ! depending on input do the specific loops
        if (present(range_flag)) then
            if (range_flag) then
                ASSERT(present(orb_res))
                ASSERT(orb_res /= 0)

                ! the case that a whole range is forbidden due to GUGA restrictions
                ! is only in the case, when orbital (a) is on of the original I,J
                ! orbitals, so in this case a will not be chosen due to the
                ! GUGA restrictions already

                ! the direction of the range is indicated through the sign of the
                ! index restriciton
                if (orb_res < 0) then
                    ! only orbitals below orb_res 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
                            cum_sum = cum_sum + &
                                      get_guga_integral_contrib(occ_orbs, orb_a, orb)
                        end if

                        cum_arr(i) = cum_sum
                    end do
                else
                    ! only orbitals above orb_res are 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
                            cum_sum = cum_sum + &
                                      get_guga_integral_contrib(occ_orbs, orb_a, orb)
                        end if

                        cum_arr(i) = cum_sum
                    end do
                end if
            end if
        else
            if (present(orb_res)) then
                ! should i assert here, that orb_res should not be 0?
                ! otherwise it would be stupid to input..
                ASSERT(orb_res /= 0)
                ASSERT(orb_res > 0)
                ! but now i have to include that orb (a) might be off-limits!
                if (tSingle) then
                    ! then orb_a is also off-limits!
                    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_a) then
                            cum_sum = cum_sum + &
                                      get_guga_integral_contrib(occ_orbs, orb_a, orb)
                        end if

                        cum_arr(i) = cum_sum
                    end do
                else
                    ! orb a is not off-limits
                    do i = 1, nOrbs
                        orb = sym_label_list_spat(label_index + i - 1)

                        if (csf_i%stepvector(orb) /= 3 .and. orb /= orb_res) then
                            cum_sum = cum_sum + &
                                      get_guga_integral_contrib(occ_orbs, orb_a, orb)
                        end if

                        cum_arr(i) = cum_sum

                    end do
                end if
            else
                ! no guga restrictions, only have to check if orb a was single
                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_a) then
                            cum_sum = cum_sum + &
                                      get_guga_integral_contrib(occ_orbs, orb_a, orb)
                        end if

                        cum_arr(i) = cum_sum
                    end do
                else
                    ! no restrictions except double occuations
                    do i = 1, nOrbs
                        orb = sym_label_list_spat(label_index + i - 1)

                        if (csf_i%stepvector(orb) /= 3) then
                            cum_sum = cum_sum + &
                                      get_guga_integral_contrib(occ_orbs, orb_a, orb)
                        end if
                        cum_arr(i) = cum_sum
                    end do
                end if
            end if
        end if

        if (near_zero(cum_sum)) then
            orb_b = 0
            return
        end if

        r = genrand_real2_dSFMT() * cum_sum
        orb_index = binary_search_first_ge(cum_arr, r)

        orb_b = sym_label_list_spat(label_index + orb_index - 1)

        if (orb_index == 1) then
            int_contrib = cum_arr(orb_index)
        else
            int_contrib = cum_arr(orb_index) - cum_arr(orb_index - 1)
        end if

    end subroutine pick_b_orb_guga_mol