Prediction of aqueous solubility of a strongly soluble solute from molecular simulation

The prediction of solubilities of compounds by means of molecular simulation has been receiving increasing attention due to the key role played by solubility in countless applications. We have predicted the aqueous solubility of urea at 300 K from chemical potential calculations for two urea model combinations: ozpinar/TIP3P and Holzl/(TIP4P/2005). The methodology assumes that the intramolecular contribution of the urea molecule to the chemical potentials is identical in the crystal and in solution and, hence, cancels out. In parallel to the chemical potential calculations, we also performed direct coexistence simulations of a urea crystal slab in contact with urea-water solutions with the aim to identify upper and lower bounds to the solubility value using an independent route. The chemical potential approach yielded similar solubilities for both urea models, despite the actual chemical potential values showing a significant dependence on the force field. The predicted solubilities for the two models were 0.013-0.018 (ozpinar) and 0.008-0.012 (Holzl) mole fraction, which are an order of magnitude lower than the experimental solubility that lies in a range of 0.125-0.216 mole fraction. The direct coexistence solubility bounds were relatively wide and did not encompass the chemical potential based solubilities, although the latter were close to the lower bound values.