My thought for an improved implementation was to throw the points into a
voxel grid and then have an iterator which traverses the voxel grid and
generates pairs on the fly (rather than generating and storing them).
Well, that sounds essentially like a cell-based nonbonded list, which is
the most common good way of implementing a nonbonded list (and is what
Modeller does) except for the bit about generating pairs on the fly.
How far the iterator searches for the second point for each voxel would
be determined by the distance cutoff. It seems easy enough to implement.
The trick is choosing the voxel size sensibly, akin to choosing a hash
bin size. This is a bit trickier if you have multiple iterators with
different cutoffs.
BTW, another thing to be aware of is that most MM force fields
(including CHARMM) exclude not just 1-2 pairs (bonds) but also 1-3 and
1-4. But I guess it would be pretty straightforward to exclude those
pairs too by following the bond graph.