Abstract We present the UNIT N-body cosmological simulations project, designed to provide precise predictions for non-linear statistics of the galaxy distribution. We focus on characterizing statistics relevant to emission line and luminous red galaxies in the current and upcoming generation of galaxy surveys. We use a suite of precise particle mesh simulations (fastpm) as well as with full N-body calculations with a mass resolution of ${\sim } 1.2\times 10^9\, h^{-1}$M⊙ to investigate the recently suggested technique of Angulo and Pontzen to suppress the variance of cosmological simulations. We study redshift-space distortions, cosmic voids, higher order statistics from z = 2 down to 0. We find that both two- and three-point statistics are unbiased. Over the scales of interest for baryon acoustic oscillations and redshift-space distortions, we find that the variance is greatly reduced in the two-point statistics and in the cross-correlation between haloes and cosmic voids, but is not reduced significantly for the three-point statistics. We demonstrate that the accuracy of the two-point correlation function for a galaxy survey with effective volume of 20 (h−1Gpc)3 is improved by about a factor of 40, indicating that two pairs of simulations with a volume of 1 (h−1Gpc)3 lead to the equivalent variance of ∼150 such simulations. The N-body simulations presented here thus provide an effective survey volume of about seven times the effective survey volume of Dark Energy Spectroscopic Instrument or Euclid. The data from this project, including dark matter fields, halo catalogues, and their clustering statistics, are publicly available.