One of the most challenging problems in computer-aided drug design is deducing one or more models of a binding site solely from the known chemical structures and measured binding affinities of a few small molecules. In particular, the X-ray crystal structure of the receptor is not known. The main difficulty is deciding how the different ligands are supposed to bind to the common site, especially when they differ substantially in chemical structure. The usual approach is to first guess or use some separate algorithm to superimpose at least the active compounds on one another and then adjust other parameters associated with the model to fit the observed affinities without further altering the proposed alignment. We have explored a long series of algorithms that instead consider the many ways each ligand could fit in the site, and then adjust the site parameters until the best way each molecule can fit has a calculated binding affinity that agrees with experiment. Thus each of the different binding modes of each molecule is some sort of constraint on the geometry and energetic features of the site model. Here we show how to unite the geometric and energetic constraints of the many different binding modes of conformationally flexible ligands in a common framework of mixed integer programming. When applied to the standard test case of corticosteroid binding globulin, even a single compound as the training set can produce a variety of site models having considerable predictive power.