# gen_all_triples_k_space Subroutine

## public subroutine gen_all_triples_k_space(nI, n_excits, det_list)

### Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nI(nel)
integer, intent(out) :: n_excits
integer(kind=n_int), intent(out), allocatable :: det_list(:,:)

## Source Code

    subroutine gen_all_triples_k_space(nI, n_excits, det_list)
integer, intent(in) :: nI(nel)
integer, intent(out) :: n_excits
integer(n_int), intent(out), allocatable :: det_list(:, :)

integer :: i, j, k, a, b, c, spin, src(3), ex(2, 3), n_bound, pos
integer(n_int) :: ilut(0:niftot), ilutJ(0:NIfTot)
integer(n_int), allocatable :: temp_list(:, :)
real(dp) :: elem

! write a routine to create all the triple excitations for the
! transcorrelated k-space hubbard model. for a given determinant
call EncodeBitDet(nI, ilut)

! do it in the most naive way for now..
! i need to loop over all different electron combinations, with the
! restriction that they are not all parallel spin (due to the form
! of the transcorrelated hamiltonian)
n_excits = 1

! todo: an estimate for the upper bound of number of triple excitations..
! this gets too high for big lattices..

! i think a more correct estimat is:
n_bound = int(nel * (nel - 1) * (nel - 2) * (nbasis - nel) * (nbasis - nel - 1) / 8)

allocate(temp_list(0:niftot, n_bound))
temp_list = 0_n_int

do i = 1, nel - 2
do j = i + 1, nel - 1
do k = j + 1, nel
! assert that the they are not all spin-parallel:
src = nI([i, j, k])
spin = sum(get_spin_pn(src))
ex(1, :) = src

if (spin == -1 .or. spin == 1) then
do a = 1, nbasis - 1
if (IsNotOcc(ilut, a)) then
do b = a + 1, nbasis
if (IsNotOcc(ilut, b)) then

! this gives me the correct spin already
! to have the same as the electrons!
c = get_orb_from_kpoints_three(src, a, b)

! i also have to check if no orbital is
! picked twice
if (IsNotOcc(ilut, c) .and. c /= a .and. c /= b) then
! this should be a valid excitation or?

ex(2, :) = [a, b, c]
! should i check the matrix element too?
! and also be sure that the momentum fits
elem = abs(get_helement_lattice(nI, 3, ex, .false.))
if (elem > EPS) then
ilutJ = make_ilutJ(ilut, ex, 3)

! and to be save, search if we have
pos = binary_search_ilut(temp_list(:, 1:n_excits), ilutJ, nifd + 1)

if (pos < 0) then
! new excitation
temp_list(:, n_excits) = ilutJ
n_excits = n_excits + 1
end if
end if
end if
end if
end do
end if
end do
end if
end do
end do
end do

n_excits = n_excits - 1
allocate(det_list(0:NIfTot, n_excits), source=temp_list(:, 1:n_excits))

end subroutine gen_all_triples_k_space