Automatic search of cavities and binding mode analysis between a ligand and a 3D protein receptor are challenging problems in drug design or repositioning. We propose a solution based on a shape theory theorem for an invariant coupled system of ligand–protein. The theorem provides a matrix representation with the exact formulas to be implemented in an algorithm. The method involves the following results: (1) exact formulae for the shape coordinates of a located-rotated invariant coupled system; (2) a parameterized search based on a suitable domain of van der Waals radii; (3) a scoring function for the discrimination of sites by measuring the distance between two invariant coupled systems including the atomic mass; (4) a matrix representation of the Lennard-Jones potential type 6–12 and 6–10 as the punctuation function of the algorithm for a molecular docking; and (5) the optimal molecular docking as a solution of an optimization problem based on the exploration of an exhaustive set of rotations. We apply the method in the xanthine oxidase protein with the following ligands: hypoxanthine, febuxostat, and chlorogenic acid. The results show automatic cavity detection and molecular docking not assisted by experts with meaningful amino acid interactions. The method finds better affinities than the expert software for known published cavities.