pick_three_opp_elecs Subroutine

    subroutine pick_three_opp_elecs(nI, elecs, p_elec, opt_sum_ms)
        integer, intent(in) :: nI(nel)
        integer, intent(out) :: elecs(3)
        real(dp), intent(out) :: p_elec
        integer, intent(out), optional :: opt_sum_ms
#ifdef DEBUG_
        character(*), parameter :: this_routine = "pick_three_opp_elecs"
#endif
        integer :: sum_ms

        ! this can be done way more effective than this simple implementation:
        first: 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

            if (get_ispn(nI(elecs(1:2))) /= 2) then
                ! if the already picked electrons have same spin
                ! take opposite spin
                do
                    elecs(3) = 1 + int(genrand_real2_dsfmt() * nel)

                    if (elecs(1) /= elecs(3) .and. elecs(2) /= elecs(3) .and. &
                        .not. same_spin(nI(elecs(1)), nI(elecs(3)))) then
                        exit first
                    end if
                end do
            else
                ! here we can pick any spin electron

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

                    if (elecs(1) /= elecs(3) .and. elecs(2) /= elecs(3)) then
                        exit first
                    end if
                end do
            end if
        end do first

        ! i have to be careful how i could have gotten the electrons..
        ! because the order does not matter and there is no restriction set
        ! output an ordered form
        elecs = sort_unique(elecs)

        ! the first two electrons are picked totally randon, independent of
        ! spin just with the restriction i /= j
        ! so p(i) = 1/nel, p(j|i) = 1/(nel-1)
        ! the third electron has a spin restriction..
        ! if spin(i) = spin(j) then spin(k) must be -spin(i) giving
        ! p(k|ij) = 1/n_opp_spin
        ! if spin(i) /= spin(j) then k is freely except /= i,j
        ! and since the order is not important the total probability is
        ! p(i) * p(i|j) * [p(k|ij) + p(i|jk) + p(j|ik)]
        ! which translates to
        ! 1/nel * 1/(nel-1) * [1/n_opp + 2/(nel-2)]
        ! so depending on the total mz, we get

        sum_ms = sum(get_spin_pn(nI(elecs)))
        ASSERT(sum_ms == 1 .or. sum_ms == -1)
        if (sum_ms == 1) then
            ! then we have 2 alpha and 2 beta electrons
            p_elec = 2.0_dp / real(nel * (nel - 1), dp) * &
                     (1.0_dp / real(nOccBeta, dp) + 2.0_dp / real(nel - 2, dp))
        else if (sum_ms == -1) then
            ! then we have 2 beta electrons and 1 alpha
            p_elec = 2.0_dp / real(nel * (nel - 1), dp) * &
                     (1.0_dp / real(nOccAlpha, dp) + 2.0_dp / real(nel - 2, dp))
        end if

        ! adjust the edge cases
!         if (nOccAlpha == 1 .or. nOccBeta == 1) p_elec = 2.0_dp * p_elec

        if (present(opt_sum_ms)) opt_sum_ms = sum_ms

    end subroutine pick_three_opp_elecs