A streamlined molecular-dynamics workflow for computing solubilities of molecular and ionic crystals

TL;DR

This paper presents a simplified and robust molecular-dynamics workflow for accurately computing the solubility of molecular and ionic crystals, addressing challenges of convergence and slow equilibration in traditional methods.

Contribution

The authors introduce a streamlined workflow that optimizes each step for calculating crystal solubility, improving robustness and ease over existing state-of-the-art methods.

Findings

Predicted solubilities are sensitive to potential energy surface choices.

Harmonic approximation often fails at or above room temperature.

Ideal solution assumption is less valid for highly soluble substances.

Abstract

Computing the solubility of crystals in a solvent using atomistic simulations is notoriously challenging due to the complexities and convergence issues associated with free-energy methods, as well as the slow equilibration in direct-coexistence simulations. This paper introduces a molecular-dynamics workflow that simplifies and robustly computes the solubility of molecular or ionic crystals. This method is considerably more straightforward than the state-of-the-art, as we have streamlined and optimised each step of the process. Specifically, we calculate the chemical potential of the crystal using the gas-phase molecule as a reference state, and employ the S0 method to determine the concentration dependence of the chemical potential of the solute. We use this workflow to predict the solubilities of sodium chloride in water, urea polymorphs in water, and paracetamol polymorphs in both water and ethanol. Our findings indicate that the predicted solubility is sensitive to the chosen potential energy surface. Furthermore, we note that the harmonic approximation often fails for both molecular crystals and gas molecules at or above room temperature, and that the assumption of an ideal solution becomes less valid for highly soluble substances.