Microsecond Molecular Dynamics Simulations of the Kinetic

ARTICLE pubs.acs.org/Langmuir

Microsecond Molecular Dynamics Simulations of the Kinetic Pathways of Gas Hydrate Formation from Solid Surfaces Dongsheng Bai,† Guangjin Chen,*,‡ Xianren Zhang,*,† and Wenchuan Wang† †

Division of Molecular and Materials Simulation, Key Lab for Nanomaterials, Ministry of Education, Beijing University of Chemical Technology, Beijing 100029, China ‡ High Pressure Fluid Phase Behavior & Property Research Laboratory, School of Chemical Engineering, University of Petroleum, Beijing 102200, China ABSTRACT: In this paper, we report microsecond molecular dynamics simulations of the kinetic pathway of CO2 hydrate formation triggered by hydroxylated silica surfaces. Our simulation results show that the nucleation of the CO2 hydrate is a three-stage process. First, an icelike layer is formed closest to the substrates on the nanosecond scale. Then, on the submicrosecond timescale, a thin layer with intermediate structure is induced to compensate for the structure mismatch between the icelike layer and the final stable CO2 hydrate. Finally, on the microsecond timescale, the nucleation of the first CO2 hydrate motif layer is generated from the intermediate structure that acts as nucleation seeds. We also address the effects of the distance between two surfaces.

1. INTRODUCTION Gas hydrates are clathrate hydrates composed of gas and liquid water molecules1 in which water molecules form cagelike cavities to encapsulate gas molecules. In the unit cell of the sI hydrate, there are typically 512 cages (small) and 51262 cages (large). A 512 cage is formed by 12 pentagonal rings, and a 51262 cage is formed from 12 pentagonal rings and 2 hexagonal rings. Recently, gas hydrates have been extensively studied because they would be essential for future energy resource and CO2 storage. As a new energy resource, natural gas (mainly CH4) hydrates are a geologically important class of materials and are abundant in the arctic region and marine sediment. For hydrates in nature, the hydrate formation zone is an area of sedimentary rock or unconsolidated clay saturated with water and gas in which the rocks and clays consist of grains of different sizes and shapes. Different from bulk systems, the presence of solid surfaces or rock pores may play an important role in the formation of hydrates. The interaction between gas hydrates and the surrounding sediment grains determines the basic physical properties of the sediment hydrate matrix.1,2 Nevertheless, a detailed knowledge of the interaction on a microscopic level is still not clear. For example, the formation of hydrates was found to consolidate rock grains, whereas the mechanism is not well known because of the poorly understood sediment hydrate interaction. Moreover, many techniques for gas hydrate detection and qualification, which are based on the acoustic properties r 2011 American Chemical Society

of the gas hydrate-bearing sediments, are highly dependent on the sediment hydrate interaction. Cha et al.3 observed that the formation rate of the CH4 hydrate is greatly enhanced by the presence of a bentonite surface. They argued that the solid surface provides a nucleation site to which the hydrate crystal can bind during its initial formation stage. Riestenberg et al.4 indicated that the presence of sediments alters the equilibrium boundary for phase transformation. One of our authors and his co-workers also demonstrated that the formation of a CH4 hydrate could be enhanced by immersing activated carbon in water.5 In general, variations in macroscopic behavior, as in this experimental research, demonstrate that the appearance of the solid surfaces could strongly affect the formation of a hydrate. This is because in the range of fluid rock interactions the local structuring of water molecules and molecule adsorption induced by solid surfaces may alter the processes of the nucleation and growth of the hydrate. However, the pathways and mechanisms (i.e., how solid surfaces trigger the nucleation of gas hydrates from a metastable system) are not yet well understood. Carbon dioxide, a greenhouse gas that is responsible for about 64% of the enhanced “greenhouse effect” of all of the greenhouse Received: December 24, 2010 Revised: February 25, 2011 Published: April 12, 2011 5961

dx.doi.org/10.1021/la105088b | Langmuir 2011, 27, 5961–5967

Langmuir

ARTICLE

Table 1. Interaction Parameters for the Model Systems in This Work

Figure 1. (a) Schematic representation of the simulation system. In this panel, 1.245 nm is the length of a unit cell of the sI hydrate. (b) Snapshot of the final configuration of the system. In this panel, only H2O and SiO2 molecules are shown; CO2 molecules are omitted for clarity. Red, purple, and white spheres indicate O, Si, and H atoms, respectively. Dashed lines are used to represent hydrogen bonds. According to the difference in the local structuring of water molecules, several parts were identified along the z direction, which are denoted by regions A D.

gases, is another type of gas molecule that can form hydrates. For environmental concerns, some authors6,7 proposed a perspective to replace CH4 from the natural gas hydrate with CO2. This idea is of high importance because it relates to both energy production and CO2 sequestration. Liquid CO2 is also proposed to deposit into deep sea and marine sediments in a hydrate form.8,9 In these cases, hydrate formation may depend on the presence of an interface, even though for the storage of CO2 in geological media10,11 the risks of hydrate crystallization in reservoir pores appear when CO2 is injected into porous rock deeply below ground12 because CO2 will expand into depressurized gas reservoirs and then cool down13 because of a large Joule Thomson effect. This phenomenon would block the injection process and affect production safety severely. In general, the formation of gas hydrate on the solid surfaces or in porous media is related to scientific, energy, and safety concerns, climate change, future energy resources, and production safety. In this respect, a detailed knowledge of the mechanism of hydrate formation affected by solid surfaces remains an attractive topic. Computer simulation has been proven to be a powerful tool for extracting detailed information about structural and dynamical properties on a molecular level and extensively used to study the formation of hydrates.14 22 In particular, Walsh et al.22 recently reported direct molecular dynamics simulations of the spontaneous nucleation and growth of methane hydrate. They extended their simulations to the timescale of several microseconds to catch nucleation events with an induction time longer than typical MD trajectories on the nanosecond scale. In this work, molecular dynamics method in a microsecond domain is applied to investigate the kinetic pathways of the nucleation of CO2 hydrates triggered by solid surfaces.

2. MODEL AND SIMULATION METHOD Our MD simulations were performed by using LAMMPS,23 an open source program for massively parallel simulations. To obtain the initial configuration, first we set 800 CO2 molecules and 4600 H2O molecules randomly into a simulation box with the size of 6.22 nm  6.22 nm  4.98 nm (5  5  4 unit cells of sI hydrate), as is shown in Figure 1a. Note that the number of H2O and CO2 molecules can exactly form 100 unit cells (5  5  4) of the sI hydrate with 600 large cages and 200 small cages in which CO2 molecules were assumed to fill in both large and

small cages, although it is somewhat difficult to fill small cages with real CO2 hydrates. Then, in the z direction of the box, two silica layers, where the surfaces were modified by attaching hydrogen ( OH model) to the nonbridging oxygen atoms, were added at the top and bottom of the box, respectively (Figure 1). In our simulations, the extensively studied SPC/E water model24 was employed and the rigidity of H2O was restricted with the SHAKE algorithm.25 CO2 molecules were represented by the EPM2 model.26 The EPM2 model for CO2 molecules has three Lennard-Jones sites with charges centered at each atom and with rigid bond lengths (1.163 Å) but a harmonic bond angle (1/2)kθ(θ θ0)2, where kθ is 1236 kJ/(mol rad2) and θ0 is 180
. The hydroxylated silica model27 was adopted for the SiO2 layers. The Lennard-Jones interaction parameters for H2O SiO2 and CO2 SiO2 were obtained by the Lorentz Berthelot mixing rule, and the parameters for H2O CO2 were taken from ref 28. The parameters for the interaction potentials used in this work are summarized in Table 1. A cutoff radius of 12.0 Å was used for the shortranged interactions, and the long-ranged interactions were evaluated by using the pppm algorithm.29 The periodic boundary conditions were imposed in all three Cartesian directions. Constant temperature and pressure during the NpT runs were maintained at 275 K and 250 bar by using the Nose Hoover algorithm.30 32 First, an NpT relaxation process of 2 ns was performed at 275 K and 250 bar to eliminate the effect of the initial configuration. Then, simulations of the hydrate formation close to the solid surfaces on the timescale of 3 μs were performed with a time step of 2 fs. In all simulations, the positions of SiO2 molecules were fixed. To reveal the structure transition of the CO2 H2O substrate system, we monitored the evolution of a four-body structural order parameter, F4j. The order parameter is a function of the torsional angle between oxygen atoms within 3 Å and the outermost hydrogen atoms in the H2O H2O pair,14,15 which is defined as F4j = (1/n)∑ni = 1cos 3ji. Here, ji is the torsional angle of the ith H2O H2O pair, and n is the total number of H2O H2O pairs. Note that in the bulk system the average values of F4j for hydrate, liquid water, and ice are 0.7, 0.04, and 0.4, respectively.19 5962

dx.doi.org/10.1021/la105088b |Langmuir 2011, 27, 5961–5967

Langmuir

Figure 2. Evolution of F4j for different regions. The evolution of the order parameter for different regions is shown in different colors: region A (black), region B (red), region C (green), and region D (blue). Note that for bulk systems the average values of F4j for hydrate, liquid water, and ice are 0.7, 0.04, and 0.4, respectively.19

Figure 3. Evolution of CO2 densities in different regions. The evolution of the density in different regions is shown in different colors: region A (black), region B (red), region C (green), and region D (blue).

3. RESULTS AND DISCUSSIONS 3.1. Nucleation and Growth of the CO2 Hydrate. Figure 1 shows the final configuration of the system, which is a typical snapshot of the equilibrium configurations. After reaching equilibrium, all of the order parameters of the system will level out versus time. According to the structure transition and nucleation mechanism, the simulation box was divided into four parts along the z direction, which are denoted by regions A D. Our simulation results indicate that the formation of the CO2 sI hydrate can be nucleated from a substrate in a step-by-step manner. First, an icelike layer is formed on the SiO2 substrate (region A). Next to region A and moving toward the center of the system, a thin layer of solid (region B) nucleates with an intermediate structure that links it to region C, which forms (on a microsecond timescale) a CO2 hydratelike structure (CO2 hydrate motif). Region D is a disordered CO2-rich fluid phase. The time evolution of F4j for different regions is given in Figure 2. The order parameter for the region closest to the SiO2 surfaces (region A) rapidly decreases to its equilibrium value of 0.27 at ∼0.13 μs and then fluctuates slightly around this value. This observation indicates that the water molecules are rapidly organized into an icelike structure in region A, whereas CO2 molecules in this

ARTICLE

Figure 4. Evolution of the number of rings in different regions. Numbers 4, 5, 6, and 7 in each panel denote the four-, five-, six- and seven-membered rings, respectively, and symbol Si denotes the ring including Si atoms.

Figure 5. Evolution of the number of cages in different regions. Symbols D and T denote the 512 and 51262 cages, respectively.

region are forced to diffuse into region B (Figure 3) before the complete formation of the icelike structure. As shown in Figure 3, the density of CO2 in region A decreases from the bulk density (∼4.15 nm 3) to nearly 0 within ∼0.07 μs, and the icelike structure is completely formed at ∼0.13 μs (Figure 2). The depletion of CO2 molecules from region A is ascribed to the interaction difference between SiO2 H2O and SiO2 CO2 as well as the structuring of water molecules in this region. The intermolecular interaction of SiO2 H2O is stronger than that of SiO2 CO2 because H2O can form hydrogen bonds with the hydroxylated SiO2 surface. In other words, the hydroxylated silica is essentially hydrophilic. Moreover, the expulsion of CO2 molecules in this layer could also be a result of the icelike structure having no affinity to occlude CO2. Region B plays a key role in the nucleation of a CO2 hydrate from the solid surfaces because it serves as an intermediate structure to reduce the mismatch between the icelike structure in region A and the final stable structure of the hydrate motif. Importantly, the region also acts as nucleation seeds for the further nucleation process. To characterize the structure of region B, we calculated the number of rings and that of cages,33,34 which are shown in Figures 4 and 5, respectively. Figure 4 shows 5963

dx.doi.org/10.1021/la105088b |Langmuir 2011, 27, 5961–5967

Langmuir

ARTICLE

Figure 6. Typical equilibrium structures for different regions. From left to right, the pictures correspond to the snapshot of regions A C. The molecules are colored as described in Figure 1. Note that in these pictures only the corresponding structures are represented by the ball-and-stick model, and other parts are shown in the wire-frame model.

Figure 7. Mean square displacement (MSD) of water molecules in different regions. The MSD is averaged over 2.0 ns near the end of the simulation. The average diffusivity of each region was calculated by using D = limtf¥(1/6t)Ær2(t)æ, and the calculated diffusion coefficients are DA = 0.47  10 9 m2/s, DB = 0.68  10 9 m2/s, DC = 0.55  10 9 m2/ s, and DD = 7.67  10 9 m2/s, respectively, for regions A D.

that after the formation of the icelike layer, the pentagonal rings in region B emerge and gradually grow in number. Unlike region A, in which hexagonal rings dominate the structure, the number of pentagonal rings in region B significantly exceeds that of hexagonal rings (Figure 4b). In region B, although there are no 512 and 51262 cages formed (Figure 5), several semi-512 cagelike structures have been found (Figure 6). The semi-512 cage is not an enclosed structure formed by several pentagonal rings (usually five, six, or seven pentagonal rings). It is an intermediate structure between independent rings and fully developed cages. The distributions of CO2 molecules in the two regions are also different. Although CO2 molecules are thoroughly excluded in region A, the density of CO2 of ∼1.3 nm 3 holds in region B after the completion of the structure transition (Figure 3). The snapshots (e.g., Figure 6b) also indicate that CO2 molecules can enter into the space encaged by some five- and six-membered rings. Figures 2 and 3 show that the structure transition in region B is a two-step process. First, the density of CO2 in this region decreased to ∼2 nm 3 at ∼0.4 μs (Figure 3). Then, the structure of H2O changed subsequently and the structure transition was completed at nearly 0.8 μs (Figure 2). In addition, the mean square displacement (MSD) of H2O molecules (Figure 7) demonstrates that the structure of region B is solidlike. The structure in region B serves as an intermediate to mediate the structure mismatch between the icelike structure near the

substrates (region A) and the sI hydrate motif structure far from the substrate (region C). This point is also confirmed by the appearance of a large number of hexagonal rings in region B (Figure 4b). It is observed that the intermediate structure in region B also serves as nucleation seeds for the further nucleation process. Figure 2 shows that the water structure of region C gradually changes from a liquid state to a hydrate state after the structure formation in region B (∼0.8 μs). The thickness of region C along the z direction is found to match roughly the size of the unit cell of the sI structure, and the induction time of CO2 hydrate formation in this region is about 1 μs (from ∼0.8 to ∼1.8 μs). The density profile for CO2 molecules in Figure 3 shows that a rather high density in the region is observed at the beginning of the crystal growth, ∼0.8 μs, because of the depletion of CO2 from regions A and B. Because H2O molecules in region C are organized into an ordered phase rather slowly, the CO2 density decreases gradually. After ∼1.8 μs, the hydrate motif in region C appears, and the CO2 density reaches ∼3 nm 3, which is very close to the calculated value for a perfect CO2 sI hydrate crystal, 3.11 nm 3. Figure 6c shows a typical structure of 51262 cages after the formation of the hydrate in region C. The number of cages in region C (Figure 5) also proves that the sI hydrate motif is formed at ∼1.8 μs. In our work, the ratio of 512 cages to 51262 cages after the formation of the sI hydrate motif (Figure 5) is larger than that for a perfect sI hydrate because 512 cages are more stable than 51262 cages22,34 as a result of the unfavorable angles for the planar six-membered rings in 51262 cages.1,34 However, the ratio of 512 cages to 51262 cages (Figure 5) in our work is somewhat smaller than the ratio for the homogeneous formation of the CH4 hydrate reported by Walsh et al.22 The difference can be understood as follows. Different from CH4 molecules, the size of CO2 molecules does not fit well in small 512 cages, and thus CO2 can form only the sI hydrate (but CH4 can form both sI and sII hydrates). Thus, the 512 cages formed in our system are mostly empty cages, which decrease their stability. In other words, the 512 cages are formed with more difficulty in our confined system with CO2 molecules than in the bulk system with CH4 molecules, and as a result, the ratio of 512 cages to 51262 cages decreases. Unlike region C, the high density of CO2 molecules in region D (nearly ∼11 nm 3) prohibits the formation of hydrates. Figure 5 shows that in region D the number of 51262 cages has a stronger fluctuation than that of 512 cages. It can be used to explain the strong oscillation of the F4j order parameter in this 5964

dx.doi.org/10.1021/la105088b |Langmuir 2011, 27, 5961–5967

Langmuir

ARTICLE

Figure 8. Evolution of F4j for different simulation systems. In each panel, the evolution of the order parameter for different regions is shown in different colors: region A (black), region B (red), region C (green), and region D (blue). The intersubstrate distance shown is scaled by the size of the unit cell (1.245 nm). With the increase in the intersubstrate distance, the panels are denoted as a i sequentially.

Figure 9. Mean square displacement (MSD) of water molecules in region C for the systems with small intersubstrate distances. The MSD is averaged over 2.0 ns near the end of the simulation.

region (Figure 2). In other words, the oscillation of F4j in region D is mainly caused by the fluctuation of the number of large cages. The mean square displacement (Figure 7) and the diffusivity of region D (7.67  10 9 m2/s) also indicate that the molecules in this region are in the liquid state. Thus, we describe region D as a liquid layer with rings and cages. Because of the space confinement in the z direction, no phase separation (the formation of CO2 bubbles) was observed, although the density of CO2 is rather high in this region. 3.2. Effects of Pore Size. To study the effects of the distance between two silica layers, we performed some additional simulations. In the simulations, different distances between two silica layers were chosen from 2.0 to 6.0 unit cells with an increment of

0.5 unit cell, respectively, while other thermodynamic conditions remain the same as those used in our previous simulations. Again, we used the order parameter F4j to characterize the state of the system and its structure evolution. Figure 8 shows the evolution of the order parameters for different simulation systems. In the Figure, all systems exhibit the same pathway for the nucleation of the CO2 hydrate: first, an icelike layer formed in the region closest to the substrate; then an intermediate structure formed to compensate for the structure mismatch between the icelike layers closest to the substrate and the final bulk CO2 sI hydrate motif; and finally, on the microsecond timescale, the nucleation of the hydrate layer occurs on the intermediate structure, which acts as nucleation seeds. The similarity of these cases demonstrates that the nucleation mechanism and its kinetic pathways for the CO2 H2O substrate systems are independent of the space between the two silica substrates in our simulations. The appearance of the icelike layers occurs at nearly the same time (∼0.15 μs) for different systems, basically being independent of the intersubstrate distance. This is because the formation of the first structured icelike layers is solely determined by the hydrophilic nature of hydroxylated silica. Similar to the formation of the icelike layers, the subsequent formation of the intermediate structures in region B takes place at the same time and nearly the same distance from the substrates, again being independent of the pore size except for very narrow pores. Note that when the intersubstrate distance is smaller than 3.5 unit cells in the z direction the size effect becomes significant and can no longer be neglected. In these cases, the F4j parameter (Figure 8a c) demonstrates that the thin layer with intermediate structure in region B is closer to the hydrate structure. However, for region C in the system with a small intersubstrate distance, both 5965

dx.doi.org/10.1021/la105088b |Langmuir 2011, 27, 5961–5967

Langmuir

ARTICLE

nucleation mechanism and its pathways with respect to the nucleation of the CO2 hydrate from a CO2 H2O substrate system are nearly independent of the space between the two silica substrates except for very narrow pores. The icelike layer for the region closest to the substrate and the thin layer with intermediate structure are found to be independent of the intersubstrate distance. However, the thickness of the first hydrate layer seem to be mildly affected by pore size in the size range we studied, and a layer-by-layer mechanism for the growth of the CO2 hydrate was observed. For the system with a very small intersubstrate distance, no CO2 hydrate can be formed because of the geometric frustration induced by the small system size.

’ AUTHOR INFORMATION Figure 10. Evolution of F4j for layer L1 and that for layer L2 for the system with an intersubstrate distance of six unit cells. In the system, region C contains two layers: layers L1 and L2, each having a thickness of one unit cell. As a comparison, the evolution of F4j in region C is also given. (See also Figure 8.)

the F4j parameter (Figure 8) and the mean square displacement of water molecules (Figure 9) show that the structure in region C is liquidlike, similar to that in region D for the system with a larger intersubstrate distance. Thus, the final CO2 sI hydrate motif disappears for the system with a small intersubstrate distance. This observation indicates that the hydrate structure may be geometrically frustrated by the small system size. For the systems with intersubstrate distances larger than 5.0 unit cells (Figure 8h,i), the inflection in the growth of region C is caused by the layer-by-layer growth of the CO2 hydrate. As an example, we calculated the F4j parameter for the system with an intersubstrate of six unit cells in which region C contains two layers: layer L1 (close to region B) and layer L2 (close to region D). The results are shown in Figure 10. The final configuration for the system also indicates that region C contains two layers of the CO2 hydrate. Figure 10 indicates that in this region the inflection of F4j at ∼1.6 μs in these cases represents an induction period for the second layer. From the evolution of the F4j order parameter, we can clearly obtain a layerby-layer mechanism for the growth of the CO2 hydrate, which is similar to that reported by Jacobson.34 Hence, the hydrate motif layer in region C, which begins to form on the preexisting intermediate structure, seems to be mildly affected by the intersubstrate distance. In other words, the size of region C increases roughly with the system size.

4. CONCLUSIONS In this article, systematic molecular dynamics simulations on the microsecond timescale were performed to illustrate the kinetic pathway of the formation of the CO2 hydrate on hydroxylated SiO2 surfaces. From our simulations, it is found that under the conditions where the CO2 sI hydrate phase is thermodynamically stable, the nucleation of the CO2 hydrate from a solid surface is, in fact, a three-stage process. First, an icelike layer is formed closest to the substrate on a nanosecond timescale. Then, on a submicrosecond timescale, a thin layer with intermediate structure is subsequently formed on the icelike layer to compensate for the structure mismatch between the icelike layers and the final hydrate motif. Finally, on the microsecond timescale, the CO2 hydrate motif layer is nucleated from the intermediate structure that acts as nucleation seeds. The effects of the intersubstrate distance are also addressed here. The simulation results of different systems indicate that the

Corresponding Author

*E-mail: [email protected]; [email protected].

’ ACKNOWLEDGMENT The work is supported by National Natural Science Foundation of China (nos. 20736005, 20876004, and 20925623). Generous allocations of computer time by the Chemical Grid Project of BUCT are acknowledged. ’ REFERENCES (1) Sloan, E. D.; Koh, C. A. Clathrate Hydrates of Natural Gases, 3rd ed.; CRC Press: Boca Raton, FL, 2008. (2) Lee, M. W.; Collett, T. S. Geophysics 2001, 66, 763. (3) Cha, S. B.; Ouar, H.; Wildemann, T. R.; Sloan, E. D. J. Phys. Chem. 1988, 92, 6492. (4) Riestenberg, D.; West, O.; Lee, S.; MaCallum, S.; Phelps, T. J. Mar. Geol. 2003, 198, 181. (5) Yan, L.; Chen, G.; Pang, W.; Liu, J. J. Phys. Chem. B 2005, 109, 6025. (6) Ohgaki, K.; Takano, K.; Sangawa, H.; Sangawa, H.; Matsubara, T.; Nakano, S. J. Chem. Eng. Jpn. 1996, 29, 478. (7) Warzinski, R. P.; Holder, G. D. Energy Fuels 1998, 12, 189. (8) Harrison, W. J.; Wendlandt, R. F.; Sloan, E. D. Appl. Geochem. 1995, 10, 461. (9) Brewer, P. G.; Orr, F. M., Jr.; Friederich, G.; Kvenvolden, K. A.; Orange, D. L. Energy Fuels 1998, 12, 183. (10) Hendriks, C. A.; Blok, K. Energy Convers. Manage. 1993, 34, 949. (11) Bachu, S. Energy Convers. Manage. 2002, 43, 87. (12) Haszeldine, R. S. Science 2009, 325, 1647. (13) Hughes, D. Carbon Capture J. 2008, 2, 2. (14) Rodger, P. M.; Forester, T. R.; Smith, W. Fluid Phase Equilib. 1996, 116, 326. (15) Baez, L. A.; Clancy, P. Ann. N.Y. Acad. Sci. 1994, 715, 177. (16) Moon, C.; Taylor, P. C.; Rodger, P. M. J. Am. Chem. Soc. 2003, 125, 4706. (17) Nada, H. J. Phys. Chem. 2006, 110, 16526. (18) Vatamanu, J.; Kusalik, P. G. J. Phys. Chem. 2006, 110, 15896. (19) Moon, C.; Taylor, P. C.; Rodger, P. M. Faraday Discuss. 2007, 136, 367. (20) Guo, G.; Zhang, Y.; Li, M.; Wu, C. J. Chem. Phys. 2008, 128, 194504. (21) Zhang, J.; Hawtin, R. W.; Yang, Y.; Nakagava, E.; Rivero, M.; Choi, S. K.; Rodger, P. M. J. Phys. Chem. B 2008, 112, 10608. (22) Walsh, M. R.; Koh, C. A.; Sloan, E. D.; Sum, A. K.; Wu, D. T. Science 2009, 326, 1095. (23) Plimpton, S. J. J. Comput. Phys. 1995, 117, 1. (24) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. 5966

dx.doi.org/10.1021/la105088b |Langmuir 2011, 27, 5961–5967

Langmuir

ARTICLE

(25) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (26) Harris, J. G.; Yung, K. H. J. Phys. Chem. 1995, 99, 12021. (27) Lopes, P. E. M.; Murashov, V.; Tazi, M.; Demchuk, E.; MacKerell, A. D. J. Phys. Chem. B 2006, 110, 2782. (28) Duan, Z. H.; Zhang, Z. G. Geochim. Cosmochim. Acta 2006, 70, 2311. (29) Hockney, R. W.; Eastwood, J. W. Computer Simulation Using Particles; Taylor & Francis: New York, 1989. (30) Hoover, W. G. Phys. Rev. A 1985, 31, 1695. (31) Hoover, W. G. Phys. Rev. A 1986, 34, 2499. (32) Melchionna, S.; Ciccotti, G.; Holian, B. L. Mol. Phys. 1993, 78, 533. (33) Matsumoto, M.; Baba, A.; Ohmine, I. J. Chem. Phys. 2007, 127, 134504. (34) Jacobson, L. C.; Hujo, W.; Molinero, V. J. Phys. Chem. B 2009, 113, 10298.

5967

dx.doi.org/10.1021/la105088b |Langmuir 2011, 27, 5961–5967