A method for predicting the effect of solvent on the morphology of organic crystals is presented, providing an efficient screening tool for identifying ideal crystallization solvents. The solvent effect is estimated by the computation of chemical potentials and activity coefficients of crystal surfaces using a first principles-based statistical thermodynamics approach. Density functional theory and COSMO-RS are utilized to determine the activity coefficients of the crystal growth faces of a selection of active pharmaceutical ingredients (APIs) in solvents across a broad range of polarities. The ability of COSMO-RS to predict and quantify the effects of solvent on crystal growth and morphology is assessed using hierarchical clustering to classify the solvents according to their overall interaction strength with the crystal faces. The COSMO-RS approach allows for a physical interpretation of the predictions in terms of surface polarity and is confirmed by comparison to published experimental data. Herein a methodology is reported for automated computation of the activity coefficients of all solvent-surface pairs directly from the drug crystal structure. The procedure goes beyond the traditional trial-and-error solvent selection process and has the potential to be used as a rapid computational screening tool in pharmaceutical drug development.