The application of the zeroth-order regular relativistic approximation (ZORA) for molecular density functional calculations is investigated. By introducing a model potential to construct the kinetic energy operator, stationarity of the energy with respect to orbital variations is gained and most problems connected with gauge dependence of the regular approximation are eliminated. The formulation of a geometry gradient is greatly facilitated using this formalism. Calculations for the coinage metal hydrides (Cu-2, Ag-2, Au-2) as well as for the homonuclear (Cu-2, Ag-2, Au-2 and heteronuclear (CuAg, CuAu, AgAu) diatomics show that the results of ZORA calculations within the electrostatic shift approximation, as introduced by van Lenthe and co-workers, can be duplicated using the simpler scheme proposed in this work. Results for the coinage metal fluorides (CuF, AgF, AuF) and chlorides (CuCl, AgC1, AuCl) are presented as well. First-order relativistic calculations have been performed for all systems to assess the applicability of leading-order relativistic perturbation theory.