In this paper we describe the implementation of a very fast molecular dynamic method to realistically handle electrostatic interactions in simulations of solvated proteins. Our scheme is based on a recently proposed reversible multiple time step (r-RESPA) algorithm and a new modification of the particle mesh Ewald method. While the latter technique provides a fast and accurate representation of the Coulombic interactions for infinite systems, the r-RESPA algorithm exploits a separation of the particle force into components with increasingly longer time scales corresponding to contributions from short-, medium-, and long-range radial shells, By combining the two techniques we are able to reduce considerably the computational cost of molecular dynamics simulation of large biomolecular system without affecting energy conservation and dynamical properties. With respect to single time step simulations employing standard Ewald summation and rigid bond constraints, we obtain a speed-up of about 1 order of magnitude, Finally, our method is about 2.5 times as fast as simulations making use of spherical cutoffs. Since the majority of today biomolecular simulations use spherical cutoffs, we expect that our algorithm will find general applications in the field.