A new method for the rapid calculation of complexation free energy surfaces (FES) in solution is presented. The Monte Carlo (MC) simulation technique for solute-solute statistics is combined with the reference interaction site model (RISM) integral equation theory for solute-solvent and pure solvent statistics. Statistical-mechanical perturbation theory is used to determine complexation potentials of mean force and, from these, standard free energies of binding. The hybrid RISM/MC method is applied to calculate the selectivity of rigid 18-crown-6 for alkali metal ions in liquid methanol and acetonitrile. A full atom model of the host in its D-3h conformation is used for the calculations. The difference of binding free energies between the various ions is found to be in excellent agreement with the experimental data in the case of methanol, less accurate for acetonitrile. Additional information about the structure of the FES is extracted by scanning its topography. Two-dimensional slices through the FES are presented to illustrate the solvation structure of the complexes. By connecting local minima and saddle points, minimum energy pathways are revealed, leading to the globally optimal complex structure.