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

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