Simulations of Ice Nucleation by Model AgI Disks and Plates - The

Feb 15, 2016 - ... Christoph Dellago , Laura Felgitsch , Janine Fröhlich-Nowoisky , Hartmut Herrmann , Swetlana Jungblut , Zamin Kanji , Georg Menzl ...
0 downloads 0 Views 4MB Size
Download PDF
Article pubs.acs.org/JPCB

Simulations of Ice Nucleation by Model AgI Disks and Plates Stephen A. Zielke, Allan K. Bertram, and G. N. Patey* Department of Chemistry, University of British Columbia, Vancouver, British Columbia V6T 1Z1, Canada S Supporting Information *

ABSTRACT: Silver iodide is one of the most effective ice nuclei known. We use molecular dynamics simulations to investigate ice nucleation by AgI disks and plates with radii ranging from 1.15 to 2.99 nm. It is shown that disks and plates in this size range are effective ice nuclei, nucleating bulk ice at temperatures as warm as 14 K below the equilibrium freezing temperature, on simulation time scales (up to a few hundred nanoseconds). Ice nucleated on the Ag exposed surface of AgI disks and plates. Shortly after supercooling an ice cluster forms on the AgI surface. The AgI-stabilized ice cluster fluctuates in size as time progresses, but, once formed, it is constantly present. Eventually, depending on the disk or plate size and the degree of supercooling, a cluster fluctuation achieves critical size, and ice nucleates and rapidly grows to fill the simulation cell. Larger AgI disks and plates support larger ice clusters and hence can nucleate ice at warmer temperatures. This work may be useful for understanding the mechanism of ice nucleation on nanoparticles and active sites of larger atmospheric particles.

I. INTRODUCTION Ice formation can occur by both homogeneous and heterogeneous nucleation.1 Homogeneous nucleation refers to ice nucleation in the absence of a foreign surface, while heterogeneous nucleation refers to ice nucleation at a surface. Heterogeneous nucleation of ice occurs at temperatures that are much warmer than the homogeneous nucleation temperature of ice2 (approximately −38 °C in cloud droplets). Both homogeneous and heterogeneous nucleation require the formation of a critical cluster of ice. Clusters larger than the critical size grow into bulk ice, whereas subcritical clusters may either grow to critical size or melt.1 In heterogeneous nucleation, a surface facilitates ice nucleation by reducing the free-energy barrier to critical cluster formation. Experimentally, a variety of materials have been shown to be effective heterogeneous ice nuclei.2 Computer simulations are rapidly emerging as a valuable tool in the ongoing effort to understand the microscopic mechanisms of ice nucleation. Homogeneous nucleation has been investigated with atomistic water models by immersing “prefabricated” ice clusters (potential critical nuclei) into supercooled water,3−6 and by direct simulation of spontaneous nucleation for the coarse grained mW water model.7 These studies highlight that it is difficult to create critical clusters homogeneously. At higher temperatures, the probability of nucleation is extremely low, and at lower temperatures the slow dynamics of atomistic water models inhibits the observation of ice nucleation and growth on currently practical simulation time scales. If a good ice nucleus (IN) can be found, then heterogeneous nucleation poses an easier task for direct simulation, and spontaneous nucleation followed by ice growth has been observed on infinite sheets and disks of graphene8 and on infinite sheets of silver iodide.9,10 It has also been shown that ice nucleation can be influenced by sheets of kaolinite.11 Cox et al.12,13 have recently reported simulations of ice nucleation on nanoparticles and disks employing the mW water model.7 They © 2016 American Chemical Society

found that the ability of a surface to organize water molecules into an ice-like configuration, controlled by the surface morphology and the water−surface interaction energy, can have a large effect on the ice nucleation rate, giving significant increases relative to the homogeneous ice nucleation rate in most cases. Silver iodide is a very effective ice nucleating agent,14 nucleating ice at temperatures as warm as −3 °C.2 Because of its efficiency for nucleating ice, AgI has been employed as a cloud seeding agent.15,16 Two polymorphs of silver iodide can exist under atmospheric conditions. The β phase is the most stable and has a hexagonal unit cell.17,18 The γ phase is metastable and has a cubic unit cell.17 Most simulation work has focused on the β-AgI(0001) face due to its hexagonal symmetry.19−28 Ice is similar to AgI in that a stable hexagonal form (Ih) and a metastable cubic form (Ic) can occur under atmospheric conditions.29 Both ice structures consist of layers of water molecules arranged in hexagonal, chair-conformed rings that provide a good match to the silver exposed β-AgI(0001) face.9 Ice structures Ih and Ic represent two different ways in which these layers can be stacked. The common layer forms the basal face (0001) of Ih and the (111) face of Ic. Because of this common layer, the occurrence of random layers of Ih and Ic in a single ice crystal is observed in both simulations30 and experiments,31,32 and is referred to as stacking disordered ice, Isd33 Previous simulation investigations of ice nucleation on infinite slabs of silver iodide9,10 have shown that ice nucleates readily on the silver exposed β-AgI(0001) face but not on the corresponding iodide exposed face. These simulations suggest that ice is not so much nucleating on favorable AgI surfaces, as Received: July 9, 2015 Revised: January 16, 2016 Published: February 15, 2016 2291

DOI: 10.1021/acs.jpcb.5b06605 J. Phys. Chem. B 2016, 120, 2291−2299

Article

The Journal of Physical Chemistry B

ions. For the silver-exposed β-AgI(0001) surface considered here we found that a wide range of potentials corresponding to partial charges ranging from 0 to ±0.8e nucleated ice efficiently. Figure S9 of the Supporting Information shows that when the induction term is included in the Hale-Keifer force field the potential continues to lie within (or very nearly so) the range of potentials for which we observe ice nucleation. For the AgI disks considered in the present paper, we also demonstrate that disks of uncharged “ions” nucleate ice just as effectively as those with partial charges of ±0.6e, the charges used in the HaleKeifer force field (see Figure S7). Our observations indicate that the ice nucleating ability of the model AgI surfaces is not strongly dependent on details of the electrostatic water−AgI interactions. As discussed in detail in ref 9, the atomistic geometry of the surface is the important factor. Finite-size disks and plates of different shape are “carved” from four layer thick sheets of β-AgI,18 such that a (0001) face forms the top (silver exposed) and bottom (iodide exposed) of each disk or plate. Top views of the disks and plates considered are shown in Figure 1. Note that all disks and plates are four layers thick (see Figure 2 for a side view) and are electrically

it is simply growing from surfaces with geometric features that closely resemble faces of ice. In this paper, we investigate ice nucleation on AgI disks and plates with nanometer dimensions. This investigation provides insight into the mechanism of ice nucleation on nanoparticles. In addition our simulations allow us to directly explore the size, nature, and relevance of critical clusters during heterogeneous nucleation. We note that such free-standing AgI disks and plates may not exist in real physical situations; however, it may be possible for a larger silver iodide particle to have small patches that are of similar size and structure to those studied here. Hence, our studies may provide insight into the ice nucleation properties of larger silver iodide particles. More generally, atmospheric particles, which typically range from 10 nm to 10 μm in size,1 may present many different crystal faces. As a result the surface of atmospheric particles may also contain crystalline patches with nanometer dimensions similar to those studied here. For convenience we refer to these patches as active sites. Our studies may also be useful for understanding the mechanism of ice nucleation on active sites of this type. We observe that ice nucleation on disks and plates with nanometer dimensions occurs in two stages. Initially, several layers of ice form on the disk or plate surface. These layers fluctuate growing and shrinking as the simulation proceeds, but the surface ice cluster never melts, and the number of “permanently” frozen water molecules can represent a significant fraction of the required critical nucleus. For a particular disk or plate, if the temperature is low enough, an ice cluster fluctuation eventually matches the critical size, and bulk ice nucleates and rapidly grows to fill the simulation cell. As the temperature is increased, larger disks are needed to create the larger critical clusters that are required to achieve nucleation. Wherever it is possible and meaningful, we compare our estimates of critical cluster size with theory, previous simulations, and experiment. The remainder of this paper has three sections. The models and simulation methods are described in Section II, results are given and discussed in Section III, and our conclusions are summarized in Section IV.

II. SIMULATION METHODS Molecular dynamics simulations were carried out (usually under NVT conditions) employing the GROMACS 4.5.5 and 4.6.2 software packages.34 The leapfrog integrator was used to propagate the dynamics with a time step of 2 fs. Coulombic interactions were handled using the smooth particle mesh Ewald method,35 and the temperature was controlled using the Nosé−Hoover thermostat.36,37 Most simulations employed the six-site model38 (normal melting point, 289 ± 3 K39), but the TIP4P/2005 model40 (normal melting point, 250.5 ± 3 K41) is also considered for comparison with previous work.3 Both are rigid water models, and in the simulations bond lengths and angles are constrained using the LINCS algorithm.42 As in our previous simulations,9 silver iodide is modeled as a rigid lattice (possible surface reconstruction43 is not considered) employing the force field developed by Hale and Keifer,20 neglecting the polarization terms. Calculations and a discussion of the polarization contribution to the silver−water and iodide−water interactions are given in the Supporting Information (see Section A and Figures S8 and S9). In our previous work9 using infinite sheets of AgI, we tested the sensitivity of ice nucleation to different interaction potentials by varying the partial charges on the silver and iodide

Figure 1. Top view of the AgI disks and plates considered. Silver and green spheres represent silver and iodide ions, respectively. 2292

DOI: 10.1021/acs.jpcb.5b06605 J. Phys. Chem. B 2016, 120, 2291−2299

Article

The Journal of Physical Chemistry B

number of water molecules included for the different systems are given in Table 1. As previously noted, most simulations were carried out under NVT conditions with water densities lying in the range 0.94 to 0.96 g/mL (Table S1); however, to check that our results are not strongly influenced by the relatively low densities (closer to ice than to liquid water) or the constant volume condition, simulations for systems A, B, C, and D (Table S2) were carried out under NPT conditions with the pressure controlled at 1 bar using the Parrinello−Rahman barostat. In these NPT simulations, ice nucleated on the surfaces via the same mechanism and in simulation times comparable to the NVT results (see Figure S7). Thus, at least for the conditions we consider, the simulation results are not sensitive to the ensemble employed. The configuration labeled C2 contained two “mirrored” plates located 10 nm apart in the simulation cell. This system was used to check that spurious electric fields were not influencing our single disk or plate simulations, as can be the case for slabs of charged particles that are periodically infinite in two dimensions, but of finite thickness.44 (Also see Figure S6 and the discussion in the Supporting Information.) As previously noted, we also carried out simulations (for systems C and D) with no charges on the sliver and iodide ions. Ice nucleation occurred in these simulations, just as it does when the full force field (with charges of ±0.6e) is used (see Figure S7). This provides a further demonstration that surface fields are not strongly influencing our results. All simulations investigating nucleation were begun with positional coordinates obtained at 0.5 ns time intervals from an equilibrium simulation trajectory at 300 K. Production runs were initiated at 300 K with velocities drawn randomly from a Maxwell−Boltzmann distribution and were relaxed at 300 K for 1 ns. Following this, the systems were cooled to the desired temperature over a period of 1 ns and held at that temperature for the duration of the run. For the six-site model, at least two runs were performed at the highest temperature where nucleation was observed. Two additional runs, at a temperature 1 K warmer than the ice nucleation temperature, were carried out for times comparable to those where ice nucleation occurred. This was done to increase confidence that we had achieved the highest possible temperature where a particular AgI disk or plate could nucleate ice, at least on simulation time scales. Complete lists of the simulations performed are given in Tables S2 and S3. Simulations were analyzed for the presence of ice using the CHILL algorithm45 with a slightly relaxed definition to allow for the reduced tetrahedral nature of the water models used here compared with the mW model46 used in the original work. The output from the CHILL algorithm was then analyzed with the DBScan algorithm47 to find the largest cluster. DBScan identifies clusters based on groups of ice molecules that have a minimum number of ice molecules as nearest neighbors. For DBScan we use a near-neighbor cutoff of 3.8 Å and a minimum of four ice-molecule neighbors. Unless otherwise stated, all temperatures are given as the degree of supercooling, defined as ΔTs = Tm − T, where Tm is the normal melting point of bulk ice for a particular water model and T is the temperature of the simulation.

Figure 2. Series of snapshots of an ice cluster nucleating bulk ice on AgI disk D. At each time, side and top views of the ice cluster are shown. Silver and iodide ions are shown in silver and green, respectively. In the first and third columns, black, blue, and orange represent hydrogen atoms, the oxygen atoms of Ih, and the oxygen atoms of Ic, respectively. In the second and fourth columns, yellow and purple represent the oxygen atoms of ice inside and outside the cylindrical volume defined by the disk radius, respectively. Liquid water molecules are not shown in the Figure.

neutral containing equal numbers of Ag and I ions. Disks and plates labeled A−E were obtained by retaining all atoms within a specified circumscribing radius, which naturally leads to hexagonal shapes for the smaller plates (A−C) and to roughly circular shapes for the larger disks (D,E), as can be seen in Figure 1. The rectangular plate (labeled Dr) has the same number of atoms as disk D and is considered to probe for possible effects of shape on the nucleation properties. The circumscribing radii of the disks and plates together with the rectangular dimensions of plate Dr are given in Table 1. Table 1. Summary of the Systems Investigated with the SixSite Model as Well as Key Results system A B C C2 D Dr E

disk radius (nm)

number of waters

X/Y/Z dimensions of cell (nm)

min. observed nucleation ΔTs (K)

1.15 1.61 2.07 2.07 (2 plates) 2.53 4.507 × 4.229 (rectangle) 2.99

5350 7620 7450 17650 11870 11870

5/5/7 6/6/7 6/6/7 6/6/16.44 7/7/8.222 7/7/8.22

35 26 22 21 16 16

17150

8/8/9

14

In all simulations, with one exception (discussed later), the central simulation cell contained a single AgI disk or plate, centered in the XY plane with the (0001) face perpendicular to the Z axis. Periodic boundary conditions are applied in all directions so the actual location of the disk or plate along the Z axis has no influence on the physical behavior of the system; however, for purposes of analysis and depiction the disk or plate is located near the bottom of the cell with the silver exposed surface on which ice nucleates directed upward in the positive Z direction (see Figure 2). The cell dimensions and the

III. RESULTS AND DISCUSSION In our simulations ice clusters form only on the silver exposed face of the AgI disks and plates, consistent with previous work 2293

DOI: 10.1021/acs.jpcb.5b06605 J. Phys. Chem. B 2016, 120, 2291−2299

Article

The Journal of Physical Chemistry B showing that this face nucleates ice very effectively.9,10 If the surface-induced ice cluster achieves a critical size (dependent on ΔTs), then ice nucleates and grows throughout the simulation cell. Several configurational snapshots of an ice cluster on the silver exposed surface of disk D nucleating bulk ice at ΔTs = 16 K are shown in Figure 2. We note that at previous times (30− 91 ns in Figure 2) the ice cluster is contained within the cylindrical boundary of the disk and is roughly hemispherical in shape. At 91 ns the cluster has achieved critical size and begins to extend beyond the cylindrical disk boundary (see later for further discussion on this point), and it then grows primarily horizontally eventually filling the simulation cell. As identified by the CHILL algorithm, the ice on this disk (and the other disk and plates considered) is a mixture of Ih and Ic, or Isd.33 The c axis of Ih is perpendicular to the disk surface. Thus, the initial ice growth is on the basal face as the cluster grows taller. Later, horizontal growth occurs on what would be prism faces if the ice was completely Ih. A video of this simulation is available in the Supporting Information. For each system considered, the smallest ΔTs for which bulk ice nucleation was observed is given in Table 1. For the disks and plates (Figure 1A−D) ΔTs is plotted versus the disk or plate radius in Figure 3 and compared with previous results. The black horizontal line in Figure 3 is the minimum observed ΔTs (8 K) required to nucleate ice on infinite silver iodide

slabs.9 The trend in our results indicates that the infinite surface limit is being approached as the disk or plate radius increases. We note (Table 1) that the rectangular surface (Dr) nucleated ice at the same ΔTs as the disk D, suggesting that the particular shape of the exposed AgI surface is not a major factor influencing ice nucleation. Also, the mirrored plate system C2 nucleated ice only one degree warmer than the single disk or plate system, indicating that the small electric fields and resulting polarization of the water produced by single disks or plates (see Figure S6) do not significantly influence the nucleation temperature. It is possible that the one degree difference found for the C2 case is due to the ice-like clusters on both plates interacting, making it easier to achieve critical cluster size. Additionally, note that the C and D surfaces nucleate ice at the same temperatures when the silver and iodide ions are completely uncharged, and hence no electric field is present (Figure S7). The red points and red horizontal line in Figure 3 are the results of Lupi et al.8 for ice nucleation on graphene obtained using the mW water model.46 In that system, ice nucleation also begins with ice-like clusters forming on the disks, and the trend observed is similar to our results; however, we note that for the mW model the infinite surface limit is achieved at a radius of ∼2 nm, whereas in our case the approach to the infinite limit appears to be asymptotic. In Figure 3 we also compare our results with predictions for homogeneous nucleation of an infinitely long cylindrical nucleus,5,6 homogeneous nucleation of a spherical critical nucleus,3,4 and predictions based on the Gibbs−Thomson (GT) equation, which relates ΔTs to the size of a critical ice cluster or nucleus. A second-order expansion of the GT equation can be written in the form48 ⎛ κT ⎛T ⎞ ΔcTm ⎞ Tm 1 ⎟ (γK )2 ΔTs = ⎜ m ⎟γK + ⎜ m − χ − 2 2ρq2 ⎠ ρq ⎝ ρq ⎠ ⎝ ρq

(1)

where (assuming homogeneous nucleation) q is the heat of melting, ρ is the density of ice, γ is the interfacial tension between ice and water, χ is the compressibility of ice, κ is the thermal expansivity, and Δc is the difference between the isobaric specific heats of water and ice. In eq 1, K depends on the geometry of the ice cluster and is 2/r for a sphere of radius r and 1/r for a cylinder, where r is the radius of the base. In both cases, eq 1 reduces to the form ΔTs = Figure 3. ΔTs required to nucleate ice in 20−200 ns versus the radius of the nucleating disk or plate compared with previous simulation results and with the GT equation based on the properties of real water. Heterogeneous nucleation results for AgI disks or plates are shown in black (see legend) for the six-site and TIP4P/2005 models. The solid black line indicates the lowest supercooling required to nucleate ice on infinite AgI slabs, and the dotted black curve represents a fit of our sixsite model results to the GT equation. Red downward triangles and the red line show heterogeneous nucleation results8 for the mW model on graphene disks and infinite slabs, respectively. Blue triangles are homogeneous nucleation results for the TIP4P-E model obtained employing infinite cylindrical ice clusters.6 The GT equation for cylinders48 is also shown (dashed blue curve). Homogeneous nucleation results employing spherical ice clusters are in green (see legend) for the TIP4P/2005,3 TIP4P/Ice,3 and TIP4P4 models. The real-water GT equation48 for spherical nuclei is shown as a green dashed curve.

C1 C − 22 r r

(2)

where C1 and C2 are constants that depend on the geometry. Not all of the parameters determining these constants are known for the water models considered here, but the experimental values for water yield48 C1 = 519 K Å, C2 = 829 K Å2 for spherical clusters and C1 = 259 K Å, C2 = 207 K Å2 for the cylindrical case. Figure 3 suggests that the radius of the disk or plate required for heterogeneous nucleation on AgI at a given supercooling agrees well with the critical radius for homogeneous nucleation of spheres at the same supercooling. A similar trend has recently been observed for heterogeneous nucleation on watersoluble macromolecules.49 In contrast, the radius of the disk or plate required for heterogeneous nucleation on AgI is larger than the critical radius for homogeneous nucleation of infinite cylinders at the same supercooling. 2294

DOI: 10.1021/acs.jpcb.5b06605 J. Phys. Chem. B 2016, 120, 2291−2299

Article

The Journal of Physical Chemistry B

Figure 4. Plots showing the number of water molecules identified as ice-like growing on each disk or plate at the minimum ΔTs, where ice nucleation was observed, and at a temperature one degree warmer. The ΔTs values are given in the legends. Our estimates of the minimum and maximum limits on the critical cluster size (see text) are shown as horizontal lines and open circles, respectively.

Three GT curves are given in Figure 3. The curves for spherical and cylindrical clusters were obtained using the experimental parameters for water. We note that the results of Pereyra et al.6 lie close to the GT equation for cylindrical nuclei. In Figure 3, we also include the curve for the parameters C1 = 450 K Å, and C2 = 537 K Å2, which represent the best leastsquares fit of eq 2 to our simulation results. We see that eq 2 does give a very good fit to our results. Thus, the twoparameter GT equation strictly applicable to homogeneous nucleation gives a reasonably good representation of heterogeneous nucleation by AgI disks or plates, likely indicating that the ice-liquid interfacial tension remains a dominant factor in determining critical cluster size. Figure 4 is a set of plots showing the number of water molecules identified as ice-like forming clusters on various AgI disks and plates as time progresses. Additional plots for the C2 case and for simulations with the TIP4P/2005 model are shown in Figures S1 and S2, respectively. As previously described, the CHILL45 algorithm was used to identify which water molecules are ice-like, and the DBScan47 algorithm was employed to determine which ice-like molecules are part of the cluster on the disk or plate surface. Simulation results for two ΔTs values are shown in Figure 4; one corresponds to the smallest ΔTs for which ice nucleation was observed, and the other corresponds to simulations one degree warmer, where ice did not nucleate. We note that in simulations where ice nucleation occurred, it did so in 100−200 ns after supercooling. Of course, it is always possible that ice nucleation would occur at higher temperatures (smaller ΔTs) with longer simulation times, but our simulations represent what is presently practical. For ice to nucleate and grow in our simulations, a critical ice cluster size (critical nucleus) must be achieved. We estimate upper and lower limits to the number of water molecules in the critical cluster. In Figure 4, lower limits are shown as black horizontal lines. These estimates were obtained from

simulations that did not nucleate ice and that were one degree warmer than simulations that did. The number of water molecules in the ice cluster was recorded every nanosecond (Figure 4), and the average number of water molecules, together with the associated standard deviation, was calculated using the results for times >20 ns. The lower limit on cluster size was taken to be two standard deviations above the average size obtained in these simulations. Upper limits on the number of water molecules in the critical cluster were estimated for each simulation that nucleated ice and were taken to be the size at which the cluster begins to extend beyond a cylindrical volume defined by the radius of the disk or plate, or, more precisely, when the number of ice-like molecules outside the cylinder exceeds 10. Our justification for this somewhat arbitrary selection of an upper limit is based on Figure 5, where it can be seen that lateral cluster growth beyond a disk-defined cylinder closely coincides with nucleation and irreversible growth of bulk ice. For comparison, the disk D results shown in Figure 5 are from the same simulation depicted in Figure 2. In Figure 5, the number of water molecules in the ice cluster growing on the disk or plate is plotted as a function of time. The ice cluster is further resolved by counting the number of molecules inside and outside the cylinder, with a base defined by the radius of the circle circumscribing the disk or plate. Similar plots for other ice-nucleating simulations employing the six-site model are given in Figure S3, and results for the TIP4P/ 2005 model are shown in Figure S4. The plots shown in Figures 4 and 5 and Figures S3 and S4 suggest that ice nucleation and growth on AgI disks and plates occur in two stages. (Note that the plots obtained in NPT simulations (Figure S7) show very similar behavior.) There is an initial period, varying in length depending on the disk or plate size, where the ice cluster slowly grows, sometimes fluctuating in size, within the confines of the cylindrical space determined by the disk or plate radius, as previously described. This is followed by rapid growth once a critical cluster size is 2295

DOI: 10.1021/acs.jpcb.5b06605 J. Phys. Chem. B 2016, 120, 2291−2299

Article

The Journal of Physical Chemistry B

Figure 5. Number of ice-like molecules as a function of time for four disks and plates. Curves are included for the total number of ice molecules (black) and for the numbers inside (red) and outside (green) the cylindrical volume defined by the disk or plate radius. In each panel, a circle denotes the point identified as the upper limit on the critical cluster size.

reached. Much of the rapid growth occurs horizontally (see the green line in Figure 5) on what would be the prism face of pure Ih, where growth is known to be faster than on the basal face50 (perpendicular to the surface of the AgI disk or plate). There is an apparent correlation between when the cluster begins to extend beyond the disk-based cylinder and the onset of rapid growth. Therefore, we define (admittedly somewhat arbitrarily) the transition from the first stage of growth to the second as when the ice cluster has at least ten molecules beyond the edge of the disk or plate, and these critical clusters are marked with circles in Figures 4 and 5. The number of water molecules in the critical cluster are plotted as vertical bars joining the lower and upper limit estimates and are shown as functions of ΔTs in Figure 6. The CHILL algorithm is not efficient at recognizing ice-like molecules that are immediately adjacent to a surface; therefore. we add additional molecules to the cluster numbers reported by CHILL and DBScan. From snapshots of water on the surfaces, we estimate the corrections to be 21, 50, 75, 114, 131, and 173 for disk or plate A, B, C, D, Dr, and E, respectively. In Figure 6A, we also included the number of water molecules in the critical cluster predicted by the GT equation for homogeneous nucleation of a sphere using the parameters of real water. GT estimates for the number of ice molecules in spherical clusters are obtained assuming an ice density of 917 kg/m3. We note that our results do follow the trend predicted by the GT equation, but our estimates lie consistently below the GT curve. The somewhat smaller critical clusters we observe possibly reflect the influence of heterogeneous nucleation by AgI. In panel A we also compare our estimates with previous simulations of the critical cluster size for homogeneous nucleation. The results of Sanz et al.3 employed the TIP4P/ 2005 40 and TIP4P/ice 51 models and come from an investigation of homogeneous ice nucleation using spherical Ih clusters as trial nuclei. Espinosa et al.4 performed similar studies for the TIP4P52 and mW46 water models, but we only show the results for TIP4P because the results for mW are

Figure 6. Comparison of critical cluster results from our simulations with previous experimental and theoretical work. (A) Comparison with previous simulations. Present results are shown as bars joining the upper and lower estimates for the critical cluster and are labeled indicating the disk or plate employed. Results for the six-site and TIP4P/2005 models are shown in black and pink, respectively. The green dotted curve represents the real-water GT equation.48 Green points (see legend) are homogeneous nucleation results for the TIP4P/2005,3 TIP4P/Ice,3 and TIP4P4 models. Blue triangles are homogeneous nucleation results for the six-site model obtained in an applied electric field.53 (B) Comparison with experimental estimates obtained using microemulsions54 for homogeneous nucleation (purple squares, see legend) and heterogeneous nucleation (orange downward triangles). The open purple squares are homogeneous nucleation results obtained employing levitated droplets.55 Other symbols are as in panel A. (C) Simulation (present work) and experimental results plotted versus the scaled temperature τ. The symbols and curves are as in panels A and B.

similar to those shown for TIP4P/ice. Yan et al.53 investigated ice nucleation by electric fields using the six-site model of water. Electric fields raise the melting temperature of water and thus allow homogeneous nucleation to be observed in direct simulations employing atomistic models but do not have a large effect on the size of the critical nucleus for a given ΔTs.53 The results of Yan et al.53 (six-site model) are in good accord with the present estimates, which is perhaps expected given that the same water model is employed in both cases. The results of Sanz et al.3 (TIP4P/20005 and TIP4P/ice) and Espinosa et al.4 (TIP4P) are higher than our estimates for heterogeneous nucleation on AgI surfaces. This difference may partially be due to the different models considered or to the fact that Sanz et 2296

DOI: 10.1021/acs.jpcb.5b06605 J. Phys. Chem. B 2016, 120, 2291−2299

Article

The Journal of Physical Chemistry B al.3 and Espinosa et al.4 used spherical ice clusters that will present a variety of crystal faces for nucleation. Different faces will grow differently and may impede nucleation, thus requiring larger sizes. In general, panel A suggests that the number of water molecules in the critical cluster for heterogeneous nucleation on AgI is within a factor of 4 of the number of water molecules in the critical cluster required for homogeneous nucleation of spheres at the same supercooling. To test for possible water model sensitivity, we performed simulations with disks D and E using the TIP4P/2005 model considered by Sanz et al.3 We find that the TIP4P/2005 model does require a larger ΔTs (compared with the six-site model) to nucleate on disks D and E, while giving critical clusters similar in size to those found for the six-site model (see Figure 6A,C). Thus, the particular water models employed to study ice nucleation are also a factor in determining the critical cluster size. Some experimental results for the critical cluster of ice obtained under different conditions are compared with our estimates in Figure 6B. Liu et al.54 investigated both homogeneous and heterogeneous nucleation using water in oil microemulsions. To investigate heterogeneous nucleation they incorporated heptacosanol as the ice nucleating agent into the emulsions. In addition, Kramer et al.55 investigated homogeneous nucleation using levitated single water droplets. Our results show reasonable agreement with the homogeneous studies, which is consistent with the conclusions reached from Figure 6A. For the case of heterogeneous nucleation on heptacosanol, our results deviate from the measurements at warmer temperatures, which is discussed later. Malaspina et al.56 have suggested that another way of comparing different water models is to plot them versus a scaled temperature, τ, as is done in Figure 6C. The scaled temperature is defined as τ = (T − Tw)/(Tm − Tw), where T is the nucleation temperature, Tm is the melting temperature of ice for a particular model, and Tw is the temperature of the Widom line, defined by heat capacity maxima. The Widom line temperatures at 1 atm were taken from Malaspina et al.56 to be 249 K for six-site, 229 K for TIP4P/2005, and 235 K for real water. Widom line temperatures are not known for the TIP4P and TIP4P/Ice models, and these models are not included in panel C. The melting temperatures used are 289 K for six-site39 and 250.5 K for TIP4P/2005.41 The scaled temperature is an attempt to account for dynamical differences between real water and water models and also between different water models that often have significantly different melting points56 and dynamical properties. We note that on the τ scale the TIP4P/2005 results (both those of Sanz et al.3 and those we obtain) still differ significantly from the estimates for the six-site model. Thus, model differences that are not captured by τ clearly persist and should be kept in mind in ice nucleation studies of water models. From Figure 6B,C, we see that our model estimates and the experimental results for heterogeneous nucleation on heptacosanol are in reasonable agreement at lower temperatures but deviate at the higher temperatures. At higher temperatures, the experiments report smaller amounts of water in critical clusters and appear to follow a different trend. The higher temperature results represent the onset of heterogeneous nucleation using the IN heptacosanol. Our results are also for heterogeneous nucleation but appear to more closely follow the trend of homogeneous nucleation.

The apparent similarity of heterogeneous nucleation by AgI disks to the homogeneous case in terms of the critical cluster size may reflect the fact that AgI is a close ice mimic.9 In other words, from the point of view of ice nucleation and growth, there may be a strong resemblance between a AgI disk or plate and a corresponding piece of ice, the crucial difference being that the AgI disk or plate obviously cannot melt. Some additional insight can be gained by examining the sizes of subcritical clusters that never achieve nucleation. We collect samples of these “failed” clusters using the local maxima in Figure 4 and in similar plots (not shown) at other temperatures. An example plot illustrating how we identified failed clusters for plate C at ΔTs = 19 K is given in the Supporting Information (Figure S5). Note that only the larger failed clusters are selected. The number of water molecules in failed clusters for different disks or plates (denoted with different colors) is plotted versus ΔTs in Figure 7. The numbers of water molecules in critical

Figure 7. Sizes of subcritical clusters (six-site model) that failed to nucleate for different disks and plates as functions of ΔTs. Critical cluster size ranges for the smallest ΔTs where nucleation occurred are shown for each disk or plate as bars joining estimates of upper and lower limits.

nuclei are included for comparison. Two trends are immediately obvious from Figure 7. First, as one would expect, larger disks and plates create larger clusters. Second, for a given disk or plate, the number of water molecules in the largest failed cluster decreases with increasing temperature. This second observation indicates an additional hurdle to ice nucleation at warmer temperatures. Not only are larger clusters required to nucleate ice at warmer temperatures but also it becomes harder to create larger clusters as the temperature is increased. Therefore, in our simulations the warmest temperature where a particular AgI disk or plate can nucleate bulk ice is the temperature where the required number of water molecules in the critical cluster just matches the largest cluster that can form on the disk or plate.

IV. CONCLUSIONS We have used molecular dynamics simulations to investigate ice nucleation by AgI disks or plates with circumscribing radii ranging from 1.15 to 2.99 nm. The six-site and TIP4P/2005 water models were considered, and a similar ice nucleation process was observed in both cases. The disk or plate size required to nucleate ice on simulation time scales increases as the degree of supercooling ΔTs decreases. For the six-site model the largest disk considered (radius of 2.99 nm) nucleated 2297

DOI: 10.1021/acs.jpcb.5b06605 J. Phys. Chem. B 2016, 120, 2291−2299

Article

The Journal of Physical Chemistry B ice at ΔTs = 14 K, compared with 8 K for an infinite AgI sheet.9 Our results demonstrate that relatively small patches of defectfree AgI surface could serve as efficient “active sites” for ice nucleation. Ice nucleation by nanoscale AgI disks and plates occurs in two stages. Initially, upon supercooling, ice clusters roughly hemispherical in shape form on the silver exposed face of the AgI disk or plate. These clusters fluctuate in size but never melt away until (if the disk size and ΔTs are suitably matched) the cluster achieves critical size, allowing bulk ice to nucleate and rapidly grow to fill the simulation cell. For a given ΔTs, our estimates of the critical cluster size are roughly comparable to those obtained from previous simulation studies of homogeneous nucleation. Additionally, the critical cluster size as a function of ΔTs is well represented by a second-order (twoparameter) Gibbs−Thomson equation, as one would expect for homogeneous nucleation. Thus, it appears that the main influence of the AgI disk or plate is to create and stabilize an ice cluster of some minimum size (dependent on ΔTs and the radius of the AgI disks or plate) that is always present in the system. This constantly present ice cluster reduces the freeenergy barrier that must be overcome to achieve critical cluster size and hence ice nucleation.

■

(3) Sanz, E.; Vega, C.; Espinosa, J. R.; Caballero-Bernal, R.; Abascal, J. L. F.; Valeriani, C. Homogeneous Ice Nucleation at Moderate Supercooling from Molecular Simulation. J. Am. Chem. Soc. 2013, 135, 15008−15017. (4) Espinosa, J. R.; Sanz, E.; Valeriani, C.; Vega, C. Homogeneous Ice Nucleation Evaluated for Several Water Models. J. Chem. Phys. 2014, 141, 18C529-1−18C529-14. (5) Pereyra, R. G.; Carignano, M. A. Ice Nanocolumns: A Molecular Dynamics Study. J. Phys. Chem. C 2009, 113, 12699−12705. (6) Pereyra, R. G.; Szleifer, I.; Carignano, M. A. Temperature Dependence of Ice Critical Nucleus Size. J. Chem. Phys. 2011, 135, 034508-1−034508-5. (7) Moore, E. B.; Molinero, V. Structural Transformation in Supercooled Water Controls the Crystallization Rate of Ice. Nature 2011, 479, 506−509. (8) Lupi, L.; Hudait, A.; Molinero, V. Heterogeneous Nucleation of Ice on Carbon Surfaces. J. Am. Chem. Soc. 2014, 136, 3156−3164. (9) Zielke, S. A.; Bertram, A. K.; Patey, G. N. A Molecular Mechanism of Ice Nucleation on Model AgI Surfaces. J. Phys. Chem. B 2015, 119, 9049−9055. (10) Fraux, G.; Doye, J. P. K. Note: Heterogeneous Ice Nucleation on Silver-Iodide-Like Surfaces. J. Chem. Phys. 2014, 141, 216101-1− 216101-2. (11) Cox, S. J.; Raza, Z.; Kathmann, S. M.; Slater, B.; Michaelides, A. The Microscopic Features of Heterogenous Ice Nucleation May Affect the Macroscopic Morphology of Atmospheric Ice Crystals. Faraday Discuss. 2014, 167, 389−403. (12) Cox, S. J.; Kathmann, S. M.; Slater, B.; Michaelides, A. Molecular Simulations of Heterogeneous Ice Nucleation. I. Controlling Ice Nucleation Through Surface Hydrophilicity. J. Chem. Phys. 2015, 142, 184704-1−184704-5. (13) Cox, S. J.; Kathmann, S. M.; Slater, B.; Michaelides, A. Molecular Simulations of Heterogeneous Ice Nucleation. II. Peeling Back the Layers. J. Chem. Phys. 2015, 142, 184705-1−184705-8. (14) Vonnegut, B. The Nucleation of Ice Formation by Silver Iodide. J. Appl. Phys. 1947, 18, 593−595. (15) Bruintjes, R. T. A Review of Cloud Seeding Experiments to Enhance Precipitation and Some New Prospects. Bull. Am. Meteorol. Soc. 1999, 80, 805−820. (16) Breed, D.; Rasmussen, R.; Weeks, C.; Boe, B.; Deshler, T. Evaluating Winter Orogrophic Cloud Seeding: Design of the Wyoming Weather Modification Pilot Project (WWMPP). J. Appl. Meteorol. Clim. 2014, 53, 282−299. (17) Burley, G. Polymorphism of Silver Iodide. Am. Mineral. 1963, 48, 1266−1276. (18) Burley, G. Structure of Hexagonal Silver Iodide. J. Chem. Phys. 1963, 38, 2807−2812. (19) Fukata, N.; Paik, Y. Water Adsorption and Ice Nucleation on Silver Iodide Surfaces. J. Appl. Phys. 1973, 44, 1092−1100. (20) Hale, B. N.; Kiefer, J. Studies of H2O on β-AgI Surfaces: An Effective Pair Potential Model. J. Chem. Phys. 1980, 73, 923−933. (21) Hale, B. N.; Kiefer, J.; Terrazas, S.; Ward, R. C. Theoretical Studies of Water Adsorbed on Silver Iodide. J. Phys. Chem. 1980, 84, 1473−1479. (22) Shevkunov, S. V. Adsorption of Water Vapor on the AgI Surface: A Computer Experiment. Russ. J. Gen. Chem. 2005, 75, 1632− 1642. (23) Shevkunov, S. V. Layer-by-Layer Adsorption of Water Molecules on the Surface of a Silver Iodide Crystal. Russ. J. Phys. Ch. 2006, 80, 769−775. (24) Shevkunov, S. V. Clusterization of Water Molecules on Crystalline β-AgI Surface. Computer Experiment. Colloid J. 2006, 68, 632−643. (25) Taylor, J. H.; Hale, B. N. Monte Carlo Simulations of Water-Ice Layers on a Model Silver Iodide Substrate: A Comparison with Bulk Ice Systems. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 47, 9732−9741.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b06605. Table of bulk water densities; tables of simulations; plots of ice cluster size versus time for different systems; additional cylinder plots for both water models; an example of failed cluster selection; polarization plots; plots of ice cluster size versus time for models with zero charge and for NPT simulations; and plots comparing force fields. (PDF) Video of simulation of ice nucleation by AgI disks and plates. (MPG)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The financial support of the Natural Science and Engineering Research Council of Canada is gratefully acknowledged. This research has been enabled by the use of WestGrid and Compute/Calcul Canada computing resources, which are funded in part by the Canada Foundation for Innovation, Alberta Innovation and Science, BC Advanced Education, and the participating research institutions. WestGrid and Compute/ Calcul Canada equipment is provided by IBM, Hewlett-Packard and SGI.

■

REFERENCES

(1) Seinfeld, J. H.; Pandis, S. N. Atmospheric Chemistry and Physics: From Air Pollution to Climate Change, 2nd ed; Wiley-Interscience: Hoboken, NJ, 2006. (2) Pruppacher, H. R.; Klett, J. D. Microphysics of Clounds and Precipitation; Springer: New York, 2010. 2298

DOI: 10.1021/acs.jpcb.5b06605 J. Phys. Chem. B 2016, 120, 2291−2299

Article

The Journal of Physical Chemistry B (26) Volkov, V. B.; Zhogolev, D. A.; Bakhanova, R. A.; Kiselev, V. J. A Quantum-Chemical Study of the Interaction of the Water Molecule with Silver Iodide Surface Clusters. Chem. Phys. Lett. 1976, 43, 45−48. (27) Ward, R. C.; Holdman, J. M.; Hale, B. N. Monte Carlo Studies of Water Monolayer Clusters on Substrates: Hexagonal AgI. J. Chem. Phys. 1982, 77, 3198−3202. (28) Ward, R. C.; Hale, B. N.; Terrazas, S. A Study of the Critical Cluster Size for Water Monolayer Clusters on a Model AgI Basal Substrate. J. Chem. Phys. 1983, 78, 420−423. (29) Malenkov, G. Liquid Water and Ices: Understanding the Structure and Physical Properties. J. Phys.: Condens. Matter 2009, 21, 283101-1−283101-35. (30) Carignano, M. A. Formation of Stacking Faults During Ice Growth on Hexagonal and Cubic Substrates. J. Phys. Chem. C 2007, 111, 501−504. (31) Malkin, T. L.; Murray, B. J.; Brukhno, A. V.; Anwar, J.; Salzmann, C. G. Structure of Ice Crystallized from Supercooled Water. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 1041−1045. (32) Kuhs, W. F.; Sippel, C.; Falenty, A.; Hansen, T. C. Extent and Relevance of Stacking Disorder in “Ice Ic. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 21259−21264. (33) Malkin, T. L.; Murray, B. J.; Salzmann, C. G.; Molinero, V.; Pickering, S. J.; Whale, T. F. Stacking Disorder in Ice I. Phys. Chem. Chem. Phys. 2015, 17, 60−76. (34) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D. GROMACS 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845− 854. (35) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (36) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519. (37) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (38) Nada, H.; van der Eerden, J. P. J. M. An Intermolecular Potential Model for the Simulation of Ice and Water Near the Melting Point: A Six-Site Model of H2O. J. Chem. Phys. 2003, 118, 7401−7413. (39) Abascal, J. L. F.; Fernández, R. G.; Vega, C.; Carignano, M. A. The Melting Temperature of the Six Site Potential Model of Water. J. Chem. Phys. 2006, 125, 166101-1−166101-2. (40) Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505-1−234505-12. (41) Fernández, R. G.; Abascal, J. L. F.; Vega, C. The Melting Point of Ice Ih for Common Water Models Calulated from Direct Coexistance of the Solid-Liquid Interface. J. Chem. Phys. 2006, 124, 144506-1−144506-11. (42) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Slover for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463−1472. (43) Tasker, P. W. The Stability of Ionic Crystal Surfaces. J. Phys. C: Solid State Phys. 1979, 12, 4977−4984. (44) Croteau, T.; Bertram, A. K.; Patey, G. N. Simulation of Water Adsorption on Kaolinite Under Atmospheric Conditions. J. Phys. Chem. A 2009, 113, 7826−7833. (45) Moore, E. B.; de la Llave, E.; Welke, K.; Scherlis, D. A.; Molinero, V. Freezing, Melting and Structure of Ice in a Hydrophilic Nanopore. Phys. Chem. Chem. Phys. 2010, 12, 4124−4134. (46) Molinero, V.; Moore, E. B. Water Modeled As an Intermediate Element Between Carbon and Silicon. J. Phys. Chem. B 2009, 113, 4008−4016. (47) Ester, M.; Kriegel, H.-P.; Sander, J.; Xu, X. A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise. Proc. 2nd Int. Conf. Knowl. Discovery Data Mining 1996, 226− 231.

(48) Mori, A.; Maruyama, M.; Furukawa, Y. Second-Order Expansion of Gibbs-Thomson Equation and Melting Point Depression of Ice Crystallite. J. Phys. Soc. Jpn. 1996, 65, 2742−2744. (49) Pummer, B. G.; Budke, C.; Augustin-Bauditz, S.; Niedermeier, D.; Felgitsch, L.; Kampf, C. J.; Huber, R. G.; Liedl, K.; Loerting, T.; Moschen, T.; et al. Ice Nucleation by Water-Soluble Macromolecules. Atmos. Chem. Phys. 2015, 15, 4077−4091. (50) Seo, M.; Jang, E.; Kim, K.; Choi, S.; Kim, J. S. Understanding Anisotropic Growth Behavior of Hexagonal Ice on a Molecular Scale: A Molecular Dynamics Simulation Study. J. Chem. Phys. 2012, 137, 154503-1−154503-8. (51) Abascal, J. L. F.; Sanz, E.; Garciá Fernández, R. G.; Vega, C. A Potential Model for the Study of Ices and Amorphous Water: TIP4P/ Ice. J. Chem. Phys. 2005, 122, 234511-1−234511-9. (52) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (53) Yan, J. Y.; Overduin, S. D.; Patey, G. N. Understanding Electrofreezing in Water Simulations. J. Chem. Phys. 2014, 141, 074501-1−074501-9. (54) Liu, J.; Nicholson, C. E.; Cooper, S. J. Direct Measurment of Critical Nucleus Size in Confined Volumes. Langmuir 2007, 23, 7286− 7292. (55) Krämer, B.; Hübner, O.; Vortisch, H.; Wöste, L.; Leisner, T.; Schwell, M.; Rühl, E.; Baumgärtel, H. Homogeneous Nucleation Rates of Supercooled Water Measured in Single Levitated Microdroplets. J. Chem. Phys. 1999, 111, 6521−6527. (56) Malaspina, D. C.; di Lorenzo, A. J. B.; Pereyra, R. G.; Szleifer, I.; Carignano, M. A. The Water Supercooled Regime as Described by Four Common Water Models. J. Chem. Phys. 2013, 139, 024506-1− 024506-8.

2299

DOI: 10.1021/acs.jpcb.5b06605 J. Phys. Chem. B 2016, 120, 2291−2299