We report a multistage lattice energy minimization methodology for generating stable packing arrangements of cocrystals containing flexible molecules. In the first approximation, the intermolecular electrostatic interactions are modeled with atomic charges and the molecular deformation energy is interpolated over a set of precomputed quantum mechanical values. At subsequent stages, the accuracy is improved by first using analytically rotated and then conformation-dependent multipole moments, computed from the isolated-molecule charge density, and "on-the-fly" quantum mechanical calculations to compute the intramolecular deformation energy. This multistage approach increases the efficiency of the search and establishes the molecule-dependent error due to the atomic charge representation of the charge density and the neglect of the conformational dependence of atomic multipole moments. The methodology is used to study the lattice energy landscapes of the cocrystals of 4-aminobenzoic acid with 2,2'-bipyridine and 4-nitrophenylacetic acid, as well as the single-component crystals. All single-component, experimentally determined crystal structures within the scope of the search were found at, or very close to, the global minimum. The experimental cocrystal with 2,2'-bipyridine is also predicted to be among the most stable packing arrangements. On the contrary, the lattice energy landscape of the cocrystal with 4-nitrophenylacetic acid contains several low energy structures that are more stable than the experimentally observed form and have different hydrogen bonding motifs. Overall, the methodology can provide worthwhile crystal energy landscapes for multicomponent organic solids and thereby contribute to understanding cocrystal formation.