The calculation of thermochemical data requires accurate molecular energies. When such high accuracy is needed, often a method such as G1, G2, or G3 is used to calculate the energy. These methods rely upon the standard harmonic normal mode analysis to calculate the vibrational and rotational contributions to the energy. We present a method for going beyond the harmonic analysis, which uses path integral Monte Carlo to calculate the vibrational and rotational contributions. Anharmonic effects are found to be as large as 2.5 kcal/mol for the molecules studied. Analytical methods for determining an optimal path discretization are presented. A novel potential energy caching scheme, which greatly improves computational efficiency, is also presented.