  We introduce a new path integral Monte Carlo method for investigating
nonadiabatic systems in thermal equilibrium and demonstrate an approach to
reducing stochastic error. We derive a general path integral expression for the
partition function in a product basis of continuous nuclear and discrete
electronic degrees of freedom without the use of any mapping schemes. We
separate our Hamiltonian into a harmonic portion and a coupling portion; the
partition function can then be calculated as the product of a Monte Carlo
estimator (of the coupling contribution to the partition function) and a
normalization factor (that is evaluated analytically). A Gaussian mixture model
is used to evaluate the Monte Carlo estimator in a computationally efficient
manner. Using two model systems, we demonstrate our approach to reduce the
stochastic error associated with the Monte Carlo estimator. We show that the
selection of the harmonic oscillators comprising the sampling distribution
directly affects the efficiency of the method. Our results demonstrate that our
path integral Monte Carlo method's deviation from exact Trotter calculations is
dominated by the choice of the sampling distribution. By improving the sampling
distribution, we can drastically reduce the stochastic error leading to lower
computational cost.
