  We present a retuning of the IllustrisTNG baryonic physics model which can be
used to run large-box realistic cosmological simulations with a lower
resolution. This new model employs a lowered gas density threshold for star
formation and reduced energy releases by stellar and black hole feedback. These
changes ensure that our simulations can produce sufficient star formation to
closely match the observed stellar and gas properties of galaxies and galaxy
clusters, despite having $\sim$160 times lower mass resolution than the
simulations used to tune the fiducial IllustrisTNG model. Using the retuned
model, we have simulated Hu-Sawicki $f(R)$ gravity within a $301.75h^{-1}{\rm
Mpc}$ box. This is, to date, the largest simulation that incorporates both
screened modified gravity and full baryonic physics, offering a large sample
($\sim$500) of galaxy clusters and $\sim$8000 galaxy groups. We have reanalysed
the effects of the $f(R)$ fifth force on the scaling relations between the
cluster mass and four observable proxies: the mass-weighted gas temperature,
the Compton $Y$-parameter of the thermal Sunyaev-Zel'dovich effect, the X-ray
analogue of the $Y$-parameter, and the X-ray luminosity. We show that a set of
mappings between the $f(R)$ scaling relations and their $\Lambda$CDM
counterpart, which have been tested in a previous work using a much smaller
cosmological volume, are accurate to within a few percent for the
$Y$-parameters and $\lesssim7\%$ for the gas temperature for cluster-sized
haloes ($10^{14}M_{\odot}\lesssim M_{500}\lesssim10^{15}M_{\odot}$). These
mappings will be important for unbiased constraints of gravity using the data
from ongoing and upcoming cluster surveys.
