Here we suggest a mixed computational and empirical approach serving to optimize the parameters of complex and polarizable molecular mechanics (PMM) models for complicated liquids. The computational part of the parameter optimization relies on hybrid calculations combining density functional theory (DFT) for a solute molecule with a PMM treatment of its solvent environment at well-defined thermodynamic conditions. As an application we have developed PMM models for water featuring ν = 3, 4, and 5 points of force action, a Gaussian inducible dipole and a Buckingham potential at the oxygen, the experimental liquid phase geometry, the experimental gas phase polarizability αexpg = 1.47 Å 3, and, for ν = 4 and 5, the gas phase value μexpg = 1.855 D for the static dipole moment. The widths of the Gaussian dipoles and, for ν = 4 and 5, also the electrostatic geometries of these so-called TLνP models are derived from self-consistent DFT/PMM calculations, and the parameters of the Buckingham potentials (and the static TL3P dipole moment) are estimated from molecular dynamics (MD) simulations. The high quality of the resulting models is demonstrated for the observables targeted during optimization (potential energy per molecule, pressure, radial distribution functions) and a series of predicted properties (quadrupole moments, density at constant pressure, dielectric constant, diffusivity, viscosity, compressibility, heat capacity) at certain standard conditions. Remaining deficiencies and possible ways for their removal are discussed.