In this paper we consider a general class of second order stochastic partial differential equations on $\mathbb{R}^d$ driven by a Gaussian noise which is white in time and it has a homogeneous spatial covariance. Using the techniques of Malliavin calculus we derive the smoothness of the density of the solution at a fixed number of points $(t,x_1), \dots, (t,x_n)$, $t>0$, assuming some suitable regularity and non degeneracy assumptions. We also prove that the density is strictly positive in the interior of the support of the law.