We study the problem of representing the molecular charge distribution in a convenient way for practical applications and we propose, instead of a single representation, a flexible procedure for building approximations with an arbitrary level of accuracy as concerns the long-range part of the electrostatic potential. We first discuss the splitting of the total electrostatic potential into a multipolar part (long-range) and a penetration part (shortrange) in connection with the usual one-center multipole expansion: at large enough distances, this expansion precisely converges towards the multipolar part. However, this representation is not practically efficient as soon as the molecule departs from a spherical shape, and we therefore consider a so-called multicenter multipole representation. In the MO-LCAO framework, the use of a basis of Gaussian atomic orbitals {χα} generates such a representation in a natural way: indeed, each elementary distribution χ*αχβ then reduces to a one-center distribution with a finite number of multipoles, hence an exact multipolar representation (corresponding to the approximate MO-LCGO wave function) involving a finite (albeit large, in general) number of centers, each one bearing a finite (small) number of multipoles. From this representation used as a reference, we then generate simplified multicenter multipole representations through a systematic procedure of reduction of the number of centers (having chosen a reduced set of new centers, we represent the multipoles of each suppressed old center through suitable new multipoles located on a few new centers closest to the old one under consideration). The number and locations of the new centers are arbitrary, hence the great flexibility of the method. We checked the accuracy of this procedure for molecules of different polarity (water, benzene, adenine), and we noticeably found that a representation involved as centers, the atoms plus one point per chemical bond exhibited quite a satisfactory accuracy. In the conclusion, we briefly discuss the possibility of extending the method to other kinds of basis functions (e.g., Slater orbitals), in connection with the problem of representing the charge distribution associated with the exact molecular wave function, we also discuss the choice of the number centers according to the desired accuracy and we briefly describe a way of applying the method to the problem of representing a large molecule in terms of molecular fragments.