The principles behind the computation of protein - ligand binding free energies by Monte Carlo simulation in the Grand Canonical Ensemble are described in detail, and two variations of the calculation are presented. The computation can be performed by bathing a protein binding site either with ligand images that interact with each other or with ligand images that pass through each other. The second method is theoretically more rigorous, but we show that both methods lead to the same result, and there are distinct numeric advantages to using ligand images that interact with each other. The Grand Canonical simulation provides gas-phase binding free energies that can be converted to aqueous energies by generalized Born-surface area (GB/SA) solvation calculations to provide values that agree with experiment within ±1.5 kcal/mol. However, the accuracy of these simple solvation calculations is a major limiting factor in the accuracy of the overall binding free-energy computation. The Grand Canonical simulation has several characteristics beneficial to free-energy calculations. One is that the number of parameters that must be set for the simulation is small and can be determined objectively, making the outcome more deterministic, with respect to choice of input conditions, as compared to perturbation methods. Second, the simulation is free from assumptions about the starting pose or nature of the binding site. A final benefit is that binding free energies are a direct outcome of the simulation, and little processing is required to determine them.