Mössbauer spectroscopy is an indispensable technique and analytical tool in iron coordination chemistry. The linear correlation between the electron density at the nucleus ("contact density") and experimental isomer shifts has been used to link calculated contact densities to experimental isomer shifts. Here we have investigated for the first time relativistic methods of systematically increasing sophistication, including the eXact 2-Component (X2C) Hamiltonian and a finite-nucleus model, for the calculation of isomer shifts for iron compounds.While being of similar accuracy as the full four-component treatment, X2C calculations are far more efficient. We find that effects from spin orbit coupling can safely be neglected, leading to further speed up. Linear correlation plots using effective densities rather than contact densities versus experimental isomer shift leads to a correlation constant a = −0.32 a −3 0 mm s −1 (PBE functional) which is in close agreement to experimental findings. Isomer shifts of similar quality can thus be obtained both with and without fitting, which is not the case if one pursues a priori a non-relativistic model approach. As an application for a biologically relevant system, we have studied three recently proposed [Fe]-hydrogenase intermediates. The structures for these intermediates were extracted from QM/MM calculations using large QM regions surrounded by the full enzyme and a solvation shell of water molecules. We show that a comparison between calculated and experimentally observed isomer shifts can be used to discriminate between the different intermediates, whereas calculated atomic charges do not necessarily correlate with Mössbauer isomer shifts. Detailed analysis shows that the difference in isomer shifts between two intermediates is due to an overlap effect.