Gravity simulation using Numpy and Pygame
The positions of the bodies are modified in-place, while the loop is running. It means that in your universe the action is not equal to reaction. Your simulation will exhibit some nonrealistic behaviors, like drift of the center of the mass.
Euler integration method is not the most accurate. Consider Runge-Kutta. At least, monitor the motion invariants (total energy, momentum, and angular momentum of the ensemble), and do corrective actions when they start diverging.
After f = m[i] * ... and f_total / m[i] the mass cancels. Multiplication and division could be safely omitted.
Along the same line, scale the masses by G once, before the simulation begins.