We study existence and regularity of the density for the solution $u(t,x)$ (with fixed $t > 0$ and $x \in D$) of the heat equation in a bounded domain $D \subset \mathbb R^d$ driven by a stochastic inhomogeneous Neumann boundary condition with stochastic term. The stochastic perturbation is given by a fractional Brownian motion process. Under suitable regularity assumptions on the coefficients, by means of tools from the Malliavin calculus, we prove that the law of the solution has a smooth density with respect to the Lebesgue measure in $\mathbb R$.