501
2007, 111, 501-504 Published on Web 12/19/2006
Formation of Stacking Faults during Ice Growth on Hexagonal and Cubic Substrates Marcelo A. Carignano* Department of Chemistry, Purdue UniVersity, 560 OVal DriVe, West Lafayette, Indiana 47907-1393 ReceiVed: NoVember 8, 2006; In Final Form: NoVember 30, 2006
We present a molecular dynamics simulation study of ice growth from supercooled water in contact with: (i) the basal plane of hexagonal ice and (ii) the (111) plane of cubic ice. We observe that, regardless of the substrate, a variable number of stacking faults appear during the crystallization. The simulations also show that the formation of cubic ice is more likely than the formation of hexagonal ice, in agreement with recent experimental findings and the Ostwald step rule.
Hexagonal ice (ice Ih) is the stable phase of ice upon freezing at ambient pressure or below.1 Cubic ice (ice Ic) instead, is a metastable phase that transforms irreversibly into ice Ih.2 Recent experiments showed that ice Ic is the dominant crystalline phase that nucleates when pure water droplets freeze homogeneously at temperatures as high as 235 K.3,4 The fraction of ice Ic in the frozen droplets depends on the temperature distance to the melting point. A similar result has been observed on water freezing in mesopores.5 In both situations, namely in water droplets and confined water, the X-ray diffraction patterns obtained after freezing correspond to a mixture of ice Ic and Ih. The picture that emerges from the experimental evidence is that a mixture of ice Ic and Ih forms naturally, suggesting that ice Ic may grow over an initial nucleus with hexagonal structure, and ice Ih may grow over nucleus with cubic structure. This scenario is investigated in this work by molecular dynamics simulations. The interest in molecular dynamics simulations of ice/water systems has been in augment in the recent years, fueled by the ever increasing speed of computers that now allow direct observation of nucleation and ice growth. In spite of this advance, computer resources are still the limiting factor at the time of selecting the water model, system size and running parameters. The great majority of the existing water models were developed to reproduce known properties of the liquid phase at ambient conditions. Yet, these models were used to investigate crystallization properties. For example, the TIP4P model was used by Matsumoto et al.6 to show the nucleation of ice after 290 ns of simulation in a system with 512 water molecules. The TIP4P model was also used by Nada and Furukawa7 to investigate the growth mechanisms on the different crystal faces of ice. Vrbka and Jungwirth8 used the SPC/E model to investigate ice growth from pure water and brine. Their system, with 720 water molecules, demanded around 300 ns of simulation to crystallize a few layers of ice. Two water models especially designed for ice/water systems were recently published: the six-site9 and TIP4P/Ice10 models. The six-site model has received more attention than the TIP4P/Ice and is the model that we use in this work. The main advantage of the six-site * To whom correspondence should be addressed. E-mail:
[email protected].
10.1021/jp067388q CCC: $37.00
model over previous models is that it clearly displays the growth of several layers of ice when in contact with supercooled water in a time that is fast enough to allow for detailed studies with large systems.11-14 Moreover, using this model, homogeneous nucleation has been observed in a pure system with a free interface.15 The growth rate is about 1 order of magnitude faster for the six-site model than for SPC/E8,12 and is consistent with experimental freezing rates at supercooling larger than 10 °C.16 The atomistic mechanism of growth of the six-site model has been shown to be the same as the one observed for TIP4P.7,12,13 The melting temperature of the six-site model, initially estimated at 271 ( 9 K,9 was recently recalculated using Hamiltonian Gibbs-Duhem integration and direct coexistence of the solid and liquid phases.17 The new estimation for the melting temperature is 289 K, which is certainly high and asks for reparametrization of the model. Nevertheless, this put the model in a similar territory as other models with the advantage of displaying a faster crystallization rate. As a reference, let us mention that the melting temperatures of TIP4P and SPC/E are 232 and 215 K, respectively, although the exact value depends on the details of the determination method.18 Ice Ih has a rhombic unit cell with four molecules, and ice Ic has a cubic unit cell with eight water molecules.1 Despite different crystal structures, there is a common feature among the two ices. The basal planes of ice Ih consist of two interpenetrating hexagonal double layers stacked following the ABABAB... sequence. Ice Ic consists of a stack of the same hexagonal double layers but following the ABCABC... sequence. In ice Ic, these layers coincide with the (111) family of planes. Then the basal planes of ice Ih and the (111) planes of ice Ic are identical and perfectly coherent, and therefore, it is possible to join the two ice forms along their hexagonal double layers with a very small interfacial energy.2 The enthalpy change of the irreversible transformation from ice Ic to Ih at about 200 K was experimentally determined to be 37 J/mol.19 The free energy difference between the two phases was theoretically estimated to increase from 107 J/mol at 53 K to 128 J/mol at 273 K.20 These very small values suggest that complete relaxation of a sample of mixed ice Ic and Ih to the stable phase should take © 2007 American Chemical Society
502 J. Phys. Chem. C, Vol. 111, No. 2, 2007
Letters TABLE 1: Freezing Sequence for the Six Simulated Trajectoriesa trajectory
final structure
hex-1 hex-2 hex-3 cubic-1 cubic-2 cubic-3
qh7k3h1k3q qh7k3h2k1h1q qh7k3h1k3q qk8h2k1h2k1q qk14q qk10h1k3q
a q stands for a quasi liquid layer, and h and k represent a layer in hexagonal or cubic arrangement, respectively.
Figure 1. Initial conformation for the simulations with hexagonal (top) and cubic substrate (bottom). The view is a projection on the (100) plane for the first, and on the (11h0) plane for the later. The lines indicate the different stacking arrangement in the two phases, as it crosses the hexagonal planes through the (projected) longer bonds. The zigzag pattern corresponds to hexagonal ice, and the straight line corresponds to cubic ice.
a long time compared with the time of the crystallization4 and impossible to be observed in computer simulations with today’s technology. The low interfacial energy and the coherence of the two ice types suggest that, if ice is growing on top of the hexagonal planes, there should be possible to form either one of the two structures. Namely, on an existing ice substrate having a pattern ABABAB..., the next plane could be in the A arrangement continuing the hexagonal ice or in the C arrangement as if it were cubic ice. In the same way, on a pure cubic substrate (ABCABC...) a new plane can form continuing the cubic sequence (A), or shifted to the B arrangement as if it were hexagonal ice (BCBC...). Any alteration of the order of one or the other phase is referred to as a stacking fault. In this paper, we present a molecular dynamics simulation study of ice growth from supercooled water on two different substrates. First, a block of 8 ice double layers (768 water molecules) in hexagonal arrangement (ABABABAB) with its basal plane in contact with 768 water molecules in liquid state. Second, a system of 8 ice double layers cubic arrangement (ABCABCAB) with its (111) surface in contact with liquid water. This initial setting is displayed in Figure 1. To analyze the structure of the final conformation once the system crystallized, it is convenient to use a different notation to indicate the number of layers in cubic or hexagonal arrangement.21 For example, the initial conformations displayed in Figure 1 are characterized by h8l (k8l) for ice Ih (ice Ic), where h8 (k8) means a sequence of 8 hexagonal (cubic) layers and l stands for the liquid layer. The initial coordinates for the proton disordered ice block were generated using the algorithm proposed by Sandler et al.22 Periodic boundary conditions are applied on the xyz direction, and the simulations are performed at constant temperature (Nose-Hoover thermostat) and semi-isotropic constant pressure
Figure 2. Structure after 15 ns of simulation for hexagonal (top) and cubic (bottom) substrates. There are statistically identical quasi liquid layer on all of the surfaces, regardless of the ice type. The ice formed during the simulation is a mixture of cubic and hexagonal ice, regardless of the substrate.
(Parrinello-Rahman barostat). The size of the simulation box in the z direction was kept at a constant value of 12 nm during the simulations. This creates an ice-vacuum (water-vacuum) interface on the left (right) side of the initial system displayed in Figure 1. Once the simulation started, the vacuum region should be regarded as the vapor in equilibrium with the system. This general setting is similar to the one used in ref 12. The advantage of having the interface with vacuum (or vapor) instead of using full contact with the periodic boundary conditions is that in this case we have only one ice growth front, which eliminate the structural constraints created when the system is freezing with two approaching growth fronts. Besides, it allows for the study of the quasi liquid layer that forms on the surface of ice at temperatures close to melting. All of the results presented here correspond to simulations performed at 275 K, which represent a supercooling of 14 °C. Three simulations of 15 ns were performed for each system. The stacking sequences of the final configurations are summarized in Table 1, and the final configuration for one run for each system is depicted in Figure 2. The main result is that regardless of the initial substrate, stacking faults appear during the crystallization. Moreover, sections of cubic ice have grown on the hexagonal substrate, and sections of hexagonal ice have grown on the cubic substrate. By looking at the stacking sequence of the final structures, we see that there are more layers that have formed following the cubic arrangement than layers in the hexagonal arrangement. Even though the size of our study
Letters
J. Phys. Chem. C, Vol. 111, No. 2, 2007 503
Figure 4. Potential energy as a function of time for a system with hexagonal (black, solid lines) and cubic (red, dashed lines) substrates. The horizontal lines indicate the energy before and after the rapid crystallization of a complete hexagonal plane.
Figure 3. Snapshots of two ice layers that crystallized in cubic (top) and hexagonal (bottom) sequence. The left panels are projections on the same plane as in Figure 1 and the right panels are projections on the xy plane.
is very limited, it suggests that ice Ic is the most likely structure to form upon freezing. Out of 42 ice layers created during the growth simulation, 32 of them followed a cubic sequence. This means that 76% of the new ice grew in cubic structure. The stability of the conformations obtained after 15 ns was checked by extending the simulation (for one case) to a total of 100 ns. No modification of the structure was observed, except for the melting and refreezing of the external ice layer (see below). The fact that the system prefers to form the metastable phase instead of the stable phase should not be surprising, since it is a rather common phenomena observed in many systems and referred to as the Ostwald step rule.23 In particular, for ice crystallization, our results are in agreement with the recent experiments of Murray et al.4 In Figure 3, we show a detail of two layers in the cubic sequence formed in the system with hexagonal substrate (top) and two layer in hexagonal sequence that crystallized on the system with the cubic substrate (bottom). The pictures on the left are projections similar to the ones displayed on Figure 2, and the pictures on the right are projections (of the same atoms displayed on the left pictures) on the xy plane which is the (111) plane for ice Ic and the basal plane for ice Ih. The snapshots display the characteristics of the two structures. In the cubic case (top), the hexagonal patterns are shifted so that the oxygen atoms of one hexagonal plane are at the center of the hexagons of the adjacent hexagonal sublayer. In the hexagonal case, the hexagonal planes are one on top of each other allowing seeing through the hexagonal channels created by the oxygen atoms. The relative orientation of the “longer” bonds in the projections displayed on the left helps to identify the structure. In the cubic case, the long bonds appear to be parallel. In the hexagonal case, the long bonds form an angle of approximately 20°.
In Figure 4, we display the potential energy as a function of time for the two studied systems. For the sake of clarity, we show only one representative case for each system. The final energy of the system is very similar among all of the studied systems, and the difference depends on the particular freezing sequence of each trajectory. For the cases displayed in Figure 4, the energy difference is 0.1 kJ/mol and is the largest one observed for all of the systems. The mechanism of freezing on the basal plane of hexagonal ice has been identified as a layerby-layer process.7 This fact is reflected in the energy vs time plot by the sharp steps occurring when a layer crystallizes. The size of these steps is 0.5 kJ/mol and the same for the two structure types within the statistical error. After 11 ns of simulation the system has frozen seven new ice layers and melted one layer on the (initial) ice/vacuum interface. Namely, at 11 ns, the system consists of an ice slab wetted at both sides by quasi liquid layers. The increase of the potential energy of the cubic system after 12 ns of simulation corresponds to the melting of one of the external ice layers. The melted layer refreezes after approximately 1 ns. Therefore, the quasi liquid layer is thicker on one side of the system for a short period of time. This behavior has been observed in several cases, independently of the initial substrate and confirming the dynamic nature of the quasi liquid layers. The size of the quasi liquid layer is the same for both structures and consistent with our previous study of pure hexagonal ice.12 In summary, we report molecular dynamics simulations that show that stacking faults naturally appear during freezing occurring on the basal plane of ice Ih and on the (111) plane of ice Ic. The freezing mechanism is the same independently of the freezing sequence and can be described as a layer by layer process. The rapid completion of each layer is reflected by a sharp jump in the potential energy, as the system decreases in enthalpy upon freezing. Our results are in agreement with recent experimental findings that showed that ice Ic is the dominant crystalline phase immediately after freezing. References and Notes (1) Hobbs, P. V. Ice Physics; Clarendon Press: Oxford, 1974. (2) Johari, G. P. Phil. Mag. B. 1998, 78, 375-383. (3) Murray, B. J.; Knopf, D. A.; Bertram, A. K. Nature 2005, 434, 202-205. (4) Murray, B. J.; Bertram, A. K. Phys. Chem. Chem. Phys. 2006, 8, 186-192. (5) Morishige, K.; Uematsu, H. J. Chem. Phys. 2005, 122, 044711. (6) Matsumoto, M.; Saito, S.; Ohmine, I. Nature 2002, 416, 409-413. (7) Nada., H.; Furukawa, Y. J. Cryst. Growth 1996, 169, 587-597. (8) Vrbka, L.; Jungwirth, P. Phys. ReV. Lett. 2005, 95, 148501.
504 J. Phys. Chem. C, Vol. 111, No. 2, 2007 (9) Nada, H.; van der Eerden, J. P. J. M. J. Chem. Phys. 2003, 118, 7401-7413. (10) Abascal, J. L. F.; Sanz, E.; Ferna´ndez, R. G.; Vega, C. J. Chem. Phys. 2005, 122 (23), 234511. (11) Nada, H.; van der Eerden, J. P.; Furukawa, Y. J. Cryst. Growth 2004, 266, 297-302. (12) Carignano, M.; Shepson, P.; Szleifer, I. Mol. Phys. 2005, 103, 2957-2967. (13) Nada, H.; Furukawa, Y. J. Cryst. Growth 2005, 283, 242-256. (14) Carignano, M.; Baskaran, E.; Shepson, P.; Szleifer, I. Ann. Glaciol. 2006, 44, in press. (15) Vrbka, L.; Jungwirth, P. J. Phys. Chem. B 2006, 110, 18126-18129. (16) Pruppacher, H. R. J. Chem. Phys. 1967, 47, 1807-1813.
Letters (17) Abascal, J. L. F.; Ferna´ndez, R. G.; Vega, C.; Carignano, M. A. J. Chem. Phys. 2006, 125, 166101. (18) Vega, C.; Sanz, E.; Abascal, J. L. F. J. Chem. Phys. 2005, 122, 114507. (19) Yomamuro, O.; Oguni, M.; Matsuo, T.; Suga, H. J. Phys. Chem. Solids 1987, 48, 935-942. (20) Tanaka, H. J. Chem. Phys. 1998, 108, 4887-4893. (21) Sebastian, M. T.; Krishna, P. Random, Non-Random and Periodic Faulting in Crystals; Gordon and Beach Science Publishers: Yverdon, Switzerland, 1994. (22) Sandler, P.; oh Jung, J.; Szczesniak, M. M.; Buch, V. J. Chem. Phys. 1994, 101, 1378-1392. (23) Ostwald, W. Z. Phys. Chem. 1879, 22, 289-330.
