Here is the explanation of the connection between positive definite functions and energy minimization on the sphere, that I promised earlier, as well as an alternative solution to the posted question, which extends to all dimensions. I'll try to keep preliminaries to a minimum.
Let $F:[-1,1]\rightarrow \mathbb R$ be a continuous function. We say that the function $F$ is positive definite on the sphere $S^{d-1}$ if for any finite subset $Z = \{ z_1, \dots, z_N\} \subset S^{d-1}$ the matrix $\big[ F ( \langle z_i, z_j \rangle ) \big]_{i,j=1}^N$ is positive semidefinite.
Let $\{P_n^d\}_{n=0}^\infty$ denote the sequence of Gegenbauer (ultraspherical) polynomials, i.e. polynomials on $[-1,1]$, orthogonal with respect to the weight $(1-t^2)^{\frac{d-3}2}$ (there are various normalizations and notations in the literature -- this is not the most common one, but it is simpler for our purposes). Notice that the weight is the one that arises when integrating a zonal function on $S^{d-1}$, i.e.
\begin{equation}\label{e.int}
\int\limits_{S^{d-1}} f ( \langle p, x \rangle ) d\sigma (x) = \frac{\omega_{d-2}}{\omega_{d-1}} \int_{-1}^1 f (t ) \big(1-t^2\big)^{\frac{d-3}2}\, dt,
\end{equation}
where $\sigma$ is the normalized uniform surface measure on the sphere $S^{d-1}$, $p \in S^{d-1}$ is an arbitrary pole, and $\omega_{d-1}$ is the Lebesgue surface measure of $S^{d-1}$.
When $d=2$ these are Chebyshev polynomials (of the first kind), and for $d=3$ these are the Legendre polynomials. Since this is an orthogonal sequence, for any function $F$ as above, one can consider its Gegenbauer expansion $\sum\limits_{n=0}^\infty \widetilde{F}_n P_n^d (t)$. We take $P_0^d (t) \equiv 1$, and one can see from the first equation, that
\begin{equation}
\widetilde{F}_0 = \int\limits_{S^{d-1}} F ( \langle p, x \rangle ) d\sigma (x) = \int\limits_{S^{d-1}} \int\limits_{S^{d-1}} F ( \langle x, y \rangle ) d\sigma (x) d\sigma (y),
\end{equation}
where the second identity follows by rotational invariance (the expression in the middle is independent of $p$). The aforementioned theorem of Schoenberg states that $F$ is positive definite on the sphere if and only if $\widetilde{F}_n \ge 0$ for all $n \ge 0$.
In some special cases there is an easier way to check for positive definiteness. Assume $F$ is also analytic in $(-1,1)$ and in its Taylor series has non-negative coefficients, then $F$ is positive definite on $S^{d-1}$ for any $d\ge 2$. This can be proved in a couple of ways: either by integration by parts and using Rodrigues formula, or alternatively as follows. Schur's theorem says that if matrices $A$ and $B$ are positive semidefinite, then their Hadamard (elementwise) product $A\circ B = [ a_{ij} b_{ij} ]$ is also positive semidefinite. Since the Gram matrix $\big[ \langle z_i, z_j \rangle \big]_{i,j=1}^N$ is always positive semidefinite, by Schur's theorem so are its Hadamard powers $\big[ \langle z_i, z_j \rangle ^k \big]_{i,j=1}^N$. And since all coefficients in the Taylor series of $F$ are non-negative, then the matrix $\big[ F ( \langle z_i, z_j \rangle ) \big]_{i,j=1}^N$ is also positive-semidefinite.
Now define the energy
\begin{equation}
I_F (\mu) = \int\limits_{S^{d-1}} \int\limits_{S^{d-1}} F ( \langle x, y \rangle ) d\mu (x) d\mu (y),
\end{equation}
and assume we want to minimize it over the class $\mathcal M$ of probability measures on the sphere $S^{d-1}$.
Claim: if $F$ is positive definite on $S^{d-1}$, up to an additive constant, then the uniform distribution $\sigma$ is a minimizer of $I_F$ (this is sort of folklore).
Here is one proof. Assume that $F$ is positive definite on $S^{d-1}$, up to an additive constant, i.e. $\widetilde{F}_n \ge 0$ for all $n \ge 1$. Recall that $\widetilde{F}_0 = I_F (\sigma)$, as discussed above. We shall use the ``addition formula" for spherical harmonics. Let $M_{n,d}$ be the dimension of the space of spherical harmonics on $S^{d-1}$ of degree $n$, and let $\{ H_{n,k} \}_{k=1}^{M_{n,d}}$ be an orthonormal basis of this space (with respect to $\sigma$). Then for any $x$, $y \in S^{d-1}$
\begin{equation}
\sum_{k=1}^{M_{n,d}} H_{n,k} (x) H_{n,k} (y) = P_n^d ( \langle x, y \rangle),
\end{equation}
where $P_n^d$ is the $n^{th}$ Gegenbauer polynomial. (In $d=2$, i.e. on the circle, this is equivalent to the formula $\cos (nx ) \cos (ny) + \sin (ny) \sin (ny) = \cos n (x-y)$, since spherical harmonics can be associated with the Fourier basis, and Chebyshev polynomials satisfy $T_n (\cos \theta) = \cos n \theta$.) We than have, using the addition formula,
\begin{align*}
I_F (\mu ) & = \widetilde{F}_0 + \sum_{n\ge1} \widetilde{F}_n \int\limits_{S^{d-1}} \int\limits_{S^{d-1}} P_n^d ( \langle x, y \rangle ) d\mu (x) d\mu (y) \\
&= I_F (\sigma) + \sum_{n\ge1} \widetilde{F}_n \sum_{k=1}^{M_{n,d}} \bigg( \int\limits_{S^{d-1}} H_{n,k} (x) d\mu (x) \bigg)^2 \ge I_F (\sigma),
\end{align*}
i.e. $\sigma$ is a minimizer. (Caution: I have used Fubini above, i.e. interchanged integration and summation, without justifying it -- but in fact there is an interesting self-improving property here: for positive definite $F$, its Gegenbauer series converges absolutely, this could be proved either using Mercer's theorem from spectral theory or directly using properties of Gegenbauer polynomials.)
This prove is the most natural, but for completeness I include another, completely elementary, proof of the Claim. Recall that positive definiteness means that for all $Z = \{ z_1, \dots, z_N\} \subset S^{d-1}$ and all $c_1, \dots, c_N \in \mathbb R$ we have
\begin{equation*}
I_F \big( \sum_{i} c_i \delta_{z_i} \big) = \sum_{i,j} F ( \langle z_i, z_j \rangle ) c_i c_j \ge 0,
\end{equation*}
(it is easy to see that the expressions above are equal). Since linear combinations of Dirac delta masses are weak-$*$ dense in the space of Borel measures, this condition is equivalent to saying that $I_F (\nu) \ge 0$ for all signed Borel measures. Then observe the following, for any $\mu \in \mathcal M$
\begin{align*}
0 &\le I_F ( \mu - \sigma) = I_F (\mu) + I_F (\sigma) - 2 \int\limits_{S^{d-1}} \int\limits_{S^{d-1}} F ( \langle x, y \rangle ) d\sigma (x) d\mu (y)\\
& = I_F (\mu) + I_F (\sigma) - 2 \int\limits_{S^{d-1}} I_F (\sigma) d\mu (y) = I_F (\mu) - I_F (\sigma) ,
\end{align*}
hence $\sigma$ is a minimizer. In the last line I used the fact that the inner integral in $\sigma$ is independent of $y$ and equals $I_F (\sigma)$.
Finally, returning to our problem, our function is $| \sin \angle(x,y) | = \sqrt{1-\langle x,y \rangle^2}$, i.e. the underlying function $F(t) = (1-t^2)^{1/2}$. We can write the Taylor series
\begin{equation}
(1-t^2)^{1/2} = \sum_{n=0}^\infty (-1)^n {{1/2} \choose n} t^{2n},
\end{equation}
and it is easy to see that all the coefficients (except $n=0$) are non-positive. Since we are maximizing (instead of minimizing), it follows from the discussion above that $\sigma$ is a maximizer in all dimensions $d\ge 2$.