A general quantum-mechanical method for computing kinetic isotope effects is presented. The method is based on the quantum-instanton approximation for the rate constant and on the path-integral Metropolis–Monte Carlo evaluation of the Boltzmann operator matrix elements. It computes the kinetic isotope effect directly, using a thermodynamic integration with respect to the mass of the isotope, thus avoiding the more computationally expensive process of computing the individual rate constants. The method should be more accurate than variational transition-state theories or the semiclassical instanton method since it does not assume a single tunneling path and does not use a semiclassical approximation of the Boltzmann operator. While the general Monte Carlo implementation makes the method accessible to systems with a large number of atoms, we present numerical results for the Eckart barrier and for the collinear and full three-dimensional isotope variants of the hydrogen exchange reaction H+H2H2+H. In all seven test cases, for temperatures between 250 and 600 K, the error of the quantum instanton approximation for the kinetic isotope effects is less than ~10%.