I think the mistake you are making is confusing a path and its image.
To make myself clearer, let's have a look at the proof that a minimizer of the energy is a geodesic. Let $\gamma : [a,b] \longrightarrow S $ be a smooth path
$$ \int{||\gamma'||} \leq (b-a)^2\int{||\gamma||^2} $$ and equality holds iff $||\gamma'||$ is a constant function.We can suppose $b-a=1$ to have $L(\gamma) \leq E(\gamma)$ without loss of generality. Since making a reparametrisation of $\gamma$ doesn't change the value of the lenght(exercise), take a minimizer of the lenght $\gamma$ and let $\tilde{\gamma}$ its reparametrisation at constant speed $L(\gamma)$
$$L(\gamma) = L(\tilde{\gamma}) \leq E(\gamma) $$
since $\gamma$ minimizes $L(\tilde{\gamma}) = E (\tilde{\gamma}) \geq E(\gamma)$
then $L(\gamma) = E(\tilde{\gamma})$ and this proves that a minimizer of the lenght is a minimizer of the energy iff and only if it is parametrised at constant speed ! Especially energy is not invariant by reparametrisation
But the parametrisation by $\theta$ is not at constant speed, and I think that's were you are wrong
You appear dazed and confused, so I will address some of your confusion.
Coordinate bases vs. orthonormal frames:
First of all, I don't know what you mean by $\hat{r}$ and $\hat{\theta}$, but in standard notation, this refers to unit vectors.
Generally speaking, coordinate bases and general frames behave differently.
If you take the plane $\mathbb{R}^2$ with standard cartesian coordinates $(x,y)$, and you make the transforms $$ x(r,\theta)=r\cos\theta \\ y(r,\theta)=r\sin\theta, $$ then the coordinate bases are given by $$ \mathbf{g}_1=\frac{\partial\mathbf{r}}{\partial r}\ \text{and}\ \mathbf{g}_2=\frac{\partial\mathbf{r}}{\partial\theta}. $$
Here $\mathbf{r}=(x,y)=x\hat{e}_1+y\hat{e}_2$, where $\hat{e}_i$ are the cartesian basis vectors.
In differential geometry notation you might see them as $\mathbf{g}_1=\partial/\partial r$ and $\mathbf{g}_2=\partial/\partial\theta$, because on a differentiable manifold, there are no "position vectors" like $\mathbf{r}$, so vectors have to be identified with differential operators.
Now, in vector calculus, it is common to use an orthonormal frame instead of the coordinate basis vector fields. The coordinate basis vectors for the polar coordinates are already orthogonal, and the $r$-vector is also unit length everywhere, but the $\theta$-vector is longer the farther you are from the origin, so we usually norm them. We have $\hat{r}=\partial\mathbf{r}/\partial r$ and $\hat{\theta}=(1/r)\partial\mathbf{r}/\partial\theta$, then these vectors form an orthonormal basis everywhere they are defined.
The components of the metric tensor are defined by the inner products of the basis vector fields (the metric tensor is abstractly the same as the inner product, but in terms of components, this is how you define it), so the metric tensor is $g_{ij}=\delta_{ij}$ in the orthonormal frame, but for polar coordinates, it is given by $$ (g_{ij})=\left(\begin{matrix}1 && 0 \\ 0 && r^2\end{matrix}\right). $$ So if you wanna calculate inner products via this form of the metric, you need to use coordinate bases, not the orthonormal frame.
Position vectors?!
What usually neither physicist literature or vector calculus literature will tell you in general, is that in vector calculus/multivariable calculus/physics, there are usually three kinds of vectors. These are:
Position/separation vectors: Arrows that point between two points. Position vectors point from the origin to another point. Separation vectors between arbitrary points.
"Tangent vectors": These vectors are firmly located in one point. They have finite length (norm), but they don't point between two points at all. If you think about it, it makes sense. Consider an electric field $\mathbf{E}(\mathbf{r})$. This associates to a point (represented by the position vector $\mathbf{r}$) a vector $\mathbf{E}(\mathbf{r})$, representing the strength of the electric field there. It is a vector because the electric field has a direction too. But ask, "to which point does that vector point to, from $\mathbf{r}$?". The question doesn't make sense. The vector simply sits in the point $\mathbf{r}$, and it pertains to that point and that point only.
"Moment" vectors: These vectors are not located anywhere. Stuff like, angular momentum of a rigid body, or total momentum of a system of multiple particles. It is a quantity with direction and magnitude, but the "location" of this quantity is ill-defined.
Now comes the important point. In differential geometry, only tangent vectors make sense, because in differential geometry, you consider spaces more general than vector/affine spaces, and the first and third vector concepts are only meaningful in a vector/affine space.
But even if you stay in euclidean spaces, but use curved coordinates to describe stuff, you have to, well, essentially use differential geometry as well.
Why is it? A coordinate system on, say, $\mathbb{R}^2$ is a map $$ \psi:\mathcal{U}\rightarrow \mathbb{R}^2, $$ where $\mathcal{U}$ is some open set of $\mathbb{R}^2$. For example the polar coordinate map is $$ \psi(x,y)=\left(\sqrt{x^2+y^2},\text{arctan}\left(\frac{y}{x}\right)\right), $$ although one usually gives the inverse of this map.
Now, if $\mathbf{r}_1=(x_1,y_1)$ is a point, and $\mathbf{r}_2=(x_2,y_2)$ is another one, then $\mathbf{r}=\alpha\mathbf{r}_1+\beta\mathbf{r}_2$ is also a point, given by the vector sum. What happens when we plug this into the coordinate change map?
$$ \psi(\mathbf{r})=\left(\sqrt{(\alpha x_1+\beta x_2)^2+(\alpha y_1+\beta y_2)^2},\text{arctan}\frac{\alpha y_1+\beta y_2}{\alpha x_1+\beta x_2}\right). $$
This is certainly not equal to $\alpha\psi(\mathbf{r}_1)+\beta\psi(\mathbf{r}_2)$, which is not surprising, since $\psi$ is not linear. But because the coordinate change map is not linear, vector addition and scalar multiplication gets messed up in the new coordinates.
This is why, if you use curvy coordinates, you cannot use "position vectors", and the distinction between points and coordinates must be kept clear.
In the standard coordinates of $\mathbb{R}^2$, or in any other coordinate system, which is related to the standard coordinates via a linear transformation, the addition and scalar multiples of points as vectors makes sense and is consistent (precisely because linear maps respect the vector space structure), so you can use them. But not in curvy coordinates, sorry.
Best Answer
To calculate the distance between two points, first you have to choose a path between those two, say $\gamma (s)$. If you want the shortest distance, as you said, solve the geodesic equation $\nabla_{\dot{\gamma}} \dot{\gamma} = 0$ to obtain function $\phi = \phi(\theta)$ and then calculate \begin{align} \int_{s_1}^{s_2} ds &= \int_{s_1}^{s_2} \sqrt{g_{\mu \nu}dx^{\mu} dx^{\nu}} \\ &=\int_{s_1}^{s_2} \sqrt{R^2 d\theta^2 + R^2 \sin^2 \theta d\phi^2} \\ &= R \int_{\theta_1}^{\theta_2} \sqrt{1+\sin^2 \theta \Big(\frac{d\phi}{d\theta}\Big)^2} d\theta \end{align} However, to show that the great circle is the geodesic on $\mathbb{S}^2$ you don't need to calculate the distance. Just solve geodesic eq.