Adsorption complexes of palladium atoms on F-s, F-s(+), F-s(2+), and O2- centers of MgO(001) surface have been investigated with a gradient-corrected (Becke-Perdew) density functional method applied to embedded cluster models. This study presents the first application of a self-consistent hybrid quantum mechanical/molecular mechanical embedding approach where the defect-induced distortions are treated variationally and the environment is allowed to react on perturbations of a reference configuration describing the regular surface. The cluster models are embedded in an elastic polarizable environment which is described at the atomistic level using a shell model treatment of ionic polarizabilities. The frontier region that separates the quantum mechanical cluster and the classical environment is represented by pseudopotential centers without basis functions. Accounting in this way for the relaxation of the electronic structure of the adsorption complex results in energy corrections of 1.9 and 5.3 eV for electron affinities of the charged defects F-s(+) and F-s(2+), respectively, as compared to models with a bulk-terminated geometry. The relaxation increases the stability of the adsorption complex Pd/F-s by 0.4 eV and decreases the stability of the complex Pd/F-s(2+) by 1.0 eV, but it only weakly affects the binding energy of Pd/F-s(+). The calculations provide no indication that the metal species is oxidized, not even for the most electron deficient complex Pd/F-s(2+). The binding energy of the complex Pd/O2- is calculated at -1.4 eV, that of the complex Pd/F-s(2+) at -1.3 eV. The complexes Pd/F-s and Pd/F-s(+) exhibit notably higher binding energies, -2.5 and -4.0 eV, respectively; in these complexes, a covalent polar adsorption bond is formed, accompanied by donation of electronic density to the Pd 5s orbital.