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

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