create_hilbert_space_kspace Subroutine

    subroutine create_hilbert_space_kspace(n_alpha, n_beta, n_orbs, nI, &
                                           n_states, state_list_ni, state_list_ilut)
        ! new functionality to create all states of the hilber space in the
        ! k-space hubbard model. reuse stuff from the real-space and apply
        ! the momentum conservation criteria!
        integer, intent(in) :: n_alpha, n_beta, n_orbs, nI(nel)
        integer, intent(out) :: n_states
        integer, intent(out), allocatable :: state_list_ni(:, :)
        integer(n_int), intent(out), allocatable :: state_list_ilut(:, :)
#ifdef DEBUG_
        character(*), parameter :: this_routine = "create_hilbert_space_kspace"
#endif

        integer(n_int), allocatable :: all_dets(:), temp_list_ilut(:, :)
        integer :: i, j, nJ(nel), n_allowed, k_vec_in(3), k_vec(3)
        type(Symmetry) :: sym_prod, sym_in
        integer, allocatable :: temp_list_ni(:, :)
        integer(n_int) :: temp_ilut(0:niftot)
        integer :: save_states

        ! this creates all possible basis states:
        all_dets = create_all_dets(n_orbs, n_alpha, n_beta)

        sym_in = G1(nI(1))%Sym
        k_vec_in = G1(nI(1))%k
        do i = 2, nel
            sym_in = SymTable(sym_in%s, G1(nI(i))%Sym%s)
            ! also add the k-vecs in a way so they do not leave the
            ! first BZ
            k_vec_in = lat%add_k_vec(k_vec_in, G1(nI(i))%k)

            ! a sanity check to see if it remained.. can be removed later!

        end do

        allocate(temp_list_ilut(0:niftot, size(all_dets)))
        allocate(temp_list_ni(nel, size(all_dets)))

        temp_list_ilut = 0_n_int
        temp_list_ni = 0

        ! now i want to loop over it to check the total momentum
        ! how do i most efficiently determine the symmetry?
        ! since adding the k-vector is not an option or? how far can the
        ! momentum then go in terms of Brillouin zones? can i efficiently
        ! map it back? or i just use the symmetry table, which has to be
        ! setup beforehand!
        n_allowed = 0

        do i = 1, size(all_dets)
            temp_ilut = all_dets(i)
            call decode_bit_det(nJ, temp_ilut)
            sym_prod = G1(nJ(1))%sym
            k_vec = G1(nJ(1))%k
            do j = 2, nel
                ! and i have to decode into the nI representation..
                ! oh god this is inefficient..
                sym_prod = symtable(sym_prod%s, G1(nJ(j))%Sym%s)
                k_vec = lat%add_k_vec(k_vec, G1(nJ(j))%k)

            end do
            if (sym_prod%s == sym_in%s) then
                ! if the symmetry fits:
                ASSERT(all(k_vec == k_vec_in))
                n_allowed = n_allowed + 1
                temp_list_ilut(:, n_allowed) = temp_ilut
                temp_list_ni(:, n_allowed) = nJ
            else
                ASSERT(.not. all(k_vec == k_vec_in))
            end if
        end do

        n_states = n_allowed

        allocate(state_list_ni(nel, n_allowed), source=temp_list_ni(:, 1:n_allowed))
        allocate(state_list_ilut(0:niftot, n_allowed), source=temp_list_ilut(:, 1:n_allowed))

        if (tHPHF) then
            ! for hphfs we need to purify the hilbert space!
            save_states = n_states
            if (allocated(temp_list_ilut)) deallocate(temp_list_ilut)

            allocate(temp_list_ilut(0:NIfTot, n_states), source=state_list_ilut)

            call spin_purify(save_states, temp_list_ilut, n_states, state_list_ilut)

        end if

    end subroutine create_hilbert_space_kspace