subroutine gen_double_back_spawn_ueg(nI, ilutI, nJ, ilutJ, tParity, ExcitMat, &
pgen, part_type)
integer, intent(in) :: nI(nel)
integer(n_int), intent(in) :: ilutI(0:niftot)
integer, intent(out) :: nJ(nel), ExcitMat(2, maxExcit)
integer(n_int), intent(out) :: ilutJ(0:niftot)
logical, intent(out) :: tParity
real(dp), intent(out) :: pgen
integer, optional :: part_type
integer :: elec_i, elec_j, iSpn, orb_b, src(2), &
loc, temp_part_type, orb_a
real(dp) :: x, pAIJ, dummy, mult
logical :: tAllowedExcit, t_temp_back_spawn
t_temp_back_spawn = .false.
if (present(part_type)) then
temp_part_type = part_type
else
temp_part_type = 1
end if
! i know here that back-spawn is on and this is a non-inititator:
! so for the beginning implementation we pick the two electrons
! randomly
! i should make this to a function below: maybe with the hubbard
! flag as additional input to provide us with two independent
! electrons in any order possible!
! i do not need two do loops in the case of the UEG
elec_i = 1 + int(genrand_real2_dsfmt() * nel)
do
elec_j = 1 + int(genrand_real2_dsfmt() * nel)
if (elec_j /= elec_i) exit
end do
src = nI([elec_i, elec_j])
iSpn = get_ispn(src)
! i need to call nI(elecs)
loc = check_electron_location(src, 2, temp_part_type)
! wait a minute.. thats incorrect or?
! if we have both in the occupied manifold we want to restrict
! our orbital choice..
! i can decide based on the occ_virt_level or?
if ((loc == 2) .or. (loc == 1 .and. occ_virt_level /= -1) .or. &
(loc == 0 .and. occ_virt_level >= 1)) then
t_temp_back_spawn = .true.
! i think i need a new one for the ueg also.. since there is not
! such a spin restriction as in the hubbard model
! i think just using the "normal" occupied picker should be fine
! how does this compile even??
! no.. write a new one without the spin-restriction
call pick_occupied_orbital_ueg(ilutI, src, ispn, temp_part_type, pAIJ, &
dummy, orb_a)
if (orb_a == 0) then
nJ(1) = 0
pgen = 0.0_dp
return
end if
else
! otherwise pick freely..
do
x = genrand_real2_dSFMT()
if (iSpn == 2) then
orb_a = int(nBasis * x) + 1
pAIJ = 1.0_dp / real(nBasis - nel, dp)
else
orb_a = 2 * (INT(nBasis / 2 * x) + 1) & ! 2*(a number between 1 and nbasis/2) gives the alpha spin
& - (1 - (iSpn / 3)) ! alpha numbered even, iSpn/3 returns 1 for alpha/alpha, 0 for beta/beta
if (iSpn == 1) then
paij = 1.0_dp / (nbasis / 2 - noccbeta)
else !(ispn = 3)..
paij = 1.0_dp / (nbasis / 2 - noccalpha)
end if
end if
if (IsNotOcc(ilutI, orb_a)) exit
end do
end if
! i guess this is not necessary, but anyway..
IF ((orb_a < 0) .or. (orb_a > nBasis)) THEN
CALL Stop_All("CreateDoubExcitLattice", "Incorrect basis function generated")
end if
tAllowedExcit = is_allowed_ueg_k_vector(src(1), src(2), orb_a)
IF (.not. tAllowedExcit) THEN
nJ(1) = 0
RETURN
end if
orb_b = get_orb_from_kpoints(src(1), src(2), orb_a)
IF (orb_b == -1 .or. orb_a == orb_b) THEN
nJ(1) = 0
pgen = 0.0_dp
RETURN
end if
! Is b occupied?
if (.not. tAllowedExcit .or. IsOcc(ilutI, orb_b)) then
nj(1) = 0
pgen = 0.0_dp
return
end if
! Find the new determinant
call make_double(nI, nJ, elec_i, elec_j, orb_a, &
orb_b, ExcitMat, tParity)
ilutJ = make_ilutJ(ilutI, ExcitMat, 2)
! we knoe it is back-spawn + ueg remember! so paij is already above
!calculate generation probabilities
! note, p(b|ij)=p(a|ij) for this system
! i have to be careful.. only if b is also in the occupied manifold
! if we forced the pick.. then it is * 2
mult = 2.0_dp
if (t_temp_back_spawn .and. (.not. is_in_ref(orb_b, part_type))) then
mult = 1.0_dp
end if
pgen = 2.0_dp * mult * paij / real(nel * (nel - 1), dp)
end subroutine gen_double_back_spawn_ueg