Anisotropy in Growth Kinetics of Tetrahydrofuran Clathrate Hydrate: A

Mar 19, 2009 - The growth kinetics of a tetrahydrofuran (THF) clathrate hydrate at the interface between the clathrate and an aqueous THF solution wer...
0 downloads 0 Views 2MB Size
Download PDF
4790

J. Phys. Chem. B 2009, 113, 4790–4798

Anisotropy in Growth Kinetics of Tetrahydrofuran Clathrate Hydrate: A Molecular Dynamics Study Hiroki Nada* National Institute of AdVanced Industrial Science and Technology (AIST), 16-1 Onogawa, Tsukuba 305-8569, Japan ReceiVed: NoVember 14, 2008; ReVised Manuscript ReceiVed: February 10, 2009

The growth kinetics of a tetrahydrofuran (THF) clathrate hydrate at the interface between the clathrate and an aqueous THF solution were investigated by means of a molecular dynamic simulation. The simulation was carried out for the interface of both the {100} and {111} planes of the THF clathrate. The simulation indicated the same anisotropic growth as that observed in real systems: the growth of the THF clathrate was much slower at the {111} interface than at the {100} interface. When the THF clathrate grew, THF molecules that were dissolved in the solution first were arranged at both large and small cage sites on the interface. Subsequently, the formation of cages by H2O molecules occurred in regions surrounded or sandwiched by those arranged THF molecules. As the formation of cages progressed, the THF molecules that had once been arranged at small cage sites gradually moved away from the sites, and finally the structure of the clathrate was completely formed. Simulation results strongly suggested that the rate-determining process for clathrate growth was the rearrangement of THF molecules at the interface from a disordered state to a state in which THF molecules were ideally arranged at large cage sites only. This rearrangement occurred much more slowly at the {111} interface than at the {100} interface, owing to the formation of a modified structure in which large and small cages were formed at opposite positions of the {111} interface. The anisotropic growth kinetics of the THF clathrate, which were obtained in this study, are consistent with the fact that growth shapes of THF clathrates in real systems are octahedral with flat {111} planes. 1. Introduction Clathrate hydrates (hereafter, clathrates) are crystalline forms of water that contain many guest molecule inclusions (hereafter, guests).1 Each guest in a clathrate is included in a cage formed by water molecules. Many types of gaseous molecules can be guests of clathrates. Therefore, clathrates can be used to store fuel gases, such as CH4 and H2, and to separate and remove greenhouse gases, such as CO2, from mixtures of gas species released from factories.2,3 Thus, clathrates have attracted a great deal of attention from researchers who work in the fields of physical chemistry, chemical engineering, energy technology, and environmental technology. A problem associated with the aforementioned practical uses of gas clathrates is that gas clathrates are not thermodynamically stable at 1 atm pressure near 0 °C, and, therefore, the growth of gas clathrates requires applying high pressure to a twocomponent system of water and gas, externally cooling the system, or both. Unfortunately, applying high pressure involves risks of explosions, and cooling consumes large amounts of electricity, which is unfavorable because such consumption leads to release of CO2 greenhouse gases. Therefore, a method by which to capture gas molecules into clathrates without applying high pressure or external cooling is highly desired. Clathrate guests are not restricted to gas molecules. For example, tetrahydrofuran (THF) molecules can also be guests of clathrates.4 THF is an organic solvent that is widely used in scientific and industrial works. The melting point and boiling point of THF at 1 atm are 164.75 and 339.15 K, respectively. Unlike gas molecules, THF molecules are soluble in liquid water * To whom correspondence should be addressed. Telephone: +81-29861-8231. Fax: +81-29-861-8722. E-mail: [email protected].

at an arbitrary concentration. Moreover, in contrast to gas clathrates, THF clathrates grow from an aqueous THF solution at temperatures above 0 °C, even at 1 atm (the melting point, Tm, of THF clathrates at 1 atm is 277.55 K).5 Thus, neither high pressure nor cooling below 0 °C is required for the growth of THF clathrates. THF clathrates are structure II clathrates consisting of large (hexakaidecahedra consisting of 12 pentagons and four hexagons, 51264) and small (pentagonal dodecahedra, 512) cages.4 Experimental results indicate that THF molecules are included as guests in the large cages only.4 However, when THF clathrates are grown from a three-component system of water, THF, and a certain species of gas, gas molecules can also be included as guests in vacant cages of the grown THF clathrates, provided that the sizes of the gas molecules are appropriate for inclusion in the vacant cages.6-12 Notably, in such instances, gas molecules are captured in THF clathrates at much lower pressures or much higher temperatures than those required for gas molecule inclusion into pure gas clathrates. Therefore, THF clathrates have received much attention for use as materials for storage of fuel gases and for separation and removal of greenhouse gases.6-11 For such practical applications, knowledge on the growth kinetics of THF clathrates is of particular importance. So far, several experimental studies of the growth of THF clathrates have been carried out.13-17 These studies indicate that the growth shape of THF clathrates is octahedral with flat {111} planes, meaning that the {111} planes have the slowest growth velocity of all crystallographic planes. The difference in the growth velocity of THF clathrates among crystallographic planes originates from the anisotropy in the growth kinetics of the

10.1021/jp810041t CCC: $40.75  2009 American Chemical Society Published on Web 03/19/2009

Growth Kinetics of Tetrahydrofuran Clathrate Hydrate planes. However, for all crystallographic planes of THF clathrates, the growth kinetics remain unclear because elucidation of these growth kinetics by experimental means is quite difficult. At present, computer simulations such as molecular dynamics (MD) are the only methods that can be used to analyze the growth kinetics of THF clathrates at the molecular level.18 So far, several MD simulation studies of the growth kinetics of CH4 clathrates have been carried out.19-25 However, the structure of CH4 clathrates differs from that of THF clathrates in that CH4 clathrates are structure I, in which CH4 molecules are included as guests in both large (tetrakaidecahedra consisting of 12 pentagons and two hexagons) and small (512) cages. Moreover, since the solubility of CH4 molecules in liquid water is very low, CH4 clathrates grow from a dilute aqueous CH4 solution. However, THF clathrates normally grow from a concentrated aqueous THF solution. Therefore, the growth kinetics of THF clathrates, which may be different from those of CH4 clathrates, should be independently studied by means of MD simulations. To our knowledge, no one has undertaken such studies to date. In this study, we carried out MD simulations to elucidate the growth kinetics on the {100} and {111} planes of a THF clathrate. The simulations indicated that the growth of a THF clathrate was much slower on the {111} plane than on the {100} plane, in agreement with experimental results.13-17 In this article, we discuss the anisotropy in the molecular-scale growth kinetics between the planes, as well as the relationship between this anisotropy and the growth shapes of THF clathrates in real systems. 2. Simulation Methods 2.1. Potential Models. The intermolecular interaction between a pair of H2O molecules was estimated with a six-site model,26 which is a rigid H2O model developed for simulations of ice and water near the real Tm of ice at 1 atm. In this model, a H2O molecule has six interaction sites: an O atom, two H atoms, two lone-pair sites (L), and a site M located on the bisector of ∠HOH. A positive charge is placed on each H atom, and a negative charge on each L and M site. The interaction is calculated as the sum of the Coulomb potentials between the charges plus the sum of the Lennard-Jones (LJ) potentials between the atoms. The six-site model reproduces the growth of an ice crystal from water near the real Tm of ice at 1 atm in MD simulations.27-30 Moreover, using the six-site model, we observed the growth of a CH4 clathrate at 298 K and 100 MPa in a previous MD simulation.24 Therefore, we believe that the sixsite model is also suitable for MD simulations of the growth of a THF clathrate. For THF molecules, a five-site model proposed by Chandrasekhar and Jorgensen was used.31 In this model, each methyl group is represented by a single interaction site. Thus, this model is a simplification of a real THF molecule. Nevertheless, this model satisfactorily reproduces the structural properties of liquid THF.31 Moreover, as we will describe below, the Tm of a THF clathrate in this THF model with the six-site model proved to be higher than the Tm of ice in the six-site model, as is the case in real systems. Therefore, this simplified THF model is believed to be sufficient for qualitatively understanding the growth kinetics of a THF clathrate. For quantitatively reproducing the growth of a THF clathrate, which was not the purpose of this study though, it would be better to use a fully atomistic model.9 The interactions between a pair of THF molecules and between a pair of H2O and THF molecules were also calculated

J. Phys. Chem. B, Vol. 113, No. 14, 2009 4791

Figure 1. Molecular geometries in the six-site H2O and five-site THF models. Large gray and small black circles are the LJ and Coulomb interaction sites, respectively. Large gray circles in which a small black circle is drawn indicate sites on which both the LJ and the Coulomb interactions act. The charge amounts for each site on which the Coulomb interaction acts are shown in parentheses. The values of ε/kB and σ for the LJ interactions for each model are given in Table 1. For a H2O molecule, the values of O-H distance, O-M distance, O-L distance, ∠HOH, and ∠LOL are 0.098 nm, 0.023 nm, 0.08892 nm, 108°, and 111°, respectively. For a THF molecule, the values of O-C distance, C-C distance, ∠COC, and ∠OCC are 0.1411 nm, 0.1529 nm, 111°, and 109.4°, respectively.

TABLE 1: Values of ε/kB and σ for the LJ Interaction Sites of the Six-Site H2O and Five-Site THF Moleculesa molecules

sites

σ (nm)

ε/kB (K)

H2O-H2O

Ow-Ow Ow-H H-H OT-OT OT-CH2 CH2-CH2 Ow-OT Ow-CH2 H-OT H-CH2

0.3115 0.1894 0.0673 0.3000 0.3453 0.3905 0.3058 0.3510 0.1837 0.2289

85.9766 34.5471 13.8817 85.4665 71.2462 59.3920 85.7212 71.4585 34.4445 28.7134

THF-THF H2O-THF

a Ow and OT represent O of H2O and THF molecules, respectively.

as the sum of the Coulomb potentials plus the sum of the LJ potentials. Parameters for the LJ potentials between H2O and THFmoleculesweredeterminedusingthestandardLorentz-Berthelot rules. Illustrations of molecular geometries and the values of the LJ parameters in the models are given in Figure 1 and Table 1, respectively. Notably, the parameters of the six-site model were optimized by assuming the truncation of intermolecular interactions at intermolecular distances around 1 nm.26 Therefore, in this study, as with previous studies,24,27-29 all intermolecular interactions were smoothly truncated at intermolecular distances from 0.95 to 1 nm by means of a switching function.32 In a previous study, the anisotropic growth kinetics of ice from water, which explained a growth shape of ice in real systems qualitatively, were successfully obtained using the six-site model with the truncation of the interactions.28 Therefore, we believe that the anisotropic growth kinetics of a THF clathrate can also be qualitatively elucidated even if the truncation of the interactions is used. 2.2. Simulation System. We prepared two simulation systems (Figure 2a). One system simulated growth on the {100} planes of a THF clathrate ({100} system), and the other simulated growth on the {111} planes ({111} system). Each system was a rectangular parallelepiped in which a THF clathrate was sandwiched between two stoichiometric solution phases (i.e., solution phases that have the same concentration of THF molecules as that in the clathrate). Each system contained two interfaces between the clathrate and the solution phase (hereafter, these interfaces are referred to as the clathrate

4792 J. Phys. Chem. B, Vol. 113, No. 14, 2009

Figure 2. (a) Illustration of the simulation system. The dotted lines show the clathrate interfaces. (b) Illustrations of large (51264) and small (512) cages in the clathrate. The solid lines show hydrogen bonds between H2O molecules. The large cage contains a single THF molecule, which is illustrated by red and orange spheres that correspond to methyl groups and the oxygen atom, respectively. (c) Illustrations of the clathrates in the {100} and {111} systems. The blue-colored open circles show the positions of empty small cage sites. The structures of a {100} layer and of the R and β layers of a {111} layer, which are projected onto the x-y plane, are also illustrated.

interfaces). Periodic boundary conditions were imposed in all x, y, and z directions of the system. At the beginning of the simulation, all H2O molecules in the clathrate were arranged on the ideal lattice sites of the structure II clathrate.4 The structure II clathrate contains eight large (51264) and 16 small (512) cages in its unit cell (Figure 2b). Because experimental results indicate that only large cages are occupied by a THF molecule in real THF clathrates,4 we positioned a THF molecule at the center of each large cage only. The clathrate for the {100} system consisted of 1088 H2O and 64 THF molecules, and that for the {111} system consisted of 1224 H2O and 72 THF molecules. For both systems, each of the two solution phases in the system consisted of the same number of H2O and THF molecules as those in the clathrate. Thus, the total numbers of molecules included in the {100} and {111} systems were 3456 and 3888, respectively. The sizes of the systems used in this study might not be sufficient for quantitative elucidation of the growth of the clathrate. However, a MD simulation study for the growth of ice from water using the six-site model by Carginano et al. indicated that the structure of an ice-water interface for a 6144 H2O system was the same as that for a 2536 H2O system.30 Moreover, the anisotropic growth kinetics of ice that they obtained for the 2536 H2O system was qualitatively the same as that we obtained for a 1440 H2O system in a previous study,28,30 suggesting that those sizes were sufficient for qualitative elucidation of the anisotropic growth kinetics of ice. Therefore, we believe that the sizes of the systems used in this study were also sufficient, at least, for qualitative elucidation of the anisotropic growth kinetics of a THF clathrate. In order to create each system including two clathrate interfaces, we used the same method as used in previous studies.27,28 First, a MD simulation of the clathrate was carried out to obtain its equilibrium structure at 285 K and 1 atm. Second, an equilibrium structure of the solution phase at 285 K and 1 atm was obtained from a MD simulation in which the dimensions of the system in the x and y directions were fixed at the equilibrium values for the clathrate (3.45 × 3.45 nm2 for the {100} system and 4.23 × 3.66 nm2 for the {111} system). Third, the clathrate was sandwiched with two copies of the solution phase to create two clathrate interfaces. Finally, in order to relax the structures of the clathrate interfaces, we carried out a short Monte Carlo (MC) simulation at 285 K for each system

Nada using a standard Metropolis sampling method for the NVT ensemble. The number of trial displacements of molecules in the MC simulation was 3.456 × 106 for the {100} system and 3.888 × 106 for the {111} system. By using the MC simulation, disruption of the structure of the clathrate due to extremely large intermolecular forces at the clathrate interfaces, which might occur if a MD simulation was used for the relaxation, was successfully avoided. We checked that the growth of the clathrate did not occur during the MC simulation. The equilibrium structure of the clathrate in each system before the relaxation is illustrated in Figure 2c. For the {100} system, eight of a unit layer, which has a thickness of approximately 0.43 nm (see the region enclosed by solid lines for the {100} system in Figure 2c) are stacked in the z direction. Each unit layer includes eight large and 16 empty small cage sites. Hereafter, those unit layers are referred to as {100} layers. The arrangement of THF molecules in a {100} layer, which is projected onto the x-y plane, is shown in Figure 2c. For the {111} system, three of a unit layer, which has a thickness of approximately 1.0 nm (see the region enclosed by solid lines for the {111} system in Figure 2c), are stacked in the z direction. Each unit layer includes 24 large and 48 empty small cage sites. Hereafter, those unit layers are referred to as {111} layers. Moreover, we also divided each {111} layer into three thin layers: a layer containing 12 large cage sites, a layer containing 12 large and 12 empty small cage sites, and a layer containing 36 empty small cage sites. Hereafter, those divided layers are referred to as the R, β, and γ layers, respectively.33 The arrangements of THF molecules in the R and β layers, which are projected onto the x-y plane, are shown in Figure 2c. (The γ layer is not shown.) 2.3. MD Simulation. MD simulations were carried out for each system by means of a leapfrog algorithm proposed by Fincham34 with a time step of 1 fs. The total run was 20 ns. Temperature and pressure were maintained at 285 K and 1 atm, respectively, by means of a method proposed by Berendsen et al.35 The thermal and pressure bath constants were set to 0.1 and 2.0 ps, respectively.35 During the simulation, the pressure was maintained by changing the dimension of the systems in the z direction only. The dimensions of each system in the x and y directions were fixed at their equilibrium values, 3.45 × 3.45 nm2 for the {100} system and 4.23 × 3.66 nm2 for the {111} system. These equilibrium values were determined by a separate MD simulation of a bulk THF clathrate at 285 K and 1 atm. Notably, the Tm of the clathrate at 1 atm in the potential models used here was higher than the real Tm of 277.6 K,5 because we observed clathrate growth at 285 K. However, at 295 K, we did not observe clathrate growth. At T > 300 K, we observed melting of the clathrate clearly. Therefore, we assumed that the Tm of the clathrate at 1 atm in this simulation was 295 K. Hence, the simulation temperature of 285 K corresponded to 0.966Tm, which is equivalent to 268 K in real systems. Extensive free energy calculations for both the clathrate and the solution are needed to determine the Tm of the clathrate precisely. 3. Results and Analysis 3.1. Anisotropic Growth Velocity of the THF Clathrate. Figure 3a shows the change in potential energy (U), ∆U ) U(t) - U(t)0), as a function of time, t, for both systems. A decrease in ∆U indicates the growth of the clathrate, and the slope of the ∆U-t relationship is proportional to the growth velocity of the clathrate. Obviously, the growth velocity was much larger

Growth Kinetics of Tetrahydrofuran Clathrate Hydrate

J. Phys. Chem. B, Vol. 113, No. 14, 2009 4793

Figure 4. (a) Trajectories of THF molecules for periods of 9.5-10.0, 11.5-12.0, and 13.5-14.0 ns, and snapshots of H2O molecules in region B for the {100} system. The snapshots were created from the time-averaged coordinates of H2O molecules over periods of 9.9-10.0, 11.9-12.0, and 13.9-14.0 ns. (b) The illustration at the left shows a potential energy surface created by the ideal structure of a {100} layer in region A. The illustration at the right shows the ideal positions of large (red) and small (blue) cage sites in region B.

Figure 3. (a) Change in potential energy (U), ∆U ) U(t) - U(t)0), of the clathrate systems as a function of time, t. (b) Snapshots of THF and H2O molecules near one of the clathrate interfaces in each system. The snapshots were created from the time-averaged coordinates of H2O molecules over periods of 0.9-1.0, 9.9-10.0, and 19.9-20.0 ns. The dotted lines show the initial positions of the clathrate interfaces. The blue molecules represent sl-H2O molecules.

for the {100} system than that for the {111} system. This result qualitatively agrees with experimental results.13-17 Figure 3b shows the time sequence of snapshots of molecules near one of the clathrate interfaces for both systems. The snapshots were created from the time-averaged coordinates of molecules over 0.9-1.0, 9.9-10.0, and 19.9-20.0 ns, respectively. The dotted lines show the initial positions of the clathrate interfaces. In the snapshots, the blue molecules represent “solidlike” H2O molecules (sl-H2O molecules), which were defined as H2O molecules that had four nearest neighboring H2O molecules and were connected by a hydrogen bond with each of the nearest neighboring H2O molecules.26,28 The sl-H2O molecules were analyzed from the time-averaged coordinates that were used for creating the snapshots. Growth of about three {100} layers occurred at the initial position of the clathrate interface for the {100} system (regions A, B, and C in Figure 3b). The thickness of each of regions A, B, and C was the same as the thickness of a {100} layer. However, for the {111} system, growth of only a single, incomplete {111} layer was observed at the initial position of the clathrate interface. Growth of about three {100} layers for the {100} system and growth of only a single, incomplete {111} layer for the {111} system were also observed at the initial positions of the other clathrate interfaces. Therefore, in this report, we present the results of growth kinetics analysis for only one clathrate interface for each system, i.e., the clathrate interfaces shown in Figure 3b. Notably, we carried out an additional MD simulation starting with different initial molecular velocities for each system, and

we observed qualitatively the same anisotropic growth kinetics as described below. Unfortunately, the growth observed for the {111} system in this additional simulation was also insufficient to allow evaluation of the growth velocity. A MD simulation with a much longer run than in this study is required to evaluate the growth velocity. 3.2. Growth Kinetics on the {100} Plane. First, we will show the results for the {100} system. Figure 4a shows the time sequence of the trajectories of THF molecules and the snapshots of H2O molecules in region B for the {100} system. Interestingly, THF molecules were arranged not only at large cage sites but also at small cage sites, even though all small cage sites in bulk THF clathrates are known to be empty. Moreover, the THF molecules were stably arranged at those cage sites, even for cases in which the surrounding H2O molecules had not yet formed cages. Figure 4b shows a potential energy surface created by an ideal structure of region A. The potential energy surface was determined in the following way. For a particle that was located on the ideal structure of region A, the minimum value of U between the ideal structure of region A and the particle was searched for by changing the z component of the particle while fixing its x and y components. This procedure was repeated for every x and y component of the particle. U was calculated from the LJ potential, and the LJ parameters of the particle were set to be the same as those of the methyl groups of a THF molecule. The purpose of determining the potential energy surface was to infer the reason for the stable arrangement of THF molecules at cage sites prior to the formation of cages by the surrounding H2O molecules. Of course, this potential energy surface determined for a LJ particle is not exactly the same as that determined for a THF molecule. However, we analyzed U of the LJ and Coulomb terms, ULJ and UCoul, acting on a THF molecule in a bulk THF clathrate, and ULJ () -3.1 kJ/mol) proved to be much lower than UCoul () -0.3 kJ/mol). This result suggests that the stability of the arranged THF molecules at cage sites is mainly dominated by the LJ potentials. Moreover, we checked that the energetically stable positions of a THF molecule on the fixed ideal structure of region A, which were searched for by MC simulations at 0 K, were consistent with those predicted from the potential energy surface for a LJ

4794 J. Phys. Chem. B, Vol. 113, No. 14, 2009

Nada

Figure 6. Distribution function, P(θ), of the angle θ between a vector normal to the THF molecular plane (n) and the z axis for THF molecules that were arranged at cage sites on the clathrate interface for the {100} system. The dashed line shows P(θ) ) 1/2 sin θ, which is P(θ) for random orientations of THF molecules.

Figure 5. Numbers of THF molecules that were arranged at large and small cage sites, NgL and NgS, and fraction of sl-H2O molecules, fsl, as a function of time in region B for the {100} system.

particle. Thus, we confirmed that the potential energy surface for a LJ particle was sufficient to infer the reason for the stable arrangement of THF molecules at cage sites. Local minima appeared at both large and small cage sites, though the minimum energy values at small cage sites were lower than those at large cage sites. This appearance of local minima at both large and small cage sites explains why THF molecules were arranged at both large and small cage sites prior to the formation of cages by the surrounding H2O molecules. Notably, three pentagons constituting each of the small cages in region B were contained in the ideal structure of region A, whereas, for each of the large cage sites in region B, only two pentagons were contained in the ideal structure of region A. The lower minimum energy values at small cage sites relative to the values at large cage sites originated from this larger number of pentagons that were included in the ideal structure of region A for each of the small cage sites in region B than the number for each of the large cage sites in region B. Figure 5 shows the time sequences of the numbers of THF molecules that were arranged at large and small cages (NgL and NgS, respectively) and the fraction of the sl-H2O molecules, fsl, in region B. The result clearly indicates that when fsl increased, thus marking the onset of cage formation, about 50% of the large and small cage sites had already been occupied by a THF molecule. The decrease in NgS with increasing fsl, which occurred after 10 ns, indicates that as the formation of cages progressed, the arrangement of the THF molecules at small cage sites was gradually destabilized, and the THF molecules finally moved away from the sites. Similar time evolutions of NgL, NgS, and fsl were also observed for growth in regions A and C. The orientations of THF molecules at large and small cage sites at the clathrate interface reflected the differences in the stability of the THF molecules’ arrangement between the two types of cage sites during cage formation. Figure 6 shows the distribution function, P(θ), of the angle θ between a vector normal to the molecular plane of a THF molecule (n) and the z axis for THF molecules that were arranged at cage sites at the clathrate interface for the {100} system.36 The dashed line shows P(θ) ) 1/2 sin θ, which is P(θ) for random orientations

Figure 7. Number density profiles, Fz, of H2O and THF molecules and fsl along the z direction for the {100} system. Fz and fsl were analyzed from simulation data for a period of 9.5-10.0 ns. The origin of the z component corresponds to the position of the clathrate interface at the beginning of the simulation.

of THF molecules. P(θ) was measured for THF molecules that were arranged at cage sites in region B for a period of 10-13 ns. Obviously, P(θ) for THF molecules at large cage sites indicates a random distribution of orientations, in agreement with previous experimental results showing that the orientations of THF molecules are isotropically distributed in large cages of bulk THF clathrates.37 However, THF molecules at small cage sites tended to take orientations in which the molecular planes became parallel to the z axis (i.e., θ tended to be 90°). Such orientations suggest that only one of the interaction sites (four methyl groups and an oxygen atom) of each THF molecule was located at a small cage site at the clathrate interface, whereas the other four interaction sites tended to be exposed to the solution. These orientations reflected the destabilization of the arrangement of THF molecules at small cage sites owing to the formation of small cages around the sites: since small cages were too small to include THF molecules, as the formation of small cages progressed, the THF molecules changed their orientations so as to avoid repulsive interactions with the formed cages. Notably, the arrangement of THF molecules at cage sites at the clathrate interface was not sufficient for the formation of entire cages around the molecules. As can be seen in Figure 3b, most of the sl-H2O molecules appeared in regions surrounded or sandwiched by THF molecules that had already been arranged at cage sites. The formation of cages mainly in regions surrounded or sandwiched by arranged THF molecules is also seen in Figure 7, which shows the number density profiles (Fz) of H2O and THF molecules and the profile of fsl in the z direction for the period of 9.5-10 ns for the {100} system. The origin of the z component corresponds to the position of the clathrate interface at the beginning of the simulation. H2O molecular layers are observed only inside of the outermost THF molecular layer, which is included in region B. Moreover, the values of

Growth Kinetics of Tetrahydrofuran Clathrate Hydrate

Figure 8. (a) Trajectories of THF molecules for periods of 5.5-6.0, 11.5-12.0, and 17.5-18.0 ns, and snapshots of H2O molecules in the R and β layers of a {111} layer at the initial position of the clathrate interface for the {111} system. The snapshots were created from the time-averaged coordinates of H2O molecules over periods of 5.9-6.0, 11.9-12.0, and 17.9-18.0 ns. (b) The illustration at left shows a potential energy surface created by the ideal structure of the R layer. The illustration at right shows the ideal positions of large (red) and small (blue) cage sites in the β layer.

fsl in the region outside of the outermost THF molecular layer are much closer to the fsl value for the solution (0) than that for the clathrate (1). Strictly speaking, partial structures of cages transiently appeared in regions outside the outermost THF molecular layer. However, those structures were not stably formed, and thus those structures were not regarded as parts of the clathrate. Therefore, we conclude that the formation of cages occurred only in regions surrounded or sandwiched by THF molecules that were arranged at cage sites on the clathrate interface. Similar results for the formation of cages were also obtained in a previous study of the growth of gas clathrates.24 In summary, the growth kinetics were as follows. First, THF molecules were arranged at both large and small cage sites on the clathrate interface, according to the local minima of a potential energy surface created by the underlying {100} layer. Subsequently, the formation of cages by H2O molecules occurred in regions surrounded or sandwiched by arranged THF molecules. As the formation of cages progressed, THF molecules that were once arranged at small cage sites gradually moved away from the sites. Finally, when all large cage sites were ideally occupied by a THF molecule, the formation of a {100} layer was completed. 3.3. Growth Kinetics on the {111} Plane. Next, we present the results for the {111} system. During the growth of the clathrate for the {111} system, the arrangement of THF molecules at cage sites at the clathrate interface prior to the formation of cages by the surrounding H2O molecules was observed. The formation of cages in regions surrounded or sandwiched by the arranged THF molecules was observed, as well. However, the growth kinetics of the two systems differed remarkably, which explains the anisotropic growth velocity between the systems. Figure 8a shows the time sequence of the trajectories of THF molecules and the snapshots of the sl-H2O molecules in the R and β layers at the initial position of the clathrate interface. In the β layer, the arrangement of THF molecules at both large and small cage sites prior to the formation of cages is clearly seen. This result also is attributable to the appearance of local minima of a potential energy surface created by the underlying R layer (Figure 8b). As can be seen in Figure 8a, although the

J. Phys. Chem. B, Vol. 113, No. 14, 2009 4795

Figure 9. Time sequences of NgL, NgS, and fsl in the β layer for the {111} system.

Figure 10. Snapshots of THF and H2O molecules near the other clathrate interface in the {111} system. The snapshots were created from the time-averaged coordinates of H2O molecules over periods of 3.9-4.0, 6.9-7.0, and 9.9-10.0 ns.

R layer had almost completely grown before 6 ns (strictly speaking, the growth progressed quickly and had been almost completed at 3 ns), the growth of the β layer after 6 ns was very slow. This very slow growth of the β layer after 6 ns is also seen in Figure 9, which shows that NgL, NgS, and fsl for the β layer became almost constant after 6 ns. Notably, as can be seen in Figure 9, fsl for the β layer increased immediately after the beginning of the simulation, suggesting that the β layer started growing immediately after the beginning of the simulation. This result would reflect that since the growth of the R layer occurred immediately after the beginning of the simulation, the subsequent growth of the β layer also occurred on the grown parts of the R layer immediately, before the growth of the R layer was completed. For the other clathrate interface, the growth of both the γ and R layers was finished within a period of 10 ns from the beginning of the simulation (Figure 10). However, the subsequent growth of the β layer was very slow as well. Consequently, we conclude that the slow growth of the clathrate for the {111} system was due to the very slow growth of the β layer. Therefore, elucidation of the growth kinetics of the β layer is key to understanding why the clathrate grew much more slowly for the {111} system than for the {100} system. As can be seen in Figure 2c, in the ideal structure of the β layer, hydrogen bonds between H2O molecules create a triangle lattice in which THF molecules are arranged at the triangle lattice sites (large cage sites). However, in some regions of the β layer, the triangle lattice appeared to be shifted from its ideal

4796 J. Phys. Chem. B, Vol. 113, No. 14, 2009

Figure 11. P(θ) for THF molecules that were arranged at cage sites in the β layer for the {111} system. The definition of the angle θ is the same as that given in Figure 6. The dashed line shows P(θ) ) 1/2 sin θ, which is P(θ) for random orientations of THF molecules.

position (see the regions enclosed by dotted lines in Figure 8a). This result suggests that a modified structure was formed in those regions. In the triangle lattice of those regions, THF molecules were arranged at small cage sites. Therefore, in those regions, large and small cages were formed at opposite positions: H2O molecules near small cage sites were arranged to form large cages, and H2O molecules near the adjacent large cage sites were arranged to form small cages. The formation of this modified structure was also reflected in P(θ). Figure 11 shows P(θ) for THF molecules that were arranged at cage sites of the β layer in the {111} system for a period of 10-13 ns. Unlike P(θ) for the {100} system, P(θ) for large cage sites was almost the same as that for small cage sites. The P(θ) values for both small and large cage sites were intermediate between the P(θ) for the random orientation and that for orientations in which THF molecular planes became parallel to the z axis. Notably, during the period of 10-13 ns, fsl was larger for the β layer in the {111} system than for region B in the {100} system (see Figures 5 and 9). Therefore, if the formation of small cages occurred around all small cage sites, the orientations of THF molecules at small cage sites for the β layer in the {111} system would have a stronger tendency to take orientations in which the molecular planes become parallel to the z axis than those for region B in the {100} system. Nevertheless, the tendency did not appear strongly in the P(θ) for small cage sites of the β layer in the {111} system. Therefore, the results of P(θ)s for the {111} system are attributable to the fact that some THF molecules at small cage sites assumed random orientations owing to the formation of large cages by the surrounding H2O molecules, and that some THF molecules at large cage sites assumed orientations in which the molecular planes became parallel to the z axis owing to the formation of small cages by the surrounding H2O molecules. In this simulation, the modified structure was not stably formed over the whole region of the β layer. Moreover, a transition from the modified structure to the ideal structure was frequently observed. Therefore, the modified structure must have been thermodynamically less stable than the ideal structure. Unfortunately, simulation data was insufficient for evaluating the stability of the modified structure quantitatively. In order to evaluate the stability, it would be required to determine the time scale of the transition from the modified structure to the ideal structure by performing a MD simulation with a much longer run than in this study. Nevertheless, as we discuss in the next section, the transient formation of the modified structure for the {111} system explains why the growth of the clathrate was much slower for the {111} system than for the {100} system.

Nada

Figure 12. Illustrations of a large cage site of the R layer and adjacent large and small cage sites of the β layer. The red and blue circles show positions of large and small cage sites, respectively. Polygons constituting cages around those cage sites are also partly illustrated.

4. Discussion 4.1. Formation of a Modified Structure on the {111} Plane. Here, we discuss why the modified structure was formed only in the β layer for the {111} system. Figure 12 shows schematic illustrations of adjacent large and small cage sites of the β layer and a large cage site of the underlying R layer. Polygons constituting cages around those cage sites are also partly illustrated. The same structure of cage polygons, i.e., a group of three pentagons, appears under both the large and the small cage sites of the β layer. Therefore, when the underlying R layer grows, the large and small cage sites of the β layer are indistinguishable, meaning that if a THF molecule is arranged at a small cage site of the β layer, the surrounding H2O molecules may form a large cage instead of a small cage. If a large cage is formed around the small cage site of the β layer, H2O molecules around the adjacent large cage site of the β layer consequently form a small cage. Thus, the identical structures of cage polygons under large and small cage sites of the β layer explain the formation of the modified structure in the β layer. In addition, we note that if the modified structure was stably formed over the β layer, the clathrate might continue to grow with the formation of a stacking fault. The ideal stacking sequence of {111} layers is...ABCABC.... If each {111} layer is divided into R, β, and γ layers, then each ABC unit is written as A1A2A3B1B2B3C1C2C3, where the subscripts 1, 2, and 3 represent the R, β, and γ layers, respectively. If we assume that the structure of the R layer that grew on the initial position of the clathrate interface in this simulation corresponded to that of A1, then the ideal structure of the β layer that grew subsequently should correspond to that of A2. However, the modified structure corresponded to the structure of B2, but not of A2. Thus, if the modified structure was stably formed over the β layer, the clathrate might continue to grow with the stacking sequence of A1B2B3C1. In this simulation, the growth of a clathrate with a stacking fault was not observed. However, the transient formation of the modified structure in this simulation implies that stacking faults might accidentally form during clathrate growth on {111} planes in real systems. More detailed studies are required to investigate the possibility of the formation of the stacking fault. 4.2. Rate-Determining Process. To explain why the growth velocity of {111} planes of clathrates is the smallest of all crystallographic planes, two hypotheses have been proposed for the rate-determining process of clathrate growth. Smelik and King hypothesized that the formation of small cages is the ratedetermining process, and therefore, the growth of {111} planes, where the concentration of small cages can be higher than those

Growth Kinetics of Tetrahydrofuran Clathrate Hydrate of other planes, is the slowest.14 In contrast, Larsen et al. hypothesized that the formation of hexagons constituting large cages is the rate-determining process, and therefore, the growth of {111} planes, in which the hexagons lie, is the slowest.15 In both of the above two hypotheses, the formation of cages by H2O molecules surrounding cage sites was considered to be the rate-determining process. However, in our simulation, the formation of cages occurred in regions surrounded or sandwiched by THF molecules but did not occur in the regions surrounding cage sites. Moreover, as can be seen in Figure 5, the organization of H2O molecules for the formation of cages progressed more quickly than did the arrangement of THF molecules at cage sites (i.e., fsl drastically increased at a period of 10-15 ns, whereas NgL gradually increased at a period of 4-14 ns). Thus, because our results suggest that the formation of cages by H2O molecules was not the rate-determining process, our results do not agree with the above-mentioned hypotheses. Notably, our results suggest that THF molecules at the clathrate interface rearranging from a disordered state to the ideal state, in which THF molecules are arranged at large cage sites only, is the rate-determining process. In this simulation, the rearrangement occurred much more slowly for the {111} system than for the {100} system, because the formation of the modified structure occurred in the {111} system. This potential rate-determining process explains why the growth of the clathrate was much slower for the {111} system than for the {100} system. Moreover, we considered the structures of other crystallographic planes of the clathrate, and we determined that the formation of the modified structure is not expected to occur on other crystallographic planes. Therefore, the rearrangement of the THF molecules on {111} planes is expected to be slower than the rearrangements on other crystallographic planes of the clathrate. This expectation is consistent with the fact that growth shapes of real THF clathrates; that is, the growth shapes are normally regular octahedral with flat {111} planes, and even if irregular shapes are formed, the growth shapes always have flat {111} planes.14,15 MD simulations for the growth kinetics of other crystallographic planes are required to confirm this expectation. It should be emphasized that this study is the first to suggest a rate-determining process on the basis of results derived from investigating the anisotropy in the molecular-scale growth kinetics of a clathrate. In our simulation, the clathrate was grown from a concentrated THF solution, as in real systems. Therefore, we believe that the rate-determining process suggested by our simulation results applies to clathrate growth in real systems. However, if a clathrate is grown from a dilute THF solution, the rate-determining process might be the transfer of THF molecules, as has been observed in the case of the growth of gas clathrates.24 In the case of structure I clathrates, growth shapes tend to be rhombic dodecahedral with flat {110} planes (the same shapes as common garnet),14,15 meaning that the {110} planes have the slowest growth velocity of all crystallographic planes. We considered whether growth shapes of structure I clathrates can also be explained by the formation of the modified structure, and we determined that the formation of the modified structure is not expected to occur on {110} planes of structure I clathrates. Normally, the growth shape of a crystal reflects the anisotropy in the growth kinetics among all of the crystal’s crystallographic planes. The anisotropic growth kinetics depend on the structure of the crystal. Therefore, it is not at all strange that the anisotropic growth kinetics of structure I clathrates differ from

J. Phys. Chem. B, Vol. 113, No. 14, 2009 4797 those of structure II clathrates. Further studies are required to elucidate the anisotropic growth kinetics of structure I clathrates, as well as the relationship between kinetics and growth shapes. 4.3. Cage Formation Process. The present result that the formation of cages by H2O molecules occurred in regions surrounded or sandwiched by THF molecules was essentially the same as the results obtained in a previous study of the growth of gas clathrates.24 In the case of the growth of gas clathrates, gas molecules act as hydrophobic walls to induce the clustering of confined H2O molecules in spaces surrounded by the gas molecules. In other words, since H2O molecules do not make hydrogen bonds with gas molecules, the confined H2O molecules form an ordered structure themselves so as to fit the spaces. THF molecules interact with H2O molecules more strongly than do gas molecules, because both LJ and Coulomb interactions occur between H2O and THF molecules, whereas only LJ interactions occur between H2O and gas molecules. However, like gas molecules, THF molecules also do not form hydrogen bonds with H2O molecules. Therefore, like gas molecules, THF molecules also have hydrophobic properties rather than hydrophilic properties. Thus, we conclude that the formation mechanism of cages is essentially the same for both THF and gas clathrates. All molecular species that can be included in clathrates as guests have hydrophobic properties rather than hydrophilic properties. Therefore, the formation of cages in all clathrates is expected to essentially involve the clustering of confined H2O molecules in spaces surrounded by guest molecules. In addition, the arrangement of THF molecules at small cage sites might accelerate the formation of empty cages (i.e., small cages of the clathrate). As a result of the arrangement of THF molecules at small cage sites, H2O molecules surrounding small cage sites were partly surrounded or sandwiched by THF molecules, meaning that those surrounded or sandwiched H2O molecules were advantageous for the formation of small cages. Thus, the formation of small cages occurred immediately after the arrangement of THF molecules at small cage sites, and finally those small cages became empty. This process for the formation of empty cages was clearly observed for the growth of {100} layers and also for the growth of the γ layer of a {111} layer. Consequently, the observed arrangement of THF molecules at small cage sites on the clathrate interface may be a key process not only for the formation of the octahedral growth shape of the clathrate but also for the formation of empty cages. 5. Conclusions The molecular-scale growth kinetics of THF clathrates on {100} and {111} planes were investigated by means of MD simulations. The anisotropy in the growth velocity between the systems, which was obtained in the simulation, qualitatively agreed with experimental results in that the growth of the clathrate was much slower for the {111} system than for the {100} system.13-17 For both {100} and {111} systems, when the clathrate grew at the clathrate interface, THF molecules were arranged at both large and small cage sites. Subsequently, the formation of cages by H2O molecules occurred in regions surrounded or sandwiched by THF molecules. As the formation of cages progressed, THF molecules that had once been arranged at small cage sites gradually moved away from the sites. Finally, when THF molecules were ideally arranged at large cage sites only, the structure of the clathrate was completely formed. The arrangement of THF molecules at cage sites prior to the formation of cages by H2O molecules was attributable to the appearance of local minima of a potential energy surface created

4798 J. Phys. Chem. B, Vol. 113, No. 14, 2009 by the underlying clathrate. The appearance of those local minima depends on the structure of the clathrate but is expected to be independent of guest species. Moreover, the observed formation of cages by H2O molecules in regions surrounded or sandwiched by THF molecules was also observed in a previous study for the growth of gas clathrates.24 This phenomenon arises from the hydrophobic properties of clathrate guest molecules; consequently, the observed arrangement of THF molecules at cage sites prior to the formation of H2O cages and the formation of cages in regions surrounded or sandwiched by guest molecules are expected to be similar for the growth of all clathrates. The present results suggest that the rate-determining process for the growth of the clathrate was the rearrangement of THF molecules at the clathrate interface from a disordered state to a state in which THF molecules were ideally arranged at all large cage sites only. This rate-determining process differs from that observed for gas clathrates, in which the transfer of gas molecules is rate-limiting.24 The formation of the modified structure in which large and small cage sites were formed at opposite positions was observed in the {111} system only. Owing to this formation of the modified structure, the rearrangement of THF molecules occurred much more slowly for the {111} system than for the {100} system. The formation of this modified structure is not expected to occur on other crystallographic planes of the clathrates. Therefore, the rearrangement of the THF molecules on {111} planes is expected to be slower than the rearrangements on other crystallographic planes of the clathrate, consistent with the fact that growth shapes of real clathrates are octahedral with flat {111} planes. Notably, the modified structure was formed only when THF molecules were arranged at small cage sites on the clathrate interface. Therefore, the observed arrangement of THF molecules not only at large cage sites but also at small cage sites at the clathrate interface at the beginning of growth is a particularly important process to form octahedral growth shapes with flat {111} planes in real systems. In conclusion, the anisotropy in the molecular-scale growth kinetics of the clathrate, which explains the growth shapes of the clathrates in real systems, was elucidated in this simulation study for the first time. Future studies will involve simulation of the growth of THF-gas binary clathrates. The results obtained by these future studies will be helpful for the capture of gas molecules in clathrates for technological applications. Moreover, simulation of clathrate growth inhibition by kinetic inhibitors, such as polyvinylpyrrolidone38 and antifreeze protein,39 will also be of interest. Growth inhibition is known to occur on {111} planes of clathrates.13,15,40,41 The molecular-scale growth kinetics of THF clathrates on {111} planes, which were elucidated in this study, will be helpful for understanding this growth inhibition mechanism. Acknowledgment. The author thanks Dr. Tsutomu Uchida (Graduate School of Engineering, Hokkaido University) for helpful comments. This work was part of a Joint Research Program of the Institute of Low Temperature Science, Hokkaido University, Japan. Some of the computations in this work were conducted by using the facilities of the Supercomputer Center at the Institute of Solid State Physics, University of Tokyo, Japan, and by using a supercomputer at the Japan Aerospace Exploration Agency (JAXA).

Nada References and Notes (1) Sloan, E. D., Koh, C. A. Clathrate Hydrates of Natural Gases, 3rd ed.; CRC Press: Boca Raton, FL, 2008. (2) Gudmundsson, J. S.; Mork, M.; Graff, O. F. Proc. 4th Int. Conf. Gas Hydrates 2002, 2, 997. (3) Moretti, E. C.; Mukhopadhyay, N. Chem. Eng. Prog. 1993, 89, 20. (4) Mak, T. C. W.; McMullan, R. K. J. Chem. Phys. 1965, 42, 2732. (5) Gough, S. R.; Davidson, D. W. Can. J. Chem. 1971, 49, 2691. (6) Yamamoto, O.; Matsuo, T.; Suga, H.; David, W. I. F.; Ibberson, R. M.; Leadbetter, A. J. Physica B 1995, 213-214, 405. (7) Florusse, L. J.; Peters, C. J.; Schoonman, J.; Heter, K. C.; Koh, C. A.; Dec, S. F.; Marsh, K. N.; Sloan, E. D. Science 2004, 306, 469. (8) Lee, H.; Lee, W.; Kim, J. Y.; Park, D.; Seo, Y.-T.; Zeng, H.; Moudrakovski, I. L.; Ratcliffe, C. I.; Ripmeester, J. A. Nature (London) 2005, 434, 743. (9) Alavi, S.; Ripmeester, J. A.; Klug, D. D. J. Chem. Phys. 2006, 124, 014704. (10) Kim, D.-Y.; Park, Y.; Lee, H. Catal. Today 2007, 120, 257. (11) Talyzin, A. Int. J. Hydrogen Energy 2008, 33, 111. (12) In the case that the gas is hydrogen, H2 molecules are included in both large and small cages of THF clathrates. (13) Makogon, T. Y.; Larsen, R.; Knight, C. A.; Sloan, E. D., Jr. J. Cryst. Growth 1997, 179, 258. (14) Smelik, E. A.; King, H. E., Jr. Am. Mineral. 1997, 82, 88. (15) Larsen, R.; Knight, C. A.; Sloan, E. D., Jr. Fluid Phase Equilib. 1998, 150-151, 353. (16) Iida, T.; Mori, H.; Mochizuki, T.; Mori, Y. H. Chem. Eng. Sci. 2001, 56, 4747. (17) Nagashima, K.; Yamamoto, Y.; Takahashi, M.; Komai, T. Fluid Phase Equilib. 2003, 214, 11. (18) Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: 1996. (19) Ba´ez, L. A.; Clancy, P. Ann. N. Y. Acad. Sci. 1994, 715, 177. (20) Westacott, R. E.; Rodger, P. M. Chem. Phys. Lett. 1996, 262, 47. (21) Rodger, P. M. In Gas Hydrate: Challenges for the Future; Holder, G. D., Bishnoi, P. R., Eds.; Ann. N. Y. Acad. Sci. 2000, 912, 474. (22) Moon, C.; Taylor, P. C.; Rodger, P. M. J. Am. Chem. Soc. 2003, 125, 4706. (23) English, N. J.; MacElroy, M. D. J. Chem. Phys. 2004, 120, 10247. (24) Nada, H. J. Phys. Chem. B 2006, 110, 16524. (25) Hawtin, R. W.; Quigley, D.; Rodger, P. M. Phys. Chem. Chem. Phys. 2008, 10, 4853. (26) Nada, H.; van der Eerden, J. P. J. Chem. Phys. 2003, 118, 7401. (27) Nada, H.; van der Eerden, J. P.; Furukawa, Y. J. Cryst. Growth 2004, 266, 297. (28) Nada, H.; Furukawa, Y. J. Cryst. Growth 2005, 283, 242. (29) Nada, H.; Furukawa, Y. J. Phys. Chem. B 2008, 112, 7111. (30) Carignano, M. A.; Shepson, P. B.; Szleifer, I. Mol. Phys. 2005, 103, 2957. (31) Chandrasekhar, J.; Jorgensen, W. L. J. Chem. Phys. 1982, 77, 5073. (32) Vlot, M. J.; Huinink, J.; van der Eerden, J. P. J. Chem. Phys. 1999, 110, 55. (33) Strictly speaking, the structures of the R and β layers were the same. Both layers contained 12 THF molecules. In this study, 12 empty small cages, which were placed between THF molecular planes of the R and β layers, were defined to be included in the β layer. The R layer was defined as a thin layer containing 12 THF molecules atop the γ layer, and the β layer was defined as a thin layer containing 12 THF molecules and 12 empty small cage sites atop the R layer. According to those definitions, growth of those divided layers was defined to occur in sequence of the R, β and γ layers on the upper clathrate interface in Figure 2c and in sequence of γ, R, and β layers on the lower clathrate interface in Figure 2c. (34) Fincham, D. Mol. Simul. 1992, 8, 165. (35) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (36) Karim, O. A.; Haymet, A. D. J. J. Chem. Phys. 1988, 89, 6889. (37) Jacob, D. M.; Zeidler, M. D.; Kanert, O. J. Phys. Chem. A 1997, 101, 5241. (38) Lederhos, J. P.; Long, J. P.; Sum, A.; Christiansen, R. L.; Sloan, E. D. Chem. Eng. Sci. 1996, 51, 1221. (39) Yeh, Y.; Feeney, R. E. Chem. ReV. 1996, 96, 601. (40) Carver, T. J.; Drew, M. G. B.; Rodger, P. M. J. Chem. Soc., Faraday Trans. 1996, 92, 5029. (41) Zeng, H.; Wilson, L. D.; Walker, V. K.; Ripmeester, J. A. J. Am. Chem. Soc. 2006, 128, 2844.

JP810041T