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}$.

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

For any non-negative test function $\Psi\in\mathscr{D}(\text{B})$ we have \begin{align*} \int (\Phi\star f)\Delta \Psi = \int (\Delta \Phi \star f)\Psi \geq \int (\Delta \Phi \star f_\varepsilon)\Psi = \int (\Phi\star\Delta f_\varepsilon)\Psi, \end{align*} where we used that $\Delta \Phi\geq 0$ on $2\text{B}$. This inequality means $\Delta(\Phi\star f) \geq \Phi\star \Delta f_\varepsilon$ in $\mathscr{D}'(\text{B})$ and since both functions are smooth, we have the same inequality pointwise in $\text{B}$. Now replace $\Phi$ by $\Phi_\delta:=\Phi(x/\delta)/\delta^d$. By standard properties of convolution you have a decreasing family $\delta_\epsilon$ 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.

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.