Article pubs.acs.org/JPCC
Voronoi Tessellation Analysis of Clathrate Hydrates Somendra N. Chakraborty,† Eric M. Grzelak,† Brian C. Barnes,† David T. Wu,*,†,‡ and Amadeu K. Sum*,† †
Center for Hydrate Research, Chemical & Biological Engineering Department and ‡Chemistry Department, Colorado School of Mines, Golden, Colorado, United States S Supporting Information *
ABSTRACT: Molecular simulation of clathrate hydrate has provided significant advancements in our understanding of hydrate properties and formation. In this work, we report the application of Voronoi tessellation to characterize the structuring of water and guest molecules forming hydrates. Tessellation of perfect sI and sII hydrate reveals positions of Voronoi vertices similar to the oxygen atoms of enclathrating water molecules. Applying tessellation to a simulation trajectory of hydrate formation, and using a further selection criteria based on polyhedra volume and coordination number, we identify numbers and types of cagelike polyhedra. Voronoi analysis of this type results in similar numbers of identified cages but with differing topologies. However, once nearest neighbor methanes are also enclathrated, the topologies of the Voronoi polyhedra approach that of the actual water cages. Since only methane coordinates are required, Voronoi tessellation is a fast and simple tool that can be used as an order parameter to identify the structuring of molecules when studying hydrates in simulations.
■
INTRODUCTION Clathrate hydrates are crystalline compounds of water and guest molecules that are usually formed under conditions of high pressure and low temperature. The water molecules in the hydrate framework are connected by hydrogen bonds, similar to those in ice, forming cages that enclose the guest molecules, such as methane, ethane, carbon dioxide, hydrogen, or other small molecules. Hydrates are commonly found in one of two crystal structures depending on the guest molecule being enclathrated. Structure I (sI) hydrate is a cubic crystal (a = 12.1 Å) with the unit cell composed of 46 water molecules forming eight cages: two pentagonal dodecahedral (512) cages, formed by 20 water molecules, and six tetrakaidecahedral (51262) cages, formed by 24 water molecules.1 Structure II (sII) is also cubic (a = 17.0 Å) and has 136 water molecules in the unit cell. It is composed of 16 512 cages and 8 hexakaidecahedral (51264) cages, formed by 28 water molecules.1 Molecular simulations of hydrate formation have been a challenging and active area of research. Molecular level understanding of hydrate formation is a fundamentally complex and important area that has implications for every area involving hydrates, including the formation and prevention of hydrates in oil/gas pipelines, the formation and dissociation of hydrates in nature, and the storage and transportation of gases. To date, there have been many computational studies aimed at understanding hydrate formation.2−22 The formation of hydrate from simulations, in general, has been described by order parameters (OP),7,23−25 which are a quantitative measure of the degree of order in the system. Identifying an appropriate OP is often a difficult task as it should capture the important events before, during, and after the formation of hydrates. For hydrates, the task of identifying a © 2012 American Chemical Society
good OP is perhaps even more challenging since formation might proceed by multiple steps and might involve cooperative participation of both water and guest molecules. Often, combinations of multiple order parameters are used together to better characterize the development of order. There have been different OPs used in the study of hydrate formation to characterize the order of water and guest molecules. One class of these OPs is based on the geometric arrangement of water molecules, such as F3 (measure of deviation from tetrahedrality in water molecules),6 F4 (measure of torsion angle between hydrogen bonded water molecules),26 and Steinhardt OP (measure of orientation of bonds).7 While useful, these OPs are typically applied globally and only capture the geometry of a few neighboring atoms if applied locally and may not be well-suited to isolate cage-scale events in the early stages of hydrate formation. Consequently, the other type of commonly used OP is based on the water cages formed during hydrate formation. Appearance of stable cages is the most obvious observable of hydrate formation.3,19,15 In our initial studies of the hydrate formation process,3−5 it is observed that when cages start appearing, nucleation initiates and thereafter the solid structure spontaneously grows. In this process, both the water and guest molecules undergo significant rearrangement to form specific cages. The nature of the methane−water cooperative ordering that leads to cages and growth thereafter is a challenge to quantify. Here, we examine a different approach that could be applied in the study of hydrate formation. We tessellate space on the Received: May 11, 2012 Revised: August 22, 2012 Published: August 24, 2012 20040
dx.doi.org/10.1021/jp304612f | J. Phys. Chem. C 2012, 116, 20040−20046
The Journal of Physical Chemistry C
Article
Figure 1. Voronoi tessellation of structure I and II clathrate hydrates. The blue, green, and yellow polyhedra are representative of the 512, 51262, and 51264 water cages, respectively. Shown are the polyhedra resulting from the tessellation of perfect (methane located at the center of the water cages at idealized lattice locations) and equilibrated (methane displaced from the center of the water cages because of thermal fluctuations) crystals.
basis of the positions of the methane molecules. Voronoi tessellation generates polyhedra, whose topological and metric properties can provide a quantification of local ordering and can thus serve as an order parameter. The goal of this study is to present the Voronoi polyhedra (VP) analysis method and to demonstrate how it may be used to characterize the ordering and structuring in clathrate hydrates. It is also important to note that our analysis is based on the coordinates of methane molecules only, which makes it much faster to perform than other analysis methods involving water molecules, such as all the current order parameters described above. As such, the VP order parameters could be useful alone or complementing existing order parameters.
methane molecules. We also compare the results of the VP analysis to the seven types of water cages which form during hydrate nucleation.4 Hereafter, the term cage will refer to such cages identified on the basis of the water molecules, whereas the term polyhedron will refer to the polyhedra obtained from the VP analysis on the basis of the methane molecules. The VP analysis is performed using the open source Voro++ software.40 For the identification of the water cages, we use an independent code reported in detail elsewhere;4 the code was developed to identify seven types of hydrate cages on the basis of the positions and orientations of the water molecules.
METHODOLOGY The VP of a given point in a set of points is the locus of all points that are closer to this point than to any other point in the set.27,28 Voronoi tessellation has been applied in many other areas to study the local environment of ordering systems, including water, Lennard-Jones solids, glass-forming liquids, colloids, and hydrated phospholipid membranes.29−37 Important metric properties of the VP include the polyhedron volume, face areas, and edge lengths. Topological properties of the polyhedra include the numbers of faces, edges, and vertices. In this study, we monitor such properties from the VP analysis on the basis only of the methane molecules to gain insight into the mechanism of hydrate formation and structuring. As a reference system, we first collect the topological and metric properties of the VP analysis for perfect sI and sII hydrate lattices. The Voronoi tessellation is performed for methane molecules placed at the center of each cage. A pure methane sII hydrate is not observed in experiment, but this idealized structure is considered here to illustrate that the tessellation can treat and identify structures on the basis only of the coordinates of the guest molecules. A comparison of the polyhedra generated from the VP analysis is then made with the cages actually formed by the water molecules in sI and sII hydrates. The coordinates of the oxygen atoms in the water molecules in the unit cell for the perfect sI and sII were those determined by X-ray crystallography.38,39 We next apply the VP analysis to a molecular dynamics trajectory of water and methane that we simulated previously for studying the nucleation and growth of hydrates (details described elsewhere4). This particular nucleating trajectory is from a canonical ensemble simulation at 245 K density-equilibrated beforehand at 2000 bar with a total simulation time of 790 ns. Nucleation (defined operationally here as the appearance of the first cage) occurred at about 90 ns from the start of the simulation. The system contained 11 776 water and 2048
VP Analysis of Perfect and Equilibrated sI and sII Crystals. Figure 1 shows the tessellated sI and sII hydrate crystals. It is clear from the figure that the shapes of the polyhedra closely resemble that of the actual cages found in the hydrate structures as expected since the Voronoi vertices are the Delaunay dual of the methane positions. For clarity, we show a portion of sI and sII tessellated crystal that has distinct cages/polyhedra. Here, we show the 512 and 51262 polyhedra of the sI hydrate and the 512 and 51264 polyhedra of the sII hydrate. The naming of polyhedra follows the convention used for cages specifying the number and type of faces present. Thus, the 512 polyhedra comprise 12 pentagonal faces, the 51262 polyhedra comprise 12 pentagonal faces and 2 hexagonal faces, and the 51264 polyhedra comprise 12 pentagonal and 4 hexagonal faces. Figure 2 superimposes images of the individual polyhedra and cage types present in the perfect sI and sII hydrate lattice. The vertices of the polyhedra closely match the positions of the oxygen atoms of the water molecules making up the cages in each case. The VP analysis identifies the same number and
■
■
RESULTS AND DISCUSSION
Figure 2. Overlap of polyhedra and water cages found in structure I and II hydrates. The oxygen atom of the water molecules and the outline of the water cages are shown in red. The polyhedra are equivalent to the 512 (blue), 51262 (green), and 51264 (orange) cages. For clarity, hydrogen atoms of the water molecules are not shown. Polyhedra are similarly colored as in Figure 1. 20041
dx.doi.org/10.1021/jp304612f | J. Phys. Chem. C 2012, 116, 20040−20046
The Journal of Physical Chemistry C
Article
Figure 3 shows the distribution of Voronoi polyhedra volumes at different time steps of the nucleating methane
types of polyhedra as cages found in the sI and sII hydrate unit cells, that is, 2 512 and 6 51262 in sI and 16 512 and 8 51264 in sII. Table 1 reports the topological and metric properties of the three different types of polyhedra found in the perfect sI and sII Table 1. Properties of Voronoi Polyhedra Identified in Perfect Structure I and II Hydratesa polyhedron 12
5 (sI) 512 (sII) 51262 51264
nedges
nvertices
nfaces
S (Å2)
V (Å3)
η
30 30 36 42
20 20 24 28
12 12 14 16
184.82 183.34 187.37 203.91
204.52 200.39 211.06 246.44
1.335 1.338 1.306 1.235
S is the surface area; V is the volume. Asphericity parameter η = S3/ (36πV2); η = 1 for a sphere.
a
crystals. The number of faces, vertices, and edges of each type of polyhedron matches those of the actual cages in the hydrate structures. Furthermore, the polyhedral volume and surface area increase in the order 512 < 51262 < 51264. The asphericity parameter, which can also be used to characterize the type of cage/polyhedron, ranges from 1.2 to 1.4. Table 1 shows that these properties of the polyhedra are equivalent to the cages and can thus be used to identify hydrate cages and structures. We tested the sensitivity of the polyhedral topology in sI and sII to random displacements of methane molecules in the cages in two different ways. First, we performed VP analysis on a 100 ps simulated trajectory of equilibrated sI and sII crystals at 300 K. Figure 1 also shows the shapes of these polyhedra at the last time step of the simulation. We observe here that for both sI and sII, the topologies of the polyhedra are preserved. Second, we performed VP analysis after randomly displacing methane molecules from their lattice positions inside the cages (see Figure S1 in the Supporting Information). We find that for maximum displacements greater than 2.0 Å from the center of each cage irregular polyhedra were generated containing three-, four-, five-, six-, and seven-sided faces. VP Analysis of a Nucleating Trajectory. Having demonstrated that the VP analysis gives polyhedra whose properties match closely with that of the actual cages in the sI and sII crystals, we then applied it to the analysis of a hydrate nucleation trajectory. In this trajectory, initially metastable phase separated methane and water molecules start intermixing and eventually nucleate to form hydrates with the first hydrate cage appearing around 90 ns. Tessellation of the methane molecules at each configuration separated by 1 ns is then performed to determine the polyhedra containing each methane molecule over the course of the nucleation. Because we are only interested in the polyhedra for the methane molecules that experience cagelike environments in the nucleation and growth process, we applied two criteria to select for the relevant cagelike polyhedra: (1) a volume range (between 170 and 300 Å3) and (2) a coordination number within a cutoff distance (minimum of 12 methane molecules within 10 Å). In the second criterion, the choice of 12 methane molecules within the cutoff is based on the expectation of identifying the smallest cages in sI and sII structures. The possibility of a cagelike polyhedron from fewer methane molecules within this cutoff is of course possible, however, the polyhedron generated may not result in a stable structure/cage. Moreover, the methane−methane radial distribution function in water shows that very few methane molecules are found beyond 10 Å because of the low solubility of methane in water.
Figure 3. Distribution of Voronoi polyhedra volume at different times in a molecular dynamics (MD) simulation trajectory for methane + water. Hydrate nucleates at about 90 ns. The inset shows the total distribution. Lines correspond to 0 (black), 200 (red), 400 (green), and 600 (blue) ns.
trajectory. The sharp peaks around 65 Å3 correspond to the volume of the VP of methane molecules in the gas phase. The smaller peak (blue) near 220 Å3 corresponds to the polyhedral volume obtained from tessellating methane molecules in clathrate hydrate cagelike environments. These peaks are expanded in the larger plot (total distribution is shown in the inset). With time, the peaks around 65 Å3 decrease in intensity and those around 220 Å3 increase. Along with these peaks, the distribution also has a decaying tail from 300 to 800 Å3. We notice that this tail of the volume distribution is associated with the polyhedra of methane molecules that are surrounded by water molecules and methane molecules that are beyond 10 Å. We also analyzed five additional nucleating trajectories with similar simulation conditions which had previously been reported by our group (two canonical ensemble and three isobaric−isothermal ensemble trajectories). The results for quantities such as the distribution of Voronoi polyhedra volumes in those trajectories varied with the amount of clathrate formed in the system similarly to the results presented here. As more methane molecules were enclathrated, the total volume of the Voronoi polyhedra in the system increased and at volume ranges associated with polyhedra seen in sI and sII clathrate structures (roughly 180−260 Å3) supporting the generality of our findings. Using only the first criterion on volume to select the relevant methane molecules, we obtain a large number of polyhedra (see top left plot in Figure 4), which include all polyhedra irrespective of the topology, with volume in a range close to that of the cages of interest (see Table 1). The volume range used is larger than the known volumes of the perfect cage/ polyhedra in order to account for fluctuations. This volume criterion acts to select polyhedra for methane molecules that are solvated (dissolved in bulk water) and excludes those in the bulk methane phase. The number of such polyhedra identified is much larger than the number of actual water cages, which is expected since we are considering any possible topology, 20042
dx.doi.org/10.1021/jp304612f | J. Phys. Chem. C 2012, 116, 20040−20046
The Journal of Physical Chemistry C
Article
polyhedron corresponding to a 512 cage, methane molecules must be suitably surrounding the central methane molecule (with a local environment similar to the 512 cage in a crystal). Interestingly, this happens most often when the surrounding methane molecules are themselves also in cagelike polyhedra, that is, methane molecules which are deep inside a growing hydrate structure are more likely to have a polyhedron that is cagelike. The topological information and discrete counts of different stuctures provided by Voronoi polyhedra analysis may be contrasted with, for example, carbon−carbon radial distribution functions as shown in Figure 5. In that figure, the enclathration
Figure 4. Number of Voronoi polyhedra (medium brown bars) and water cages (dark maroon bars) identified at different times during an MD simulation trajectory of hydrate forming from methane + water. The light gray bars in the first subplot correspond to the number of polyhedra when only a volume criterion is applied. All the medium brown bars correspond to the number of polyhedra when both volume and coordination number criteria are applied. Figure 5. Carbon−carbon radial distribution functions (RDFs) for an MD simulation trajectory of hydrate forming from methane + water. The simulation duration was 790 ns, and the eight RDFs were generated from data in sequential 100 ns bins along the trajectory (90 ns for the final bin). The RDF formed from data in the first 100 ns is represented by a red line and in the final 90 ns is represented by a thick black line. There is a smooth progression in changes to the RDFs from start to finish.
including those that may be noncagelike. However, since we are mostly interested in the polyhedra for methane molecules in cagelike environments, we also filter the polyhedra on the basis of the second criterion. The second criterion significantly restricts the polyhedra identified to those that more closely resemble the water cages. As shown in Figure 4, the total number of polyhedra obtained in this way mirror those of the actual water cages. This second criterion also eliminated polyhedra for methane molecules near the methane−water interface. The choice of at least 12 methane molecules is made because this corresponds to the minimum coordination number of a methane molecule in a hydrate cage, specifically, the 512 cage. Figure 4 also shows the number of cages and polyhedra identified for seven important cage types4 at different time steps of the hydrate simulation trajectory (the topology may be the same but the polyhedral vertices may not necessarily align with the water positions of the cage). The number of 512 polyhedra is the highest among all the polyhedra. However, the actual number of 512 polyhedra is smaller than for the 512 cages (in the perfect sI and sII hydrates, the number and type of polyhedra exactly matched the cages). This difference in number arises because not all 512 cages are identified as 512 polyhedra but instead as 3m4n5o6p7q (with possibly nonzero m, n, o, p, q) and 4n5o6p (n > 1, o < 10, and p > 1) polyhedra. Likewise, 51262, 51263, and 51264 cages can sometimes correspond to different polyhedra with topology having four-sided faces. As a result, the total number for each of these polyhedra is always lower in number than the corresponding cage counterpart. This is also part of the reason why the numbers of the 4151062, 4151063, and 4151064 polyhedra are higher than those of the corresponding cage (these polyhedra with four-sided faces are identified in place of the 5126p polyhedra). We also find that to observe a 512
of gas is inferred by the decrease in size of the 4 Å nearestneighbor peak for methanes in the gas phase and the gradual change of peak size and location in the 6.5−7.5 Å range. The methane−water slab interface no longer exists at the end of this trajectory, but a smaller cylindrically shaped phase of methane gas is present. The presence of both bulk methane gas and amorphously structured clathrate hydrate makes determination of structure from pair correlation functions difficult. The C−C radial distribution function is most useful as an indicator of structure when the final system is in a nearly perfect crystal. However, the Voronoi analysis is quite useful for locally quantifying the topology of domains within a system showing which type of polyhedra are being formed as the hydrate growth proceeds. That type of information is not easily available from the analysis of pair correlation functions. This supports the usefulness of Voronoi tessellation analysis in systems which are not entirely well-ordered. When visualized, the Voronoi tessellation directly and efficiently provides anisotropic information that is not readily extracted from isotropic radial distribution analysis and other common order parameters. Results for coordination number as determined by integration of radial distribution functions and other common order parameters are also presented in the Supporting Information (Figures S2−S4). While order parameters such as the q6 or F4 parameters may be useful for characterizing the 20043
dx.doi.org/10.1021/jp304612f | J. Phys. Chem. C 2012, 116, 20040−20046
The Journal of Physical Chemistry C
Article
general state of a system,3,5 they do not offer the richness of geometric or topological information provided by Voronoi polyhedra analysis as shown in Figures 3 and 4. A more detailed comparison was made between the topology of cages and polyhedra for individual methane molecules at a specific time. For the cage and polyhedron comparison, we selected the methane molecules occupying a particular cage type and identified the corresponding polyhedron topology (see Table 2). Our analysis suggests that almost 50% of the 512
local ordering of methane molecules clarifying the scope of the approach for one-to-one comparisons of cages and polyhedra. To further illustrate the usefulness of the VP analysis, we capture the evolution of the Voronoi polyhedron of a methane molecule that ultimately becomes enclathrated by a water cage. The coordination number and the polyhedron properties of a particular methane are shown in Figure 6. In the top left plot,
Table 2. Cages to Polyhedra Comparison for Individual Methane Molecules at 600 nsa 512 51262 51263 51264 4151062 4151063 4151064
cages
polyhedra
511
510
59
58
57
56
55
54
107 78 21 6 21 6 6
51 24 5 0 5 0 0
1 2 0 1 0 0 0
11 18 0 3 0 0 0
9 8 2 0 3 0 1
12 12 4 2 5 2 1
6 5 2 0 1 0 2
7 5 3 0 1 1 1
1 1 4 0 0 1 0
6 0 0 0 2 1 0
a
The data in the table show the counts of different polyhedra for each cage type. The second and third columns compare the numbers between same types of cages and polyhedra. As an example, the information in the table should read as follows: 107 methane molecules are enclathrated in 512 cages of which 51 are 512 polyhedra and the remaining polyhedra are 3m4n511−46p7q types (see the Supporting Information for a description of the polyhedra types).
Figure 6. Voronoi polyhedron properties and coordination number for an individual methane molecule. The evolution of the Voronoi polyhedron around this methane molecule is depicted in Figure 7.
cages are identified as 512 polyhedra, and this percentage decreases for 51262 and becomes near zero for the other cage types. The polyhedra with topology differing from the corresponding cage had various topologies represented as 3m4n5o6p7q. Similarly to the polyhedron and cage comparison, we selected the methane molecules with a particular polyhedron topology and identified the actual cage type (see Table 3). Here, we see that for the most common types of polyhedra (512, 51262, 51263, 51264), the VP analysis correctly identified the cages more than ∼80% of the time. As also seen from this comparison, a large number of unusual polyhedra are often identified in place of the actual cages. These observations demonstrate the sensitivity of the polyhedra topology to the
the methane−methane (red) and methane−water (black) coordination numbers within a 10 Å and 5 Å spherical cutoff stabilize to 12 and 20, respectively. These spherical cutoffs are chosen to correspond to the first minimum in the methane− methane (in water) and methane−water radial distribution functions. In the top right plot, the number of four, five, and six faces are shown in black, red, and blue, respectively. The methane environments identified by both of these properties (coordination number and topology of the polyhedron) are similar as they stabilize concurrently. Moreover, the fluctuations
Table 3. Comparison of Polyhedra to Cages at 600 nsa 512 51262 51263 51264 4151062 4151063 4151064 425861 425862 425863 425864 3m4n5o6p7q
polyhedra
512
51262
51263
51264
4151062
4151063
4151064
67 30 7 0 44 15 6 11 18 10 13 76
51 0 0 0 9 2 0 3 3 0 0 25
0 24 2 0 13 0 1 0 4 2 1 20
0 1 5 0 0 0 0 0 0 0 1 6
0 0 0 0 1 2 0 0 0 0 0 1
3 0 0 0 5 0 1 1 0 0 1 5
0 0 0 0 1 0 0 0 0 0 0 4
0 0 0 0 1 0 0 0 0 0 1 2
The first and second columns correspond to the polyhedra identified and the corresponding number, respectively. The remaining columns show the number of actual cages matched for each polyhedron type. The last row has entries for triangular-faced polyhedra with enclathrated methane molecules and their cage count.
a
20044
dx.doi.org/10.1021/jp304612f | J. Phys. Chem. C 2012, 116, 20040−20046
The Journal of Physical Chemistry C
Article
Figure 7. Snapshot of the structuring around the same (topology and coordination numbers previously plotted in Figure 6) methane molecule (green) surrounded by water molecules (red) and Voronoi vertices (joined by blue cylinders) at progressing times in a methane + water MD simulation trajectory. The topology and volume of the formed polyhedron are also shown. Methane molecules in the vicinity are not shown. For clarity, hydrogen atoms in the water molecules are not shown.
with time. This implies that with time, dissolved methane molecules cluster with water molecules with Voronoi polyhedra similar to the water cages. In some cases, the topology of the polyhedra perfectly matches that of its corresponding cage. From this study, we can conclude that Voronoi polyhedra analysis provides useful and convenient order parameters to describe the cooperative ordering of methane and water molecules in the formation of clathrate hydrates. One of the major advantages of the VP analysis is that it only requires the position of the methane molecules, which consequently results in a much faster analysis than cage counting. Cage-counting algorithms typically loop multiple times over pairs of water molecules to obtain the edges, polygons, then polyhedra, whereas three-dimensional Voronoi tessellation algorithms scan progressively through the guest positions running in O(N ln N). The number of methane molecules in hydrates is ∼1/6 the number of water molecules, which also helps give tessellation an advantage in speed. The VP analysis was found to be about 2 orders of magnitude faster than an existing cage count code.4 While the polyhedron type obtained does not always match with the cage type, the tessellated polyhedra do capture the cooperative ordering of methane and water molecules. One of the shortcomings of the VP analysis for identification of cage types in clathrate hydrates is the necessity for the methane molecules to cluster sufficiently to provide a local environment whose Voronoi polyhedra are equivalent to the cages. This condition also excludes identification of empty or doubly occupied cages from the VP analysis. Because the VP analysis does not always reliably generate the cages, its main strength as an order parameter may be to efficiently screen for the ordered environment between water and methane molecules, so that other local structure analyses can be applied efficiently. Thus, Voronoi polyhedra analysis may contribute to a more complete and detailed mechanistic understanding of hydrate structure formation.
in all these properties stabilize to specific values indicating that the methane molecule is in a polyhedron equivalent to a hydrate cage, which for this methane is a 512 polyhedron. Figure 6 also explains the journey of this methane from the gas phase to enclathration. The methane is initially in the gas phase (small polyhedron volume), becomes solvated (large polyhedron volume), and finally becomes enclathrated (polyhedron volume similar to a cage). This process is displayed visually in Figure 7, which shows a sequence of snapshots of the polyhedron and nearest water molecules (within 5 Å) to the selected methane molecule. The methane molecule (light green sphere) is initially surrounded by a few water molecules (red spheres). To make the polyhedron visible, we connect the vertices with thick blue lines and color the faces in blue. In this example, at zero ns, the first solvation shell is shared by both water and methane molecules as the central methane is located at the interface of bulk water and bulk methane phases resulting in an irregular polyhedron. At 94 ns, the central methane molecule is only solvated by water molecules (no other methane molecules are found in the first solvation shell), and thus the polyhedron generated is large in volume. At 216 and 282 ns, a polyhedron that more closely resembles a cagelike environment appears, and by 302 ns, the locations of the polyhedron vertices and the water molecules nearly overlap. At this configuration, this particular methane molecule is surrounded by other methane molecules that are suitably located in the near vicinity; therefore, the tessellation leads to a 5 12 polyhedron corresponding to a 512 cage.
■
CONCLUSION By performing three-dimensional Voronoi tessellation of methane molecules, we obtain useful information about the ordering of methane and water that is closely related to the water cages formed in clathrate hydrates. For perfect sI and sII hydrates, tessellation using only the positions of the methane molecules generates Voronoi polyhedra that closely match the actual water cages. When we extend this analysis to investigate a hydrate formation trajectory, we note that a one-to-one correspondence between the number and type of polyhedra and cages is not necessarily obtained. Such a one-to-one correspondence is recovered when the surrounding (12 or more) methane molecules are in place around a particular methane molecule. One consequence is that the VP analysis is sensitive to structure formation after several water cages have formed. For a perfect crystal, however, all the cages are occupied, resulting in Voronoi polyhedra that resemble the cages. In monitoring the polyhedra formed in the MD trajectory, we do see an increase in the 512 and other polyhedra
■
ASSOCIATED CONTENT
S Supporting Information *
VP properties of all the methanes of sI and sII crystal, polyhedra detected from VP analysis of the nucleating trajectory at different time steps, plots of several order parameters for the nucleating trajectory and a specific methane within the trajectory, plot of selected VP properties (volume and number of faces) after tessellation of randomly displaced methane molecules from the center of the cages. This material is available free of charge via the Internet at http://pubs.acs. org/. 20045
dx.doi.org/10.1021/jp304612f | J. Phys. Chem. C 2012, 116, 20040−20046
The Journal of Physical Chemistry C
■
Article
(28) Finney, J. L. Proc. R. Soc. Lond. 1970, 319, 479−493. (29) Swope, W. C.; Andersen, H. C. Phys. Rev. B. 1990, 41, 7042− 7054. (30) Ruocco, G.; Sampoli, M.; Valluri, R. J. Chem. Phys. 1992, 96, 6167−6176. (31) Idrissi, A.; Vyalov, I.; Kiselev, M.; Fedorov, M.; Jedlovszky, P. J. Phys. Chem. B 2011, 115, 9646−9652. (32) Starr, F. W.; Sastry, S.; Douglas, J. F.; Glotzer, S. C. Phys. Rev. Lett. 2002, 89, 125501. (33) Alinchenko, M. G.; Voloshin, V. P.; Medvedev, N. N.; Mezei, M.; Partay, L.; Jedlovszky, P. J. Phys. Chem. B. 2005, 109, 16490− 16502. (34) Jedlovsky, P. J. Chem. Phys. 2000, 113, 9113−9121. (35) Ruocco, G.; Sampoli, M.; Torcinia, A.; Vallauri, R. J. Chem. Phys. 1993, 99, 8095−8104. (36) Kondo, T.; Tsumuraya, K. J. Chem. Phys. 1991, 94, 8220−8226. (37) Chakraborty, D.; Chandra, A. J. Mol. Liq. 2011, 163, 1−6. (38) Mak, T. C. W.; McMullan, R. K. J. Chem. Phys. 1965, 42, 2732− 2737. (39) McMullan, R. K.; Jeffrey, G. A. J. Chem. Phys. 1965, 42, 2725− 2732. (40) Rycroft, C. H. Multiscale Modeling in Granular Flow. Ph.D. Thesis, Massachusetts Institute of Technology, 2007.
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected] (D.T.W.);
[email protected] (A.K.S.). Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS We acknowledge funding from the U.S. National Science Foundation (CBET-0933856 and CHE-1125235). The simulations were performed on HPC facilities at the Golden Energy Computing Organization at the Colorado School of Mines using resources acquired with financial assistance from the NSF and the National Renewable Energy Laboratory.
■
REFERENCES
(1) Sloan, E. D.; Koh, C. A. Clathrate Hydrates of Natural Gases, 3rd ed.; CRC Press: Boca Raton, FL, 2008; pp 45−111. (2) Walsh, M. R. Methane Hydrate Nucleation Rates and Mechanisms from MD Simulations. Ph.D. thesis, Colorado School of Mines, 2011. (3) Walsh, M. R.; Koh, C. A.; Sloan, E. D.; Sum, A. K.; Wu, D. T. Science 2009, 326, 1095−1098. (4) Walsh, M. R.; Rainey, J. D.; Lafond, P. G.; Park, D.-H.; Beckham, G. T.; Jones, M. D.; Lee, K.-H.; Koh, C. A.; Sloan, E. D.; Wu, D. T.; et al. Phys. Chem. Chem. Phys. 2011, 13, 19951−19959. (5) Walsh, M. R.; Beckham, G. T.; Koh, C. A.; Sloan, E. D.; Wu, D. T.; Sum, A. K. J. Phys. Chem. C 2011, 115, 21241−21248. (6) Baez, L. A.; Clancy, P. Ann. N.Y. Acad. Sci. 1994, 715, 177−186. (7) Radhakrishnan, R.; Trout, B. L. J. Chem. Phys. 2002, 117, 1786− 1796. (8) Moon, C.; Taylor, P. C.; Rodger, P. M. J. Am. Chem. Soc. 2003, 125, 4706−4707. (9) Moon, C.; Hawtin, R. W.; Rodger, P. M. Faraday Discuss. 2007, 136, 367−382. (10) Guo, G. J.; Zhang, Y. G.; Liu, H. J. Phys. Chem. C 2007, 111, 2595−2606. (11) Hawtin, R. W.; Quigley, D.; Rodger, P. M. Phys. Chem. Chem. Phys. 2008, 10, 4853−4864. (12) Zhang, J.; Hawtin, R. W.; Yang, Y.; Nakagava, E.; Rivero, M.; Choi, S. K.; Rodger, P. M. J. Phys. Chem. B 2008, 112, 10608−10618. (13) Guo, G. J.; Li, M.; Zhang, Y. G.; Wu, C. H. Phys. Chem. Chem. Phys. 2009, 11, 10427−10437. (14) Conde, M. M.; Vega, C. J. Chem. Phys. 2010, 133, 064507. (15) Jacobson, L. C.; Molinero, V. J. Phys. Chem. B 2010, 114, 7302− 7311. (16) Jacobson, L. C.; Hujo, W.; Molinero, V. J. Am. Chem. Soc. 2010, 132, 11806−11811. (17) Liang, S.; Kusalik, P. G. Chem. Phys. Lett. 2010, 494, 123−133. (18) Tung, Y. T.; Chen, L. J.; Chen, Y. P.; Lin, S. T. J. Phys. Chem. B 2010, 114, 10804−10813. (19) Vatamanu, J.; Kusalik, P. G. Phys. Chem. Chem. Phys. 2010, 12, 15065−15072. (20) Jensen, L.; Thomsen, K.; von Solms, N.; Wierzchowski, S.; Walsh, M. R.; Koh, C. A.; Sloan, E. D.; Wu, D. T.; Sum, A. K. J. Phys. Chem. B 2010, 114, 5775−5782. (21) Bai, D.; Chen, G.; Zhang, X.; Wang, W. Langmuir 2011, 27, 5961−5967. (22) Lyakhov, G. A.; Mazo, D. M. Phys. Wave Phenom. 2009, 17, 289−300. (23) van Duijneveldt, J. S.; Frenkel, D. J. Chem. Phys. 1992, 96, 4655−4668. (24) Trudu, F.; Donadio, D.; Parrinello, M. Phys. Rev. Lett. 2006, 97, 105701. (25) Quigley, D.; Rodger, P. M. J. Chem. Phys. 2008, 128, 154518. (26) Rodger, P. M.; Smith, W.; Forester, T. R. Fluid Phase Equilib. 1996, 116, 326−332. (27) Voronoi, G. J. Reine Angew. Math. 1908, 133, 97−178. 20046
dx.doi.org/10.1021/jp304612f | J. Phys. Chem. C 2012, 116, 20040−20046
