We present a gauge-origin independent method for the calculation of nuclear magnetic shielding tensors of molecules in a structured and polarizable environment. The method is based on a combination of density functional theory (DFT) or Hartree-Fock wave functions with molecular mechanics. The method is unique in the sense that it includes three important properties that need to be fulfilled in accurate calculations of nuclear magnetic shielding constants: (i) the model includes electron correlation effects, (ii) the model uses gauge-including atomic orbitals to give gauge-origin independent results, and (iii) the effect of the environment is treated self-consistently using a discrete reaction-field methodology. The authors present sample calculations of the isotropic nuclear magnetic shielding constants of liquid water based on a large number of solute-solvent configurations derived from molecular dynamics simulations employing potentials which treat solvent polarization either explicitly or implicitly. For both the (17)O and (1)H isotropic shielding constants the best predicted results compare fairly well with the experimental data, i.e., they reproduce the experimental solvent shifts to within 4 ppm for the (17)O shielding and 1 ppm for the (1)H shielding.