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

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