  A large fraction of cosmological information on dark energy and gravity is
encoded in the nonlinear regime. Precision cosmology thus requires precision
modeling of nonlinearities in general dark energy and modified gravity models.
We modify the Gadget-2 code and run a series of N-body simulations on modified
gravity cosmology to study the nonlinearities. The modified gravity model that
we investigate in the present paper is characterized by a single parameter
\zeta, which determines the enhancement of particle acceleration with respect
to general relativity (GR), given the identical mass distribution (\zeta = 1 in
GR). The first nonlinear statistics we investigate is the nonlinear matter
power spectrum at k < 3h/Mpc, which is the relevant range for robust weak
lensing power spectrum modeling at l < 2000. In this study, we focus on the
relative difference in the nonlinear power spectra at corresponding redshifts
where different gravity models have the same linear power spectra. This
particular statistics highlights the imprint of modified gravity in the
nonlinear regime and the importance to include the nonlinear regime in testing
GR. By design, it is less susceptible to the sample variance and numerical
artifacts. We adopt a mass assignment method based on wavelet to improve the
power spectrum measurement. We run a series of tests to determine the suitable
simulation specifications (particle number, box size and initial redshift). We
find that, the nonlinear power spectra can differ by ~30% for 10% deviation
from GR (|\zeta-1| = 0.1) where the rms density fluctuations reach 10. This
large difference, on one hand, shows the richness of information on gravity in
the corresponding scales, and on the other hand, invalidates simple
extrapolations of some existing fitting formulae to modified gravity cosmology.
