Molecular systems are quantum systems, but the complete characterization of molecular motions within a fully quantum framework might appear to be an unfeasible task because it would require that the actual nuclear positions are established at any time. One would like to use a quantum molecular trajectory that defines the instantaneous nuclear positions and satisfies the predictions of quantum mechanics in terms of its statistical properties. Even though it can be proven that the single Bohm trajectory provides a representation of the quantum molecular trajectory, this solves the issue only on a theoretical ground: exact solutions of the Schrödinger-Bohm dynamical system are extremely computationally demanding. Therefore, we derive a stochastic equation of Smoluchowski type from the Schrödinger-Bohm dynamics, through projection operator techniques, in order to characterize the molecular motions of open quantum systems. The main quantum features of the motions emerge from the equilibrium distribution, i.e., the wave function's squared modulus integrated on the environment degrees of freedom. Furthermore, we verify the accuracy of the stochastic equation by comparing its predictions with those of the deterministic dynamics for a model system of six interacting harmonic oscillators. The indisputable advantage of this full quantum mechanical approach is that of representing the molecular dynamics, which controls important phenomena like vibrational relaxation, conformational transitions and activated processes, in a self consistent way and at the low computational cost of solving simple stochastic equations.