Withtout loss of generality, we can assume $p=0$ and $f(0)=0$. 

Also, the problem is purely local : we can assume that $f$ and all the functions $f_\epsilon$ are compactly supported in the unit ball $\text{B}$, which contains also all the neighborhoods $U_\epsilon$. You can also assume that for all $\epsilon>0$ the estimate $f\leq f_\epsilon$ is satisfied on $\text{B}$.

Now pick a smooth function $\Phi$ subharmonic on $2 \text{B}$, integrable over $\mathbf{R}^d$ with total mass equal $1$.

Fix $\epsilon$ and take $f_\epsilon$ given by your statement. Since $f$ and $f_\epsilon$ vanish outside $\text{B}$, we claim $\Delta \Phi\star (f-f_\epsilon)\geq 0$ on $\text{B}$, because $\Delta \Phi\geq 0$ on $2\text{B}$. This means in particular that $\Delta(\Phi\star f) \geq \Phi\star \Delta f_\varepsilon$ is true at the origin. Now replace $\Phi$ by $\Phi_\delta:=\Phi(x/\delta)/\delta^d$. By standard properties of convolution you have a decreasing family $\delta_\epsilon\rightarrow 0$ such that $\Phi_{\delta_\epsilon} \star \Delta f_\epsilon(0)\geq-2\epsilon$.

The corresponding sequence $g_\epsilon:=\Phi_{\delta_{\epsilon/2}}\star f$ does the job.