  We expand iterative numerically-exact influence functional path-integral
tools and present a method capable of following the nonequilibrium time
evolution of subsystems coupled to multiple bosonic and fermionic reservoirs
simultaneously. Using this method, we study the real-time dynamics of charge
transfer and vibrational mode excitation in an electron conducting molecular
junction. We focus on nonequilibrium vibrational effects, particularly, the
development of vibrational instability in a current-rectifying junction. Our
simulations are performed by assuming large molecular vibrational anharmonicity
(or low temperature). This allows us to truncate the molecular vibrational mode
to include only a two-state system. Exact numerical results are compared to
perturbative Master equation calculations demonstrating an excellent agreement
in the weak electron-phonon coupling regime. Significant deviations take place
only at strong coupling. Our simulations allow us to quantify the contribution
of different transport mechanisms, coherent dynamics and inelastic transport,
in the overall charge current. This is done by studying two model variants: The
first admits inelastic electron transmission only, while the second one allows
for both coherent and incoherent pathways.
