Thank you for helping me out.
Yes, by replacing the expression to index_add_, the results are more or less the same. Although there is some round-off differences between CPU and GPU, these errors are of the order 1e-8, so that’s the behavior that I expected.
For completeness: I had to replace the command
tuss[j,:,self.coordlist[:]] += x[j,:,:]

with:
tuss[j].index_add_(1, self.coordlist, x[j])

I’ve also noticed that the old version yielded wrong results on the CPU as well.