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

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