pick_second_occupied_orbital Subroutine

    subroutine pick_second_occupied_orbital(ilutI, cc_b, orb_a, ispn, &
                                            part_type, cpt, cum_sum, orb, calc_pgen)
        ! routine which picks second orbital from the occupied manifold for
        ! a double excitation. this function gets called if we have picked
        ! two electrons also from the occupied manifold in the flex version
        ! of the back-spawning method. to ensure we keep the excitation
        ! level the same but also increase the flexibility of the method
        ! this now has to take symmetries into account, which makes it a
        ! bit more complicated
        integer, intent(in) :: cc_b, orb_a, ispn, part_type
        integer(n_int), intent(in) :: ilutI(0:niftot)
        real(dp), intent(out) :: cpt, cum_sum
        integer, intent(out) :: orb
        logical, intent(in), optional :: calc_pgen
        character(*), parameter :: this_routine = "pick_second_occupied_orbital"

        integer :: label_index, norb, sym_orbs(OrbClassCount(cc_b))
        integer :: i, n_valid, occ_orbs(nel), j, ind, orb_b
        ! i need to take symmetry and spin into account for the valid
        ! "occupied" orbitals.
        ! because we have picked the first indepenent of spin and symmetry

        ! i could compare the reference det and the symmetry allowed list
        ! of orbitals
        label_index = SymLabelCounts2(1, cc_b)
        norb = OrbClassCount(cc_b)

        ! create the symmetry allowed orbital list
        do i = 1, norb
            sym_orbs(i) = SymLabelList2(label_index + i - 1)
        end do

        j = 1
        occ_orbs = -1
        orb = -1

        ! check which occupied orbitals fit all the restrictions:
        ! or i guess this is already covered in the symlabel list!
        ! check that!
        if (ispn == 2) then
            ! then we want the opposite spin of orb_a!
            do i = 1, nel
                orb_b = projedet(i, part_type_to_run(part_type))
                ! check if in occupied manifold
                if (IsNotOcc(ilutI, orb_b)) then
                    ! check if symmetry fits
                    if (any(orb_b == sym_orbs)) then
                        ! and check if spin is opposit
!                         if (.not. (is_beta(orb_a) .eqv. is_beta(orb_b))) then
                        if (.not. same_spin(orb_a, orb_b)) then
                            occ_orbs(j) = orb_b
                            j = j + 1
                        end if
                    end if
                end if
            end do
        else
            ! otherwise we want the same spin but have to ensure it is not
            ! already picked orbital (a)
            do i = 1, nel
                orb_b = projedet(i, part_type_to_run(part_type))
                if (IsNotOcc(ilutI, orb_b)) then
                    if (any(orb_b == sym_orbs)) then
                        if (same_spin(orb_a, orb_b) .and. (orb_a /= orb_b)) then
                            occ_orbs(j) = orb_b
                            j = j + 1
                        end if
                    end if
                end if
            end do
        end if

        n_valid = j - 1

        if (n_valid == 0) then
            orb = 0
            cpt = 0.0_dp
            ! ok here i decide to output it to 1.0.. so do it in all the
            ! other routines too..
            cum_sum = 1.0_dp
            return
        end if

        ASSERT(n_valid > 0)
        cpt = 1.0_dp / real(n_valid, dp)
        cum_sum = 1.0_dp

        if (present(calc_pgen)) then
            if (calc_pgen) return
        end if

        ind = 1 + int(genrand_real2_dSFMT() * n_valid)

        orb = occ_orbs(ind)

        ASSERT(orb > 0)
        ASSERT(orb <= nbasis)

    end subroutine pick_second_occupied_orbital