quantum mechanics - 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?