Complex values for the dispersion relation obtained through an $s$-band only tight binding model for diamond cubic crystal

Any given atom in a diamond cubic lattice (Like Si or Ge) has four nearest neighbours at at a distance $$\sqrt{3}a/4$$, being $$a$$ the lattice constant. The translation vectors to these neighbours can be chosen as

$$$$(1,1,1)\frac{a}{4}, {\phantom{x}} (1,-1,-1)\frac{a}{4}, {\phantom{x}} (-1,1,-1)\frac{a}{4}, {\phantom{x}} (-1,-1,1)\frac{a}{4},$$$$

To obtain the vectors above I consider the left face centred site and considered the x-axis positive to the right, the y-axis positive upwards and the z-axis coming out of the screen

Now, the s-band only dispersion relation given by the tight-binding method is

$$$$E(\vec{k}) = e_{s} + \Sigma_{\vec{\tau}}\gamma(|\tau|)e^{i\vec{k}\cdot\vec{\tau}}$$$$

where $$\tau$$ is the vectors described above to the nearest neighbours. This summation will result in

$$$$E(\vec{k}) = e_{s} +\gamma(|\tau|)(e^{i(k_{x}+k_{y}+k_{z})a/4}+e^{i(k_{x}-k_{y}-k_{z})a/4}+e^{i(-k_{x}+k_{y}-k_{z})a/4}+e^{i(-k_{x}-k_{y}+k_{z})a/4})$$$$

My issue is the following: When you do the same procedure for the face-centred cubic crystal you have 12 nearest neighbours, but the summation can be simplified to an explicitly real expression, depending only on combinations of cosines. In this case the expression seems to always have an imaginary component. Where have I made a mistake?

edit.: the expression multiplying the $$\gamma$$ factor can be reduced to $$2\cos(k_x\frac{a}{4})\cos(k_y\frac{a}{4})\cos(k_z\frac{a}{4})-2i\sin(k_x\frac{a}{4})\sin(k_y\frac{a}{4})\sin(k_z\frac{a}{4}).$$

The formula you have given works only if there is a single orbital in the basis of the lattice. Diamond has zinc blende structure, meaning a face-centred cubic (FCC) lattice with a basis of atoms at $$[0,0,0]$$ and $$\left[\frac{1}{4},\frac{1}{4},\frac{1}{4}\right]$$ (in Cartesian coordinates and units of the lattice constant). Thus there are two atoms in the basis, and the strategy is instead to solve the matrix equation $$E(\boldsymbol{k})c_i=\sum_j c_j \sum_{\boldsymbol{R}}\gamma_{ij} \,e^{i\boldsymbol{k}\cdot\boldsymbol{R}}$$ where the sum is over atoms in the basis $$i$$ and lattice vectors $$\boldsymbol{R}$$, and $$\gamma_{ij}$$ is the coupling matrix element between orbital $$i$$ and $$j$$ (any solid state physics book should derive something like this). The tight binding approximation then assumes that only the nearest neighbour couplings are non-negligible, in our case given by a constant $$\gamma$$.
Importantly, note that we use the lattice vectors whose basis contains the nearest neighbours, not the coordinates of the nearest neighbours themselves. The nearest neighbours to the atom at $$[0,0,0]$$ are correctly as you listed, but we must find the lattice vectors to which they belong. In the FCC lattice, these are $$[0,0,0],\quad \left[0,-\frac{1}{2},-\frac{1}{2}\right], \quad \left[-\frac{1}{2},0,-\frac{1}{2}\right], \quad \left[-\frac{1}{2},-\frac{1}{2},0\right].$$ In other words, the above lattice vectors plus $$\left[\frac{1}{4},\frac{1}{4},\frac{1}{4}\right]$$ give back the coordinates you listed. Similarly, we must find the lattice vectors to which the nearest neighbours of the atom at $$\left[\frac{1}{4},\frac{1}{4},\frac{1}{4}\right]$$ belong; a bit of thought shows that these are minus the above.
Substituting these into the matrix equation shows that the energies of the bands are given by the eigenvalues of the matrix $$\gamma\begin{pmatrix} 0 & 1 + e^{-i(k_y+k_z)a/2} + e^{-i(k_x+k_z)a/2}+e^{-i(k_x+k_y)a/2} \\ 1 + e^{i(k_y+k_z)a/2} + e^{i(k_x+k_z)a/2}+e^{i(k_x+k_y)a/2} & 0 \end{pmatrix}$$ where we have ignored any diagonal term (just giving an overall constant). I won't do this, but note that it is Hermitian so the eigenvalues must be real. Note also that there are two different bands deriving from the two atoms in the basis.