9106
J. Phys. Chem. C 2008, 112, 9106–9113
Interactions between Structure H Hydrate Formers and Water Molecules Robin Susilo,†,‡ Saman Alavi,*,† Steven Lang,† John Ripmeester,† and Peter Englezos‡ Steacie Institute for Molecular Sciences, National Research Council Canada, Ottawa, ON, K1A 0R6, Canada and Department of Chemical & Biological Engineering, UniVersity of British Columbia, VancouVer, BC, V6T 1Z3, Canada ReceiVed: January 23, 2008; ReVised Manuscript ReceiVed: March 18, 2008
Studies on methane storage in hexagonal structure H (sH) hydrates reveal that hydrate stability conditions, formation rates, and cage occupancies differ for the three large molecule guest substances (LMGSs) studied: tert-butyl methyl ether (TBME), neo-hexane (NH), and methyl-cyclohexane (MCH). In this work, the structures of aqueous solutions with the LMGS are examined with molecular dynamics to gain insights into the interactions between water and LMGS. Water molecules form solvation cavities that follow the contours of the LMGS molecules and have diameters between 4 and 7 Å. The cavity size increases from TBME < NH < MCH, and the nearest water molecules to the LMGS have a local density that is higher by a factor of 2 than that of the bulk water phase. Water molecules interact strongly with the ether oxygen of TBME but not with other atoms in the LMGS molecules. This strong attraction gives TBME a greater solubility in water and better wetting on ice than NH or MCH. Moreover, TBME also has a larger diffusivity than NH in water. The greater contact between TBME and the host water molecules may cause the faster rate of initial hydrate growth from solid ice particles (fixed bed) or liquid water with mixing for this guest species. However, the strong attraction of the TBME molecule with water may also slow down the subsequent crystal growth and may perhaps limit the occupancy of the other guests (methane) in the sH clathrate. The volume change due to solvation is negative for all three LMGSs. This suggests a more compact hydrogen bond network in the hydration shell that may facilitate the formation of clathrates when a hydrate seed is available and when mass and heat transfer resistances are negligible. For example, this may occur during hydrate formation from a fixed ice bed with temperature ramping. Hence, the clathrate growth with hydrophobic guests (NH and MCH) proceeds rapidly toward completion as the temperature is increased above the ice melting point. On the other hand, the formation of hydrate cages for the TBME system is disrupted because of the strong water-solute interactions. Consequently, hydrate formation from ice with thermal ramping with TBME is slower. 1. Introduction Clathrate hydrates, also known as gas hydrates, are naturally occurring inclusion compounds that can be synthesized under suitable temperature and pressure conditions.1 Clathrates are crystalline materials with a physical appearance similar to ice that can store large amounts of gas inside their microscopic lattice cavities. When natural gas components are the guest molecules associated with the clathrate hydrates, these materials sometimes are known as flammable ice, as the gas released from decomposing hydrate can be ignited.2 Interest in natural hydrate research has grown over the past decade mainly due to the potential for use of hydrates in energy and environmental applications.3 The utilization of synthetic clathrates for their gas storage potential has been recognized for several prospective applications that include natural gas transport,4–6 hydrogen storage for fuel cells,7–9 and carbon capture/sequestration.10–16 There are three well-known clathrate hydrate structures, namely, cubic structure I (sI), cubic structure II (sII), and hexagonal structure H (sH). Recently, attention has been directed toward sH clathrates because of their ability to maintain relatively high gas content and their greater stability, for example, in methane storage and transport, as compared to sI and sII hydrates.17 All three hydrate structures are known to * Corresponding author e-mail:
[email protected]. † National Research Council of Canada. ‡ University of British Columbia.
naturally occur, with sH clathrates having been reported only recently.18 Structure H clathrates have the largest cage among all hydrate structures, an icosahedron, and they require two different sizes of guest molecules to stabilize the crystal.19 The large cages accommodate large guest molecules, whereas the smaller cages are filled with smaller molecules such as methane, xenon, or carbon dioxide, which by themselves, are known as typical sI clathrate formers. In the presence of these small molecules, the sH clathrates form at lower pressures than the corresponding sI hydrates. The smaller molecules (e.g., methane) may also occupy the large cages of sH hydrate with multiple occupancy, although this has been reported only at extremely high pressures.20–22 The formation of sH clathrates involves at least three different molecules in three separate phases, which adds to the complexity of these systems. The sH clathrate formers, also known as large molecule guest substances (LMGS), are generally poorly miscible with water, and the smaller guest molecules are commonly present as a gas phase. As a result, a gas phase, a nonaqueous liquid (LMGS) phase, and an aqueous phase have to be in contact in order for the sH clathrate to form. Earlier studies on methane storage in sH hydrates at our laboratories found that the formation rates and methane occupancies differ among the three LMGS studied: tert-butyl methyl ether (TBME), neo-hexane (NH), and methyl-cyclohexane (MCH). The initial rate of clathrate formation grown from solid ice powder (without stirring)22,23 and liquid water with
10.1021/jp8006848 CCC: $40.75 2008 American Chemical Society Published on Web 05/21/2008
Structure H Hydrate Formers and Water Molecules efficient mixing by agitation24 were found to increase in the order MCH < NH < TBME. The system with TBME showed the fastest initial clathrate growth rate, which is consistent with the fact that TBME is much more soluble in water than NH and MCH25 and that TBME wets ice much better than the other two liquids.22 Surprisingly, the hydrate composition determined by solid-state NMR spectroscopy and the gas content measured by decomposing the clathrate revealed that the methane content with TBME was the smallest, followed by NH and MCH.27 The methane pressure required to form a stable sH clathrate with TBME is also higher than with NH or MCH.27,28 It was also observed during hydrate formation from ice powder that as the temperature was ramped above the ice-point, the system with TBME showed a slower transformation, whereas the other two systems proceeded quickly toward full clathrate conversion.23 Consequently, TBME is not the best LMGS because the total time required to convert ice/water fully into clathrate is actually longer, and the methane content is also less than for the NH/MCH systems. The selection of an appropriate LMGS is, therefore, essential to ensure that the clathrate fits the requirements for gas storage applications, which are (a) a high gas content and (b) a fast hydrate formation rate. The above experimental observations and measurements motivated us to study the aqueous solvation of the LMGS at the molecular level. More specifically, the spatial distributions and orientations of water molecules surrounding the LMGS molecule were determined, and the volume of solvation, solution density, and solute diffusivity in water are also reported. The results are related to the solubility and wetting properties of LMGS and rates of hydrate formation.
J. Phys. Chem. C, Vol. 112, No. 24, 2008 9107 oxygen atoms and the protons of water. Consequently, it is more difficult to implement those polarizable water potentials in MD and hence the fixed-charge model, which is good for qualitative predictions, is used for this first study. All LMGS molecules are considered flexible with the initial structure determined by energy optimization using density functional theory at the B3LYP/6-311++G(d,p) level with the GAUSSIAN 98 suite of programs.42 The electrostatic charges are estimated from the electrostatic potential grid (CHELPG)43 method as implemented in the GAUSSIAN program, and the van der Waals interactions of atoms and intramolecular force constants are obtained from the general AMBER force field (GAFF).44 The point charges, van der Waals Lennard-Jones intermolecular potential parameters, and atomic labels from the AMBER force field are provided in our previous work.45 Intramolecular force constants for the bond stretches, the angle bending, and the dihedral interactions are taken directly from the GAFF tables. The van der Waals interactions among guest-guest and guest-host molecules are based on Lennard-Jones (12-6) potential. Coulombic interactions are used to model the electrostatic intermolecular interactions between point charges qi and qj located on the atomic nuclei i and j. Long-range electrostatic interactions were calculated using the Ewald summation method33,34 with a precision of 1 × 10-6, and all intermolecular interactions in the simulation box were calculated within a cutoff distance of Rcutoff ) 13.0 Å. The standard combining rules, ij ) (ii ij)1/2and σij ) (σii+ σjj)/2 are used for Lennard-Jones potential parameters between unlike atomtype force centers i and j. The intermolecular potential is given by eq 1,
2. Molecular Dynamics Methodology The DL_POLY molecular dynamics (MD) simulation package version 2.14 was employed in this study.29 The simulations were performed at 274 K and 2 MPa with the Nose´-Hoover thermostat-barostat algorithm30,31 and the modification of Melchionna et al.32 with thermostat and barostat relaxation times of 0.5 and 2.0 ps, respectively. The equations of motion were solved using the Verlet leapfrog algorithm with a time step of 1 fs.33,34 Energy conservation and configurational energy convergence criteria were used to test the appropriateness of this time step. The aqueous solution structures were determined at constant temperature and pressure (NPT) and the mean square displacements (MSD) were obtained at constant volume and energy (NVE) simulations. The extended simple point charge (SPC/E) model is used to describe the intermolecular potential for water.35 Water molecules are considered rigid and are treated with the SHAKE34 algorithm. The fixed-charge SPC/E model has been used in clathrate simulations and gives qualitatively and semiquantitatively correct predictions regarding the behaviors of these inclusion compounds. The use of polarizable potential models for water36,37 leads to a more accurate representation of molecular dipole and charge polarization enhancements in the methane clathrate. The calculated vibrational spectrum of methane hydrate in the water libration region was also found to depend sensitively on the nature of the water potential used in the classical force field. The AMOEBA,38 COS/G2,39 SWM4NDP,40 and GCPM41 are among polarizable water models that have been used to predict the properties of solvation and hydrogen bonding in aqueous solutions more accurately. By allowing a degree of charge mobility around the atoms, polarizable potentials provide a more accurate representation of the interactions between the nonbonded electron pairs of the
V (inter) )
∑∑ i)1 j>i
{ [( ) ( ) ] 4εij
σij rij
12
-
σij rij
6
+
qiqj 4πε0rij
}
(1)
The simulations for aqueous solution were performed for a minimum total time of 250 ps, and the initial 50 ps was used to equilibrate the system. The solution consists of 918 water molecules and a single solute (LMGS) molecule with periodic boundary and initial dimension of 33.7 × 33.7 × 28.0 Å3. Snapshots of the solution are collected at every 0.2 ps interval and 1000 snapshots of the solutions were used to acquire structural data for analysis. The spatial data from each snapshot is projected onto the reference coordinate systems shown in Figure 1. The C1-O-C2 atoms for TBME, C1-C2-C3 atoms for MCH, and C6-C5-C4 atoms for NH atoms are selected to define the reference axes. The two end carbon atoms define the Y-axis in each case. The X-axis is drawn perpendicular from the Y-axis through the central atom (oxygen atom from the TBME molecule or analogous carbon atoms in NH and MCH). The Z-axis faces outward perpendicular to the XY-plane. The reference axes are chosen to arrange analogous structural features of the solutes in a similar manner. A rotation and translation matrix based on the reference LMGS molecules is used to transform the coordinates of the water molecules in the simulation. The LMGS studied are considered as flexible molecules, and hence, vibrational motions are allowed for the atoms that are not used in defining the coordinate axes. The water spatial distributions are obtained by analyzing all the 1000 snapshots after transforming the coordinates to the reference system based on the LMGS molecules. The water molecule positions are stored in 3-dimensional bins 0.5 Å per side. The atomic (hydrogen or oxygen) densities in each box are normalized to a pure bulk water simulation. Water molecular
9108 J. Phys. Chem. C, Vol. 112, No. 24, 2008
Susilo et al.
Figure 1. Definitions of the coordinate systems with the LMGS molecule as the reference (top panels). The atoms selected represent similar geometries. The colors for the carbon, oxygen, and hydrogen atoms are cyan, red, and white, respectively. Thermal vibrations are seen for the free atoms of the LMGS molecules in the lower panels.
Figure 2. The definition of the orientation of water molecule relative to the reference atom (RA) as the origin. The reference atoms are O (TBME), C5 (NH) and C2 (MCH), see Figure 1. The water molecule is represented in red. The distance between the reference atom and the oxygen-water is designated as R. The orientation of the water molecule is set as the direction of the moment dipole of water (µ). The angular orientation of water with respect to the reference frame is given in terms of the cosine of the angle ω.
distributions and orientations up to 10 Å away from the solute are analyzed. Figure 2 illustrates the definitions of the water orientation (µ) and radial distribution (R) vectors and the angle (ω) between these two vectors. The water orientation vector µ is drawn from the oxygen atom along the bisector of the HOH angle. The angle ω determines the angular orientations of water molecules with respect to the referenced atom. The water orientation distribution is also normalized relative to the random orientation of water molecules in the bulk phase. 3. Results and Discussion 3.1. Spatial Distributions of Water Molecules near the Solutes. The radial distribution functions (RDFs) of the reference-atom from the LMGS (O-TBME, C5-NH, and C2-MCH)
with oxygen-water (OW) and hydrogen-water (HW) are plotted in Figure 3. For comparison, the corresponding RDFs of pure water oxygen (OW) are also plotted. The HW-OW and OWOW interactions clearly indicate the presence of hydration shells. The strong first OW-OW shell appears at less than 3.3 Å and the first HW-OW shell at less than 2.5 Å, which are typical hydrogen bond distances between the reference water molecule and the nearest-neighbor shell of waters. The second hydration shell is centered at ∼4.5 Å for OW-OW atoms (or ∼3.8 Å for OW-HW atoms) and is broader and less structured than the first shell. The water molecules in the second hydration shell are mostly those that are hydrogen bonded to the first hydration shell. A weak third hydration shell is also visible between 5.7 and 7.7 Å (OW-OW). A strong O-HW peak in the first hydration shell below ∼3 Å is also seen in the TBME solution, indicating the presence of hydrogen bonding between the water and oxygen from TBME, but this peak is absent for the other two solutes. The TBME RDFs also have broad peaks in the range of 4-6 Å, which represents the solvation sphere of the hydrophobic parts (methyl groups) of the molecule. The RDFs of the NH and MCH solutions have broad peaks in the range of 4-6 Å, representing the solvation sphere of the hydrophobic molecules. There are only weak longer range peaks for these solutes. The RDF plots provide the radial distribution of water molecules around the solute. To gain more insight on the solvation process and its relationship with clathrate formation and stability, information on the water spatial distribution and orientation with respect to the solutes are needed.
Structure H Hydrate Formers and Water Molecules
J. Phys. Chem. C, Vol. 112, No. 24, 2008 9109
Figure 3. Radial distribution function (RDF) of solute reference atoms (O-water, O-TBME, C5-NH, and C2-MCH) with water-oxygen OW (top) and water-hydrogen HW (bottom).
The spatial distributions of water oxygen and hydrogen atoms near the solutes in the XY-plane at Z ) 0 are shown in Figure 4. The oxygen and hydrogen density maps are similar, and the water spatial distribution function near the solute shows a local density that is approximately twice that of the bulk water phase. The solute hydrophobic groups weakly interact with the water molecules. Consequently, the water molecules in the hydration sphere attract each other more strongly so as to minimize the energy loss due to hydrophobic moiety–water interactions. This high density water in the hydration shell was also reported for water-methane,46 water-cyclohexane, and water-benzene systems.47 The thickness of this high density hydration shell is less than ∼1 Å. In Figure 4 it is also seen that the hydrogen bonding between the O-TBME and HW atoms distorts the hydration shell for the TBME system. Figure 5 illustrates the three-dimensional spatial distribution of water oxygen atoms surrounding the reference atom at 2.7 Å (green), 4 Å (white), and 5 Å (gray). The shortest distance (2.7 Å) is characteristic of hydrogen-bonding and corresponds to water-water and water- TBME interactions. The hydrogen bonds of water are present as two caps above the hydrogen atoms of the reference water and a wide distribution at the oxygen site (right-hand side) of the reference water.48 The latter is also observed in the water-TBME system but not with the NH or MCH system. The water molecule donates a proton to the oxygen atom of TBME. However, the distributions of H-donor water molecules are not as widely spread as the water-water system due to the hydrophobic nature of the methyl group. At a distance 4 Å away from the reference atom, the water molecules begin to form a shell that encapsulates the hydrophobic parts of the LMGS. The LMGS molecules are almost completely enclosed by a water molecule network at 5 Å except for the hydrophobic part on the lower left for both TBME and NH molecule and on the lower right for the MCH molecule.
Figure 4. Spatial contours for OW and HW distributions around the solute molecules in the XY-plane at Z ) 0 with respect to bulk water. The water density near the solutes is approximately twice that of bulk water.
The cavity size increases upon going from TBME to NH and MCH. The TBME and NH molecule are rather similar in profile and molecular weight and, hence, a comparison between the two is particularly interesting. The only difference between the two is that the oxygen atom in TBME is replaced by a CH2 group in NH. It is clear from Figure 5 that water encloses the TBME molecule to a greater extent than the NH molecule at the same distance from the reference atom due to the stronger attraction between TBME and water molecules. 3.2. Spatial Orientation of Water Molecules near LMGS Molecules. The orientational density maps of water around the solutes are shown in Figure 6, where the distribution of the cos(ω) for the water molecules is plotted as a function of the radial distance R of the water molecules away from the reference atom. The orientational density is normalized to the orientational correlations in the pure water system. In the solvent, at a separation of ∼2.5 Å, water molecules have their preferential orientational densities approximately 18 times larger than the random distribution in the bulk phase. The water molecules that accept a hydrogen from the reference water mostly face away from the reference atom (cos ω ≈ -1). The proton-donating water molecules preferably face the reference atom with angles between 45-55° (cos ω ≈ +0.6 to +0.7). Water molecules in the second hydration shell (∼4 Å) are mostly oriented either facing toward (H-donor, cos ω ≈ +1), away from (H-acceptor, cos ω ≈ -1), or tangential to (cos ω ≈ 0) the reference atom. The water molecules in the bulk phase (R > 6 Å) have no preferred orientation. The proton-donor water molecules at ∼3 Å that are hydrogenbonded to the oxygen atom of TBME face toward the reference
9110 J. Phys. Chem. C, Vol. 112, No. 24, 2008
Susilo et al.
Figure 5. Distribution of the water molecules with respect to the solute at selected distances. Strong short-ranged hydrogen bonds for water-water and TBME-water are observed.
TABLE 1: Volumes of Water, LMGS, and Water-LMGS Solutions component/number of molecule
volume (Å3)
Fcalc.a (g/cm3)
water (918) TBME (1)c NH (1)c MCH (1)c water (918) + TBME (1) water (918) + NH (1) water (918) + MCH (1)
27196 ( 208 197 ( 7 204 ( 5 210 ( 5 27390 ( 196 27364 ( 192 37375 ( 197
1.009 ( 0.008 0.742 ( 0.027 0.700 ( 0.017 0.776 ( 0.020 1.007 ( 0.007 1.008 ( 0.007 1.008 ( 0.007
Fexpt.b (g/cm3) 1.000 0.741 0.644 0.769
a Density values from the simulation b Experimental density values from ref.45 c From separate simulations of 27 LMGS molecules.
Figure 6. Orientational distributions of water with respect to a reference atom from OW, O-TBME, C5-NH and C2-MCH. The distance and angle variables are defined in Figure 2. At short distances, there is a strong preferred orientation for each system. Higher orientation density is also observed due to the higher density of water.
atom with angles ω ≈ 45-55°. The water orientation density is a factor of 6 higher than the bulk water phase but is three times lower than the water-water system. The water molecules up to ∼6 Å from TBME are preferentially oriented away from the reference atom by 90-125° (cos ω ≈ 0 to -0.6) due to the hydrophobic nature of the methyl groups. NH and MCH have
orientational densities similar to the hydrophobic part of TBME but in a broader range, from 50-135°. There are no orientational correlations for water molecules that are further than 7 Å from the reference atoms in all LMGS. 3.3. Solvation Volumes. The volumes of pure solvent (water), pure solutes (LMGSs), and the LMGS-water solutions are summarized in Table 1. The volumes reported correspond to the number of molecules indicated in the bracket. The volume change differences between the solution volume and the sum of the separate solute and solvent volumes are generally small negative values for all systems. Although the computational quantification of molar volume change is not possible due to the small number of solute molecules present in the solutions studied, the decrease in molar volume for hydrophobic solutes such as methane, ethane, propane, and benzene has been reported in several experimental49 and theoretical50 studies. The
Structure H Hydrate Formers and Water Molecules
Figure 7. Mean square displacement (MSD) of TBME and NH in water over a time range of 60 ps.
volume reduction can be attributed to stronger hydrogen bonds between the solvating water molecules in order to minimize the hydrophobic water-guest interactions. Hence, the water molecules are expected to be more structured when solvating the nonpolar solutes. The density of the pure solute and solvent from the simulation are comparable to literature values.51 The density of all solutions is slightly lower than the water density due to the presence of the LMGS molecule. It should be noted that the magnitudes of the volume reductions in the solutions given in Table 1 are semiquantitative. For better quantitative estimates for the volume of solution, more accurate polarizable models for water (such as TIP4P-FQ52) and the LGMS (CHARMM-FQ53) should be used. The use of polarizable models would additionally lead to a more accurate representation of hydrogen bonding in the solvation shell of the LGMS. 3.4. Solute Diffusivity in Aqueous Solution. The dynamics of LMGS molecules as solutes in water can be observed from the mean square displacement (MSD) plot, which determines the diffusion coefficient. Figure 7 shows the time variation of the MSD of the TBME and NH molecules. The MSD profiles of TBME and NH show similar structures but differ in quantitative features. For time scales
