Something like this should give you what you want.

radius /= gl_ProjectionMatrix[2][3] * position.z + gl_ProjectionMatrix[3][3];

Basically you want the final radius in screen space (what the shader expect) but have only the world space radius. The code above should do the conversion (I haven't tested).

Remember to test in persp and ortho views.

But I'm not sure about replacing sphere by points though.
Points will have depth problem on viewfacing surfaces and will be cut in half most of the time.
So either use points and disable depth test, and let point leak bellow surfaces.
Or use instancing like I did in clay engine to instance a sphere batch. (might be less trouble than all of this)