Is there an easier way to show that a tetrahedron is optimal? Basically, I'm trying to address a problem that somewhat mirrors electron geometry, and I'm phrasing it like this:
Consider four points in three dimensional space ($P_1,P_2,P_3,P_4$) such that each $P_i$ is on the unit sphere. Define the disturbance of the four points as:
$$ D(P_1,P_2,P_3,P_4)=\sum_{1 \leq i < j \leq 4} \frac{1}{|P_i-P_j|^2}$$
Where $|P_i-P_j|$ is the euclidean distance between $P_i$ and $P_j$. My goal is to minimize the disturbance of the four points, and show that the four points which minimize the disturbance form a tetrahedron.
Since we are only considering the shape (up to rotation), we can assume w.l.o.g. that $P_1=(1,0,0)$, and that $P_2$ also lies on the plane $z=0$. To do this, I've decided to parametrize these points as:
$$ P_2=(\cos(\alpha),\sin(\alpha),0)$$
$$ P_3=(\sin(\phi_1)\cos(\theta_1),\sin(\phi_1)\sin(\theta_1),\cos(\phi_1))$$
$$ P_3=(\sin(\phi_2)\cos(\theta_2),\sin(\phi_2)\sin(\theta_2),\cos(\phi_2))$$
And in this sense, the disturbance is a function from $\mathbb{R}^5$ to $\mathbb{R}$. I took the partial derivatives with respect to each variable and set them equal to zero, but it just all gives an incredibly nasty system of equations that I don't know how to simplify. My goal is to show that they all must form a tetrahedron but I really don't know how to do that, and I am hoping there is an easier method.
 A: This is really a comment rather than an answer, but unfortunately I can't comment yet, so.. And I apologise if it is not helpful, I'm not a mathematician, but the question seems interesting!
This question seems reminiscent of the Coulomb problem physics and chemistry, where interaction between N bodies is considered as a sum of pairwise interactions between single bodies. 
For the basic form of pairwise interaction, the Coulomb kernel is
$$K(\textbf{r}_i,\textbf{r}_j) = \frac{1}{|\textbf{r}_i-\textbf{r}_j|} $$
In physics this is intimately related to the Coulomb potential, the derivative of which gives the inverse square term you've written, which gives forces.
I think the general way to get something workable in this kind of problem is to do a Taylor expansion in 3D, then express in terms of the exponential planewave operator, then in terms of solid harmonics for which rules sum exist and expressions are tabulated.
The standard result for the Coulomb kernel written in terms of spherical harmonics is
$$K(\textbf{r}_i,\textbf{r}_j) = \sum^\infty_{l=0}\sum^l_{m=-l}\frac{4\pi}{2l+1}Y_{l,m}(\mathbf{r}_1)Y^*_{l,m}(\mathbf{r}_2)\frac{{\mathrm{min}(r_i,r_j)}^{l+1}}{\mathrm{min}(r_i,r_j)^l}$$
where the Ylm expressions can be looked up in tables.
I imagine for four unit point charges you consider a density distribution of four Dirac deltas
$$\rho(\mathbf{r})=\sum^4_{i=1}\delta(\mathbf{r}-\mathbf{r}_i)$$
evaluate a Hartree-type interaction
$$H(\mathbf{r},\mathbf{r'})=\int\int\rho(\mathbf{r})K(\textbf{r}_i,\textbf{r}_j)\rho(\mathbf{r'})drdr'$$
Hopefully that simplifies down using the properties of integrated Dirac deltas, and you can get an expression which you can minimise by the tetrahedron points.
