Electron groups provide a natural way to introduce local concepts into quantum chemistry, and the wave functions based on the group products can be considered as a framework for constructing efficient computational methods in terms of "observable" parts of molecular systems. The elements of the group wave functions (electronic structure variables) can be optimized by requiring the number of operations proportional to the size of the molecule. This directly leads to computational methods linearly scaling for large molecular systems. In the present work we consider a particular case of such a wave function implemented for the semiempirical NDDO Hamiltonian. The electron groups are expressed in terms of optimized atomic (hybrid) orbitals with chemical bonds described by geminals and the delocalized groups described by Slater determinants (with or without spin restriction). This scheme is very fast by itself but its speed is considerably limited by the computations of the interatomic Coulomb interactions. Here we develop a consistent method based on group functions which uses the multipole scheme for interatomic interactions. The explicit usage of the atomic multipoles makes the method extremely fast, although the numerical efficiency is largely achieved due to the local character of the electron groups involved. We discuss numerical characteristics of the new method as well as its possible parametrization. We apply this method to study dodecahedral water clusters with hydrogen fluoride substitution and base the analysis on the exhaustive calculation of all symmetry-independent hydrogen-bond networks.