An accurate and numerically efficient treatment of electrostatics is essential for biomolecular simulations, in particular, when a smooth interface to quantum chemical descriptions is needed. Force field used in classical biomolecular simulation codes such as AMBER and CHARMM assign "partial charges" to every atom in a simulation in order to model the interatomic electrostatic forces. The respective charge values are obtained via least-squares fitting to the Coulombic potential produced by quantum chemical procedures Unfortunately, the fitting procedure for large, conformationally flexible molecules is under-determined, which is a major source of errors. There are two main problems associated with the treatment of classical electrostatics: (i) how does one eliminate artifacts associated with the point charges as used in force fields, and thereby improve the electrostatic potentials in a physically meaningful way?; (ii) how does one efficiently simulate the very costly long-range electrostatic interactions? Here, we present results on a recently developed distributed multipole method and discuss the importance of this method for large scale biomolecular simulations.