Derivation of Matrix Components of Hamiltonian in Tight Binding Method Im currently struggling with the description of the tight binding method in the original paper by Slater and Koster from 1954 (where a free version of the paper can be found under this link).
In principle I understand the concept of summing over the nearest neighbours and building up the Hamiltonian from these sums.
However in Table II of the original paper, Slater and Koster give these Hamiltonian elements as sums of three center integrals like ($E_{s,x}(p,q,r), E_{s,x^2-y^2}(p,q,r)$) for the simple cubic crystal as an example.
I understand how for instance the element $(s,s)$ in Table II gets constructed, however the element $(s,x)$ is unclear. It seems like that depending on the symmetry of the two orbitals centered on the i and j atom position just some three center integrals get selected and others do not enter the sum. Could anybody explain the general method of building up the Table II?
UPDATE: To be more clear and not to mix up the three-center approximation with the two-center approximation I'm listing two matrix elements of the Table II as an exampe:
$(s/x)=2\mathbb{i}E_{s,x}(100)\sin{\xi}+4\mathbb{i}E_{s,x}(110)(\sin{\xi}\cos{\eta}+\sin{\xi}\cos{\zeta})+8\mathbb{i}E_{s,x}(111)\sin{\xi}\cos{\eta}\cos{\zeta}$
As a comparison the $(s,s)$ element shows all nearest neighbor elements of the Bloch summation
$(s/s)=E_{s,s}(000)+2E_{s,s}(100)(\cos{\xi}+\cos{\eta}+\cos{\zeta})+4E_{s,s}(110)(\cos{\xi}\cos{\eta}+\cos{\xi}\cos{\zeta}+\cos{\eta}\cos{\zeta})+8E_{s,s}(111)\cos{\xi}\cos{\eta}\cos{\zeta}$
The Integrals $E_{s,s}(p,q,r)$ and $E_{s,x}(p,q,r)$ in above elements are three center integrals which can be substituted with the expressions in Table I of the reference to yield the two-center integral Hamiltonian components.
FURTHER UPDATE:
The original paper states that the expressions in Table II have been built up by minimizing the number of three center integrals by using symmetry considerations. One such relation I found already namely:
$E_{s,x}(-1,0,0)=-E_{s,x}(1,0,0)$ and the equivalente relations for the q and r coordinates which result in the first term of (s/x) as
$2\mathbb{i}\sin{\xi}E_{s,x}(100)$ as expected.
However for the $E_{s,x}(110)$ I'm already lost again....
FURTHER UPDATE:
The nice answer by Timok shows how to resolve the problem. There is additional fineprint to the topic.
By using tesseral functions (which can be used the basis functions of the Slater Coster approach and are in fact the real representation of the spherical harmonics. Refer also to this Wikipedia article) one can inspect directly the symmetry relations for the basis functions.
E.g. the tesseral function for the xy orbital is
$\frac{1}{4}\sqrt{\frac{15}{\pi}}\Im(e^{2 \mathbb{i} \phi}\sin^2{\theta})$
which is in cartesian coordinates
$\frac{1}{2}\sqrt{\frac{15}{\pi}} x y$.
Thus we can see immediately that by changing the sign of the coordinates x y z in any possible permutation we get the symmetry relation:
$E_{s,xy}(x,y,z)=Sign(x)Sign(y)E_{s,xy}(|x|,|y|,|z|)$
FURTHER UPDATE:
The usual definition of the tesseral functions for the slater coster integrals are nicely working except the $d_{3z^2-r^2}$ orbital which results in non-orthogonality of the matrix element $<s|d_{3z^2-r^2}>$ when using the tesseral definition with $l=2$ and $m=0$ giving $\sqrt{\frac{5}{16\pi}}(3\cos^2\theta-1)$. However this non-orthogonality can be corrected if one uses a modified term for $d_{3z^2-r^2}=\frac{3}{8}\sqrt{\frac{5}{\pi}}\cos{2\theta}$.
Now the big question remains: Which term has been used by Slater and Coster? All newer publications have been using the standard tesseral definition which is not completely orthogonal!
Kind regards,
Rainer 
 A: An important discrimination in the integrals is by looking at the zero values of the $E's$ in the table $I$ for particular values of the $(l,m,n)$, because the $(l,m,n)$  of table $I$ are proportionnal to the $(p,q,r)$ of table $II$ (see page $1504$ of the ref): 
a) if you compare $(s/s)$ and $(s,x)$ in table $II$, the difference is that you have not the integral $E_{s,x}(p=0,q=0,r=0)$ in the development of $(s,x)$
But, looking at the table $I$, we  have $E_{s,x}(l,m,n) = l ~(sp\sigma)$, so  $E_{s,x}(p=0,q=0,r=0) = E_{s,x}(l=0,m=0,n=0)=0$.
b) If you look at $(x,y)$, you have not $E_{x,y}(p=0,q=0,r=0)$ , $E_{x,y}(p=1,q=0,r=0)$, $E_{x,y}(p=0,q=1,r=0)$ terms. But we see in table $I$, that  : $E_{x,y}(l,m,n) = lm ((pp\sigma) - (pp\pi))$.
So, for instance you have $E_{x,y}(p=1,q=0,r=0) = E_{x,y}(l=\lambda,m=0,n=0)=0$, and idem for the other terms.
This is not the full story, of course...
[EDIT]
I am using the ideas given in this paper, (see for instance, example $2$, Lithium).
The tight-binding wave-function : 
$$\psi_{j, \vec k}(\vec r) = \sum_{p,q,r} e^{i(p \vec k.\vec a_1 + q \vec k.\vec a_2 +r \vec k.\vec a_3)} \phi_j (\vec r - p \vec a_1 - q \vec a_2 - r \vec a_3)$$
where $i,j$ are orbital indices.
We have $ E_{ij} = \langle \phi_j|\hat H| \psi_{i}\rangle $
In the tight binding approximation, only the on-site and nearest neighbor matrix elements are retained.
$E_{ij} = \sum_{(on-site ~+~ nearest ~~neighbours)~~p,q,r} e^{i(p \vec k.\vec a_1 + q \vec k.\vec a_2 +r \vec k.\vec a_3)} \langle \phi_j| \hat H|\phi_i(p,q,r) \rangle \tag{1}$
(s,x) and its symmetries
Now, considering $(s,x)$, we are going to use symmetries, that is : 
$\langle \phi_s| \hat H|\phi_x(\epsilon_1 p, \epsilon_2  q, \epsilon_3 r) \rangle = (-1)^{\epsilon_1} \langle \phi_s| \hat H|\phi_x( p,  q,  r) \rangle \tag{2a}$ 
$\langle \phi_s| \hat H|\phi_x( p,  r,  q) \rangle =  \langle \phi_s| \hat H|\phi_x( p,  q,  r) \rangle \tag{2b}$ 
$\langle \phi_s| \hat H|\phi_x( 0,  q,  r) \rangle =  0\rangle \tag{2c}$ 
Note that $(2c)$ is a consequence of $(2a)$
Nearest neighbours
We only consider $p,q,r = (1,0,0),(-1,0,0)$. Call $E_{s,x}(1,0,0) = \langle \phi_s| \hat H|\phi_x(1,0, 0) \rangle$. From formula $(1)$ and the symmetries relations $(2)$, we get : 
$(s,x)_1 = e^{ik_1 a_1} E_{s,x}(1,0,0) + e^{-ik_1 a_1} E_{s,x}(-1,0,0) \\ = 
E_{s,x}(1,0,0) (e^{ik_1 a_1} - e^{-ik_1 a_1}) \\ = 2i E_{s,x}(1,0,0) sin (k_1 a_1) \\ 
= 2i E_{s,x}(1,0,0) sin (\xi)$
Second Nearest neighbours
We only consider $(p,q,r) = {(1,1,0), (1,-1,0), (-1,1,0), (-1,-1,0)} + {q \leftrightarrow r}$. So, we have : 
$(s,x)_2 = e^{ik_1 a_1 +ik_2 a_2} E_{s,x}(1,1,0) + e^{ik_1 a_1 - ik_2 a_2} E_{s,x}(1,-1,0) + e^{-ik_1 a_1 +ik_2 a_2} E_{s,x}(-1,1,0) + e^{-ik_1 a_1 - ik_2 a_2} E_{s,x}(-1,-1,0) + \{q \leftrightarrow r, k_2 a_2 \leftrightarrow k_3 a_3 \}$ 
From formula $(1)$ and the symmetries relations $(2)$, we get : 
$(s,x)_2 = E_{s,x}(1,1,0)[e^{ik_1 a_1 +ik_2 a_2} + e^{ik_1 a_1 - ik_2 a_2}  - e^{-ik_1 a_1 +ik_2 a_2} - e^{-ik_1 a_1 - ik_2 a_2}] + \{ k_2 a_2 \leftrightarrow k_3 a_3 \} \\ 
= 4i E_{s,x}(1,1,0) (\sin k_1 a_1 \cos k_2 a_2 + \{ k_2 a_2 \leftrightarrow k_3 a_3 \}) \\
= 4i E_{s,x}(1,1,0) (\sin \xi \cos \eta + \sin \xi  \cos \zeta)$
