Classical molecular dynamics simulations have been the preferred method to cope with the characteristic sizes and time scales of complex life-science systems. However, while classical methods have well known limitations, such as that their accuracy strongly depends on empirical tuning, the practical use of far more accurate methods that rely on quantum Hamiltonians, has been limited by the current efficiency and scalability of sparse matrix-matrix multiplication algorithms used in the self-consistent field equations. In this paper, we show unprecedented massive scalability of a recently presented method, called MPSM3, for sparse matrix-matrix multiplication. The algorithmic basis of the method was presented in a recent publication, while here we describe the algorithmic enhancements that allow us to claim at least one order of magnitude improvement in scalability and time to solution over the state of the art (original MPSM3). We achieve a time to solution for the multiplication of density matrices within the self-consistent field scheme that is approaching the time needed to evaluate energy and forces with classical force-field methods and that is independent from the system size, provided proportional resources. This latest development renders the application of entirely quantum Hamiltonians to systems of several millions of atoms for extended molecular dynamics investigations feasible.