# gen_parallel_double_hubbard Subroutine

## public subroutine gen_parallel_double_hubbard(nI, ilutI, nJ, ilutJ, ex, tParity, pgen)

### Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nI(nel)
integer(kind=n_int), intent(in) :: ilutI(0:niftot)
integer, intent(out) :: nJ(nel)
integer(kind=n_int), intent(out) :: ilutJ(0:niftot)
integer, intent(out) :: ex(2,2)
logical, intent(out) :: tParity
real(kind=dp), intent(out) :: pgen

## Source Code

    subroutine gen_parallel_double_hubbard(nI, ilutI, nJ, ilutJ, ex, tParity, pgen)
integer, intent(in) :: nI(nel)
integer(n_int), intent(in) :: ilutI(0:niftot)
integer, intent(out) :: nJ(nel), ex(2, 2)
integer(n_int), intent(out) :: ilutJ(0:niftot)
logical, intent(out) :: tParity
real(dp), intent(out) :: pgen
real(dp) :: p_elec, p_orb
integer :: elecs(2), orbs(2), src(2)

! in the transcorrelated case we have to decide
! i first have to choose an electron pair (ij) at random
! but with the condition that they have to have opposite spin!
! this is the only difference: i pick two spin-parallel electrons..
! the rest stays the same.. i just have to adjust the
! get_orb_from_kpoints routine
! and the matrix element calculation
call pick_spin_par_elecs(nI, elecs, p_elec)

src = nI(elecs)

! i realised i could reuse the already implemented orbital picker,
! but the other one does not use the fact that we know that we
! have parallel spin in this case! so implement a parallel
! spin-orbital picker!! (also better for matrix elements.. so we
! must not check if the spins fit, if we only take the correct ones!)
call pick_ab_orbitals_par_hubbard(nI, ilutI, src, orbs, p_orb)

if (orbs(1) == ABORT_EXCITATION) then
nJ(1) = ABORT_EXCITATION
pgen = 0.0_dp
return
end if

! and make the excitation
call make_double(nI, nJ, elecs(1), elecs(2), orbs(1), orbs(2), ex, tParity)

ilutJ = make_ilutJ(ilutI, ex, 2)

! i think in both the electrons and the orbitals i have twice the
! probability to pick them
! already modified in the orbital picker!
! i am not super sure about a factor of 2 here..
pgen = p_elec * p_orb

end subroutine gen_parallel_double_hubbard