Using Lennard-Jones potential to generate $n$ oxygen atom coordinates that are evenly distributed on a sphere

Question

I am a student majoring in Chemistry. In my daily life, I often write input files of computational chemistry software. Most of the input files contain the coordinates block as shown below. It is quite easy to write the coordinates manually if the total number of atoms is not big. However, if the molecule has $1000$ or even $10000$ atoms, inputting their coordinates is not a trivial job, and that is why I am here.

Elements      x             y           z
O           1.084802    -0.240802   2.133771
H           -0.774428   0.497993    1.88079
O           2.419005    0.899744    0.191657
C           0.739285    0.563107    -1.493075

Recently, I want to write a script that can generate $n(n>100)$ oxygen atom coordinates automatically. These oxygen atoms are evenly distributed on a sphere and form a cavity, like C60. First, I thought of the Lennard-Jones equation: \begin{equation}\tag{1} V(r) = 4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12} -\left(\frac{\sigma}{r}\right)^{6}\right] \end{equation} Briefly, the potential(V) between two atoms is a function of the distance(r) between them. Here, the constants $\epsilon$ and $\sigma$ are related to the element type. According to the function diagram, when $r$ is not too large and not too small, the potential can be minimized. In other words, $V(r)$ has a global minimum. These must be very simple to understand for you guys. ; ) Because of the following two conditions, the equation can be simplified.

1. Because I want to restrain the distance between oxygen atoms to 3 angstrom, using the $\frac{\partial{V}}{\partial{r_{min}}}=0$ and $r_{min}=3$, I can determine the $\sigma$ is about $2.673$.
2. Meanwhile, the "amplitude" of the equation is not important. That is I can set $4\epsilon = 1$.
   \begin{equation}\tag{2} V(r) = \left(\frac{2.673}{r}\right)^{12} -\left(\frac{2.673}{r}\right)^{6} \end{equation}

After getting equation 2, I construct a cost function: $$\tag{3} \Theta = \sum_{\{ij\}\in set}V(r_{ij}) $$ $\{ij\}\in set$ means that $r_{ij}$ are all the possible connections between oxygen atoms. Therefore, an adjacency matrix whose diagonal is 0 and other elements are 1 can be constructed.
The workflow of my script is as follows. First, the coordinates of all oxygen atoms are randomly initialized in range $(0, 1)$. Second, the gradient descendent method was employed to minimize $\Theta$. \begin{equation}\tag{4} \begin{split} \frac{\partial{V}}{\partial{x_i}} &= \left(x_j - x_i\right)\left(\frac{12\sigma^{12}}{r_{ij}^{14}} - \frac{6\sigma^{6}}{r_{ij}^{8}}\right) \\ \frac{\partial{V}}{\partial{y_i}} &= \left(y_j - y_i\right)\left(\frac{12\sigma^{12}}{r_{ij}^{14}} - \frac{6\sigma^{6}}{r_{ij}^{8}}\right) \\ \frac{\partial{V}}{\partial{z_i}} &= \left(z_j - z_i\right)\left(\frac{12\sigma^{12}}{r_{ij}^{14}} - \frac{6\sigma^{6}}{r_{ij}^{8}}\right) \end{split} \end{equation} During the minimization, I ran into two problems:

1. Gradient explosion. When $r$ is small, the gradient can be quite large, say $10^{12}$, and a reasonable $\alpha$ value is difficult to set.
2. Gradient vanishing. When $r$ is large, the gradient can be quite small, say $10^{-12}$, and a reasonable $\alpha$ value is difficult to set.

My question is :

1. Is it possible to set an appropriate $\alpha$?
2. Is my cost function appropriate?

Any insights would be appreciated. ; )