We develop here a new method to fit the molecular electrostatic potentials obtained in quantum mechanical calculations to a set of classical electrostatic multipoles, usually point charges located at atomic positions. We define an object function of fitting as an integration of the difference of electrostatic potentials in the entire 3-dimensional physical space. The object function is thus rotationally invariant with respect to the molecular orientation and varies smoothly with respect to molecular geometric fluctuations. Compared with commonly employed methods such as the Merz-Singh-Kollman and CHELPG schemes, this new method, while possessing comparable accuracy, shows greatly improved numerical stability with respect to the molecular positions and geometries. The method can be used in the fitting of electrostatic potentials for the molecular mechanics force fields and also can be applied to the calculation of electrostatic polarizabilites of molecular or atomic systems.