A quasi-restricted orbital (QRO) approach for the calculation of the spin-orbit term of zero-field splitting tensors (D-SO tensors) by means of density functional theory (DFT) importantly features in the fact that it is free from spin contamination problems because it uses spin eigenfunctions for the zeroth order wave functions. In 2011, however, Schmitt and co-workers pointed out that in the originally proposed QRO working equation some possible excitations were not included in their sum-over-states procedure, which causes spurious D-SO contributions from closed shell subsystems located far from the magnetic molecule under study. We have revisited the derivation of the QRO working equation and modified it, making it include all possible types of excitations in the sum-over-states procedure. We have found that the spurious D-SO contribution can be eliminated by taking into account contributions from all possible types of singly excited configuration state functions. We have also found that only the SOMO(alpha) -> SOMO(beta) excited configurations have nonzero contributions to the D-SO tensors as long as alpha and beta spin orbitals have the same spatial distributions and orbital energies. For the D-SO tensor calculations, by using a ground state wave function free from spin contamination, we propose a natural orbital-based Pederson-Khanna (NOB-PK) method, which utilizes the single determinant wave function consisting of natural orbitals in conjunction with the Pederson-Khanna (PK) type perturbation treatment. Some relevant calculations revealed that the NOB-PK method can afford more accurate D-SO tensors than the conventional PK method as well as the QRO approach in Mn-II complexes and Re-IV-based single molecule magnets.