  In a previous paper we introduced a new method for simulating collisional
gravitational $N$-body systems with linear time scaling on $N$, based on the
Multi-Particle Collision (MPC) approach. This allows us to simulate globular
clusters with a realistic number of stellar particles in a matter of hours on a
typical workstation. We evolve star clusters containing up to $10^6$ stars to
core collapse and beyond. We quantify several aspects of core collapse over
multiple realizations and different parameters, while always resolving the
cluster core with a realistic number of particles. We run a large set of N-body
simulations with our new code. The cluster mass function is a power-law with no
stellar evolution, allowing us to clearly measure the effects of the mass
spectrum. Leading up to core collapse, we find a power-law relation between the
size of the core and the time left to core collapse. Our simulations thus
confirm the theoretical self-similar contraction picture but with a dependence
on the slope of the mass function. The time of core collapse has a
non-monotonic dependence on the slope, which is well fit by a parabola. This
holds also for the depth of core collapse and for the dynamical friction
timescale of heavy particles. Cluster density profiles at core collapse show a
broken power law structure, suggesting that central cusps are a genuine feature
of collapsed cores. The core bounces back after collapse, and the inner density
slope evolves to an asymptotic value. The presence of an intermediate-mass
black hole inhibits core collapse. We confirm and expand on several predictions
of star cluster evolution before, during, and after core collapse. Such
predictions were based on theoretical calculations or small-size direct
$N$-body simulations. Here we put them to the test on MPC simulations with a
much larger number of particles, allowing us to resolve the collapsing core.
