gen_parallel_double_hubbard Subroutine

    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