A method for the generation of highly accurate, nearly-exact, full-dimensional interaction energy surfaces for weakly interacting subsystems is proposed. The method is based on the local expansion of the exact interaction energy surface in the Taylor series with respect to intramolecular coordinates. It is shown that without any significant loss of accuracy this expansion can be limited to a few low-order terms. This leads to significant savings in computations of the full-dimensional interaction energy surfaces. Also a method for the direct calculation of the interaction energy surface of reduced dimensionality, corresponding to averaging over the intramolecular vibrations, without explicit knowledge of the full-dimensional surface, is presented. The main ideas and computational features of the proposed scheme are comprehensively tested for the Ar-HF system. (C) 2004 American Institute of Physics.