  The thermal lines method for the evaluation of vibrational expectation values
of electronic observables [B. Monserrat, Phys. Rev. B 93, 014302 (2016)] was
recently proposed as a physically motivated approximation offering balance
between the accuracy of direct Monte Carlo integration and the low
computational cost of using local quadratic approximations. In this paper we
reformulate thermal lines as a stochastic implementation of quadrature grid
integration, analyze the analytical form of its bias, and extend the method to
multiple point quadrature grids applicable to any factorizable harmonic or
anharmonic nuclear wave function. The bias incurred by thermal lines is found
to depend on the local form of the expectation value, and we demonstrate that
the use of finer quadrature grids along selected modes can eliminate this bias,
while still offering a ~30% lower computational cost than direct Monte Carlo
integration in our tests.
