A methodology is suggested for modeling the phase equilibria of complex chemical mixtures with an equation of state (EoS) for the case where only limited experimental data exist. The complex hydrogen bonding behavior is explicitly accounted for and the corresponding parameters are adopted from simpler molecules of similar chemical structure and/or are fitted to Hansen's partial solubility parameters. The methodology is applied to modeling the solubility of three pharmaceuticals, namely acetanilide, phenacetin, and paracetamol, using the nonrandom hydrogen bonding (NRHB) EoS. In all cases, accurate correlations were obtained. The prediction ability of the approach was evaluated against predictions from the COSMO-RS model. A thorough discussion is made for the appropriate modeling of solid solubility considering the effect of the difference of the heat capacities of the solute in liquid and solid state, Delta Cp = Cp-l - Cps, in the determination of solid chemical potential and, also, of the polymorphism of drugs. (c) 2009 American Institute of Chemical Engineers AIChE J, 55: 756-770, 2009