  Primordial non-Gaussianity (PNG) is one of the most powerful probes of the
early Universe and measurements of the large scale structure of the Universe
have the potential to transform our understanding of this area. However
relating measurements of the late time Universe to the primordial perturbations
is challenging due to the non-linear processes that govern the evolution of the
Universe. To help address this issue we release a large suite of N-body
simulations containing four types of PNG: \textsc{quijote-png}. These
simulations were designed to augment the \textsc{quijote} suite of simulations
that explored the impact of various cosmological parameters on large scale
structure observables. Using these simulations we investigate how much
information on PNG can be extracted by extending power spectrum and bispectrum
measurements beyond the perturbative regime at $z=0.0$. This is the first joint
analysis of the PNG and cosmological information content accessible with power
spectrum and bispectrum measurements of the non-linear scales. We find that the
constraining power improves significantly up to $k_\mathrm{max}\approx 0.3
h/{\rm Mpc}$, with diminishing returns beyond as the statistical probes
signal-to-noise ratios saturate. This saturation emphasizes the importance of
accurately modelling all the contributions to the covariance matrix. Further we
find that combining the two probes is a powerful method of breaking the
degeneracies with the $\Lambda$CDM parameters.
