How to generate Hubbard Hamiltonian? I am a beginner with Hubbard Hamiltonian and my question is very basic:
how can I generate the matrix form of the Hubbard Hamiltonian?
I know the theory but I don't know how to put it numerically.
$$\hat{H} = -\sum_{\langle ij \rangle\sigma} (\hat{c}_{i\sigma}^{\dagger}\hat{c}_{j\sigma} + H.c.) + U\sum_i \hat{n}_{i\uparrow}\hat{n}_{i\downarrow} + \sum_{i\sigma} V_i\hat{n}_{i\sigma}$$
I mean, assuming that $\hat{H} = \hat{T} + \hat{U} + \hat{V}$, how can I generate each Hamiltonian?? 
Thanks in advance.
 A: To be able to write this into a matrix form, we need a priori to escape from the many-body picture. To do so we are going to suppose that for a reason due to the physics of the problem, the spin-down fermions have no dynamics (the spin is represented by $\sigma$, we are facing a Fermi-Hubbard model). This means hat we can replace $\hat n_{\downarrow}$ by its mean value $n_{\downarrow}$.
From this we see that the total Hamiltonian $\hat H$ can be decomposed in the following way:
$$ T^{(\uparrow)} = -\sum_{<i,j>} \hat c_{i\uparrow} \hat c^\dagger_{j\uparrow} + H.c. $$
$$ V^{(\uparrow)} = \sum_i \hat n_{i\uparrow} V_i.$$
Their spin-down counterpart are simply (sorry for the redundancy):
$$ T^{(\downarrow)} = -\sum_{<i,j>} \hat c_{i\downarrow} \hat c^\dagger_{j\downarrow} + H.c. $$
$$ V^{(\downarrow)} = \sum_i \hat n_{i\downarrow} V_i $$
$$ U^{(\uparrow)} = U \sum_i \hat n_{i\uparrow} n_{i\downarrow}.$$
This allows us to separate the spin-up and spin-down Hamiltonians:
$$ \hat H^{(\uparrow)} = T^{(\uparrow)} + V^{(\uparrow)} + U^{(\uparrow)}$$
$$ \hat H^{(\downarrow)} = T^{(\downarrow)} + V^{(\downarrow)}. $$
Note that $\hat n_{i\downarrow}$ has lost its hat. By writing this we assume that there is no back-action from the spin-up fermions on the spin-down fermions. Here we just used the assumption of no dynamics for the spin-down fermions.
From this we can first diagonalize $\hat H^{(\downarrow)}$ and pick our favorite eigenstate or linear combination of eigensates, take the absolute squared value (which is the density $n_{\downarrow}$) and plug it into $\hat H^{(\uparrow)}$. We then diagonalize $\hat H^{(\uparrow)}$ if we want to know its spectrum. Or we use $\hat H^{(\uparrow)}$ to compute the time-evolution of the spin-up fermions in presence of spin-down fermions as a background.
A: As  the occupation number of each the Fermi states is 0 or 1, when you have an $N$-site Hubbard model then the Hilbert space has dimension $2^{N}$ and the  Hamiltonian matrix will   be $2^N$-by-$2^N$. So, for a rather  small model with say a 10-by-10 lattice--- i.e. 100 sites--- the matrix is  $2^{100}$-by-$2^{100}$. As $2^{10}\approx 10^3$  this is a $10^{30}$-by-$10^{30}$ matrix. Do you plan to diagonalize this matrix numerically?
