We have developed an application and implemented parallel algorithms in order to provide a computational framework suitable for massively parallel supercomputers to study the unitary dynamics of quantum systems. We use renowned parallel libraries such as PETSc/SLEPc combined with high-performance computing approaches in order to overcome the large memory requirements to be able to study systems whose Hilbert space dimension comprises over 9 billion independent quantum states. Moreover, we provide descriptions on the parallel approach used for the three most important stages of the simulation: handling the Hilbert subspace basis, constructing a matrix representation for a generic Hamiltonian operator and the time evolution of the system by means of the Krylov subspace methods. We employ our setup to study the evolution of quasidisordered and clean many-body systems, focussing on the return probability and related dynamical exponents: the large system sizes accessible provide novel insights into their thermalization properties.