When accessing a numpy multi-dimensional array with other multi-dimensional arrays of integer type the arrays used for the indices need to have the same shape. 
Numpy will happily broadcast, if possible - but for that to be possible the arrays need to have the same dimensionality, e.g. this works:
sumkma[site[j], j[:,np.newaxis]]
The np.newaxis results in j[:,np.newaxis] being two-dimensional (shape is now (20,1), whereas shape of j is one-dimensional (20,)). It now has a shape that can be broadcasted to the shape of site[j]:
>>> j.shape
(20,)
>>> site[j].shape
(20,200)
>>> j[:,np.newaxis].shape
(20,1)
The same dimensionality for the index arrays allows numpy to broadcast them to have the same shape (20,200).
https://docs.scipy.org/doc/numpy-1.13.0/user/basics.indexing.html#indexing-multi-dimensional-arrays