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

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