The self-consistent scalar-relativistic linear combination of Gaussian-type orbitals density functional (LCGTO-DF) method has been extended to calculate analytical energy gradients. The method is based on the use of a unitary second order Douglas-Kroll-Hess (DKH) transformation for decoupling large and small components of the full four-component Dirac-Kohn-Sham equation. The approximate DKH transformation most common in molecular calculations has been implemented; this variant employs nuclear potential based projectors and it leaves the electron-electron interaction untransformed. Examples are provided for the geometry optimization of a series of heavy metal systems which feature a variety of metal-ligand bonds, like Au-2, AuCl, AuH, Mo(CO)(6) and W(CO)(6) as well as the d(10) complexes [Pd(PH3)(2)O-2] and [Pt(PH3)(2)O-2]. The calculated results, obtained with several gradient-corrected exchange-correlation potentials: compare very well with experimental data and they are of similar or even better accuracy than those of other high quality relativistic calculations reported so far.