For $a,x\in\mathbf{R}^n$, let $$f_a(x) = -\exp(-\|x-a\|^2)$$ Then a linear combination of $f_a$ will satisfy these conditions, where $a$ ranges over your chosen minima. E.g., if $n=2$, $$4f_{(0,0)}+ f_{(1,0)}+ f_{(0,1)}$$ will have a global minimum at the origin and two symmetric local minima.

For $\lambda\in\mathbf{R}$, $a,x\in\mathbf{R}^n$, let $$f_{\lambda;a}(x) = -\exp(-\lambda\|x-a\|^2)$$ Then linear combinations of $f_{\lambda;a}$ can satisfy these conditions, where $a$ ranges over your chosen minima. E.g., if $n=2$, $$f_{6;(0,0)}+ f_{1;(1,0)}+ f_{1;(0,1)}$$ will have a global minimum at the origin and two symmetric local minima, as in these plots and computations.