The implementation of analytical geometry gradients within the framework of the relativistic density functional procedure described earlier allows the calculation of the geometrical structure and bond dissociation energies of polyatomic molecules. This has been done for the nine transition metal carbonyls M(CO)(n) (n=6 : M=Cr, Mo, W; n=5 : M=Fe, Ru, Os; n=4 : M=Ni, Pd, Pt). To determine the first metal-carbonyl bond dissociation energy, a complete geometry optimization of the fragments M(CO)(n-1) has been performed, and the energy differences have been corrected for the basis set superposition error (BSSE). The same procedure has been applied to the molecular structure of the nine complexes M(CO)(5)L (M=Cr, Mo, W; L=N-2, CS, NO+) and their M-L bond dissociation energies. The results are in good agreement with quasirelativistic density functional and high-level ab initio calculations. In most cases, the agreement with experimental values, where available, is good as well.