A self consistent nonlocal density functional method (NL-SCF) based on the gradient corrected exchange term by Becke [Phys. Rev. A38, 3098(1988)] and a gradient corrected correlation term by Perdew [Phys. Rev. B33, 8822(1986)] has been implemented into the LCAO-HFS program by Baerends et al. [Chem. Phys. 2, 41(1973)]. An expression for the variational potential corresponding to Becke's gradient correction has been derived. A detailed analysis is given of the change in density, dr, induced by the gradient corrections. It is shown that dr is two orders of magnitude smaller than the deformation density Dr = rmolecule -ratom. The dipole moments calculated by NL-SCF differ as a consequence only modestly from dipole moments obtained by the Local Density Approximation (LDA) in which the gradient corrections are absent. The importance of self-consistency for calculations on bond distances, bond energies and vibrational frequencies are assessed by comparing results from NL-SCF calculations with data obtained by a perturbative approach (NL-P) in which densities obtained from LDA calculation were used in the evaluation of the gradient corrections. It is further shown that results from calculations on metal-ligand bond distances and bond energies are improved considerably by nonlocal corrections.