Context <p>The ability to predict the relative binding energies of ligands to a biological receptor would be of great value in drug discovery. However, accurately calculating the predicted binding energies is limited by the high accuracy required, by the presence of multiple minima on the potential energy surface, and by issues specific to the intrinsic properties of the binding site, such as details of the geometry of the ligand–protein complex. To address these issues, a systematic analysis of potential sources of error was carried out which resulted in a few relatively small changes being made to the MOPAC program.</p> Methods <p>A set of 77 ligands was constructed for which experimentally determined IC₅₀ values were available. For each of the ligands, prediction of the protein–ligand interaction energy was carried out in two distinct stages. In the first stage, the Protein–Ligand docking program GOLD was used to generate several distinct conformations of the ligand bound to a protein. The geometries of these systems were then optimised using the MOPAC program. A comparison of the relative binding energies of the ligands with the reported IC₅₀ values showed a very poor predictive power. By partitioning the ligand set into two subsets, and eliminating six ligands that were inconsistent with the experimental results, a large increase in accuracy was obtained.</p>