Our BD- or GR-specific simulation pipeline begins from a set of input parameters θ_BD or θ_GR. First, they are passed to the Einstein-Boltzmann solver HI_CLASS, which solves the background and perturbations and outputs the linear power spectrum. Next, FML back-scales P_{HI_CLASS}(k,z) from z = 0 to z = z_init, realizes initial conditions for an N-body simulation with the random seed s and evolves it forward with the COLA method on a uniform grid, outputting the nonlinear power spectrum. Finally, FML’s first snapshot can be passed on to RAMSES, which evolves the same initial conditions with the standard N-body method on an adaptive mesh and yields an independent nonlinear power spectrum.