  The reduced density matrix of excitons coupled to a phonon bath at a finite
temperature is studied using the path integral Monte Carlo method. Appropriate
choices of estimators and importance sampling schemes are crucial to the
performance of the Monte Carlo simulation. We show that by choosing the
population-normalized estimator for the reduced density matrix, an efficient
and physically-meaningful sampling function can be obtained. In addition, the
nonadiabatic phonon probability density is obtained as a byproduct during the
sampling procedure. For importance sampling, we adopted the Metropolis-adjusted
Langevin algorithm. The analytic expression for the gradient of the target
probability density function associated with the population-normalized
estimator cannot be obtained in closed form without a matrix power series. An
approximated gradient that can be efficiently calculated is explored to achieve
better computational scaling and efficiency. Application to a simple
one-dimensional model system from the previous literature confirms the
correctness of the method developed in this manuscript. The displaced harmonic
model system within the single exciton manifold shows the numerically exact
temperature dependence of the coherence and population of the excitonic system.
The sampling scheme can be applied to an arbitrary anharmonic environment, such
as multichromophoric systems embedded in the protein complex. The result of
this study is expected to stimulate further development of real time
propagation methods that satisfy the detailed balance condition for exciton
populations.
