The grand canonical simulation algorithm is considered as a general methodology to sample the configuration of water molecules confined within protein environments. First, the probability distribution of the number of water molecules and their configuration in a region of interest for biochemical simulations, such as the active site of a protein, is derived by considering a finite subvolume in open equilibrium with a large system serving as a bulk reservoir. It is shown that the influence of the bulk reservoir can be represented as a many-body potential of mean force acting on the atoms located inside the subvolume. The grand canonical Monte Carlo (GCMC) algorithm, augmented by a number of technical advances to increase the acceptance of insertion attempts, is implemented, and tested for simple systems. In particular, the method is illustrated in the case of a pure water box with periodic boundary conditions. In addition, finite spherical systems of pure water and containing a dialanine peptide, are simulated with GCMC while the influence of the surrounding infinite bulk is incorporated using the generalized solvent boundary potential [W. Im, S. Berneche, and B. Roux, J. Chem. Phys. 114, 2924 (2001)]. As a last illustration of water confined in the interior of a protein, the hydration of the central cavity of the KcsA potassium channel is simulated. (C) 2004 American Institute of Physics.