The theory of a hybrid quantum mechanics/molecular mechanics (QM/MM) approach for gauge-origin independent calculations of the molecular magnetizability using Hartree-Fock or Density Functional Theory is presented. The method is applied to liquid water using configurations generated from classical Molecular Dynamics simulation to calculate the statistical averaged magnetizability. Based on a comparison with experimental data, treating only one water molecule quantum mechanically appears to be insufficient, while a quantum mechanical treatment of also the first solvation shell leads to good agreement between theory and experiment. This indicates that the gas-to-liquid phase shift for the molecular magnetizability is to a large extent of non-electrostatic nature. (c) 2007 Elsevier B.V. All rights reserved.