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.