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

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