Long range electrostatic forces are involved at a fundamental level in many biological phenomena. Their prohibitive computational costs often prevents their correct calculation in molecular dynamics (MD) simulations of biological molecules. In this paper we present a method to handle efficiently and exactly electrostatic interactions in MD simulations with periodic boundary conditions. Our scheme employs a multiple time step r-RESPA integration algorithm in combination with the Ewald summation technique, and is specifically targeted to simulation of large size complex molecular systems such as solvated proteins. In this approach, the force associated with each particle of the system is partitioned into four components which evolve in time with distinct and increasingly longer time scales. We found that a suitable time scale separation is achieved by subdividing the direct space nonbonded interactions, inclusive of Coulombic and van der Waals contributions, in a short, medium and long range shells. The fastest bonded forces are associated to the component with the shortest timescale, while the reciprocal space nonbonded interaction is included with the medium range direct space forces. Our method is general. It can be straightforwardly implemented for any biomolecular force fields used in condensed phase simulations and can be applied to other complex molecular systems such as molecular liquids and heterogeneous solutions (e.g., micelles, membranes, etc.). We carried out tests on solvated proteins samples containing 7 040 and 20 627 atoms. Due to the nonlinear scaling of the Ewald computational cost with system size, the performance of our r-RESPA algorithm with respect to single time step algorithms with bond constraints grows with the number of particles. We reach a relative speed-up of 4.2 for a system composed of a solvated photosynthetic reaction center.