In this work we present a new proposal to model intermolecular interactions and use it for water molecules. The parameters of the model were fitted to reproduce the single molecule's electrostatic properties, a sample of 352 points in a refined ab initio single molecule deformation potential energy surface (PES), and the theoretical limit of the dimerization energy, -20.8 kJ/mol. The model was able to reproduce a sample of 180 additional points in the single molecule deformation PES, and 736 points in a pair-interaction surface computed at the MP2/aug-cc-pVQZ(') level with the counterpoise correction. Though the model reproduced the diagonal of the polarizability tensor, it could account for only 60% of the three-body nonadditive contributions to the interaction energies in 174 trimers computed at the MP2/6-311++(2d,2p) level with full counterpoise correction, but reproduced the four-body nonadditivities in 34 tetramers computed at the same level as the trimers. The model's predictions of the structures, energies, and dipoles of small clusters resulted in good agreement with experimental data and high quality ab initio calculations. The model also reproduced the second virial coefficient of steam at various temperatures, and the structure and thermodynamical properties of liquid water. We found that the short-range water-water interactions had a critical influence on the proper performance of the model. We also found that a model based on the proper intermolecular interactions requires the inclusion of intramolecular flexibility to be adequate.