Low variance energy estimators for systems of quantum Drude oscillators: treating harmonic path integrals with large separations of time scales

In the effort to develop atomistic models capable of accurately describing nanoscale systems with complex interfaces, it has become clear that simple treatments with rigid charge distributions and dispersion coefficients selected to generate bulk properties are insufficient to predict important physical properties. The quantum Drude oscillator model, a system of one-electron pseudoatoms whose "pseudoelectrons" are harmonically bound to their respective "pseudonuclei," is capable of treating many-body polarization and dispersion interactions in molecular systems on an equal footing due to the ability of the pseudoatoms to mimic the long-range interactions that characterize real materials. Using imaginary time path integration, the Drude oscillator model can, in principle, be solved in computer operation counts that scale linearly with the number of atoms in the system. In practice, however, standard expressions for the energy and pressure, including the commonly used virial estimator, have extremely large variances that require untenably long simulation times to generate converged averages. In this paper, low-variance estimators for the internal energy are derived, in which the large zero-point energy of the oscillators does not contribute to the variance. The new estimators are applicable to any system of harmonic oscillators coupled to one another (or to the environment) via an arbitrary set of anharmonic interactions. The variance of the new estimators is found to be much smaller than standard estimators in three example problems, a one-dimensional anharmonic oscillator and quantum Drude models of the xenon dimer and solid (fcc) xenon, respectively, yielding 2-3 orders of magnitude improvement in computational efficiency.