Incorporation of relativistic effects into all-electron density functional (DF) calculations via scalar relativistic corrections to the kinetic energy and the nuclear attraction energy has become rather common. On the other hand, a self-consistent treatment of the spin-orbit interaction and relativistic corrections to the electron-electron interaction (or both simultaneously) requires substantially more computational effort. We present an extension of the Douglas-Kroll approach to the Dirac-Kohn-Sham problem that accounts for relativistic corrections to the Hartree potential and permits a self-consistent treatment of spin-orbit interaction. To construct computationally efficient approximations, we exploit the electron charge density fitting scheme with an auxiliary basis set. These approximate schemes introduce effects of the relativistic transformation of the Hartree part of the electron-electron interaction, but leave the (smaller) exchange-correlation contributions untransformed. These approximations were implemented in the parallel DF program PARAGAUSS. Quantitative effects of the new relativistic DF procedures were illustrated for the spin-orbit splittings of the Kohn-Sham levels in the Hg atom and the g-tensor shifts of NO2 where we show the value of new scheme. We also studied how properties of the diatomic molecules TlH, PbO, Pb-2, and Bi-2 change due to the improved treatment of relativistic effects. (C) 2003 American Institute of Physics.