In this work, we have analyzed the key theoretical and computational considerations for obtaining accurate chemical shifts in paramagnetic systems from QC methods. We have paid particular attention to the methods which use the correct physics in order to obtain the best possible calculated shifts. We have highlighted two important points through our choice of model transition metal complexes from the literature with a significant number of the available ex-perimental NMR and magnetic properties data. The complexes are chosen to emphasize the different contributions to the paramagnetic shifts i.e. the contact and the pseudocontact parts. Firstly, we have analyzed the predictive performances of the DFT functionals across the hierarchy and the core-augment basis sets in characterizing the experimental chemical shifts. For nuclei with significant contact shifts, the results vary considerably according to the method used. Here we reaffirm that currently the contact shifts are best reproduced using hybrid functionals, specially with range-separated exchange. In addition at the atomic level the local bonding and radial correlation effects are also important through the description of the basis sets in the better performance for a given DFT method. Secondly, we have addressed the recent controversy over the "exact" equations that should be used to describe the paramagnetic shifts in the transition metal complexes. The correct formalism is given by eq. 1 with the evaluation the χ s tensor according to eq. 9. For metal ion with significant SOC and low lying states (eg. Co 2+ ) at a ligand environment, the wrong choice of formalism can lead to significant errors in the calculated shifts due to poor description of the magnetic anisotropy and its impact on the dipolar shifts. On the other hand, for metal ions with no (eg. Mn 2+ ) or small (eg. Ni 2+ ) SOC, the discrepancy between eq. 9 and other formalisms is probably less significant. The VV formalism properly includes contributions from the low lying excited states unlike the mPV (eq. 10) and BB (eq. 11) formalisms. One significant limitation to the current state-of-the-art is that we cannot currently include the state-specific hyperfine tensor contributions from each individual KD. The use of the pointdipole approximation circumvents these shortcoming for anisotropic part for the hyperfine tensor, however for the contact interaction and indeed a proper DFT calculated spin-dipolar contribution there is no equivalent remedy.