Dynamical Properties of Hydrogen Sulphide Motion in its Clathrate

Bioprocess Engineering, University College Dublin, Belfield, Dublin 4, Ireland ..... bonds prior to breakage and possible subsequent reformation o...
0 downloads 0 Views 3MB Size
Download PDF
ARTICLE pubs.acs.org/JPCA

Dynamical Properties of Hydrogen Sulphide Motion in its Clathrate Hydrate from Ab Initio and Classical Isobaric-Isothermal Molecular Dynamics Niall J. English*,† and John S. Tse*,‡ †

The SEC Strategic Research Cluster and the Centre for Synthesis and Chemical Biology, School of Chemical and Bioprocess Engineering, University College Dublin, Belfield, Dublin 4, Ireland ‡ Department Physics and Engineering Physics, University of Saskatchewan, Saskatoon, Saskatchewan S7N 5E2, Canada ABSTRACT: Equilibrium ab initio (AI) and classical constant pressure-constant temperature molecular dynamics (MD) simulations have been performed to investigate the dynamical properties in (type I) hydrogen sulphide hydrate at 150 and 300 K and 1 bar, and also lower temperatures, with particular scrutiny of guest motion. The rattling motions of the guests in the large and small cavities were around 45 and 75-80 cm-1, respectively, from AIMD, with the corresponding classical MD modes being 10-12 cm-1 less at 150 K and around 5 cm-1 lower at 300 K. The rattling motion in the small cavity overlapped somewhat with the translational motion of the host lattice (with modes at circa 85 and 110 cm-1), due in part to a smaller cage radius and more frequent occurrences of guest-host hydrogen bonding leading to greater coupling in the motion. The experimentally determined H-S stretch and H-S-H bending frequencies, in the vicinity of 2550-2620 and 1175 cm-1 [H.R. Richardson et al., J. Chem. Phys. 1985, 83, 4387] were reproduced successfully in the AIMD simulations. Consideration of Kubic harmonics for the guest molecules from AIMD revealed that a preferred orientation of the dipole-vector (or C2-axis) exists at 150 K vis-a-vis the [100] cube axis in both the small and large cavities, but is markedly more significant for the small cavity, while there is no preferred orientation at 300 K. In comparison, classical MD did not reveal any preferred orientation at either temperature, or at 75 K (closer to the AIMD simulation at 150 K vis-a-vis that approach’s estimated melting point). Probing rotational dynamics of the guests reinforced this temperature effect, revealing more rapid rotational time scales at 300 K with faster decay times of dipole-vector (C2) and H-H-vector (C2y) being similar for each cage, at around 0.25 and 0.2 ps, respectively, versus approximately 0.45 and 0.5 ps (large) and 0.8 ps (small) at 150 K. It was found that the origin of the observed preferred orientations, especially in the small cages, at 150 K via AIMD was attributable to optimization of the dipolar interaction between the guest and outward-pointing water dipoles in the cavity, with guests “flitting” rotationally between various such configurations, forming occasionally hydrogen bonds with the host molecules.

’ INTRODUCTION Clathrate hydrates are nonstoichiometric crystalline inclusion compounds in which a water host lattice encages small guest atoms or molecules in cavities; the empty lattice is thermodynamically unstable, and its existence is due to hydrogen bond stabilization, resulting from the enclathration of the trapped solutes in its cages.1,2 There are three known common hydrate structures: (s)I, II, and H. In type I hydrate, the unit cell is formed from two small 512 pentagonal dodecahedral cavities and six slightly larger tetrakaidecahedral 51262 cages, with 46 water molecules.1,2 For type I hydrate, there are up to 8 guests per unit cell (i.e., one per cavity). Hydrogen sulphide (H2S) clathrate holds a particularly important position in hydrate research. In the natural environment, gaseous H2S hydrate is often found together with methane hydrate in oil and gas wells. Therefore, the understanding of the thermodynamic and physical properties of pure and CH4-H2S hydrates are relevant to the exploration of methane hydrate. Although H2S is immiscible with water, H2S hydrate also has a r 2011 American Chemical Society

peculiarly high melting point, larger than that of ice, indicating strong guest-host interaction not presented in other gas hydrates. Richardson et al.3 determined the H-S stretch and H-S-H bending frequencies, in the vicinity of 2550-2620 and 1175 cm-1, respectively. The purpose of this article is to investigate the dynamical properties of fully occupied (type I) hydrogen sulphide hydrates. The objective is to employ both ab initio (AI) and classical MD methods to scrutinize the guest dynamics and determine the underlying factors governing their motion.

’ METHODOLOGY The higher meting point of H2S hydrate than ordinary ice (Ih) suggests strongly that the guest-host interaction may be Special Issue: Victoria Buch Memorial Received: December 2, 2010 Revised: February 11, 2011 Published: March 10, 2011 6226

dx.doi.org/10.1021/jp111485w | J. Phys. Chem. A 2011, 115, 6226–6232

The Journal of Physical Chemistry A dominant in this clathrate. Flexible classical potentials are parametrized usually at a single temperature (often near or at ambient temperature), and therefore, this is a drawback when assessing temperature affects on intramolecular motion. Therefore, it is essential to incorporate possible quantum effects in molecular dynamics calculations (AIMD), particularly at low temperature (150 K). AIMD also serves as a useful comparison to classical MD for consideration of other dynamical properties of the guest molecules. In this study, the Perdew-Becke-Ernzerhof (PBE) functional4 was used, within the pseudopotential localized basis set approximation, as implemented in the density functional (DFT) SIESTA code.5 Isobaric-isothermal (variable-cell NPT) Born-Oppenheimer molecular dynamics (BOMD) was carried out at 1 bar with a time step of 1 fs. Split-valence plus polarization (DZP) basis sets were employed for oxygen, sulfur, and hydrogen atoms. The Coulomb potential was expanded in a plane-wave basis with an energy cutoff of 400 Ryd; a very large plane-wave expansion is essential to ensure convergence in the calculated total energy and atomic forces. For instance, the temporal values of the cell dimensions, total energy, and temperature during the MD calculations at 150 K are shown in Figure 1. Apart from the initial 2000 one-femtosecond time steps required for the equilibration of the system, all of the calculated quantities are well converged (also at 300 K after a 2 ps relaxation), indicating the appropriateness of the chosen computational parameters. For classical MD, the SPC/E water model was used for water-water interactions6 and the rigid potential of Jorgensen et al. for H2S.7 Rigid potentials were employed, as it was felt that AIMD would be superior to flexible classical potentials for gauging intramolecular motion (vide supra). Lorentz-Berthelot combining rules were used for the H2S-water Lennard-Jones (LJ) interactions.8 The intermolecular potential between the Jorgensen’s potential for H2S with the SPC/E model for water involves only the LJ term between sulfur and oxygen atoms. Usually geometric mixing rules were used in this case,7 however, it was found from test simulations that both geometric and Lorentz-Berthelot combining rules gave essentially identical results for all calculated quantities. This can be understood clearly, because the S-O LJ radii are 3.433 and 3.423 Å for Lorentz-Berthelot and geometric mixing cases, respectively. The cutoff radius for Lennard-Jones interaction parameters was 10 Å. The Ewald method was used to handle long-range electrostatic interactions.8,9 The real-space cutoff distance for the Ewald method was 11 Å, and the screening parameter and the number of wavevectors were set such that the relative error in the Ewald summation was less than 1  10-6. In practice, this led to the product of the screening parameter R and the real-space cutoff distance to be in the range of around 3.2 to 3.5, and the Ewald electrostatic energy and forces varied very little with R and increasing number of wavevectors. The velocity Verlet scheme was used for MD,8 and the RATTLE method used to impose constraints.10 For extended system dynamics in the NVT ensemble, light coupling to a Nose-Hoover thermostat was used, with a thermostat relaxation time of 1 ps.11 Where dynamics was carried out in the NPT case, light coupling was applied with Melchionna’s modified form of the Hoover barostat using isotropic cell fluctuations,12 with thermostat and barostat relaxation times of 1 and 2.5 ps, respectively. A 1 fs time step was used for all simulations and found to be highly satisfactory for energy conservation in the NVE ensemble. The starting coordinates of the oxygen atoms were those of the unit cell of sI.13 The initial unit cell length was 12.03 Å for the

ARTICLE

Figure 1. Evolution during variable-cell NPT BOMD at 150 K of (a) cell size, (b) cell angle, (c) temperature, and (d) total energy. Data is shown for the full 25 ps trajectory, but the first 2 ps was rejected for analysis, owing to equilibration. The curves in red, green and black color in (a) and (b) are the temporal variation of the lengths and angles of the model unit cell.

cubic sI system. The initial orientations of the water molecules were selected in a random manner so as to conform to the BernalFowler rules,14 and so that the total dipole moment of the system would be vanishingly small. The Rahman-Stillinger procedure was used to achieve a small total dipole moment.15 H2S molecules were placed in the cavities of each unit cell to allow for full occupation, i.e., 8 guest molecules. For AIMD, a single unit cell containing 46 water molecules was used. For classical MD, a supercell constructed from 3  3  3 replication of the unit cell was employed. Although the sampling time for rotational motion of water molecules in hydrates is much longer than the 18-23 ps of time for AIMD used in this study (vide infra),16 the different rotational water-molecule configurations in each of the cavities led to reasonable statistical averaging of guest (rotational) motion in different cage proton disordered environments. NPT relaxation was initially carried out at 1 bar at 150 and 300 K for 100 ps using MD with the classical potentials. Once the 6227

dx.doi.org/10.1021/jp111485w |J. Phys. Chem. A 2011, 115, 6226–6232

The Journal of Physical Chemistry A

ARTICLE

system volumes had reached steady values with mean pressures at their set values, 200 ps of production MD was performed. AIMD was carried out at 100, 150, and 300 K. Usually, a period of 2 ps was needed to reach equilibration (cf. Figure 1). This relaxation period was then followed by 18-23 ps of variable-cell production MD. It was found that, at 100 K, even though the lattice parameters and total energy converged readily within this time regime, the energy exchange between water and H2S was very slow. Therefore, the results at this temperature will not be discussed in detail. The slow equilibrium between the water and H2S at 100 K is probably due to the fact that PBE overestimates the water interaction as the theoretical melting point of ice Ih 417 K.17 From a simple scaling argument, the “true” temperature of this system is close to 73 K. Therefore, considering this, we also carried out classical MD simulations at 75 K, to compare more directly to the AIMD results at 150 K. The calculated mean sizes of the cubic unit cells at 100, 150, and 300 K were 11.24, 11.27, and 11.33 Å, respectively: the resultant estimated coefficient of linear thermal expansion of 4  10-5 is comparable to the experimental measured value of around 4-7  10-5 in structure I methane and ethylene oxide hydrates.18 To investigate the frequency modes responsible for rattling motion of the guest molecules in the cage, as well as that of the lattice, (normalized) velocity autocorrelation functions (VACF’s) were evaluated for each atom type i, that is, ZðtÞ ¼ Ævi ð0Þ 3 v i ðtÞæ=Ævi ð0Þ 3 vi ð0Þæ

ð1Þ

VACFs measure the degree of significance of coupling of atomic motions with themselves. Power spectra (Fourier transforms) of the VACFs reveal the frequency modes contributing to the motion. For the lattice, power spectra reveal the densities of states (DOS); the oxygen atoms in water dominate the translational motion of the lattice, while the hydrogen atoms reflect librational motion.19 The VACFs of sulfur and hydrogen atoms in the guests were also studied in the large and small cavities, to characterize H-S stretch and H-S-H bend frequency modes from AIMD, and also to study any coupling of guest rattling modes with the lattice translational DOS. To study possible preferential alignment of H2S in the cavities, the orientation distribution function f(r) is expanded in terms of Kubic harmonics Kn(r), 4πf ðrÞ ¼ 1 þ C4 K4 ðrÞ þ C6 K6 ðrÞ þ :::

ð2Þ

where r = (x,y,z) is a unit vector specifying a direction of interest, and in the case of H2S, it is the C2ν symmetry axis. K4 ðrÞ ¼ ð21=16Þ1=2 ð5Q - 3Þ

ð3Þ

K6 ðrÞ ¼ ð13=128Þ1=2 ð462R þ 21Q - 17Þ

ð4Þ

with Q = x þ y þ z and R = x y z . The expansion coefficients C4 and C6 are equal, respectively, to ÆK4æ and ÆK6æ. In this study, K4 and K6 were considered in the large and small cavities.20 A nonzero ÆK4æ denotes preferential alignment relative to the [100] cube axis. The (normalized) autocorrelation functions (ACFs) of these quantities in both cage types were also monitored. In an effort to characterize rotational motion, the ACF of the H2S dipole vector, μ, which coincides with the H2S C2ν-symmetry axis, was also computed 4

4

4

2 2 2

ZðtÞ ¼ Æμi ð0Þ 3 μi ðtÞæ=Æμi ð0Þ 3 μi ð0Þæ

ð5Þ

alongside its counterpart (i.e., the “y”-axis joining the two hydrogen atoms in H2S as the principal axis of rotation in Cy2

Figure 2. Power spectrum of oxygen VACF, the translational density of states, from AIMD at (a) 150 and (b) 300 K.

Figure 3. Power spectrum of water hydrogen atoms’ VACF, the librational density of states, from AIMD at (a) 150 and (b) 300 K.

H2S)21,22 from consideration of the more general form cRl ðtÞ ¼ ÆPl ½eRi ðtÞ 3 eRi ð0Þæ=ÆPl ½eRi ð0Þ 3 eRi ð0Þæ

ð6Þ

where Pl is a Legendre polynomial and eRi is a unit vector along molecule i’s principal axis of rotation. The decays of these ACFs may be modeled to first-order approximation as exponential, that is, A exp(-t/τ), and characteristic relaxation times τ derived to describe characteristic rotational time scales.

’ RESULTS AND DISCUSSION The translational and librational DOS of the host lattice, i.e., power spectra of the VACF’s for oxygen and hydrogen atoms in water, are shown in Figures 2 and 3 for 150 and 300 K from AIMD, respectively. The lattice translational modes occur at circa 85 and 110 cm-1, with a very marginal shift toward higher frequencies in the translational and librational DOS at 150 K; results from classical MD were similar, except the H-O-H bending and O-H stretch mode were not manifest in the librational DOS due to the use of rigid potentials, ipso facto, and qualitatively consistent with previous host-DOS results for other types of clathrate.19,20 The vibrational density of states (VDOS, power spectra) for the sulfur and hydrogen atoms in the guests are shown in Figures 4 and 5, respectively, in each cage type at both 150 and 300 K from classical MD and AIMD. For sulfur, the (classical) MD modes were 10-12 cm-1 less at 150 K and around 5 cm-1 lower at 300 K vis-a-vis AIMD, while neither MD nor AIMD spectra exhibit 6228

dx.doi.org/10.1021/jp111485w |J. Phys. Chem. A 2011, 115, 6226–6232

The Journal of Physical Chemistry A

ARTICLE

Figure 5. Power spectrum of H2S hydrogen atoms’ VACF from AIMD and classical MD at (a) 150 and (b) 300 K. Figure 4. Power spectrum of sulfur VACF from AIMD and classical MD at (a) 150 and (b) 300 K.

marked temperature dependence. Classical MD results at 75 K were similar to those at 150 K (data not shown for clarity). However, it is quite obvious that the vibrational features in the S powder patterns are more distinct at 150 K than at 300 K. As noted above, the classical potentials are rigid and therefore do not display any intramolecular features in the guest hydrogen spectrum. Considering AIMD results, the rattling motion in the small cavity (at around 80 and 100 cm-1, cf. Figure 4) overlaps somewhat with the translational motion of the host lattice (circa 85 and 110 cm-1, cf. Figure 2), while the rattling motion in the large cavity is around 45-48 cm-1. Incidentally, AIMD predicted a weak vibrational band at about 45 cm-1 at the lower energy side of the main translational acoustic maximum of the water lattice at 85 cm-1 (Figure 2), 150 K, which was missing at 300 K, is suggestive of the hybridization of water and H2S vibrational modes at low temperature. The greater apparent coupling in guest rattling motion with the host in small cavities is due primarily to a smaller cage radius and resultant closer contacts between the H2S and adjoining water molecules. To probe this further, the Luzar-Chandler geometric hydrogen bond identification criteria23 were employed to assess the probability of occurrence and distribution of persistence (survival) times of each individual H2S-water hydrogen bond on the AIMD trajectories. The H2S hydrogen atom acts occasionally as a donor to oxygen atoms in some of the adjoining water molecules of its cage; a water-H2S pair was considered to be hydrogen-bonded if the oxygen-sulfur distance is less than 3.5 Å and the S-H 3 3 3 O angle is greater than 150 (with all S-H 3 3 3 O angles classified as being between 0 and 180).23 It was found that more frequent, transient occurrences of guest-host hydrogen bonding took place within the small cavities, typically with greater persistence times vis-a-vis the large cages. The persistence times at 150 K tended to be greater than at

300 K, especially in the small cages, owing to a smaller radius and closer guest-host contact. For instance, it was found that in small cages at 150 K, the guest had hydrogen bond contacts (though never more than one at a time) with 9-10 of its 20 water-molecule neighbors for between 0.4 and 2.5% of the simulation time for each individual bond and 11% in terms of the temporal presence of any hydrogen bond at any given time; typical persistence times in the 15-80 fs range, with 30-50 fs being the most probable for survival times of individual hydrogen bonds prior to breakage and possible subsequent reformation of the order of picoseconds later. For some water molecules, hydrogen bond contact was more popular, with both guest hydrogen atoms bonded thereto at various times (i.e., two distinct, but closely related hydrogen bonds), owing to guest rotational motion and preferred orientations (vide infra). In any event, the greater probability of occurrence of hydrogen bonds in small cages (9-11% of the time) versus large cavities (2-3%) illustrates the more intimate coupling of guests’ motion with the lattice borne of a smaller cage radius and closer contact, as reflected by the overlap of rattling modes with the lattice translational DOS (cf. Figures 2 and 4). Scrutiny of the AIMD power spectrum for guest hydrogen atom H-S stretch and H-S-H bend modes (cf. Figure 5) reveals the pertinent frequencies, which are provided in Table 1, along with the experimental data of Richardson et al.3 at 10 K. It can be seen that these frequencies are essentially reproduced by AIMD, with relatively little temperature dependence; the frequencies at 150 K are very marginally higher than at 300 K. The limited temperature dependence is reassuring in terms of comparison to experimental data measured at 10 K. It can be seen that the frequencies are somewhat (20-40 cm-1) higher in the small cages, owing to more intimate contacts with water (vide supra for a discussion of more frequent, transient hydrogen bonding in small cavities) affecting intramolecular vibrations in H2S. 6229

dx.doi.org/10.1021/jp111485w |J. Phys. Chem. A 2011, 115, 6226–6232

The Journal of Physical Chemistry A

ARTICLE

Table 1. Comparison of the Frequencies (cm-1) of the S-H Stretch and H-S-H Bending Modes from AIMD (cf. Figure 5) with Experimental Data of Richardson et al. at 10 K (Ref 3)

a

cage/vibration

150 K

300 K

expt

small, S-H

2604

2582, 2618

large, S-H

2591

2574

2591, 2604, 2616 2532, 2548, 2556, 2570

small, H-S-H

1169

1136, 1176

1176a

large, H-S-H

1159

1156

1176a

Not attributed specifically to this cage type.

Figure 7. Normalized autocorrelations of K4 and K6, with the same notation as in Figure 6 at (a) 150 and (b) 300 K. The short-time exponential decay is evident in the inset using logarithmic axes; the decay times are specified in Table 2.

Figure 6. Evolution of estimates of time-averaged ÆK4æ and ÆK6æ from AIMD at (a) 150 and (b) 300 K.

To assess the possibility of preferred guest orientations, the time evolution of the running averages for the ÆK4æ and ÆK6æ Kubic harmonics is depicted for AIMD in Figure 6. As mentioned previously, a nonzero value for ÆK4æ indicates preferential alignment relative to the [100]-direction.20 It can be seen that ÆK4æ is 0.23 and 0.09 in the small and large cavities, respectively, at 150 K (cf. Figure 6a), while there is much less evidence of preferential alignment at 300 K (cf. Figure 6b). Classical MD shows markedly less preferential alignment at all temperatures (75, 150, and 300 K), with, for instance, ÆK4æ values of 0.045 and 0.01 in the small and large cavities, respectively, at 75 K (data not shown), and similarly small ÆK6æ values (like AIMD, in that respect, cf. Figure 6). ACFs of K4 and K6 are provided in Figure 7 for AIMD, and the initial decay also depicted on logarithmic axes to emphasize the approximation to exponential decay; relaxation times τ from corresponding fits are provided in Table 2, along with results for (classical) MD. In any event, the decay is very rapid, with classical MD and AIMD in accord in this respect, with τ values in the sub-0.04 ps range; this indicates that the Kubic harmonics vary very rapidly; the convergence in running-average estimates of ÆK4æ and ÆK6æ is evident from Figure 6.

Figure 8. (a) H2S dipole-vector autocorrelation functions in the large (“L”) and small (“S”) cages at 100 and 300 K. The short-time exponential decay is evident in the inset using logarithmic axes; the decay times are specified in Table 2.

Prior to further detailed discussion on the origin of this preferential alignment at 150 K, particularly in the small cages, it is useful to examine rotational motion of the guests. ACF’s of μ and Cy2 are shown in Figures 8 and 9 for AIMD, with relaxation times τ specified in Table 2 for classical MD and AIMD. The agreements on the predicted relaxation times from the classical potential and quantum calculations for the decay of K4 and K6 are very close. In comparison to 300 K, there are significant differences in the relaxation times for μ and Cy2 at 150 K with longer τs predicted from AIMD calculations. This observation indicates the importance of quantum effects in the enhancement of stronger intermolecular H2S-water interactions at low temperature (vide supra). In the various ACFs of Figures 7-9, it is worth noting that a transient peak in the 300 K AIMD results 6230

dx.doi.org/10.1021/jp111485w |J. Phys. Chem. A 2011, 115, 6226–6232

The Journal of Physical Chemistry A

ARTICLE

appears after around 2.7 ps in both cavity types; this appears to be an artifact due to constructive/destructive interference patterns at this approximate time-scale in length dilations in the variablecell NPT box, the amplitude of which is less evident at 150 K (cf. Figure 1a). τ for the dipole-vector ACF is essentially doubled from around 0.25 to 0.5 ps for AIMD in reducing temperature from 300 to 150 K, with similar values for both cage types; this observation is consistent for both classical MD and AIMD, but the classical MD results exhibit substantially less temperature dependence, increasing only from around 0.24-0.30 to 0.280.34 ps (mirroring the small Kubic harmonics values at 75, 150, and 300 K, unlike AIMD; cf. Figure 6). Given that the ACF of Cy2 reflects rotational motion about the principal axis of rotation in H2S, it is perhaps more relevant to examine this in greater detail. τ in this case is larger for motion in the small cages at 150 K, at 0.82 ps, vis-a-vis the large cage at 0.46 ps, and the AIMD result for the small cavity is substantially larger than (classical) MD, that is, 0.82 versus 0.47 ps, while that for the large cage is also somewhat larger than its classical counterpart, that is, 0.46 versus 0.34 ps (cf. Table 2). Allowing for 75 K as being a more relevant temperature for comparison with 150 K in AIMD (in terms of ratio to melting temperature, vide supra), the respective comparisons are still 0.82 versus 0.54 ps and 0.46 versus 0.36 ps. τ in the large cavity roughly doubles from 0.21 to 0.46 ps for AIMD upon reduction in temperature, while it increases by a more dramatic factor of 3.5 in the small cavity (0.23 to 0.82 ps). Again, as was the case for the dipole-vector ACF, classical MD displays a much weaker temperature dependence, with the corresponding factors of increase

Figure 9. (a) Normalized Cy2 rotational autocorrelation functions in the large (“L”) and small (“S”) cages at 100 and 300 K. The short-time exponential decay is evident in the inset using logarithmic axes; the decay times are specified in Table 2.

being 1.4 and 1.75 for the large and small cages, respectively (0.24 to 0.34 and 0.27 to 0.47 ps). Therefore, it is evident from these observations of rotational time scales and evidence preferred orientations that AIMD at lower temperatures, like 150 K (but not as low as 100 K, vide supra), is preferable to classical MD, and that classical MD at low temperature (even 75 K) does not capture the more dramatic slow-down in rotational motion (as opposed to translational rattling, cf. Figure 4) of H2S, particularly in the small cages where there is greater evidence of preferred orientation (cf. Figure 6a). To discuss the underlying origins of preferential H2S orientation at 150 K revealed by AIMD, especially in small cavities, the angular orientation of the (unit) dipole-vector, μ/|μ|, vis-avis the laboratory (X-ray crystallography) x-, y-, z-axes were studied via transformation to spherical polar coordinates, that is,

Figure 10. Representative snapshot of dipolar orientation of H2S in a small cavity at 150 K from AIMD (with the dipole vector denoted by an arrow) toward outward (dipole)-pointing water molecules to optimize the dipolar interaction (in this instance, water molecules numbered 25, 18, 7 and 9). Water-water hydrogen bonds are denoted by dotted lines, while a hydrogen bond with the H2S hydrogen as donor and the oxygen atoms in water molecule number 9 is also shown in green and is 2.15 Å long.

Table 2. Exponential Decay Times (ps) for Autocorrelation Functions of K4 and K6, the H2S Dipole-Vector and Also Cy2 from Classical MD at 75 K and from AIMD and Classical MD at 150 and 300 K (AIMD, MD; cf. Figures 7-9), with the AIMD Results Specified First in Each Case (150 and 300 K) ACF/cage

75 K (MD)

150 K

300 K

K4, small

0.036 ( 0.0048

0.039 ( 0.0052, 0.032 ( 0.0039

0.030 ( 0.0049, 0.028 ( 0.0036

large

0.019 ( 0.0032

0.020 ( 0.0034, 0.017 ( 0.0028

0.015 ( 0.0030, 0.014 ( 0.0024

K6, small

0.021 ( 0.0037

0.014 ( 0.0027, 0.017 ( 0.0033

0.010 ( 0.0023, 0.012 ( 0.0027

large

0.018 ( 0.0029

0.012 ( 0.0028, 0.016 ( 0.0032

0.009 ( 0.0018, 0.014 ( 0.0020

μ, small

0.34 ( 0.037

0.44 ( 0.062, 0.31 ( 0.051

0.24 ( 0.042, 0.28 ( 0.043

large

0.30 ( 0.029

0.52 ( 0.075, 0.27 ( 0.038

0.27 ( 0.038, 0.24 ( 0.035

Cy2, small

0.54 ( 0.076

0.82 ( 0.14, 0.47 ( 0.062

0.23 ( 0.047, 0.27 ( 0.045

large

0.36 ( 0.059

0.46 ( 0.081, 0.34 ( 0.58

0.21 ( 0.038, 0.24 ( 0.040

6231

dx.doi.org/10.1021/jp111485w |J. Phys. Chem. A 2011, 115, 6226–6232

The Journal of Physical Chemistry A θ = cos-1(μz/|μ|) and j = tan-1(μy/μx), and scrutiny of the temporal variations of j and θ was undertaken. The timeaverages in the small cavity at 150 K of |θ| and |j| were 90 and around 0, respectively, although there are in practice very large angular fluctuations, with the guest dipole-vector appearing to be oriented in certain configurations for a period of around a half-picosecond and spending intermediate times of roughly a quarter-picosecond in transition to the next dipolar alignment, that is, de facto “flitting” between these orientations. In addition to this, much computer graphics visualization was undertaken of snapshots corresponding to these more preferred H2S dipolar orientations; it was found that hydrogen bonding was present in almost 25% of these particular orientations, as compared to 11% in terms of the entire production trajectory (23 ps, i.e., 25 ps less 2 ps of relaxation). One such representative configuration is depicted in Figure 10 (which, as it happens, corresponds roughly in θ and j to their time-averages), and this example also has a guest-water hydrogen bond (shown in green). From the study of the preferred dipolar orientations (shown by an arrow in the example of Figure 10), it appears that these preferred orientations, to which the H2S guests “flit” their dipoles via rotation (the Cy2 relaxation time of 0.82 ps in the small cage overlapping with a “residence” time of about 0.5 ps and a “transition” time of 0.20.3 ps) is motivated by the optimization (maximization) of the energetically favorable dipolar interaction with water molecules consisting of outward-pointing protons from the cage, that is, crude alignment of the guest dipole-vector with those of similarly oriented water dipoles: this may be seen with water molecules numbered 25, 18, 7, and 9 in Figure 10. In such orientations, which tend to align guest dipoles with those several water molecules, guest-water hydrogen bonds tend to be more likely; for instance, the hydrogen bond in the example of Figure 10 survived for almost 70 fs before rotating away from that preferred dipolar configuration (and then flitting back several more times during the course of the simulation).

’ CONCLUSIONS Equilibrium ab initio (AI) and classical constant pressure-constant temperature molecular dynamics (MD) simulations have been performed to investigate the dynamical properties in (type I) hydrogen sulphide hydrate at 150 and 300 K and 1 bar, with particularly emphasis on guest motion. The rattling motions of the guests in the large and small cavities were around 45 and 75-80 cm-1, respectively, from AIMD, with the corresponding classical MD modes being 10-12 cm-1 less at 150 K and around 5 cm-1 lower at 300 K. The rattling motion in the small cavities overlapped somewhat with the translational motion of the host lattice (with modes at circa 85 and 110 cm-1), due to closer contact. Experimental data for the guests’ intramolecular stretch and bending modes was also replicated by AIMD. Consideration of Kubic harmonics for the guest molecules from AIMD revealed that a preferred orientation of the dipole-vector exists at 150 K in both the small and large cavities but is markedly more significant for the small cavity, while there is no preferred orientation at 300 K. In comparison, classical MD did not reveal any preferred orientation at either temperature or at 75 K where there would be a closer correspondence to the AIMD case of 150 K (vide supra); classical MD may therefore be considered less reliable, at least at 75-150 K. The Cy2 time scales of 0.5 ps (large) and 0.8 ps (small) tended to correspond with “flitting” time scales to these preferred orientations, which were attributable to optimization of the

ARTICLE

dipolar interaction between the guest and outward-pointing water dipoles in the cavity. Although the exact, quantitative results for these observations might vary with a more robust incorporation of dispersion into DFT and the shortcomings of the PBE functional, for example, for the ice Ih melting point, and more refined criteria for the existence of guest-host hydrogen bonds, it has been shown that BOMD can offer a qualitatively accurate picture of guest dynamics in H2S hydrate. It should be borne in mind that other species of hydrate, such as those of xenon and propane, also contain distinct guest-host interactions, so it may be worthwhile applying BOMD to a detailed comparison of guest motion therein, in comparison to both classical MD and available experimental data.

’ AUTHOR INFORMATION Corresponding Author

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

’ ACKNOWLEDGMENT N.J.E. and J.S.T. thank the Ireland-Canada University Foundation and the Royal Irish Academy for research visit funding and the High Performance Computing Virtual Laboratory for the provision of computational resources. ’ REFERENCES (1) Makogon, Y. F. Hydrates of Hydrocarbons; PennWell Books: Tulsa, OK, 1997. (2) Sloan, E. D.; Koh, C. A. Clathrate Hydrates of Natural Gases, 3rd rev. ed.; CRC Press, Taylor & Francis: New York, 2007. (3) Richardson, H. H.; Wooldridge, P. J.; Devlin, J. P. J. Chem. Phys. 1985, 83, 4387. (4) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865. (5) Sanchez-Portal, D.; Artacho, E.; Soler, J. M. Solid State Commun. 1995, 95, 685. (6) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (7) Jorgensen, W. L. J. Phys. Chem. 1986, 90, 6379. (8) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, 1987. (9) de Leeuw, S. W.; Perram, J. W.; Smith, E. R. Proc. R. Soc. London, Ser. A 1980, 373, 27. (10) Andersen, H. C. J. Comput. Phys. 1983, 52, 24. (11) Hoover, W. G. Phys. Rev. A 1985, 31, 1695. (12) Melchionna, S.; Ciccotti, G.; Holian, B. L. Mol. Phys. 1993, 78, 533. (13) McMullan, R. K.; Jeffrey, G. A. J. Chem. Phys. 1965, 42, 2725. (14) Bernal, J. D.; Fowler, R. H. J. Chem. Phys. 1933, 1, 515. (15) Rahman, A.; Stillinger, F. H. J. Chem. Phys. 1972, 57, 4009. (16) Lindberg, G. E.; Wang, F. J. Phys. Chem. B 2008, 112, 6436 and references therein. (17) Yoo, S.; Zeng, X. C.; Xantheas, S. S. J. Chem. Phys. 2009, 130, 221102. (18) Hester, K. C.; Huo, Z.; Ballard, A. L.; Koh, C. A.; Miller, K. T.; Sloan, E. D. J. Phys. Chem. B 2007, 111, 8830 and references therein. (19) English, N. J.; MacElroy, J. M. D. J. Comput. Chem. 2003, 24, 1569 and references therein. (20) Tse, J. S.; Klein, M. L.; McDonald, I. R. J. Chem. Phys. 1984, 81, 6146 and references therein. (21) English, N. J.; MacElroy, J. M. D. Mol. Phys. 2002, 100, 3753. (22) Impey, R. W.; Madden, P. A.; McDonald, I. R. Mol. Phys. 1982, 46, 513. (23) Luzar, A.; Chandler, D. Phys. Rev. Lett. 1996, 76, 928. 6232

dx.doi.org/10.1021/jp111485w |J. Phys. Chem. A 2011, 115, 6226–6232