Abstract
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.