  We present a molecular dynamics simulation method for the computation of the
solubility of organic crystals in solution. The solubility is calculated based
on the equilibrium free energy difference between the solvated solute and its
crystallized state at the crystal surface kink site. In order to efficiently
sample the growth and dissolution process, we have carried out well-tempered
Metadynamics simulations with a collective variable that captures the slow
degrees of freedom, namely the solute diffusion to and adsorption at the kink
site together with the desolvation of the kink site. Simulations were performed
at different solution concentrations using constant chemical potential
molecular dynamics and the solubility was identified at the concentration at
which the free energy values between the grown and dissolved kink states were
equal. The effectiveness of this method is demonstrated by its success in
reproducing the experimental trends of solubility of urea and naphthalene in a
variety of solvents.
