The problem is that the brute-force "nested loops" algoritm scales as N^2. If you have 1000 atoms, you need almost 500,000 pairwise distance calculations! The solution is to divide-and-conquer by recursively partitioning the space into smaller boxes with a small number of atoms. That way you can get an almost N^1 scaling.

Take a look at the code in Chemistry::Bond::Find for an example. I just released it yesterday (a very early version...) because I had almost the same problem (maybe I was reading your mind?). With the naive N^2 algorithm it was taking almost three minutes to find the bonds in a protein with 1717 atoms; with the recursive partitioning it took only 10 seconds. Note, however, that if your cutoff distance is relatively large compared to the total system size, you won't get such a big improvement. Since I wanted to find the bonds, cutoffs were around 2 Å.