Electron-coupled interaction between two nuclear spins IA and IB is of the form hδABIA.IB. In this paper, formulas for calculating δAB are derived by using the variational method. By using these formulas, we have calculated δAB of the HD molecule at internuclear distances of 1.3, 1.4, and 1.5 atomic units. The total wave function was chosen as a linear combination of seven independent functions of 1Σg+ symmetry, representing the unperturbed state, and several additional functions, to represent the perturbation. The average value of δAB over the zero-point vibration, <δAB>00, is found to be 37.138 cps. Further, the refinement of a part of the calculations is undertaken by use of the 11-term James and Coolidge wave function. The final result is <δAB>00=35.217 cps. The agreement with the observed value of 42.7+/-0.7 cps is satisfactory, considering that we have used only a few terms for the additional perturbative wave function and that the result is a sum of terms which cancel each other appreciably. Contributions from each perturbing Hamiltonian and also from a set of wave functions with different symmetry character were obtained separately. The various results are presented and discussed.