The whole problem becomes far simpler if you re-base your points (both those you need to rotate, and those that form stopping points) so that they are relative to your axis of rotation. Ie. move the entire field of view so that the axis unit vector lines up with the Z-axis of your re-based coordinate system.

You can then subset your stopping points against your input points by their radius.

That is, once the points are re-based, sqrt( x² + y² ) for each of the points will give you their radius of rotation. Now you only need compare any given input point against those stopping points that have that same radius (+-some tolorance).

And if once you've re-based the points, you convert them to 2D vectors, the length is the radius directly, and the difference between their angles tells you how far you would have to rotate them before they collided.

So now you don't need to rotate all the points incrementally at all, you simply find all the pairs of input/stopping points that can collide; and choose the minimum or maximum angle (depending upon which direction you are rotating them, as your "first collision point".

This is a much simpler, and far more computationally efficient way of folding your molecule.