Modeling molecular flexibility continues to be a frontier in predicting how two proteins bind. Allowing full dynamic flexibility can be excessive, yet neither traditional rigid body docking nor small-scale rearrangements of amino acid side chains can capture many biologically relevant motions. We propose an intermediate approach, which is to combine rigid optimization with discrete deformations that represent a protein's ``natural'' motions. Through exhaustive low-resolution search, hierarchical clustering, and global optimization within each energy basin, we are able to achieve accurate docking predictions. Our optimization approach employs a Convex Global Underestimation method that has been adapted for energy functions defined on $R^3 \times SO(3)$.