the same for all particles and is allowed to vary adaptively in time. In particular, according to the stability criterion for the leap-frog time integration, we adopt $\Dt=\eta/\sqrt{\max{|\nabla^2
\phi|}}$, where $\eta\lsim0.3$ is a dimensionless parameter. We found that $\eta=0.1$ assures good conservation of the total energy in the Newtonian cases (see Section \[secic\]). In the present version of the code, all the computations on the particles and the particle-mesh interpolations can be split among different processors, while the computations relative to the potential solver are not performed in parallel. The solution of equation (\[eqMOND\]) over the grid is then the bottleneck of the simulations: however, the iterative procedure on which the potential solver is based (see CLN) allows to adopt as seed solution at each time step the potential previously determined.

The MOND potential solver can also solve the Poisson equation (obtained by imposing $\mu=1$ in equation \[eqMOND\]), so Newtonian simulations can be run with the same code. We exploited this property to test the code by running several Newtonian simulations of both equilibrium distributions and collapses, comparing the results with those of simulations (starting from the same initial conditions) performed with the FVFPS treecode (Londrillo, Nipoti & Ciotti 2003; 