  We investigate the exponential $f(Q)$ symmetric teleparallel gravitation,
namely $f(Q)=Q+\alpha Q_0(1-e^{-\beta\sqrt{Q/Q_0}})$ using \texttt{ME-GADGET}
code to probe the structure formation with box sizes $L_{\mathrm{box}}=10/100$
Mpc$/h$ and middle resolution $N_p^{1/3}=512$. To reproduce viable cosmology
within the aforementioned modified gravity theory, we first perform Markov
Chain Monte Carlo (MCMC) sampling on OHD/BAO/Pantheon datasets and constrain a
parameter space. Furthermore, we also derive theoretical values for
deceleration parameter $q(z)$, statefinder pair $\{r,s\}$ and effective
gravitational constant $G_{\mathrm{eff}}$, perform $Om(z)$ diagnostics. While
carrying out N-body+SPH simulations, we derive CDM+baryons over
density/temperature/mean molecular weight fields, matter power spectrum (both
2/3D, with/without redshift space distortions), bispectrum, two-point
correlation function and halo mass function. Results for small and big
simulation box sizes are therefore properly compared, halo mass function is
related to the Seth-Tormen theoretical prediction and matter power spectrum to
the standard \texttt{CAMB} output.
