  Dirac-Frenkel variational method with Davydov D2 trial wavefunction is
extended by introducing a thermalization algorithm and applied to simulate
dynamics of a general open quantum system. The algorithm allows to control
temperature variations of a harmonic finite size bath, when in contact with the
quantum system. Thermalization of the bath vibrational modes is realised via
stochastic scatterings, implemented as a discrete-time Bernoulli process with
Poisson statistics. It controls bath temperature by steering vibrational modes'
evolution towards their canonical thermal equilibrium. Numerical analysis of
the exciton relaxation dynamics in a small molecular cluster reveals that
thermalization additionally provides significant calculation speed up due to
reduced number of vibrational modes needed to obtain the convergence.
