Atmospheric methane levels increased over geological formations in Northern Siberia following a heat wave in 2020, suggesting that permafrost thaw releases methane from reservoirs.
Image credit: Wikimedia Commons.
Edited by Alexis T. Bell, University of California, Berkeley, CA, and approved October 16, 2019 (received for review June 21, 2019)
Significance
Solution crystallization is a widely used technique to grow crystalline materials such as pharmaceuticals, agrochemicals, catalysts, semiconductors, and metal–organic frameworks. The highly probabilistic nature and complex molecular processes spanning multiple lengths and time scales have prevented full mechanistic understanding of crystallization processes for over a century. Here we show fluctuations in the solvation shell are key molecular events that explain and unify the empirical observations laid down since the inception of crystallization research. The activation barrier of the growth processes obtained from these molecular fluctuations can be used to estimate growth rates that show reasonable agreement with the experimental observations, hinting at the rise of additional and efficient strategies to control solution crystallization.
Abstract
Solution crystallization is a common technique to grow advanced, functional crystalline materials. Supersaturation, temperature, and solvent composition are known to influence the growth rates and thereby properties of crystalline materials; however, a satisfactory explanation of how these factors affect the activation barrier for growth rates has not been developed. We report here that these effects can be attributed to a previously unrecognized consequence of solvent fluctuations in the solvation shell of solute molecules attaching to the crystal surface. With increasing supersaturation, the average hydration number of the glutamic acid molecule decreases and can reach an asymptotic limit corresponding to the number of adsorption sites on the molecule. The hydration number of the glutamic acid molecule also fluctuates due to the rapid exchange of solvent in the solvation shell and local variation in the supersaturation. These rapid fluctuations allow quasi-equilibrium between fully solvated and partially desolvated states of molecules, which can be used to construct a double-well potential and thereby to identify the transition state and the required activation barrier. The partially desolvated molecules are not stable and can attach spontaneously to the crystal surface. The activation barrier versus hydration number follows the Evans–Polanyi relation. The predicted absolute growth rates of the α-glutamic acid crystal at lower supersaturations are in reasonable agreement with the experimental observations.
Crystallization is one of the most common processes used for the separation of molecules from solution and for the synthesis of crystalline materials for various applications such as catalysis, gas separation, and energy conversion (1, 2). The underlying processes of crystallization, namely nucleation and growth, govern the structure, morphology, and size of crystals, which impart specific properties to the crystalline materials (3⇓–5). Supersaturation, defined as the chemical potential of solution relative to the state when solution is in equilibrium with the crystals, is generally regarded as the driving force for crystallization processes. Several empirical and semiempirical kinetic models have been developed based on this hypothesis and validated experimentally to relate nucleation and growth rates to supersaturation (2, 6⇓–8). Examples of empirical models for nucleation rate include classical nucleation theory (9), 2-step nucleation theory (10), and probabilistic theory (9) and for growth rates are the Bravais–Friedel–Donnay–Harker (BFDH) method (11), the Hartman–Perdok theory (12⇓–14), the Burton–Cabrera–Frank (BCF) model (15), the birth and spread model (16), the bulk diffusion model (17), the thermochemical model (18), and MONTY (19). While the majority of these models rely on parameters estimated from the experimental data (16), the other models are limited to only the prediction of relative rates or growth shapes (20, 21). The molecular simulation methods such as molecular dynamics (MD) (22), kinetic Monte Carlo (kMC) (23, 24), and hybrid MD-kMC (25) have also been applied to simulate the growth of the crystal. However, limitations due to variation in supersaturation (22), absence of desolvation kinetics (22, 24), influence of diffusion-limited growth (23), and computational complexity (25) have made it difficult to obtain intrinsic growth rates from such molecular simulations. Currently, there are no theoretical models that can predict intrinsic rates of nucleation and face-specific growth from the first principles. The present state is due to a limited consideration of the underlying physics relating supersaturation to the rates of crystallization and integrating them into the existing models. The key physical processes in crystallization include self-assembly of molecules in the solution (nucleation) as well as on the facet of the crystal (growth), both of which are influenced by solvents, additives, pH, and temperature of a supersaturated solution (26, 27). The modeling of the rates of nucleation and growth requires detailed consideration of such self-assembly processes to estimate the transition state and the activation barrier.
Fig. 1 shows a typical process of solvent-directed self-assembly of small molecules on a growing facet of a crystal. This process initiates in a supersaturated solution where the solute molecules in the solution diffuse toward the crystal surface due to the chemical potential gradient. As the solvated molecules move toward the crystal surface, the intermolecular interaction between solute molecules increases, causing the solvation layer to deplete. The partially solvated (or desolvated) solute molecules then integrate spontaneously into the “kink” sites on the crystal surface by shedding the remaining solvation layer. The integration step is kinetically irreversible as crystals only grow and do not dissolve in a supersaturated solution. The partially solvated state of the molecule can be considered as a transition state, which determines the activation barrier and hence the rate of integration of solvated molecules into the kink sites of the crystal surface (28). The formation of the partially solvated state is the rate-limiting step, which is energetically uphill such that the fully and partially solvated states of the molecule are in dynamic equilibrium. Similar arguments may be applied to the process of nucleation, where the fluctuations in supersaturations (7) lead to the instantaneous formation of partially solvated molecules which aggregate spontaneously to form subcritical clusters. The desolvation of the ion/molecules has been previously identified as a governing step for crystal growth (29, 30), where the formation of the partially desolvated ion/molecule is the rate-limiting step (31). Joswiak et al. (28) have also determined the partially desolvated state of Na+/Cl− ions to predict the activation barrier for the integration of ions into the kink sites of the NaCl crystal. However, the estimation of the activation barrier for solvated molecules to integrate or attach to the crystal surface and thereby the prediction of absolute growth rates remain the outstanding issues in crystallization.
Fig. 1.
Mechanism of integration of molecules on the growing face of a crystal. The solvation shell of the glutamic acid molecule in the supersaturated solution (GLU(l)) depletes as it moves toward the crystal surface. This partially solvated molecule (GLU‡(l)) then integrates spontaneously to the kink site on the crystal surface. The potential energy diagram shows the fully solvated GLU(l) in equilibrium with the energetically uphill partially solvated GLU‡(l), which then integrates irreversibly to the kink site on the crystal surface (GLU(s)).
Two different approaches have been used to predict growth rates from molecular simulations: 1) hybrid MD-kMC simulations and 2) rare event simulations. The first method involves MD simulation of the solid–liquid interface in supersaturated solution to determine the rate of addition of solute to the solid interface and hence obtain the rate constants, which are then used as transition probabilities in kMC simulation to predict growth rates and growth shapes of crystals (22, 25). These approaches rely on either using simplified force fields for ideal particles or approximating the rate constant for nonequilibrium processes by counting events in the equilibrium simulations. The hybrid MD-kMC simulation has been implemented to predict either the growth shapes (25) or the rate constants (22), but not the growth rates. The rare event simulations use biased-MD methods by applying an external force to the simulation box and use methods such as umbrella sampling that can sample the crystallization coordinate while maintaining the information on conformational changes during the attachment of molecules and calculating the energy barrier (28). Due to the stochastic nature of crystal growth, the identification of the sampled path and the associated energy barrier that leads to the growth of the crystal is not trivial. The other rare event simulation method is forward flux sampling, which requires prior information on the terminal states such as the solvated molecule in the bulk and the desolvated molecule at the kink site to determine blind pathways. The complexity and longer computational time involved in such rare event simulations have limited its application to systems such as ionic crystals for prediction of the rate constants for ion attachments (28), clathrate hydrates to understand nucleation (32), and biomolecules to study conformational changes (33). None of these molecular simulation methods have been extended to predict intrinsic, absolute growth rates of organic crystals from first principles. Here we present an MD simulation method to estimate the energy barrier for partial desolvation of solute molecules before their integration into the kink sites. The average configuration of the solvation layer of solute molecules in a supersaturated solution was used to construct a double-well potential, similar to the method of Shustorovich (34, 35), that yields the transition state as well as the activation barrier. The growth rates were estimated using the established thermochemical models (21) that account for the mechanism of growth, energetics of kink sites, and the face-specific growth factors.
Theory
Multiscale computations were performed to obtain the activation barrier for attachment of glutamic acid molecules on the facets of the α-glutamic acid crystal and to calculate face-specific growth rates. The computation scheme was organized in 3 steps: 1) MD simulations of the supersaturated solution to characterize the solvation shell dynamics; 2) molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) calculation to obtain the binding energy of water molecules in the solvation shell, followed by the estimation of the activation barrier using the double-well approach; and 3) Hartman–Perdok simulation to obtain the slice energy of the crystal face and the estimation of absolute face-specific growth rates using mechanistic growth rate models. The overview of the computational approaches and the details of the calculations performed are provided in SI Appendix.
MD Simulations.
GROMACS code (36) was used to perform the molecular dynamics simulations. A dodecahedron box of fixed volume 216.48 nm3 was set up and centered around the origin. To study the molecular events at high supersaturation, some of the water molecules were randomly replaced with solute molecules (22). Four different systems of supersaturation ratios, 1.2, 1.5, 2, and 2.5, were set up on the basis of experimentally measured solubility of α-glutamic acid crystals in water (37⇓–39). The systems were simulated using the Optimized Potential for Liquid Simulations–All Atom (OPLS-AA) force field for 5 ns with a time step of 2 fs. The frames and energies were obtained at every 1 ps. Before the production run in the isothermal–isobaric (NPT) ensemble, systems were equilibrated for 1 ns each in canonical (NVT) and NPT ensembles, making the total simulation time equal to 7 ns (40). After the simulations, the distances between glutamic acid molecules were obtained and the radial distribution function (RDF) was plotted for molecules that spend the highest time in close vicinity to the other glutamic acid molecules. The RDF was plotted with the center of mass of the individual glutamic acid molecule as a reference point to analyze the distribution of water molecules around itself. The first peak of the RDF was used to obtain the solvation shell thickness at each supersaturation. Further, to obtain the energy of interaction between the glutamic acid molecule and each water molecule in the solvation shell, the MM-PBSA method was applied to decompose the energy contribution on a per-molecule basis using the Poisson–Boltzmann algorithm (41). The additional details of molecular simulations including setting up of the simulation box, the MD simulation procedure, calculation of the number of water molecules in the solvation shell (hydration number) and solvation shell thickness, plotting of the RDF, and calculation of binding energy of water molecules in the solvation shell are provided in SI Appendix.
Calculation of Activation Barrier and Growth Rates.
Correlating the distance between glutamic acid molecules and the hydration number allows classification of the structure of the solvated glutamic acid molecule at a certain time in the trajectory as either fully solvated or partially desolvated (SI Appendix). A frame in the trajectory was classified as “fully solvated” if the distance between one glutamic acid molecule and other glutamic acid molecule was greater than twice the solvation shell thickness. All of the other frames, where distance was shorter than twice the solvation shell thickness, were classified as “partially desolvated.” Since all of the simulations were performed without any external force, frames showing the extreme cases in each scenario were considered and were assumed to be in quasi-equilibrium. The structure with the highest hydration number of all of the fully solvated structures and the structure with the lowest hydration number of all of the partially desolvated structures were further analyzed to calculate the activation barrier for attachment. The fully solvated structure is denoted as state A of glutamic acid, the partially desolvated structure as state B, and the fully desolvated structure integrated into the lattice as state C. The transition of state B to state C is difficult to estimate using the molecular simulation approaches while keeping the supersaturation constant. However, the energy of the state C can be approximated from the lattice energy of the α-glutamic acid crystal (SI Appendix). To obtain the activation barrier for the transition of state A to state B, the interaction energy profiles or potential wells of 2 configurations representing state A and state B were obtained as a function of distance, and were superimposed to obtain the double-well potential. The potential well for the state B was separated by the amount of energy required to remove the respective number of water molecules from the state A. The intersection point of 2 potential wells determines the transition state and the height of the intersection gives the activation barrier. The schematic details of this approach are provided in SI Appendix. The intermolecular potentials of the OPLS-AA force field with Lorentz–Berthelot mixing rules were used to calculate the potential wells.
The calculation for face-specific growth rates (42) was performed in 3 steps: 1) Hartman–Perdok simulations to obtain stable flat faces (F-faces) and their attachment energies, 2) estimation of spiral morphology (step height and interstep distance) on F-faces, and 3) prediction of absolute growth rate using the mechanistic growth rate model. The details of growth rate calculations using the estimated activation barriers are described in SI Appendix.
Results and Discussion
Fig. 2A shows the variation in the RDFs of water molecules around a glutamic acid molecule with increasing supersaturation from 1.2 to 2.5. The height of the RDF peaks decreases with increasing supersaturation due to dehydration for the primary and secondary solvation shells. Such a phenomenon of dehydration of the solvation shell has been observed for simple ions, which is an effect of higher interionic interactions in concentrated solutions (43). Fig. 2B shows a similar behavior of decrease in the hydration number of glutamic acid molecules with increasing supersaturation, which is attributed to increasing interactions between glutamic acid molecules. At a sufficiently higher supersaturation, the average hydration number may reach an asymptotic limit corresponding to the number of adsorption sites on a molecule (43). The solvation shell also imposes a physical barrier for attachment of solute molecules to either the crystal face or the nucleus. This physical barrier can be anisotropic as the hydration numbers are different around each functional group of the molecule. It can be seen that the effect of supersaturation is to lower the average hydration number and thereby reduce the physical barrier for crystallization. However, the rates of crystallization are governed by fluctuations in the solvation shell that lead to desolvation and integration of molecules to the crystal surface at a fixed supersaturation.
Fig. 2.
(A) RDFs of water molecules at different supersaturations around the center of mass of a glutamic acid molecule. (B) Variation in the average hydration number of the primary solvation shell obtained from RDFs with increasing supersaturation.
Fig. 3 shows the decrease in the average hydration number (solid blue line) of a glutamic acid molecule with decreasing intermolecular distance (solid red line) between glutamic acid molecules at a supersaturation of 2.5. The method to identify the solvation shell and the dynamics of the hydration number at different supersaturations are given in SI Appendix. The lower average hydration number is observed when the relative distance between molecules is lower and the intermolecular interactions are maximum. The hydration number not only varies with the relative motion of glutamic acid molecules but also fluctuates (solid black line) due to the random motion of solvent molecules. The variation of hydration number over a longer time scale (solid blue line) is mostly due to variations in relative distances of glutamic acid molecules, which are dictated by local fluctuations in concentration or supersaturation. However, the fluctuations of hydration number at a much shorter time scale (solid black line) are due to the rapid exchange of water molecules in the solvation shell, which is governed by the bulk temperature. It can be seen that both temperature and supersaturation affect the variation in the hydration number and thereby the physical barrier for crystallization. The fully solvated and partially desolvated states of the molecule can be determined based on the highest and the lowest hydration number, respectively, as shown in Fig. 3. Since the fluctuations in hydration number due to temperature are much faster, the fully solvated and partially desolvated states can be considered in quasi-equilibrium. According to the transition state theory shown in Fig. 1, the equilibrium between fully solvated and partially desolvated molecules determines the activation barrier.
Fig. 3.
Dynamics of hydration number of the glutamic acid molecule at a supersaturation of 2.5. The running average of the hydration number (solid blue line) decreases with decrease in the intermolecular distance of the glutamic acid molecule (solid red line), which is governed by fluctuations in concentration. The instantaneous fluctuations in the hydration number (solid black line) are due to exchange of water molecules in the solvation shell, which increases with increase in temperature. The shorter time-scale dynamics of hydration number are dictated by temperature, whereas the longer time-scale dynamics are governed by local fluctuations in supersaturation.
Fig. 4A shows the energy profile for the interaction between glutamic acid molecules with increasing crystallization coordinate. The normalized crystallization coordinate spans the relative distance between molecules as they travel from the bulk solution (value 0) to the crystal face (value 1). The conformation of molecules does not change with increasing crystallization coordinate for small molecules such as glutamic acid, which was confirmed from molecular simulations. As the crystallization coordinate increases, the interaction between fully solvated molecules increases, which leads to desolvation of the primary shell and formation of partially desolvated molecules before integrating into the crystal face. The energy required to remove water molecules from the solvation shell is higher for the water molecules in the vicinity of the glutamic acid molecule (SI Appendix). The energy of transition, where the fully solvated molecules switch to a partially solvated configuration, can be obtained from the double-well diagram (SI Appendix) (34, 35). The activation energy for forming partially desolvated molecules decreases from 26 kJ⋅mol−1 to 17 kJ⋅mol−1 with increasing supersaturation from 1.2 to 2.5. Since the transition state is attained long before the integration of molecules to the crystal face, the activation energy is less sensitive to the relative orientation of molecules. The energy profiles for glutamic acid molecules oriented in [110] and [111] crystallographic directions and attaching to [110] and [111] crystal faces, respectively, are shown in SI Appendix. Since the integration of partially desolvated molecules to the crystal face is an irreversible process, it does not affect the activation barrier and the growth rate. A hypothetical pathway from the transition state to the crystal state is shown in Fig. 4A (dashed black line).
Fig. 4.
(A) Energy profiles for the integration of glutamic acid molecules on a crystal face at different supersaturations. The crystallization coordinate represents relative distances between molecules as they move from the bulk phase (value 0) to the crystal surface (value 1). The transitions of the fully solvated molecule to the partially desolvated molecule shown as solid curves determine the activation barrier. The dashed lines represent the hypothetical energy pathway for integration of partially desolvated molecules on the crystal surface. (B) A linear relationship between the activation barrier to form partially solvated molecules and the hydration number of fully solvated molecules for different supersaturations, which follows the Evans–Polanyi principle.
The heat of integration is the energy difference between the initial state (crystallization coordinate = 0.00, fully solvated molecules in solution) and the final state (crystallization coordinate = 1.00, molecules in the crystal lattice), which is also the sum of the energy of the crystal face (also referred to as a slice energy) and the hydration energy. Since the energy of the crystal face is constant, the heat of integration is proportional to the hydration energy or the hydration number. Fig. 4B shows a linear relationship between the activation barrier and the hydration number of fully solvated glutamic acid molecules at different supersaturations. This linear relationship is similar to the Evans–Polanyi relationship between the activation barrier and the heat of reaction.
Fig. 5 A and B shows the experimentally measured (37) and the predicted growth rates of [111] and [011] faces, respectively, of the α-glutamic acid crystal. The experimental growth rates obtained for the supersaturation range of 1.1 to 1.7 were extrapolated based on the fitting equations provided in ref. 37 up to the supersaturation of 2.5 to extend the comparison between simulated and experimental data. The growth rates were calculated using the supersaturation-dependent activation barriers (shown in Fig. 4A) in the established growth rate model. The predicted values of growth rates show a good agreement with the experimental growth rates for lower supersaturations. The predicted and experimental steady-state morphologies of the α-glutamic acid crystal are also in good agreement (SI Appendix). For higher supersaturation, the predicted growth rates are higher than the experimental growth rates, which can be due to either overestimation of activation barriers or change in growth mechanism from spiral to 2D nucleation with increasing supersaturation. The overestimation of the activation barriers might be due to the absence of crystal surface in the MD simulation. The higher interactions of fully solvated molecules near the crystal surface can result in higher desolvation and lower hydration number of partially desolvated molecules and thereby lower activation barrier. However, the solvation dynamics obtained from such MD simulations involving the crystal–solution interface are often affected by the continuous depletion of supersaturation (22). Fig. 5C shows the deviation of predicted growth rates from the experimental growth rates for [111] and [011] faces of the α-glutamic acid crystal. The extrapolated experimental growth rates for supersaturation of 2.5 are outliers as the rough growth is most likely the mechanism at such a high supersaturation, which can increase the predicted growth rate by an order of magnitude. Therefore, the predicted growth rates deviate significantly from the experimental growth rates at a supersaturation of 2.5.
Fig. 5.
Comparison between experimentally measured growth rates (solid triangles) and ab initio predicted growth rates (solid squares) for (A) the [111] face and (B) the [011] face of the α-glutamic acid crystal as a function of supersaturation. The extrapolated experimental growth rates are shown as open triangles (37) and the fitted simulation data are shown as a solid black line. (C) Deviation of the predicted growth rates from the experimental growth rates for a family of faces (011) and (111).
Conclusion
Supersaturation is generally regarded as a driver for crystallization. Here we show that the actual driving force for crystallization corresponds to the energy required to desolvate molecules, which in turn is a function of supersaturation. The hydration number and solvation shell can be considered as the energy barrier and physical barrier, respectively, for solute molecules to attach to the crystal surface. It is possible to control the solvation shell and hydration number of specific functional groups in a solute molecule by manipulating the solvent composition such that attachment of solute molecules in a certain crystallographic direction is favored. The average hydration number of a solute molecule in the solution can be reduced by increasing its concentration at a constant temperature or increasing the supersaturation. By reducing the average hydration number, the energy barrier to desolvate the solute molecules can be lowered.
The removal of solvent molecules from the fully solvated molecules is an energetically uphill process. When a critical number of solvent molecules are removed, the partially desolvated molecule can spontaneously desolvate and attach to a crystal surface via energetically favorable noncovalent interactions. The partially desolvated state has been previously identified as a transition state for molecular attachment processes (29, 31). The energy required to desolvate a fully solvated molecule to form a partially desolvated molecule can be considered as the activation barrier for crystal growth. This activation barrier is a function of supersaturation, temperature, and solvent composition, which can be determined from the fluctuations in the solvation shell. The fluctuations in the hydration number are driven by 1) rapid solvent exchange in the solvation shell, which is a function of temperature, and 2) local variation in concentration, which is a function of supersaturation. These fluctuations cause instantaneous formation of fully and partially solvated molecules that are in quasi-equilibrium. The transition state and activation barrier can be obtained from the construction of a double-well potential using the average configuration of fully solvated and partially desolvated molecules. The activation barrier was linearly dependent on the hydration number, which is similar to the Evans–Polanyi relationship. Since the activation barrier is directly proportional to the average hydration number and the hydration number is inversely proportional to the supersaturation, it can be seen how the activation barrier for crystal growth decreases with increasing supersaturation.
The supersaturation-dependent activation barriers can be used in the transition state theory to obtain rate constants for growth, which can be supplied to the mechanistic growth model to predict absolute growth rates. The predicted values of growth rates were in reasonable agreement with the experimental growth rates of the α-glutamic acid crystal for lower supersaturations, which confirms that the predicted supersaturation-dependent activation barriers are physically meaningful. The larger disagreement between experiment and theory at higher supersaturation can be due to either the variation in the growth mechanism or overestimation of activation barriers. The prediction of the activation barrier can be improved by considering a crystal–solution interface at a constant supersaturation in MD simulations. These results confirm that the solvent fluctuation in the solvation shell is a key process that governs the solution-based synthesis of crystalline materials.
Data Availability.
All the computational data generated and associated simulation details are provided in the main text and SI Appendix.
Acknowledgments
This material is based on the work performed in the Materials and Systems Engineering Laboratory at the University of Illinois at Chicago. A.V.D. thanks Ms. Sneha M. Dixit from Max Planck Institute of Colloids and Interfaces for stimulating discussion on molecular simulations. We are grateful to our funding source, the US National Science Foundation (Award CBET-1706921).
Footnotes
- ↵1To whom correspondence may be addressed. Email: mrsingh{at}uic.edu.
Author contributions: M.R.S. designed research; A.V.D. performed research; A.V.D. and M.R.S. analyzed data; and A.V.D. and M.R.S. wrote the paper.
The authors declare no competing interest.
This article is a PNAS Direct Submission.
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1910691116/-/DCSupplemental.
Published under the PNAS license.
References
Solvent fluctuations in the solvation shell determine the activation barrier for crystal growth rates
Article Classifications
- Physical Sciences
- Engineering
Sign up for Article Alerts
Jump to section
You May Also be Interested in
Pottery, burnt clay, and burnt flint from Jordan dated from 7752 BCE to 5069 BCE yielded a record of variations in Earth’s magnetic field strength during the Neolithic period.
Image credit: Thomas E. Levy.
Mechanics of ant tunneling could inspire the development of improved robotic mining techniques.
Image credit: Jose E. Andrade and David R. Miller (California Institute of Technology, Pasadena, CA).
When it comes to potential geoengineering initiatives, researchers and policymakers need to know what to study—and when to stop.
Image credit: Shutterstock/Venera Salman.
The gasdermin protein family has unexpectedly ancient origins.
Image credit: Alex Johnson.