ABSTRACT We introduce algorithms for black hole physics, i.e. black hole formation, accretion, and feedback, into the Numerical Investigation of a Hundred Astrophysical Objects (NIHAO) project of galaxy simulations. This enables us to study high mass, elliptical galaxies, where feedback from the central black hole is generally thought to have a significant effect on their evolution. We furthermore extend the NIHAO suite by 45 simulations that encompass z = 0 halo masses from 1 × 1012 to $4 \times 10^{13}\, \mathrm{M}_{⊙ }$, and resimulate five galaxies from the original NIHAO sample with black hole physics, which have z = 0 halo masses from 8 × 1011 to $3 \times 10^{12}\, \mathrm{M}_{⊙ }$. Now NIHAO contains 144 different galaxies and thus has the largest sample of zoom-in simulations of galaxies, spanning z = 0 halo masses from 9 × 108 to $4 \times 10^{13}\, \mathrm{M}_{⊙ }$. In this paper we focus on testing the algorithms and calibrating their free parameters against the stellar mass versus halo mass relation and the black hole mass versus stellar mass relation. We also investigate the scatter of these relations, which we find is a decreasing function with time and thus in agreement with observations. For our fiducial choice of parameters we successfully quench star formation in objects above a z = 0 halo mass of $10^{12}\, \mathrm{M}_{⊙ }$, thus transforming them into red and dead galaxies.