  Starting with the exact factorization of the molecular wavefunction, this
paper presents the results from the numerical implementation in nonadiabatic
molecular dynamics of the recently proposed bohmion method. Within the context
of quantum hydrodynamics, we introduce a regularized nuclear Bohm potential
admitting solutions comprising a train of $\delta$-functions which provide a
finite-dimensional sampling of the hydrodynamic flow paths. The bohmion method
inherits all the basic conservation laws from its underlying variational
structure and captures electronic decoherence. After reviewing the general
theory, the method is applied to the well-known Tully models, which are used
here as benchmark problems. In the present case of study, we show that the new
method accurately reproduces both electronic decoherence and nuclear population
dynamics.
