Accurately predicting the relative free energies of polymorphic molecular crystals is an important aspect of crystal structure prediction (CSP), which has found considerable utility in solid-form pharmaceutical development to assess the risk of conversion to a more stable, late-appearing polymorph. In this work, we investigated the Einstein-Debye phonon approximation to evaluate the vibrational free energy, \(F_{\text{vib}}\), in conjunction with dispersion-corrected density-functional theory. Three data sets were considered: (1) The "PV17'" benchmark of seventeen polymorph pairs exhibiting little or no conformational flexibility; (2) four large, flexible compounds that appeared in previous CSP blind tests; and (3) a new "FP10" set of 10 highly flexible drug molecules that each have two or more known polymorphs. It was found that the Einstein-Debye approximation provides a good balance of accuracy and efficiency, giving mean absolute errors of \(\le 1.4\)~kJ/mol relative to full supercell calculations of \(F_\text{vib}\), with a computational cost that is up to 4.5 times lower. Considering the magnitudes of free-energy differences between polymorphs, \(|\Delta F_{\text{vib}}|\), a very broad distribution was observed, with average values in excess of 3 kJ/mol. In CSP studies of drug-like molecules, we recommend that \(\Delta F_{\text{vib}}\) be considered for all candidate structures within at least ca.~6~kJ/mol of the global electronic-energy minimum to provide a high probability of identifying the correct free-energy minimum.