Closest Pair of Points Algorithm - Fortune and Hopcroft
One thing that you might not be aware of is that there's a 'hidden' exponential factor (it's explained on page 3). In $k$ dimensions the runtime is $O(kn \log(\log(n)) + 3^kn)$, as when you partition $k$-dimensional space into hypercubes each cube has $3^k - 1$ neighbors.
The rest of the algorithm really is 'obvious' then: you split your space up in hypercubes (which is just a fancy word for splitting up each dimension separately into intervals). Each point will then map to exactly one cube, and the rest of the algorithm follows.
E.g. suppose our interval size is $1$ and our domain is $[0, 10)$. Then in 1D we have that the point $3.4$ maps to bucket $3$ and $8.8$ maps to bucket $8$. In three dimensions if our interval size and (per-dimension) domain is the same, we have that $(2.3, 8.4, 4.4)$ maps to bucket $(2, 8, 4)$ and $(4.1, 5.2, 9.5)$ maps to $(4, 5, 9)$. So you end up with three-dimensional buckets but the algorithm hashes the buckets afterwards so that doesn't matter.

If you have anything more than two or three dimensions my suggestion would be to look for another algorithm.