  A numerically exact Monte Carlo scheme for calculation of open quantum system
dynamics is proposed and implemented. The method consists of a Monte-Carlo
summation of a perturbation expansion in terms of trajectories in Liouville
phase-space with respect to the coupling between the excited states of the
molecule. The trajectories are weighted by a complex decoherence factor based
on the second-order cumulant expansion of the environmental evolution. The
method can be used with an arbitrary environment characterized by a general
correlation function and arbitrary coupling strength. It is formally exact for
harmonic environments, and it can be used with arbitrary temperature. Time
evolution of an optically excited Frenkel exciton dimer representing a
molecular exciton interacting with a charge transfer state is calculated by the
proposed method. We calculate the evolution of the optical coherence elements
of the density matrix and linear absorption spectrum, and compare them with the
predictions of standard simulation methods.
