  Next generation galaxy surveys demand the development of massive ensembles of
galaxy mocks to model the observables and their covariances, what is
computationally prohibitive using $N$-body simulations. COLA is a novel method
designed to make this feasible by following an approximate dynamics but with up
to 3 orders of magnitude speed-ups when compared to an exact $N$-body. In this
paper we investigate the optimization of the code parameters in the compromise
between computational cost and recovered accuracy in observables such as
two-point clustering and halo abundance. We benchmark those observables with a
state-of-the-art $N$-body run, the MICE Grand Challenge simulation (MICE-GC).
We find that using 40 time steps linearly spaced since $z_i \sim 20$, and a
force mesh resolution three times finer than that of the number of particles,
yields a matter power spectrum within $1\%$ for $k \lesssim 1\,h {\rm
Mpc}^{-1}$ and a halo mass function within $5\%$ of those in the $N$-body. In
turn the halo bias is accurate within $2\%$ for $k \lesssim 0.7\,h {\rm
Mpc}^{-1}$ whereas, in redshift space, the halo monopole and quadrupole are
within $4\%$ for $k \lesssim 0.4\,h {\rm Mpc}^{-1}$. These results hold for a
broad range in redshift ($0 < z < 1$) and for all halo mass bins investigated
($M > 10^{12.5} \, h^{-1} \, {\rm M_{\odot}}$). To bring accuracy in clustering
to one percent level we study various methods that re-calibrate halo masses
and/or velocities. We thus propose an optimized choice of COLA code parameters
as a powerful tool to optimally exploit future galaxy surveys.
