# Numerical exact diagonalization of tight binding Hamiltonian

I want to exactly diagonalize the following Hamiltonian for $$10$$ number of sites and $$4$$ number of spinless fermions $$H = -t\sum_i^{L-1} \big[c_i^\dagger c_{i+1} - c_i c_{i+1}^\dagger\big] + V\sum_i^{L-1} n_i n_{i+1}$$ here $$L$$ is total number of sites, creation ($$c^\dagger$$) and annihilation ($$c$$) operators are defined as following $$c = \begin{bmatrix} 0&0\\1&0 \end{bmatrix}$$ and $$n_i = c_i^\dagger c_i$$ is number operator.

To exactly diagonalize (for simplicity let's take $$L=4$$ sites), one can expand $$H$$ as

$$H = -t\big[ c_1^\dagger \sigma_1^z \otimes c_2\otimes I_3 \otimes I_4 \\ + I_1 \otimes c_2^\dagger \sigma_2^z \otimes c_3\otimes I_4\\ + I_1 \otimes I_2 \otimes c_3^\dagger\sigma_3^z \otimes c_4 \big]+h.c.\\ +V\big[ n_1 \otimes n_2 \otimes I_3 \otimes I_4\\ +I_1 \otimes n_2 \otimes n_3 \otimes I_4\\ +I_1 \otimes I_2 \otimes n_3 \otimes n_4 \big]$$ where $$\sigma^z$$ (Pauli matrix) is just simple matrix multiplication for the sake of anti-commutation relation.

So far so good. (please correct me if I am doing anything wrong)!

Question:

I used the above method and numerically calculated the ground state and found that above method gives correct results for $$V=0$$ but when $$V\ne0$$ the results are wrong.

Eventually, I get to the point that I am not taking care of number of particles in the system. How do we numerically diagonalize a Hamiltonian matrix in the sector with chosen number of particles?