New implementations of the Levy-Leblond, zeroth-order regular approach (ZORA) and spin-free Dirac equation are presented within the framework of the four-component relativistic program system DIRAC. This implementation allows systematic incorporation of relativistic effects at different levels of theory and corresponding computational cost. One of the possibilities of the new code is to neglect the effect of spin-orbit coupling in the orbital optimization process and introduce it in a later stage of the calculation. This method is shown to be unstable despite the boundedness of the spin-orbit operator itself. (C) 2000 American Institute of Physics. [S0021-9606(00)30434-2].