A new method is proposed for the treatment of systems consisting of many interacting molecules. Properties of the individual molecules are computed and their response to the influence of external electric fields generated by the partner systems is determined. The necessary quantities are the point charges representing the electronic distribution and the pseudo‐polarization tensors for these point charges which are easily obtained with the help of optimization and iterative procedures, respectively. Thus information is then used in a mutually consistent iteration process to calculate the interaction energy. In this way, not only the pairwise additive energy terms but, in addition, the nonadditive contributions are obtained. The ρ1/3 approximation is used for the exchange energy and the dispersion energy is taken into account employing a London‐type expression, together with empirical parameters. The results of applying this approach to the sensitive systems of two and three interacting water molecules are compared with ab initio supermolecule results. Equilibrium geometries, total binding energies, and nonadditive contributions are in excellent agreement.