Electrostatic interactions play a crucial role in many biomolecular processes, including molecular recognition and binding. Biomolecular electrostatics is modulated to a large extent by the water surounding the molecules. Here, we present a novel approach to the computation of electrostatic potentials which allows the inclusion of water structure into the classical theory of continuum electrostatics. Based on our recent purely differential formulation of nonlocal electrostatics (Hildebrandt et al., Phys. Rev. Lett. {\bf 93} (2004) 108104) we have developed a new algorithm for its efficient numerical solution. The key component of this algorithm is a boundary element solver, having the same computational complexity as established boundary element methods for local continuum electrostatics. This allows, for the first time, the computation of electrostatic potentials and interactions of large biomolecular systems immersed in water including effects of the solvent's structure in a continuum description. We illustrate the applicability of our approach with two examples, the enzymes trypsin and acetylcholinesterase (AChE). The approach is applicable to all problems requiring precise prediction of electrostatic interactions in water, such as protein-ligand and protein-protein docking, folding, and chromatin regulation. Initial results indicate that this approach may shed a new light on biomolecular electrostatics and on aspects of molecular recognition that classical local electrostatics cannot reveal.