The indirect metal proton spin spin coupling constants between protons in cysteine groups and the mercury or cadmium nucleus have been calculated for a small model of Me rubredoxin complex (Me = Cd, Hg) by means of density functional theory with zeroth-order regular approximation Hamiltonian (DFT-ZORA). The calculated spin spin coupling constants, in spite of the moderate size of the model system, are in good agreement with the values measured in NMR experiment, which are in the 0.29-0.56 Hz range for the Cd complex and in the 0.57-2.20 Hz range for the Hg complex. The robustness of the chosen method has been verified by calculations with a number of different exchange-correlation functionals and basis sets. Additionally, it has been shown that the short- and long-distance metal proton coupling constants are affected mainly by the values of the metal-proton distance and the H-N-C-C dihedral angle.