# pick_three_orbs_ueg Subroutine

## public subroutine pick_three_orbs_ueg(nI, src, tgt, pgen, ms)

### Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nI(nel)
integer, intent(in) :: src(3)
integer, intent(out) :: tgt(3)
real(kind=dp), intent(inout) :: pgen
integer, intent(in) :: ms

## Source Code

subroutine pick_three_orbs_ueg(nI, src, tgt, pgen, ms)
use SymExcitDataMod, only: kPointToBasisFn
use SystemData, only: nBasisMax, tOrbECutoff, OrbECutoff, nmaxx, nmaxy, nmaxz
use sym_mod, only: mompbcsym

! picks three random unoccupied orbitals, given the occupied orbitals
integer, intent(in) :: nI(nel), ms, src(3)
integer, intent(out) :: tgt(3)
real(dp), intent(inout) :: pgen

integer :: i, msCur, msOrb, k1(3), k2(3), k3(3), p1(3), p2(3), p3(3), j
integer :: k, l
real(dp) :: testE
logical :: i_occup, is_allowed

msCur = ms
do i = 0, 1
! get the ms of this orb
! we take ms = 1 until the total leftover ms is negative
if (msCur > 0) then
msOrb = 1
else
! then we start taking ms = -1
msOrb = -1
end if
call get_rand_orb(nI, tgt, msOrb, i, pgen)
! the remaining ms
msCur = msCur - msOrb
end do

if (msCur > 0) then
msOrb = 1
else
! then we start taking ms = -1
msOrb = -1
end if

do i = 1, 3
k1(i) = G1(src(1))%k(i)
k2(i) = G1(src(2))%k(i)
k3(i) = G1(src(3))%k(i)
end do

do i = 1, 3
p1(i) = G1(tgt(1))%k(i)
p2(i) = G1(tgt(2))%k(i)
end do

p3 = k1 + k2 + k3 - p1 - p2

! Is p3 allowed by the size of the space?
testE = sum(p3**2)
if (abs(p3(1)) <= nmaxx .and. abs(p3(2)) <= nmaxy .and. &
abs(p3(3)) <= nmaxz .and. &
(.not. (tOrbECutoff .and. (testE > (OrbECutoff + 1.d-12))))) then
is_allowed = .true.
else
is_allowed = .false.
end if

if (.not. is_allowed) then
tgt(3) = 0
pgen = 0.0_dp
return
end if

i = kPointToBasisFn(p3(1), p3(2), p3(3), (msOrb + 1) / 2 + 1)

if (i > nBasis .or. i < 1) then
print *, 'bug kPointToBasisFn', p3, msOrb, i
do i = -1, 1
do j = -1, 1
do k = -1, 1
do l = -1, 1
print *, i, j, k, l, kPointToBasisFn(i, j, k, (l + 1) / 2 + 1)
end do
end do
end do
end do
call stop_all('pick_three_orbs_ueg', "kPointToBasisFn")
end if

! i occupied? check:
i_occup = .false.
if (i == tgt(1) .or. i == tgt(2)) then
i_occup = .true.
else
do j = 1, nel
if (i == nI(j)) i_occup = .true.
end do
end if

if (i_occup) then
tgt(3) = 0
pgen = 0.0_dp
return
else
tgt(3) = i
tgt = sort_unique(tgt)
end if

! adjust the probability by taking permutations into account
pgen = pgen * 2 * abs(ms)

end subroutine pick_three_orbs_ueg