Just an update here, for the register. In recent days I have improved the implementation of CellListMap.jl, which is an alternative to compute the exact same thing that was desired by the OP. It is also general in the sense that it can compute other properties. With the current version (0.4.3) I obtain a similar performance to that of halotools for dense systems, and for constant density and increasing number of particles the package scales really well (linearly) and is quite efficient:
I think now we can say that we have a good alternative for the original OP question, for sure.