  We implement an adaptation of the COLA approach, a hybrid scheme that
combines Lagrangian perturbation theory with an N-body approach, to model
non-linear collapse in chameleon and symmetron modified gravity models.
Gravitational screening is modeled effectively through the attachment of a
suppression factor to the linearized Klein-Gordon equations.
  The adapted COLA approach is benchmarked, with respect to an N-body code both
for the $\Lambda$CDM scenario and for the modified gravity theories. It is
found to perform well in the estimation of the dark matter power spectra, with
consistency of 1 % to $k\sim2.5$ h/Mpc. Redshift space distortions are shown to
be effectively modeled through a Lorentzian parameterization with a velocity
dispersion fit to the data. We find that COLA performs less well in predicting
the halo mass functions, but has consistency, within $1\sigma$ uncertainties of
our simulations, in the relative changes to the mass function induced by the
modified gravity models relative to $\Lambda$CDM.
  The results demonstrate that COLA, proposed to enable accurate and efficient,
non-linear predictions for $\Lambda$CDM, can be effectively applied to a wider
set of cosmological scenarios, with intriguing properties, for which clustering
behavior needs to be understood for upcoming surveys such as LSST, DESI, Euclid
and WFIRST.
