Molecular-dynamics simulations of urea nucleation from aqueous solution

Abstract

Despite its ubiquitous character and relevance in many branches of science and engineering, nucleation from solution remains elusive. In this framework, molecular simulations represent a powerful tool to provide insight into nucleation at the molecular scale. In this work, we combine theory and molecular simulations to describe urea nucleation from aqueous solution. Taking advantage of well-tempered metadynamics, we compute the free-energy change associated to the phase transition. We find that such a free-energy profile is characterized by significant finite-size effects that can, however, be accounted for. The description of the nucleation process emerging from our analysis differs from classical nucleation theory. Nucleation of crystal-like clusters is in fact preceded by large concentration fluctuations, indicating a predominant two-step process, whereby embryonic crystal nuclei emerge from dense, disordered urea clusters. Furthermore, in the early stages of nucleation, two different polymorphs are seen to compete.

Keywords: enhanced sampling; molecular dynamics; nucleation; solution; well-tempered metadynamics.

Fig. 1.

Total free energy $\left(\mathrm{\Delta }{G}_{\ell \to c}\right)$, bulk contribution $\left(\mathrm{\Delta }{G}_{bulk}=\mathrm{\Delta }{G}_{\ell \to c}-\sigma A\right)$, surface contribution $\left(\sigma A\right)$ obtained from Eq. 7 for the formation of clusters of up to 200 solute molecules for a system consisting of 600 total solute molecules in ${10}^4$ solvent molecules, with a thermodynamic solubility $x^{\mathrm{\ast }}=0.0425$ at $T=273$ K. The surface term $\sigma A$ has been modeled as ${\sigma }^{\prime }{N}^{2/3}$ with ${\sigma }^{\prime }=2.5$. The free-energy profile is compared with $\mathrm{\Delta }{G}_{CNT}$ at constant composition $x=x_0$.

Fig. 2.

$\mathrm{\Delta }{G}_{\ell \to c}$ (reported as $\mathrm{\Delta }G$ in figure) as a function of $N_c/N_{tot}$ and the initial supersaturation $x_0/x^{\mathrm{\ast }}$ (A) and $\log N_s$ (B). In A, $\log N_s$ has been fixed to 3.5, whereas in B, $x_0/x^{\mathrm{\ast }}=2.0$. The red contour represents the locus of the points where $d\mathrm{\Delta }{G}_{\ell \to c}/d{N}_c=0$, and the green contour represents the locus of the points where $\mathrm{\Delta }G=0$. The black dotted line represents a typical free-energy profile obtained in the finite-size case. It exhibits two stationary states, identified by the intersections of the dotted line with the red contour. The first intersection represents the local maximum corresponding to the critical nucleus condition, whereas the second intersection corresponds to the local minimum representing a metastable state in which a finite-sized crystal and the solution coexist. The red dashed line represents the threshold value of $x_0/x^{\mathrm{\ast }}$ (A) or $\log N_s$ (B) that discriminates between a region in which the free-energy profile is monotonically increasing (A) or has a local maximum and a local minimum (B and C). The green dashed line discriminates between a region of the parameter space in which the minimum corresponds to a metastable (B) or to the most stable state (C) for the system.

Fig. 3.

Existence domains of regimes A, B, and C described in Fig. 2 as a function of the initial supersaturation $x_0/x^{\mathrm{\ast }}$ and the size of the system $\log N_s$. In gray is reported the domain in which both macroscopic and confined systems do not have stationary points, i.e., when $x_0/x^{\mathrm{\ast }}<1$, and nucleation is not thermodynamically favored even in a macroscopic system. Exemplary free-energy profiles obtained in the region A, B, and C are displayed as Insets within the plot.

Fig. 4.

Comparison between the fitted analytical expression reported in Eq. 7, using the parameters reported in Table 1 and the free-energy profile obtained from WT metadynamics for the nucleation of urea in water in the limit of $N_{cluster}=1$. In these conditions, the hypotheses used to derive the free-energy profile are satisfied and the analytical expression describes consistently the free-energy profile obtained from simulations.

Fig. 5.

Size of the critical cluster $N^{\mathrm{\ast }}$ as a function of supersaturation calculated at 273 and 300 K using the parameter fitted from simulations S1, S2, and S3. The size of the critical nucleus for an infinite system at a composition equal to that of the homogeneous solution in simulations S2 and S3 is reported as a circle. An analogous point for simulation S1 cannot be identified, as the initial conditions for simulations S1 are such that $x_0/x^{\mathrm{\ast }}<1$. Nevertheless, Eq. 7 also applies in these conditions, providing estimates in good agreement with those of simulation S2.

Fig. 6.

Simulation S2. (A) Contour plot of the FES obtained in the plane $\left(n_o,n\right)$ (see text for a definition) and corrected with the term reported in Eq. 9 for a constant supersaturation $x_0/x^{\mathrm{\ast }}=2.5$. This value of supersaturation has been selected to clearly display the basin belonging to the crystalline nucleus within the $n_o$ interval sampled in simulation. (B) Representative states sampled during the nucleation process. For the sake of simplicity, urea molecules are reported simply as blue spheres. In order to facilitate the visualization of ordered arrangements within dense clusters, red connections are drawn between urea molecules falling within a cutoff distance of 0.6 nm of each other.

Fig. 7.

Crystal-like clusters formed by urea displaying two distinct internal structures. Form I (A) corresponds to the experimental urea crystal structure, whereas form II (B) exhibits an internal arrangement in which the C–O axis is oriented as in form I while the N–N axis has a different orientation. Both clusters are displayed in such a way that the C–O axis is orthogonal to the viewing plane. The crystal-like particle that has an internal structure as form I has solid/liquid interfaces that correspond to the {110} and {001} faces, typical of macroscopic crystals grown in water.

Fig. 8.

Simulation S2. (A) Contour plot of the FES obtained in the plane ($N_{\mathit{I}}$, $N_{\mathit{II}}$). (B) Free-energy profiles extracted from the FES for increasing cluster sizes and reported as relative to the pure form I nucleus. It can be observed that the free-energy difference between form I and form II decreases with cluster size, whereas for small nuclei form II is markedly favored, increasing in size form I and form II become almost equally probable.