First of all, instead of your $m(p_i)$ let's talk about any type of charge, denoted by $q_i$, because that's shorter and we physicists are lazy and also because we don't use the concept of negative mass too, eh.. lightly.
Secondly, the $|d(p_i,p_j)|^2$ in your expression seems to be wrong, if $d$ means distance then the expression for energy should have $|d(p_i,p_j)|$ without the square in the denominator.
So with those changes, we have to show that the energy
$$
E(x_1, x_2, \dots, x_n) = \sum_{\underset{\large i,j\leq n}{\small i\neq j}} \frac{q_i q_j}{|x_i-x_j|}
$$
has no stable configurations, i.e. no local minima, for variation of the positions $x_i$.
Actually the [link] you gave does not clearly state what variations are meant. It should be obvious that simultaneuosly scaling all $x_i \rightarrow \lambda x_i$ will decrease $E$ with $1/\lambda$ so that would immediately settle it.
But we can go further. Also in situations where some of the positions are constrained it is true. Even if $n-1$ positions are fixed and only one, say $x_n$, can move, $E$ still cannot have a minimum, because we can separate the terms involving $q_n$ in the sum:
$$
\sum_{\underset{\large i,j\leq n}{\small i\neq j}} \frac{q_i q_j}{|x_i-x_j|}
=\sum_{\underset{\large i,j\leq n-1}{\small i\neq j}} \frac{q_i q_j}{|x_i-x_j|}
+ q_n\ \underbrace{ \sum_{i=1}^{n-1} \frac{q_i}{|x_i-x_n|} }_{f(x_n)}
$$
so we only need to prove that the function $f(x_n)$ can have no minimum. Since all its terms $\frac{q_i}{|x_i-x_n|}$ are harmonic functions of $x_n$, i.e. solutions of the Laplace equation, the whole $f(x_n)$ must be a harmonic function, which proves it.
(For non-mathematicians: if a function $f(v)$ of a vector $v$ has a local minimum it must have a minimum in the $x,y$ and $z$ direction, so $\frac {d^2f} {dv_x^2}, \frac {d^2f} {dv_y^2}, \frac {d^2f} {dv_z^2}$ must all be positive, which implies $\nabla^2 f >0$, which would contradict $f$ being harmonic with $\nabla^2 f =0$.)
