The solubility of molecular crystals is of interest in many areas of chemistry, of which pharmaceutical applications are of particular importance. Predicting solubility using atomistic first-principle methods, that compare the chemical potential of the solid and the solvated phase, has become more common, but remains challenging due to the difficulty in modeling both the crystal form, including its polymorphs, and the interactions with the solvent. Here we aim to compute the solubilities of three active pharmaceutical ingredients of increasing size and complexity: paracetamol, carbamazepine and indomethacin. The known anhydrous polymorphs of each compound are considered in the calculation of the free energy of the solid form and the solvation is explored both in water and in some organic solvents (ethanol, methanol and acetonitrile). The compounds are categorized as poorly soluble or well-soluble based on the comparison of the solid form free energy and the solvation free energy of a single molecule at infinite dilution. Poor solubility then prompts the use of a solubility estimation based on excess free energies, while for well-soluble molecules their chemical potential has to be calculated as a function of concentration. This is done using the recently developed S⁰ method. While promising results are obtained for paracetamol, carbamazepine and indomethacin predictions systematically underestimate the solubility. This can be ascribed to incomplete descriptions of intermolecular interactions by the force field, but the order of stability of the solid forms and systematic nature of the deviations point toward additional problems originating in the accuracy of the method used to calculate the solid free energies.