As I mentioned in my comment, the tetrahedral formula is invariant under translations, so let's focus on regular tetrahedra conveniently centered at the origin.
Let $T$ be the coordinate matrix such a tetrahedron; that is, the matrix whose columns are coordinates in $\mathbb{R}^3$ of the tetrahedron's vertices. The columns of the matrix obviously sum to zero, but there's something less-obvious that we can say about the rows:
Fact: The rows of $T$ form an orthogonal set of vectors of equal magnitude, $m$.
For example (and proof-of-fact), take the tetrahedron that shares vertices with the double-unit cube, for which $m=2$:
$T = \begin{bmatrix}1&1&-1&-1\\1&-1&1&-1\\1&-1&-1&1\end{bmatrix} \hspace{0.25in}\text{so that}\hspace{0.25in} T T^\top=\begin{bmatrix}4&0&0\\0&4&0\\0&0&4\end{bmatrix}=m^2 I$
Any other origin-centered regular tetrahedron is similar to this one, so its coordinate matrix has the form $S = k Q T$ for some orthogonal matrix $Q$ and some scale factor $k$. Then
$SS^\top = (kQT)(kQT)^\top = k^2 Q T T^\top Q^\top = k^2 Q (m^2 I) Q^\top = k^2 m^2 (Q Q^\top) = k^2 m^2 I$
demonstrating that the rows of $S$ are also orthogonal and of equal magnitude. (Fact proven.)
For the general case, take $T$ as follows
$T=\begin{bmatrix}a_x&b_x&c_x&d_x\\a_y&b_y&c_y&d_y\\a_z&b_z&c_z&d_z\end{bmatrix}$
Now, consider the matrix $J := \left[1,i,0\right]$. Left-multiplying $T$ by $J$ gives $P$, the coordinate matrix (in $\mathbb{C}$) of the projection of the tetrahedron into the coordinate plane:
$P := J T = \left[a_x+i a_y, b_x+ib_y, c_x+i c_y, d_x + i d_y\right] = \left[a, b, c, d\right]$
where $a+b+c+d=0$. Observe that
$P P^\top = a^2 + b^2 + c^2 + d^2$
On the other hand,
$PP^\top = (JT)(JT)^\top = J T T^\top J^\top = m^2 J J^\top = m^2 (1 + i^2) = 0$
Therefore,
$(a+b+c+d)^2=0=4(a^2 + b^2 + c^2 + d^2)$
Note: It turns out that the Fact applies to all the Platonic solids ... and most Archimedeans ... and a great many other uniforms, including wildly self-intersecting realizations (even in many-dimensional space). The ones for which the Fact fails have slightly-deformed variants for which the Fact succeeds. (The key is that the coordinate matrices of these figures are (right-)eigenmatrices of the vertex adjacency matrix. That is, $TA=\lambda T$. For the regular tetrahedron, $\lambda=-1$; for the cube, $\lambda = 1$; for the great stellated dodecahedron, $\lambda=-\sqrt{5}$; for the small retrosnub icosicosidodecahedron, $\lambda\approx-2.980$ for a pseudo-classical variant whose pentagrammic faces have non-equilateral triangular neighbors.)
The argument of my answer works for all "Fact-compliant" origin-centered polyhedra, so that $(\sum p_i)^2 = 0 = \sum p_i^2$ for projected vertices $p_i$. Throwing in a coefficient --namely $n$, the number of vertices-- that guarantees translation-invariance, and we have
$\left( \sum p_i \right)^2 = n \sum p_i^2$