It looks like for small $n$, your assertion is true. For large $n$, there are symmetrically placed pairs on the side achieving the maximum. Please see this paper for details. (I do not have access to this paper yet.)

[Previuos comments]
The result is clear for $n=1$ from the AM-GM inequality. The case $n=2$ can be shown by complex analysis: Represent the regular triangle by the roots of unity $1,\zeta,\zeta^2$, where $\zeta =-\frac 12+\frac {\sqrt{3}}2i$. Then the required function is represented by $$p(z)=|(z-1)(z-\zeta)(z-\zeta)|=|z^3-1|.$$ By the maximum modulus principle and by symmetry we may assume that $z$ lies on the side joining $\zeta$ to $\zeta^2$, i.e. $z$ is of the form $z=-\frac 12+ti,$ where $-\frac{\sqrt 3}2\leq t\leq\frac{\sqrt 3}2$. By calculus, the maximum happens when $t=0$ and $z=-1/2$, which is the midpoint of a side of the regular triangle.

EDIT. Here is a suggestion for the general case. Taking logarithm, the problem is equivalent to finding the maximum of the function $$f=\frac 12\sum_{i=0}^n \log [(x_1-a_1^i)^2+\cdots+(x_n-a_n^i)^2],$$ where $(x_1,\cdots,x_n)$ is the coordinates for $P$, and $(a_1^i,\cdots,a_n^i)$ represents the vertices with $i=0,\cdots,n$. If $n>2$, then $\Delta f>0$, so $f$ is a subharmonic function and the maximum principle applies. This does not solve the problem yet, but it follows that the maximum occurs on the boundary faces.