A new algorithm for efficient and fully time-reversible integration of first-principles molecular dynamics based on orbital-free density functional theory (OFDFT) is presented. The algorithm adapts to this nontrivial case, the recently introduced Mass-Zero (MaZe) constrained dynamics. The formalism ensures that full adiabatic separation is enforced between nuclear and electronic degrees of freedom and, consequently, that the exact Born-Oppenheimer probability for the nuclei is sampled. Numerical integration of the MaZe dynamics combines standard molecular dynamics algorithms, e.g., Verlet or velocity Verlet, with the SHAKE method to impose the minimum conditions on the electronic degrees of freedom as a set of constraints. The developments presented in this work, which include a bespoke adaptation of the standard SHAKE algorithm, ensure that the quasilinear scaling of OFDFT is preserved by the new method for a broad range of kinetic and exchange-correlation functionals, including nonlocal ones. The efficiency and accuracy of the approach are demonstrated via calculations of static and dynamic properties of liquid sodium in the constant energy and constant temperature ensembles.