subroutine pick_virtual_electrons_double_hubbard(nI, part_type, elecs, src, &
ispn, pgen, calc_pgen)
! specific routine to pick 2 electrons in the k-space hubbard,
! since apparently it is important to allow all orderings of
! electrons possible.. although this could just be a artifact of the
! old hubbard excitation generation
integer, intent(in) :: nI(nel), part_type
integer, intent(out) :: elecs(2), src(2), ispn
real(dp), intent(out) :: pgen
logical, intent(in), optional :: calc_pgen
character(*), parameter :: this_routine = "pick_virtual_electrons_double_hubbard"
integer :: n_valid, i, j, ind_1, ind_2
integer :: virt_elecs(nel)
integer :: n_beta, n_alpha
! but it is also good to to it here so i can do it more cleanly
n_valid = 0
n_beta = 0
n_alpha = 0
! actually for the correct generation probabilities i have to count
! the number of valid alpha and beta electrons!
virt_elecs = -1
elecs = -1
src = -1
j = 1
do i = 1, nel
if (is_in_virt_mask(nI(i), part_type)) then
if (is_beta(nI(i))) n_beta = n_beta + 1
if (is_alpha(nI(i))) n_alpha = n_alpha + 1
! the electron is in the virtual of the
n_valid = n_valid + 1
virt_elecs(j) = i
j = j + 1
end if
end do
! in the hubbard case i also have to check if there is atleast on
! pair possible with opposite spin
if (n_valid < 2 .or. n_beta == 0 .or. n_alpha == 0) 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
! what is ispn on return?? argh too many uninitialized vars.
ispn = 0
return
end if
! apparently i have to have both ordering of the electrons in
! the hubbard excitation generator
! but it must be easier to do that... and more efficient
ASSERT(n_valid > 1)
pgen = 1.0_dp / real(n_alpha * n_beta, dp)
if (present(calc_pgen)) then
if (calc_pgen) return
end if
do i = 1, 1000
ind_1 = 1 + int(n_valid * genrand_real2_dSFMT())
do j = 1, 1000
ind_2 = 1 + int(n_valid * genrand_real2_dSFMT())
if (ind_1 /= ind_2) exit
end do
elecs(1) = virt_elecs(ind_1)
elecs(2) = virt_elecs(ind_2)
src = nI(elecs)
if (is_beta(src(1)) .neqv. is_beta(src(2))) then
ispn = 2
exit
end if
end do
if (i > 999 .or. j > 999) then
print *, "something went wrong, did not find two valid electrons!"
print *, "nI: ", nI
print *, "mask_virt_ni:", mask_virt_ni(:, part_type_to_run(part_type))
print *, "virt_elecs: ", virt_elecs
end if
end subroutine pick_virtual_electrons_double_hubbard