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

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