Abstract
We present a method for calculating the relativistic correction in hydrogen molecules that significantly exceeds the accuracy of all the previous literature results. This method utilizes the explicitly correlated nonadiabatic exponential wave function, and thus treats electrons and nuclei equivalently. The proposed method can be applied to any rovibrational state, including highly excited ones. The numerical precision of the relativistic correction reaches several kHz (∼10(-7) cm(-1)), which is below the best experimental accuracy.