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

-2*h + U	(I/12)*(13*I + Sqrt[3])*t	((3*I + 7*Sqrt[3])*t)/(6*Sqrt[2])	0	((-3*I + 5*Sqrt[3])*t)/(6*Sqrt[2])	0	0	0	0	0	0	0	0	((I/3)*(I + Sqrt[3])*t)/Sqrt[2]	0	-((-3*I + Sqrt[3])*t)/4
(-I/12)*(-13*I + Sqrt[3])*t	-2*h + U	0	-(t/Sqrt[6])	0	-((-3*I + Sqrt[3])*t)/(6*Sqrt[2])	0	0	0	0	((-9*I + 5*Sqrt[3])*t)/(6*Sqrt[2])	-((3*I + 5*Sqrt[3])*t)/(6*Sqrt[2])	(Sqrt[2]*(1 + I*Sqrt[3])*t)/3	0	((3*I + Sqrt[3])*t)/12	0
((-3*I + 7*Sqrt[3])*t)/(6*Sqrt[2])	0	-2*(h - U)	((3 - I*Sqrt[3])*t)/2	0	0	0	0	0	0	-t	t	((-3*I + Sqrt[3])*t)/3	0	((1 + I*Sqrt[3])*t)/(2*Sqrt[2])	0
0	-(t/Sqrt[6])	((3 + I*Sqrt[3])*t)/2	-2*(h - U)	0	0	0	0	0	0	0	0	0	(-2*t)/Sqrt[3]	0	I*Sqrt[3/2]*t
((3*I + 5*Sqrt[3])*t)/(6*Sqrt[2])	0	0	0	-2*(h - U)	(t + I*Sqrt[3]*t)/2	0	t	0	0	0	0	-((-3*I + Sqrt[3])*t)/3	0	((3 - I*Sqrt[3])*t)/(2*Sqrt[2])	0
0	-((3*I + Sqrt[3])*t)/(6*Sqrt[2])	0	0	(t - I*Sqrt[3]*t)/2	-2*(h - U)	-(Sqrt[2]*t)	0	(t + I*Sqrt[3]*t)/2	Sqrt[2]*t	0	0	0	-((3*I + Sqrt[3])*t)/3	0	((1 + I*Sqrt[3])*t)/(2*Sqrt[2])
0	0	0	0	0	-(Sqrt[2]*t)	-2*(h - U)	-(Sqrt[2]*t)	0	0	0	0	0	0	I*(I + Sqrt[3])*t	0
0	0	0	0	t	0	-(Sqrt[2]*t)	-2*(h - U)	(t + I*Sqrt[3]*t)/2	((1 + I*Sqrt[3])*t)/Sqrt[2]	0	0	0	0	0	(I*(I + Sqrt[3])*t)/Sqrt[2]
0	0	0	0	0	(t - I*Sqrt[3]*t)/2	0	(t - I*Sqrt[3]*t)/2	-2*(h - U)	0	-t	t	0	0	0	0
0	0	0	0	0	Sqrt[2]*t	0	((1 - I*Sqrt[3])*t)/Sqrt[2]	0	-2*(h - U)	0	0	0	0	0	0
0	((9*I + 5*Sqrt[3])*t)/(6*Sqrt[2])	-t	0	0	0	0	0	-t	0	-2*h + U	0	0	(2*t)/Sqrt[3]	0	((1 - I*Sqrt[3])*t)/(2*Sqrt[2])
0	((3*I - 5*Sqrt[3])*t)/(6*Sqrt[2])	t	0	0	0	0	0	t	0	0	-2*h + U	0	((3*I + Sqrt[3])*t)/3	0	((3 + I*Sqrt[3])*t)/(2*Sqrt[2])
0	(Sqrt[2]*(1 - I*Sqrt[3])*t)/3	((3*I + Sqrt[3])*t)/3	0	-((3*I + Sqrt[3])*t)/3	0	0	0	0	0	0	0	-2*h + U	(I/6)*(7*I + Sqrt[3])*t	0	0
((-I/3)*(-I + Sqrt[3])*t)/Sqrt[2]	0	0	(-2*t)/Sqrt[3]	0	-((-3*I + Sqrt[3])*t)/3	0	0	0	0	(2*t)/Sqrt[3]	((-3*I + Sqrt[3])*t)/3	(-I/6)*(-7*I + Sqrt[3])*t	-2*h + U	((3*I + Sqrt[3])*t)/(3*Sqrt[2])	0
0	((-3*I + Sqrt[3])*t)/12	((1 - I*Sqrt[3])*t)/(2*Sqrt[2])	0	((3 + I*Sqrt[3])*t)/(2*Sqrt[2])	0	(-1 - I*Sqrt[3])*t	0	0	0	0	0	0	((-3*I + Sqrt[3])*t)/(3*Sqrt[2])	-2*h + U	(I/4)*(5*I + Sqrt[3])*t
-((3*I + Sqrt[3])*t)/4	0	0	(-I)*Sqrt[3/2]*t	0	((1 - I*Sqrt[3])*t)/(2*Sqrt[2])	0	((-1 - I*Sqrt[3])*t)/Sqrt[2]	0	0	((1 + I*Sqrt[3])*t)/(2*Sqrt[2])	((3 - I*Sqrt[3])*t)/(2*Sqrt[2])	0	0	(-I/4)*(-5*I + Sqrt[3])*t	-2*h + U