For particles of uniform size, a simple uniform grid will work best. If you set the grid cell size to the particle diameter, then you only need to search the adjacent voxels for each particle. The uniform grid method is also easy to port to the GPU, unlike other acceleration structures like BVHs or Octrees.
