Potential function with a fixed set of minima and no poles. Can we find a continuous, almost-everywhere smooth, potential function with no poles, such that the set of local minima for the sum of these potentials around $n$ points are exactly those $n$ points?

We are given a set of $n$ points $(x_0, \ldots, x_{n-1})$ in $\mathbb{R}^d$. We consider the additive potential function
$$g(x) = \sum_{i=1}^n f(x, x_i)$$
as a sum of potentials $f : \mathbb{R}^d \times \mathbb{R}^d \mapsto \mathbb{R}$, where $f$ is continuous and almost everywhere smooth.
We would like $g$’s minima to be exactly the $x_i$. That is, any point $x$ where the gradient of $g$ vanishes is either one of the $x_i$, a saddle point, or a maximum.
$$\forall x, \nabla g(x) = \vec{0} \implies x \in \{x_1, \ldots, x_n\}~\vee~ \exists v, \lambda < 0, H(g)(x) v = \lambda v$$
By way of example, we can construct such a function with the gravitational potential. If we let $$f(x,y) = -\frac{1}{||x-y||},$$ I believe the property is respected. A point mass will always end up “falling” towards one of the attractors, or at least I believe that’s implied by Earnshaw’s theorem. I’m not 100% sure it works for all norms, but the Euclidian distance will do.
For this potential, the $x_i$ aren’t just local minima, they are poles. We would like to find a better-behaved function which does not exhibit poles.
One way to “cheat” is to cut-off the poles. Once we get close enough to one of the $x_i$, the influence of the other points is low enough that the potential can be approximated as just the potential coming from that $x_i$. This means that, for a given set of $x_i$, there will always be $\epsilon > 0$ s.t. $f(x,y) = -\frac{1}{\epsilon + ||x-y||}$ works as intended.
However, a third constraint is that we would like $f$ to not depend on the $x_i$. Are there such functions?
A famous non-example would be a standard multivariate Gaussian. If you sum three of those on the vertices of an equilateral triangle, counter-intuitively, the center is a stable equilibrium.
In general, I suspect that $f$ cannot be smooth at $f(x,x)$. Otherwise, $f$ would be locally quadratic and if three points in $x_i$ are sufficiently close to each other, we would expect their barycenter (w.r.t the norm induced by the Hessian of $f$) to be a minimum.
However, this does not rule out solutions where $f$ is non-differentiable around $f(x,x)$.
 A: Given how you phrased your problem, I assume that you want $f(x,y)$ to exclusively depend on $x-y$ and we can replace it by $\hat{f}(x-y)$. I strongly suggest thinking in terms of such a function as it makes everything easier and the alternative would be that the behaviour of your potential depends on absolute location.
$\hat{f}(x) := \sqrt{|x|}$ should fulfil your requirements.
It is continuous, smooth almost everywhere, and bounded.
It fulfils the necessary requirement of having a local minimum at $\vec{0}$.
As $\vec{0}$ is a pole of $∇\hat{f}$ (and $\hat{f}$ is smooth everywhere else), the local minimum at $\vec{0}$ will outcompete any smooth function added to it: $\hat{f}(x) + h(x)$ has a local minimum at $\vec{0}$ for any smooth and bounded $h$.
To show that there cannot be any further local minima, note that the divergence of the field induced by $\hat{f}$ is positive everywhere except for $\vec{0}$. Thus the divergence of the field induced by $g$ must be positive everywhere, except at any of the $x_i$. As the field induced by a smooth potential (here $g$) must have a negative divergence at a local minimum, none can exist (except where $g$ is not smooth).
