H is the Hamiltonian. Here we set µ=0. For grand-canonical calculations add -6µ to the main diagonal.

-4*h + U	-2*Sqrt[2]*t	0	0
-2*Sqrt[2]*t	-4*h + U	-2*t	0
0	-2*t	-4*h + U	-2*Sqrt[2]*t
0	0	-2*Sqrt[2]*t	-4*h