The scalar relativistic variant of the linear combination of Gaussian-type orbitals-fitting functions-density-functional (R-LCGTO-FF-DF) method is extended to a two-component scheme which permits a self-consistent treatment of the spin-orbit interaction. The method is based on the Douglas-Kroll transformation of the four-component Dirac-Kohn-Sham equation. The present implementation in the program PARAGAUSS neglects spin-orbit effects in the electron-electron interaction. This approximation is shown to be satisfactory as long as bonding is restricted to s, p, and d orbitals. The method is applied to the diatomics Au-2, Bi-2, Pb-2, PbO, and TlH using both a local density (LDA) and a gradient-corrected approximation (GGA) of the exchange-correlation functional. At the LDA level, bond lengths and vibrational frequencies are reproduced with high accuracy. For the determination of binding energies the open-shell reference atoms Au, Tl, Pb, Bi have been treated by a jj coupling approach based on a self-consistent noncollinear spin density-functional scheme and with an intermediate coupling procedure. The atomic state energies obtained with the jj coupling scheme agree well with experiment, but they are somewhat too high due to the incomplete inclusion of static correlation. Binding energies of diatomics at the GGA level are considerably improved due to the inclusion of spin-orbit interaction. The jj derived values are somewhat overestimated (by about 10%) compared to experiment, and they compare slightly worse with experiment than results based on the intermediate coupling approximation.