pick_spin_par_elecs Subroutine

    subroutine pick_spin_par_elecs(nI, elecs, p_elec, opt_ispn)
        integer, intent(in) :: nI(nel)
        integer, intent(out) :: elecs(2)
        real(dp), intent(out) :: p_elec
        integer, intent(out), optional :: opt_ispn
#ifdef DEBUG_
        character(*), parameter :: this_routine = "pick_spin_par_elecs"
#endif
        integer :: ispn

        do
            elecs(1) = 1 + int(genrand_real2_dsfmt() * nel)

            do
                elecs(2) = 1 + int(genrand_real2_dsfmt() * nel)

                if (elecs(1) /= elecs(2)) exit

            end do
            ispn = get_ispn(nI(elecs))
            if (ispn /= 2) exit

        end do

        ASSERT(ispn == 1 .or. ispn == 3)

        ! and always output an ordered form
        elecs = [minval(elecs), maxval(elecs)]

        ! i have to rework the probability here..
        if (ispn == 1) then
            p_elec = 1.0_dp / real(nOccBeta * (nOccBeta - 1), dp)
        else if (ispn == 3) then
            p_elec = 1.0_dp / real(nOccAlpha * (nOccAlpha - 1), dp)
        end if

        ! in the special case of all spin-parallel:
        if (nOccBeta == 0 .or. nOccAlpha == 0) p_elec = 2.0_dp * p_elec

        if (present(opt_ispn)) opt_ispn = ispn

    end subroutine pick_spin_par_elecs