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

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