gen_all_triples_k_space Subroutine

    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
                                                ! this excitation already..
                                                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