A method to compute probability current and its surface integral, the total flux, for systems of many particles of different masses is presented, based on transforming the wave function and its gradient onto a mass-weighted coordinate system. As a test for this methodology, it has been applied to a nontrivial 6-dimensional quantum dynamics study of a model of the operation of the proton-wire in Green Fluorescent Protein [O. Vendrell, R. Gelabert, M. Moreno, and J. M. Lluch, J. Phys. Chem. B, 112, 5500-5511 (2008)]. An adaptive Monte Carlo method is proposed, with favorable scaling properties for future applications, to solve the flux integral. Comparison of total reactive flux with the time derivative of the survival probability is satisfactory, corroborating the adequacy of the derivation. Using the new method the flux can quantitatively be divided into its positive and negative contributions, or more relevantly, into tunneling and classical parts. © 2011 American Institute of Physics.