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

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