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

2*h + U	-(Sqrt[2]*t)	Sqrt[3]*t	-2*t
-(Sqrt[2]*t)	2*h + U	0	0
Sqrt[3]*t	0	2*(h + U)	0
-2*t	0	0	2*h + U