This is the second of a series of papers discussing the possibility of separating and accurately calculating electrostatic and polarization energies in simulations using classical force fields. A method is described for determining a set of effective distributed multipoles which have significantly improved convergence properties in evaluating the electrostatic interaction energy between molecules. These fitted multipoles are derived to reproduce the electrostatic potential and its derivatives as calculated from a full distributed multipole analysis. The method is based on previous work on the determination of multipole-fitted charges (Ferenczy, G. G. J. Comput. Chem. 1991, 12, 913; Chipot et al. J. Phys. Chem. 1993, 97, 6628) and does not involve the use of a numerical grid. In applications on model systems, fitted charges and dipoles are able to reproduce both the interaction energy and the optimized geometry obtained from a full distributed multipole analysis. Potential-derived charges, however, result in significant errors when the molecules are in close proximity to each other. The method was also used to investigate a possible reason why norepinephrine has a higher affinity than epinephrine in the beta(1)-adrenergic receptor subtype, while the specificity is reversed in the beta(2)- and beta(3)-adrenergic receptor subtypes. This new method offers much potential in the development of new force fields, particularly those involving polarization through induced dipoles, because only fitted charges and dipoles are required to reproduce quantitatively electrostatic interactions.