We develop the halo model of large-scale structure as an accurate tool for probing primordial non-Gaussianity (PNG). In this study we focus on understanding the matter clustering at several redshifts in the context of PNG that is a quadratic correction to the local Gaussian potential, characterized by the parameter fNL. In our formulation of the halo model we pay special attention to the effect of halo exclusion, and show that this can potentially solve the long standing problem of excess power on large scales in this model. The halo model depends on the mass function, clustering of halo centers and the density profiles. We test these ingredients using a large ensemble of highresolution Gaussian and non-Gaussian numerical simulations, covering fNL = {0, +100, −100}. In particular, we provide a first exploration of how halo density profiles change in the presence of PNG. We find that for fNL positive/negative high mass haloes have an increased/decreased core density, so being more/less concentrated than in the Gaussian case. We also examine the halo bias and show that, if the halo model is correct, then there is a small asymmetry in the scale-dependence of the bias on very large scales, which arises because the Gaussian bias must be renormalized. We show that the matter power spectrum is modified by ∼ 2.5% and ∼ 3.5% on scales k ∼ 1.0 h Mpc −1 at z = 0.0 and z = 1.0, respectively. Our halo model calculation reproduces the absolute amplitude to within 10% and the ratio of non-Gaussian to Gaussian spectra to within 1%. We also measure the matter correlation function and find similarly good levels of agreement between the halo model and the data. We anticipate that this modeling will be useful for constraining fNL from measurements of the shear correlation function in future weak lensing surveys such as Euclid.