Letter pubs.acs.org/JPCL
Reaction Coordinate for Ice Crystallization on a Soft Surface Laura Lupi, Rebecca Hanscam, Yuqing Qiu, and Valeria Molinero* Department of Chemistry, The University of Utah, 315 South 1400 East, Salt Lake City, Utah 84112-0850, United States S Supporting Information *
ABSTRACT: The control of assembly and crystallization of molecules is becoming increasingly important in chemistry, engineering, and materials sciences. Crystallization is also central to understand natural processes that include the formation of atmospheric ice and biomineralization. Organic surfaces, biomolecules, and even liquid/vapor interfaces can promote the nucleation of crystals. These soft surfaces present significant structural fluctuations, which have been shown to strongly impact the rate of crystallization. This raises the question of whether degrees of freedom of soft surfaces play a role in the reaction coordinate for crystal nucleation. Here we use molecular simulations to investigate the mechanism of ice nucleation promoted by an alcohol monolayer. Our analysis indicates that while the flexibility of the surface strongly depresses its ice nucleation ability, it does not play a role in the coordinate that controls the transformation from liquid to ice. We find that the variable that drives the transformation is the size of the crystalline cluster, the same as that for the homogeneous crystallization. We argue that this is a general result that arises from the separation of time scales between surface fluctuations and the crossing of the transition state barrier for crystallization.
T
assumed by CNT. To the best of our knowledge, the reaction coordinate for nucleation of any crystal on a flexible, soft surface has never been determined. Organic and biological surfaces are some of the most efficient ice nucleating agents. Among nonbiological surfaces, alcohol monolayers are particularly efficient in promoting heterogeneous nucleation of ice.7,14 Recent studies reveal that curvature, roughness, and structural fluctuations in ice nucleating surfaces strongly depress their ice nucleation efficiency.6,7 This demonstrates that the nucleation barrier depends on the flexibility of the surface. However, it does not answer whether changes in the curvature of the surface could propel the transformation from liquid to ice. The participation of such degrees of freedom in the reaction coordinate of heterogeneous nucleation would add a new variable in the calculation of free energy profiles from which the nucleation barrier is computed and require a substantial modification of CNT. In this work, we use molecular simulations with the mW water model,15 rare event sampling methods,16−18 and committor analysis19 to determine the reaction coordinate for heterogeneous nucleation of ice on a monolayer of tricontanol (C30H61OH). We find that the degrees of freedom of the flexible surface modulate the pathway of nucleation but do not play a role in the crossing of the transition state. In agreement with the assumption of CNT, the size of the ice cluster is the variable that predicts the fate of the crystal embryo along the nucleation pathway. We first investigate the molecular process of ice nucleation on the alcohol monolayer at 220 K, 19 K above the homogeneous nucleation temperature of the mW water
he reaction coordinate is the key variable that controls the crossing of the transition state of a chemical process.1 Catalysts typically change not only the barrier but also the mechanism of reactions.2 Crystallization can be considered a chemical process in which the reaction coordinate is a collective variable that describes the crossing of the transition state between the liquid and crystal. The nucleation of a crystal can be “catalyzed” by a surface that promotes the transformation, a process known as heterogeneous nucleation. Classical Nucleation Theory (CNT)3 is generally used to interpret rates from laboratory measurements and molecular simulations. CNT proposes that the reaction coordinate for crystal nucleation is the size of the nascent crystallite, irrespective of whether the process is assisted by a surface or not. This assumes that the promoting surface does not change the reaction coordinate of crystallization. Experiments and simulations indicate that surfaces alter the pathway from liquid to crystal.4−10 It is an open question whether this change in the pathway can result in a change of the variable that controls the crossing of the transition state. Water freezing is one of the most ubiquitous crystallization processes. In the atmosphere, the formation of ice mostly proceeds assisted by soot, mineral dusts, and organic and biological materials.11 Recent studies of ice nucleation by rigid graphitic surfaces demonstrate that the pathway for nucleation involves the early creation of preordered ice-like regions at the surface, where the ice nuclei are born.5,6 Nevertheless, statistical analysis of nucleation trajectories reveal that the reaction coordinate for ice crystallization on graphitic surfaces is the size of the crystalline nucleus,5,12 the same as that for the homogeneous nucleation of ice.13 These results indicate that the reaction coordinate for homogeneous nucleation and heterogeneous nucleation of ice on a rigid surface is the one © XXXX American Chemical Society
Received: July 18, 2017 Accepted: August 19, 2017 Published: August 19, 2017 4201
DOI: 10.1021/acs.jpclett.7b01855 J. Phys. Chem. Lett. 2017, 8, 4201−4205
Letter
The Journal of Physical Chemistry Letters
Figure 1. Snapshots along a crystallization trajectory of water in contact with a monolayer of tricontanol at 220 K and 0 bar. Upper panels show the periodic simulation box in a direction perpendicular to the monolayer, central panels show the top view of the interfacial layer with bonds (gray sticks) between hydroxyl groups of the monolayer (red balls) and the intercalated water (gray balls), and bottom panels show the topography of the interfacial layer (gray surface). The ice clusters are shown with blue sticks, liquid water with cyan points, carbon chains with dark gray balls, and defects in the interfacial layer are shown with green sticks. We consistently find, here and in ref 7, that ice nucleation occurs in flat, defect-free regions of the surface. Note that even rigid flat alcohol monolayers display dynamic lines of defects as these originate in reordering of the intercalated water to relieve lattice mismatch between the monolayer and ice. The size of the defect-free regions decreases with increasing mobility of the monolayer, resulting in a decrease in the ice nucleation temperature.7
model.20 At this temperature, the heterogeneous nucleation barrier can be overcome through spontaneous fluctuations. Figure 1 shows snapshots along the crystallization trajectory. Water molecules at the interface (light gray balls) intercalate between the hydroxyl groups of the alcohols (red balls), forming an interfacial layer with ice-like domains.7 There is a lattice mismatch between the hydroxyl groups in the monolayer and ice: the surface of the monolayer is 10.4% more extended in the vertical direction and 4.3% compressed in the horizontal direction with respect to mW ice,7 the same as that for the tricontanol monolayer in experiments.21 Despite the large mismatch, the monolayer templates the formation of ice, which is directly hydrogen bonded to the water molecules intercalated with the hydroxyl groups of the surface.7 The lattice mismatch to ice is relieved through the formation of defects (shown with green sticks) in the interfacial layer. The defects consist of additional intercalated water molecules in the interfacial layer, aligned in the direction of the hydrogen bonds between the intercalated waters and hydroxyl groups. The lines of defects protrude above the plane of the surface, creating wrinkles. (lower panel of Figure 1). This imparts instantaneous roughness and curvature to the surface, both detrimental for ice nucleation.6,7 Fluctuations of the alcohol and intercalated water molecules result in rearrangement of the interfacial layer during the induction time preceding nucleation (Figure 1a−d). Wrinkles and flat domains are created and annihilated within a few nanoseconds, that is, in time scales that are short compared to the induction time for crystal nucleation. The flat regions in the
interfacial layer offer defect-free hydrogen bonding sites to ice and are the loci of nucleation (Figure 1e−f). In what follows, we investigate whether changes in the curvature of the surface and the size of flat regions in the interfacial layer contribute to the crossing of the barrier from liquid to ice. A single simulation trajectory, such as that shown in Figure 1, provides important physical insight on the nucleation pathway but is insufficient to resolve the reaction coordinate and structure of the transition state. To determine the reaction coordinate for heterogeneous nucleation of ice on the soft organic surface, we first collect over 16 000 simulation trajectories of ice crystallization using aimless shooting transition path sampling.16−18 The simulations were performed at 235 K (33 K above the homogeneous freezing temperature of mW20) to produce critical nuclei that contained about 500 molecules (i.e., 10% of the water molecules in the simulation cell; see Methods in the Supporting Information). We then propose several candidate reaction coordinates, including various measures of the size of the crystalline cluster, shape, number of hydrogen bonds of the ice crystal to the surface, as well as the size of ice-like ordered and flat patches in the interfacial layer. We finally rank these candidates using maximum likelihood optimization,16 according to their ability to represent the ensemble of trajectories. We find that the best reaction coordinate is the size of the crystalline cluster identified with neighbor-correlated bondorder parameters through the CHILL+ algorithm, NCHILL+22 (Table 1), the same as for the homogeneous nucleation of ice.13 The size of the crystalline cluster measured with the 4202
DOI: 10.1021/acs.jpclett.7b01855 J. Phys. Chem. Lett. 2017, 8, 4201−4205
Letter
The Journal of Physical Chemistry Letters Table 1. Maximum Likelihood Scores and Values at the Predicted Transition State Surface for Candidate Reaction Coordinatesa candidate reaction coordinate NCHILL+ (ice molecules) Nq6>0.57 (ice molecules) Nflat‑domains (int. water + C3oOH) Imax/Imin Rg (Å) NH‑bonds (hydrogen bonds)
score 0 79 450 613 644 697
transition state 429 293 418 1.3 17 132
a
The lower the score, the better that the coordinate describes the transition state (see the Methods).
noncorrelated bond order parameter q6,23 Nq6>0.57, has a statistically significantly lower likelihood of describing the ensemble of trajectories. Not all measures of size of the ice crystal are equally good reaction coordinates. The size of ice-like flat domains formed by hydroxyl groups of the monolayer and intercalated waters (Nflat‑domains) and the number of hydrogen bonds between the ice cluster and the interfacial layer (NH‑bonds) rank poorly as reaction coordinates for the nucleation of ice on the alcohol monolayer (Table 1). That is also the case for the shape of the ice cluster measured by the radius of gyration (Rg) and the ratio of moments of inertia in the plane of the surface (Imax/Imin). A combination of pairs of any of the six variables of Table 1 does not yield a statistically significant improvement of the maximum likelihood score (Supporting Information Table S1). This implies that two-dimensional order parameters are not needed to describe the nucleation process. We conclude that the degrees of freedom of the surface and the shape of the crystallite cannot predict the crossing of the barrier from liquid to ice. To provide a stringent, absolute test of the accuracy of the size of the crystallite defined by NCHILL+ as a reaction coordinate, we calculate the committor histogram19 (Figure 2c) for an ensemble of 100 transition state configurations. Figure 2a displays randomly selected configurations of the transition state ensemble, corresponding to NCHILL+* = 429. The standard deviation of the committor distribution (σ = 0.117) compares very well with the one expected for the exact transition state surface (σ = 0.11; see Methods). The quantitative agreement in the committor distributions indicates that the size of the crystalline cluster is the variable that controls the crossing of the nucleation barrier. The same level of accuracy was previously found for the size of the crystal as a reaction coordinate for homogeneous ice nucleation13 and heterogeneous nucleation of ice on rigid graphitic surfaces.5,12 We conclude that, different from most chemical reactions, the reaction coordinate for ice crystallization is not affected by the surface that promotes the transformation. Analysis of the ensemble of nucleation trajectories indicates that ice-like flat domains in the interfacial layer are the birthplace of the subcritical crystallites that successfully reach the transition state. To understand why the degrees of freedom of the surface do not play a role in the reaction coordinate, we show in Figure 3 the average time that it takes for (i) the critical nuclei to dissolve and (ii) the flat surface domains on which the critical nuclei sit to shrink as the nuclei dissolve. We note that the time to dissolve a specific nucleus is identical to the time it takes for that nucleus to climb the nucleation barrier because both the algorithms used in the simulations and the path sampling trajectories that they produce are time-reversible.24
Figure 2. (a) Representative transition state configurations for the crystallization of ice (blue) by an alcohol monolayer. The topography of the interfacial layer, composed of intercalated water and hydroxyl groups, is shown in gray. (b) Distribution of the ratio of the largest to smallest component of the moments of inertia tensor in the plane of the surface (Imax/Imin) for the transition state ensemble of the crystalline nuclei measured with the NCHILL+ reaction coordinate. (c) Committor probability histogram (red line) calculated for sizes of the crystalline cluster within five molecules of the predicted transition state, NCHILL+* = 429. The black dashed line corresponds to the binomial distribution expected for the exact transition state ensemble (see Methods).
Figure 3. Evolution of the average size of the critical ice nuclei (blue symbols) and the underlying flat domains at the nucleating surface (red symbols) on their way toward the top of the nucleation barrier, which is reached at t = 0. Sizes of nuclei and their underlying surface domains are normalized to their values at the top of the barrier. The lines are fits of the data to stretched exponentials, ⟨size(t)/size(0)⟩ = exp(−(t/τ)β), where τ = 29.3 ns and β = 0.52 for the surface (red line) and τ = 0.95 ns and β = 0.89 for the crystal nuclei (blue line). This demonstrates that the flat domains of the surface wherein the nuclei form evolve slowly compared to the time that it takes for the nuclei to cross the barrier top.
We find that the time scale of disappearance of the flat surface domains is long (tens of nanoseconds) compared to the time 4203
DOI: 10.1021/acs.jpclett.7b01855 J. Phys. Chem. Lett. 2017, 8, 4201−4205
Letter
The Journal of Physical Chemistry Letters required for the critical nuclei to relax from the transition state (1 ns). The separation of time scales explains why surface fluctuations modulate the barrier height without changing the progress coordinate. Although the size of the flat, ice-like regions of the surface is not part of the reaction coordinate, these regions nevertheless play a key role in the overall pathway and rate of ice nucleation.7 The topography of the surface constrains the size and shape of the critical nuclei (Figure 2a), which sit in flat domains of the surface delimited by the wrinkles created by the lines of defects. The critical embryos can be highly anisotropic (Figure 2a): their distribution of the ratio of largest to smallest moments of inertia is bimodal, with peaks at 1.1 and 2.2 (Figure 2b), corresponding to less and more stretched nuclei, respectively. This heterogeneity arises from the distinct configurations of the surface with flat domains large enough to host a crystallite of critical size. Sampling over a larger number of independent initial configurations for the transition path sampling procedure may result in a broad unimodal distribution instead of a bimodal one. Nevertheless, we find that the critical size does not depend on the anisotropy of the critical crystallite. This confirms that the shape of the crystallites is not relevant for the successful crossing of the nucleation barrier. When postcritical crystallites grow across lines of defects, the alignment of the hexagons in the ice layer in contact with the interface transitions from out-of-registry to in-registry with respect to the interfacial layer. Out-of-registry corresponds to cubic ice-like order at the base of the crystallite, while inregistry corresponds to hexagonal ice-like order.25 The cost of the line of defects between cubic and hexagonal order in a plane is 1 ± 0.1 pN nm−1.26 This can impose a significant cost on small crystals. For example, crystallites bisected by a line of defects are transition states for transformation between hexagonal and cubic ice clusters at the ice/vapor interface.25 Likewise, we find that in heterogeneous ice nucleation from the liquid it is more favorable for the crystal to adapt its shape to follow the lines of defects of the surface than to grow across them. We conclude that the flexibility of the surface controls the rate of crystal nucleation through the probability of building a flat domain sufficiently large to host a critical nucleus. The results of this work suggest that surfaces that evolve slowly compared to the time scales involved in the crossing of the transition state barrier do not contribute to the reaction coordinate of crystallization. We expect this to be valid for nucleation of crystals on hard or soft solid surfaces. It would also hold for crystal nucleation at liquid/liquid27 and at liquid/ gas interfaces28 if their capillary waves with length scales comparable to those of the critical nuclei evolve slowly compared to the rate of attachment of new particles to the critical nucleus. Surface modes in heterogeneous crystallization could be considered as equivalent to promoter modes in catalysis, which are important in early stages of the reaction but do not determine the crossing of the transition state.29 Our analysis suggests that despite the richness in nucleation pathways observed for heterogeneous crystal nucleation, only the evolution of the free energy as a function of the size of the crystal embryo is needed to predict crystal nucleation rates.
■
■
Models, methods, and a table with ranking of pairs of variables (PDF)
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. ORCID
Valeria Molinero: 0000-0002-8577-4675 Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS We gratefully acknowledge discussions with Baron Peters and Michael Grünwald. This work was supported by the National Science Foundation through Award CHE-1305427 “Center for Aerosols Impacts on Climate and the Environment” and CHE1358740 “Catalysis in a Collaborative Research Experience for Undergraduates Program”. We thank the Center for High Performance Computing at the University of Utah for technical support and a grant of computer time.
■
REFERENCES
(1) Peters, B. Reaction Rate Theory and Rare Events; Elsevier: Amsterdam, The Netherlands, 2017. (2) Klippenstein, S. J.; Pande, V. S.; Truhlar, D. G. Chemical Kinetics and Mechanisms of Complex Systems: A Perspective on Recent Theoretical Advances. J. Am. Chem. Soc. 2014, 136, 528−546. (3) Turnbull, D.; Vonnegut, B. Nucleation Catalysis. Ind. Eng. Chem. 1952, 44, 1292−1298. (4) Jungblut, S.; Dellago, C. On the Reaction Coordinate for Seeded Crystallisation. Mol. Phys. 2015, 113, 2735−2741. (5) Lupi, L.; Peters, B.; Molinero, V. Pre-Ordering of Interfacial Water in the Pathway of Heterogeneous Ice Nucleation Does Not Lead to a Two-Step Crystallization Mechanism. J. Chem. Phys. 2016, 145, 211910. (6) Lupi, L.; Hudait, A.; Molinero, V. Heterogeneous Nucleation of Ice on Carbon Surfaces. J. Am. Chem. Soc. 2014, 136, 3156−64. (7) Qiu, Y.; Odendahl, N.; Hudait, A.; Mason, R.; Bertram, A. K.; Paesani, F.; DeMott, P. J.; Molinero, V. Ice Nucleation Efficiency of Hydroxylated Organic Surfaces Is Controlled by Their Structural Fluctuations and Mismatch to Ice. J. Am. Chem. Soc. 2017, 139, 3052− 3064. (8) Nielsen, M. H.; Lee, J. R. I.; Hu, Q.; Yong-Jin Han, T.; De Yoreo, J. J. Structural Evolution, Formation Pathways and Energetic Controls during Template-Directed Nucleation of CaCO3. Faraday Discuss. 2012, 159, 105−121. (9) De Yoreo, J. J.; Waychunas, G. A.; Jun, Y. S.; Fernandez-Martinez, A. In Situ Investigations of Carbonate Nucleation on Mineral and Organic Surfaces. Rev. Mineral. Geochem. 2013, 77, 229−257. (10) Mosses, J.; Turton, D. A.; Lue, L.; Sefcik, J.; Wynne, K. Crystal Templating through Liquid-Liquid Phase Separation. Chem. Commun. 2015, 51, 1139−1142. (11) Murray, B. J.; O’Sullivan, D.; Atkinson, J. D.; Webb, M. E. Ice Nucleation by Particles Immersed in Supercooled Cloud Droplets. Chem. Soc. Rev. 2012, 41, 6519−6554. (12) Cabriolu, R.; Li, T. Ice Nucleation on Carbon Surface Supports the Classical Theory for Heterogeneous Nucleation. Phys. Rev. E 2015. (13) Lupi, L.; Hudait, A.; Peters, B.; Grünwald, M.; Mullen, R. G.; Nguyen, A. H.; Molinero, V., Role of Stacking Disorder in Ice Nucleation. Nature, Under Review. (14) Gavish, M.; Popovitz-Biro, R.; Lahav, M.; Leiserowitz, L. Ice Nucleation by Alcohols Arranged in Monolayers at the Surface of Water Drops. Science 1990, 250, 973−975. (15) Molinero, V.; Moore, E. B. Water Modeled as an Intermediate Element between Carbon and Silicon. J. Phys. Chem. B 2009, 113, 4008−4016.
ASSOCIATED CONTENT
S Supporting Information *
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpclett.7b01855. 4204
DOI: 10.1021/acs.jpclett.7b01855 J. Phys. Chem. Lett. 2017, 8, 4201−4205
Letter
The Journal of Physical Chemistry Letters (16) Peters, B.; Trout, B. Obtaining Reaction Coordinates by Likelihood Maximization. J. Chem. Phys. 2006, 125, 054108. (17) Peters, B.; Beckham, G. T.; Trout, B. L. Extensions to the Likelihood Maximization Approach for Finding Reaction Coordinates. J. Chem. Phys. 2007, 127, 034109. (18) Mullen, R. G.; Shea, J.-E.; Peters, B. Transmission Coefficients, Committors, and Solvent Coordinates in Ion-Pair Dissociation. J. Chem. Theory Comput. 2014, 10, 659−667. (19) Bolhuis, P. G.; Chandler, D.; Dellago, C.; Geissler, P. L. Transition Path Sampling: Throwing Ropes over Rough Mountain Passes, in the Dark. Annu. Rev. Phys. Chem. 2002, 53, 291−318. (20) Moore, E. B.; Molinero, V. Structural Transformation in Supercooled Water Controls the Crystallization Rate of Ice. Nature 2011, 479, 506−508. (21) Popovitz-Biro, R.; Wang, J. L.; Majewski, J.; Shavit, E.; Leiserowitz, L.; Lahav, M. Induced Freezing of Supercooled Water into Ice by Self-Assembled Crystalline Monolayers of Amphiphilic Alcohols at the Air-Water Interface. J. Am. Chem. Soc. 1994, 116, 1179−1191. (22) Nguyen, A. H.; Molinero, V. Identification of Clathrate Hydrates, Hexagonal Ice, Cubic Ice, and Liquid Water in Simulations: The Chill+ Algorithm. J. Phys. Chem. B 2015, 119, 9369−9376. (23) Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M. BondOrientational Order in Liquids and Glasses. Phys. Rev. B: Condens. Matter Mater. Phys. 1983, 28, 784−805. (24) Rowley, C. N.; Woo, T. K. New Shooting Algorithms for Transition Path Sampling: Centering Moves and Varied-Perturbation Sizes for Improved Sampling. J. Chem. Phys. 2009, 131, 234102. (25) Hudait, A.; Molinero, V. What Determines the Ice Polymorph in Clouds? J. Am. Chem. Soc. 2016, 138, 8958−8967. (26) Hudait, A.; Qiu, S.; Lupi, L.; Molinero, V. Free Energy Contributions and Structural Characterization of Stacking Disordered Ices. Phys. Chem. Chem. Phys. 2016, 18, 9544−9553. (27) Collier, K. N.; Brooks, S. D. Role of Organic Hydrocarbons in Atmospheric Ice Formation Via Contact Freezing. J. Phys. Chem. A 2016, 120, 10169−10180. (28) Qiu, Y.; Molinero, V. Strength of Alkane−Fluid Attraction Determines the Interfacial Orientation of Liquid Alkanes and Their Crystallization through Heterogeneous or Homogeneous Mechanisms. Crystals 2017, 7, 86. (29) Peters, B. Transition-State Theory, Dynamics, and Narrow Time Scale Separation in the Rate-Promoting Vibrations Model of Enzyme Catalysis. J. Chem. Theory Comput. 2010, 6, 1447−1454.
4205
DOI: 10.1021/acs.jpclett.7b01855 J. Phys. Chem. Lett. 2017, 8, 4201−4205
