  We investigate the viability of producing galaxy mock catalogues with
COmoving Lagrangian Acceleration (COLA) simulations in Modified Gravity (MG)
models employing the Halo Occupation Distribution (HOD) formalism. In this
work, we focus on two theories of MG: $f(R)$ gravity with the chameleon
mechanism, and a braneworld model (nDGP) that incorporates the Vainshtein
mechanism. We use a suite of full $N$-body simulations in MG as a benchmark to
test the accuracy of COLA simulations. At the level of Dark Matter (DM), we
show that COLA accurately reproduces the matter power spectrum up to $k \sim 1
h {\rm Mpc}^{-1}$, while it is less accurate in reproducing the velocity field.
To produce halo catalogues, we find that the ROCKSTAR halo-finder does not
perform well with COLA simulations. On the other hand, using a simple
Friends-of-Friends (FoF) finder and an empirical mass conversion from FoF to
spherical over-density masses, we are able to produce halo catalogues in COLA
that are in good agreement with those in $N$-body simulations. To consider the
effects of the MG fifth force on the halo profile, we derive simple fitting
formulae for the concentration-mass and the velocity dispersion-mass relations
that we calibrate using ROCKSTAR halo catalogues in $N$-body simulations. We
then use these results to extend the HOD formalism to modified gravity
simulations in COLA. We use an HOD model with five parameters that we tune to
obtain galaxy catalogues in redshift space. We find that despite the great
freedom of the HOD model, MG leaves characteristic imprints in the redshift
space power spectrum multipoles and these features are well captured by the
COLA galaxy catalogues.
