We present an O(N) multigrid-based method for the efficient calculation of the long-range electrostatic forces needed for biomolecular simulations, that is suitable for implementation on massively parallel architectures. Along general lines, the method consists of: (i) a charge assignment scheme, which both interpolates and smoothly assigns the charges onto a grid; (ii) the solution of Poisson's equation on the grid via multigrid methods; and (iii) the back interpolation of the forces and energy from the grid to the particle space. Careful approaches for the charge assignment and the force interpolation, and a Hermitian approximation of Poisson's equation on the grid allow for the generation of the high-accuracy solutions required for high-quality molecular dynamics simulations. Parallel versions of the method scale linearly with the number of particles for a fixed number of processors, and with the number of processors, for a fixed number of particles.