pick_virtual_electrons_double Subroutine

    subroutine pick_virtual_electrons_double(nI, part_type, elecs, src, ispn, &
                                             sum_ml, pgen, calc_pgen)
        ! this is the important routine!
        ! for non-initiator determinants this pick electrons only from the
        ! virtual orbitals of the reference determinant to increase the
        ! chance to de-excite and to spawn to an already occupied
        ! determinant from an non-initiator!
        integer, intent(in) :: nI(nel), part_type
        integer, intent(out) :: elecs(2), src(2), ispn, sum_ml
        real(dp), intent(out) :: pgen
        logical, intent(in), optional :: calc_pgen
        character(*), parameter :: this_routine = "pick_virtual_electrons_double"

        integer :: i, n_valid, j, ind, n_valid_pairs, ind_1, ind_2
        integer :: virt_elecs(nel)

        ! i guess for now i only want to choose uniformly from all the
        ! available electron in the virtual orbitals of the reference

        ! what do we need here?
        ! count all the electrons in the virtual of the reference, then
        ! pick two random orbitals out of those!
        ! check the routine in symrandexcit3.f90 this does the job i guess..

        n_valid = 0
        j = 1
        virt_elecs = -1
        elecs = -1
        src = -1
        ispn = -1
        sum_ml = -1

        do i = 1, nel
            if (is_in_virt_mask(nI(i), part_type)) then
                n_valid = n_valid + 1
                virt_elecs(j) = i
                j = j + 1
            end if
        end do

        if (n_valid < 2) then
            ! something went wrong
            ! in this case i have to abort as no valid double excitation
            ! could have been found
            elecs = 0
            src = 0
            pgen = 0.0_dp
            ! can i set defaults here which do not break the rest of the
            ! code?..
            iSpn = -1
            sum_ml = -1
            return
        end if

        ! determine how many valid pairs there are now
        n_valid_pairs = (n_valid * (n_valid - 1)) / 2

        ASSERT(n_valid_pairs > 0)
        ! and the pgen is now:
        pgen = 1.0_dp / real(n_valid_pairs, dp)

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

        ! and is it now enough to do is just like in the symrandexcit3 routine:
        ind = 1 + int(n_valid_pairs * genrand_real2_dSFMT())
        ind_1 = ceiling((1 + sqrt(1 + 8 * real(ind, dp))) / 2)
        ind_2 = ind - ((ind_1 - 1) * (ind_1 - 2)) / 2

        ! and retro pick the electron number from the created list?
        elecs(1) = virt_elecs(ind_1)
        elecs(2) = virt_elecs(ind_2)

        ! hm.. test this tomorrow

        ! now i have to pick two random ones from the list!
        ! all the symmetry related stuff at the end:
        src = nI(elecs)

        ispn = get_ispn(src)

        ! The Ml value is obtained from the orbitals
        sum_ml = sum(G1(src)%Ml)

    end subroutine pick_virtual_electrons_double