We developed a methodology to optimize exponential damping functions to account for charge penetration effects when computing molecular electrostatic properties using the multicentered multipolar expansion method (MME). This methodology is based in the optimization of a damping parameter set using a two-step fast local fitting procedure and the ab initio (Hartree-Fock/6-31G** and 6-31G**+) electrostatic potential calculated in a set of concentric grid of points as reference. The principal aspect of the methodology is a first local fitting step which generates a focused initial guess to improve the performance of a simplex method avoiding the use of multiple runs and the choice of initial guesses. Three different strategies for the determination of optimized damping parameters were tested in the following studies: (1) investigation of the error in the calculation of the electrostatic interaction energy for five hydrogen-bonded dimers at standard and nonstandard hydrogen-bonded geometries and at nonequilibrium geometries; (2) calculation of the electrostatic molecular properties (potential and electric field) for eight small molecular systems (methanol, ammonia, water, formamide, dichloromethane, acetone, dimethyl sulfoxide, and acetonitrile) and for the 20 amino acids. Our results show that the methodology performs well not only for small molecules but also for relatively larger molecular systems. The analysis of the distinct parameter sets associated with different optimization strategies show that (i) a specific parameter set is more suitable and more general for electrostatic interaction energy calculations. with an average absolute error of 0.46 kcal/mol at hydrogen-bond geometries; (ii) a second parameter set is more suitable for electrostatic potential and electric field calculations at and outside the van der Waals (vdW) envelope, with an average error decrease > 72% at the vdW surface. A more general amino acid damping parameter set was constructed from the original damping parameters derived for the small fragments and for the amino acids. This damping set is more insensitive to protein backbone and residue side-chain conformational changes and can be very useful for future docking and molecular dynamics protein simulations using ab initio based polarizable classical methods.