A massively parallel, direct quantum molecular dynamics method is described. The method combines a quantum trajectory (QT) representation of the nuclear wave function discretized into an ensemble of trajectories with an electronic structure (ES) description of electrons, namely using the density functional tight binding (DFTB) theory. Quantum nuclear effects are included into the dynamics of the nuclei via quantum corrections to the classical forces. To reduce computational cost and increase numerical accuracy, the quantum corrections to dynamics resulting from localization of the nuclear wave function are computed approximately and included into selected degrees of freedom representing light particles where the quantum effects are expected to be the most pronounced. A massively parallel implementation, based on the message passing interface allows for efficient simulations of ensembles of thousands of trajectories at once. The QTES-DFTB dynamics approach is employed to study the role of quantum nuclear effects on the interaction of hydrogen with a model graphene sheet, revealing that neglect of nuclear effects can lead to an overestimation of adsorption.