We combine path-integral Monte Carlo methods with high-quality potential energy surfaces to compute equilibrium isotope effects in a variety of systems relevant to 'clumped' isotope analysis and isotope geochemistry, including CO2, N2O, methane, and propane. Through a systematic study of heavy-atom isotope-exchange reactions, we quantify and analyze errors that arise in the widely used Urey model for predicting equilibrium constants of isotope-exchange reactions using reduced partition function ratios. These results illustrate that the Urey model relies on a nontrivial cancellation of errors that can shift the apparent equilibrium temperature by as much as 35 K for a given distribution of isotopologues. The calculations reported here provide the same level of precision as the best existing analytical instrumentation, resolving the relative enrichment of certain isotopologues to as little as 0.01 parts per thousand. These findings demonstrate path-integral methods to be a rigorous and viable alternative to more approximate methods for heavy-atom geochemical applications.