On the calculation of solubilities via direct coexistence simulations: Investigation of NaCl aqueous solutions and Lennard-Jones binary mixtures.

Direct coexistence molecular dynamics simulations of NaCl solutions and Lennard-Jones binary mixtures were performed to explore the origin of reported discrepancies between solubilities obtained by direct interfacial simulations and values obtained from the chemical potentials of the crystal and solution phases. We find that the key cause of these discrepancies is the use of crystal slabs of insufficient width to eliminate finite-size effects. We observe that for NaCl crystal slabs thicker than 4 nm (in the direction perpendicular to the interface), the same solubility values are obtained from the direct coexistence and chemical potential routes, namely, 3.7 ± 0.2 molal at T = 298.15 K and p = 1 bar for the JC-SPC/E model. Such finite-size effects are absent in the Lennard-Jones system and are likely caused by surface dipoles present in the salt crystals. We confirmed that μs-long molecular dynamics runs are required to obtain reliable solubility values from direct coexistence calculations, provided that the initial solution conditions are near the equilibrium solubility values; even longer runs are needed for equilibration of significantly different concentrations. We do not observe any effects of the exposed crystal face on the solubility values or equilibration times. For both the NaCl and Lennard-Jones systems, the use of a spherical crystallite embedded in the solution leads to significantly higher apparent solubility values relative to the flat-interface direct coexistence calculations and the chemical potential values. Our results have broad implications for the determination of solubilities of molecular models of ionic systems.