The interaction between nuclear spins in a molecule is exceptionally sensitive to the physics beyond the standard model. However, all present calculations of the nuclear spin-spin coupling constant J are burdened by computational difficulties, which hinders the comparison to experimental results. Here, we present a variational approach and calculate the constant J in the hydrogen molecule with the controlled numerical precision, using the adiabatic approximation. The apparent discrepancy with experimental result is removed by an analysis of nonadiabatic effects based on the experimental values of the J constant for HD, HT, and DT molecules. This study significantly improves the reliability of the NMR theory for searching new physics in the spin-spin coupling.