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)	0	-2*t	0	0	0	(2*t)/Sqrt[3]
-(Sqrt[2/3]*t)	2*h + U	Sqrt[6]*t	0	-(Sqrt[3]*t)	0	(2*t)/Sqrt[3]	0
0	Sqrt[6]*t	2*(h + U)	0	0	0	0	0
-2*t	0	0	2*(h + U)	-(Sqrt[2]*t)	0	0	0
0	-(Sqrt[3]*t)	0	-(Sqrt[2]*t)	2*(h + U)	2*t	0	0
0	0	0	0	2*t	2*(h + U)	0	0
0	(2*t)/Sqrt[3]	0	0	0	0	2*h + U	-2*Sqrt[2/3]*t
(2*t)/Sqrt[3]	0	0	0	0	0	-2*Sqrt[2/3]*t	2*h + U