# pick_three_opp_elecs Subroutine

## public subroutine pick_three_opp_elecs(nI, elecs, p_elec, opt_sum_ms)

### Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nI(nel)
integer, intent(out) :: elecs(3)
real(kind=dp), intent(out) :: p_elec
integer, intent(out), optional :: opt_sum_ms

## Source Code

    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