get_missing_elec Subroutine

    subroutine get_missing_elec(nI, elecs, nOcc, nPicked, tAlpha, pgen)
        ! after picking two electrons using the symrandexcit routines,
        ! get a third one with the missing spin
        integer, intent(in) :: nI(nel)
        integer, intent(inout) :: elecs(3) ! picked electrons, elecs(1:2) has to be assigned on entry
        integer, intent(in) :: nOcc    ! number of occupied orbs of the type of the missing one
        integer, intent(in) :: nPicked ! number of already picked elecs
        logical, intent(in) :: tAlpha  ! if the missing electron is alpha
        real(dp), intent(inout) :: pgen ! probability of generation
        real(dp) :: r
        integer :: index, i, count

        elecs(3) = 0
        if (nOcc <= nPicked) then
            ! there are not enought elecs available - invalid excitation
            return
        end if

        r = genrand_real2_dSFMT()
        ! pick a random index of an electron
        index = int((nOcc - nPicked) * r) + 1

        ! get the corresponding electron
        count = 0
        do i = 1, nel
            if (tAlpha .neqv. is_beta(nI(i))) then
                count = count + 1
                if (count == index) elecs(3) = i
            end if
        end do
        ! the picked electrons are not counted, so skip their indices

        ! if we need, skip the first index
        call skipElec(1)
        ! if we need, skip the second index
        call skipElec(2)
        ! note that elecs(2) > elecs(1), so we cannot end up on elecs(1) at the end of
        ! skipping

        ! uniformly chosen
        pgen = pgen / (nOcc - nPicked)

    contains

        subroutine skipElec(ind)
            integer, intent(in) :: ind

            ! if we need to skip an index
            if ((tAlpha .neqv. is_beta(nI(elecs(ind))))) then
                ! if we are above the index, we need to add 1 more, because we did not
                ! take elecs(ind) into account when picking elecs(3)
                if (elecs(3) >= elecs(ind)) then
                    ! jump to the next electron with the right spin
                    elecs(3) = elecs(3) + 1
                    do while (.not. (tAlpha .neqv. is_beta(nI(elecs(3)))))
                        elecs(3) = elecs(3) + 1
                    end do
                end if
            end if
        end subroutine skipElec
    end subroutine get_missing_elec