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

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