Is there an algorithm that takes atom coordinates and produces VSEPR geometries?
Quantum chemistry software usually exploits molecular symmetries to speed up computations. To do so they need some way to determine the point group for a given set of atoms and their coordinates. So this is basically what you are asking about, just that get something like $D_{4h}$ instead of "square planar", or $T_d$ instead of "tetrahedral".

I am not sure whether this is actually used by the mentioned software, but it is quite simple and efficient if you only want to test for certain point groups:


Each symmetry operation can be represented as a $3\times 3$ matrix, for example mirroring in the $xy$ plane is 


\begin{equation}\sigma_{xy} = \begin{pmatrix}
 1 &amp; 0 &amp; 0 \\
 0 &amp; 1 &amp; 0 \\
 0 &amp; 0 &amp; -1
\end{pmatrix}\end{equation}


We can represent the atom coordinates in a $3\times N$ matrix, where $N$ is the number of atoms. The molecule needs to be properly oriented (the center of mass should be in the origin of the coordinate system and the principle axes of inertia should coincide with the $x$, $y$ and $z$ axes).
For example $\ce{H2}$ (or any other homonuclear diatomic):


\begin{pmatrix}
 0 &amp; 0 \\
 0 &amp; 0 \\
 1 &amp; -1
\end{pmatrix}


To get the transformed coordinates according to a certain symmetry operation, we just apply the matrix of the symmetry operation to the coordinates matrix.


\begin{equation}\begin{pmatrix}
 1 &amp; 0 &amp; 0 \\
 0 &amp; 1 &amp; 0 \\
 0 &amp; 0 &amp; -1
\end{pmatrix}\times
\begin{pmatrix}
 0 &amp; 0 \\
 0 &amp; 0 \\
 1 &amp; -1
\end{pmatrix}=
\begin{pmatrix}
 0 &amp; 0 \\
 0 &amp; 0 \\
 1 &amp; -1
\end{pmatrix}
\end{equation}


We can then compare the resulting set of coordinates with the original one. Since for $\ce{H2}$ both atoms are the same, we can interchange the 2 columns and see, that we have identical matrices. Thus, $\sigma_{xy}$ is a valid symmetry element of our input coordinates. At this point you may implement a tolerance threshold, by allowing the norm (distance) between two "same" points being larger than 0.
The set of found symmetry operations will then directly give us the point group of that molecule (this information is available in character tables). So we can check for a certain point group by checking whether the symmetry operations of that point group can be found in our input geometry. This however means we need to define which point groups we want to test for. 


There are also flow charts (available in books about molecular symmetry or via your favorite search engine), telling you which symmetries elements to check for and then guide you to the point group.

I am not aware of a more general approach. If someone knows more, please extend.