A Monte Carlo docking procedure that combines random displacements of the substrate and protein side chains with minimization of the enzyme—substrate complex is described and applied to finding the binding mode of the blocked tetrapeptide N-acetyl-Leu-Pro-Phe-methylamide to the FK506 binding protein (FKBP). The tetrapeptide, an analog of the preferred FKBP substrate, and the FKBP binding site are flexible during the docking procedure. The twisted-imide transition-state form of the substrate is used during docking. The enzyme charges are scaled individually to account for solvent screening of specific binding site residues during the Monte Carlo sampling. To evaluate the relative binding free energies of the resulting structures, a rapid method for calculating polar and nonpolar solvation effects is introduced. Accurate electrostatic solute—solvent energies are calculated by solving the finite-difference linearized Poisson—Boltzmann equation; nonpolar contributions to the stability of the different conformers are estimated by the free energy of cavity formation, which is obtained from the molecular surface, and the solute—solvent van der Waals energy, which is calculated with a continuum approach. In the conformation of the enzyme—substrate complex with the lowest free energy, the tetrapeptide is bound as a type VIa proline turn with solvent accessible ends to permit longer polypeptide chains to act as substrates. Except for the imide carbonyl, which is involved in polar interactions with aromatic side chains of the FKBP binding site, all of the seven potential hydrogen bond donors or acceptors of the tetrapeptide are satisfied. The FKBP binding site has a similar conformation in the substrate complex as in the FKBP-FK506 cocrystal structure, except for the predicted reorientation of the Tyr 82 hydroxyl, which plays an important role in substrate binding. The present model for the FKBP—substrate complex is in agreement with the recently determined crystal structure of a cyclic peptide—FK506 hybrid bound to FKBP and supports the structure obtained previously by iterative model building. In addition, it is consistent with the observed effects of FKBP point mutations on the enzyme activity. The approach described here should be useful, in general, for the prediction of the structure of a molecule in solution or as part of a complex. It provides for the effective sampling of conformational space and for the inclusion of solvent effects.