Vibrational Modes of Hydrogen Hydrates: A First-Principles Molecular

Feb 8, 2017 - School of Chemical and Bioprocess Engineering, University College Dublin, Belfield, Dublin 4, Ireland. ‡ Istituto dei Sistemi Compless...
0 downloads 0 Views 1MB Size
Article pubs.acs.org/JPCC

Vibrational Modes of Hydrogen Hydrates: A First-Principles Molecular Dynamics and Raman Spectra Study Zdenek Futera,*,† Milva Celli,‡ Leonardo del Rosso,‡,§ Christian J. Burnham,† Lorenzo Ulivi,*,‡ and Niall J. English*,† †

School of Chemical and Bioprocess Engineering, University College Dublin, Belfield, Dublin 4, Ireland Istituto dei Sistemi Complessi (CNR-ISC), Consiglio Nazionale delle Ricerche, Via Madonna del piano 10, I-50019 Sesto Fiorentino, Italy § Dipartimento di Fisica e Astronomia, Università degli Studi di Firenze, via Sansone 1, I-50019 Sesto Fiorentino, Italy ‡

S Supporting Information *

ABSTRACT: We have employed classically propagated molecular dynamics (MD), within the framework of density functional theory (DFT), to calculate vibrational spectral band of molecular hydrogen trapped in clathrate hydrate, with large-cage occupancy from 1 to 4, at ∼260 K and ∼2 kbar. The predicted vibrations, obtained by applying a state-of-the-art generalized gradient approximation (GGA) functional with nonlocal correlation (VdW-DF), reproduce satisfactorily our own accurate Raman spectra (at the same temperature and pressure conditions). We decomposed the MD-sampled vibrational band to individual peaks and assigned them to the vibration of H2 molecules enclosed in small and large cages of SII hydrate. By summing the resulting spectral bands, we have demonstrated that the measured spectral response is a complex composition of signals originating from H2 molecules experiencing different local, intracage environments.



INTRODUCTION Clathrate hydrates are nonstoichiometric crystalline inclusion compounds in which a water host lattice encages small guest atoms or molecules in cavities,1,2 and they have potential utility as a hydrogen-storage medium.3−15 Mao et al. found that at 200 MPa hydrogen and water mixtures crystallize in the form of SII clathrate, with maximal hydrogen storage capacity of 5.3 wt % at 300 MPa and 250 K,16 which is close to the US DOE’s 2020 5.5 wt % target.16 SII clathrate’s host lattice is composed of a combination of relatively large hexadecahedral 51264 cages and smaller 512 dodecahedral cages in a ratio of 1:2. H2 molecules may occupy both types of cages, and can “hop” between cavities,4−9,11−18 with multiple occupancy of the large cage. Hydrogen transport is “anomalous” in the sense that intracage diffusivity is “superconducting” prior to developing essentially Fickian intercage migration.17 In any event, modeling of the dynamics of hydrogen clathrate hydrates has been approached in the past taking into account quantum aspects of the H2 motion, in terms of both molecular rotation and translation of the center of mass. This approach has unveiled important features, such as the discrete spectrum of the translational states, and their coupling with the rotational ones, and has been able to reproduce quantitatively, with good approximation, the spectra of the H2 molecules, measured by either Raman scattering19 and neutron inelastic scattering.20−25 This quantum-dynamics calculation is based on the firstprinciples Schrödinger’s equation. The essential drawback of © XXXX American Chemical Society

this approach is that only the simple case of a single H2 molecule in one cage can be computed rigorously at the present time. The theoretical study of the quantum dynamics of two or more molecules trapped together in one cage, as is the case in 51264 “large” cages of SII hydrogen clathrates, presents more difficulties and has instead been approached by different quantum-dynamics methods, such as diffusion Monte Carlo,26 path-integral molecular dynamics,27 and nuclear configuration interaction.28,29 In any case, only a few of these methods give access to the higher energy levels and detailed spectroscopic properties, and in all cases, the lattice of water molecules is considered rigid, a semiempirical potential model has to be used, and the H2 intramolecular vibrational states are not considered. In order to study intramolecular vibrational states in clathrates rigorously, both of the guest species and potentially of the water lattice, Raman spectroscopy represents an important tool for the study, since the guest-molecule excitations and the lattice vibration give rise to well-identifiable bands. This technique has been applied efficiently by several authors to hydrogen clathrates prepared ex situ and examined at 20−77 K.30,31 At these low temperatures, it has been possible to resolve and assign the different components of vibrational band, Received: November 2, 2016 Revised: January 25, 2017

A

DOI: 10.1021/acs.jpcc.6b11029 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

optical window (made of artificial diamond) to allow for the collection of the Raman signal in a backscattering configuration.35 The sample sits on a horizontal container, inside the cell, with the top surface exposed to laser light. The cell is connected with a closed-cycle cryo-refrigerator under vacuum. The Raman spectra were collected at subsequent times, at a constant temperature and pressure, during an overnight run. The Raman signal from hydrogen molecules in the different cages grows at the same time as ice transforms to the clathrate phase. The collection of each spectrum lasted half an hour. The second experiment studies clathrate formation from a liquid solution of water and hydrogen. The high-pressure Raman cell has been built to contain two phases, a liquid and a gaseous one, separated by gravity. The cell has an almost cubic shape and optical access through three optical artificial diamond windows on three vertical faces. The scattering angle is 90°. A steel capillary, linked to high pressure circuit, allows the entrance of bidistilled water. A second steel capillary, entering the cell from the top, allows the gas supply. The cell is described in greater detail in ref 36. Water was introduced into the cell at nearly atmospheric pressure and at ambient temperature, filling the cell so as to position the water meniscus just above the optical diamond windows. The temperature was subsequently lowered down to 280 K, the hydrogen pressure was raised to 1990 bar, and a sufficiently long time has been waited to allow for complete diffusion of hydrogen gas in the liquid sample. Spectra were collected at subsequent times during clathrate growth. The increasing intensity of the bands due to trapped hydrogen molecules testifies to the increase of volume of solid clathrate particles, initially floating in the liquid. In both experiments, the pressure needed for hydrogen clathrate synthesis was obtained by the use of a 3 kbar membrane compressor (Nova-Swiss) and measured with a nominal accuracy of 0.2% of full scale value (3.5 kbar). The Raman spectra were obtained by using the green line (λ = 514.5 nm) of an argon ion laser for excitation, typically operating at about 35 mW on the sample. The beam image was focused on the entrance slit of a Spex spectrometer (equipped with 1800 grooves/mm holographic gratings). The spectra were recorded by a thermoelectrically cooled CCD detector (Andor DU401ABV). The overall spectral resolution was 0.4 cm−1.

arising from different hosting cages, different large-cage occupation, and different spin isomers of guest molecule (ortho and para hydrogen).31,32 At higher temperatures, e.g., close to that of formation of hydrogen clathrate, the enhanced guest−host interactions complicate the band structure and prevent an immediate assignment. Theoretically, hydrogen vibrational frequencies in clathrate cages have been estimated by static calculations based on cluster models.33,34 Although the calculations predicted a redshift of H2 stretching vibration in the clathrate comparing to its gas-phase value, such calculations are limited in several ways. First, the condensed-matter environment of a clathrate cage is completely neglected in cluster models, though it might have only a minor effect on hydrogen vibrations. In addition, only optimized structures corresponding to total energy minimum of a model are analyzed, which corresponds to an unrealistic zerotemperature scenario. Although quantal effects are important at low temperature, these are not considered in such studies. Therefore, for more relevant interpretation of the spectra, molecular dynamics (MD) applied to an extended system should be used to properly sample clathrate structures at finite temperate. In the present study, we have measured vibrational-Raman spectra of hydrogen molecules in simple clathrates, during synthesis in two different conditions, namely when clathrate is forming from ice Ih, at −10 °C and 1510 bar, and when it is forming from a liquid hydrogen−water solution at −15 °C and 1990 bar; the underlying experiments are described in refs 35 and 36. Here, we combine these vibrational-Raman spectra with first-principles molecular dynamics (MD) in order to study the atomistic details of hydrogen-molecule and water-lattice intramolecular vibrational modes and characterize optimal functionals from density functional theory (DFT) to adopt. Indeed, although the MD is classical via Newton’s equations, this matters relatively little for a reasonably accurate treatment of intramolecular vibrations. However, admittedly, with classical MD, there may well be a near-constant “shift” in spectral bands vis-à-vis a more computationally prohibitive quantum treatment of stretch-vibrational dynamics (e.g., a direct solution of timedependent Schrödinger’s equation for H−H stretch). However, path-integral dynamics methods fail to account rigorously for dynamical properties and so cannot be relied upon fully to compute intramolecular vibrations. In this sense, the use of classical MD here is quite a reasonable approach, in that we do not suppose to achieve exact quantitative agreement with Raman spectra, and some relatively modest level of spectralband shift is not unacceptable. What is crucial, however, is to use an accurate functional and treatment of dispersion interactions, particularly to mimic the relatively high pressure conditions employed experimentally,36 given that available empirical potentials for guests (such as hydrogen) typically fail to adjust and polarize to their (intracage) environment and readjust their vibrational modes accordingly. In the present study, we seek to apply MD to mimic Raman-vibrational modes as accurately as possible, bearing in mind computational expediency.



COMPUTATIONAL DETAILS Born−Oppenheimer (BO) MD simulations under periodic boundary conditions (PBC) were carried out within the framework of density functional theory (DFT) in generalizedgradient approximation (GGA). We employed DRSLL exchange-correlation functional37,38 which involves van der Waals (vdW) interaction via explicit nonlocal correlation corrections. This functional has been shown recently to perform very well for liquid water.39−44 A time step of 0.2 fs was used in the canonical (NVT) ensemble, at 258 K in conjunction with a Nosé−Hoover thermostat with fictitious mass 120 Ry.fs2.45,46 Charge density in each MD step was converged to accuracy 10−4 au. The calculations were performed in DZP basis set using the SIESTA program.47 A cubic unit cell of a type II clathrate hydrate was studied comprising of 136 water molecules forming structure of 16 small cages and 8 large cages, with vanishingly small dipole.4,17,18 One hydrogen molecule resides in each small cage, while large-cage occupations from 1 to 4 were explored here. The lattice constant of the cubic system was calibrated



EXPERIMENTAL METHODS The first experiment investigates clathrate formation from solid ice exposed to hydrogen gas. In this case we have employed a specially designed cell, made of beryllium−copper alloy, which could withstand a hydrogen gas pressure up to 3 kbar at all temperatures ranging from 300 to 20 K, and equipped with one B

DOI: 10.1021/acs.jpcc.6b11029 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C carefully to 16.6 Å, which yields a pressure of about 2 kbar. Although the optimized constant is by ∼2.6% smaller than the experimental value of 17.047 Å,3 which indicates certain inaccuracies in applied DFT potential under high-pressure conditions, we employed this cell dimension in the simulation as it corresponds the desired pressure for the given state-of-theart DFT functional. For comparison, we also investigated hydrogen gas dissolved in liquid water, for which the cubic cell of 15.5 Å side was filled by 152 water molecules and 4 H2 molecules, which roughly corresponds to the 2.5% molar concentration of the experiment,36 estimated from the known hydrogen solubility in water,48 governed essentially by Henry’s law.49 All models were equilibrated carefully before the NVT production runs. The vibrational spectral band of H2 molecules was calculated by Fourier transforming the H−H bond length (dHH) autocorrelation function C(t) = ⟨dHH (t 0 )dHH (t 0 + t)⟩ accumulated from ca. 5 ps long trajectories I(ω) ∝

∫ ⟨dHH(t0)dHH(t0 + t )⟩e−iωt dt

Figure 2. Hydrogen vibrational Raman spectra recorded during hydrogen clathrate synthesis from liquid solution at −15 °C and 1990 bar. The different colors correspond to different spectra acquisitions performed during clathrate growth (cf. ref 32).

well with that measured during the synthesis from the solid phase. The relative intensity of the bands is not constant during clathrate growth, demonstrating a different growth kinetic for small and large cages.36 The calculated vibrational band of the H2 molecules trapped in both small and large (quadruply occupied) clathrate cages is shown in Figure 3, together with the spectral profile of H2 dissolved in pressurized liquid water. The clathrate band is redshifted with respect to the water one, with absolute positions of the computed vibrational peaks shifted by ca. 220 cm−1 to higher frequencies vis-à-vis experimental data. This is consistent with results published in ref 34 which were computed at the B97-D/cc-pVTZ level of theory in optimized cluster model. Similar red-shifted peak positions were obtained also in ref 33 although they were rescaled in that work to fit the experimental position. The reference value has been obtained by calculating the vibrational frequency of the H2 molecule in the gas phase with the same DRSLL method (via BO-MD). The origin of this shift is not caused by the DFT functional (vide inf ra), but rather by classical time evolution assumed, ipso facto, by BOMD. To justify the DRSLL functional, we performed a set of benchmark calculations of single-molecule H2 vibration at different levels of theory. The results are summarized in Table 1. First, we calculated the H2 frequency in vacuum using static quantum-chemistry approach based on harmonic approximation where the frequencies are obtained from Hessian matrix of total-energy second derivatives; these calculations were performed by Gaussian software.51 The expected frequency value for this kind of calculation is 4400 cm−1, which is again similarly shifted with respect to the experimentally observed H2 vibration 4160 cm−1 (disregarding the ortho- and para-H2 effects for simplicity), taking into account corrections for anharmonicity.52 This value was obtained almost exactly by coupled-cluster theory at the CCSD(T)/Aug-CC-pVTZ level. The pure Hartree−Fock (HF)/Aug-CC-pVQZ method, without dynamical electronic correlations, yields a higher frequency by ca. 180 cm−1, and this overestimation is only partly corrected by second-order perturbation theory (specifically, MP2/AugCC-pVTZ). On the other hand, pure GGA DFT functionals (all checked with Aug-CC-pVQZ basis set) underestimate the expected value by 60−90 cm−1. Clearly, this is caused by the complete lack of exact HF exchange; therefore, the hybrid B3LYP functional,53−55 where this term is partly included,

(1)

The autocorrelation function itself was computed on the basis of Wiener−Khintchine relations,50 employing the fast Fourier transform (FFT) algorithm with efficient O(n log n) scaling 1 ⟨f (t0)f (t0 + t )⟩ = |f (t ′)e−iωt ′ dt ′|2 eiωt dω (2) 2π

∫∫

In this way, the calculated vibrational spectrum is decoupled from the translational and rotational modes, and it is directly comparable with the measured Raman spectra.



RESULTS AND DISCUSSION Raman spectra recorded in the first experiment at two subsequent times after gas injection are shown in Figure 1,

Figure 1. Hydrogen vibrational Raman spectra recorded at different times (indicated in the legend) after gas insertion in the cell, during hydrogen clathrate synthesis from ice Ih at −10 °C and 1510 bar (cf. ref 31).

after the subtraction of the contribution of the gas present in the cell, as described in ref 35. A series of spectra obtained in the second experiment at different times after the onset of clathrate formation are reported in Figure 2, after the subtraction of gas or liquid contribution, as explained in ref 36. The recorded intensity, due to the H2 molecules trapped in the small and large cages of the clathrate, grows with time as the cell gets filled with the solid clathrates at the expenses of the liquid portion of the sample. The position of the bands agrees C

DOI: 10.1021/acs.jpcc.6b11029 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 3. Calculated vibration spectra of H2 in clathrate hydrate with quadruple-occupied large cages (solid red curve). The position of vibration band of H2 dissolved in water under same temperature and pressure is shown for comparison (dashed green curve).

Table 1. Calculated Vibrational Frequency [cm−1] of the H2 Molecule Using Harmonic Approximation Approach Based on Second Derivatives of Total Energy and Molecular Dynamics Approach Based on Fourier Transform of Velocity Autocorrelation Functiona harmonic approx gas phase water molecular dynamics gas phase water a

HF

MP2

4581 4557 HF

4516 4495 MP2

4571 4550

4509 4489

CCSD

PBE

BLYP

B3LYP

DRSLL VdW

4399 4380 CCSD

4311 4294 PBE

4342 4326 BLYP

4415 4397 B3LYP

− − DRSLL VdW

4387 4367

4298 4275

4326 4316

4408 4387

4382 4362

Explicit water in periodic boundary conditions was used in DRSLL-VdW case while implicit solvent (COSMO) was employed in other cases.

rather than the DFT functional or electronic-structure treatment per se. To analyze the origin of the vibrational peaks observed in the hydrate H2 spectra, we decomposed the band to contributions of H2 trapped in small and large clathrate cages, respectively. The decomposition is shown in Figure 4, where clathrates with large-cage occupancy of 1, 2, 3, and 4 are compared. The vibrational bands shown in Figure 4 were constructed from Gaussian fits of their single-molecular contributions. The original spectrum obtained directly by Fourier transforming the H−H bond-length autocorrelation function is shown in Figure S1 (Supporting Information). For all occupations, the low-frequency part of the vibrational band originates from H2, trapped in small cages while the hydrogen molecules accommodated in the large cages contribute to the higher frequency peaks. This demonstrates clearly that the molecular environment unambiguously affects the spectral pattern. The ratio of the peak areas corresponds to number of the H2 molecules located in the large/small cages, i.e., 0.5, 1, 1.5, and 2 for systems with large-cage occupations of 1, 2, 3, and 4, respectively. However, the large-cage occupancy is not the same in the whole sample in experimental measurements, especially during the clathrate growth when the cages are filled gradually by hydrogen, as is shown nicely in Figures 1 and 2. To explore this situation on the basis of calculated results, we constructed several composed spectra, which are shown Figure 5. As the fraction of large cages occupied by specific number of H2 molecules is not known experimentally, we summed the

provides the very accurate results. The effect of solvation, here simulated implicitly by COSMO (ε = 78.4),56 leads to frequency lowering by 18−20 cm−1 for all the benchmarked methods. The solvation effect on the H−H stretching frequency is therefore by 1 order of magnitude lower than the frequency change induced by considering the electron correlation in the calculations. Next, we calculated the same single-molecular H2 vibrations sampled from MD via Fourier transformation (cf. eqs 1 and 2) to obtain frequencies of vibrational modes. We performed NVE dynamics in Gaussian (for the same methods discussed in the previous paragraph) and in SIESTA (for DRSLL), where the initial temperature was set to 300 K. The results are summarized in the second part of Table 1. Comparing to the static harmonic approach, the MD frequencies are lower by ca. 10 cm−1 in all cases. This is a surprisingly small effect regarding the experimentally observed 4160 cm−1 value (240 cm−1 lower than the harmonic 4400 cm−1). All of the methods therefore predict a rather harmonic shape of the H2 potential, and the anharmonicities are practically negligible. Solvation effects, investigated by simulation of H2 in explicit liquid water environment, yield a change of ca. 20 cm−1, which is consistent with the implicit solvent models. The functional of interest, DRSLL, leads to a frequency of 4362 cm−1 in water, which is practically identical with a highly accurate CCSD value of 4367 cm−1. Therefore, we can conclude that the shift of calculated spectra is caused by classical (Verlet-type) dynamics, based on second-order approximation of Newton’s equations of motions, D

DOI: 10.1021/acs.jpcc.6b11029 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

(necessarily) limited number of molecules used in the simulation. To smooth the spectra and show the band decomposition more schematically, we fitted the small- and large-cage contributions independently by Gaussian curves and reconstructed the whole vibrational band by their sum with weights corresponding to number of H2 molecules. The fitting is shown in Figure S2. In any event, the Gaussians’ average frequencies are listed in Table 2. Their widths range from 26 to Table 2. Average Frequencies of Fitted Gaussians Fitted to BO-MD Data for Small and Large Cages (SC and LC) (Cf. Fits in Figure 4/Figure S2) as a Function of LC Occupation LC occupation

SC frequency (cm−1)

LC frequency (cm−1)

single double triple quadruple

4258.5/4250.7 4251.8/4239.5 4248.3/4231.8 4245.7/4238.5

4285.2/4275.8 4303.6/4295.8 4310.2/4301.5 4332.6/4311.5

35 cm−1, overestimating the experimental ones, which are of the order of 6−10 cm−1. From these values, we note that the average frequency of the molecules in the large cages increases with the occupation number, as expected from experimental data. We see also a quite apparent nonmonotonic and somewhat unexpected dependence on the same quantity of the frequency of the molecules in the small cages. Returning to the Gaussian fitting per se, reference to Figure S2 shows that the small- and large-cage bands overlap mainly in the singly occupied system; however, they become separated when there is more than one H2 molecule in the large cages. This approach cannot describe the triple-peak structure of the experimental spectra (due probably to the larger widths of the computed Gaussians), though, even when we constructed the composed bands (cf. Figure S3). More realistically, we fitted by Gaussian curves to all individual H2 peaks and reconstructed the spectra (cf. Figure 4) and their compositions from them (cf. Figure S4): despite this, the triple-peak structure was not revealed. This indicates that the splitting of single-molecular H2 peaks as a result of the intermolecular H2−H2 interaction within the cages cannot be neglected by single-Gaussian fitting, and the experimental spectra reflect closely these weak interactions in the clathrate structures.

Figure 4. Fits of calculated vibrational bands of SII clathrate hydrates with different occupation of large cages by H2 molecules. Small-cage H2 contributions are shown in blue, large cage H2 contributions in red, and their sum is drawn by solid black line. The small-cage and largecage bands were constructed from independent Gaussian fits of singlemolecular H2 peaks. Areas of the bands correspond to number of the H2 molecules.



CONCLUSIONS The vibrational spectra of hydrogen molecules enclosed in the clathrate hydrate were measured by Raman spectroscopy. The obtained vibrational bands were then calculated by firstprinciples MD approach using the DFT. We applied the state-of-the-art vdW-DF functional with nonlocal correlation correction (DRSLL), known to perform well for liquid water. In essence, the calculated vibrational band reproduced qualitatively the measured spectra, and the individual peaks were identified based on the structural data. The multiple-peak structure of the observed band can be explained by presence of large cages with different H2 occupations in the measured hydrate sample. The peak decomposition indicates clearly that large cages with various H2 occupations are present in the sII hydrate during experimental measurements. We showed that the H2 molecules trapped in the small cages provides the lowfrequency peak in the spectral band, while the large-cage hydrogen contributes to the high-frequency peak. The mutual ratio of the peaks is related to the amount of H2 accommodated in small and large cages.

Figure 5. Compositions of calculated vibration spectra of H2 in clathrate hydrate with different occupancies of large cages. Small-cage H2 contributions are shown in blue, large cage H2 contributions in red, and their sum is drawn by solid black line. Areas of the small-cage and large-cage H2 bands are proportional to the number of H2 molecules contained in each one.

vibrational bands of different systems with equal weights as an approximation. Nevertheless, the composed bands agree qualitatively with the experimental data. All of the calculated spectra exhibit a series of narrow peaks, which are not reproducible from one simulation run to another, due to statistical noise of finite time and size sampling. In our view, these fine structures arise primarily as a result of the E

DOI: 10.1021/acs.jpcc.6b11029 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

(12) Lee, J. W.; Yedlapalli, P.; Lee, S. Prediction of Hydrogen Hydrate Equilibrium by Integrating Ab Initio Calculations with Statistical Thermodynamics. J. Phys. Chem. B 2006, 110, 2332−2337. (13) Sluiter, M. H. F; Adachi, H.; Belosludov, R. V.; Belosludov, V. R.; Kawazoe, Y. Ab Initio Study of Hydrogen Storage in Hydrogen Hydrate Clathrate. Mater. Trans. 2004, 45, 1452−1454. (14) Chattaraj, P. K.; Bandaru, S.; Mondal, S. Hydrogen Storage in Clathrate Hydrates. J. Phys. Chem. A 2011, 115, 187−193. (15) Willow, S. Y.; Xantheas, S. S. Enhancement of Hydrogen Storage Capacity in Hydrate Lattices. Chem. Phys. Lett. 2012, 525, 13− 18. (16) Mao, W. L.; Koh, C. A.; Sloan, E. D. Clathrate Hydrates Under Pressure. Phys. Today 2007, 60, 42−47. (17) Gorman, P. D.; English, N. J.; MacElroy, J. M. D. Dynamical Cage Behaviour and Hydrogen Migration in Hydrogen and HydrogenTetrahydrofuran Clathrate Hydrates. J. Chem. Phys. 2012, 136, 044506. (18) Burnham, C. J.; English, N. J. Free-Energy Calculations of the Intergace Hopping Barriers of Hydrogen Molecules in Clathrate Hydrates. J. Phys. Chem. C 2016, 120, 16561−16567. (19) Giannasi, A.; Celli, M.; Zoppi, M.; Moraldi, M.; Ulivi, L. Experimental and Theoretical Analysis of Rotational Raman Spectrum of Hydrogen Molecules in Clathrate Hydrates. J. Chem. Phys. 2011, 135, 054506. (20) Ulivi, L.; Celli, M.; Giannasi, A.; Ramirez-Cuesta, A. J.; Bull, D. J.; Zoppi, M. Quantum Rattling of Molecular Hydrogen in Clathrate Hydrate Nanocavities. Phys. Rev. B 2007, 76, 161401. (21) Xu, M.; Ulivi, L.; Celli, M.; Colognesi, D.; Bačić, Z. Quantum Calculation of Inelastic eutNron Scattering Spectra of a Hydrogen Molecule Inside a Nanoscale Cavity Based on Rigorous Treatment of the Coupled Translation-Rotation Dynamics. Phys. Rev. B 2011, 83, 241403. (22) Xu, M.; Ulivi, L.; Celli, M.; Colognesi, D.; Bačić, Z. Rigorous Quantum Treatment of Inelastic Neutron Scattering Spectra of a Heteronuclear Diatomic Molecule in a Nanocavity: HD in the Small Cage of Structure II Clathrate Hydrate. Chem. Phys. Lett. 2013, 563, 1−8. (23) Colognesi, D.; Celli, M.; Ulivi, L.; Xu, M.; Bačić, Z. Neutron Scattering Measurements and Computation of the Quantum Dynamics of Hydrogen Molecules Trapped in the Small and Large Cages of Clathrate Hydrates. J. Phys. Chem. A 2013, 117, 7314−7326. (24) Celli, M.; Powers, A.; Colognesi, D.; Xu, M.; Bačić, Z.; Ulivi, L. Experimental Inelastic Neutron Scattering Spectrum of Hydrogen Hexagonal Clathrate-Hydrate Compared with Rigorous Quantum Simulations. J. Chem. Phys. 2013, 139, 164507. (25) Colognesi, D.; Powers, A.; Celli, M.; Xu, M.; Bačić, Z.; Ulivi, L. The HD Molecule in Small and Medium Cages of Clathrate Hydrates: Quantum Dynamics Studied by Neutron Scattering Measurements and Computation. J. Chem. Phys. 2014, 141, 134501. (26) Sebastianelli, F.; Xu, M.; Bačić, Z. Quantum Dynamics of Small H2 and D2 Clusters in the Large Cage of Structure II Clathrate Hydrate: Energetics, Occupancy, and Vibrationally Averaged Cluster Structures. J. Chem. Phys. 2008, 129, 244706. (27) Witt, A.; Sebastianelli, F.; Tuckerman, M. E.; Bačić. Path Integral Molecular Dynamics Study of Small H2 Clusters in the Large Cage of Structure II Clathrate Hydrate: Temperature Dependence of Quantum Spatial Distributions. Z. J. Phys. Chem. C 2010, 114, 20775− 20782. (28) Felker, P. M. Fully Quantal Calculation of H2 TranslationRotation States in (H2)4@51264 Clathrate SII Inclusion Compounds. J. Chem. Phys. 2013, 138, 174306. (29) Felker, P. M. Fully Quantal Calculation of H2 TranslationRotation States in (p-H2)2@51264 Clathrate Hydrate Inclusion Compounds. J. Chem. Phys. 2014, 141, 184305. (30) Strobel, T. A.; Sloan, E. D.; Koh, C. A. Raman Spectroscopic Studies of Hydrogen Clathrate Hydrates. J. Chem. Phys. 2009, 130, 014506. (31) Giannasi, A.; Celli, M.; Ulivi, L.; Zoppi, M. Low Temperature Raman Spectra of Hydrogen in Simple and Binary Clathrate Hydrates. J. Chem. Phys. 2008, 129, 084705.

Although the applied classically propagated BO-MD leads to shift of the calculated vibrational frequencies, the DFT method used, which is still relatively fast compared to more rigorous many-body approaches, has been demonstrated as an effective tool for analyzing both structural and dynamical properties of clathrate systems.



ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.6b11029. Calculated vibrational H2 spectra in clathrate with different occupancies, single-Gaussian fits of calculated small- and large-cage vibrational bands, compositions of calculated fitted spectra with different occupancies (PDF)



AUTHOR INFORMATION

Corresponding Authors

*E-mail [email protected] (N.J.E.). *E-mail [email protected] (Z.F.). *E-mail [email protected] (L.U.). ORCID

Christian J. Burnham: 0000-0001-5574-4339 Lorenzo Ulivi: 0000-0002-3541-0419 Niall J. English: 0000-0002-8460-3540 Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS Z.F., C.J.B., and N.J.E. thank the Science Foundation Ireland for funding under grant SFI 15/ERC-13142. 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: USA, 2007. (3) Mao, W. L.; Mao, H. K.; Goncharov, A. F.; Struzhkin, V. V.; Guo, Q.; Hu, J.; Shu, J.; Hemley, R. J.; Somayazulu, M.; Zhao, Y. Hydrogen Clusters in Clathrate Hydrate. Science 2002, 297, 2247−2249. (4) Cao, H.; English, N. J.; MacElroy, J. M. D. Diffusive Hydrogen Inter-Cage Migration in Hydrogen and Hydrogen-Tetrahydrofuran Clathrate Hydrates. J. Chem. Phys. 2013, 138, 094507. (5) Alavi, S.; Ripmeester, J. A. Hydrogen-Gas Migration Through Clathrate Hydrate Cages. Angew. Chem., Int. Ed. 2007, 46, 6102−6105. (6) Frankcombe, T. J.; Kroes, G.-J. Molecular Dynamics Simulations of Type-sII Hydrogen Clathrate Hydrate Close to Equilibrium Conditions. J. Phys. Chem. C 2007, 111, 13044−13052. (7) Peters, B.; Zimmermann, N. E.; Beckham, G. T.; Tester, J. W.; Trout, B. L. Path Sampling Calculation of Methane Diffusivity in Natural Gas Hydrates From a Water-Vacancy Assisted Mechanism. J. Am. Chem. Soc. 2008, 130, 17342−17350. (8) English, N. J.; MacElroy, J. M. D. Perspectives on Molecular Simulation of Clathrate Hydrates: Progress, Prospects and Challenges. Chem. Eng. Sci. 2015, 121, 133−156. (9) Alavi, S.; Ripmeester, J. A.; Klug, D. D. Molecular-Dynamics Study of Structure II Hydrogen Clathrate. J. Chem. Phys. 2005, 123, 024507. (10) Mao, W. L.; Mao, H.-K. Hydrogen Storage in Molecular Compounds. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 708−710. (11) Patchkovskii, S.; Tse, J. S. Thermodynamic Stability of Hydrogen Clathrate. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 14645−14650. F

DOI: 10.1021/acs.jpcc.6b11029 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C (32) Zaghloul, M. A. S.; Celli, M.; Salem, N. M.; Elsheikh, S. M.; Ulivi, L. High Pressure Synthesis and In Situ Raman Spectroscopy of H2 and HD Clathrate Hydrates. J. Chem. Phys. 2012, 137, 164320. (33) Wang, J.; Lu, H.; Ripmeester, J. A. Raman Spectroscopy and Cage Occupancy of Hydrogen Clathrate from First Principles Calculations. J. Am. Chem. Soc. 2009, 131, 14132−14133. (34) Ramya, K. R.; Venkatnathan, A. Vibrational Raman Spectra of Hydrogen Clathrate Hydrates from Density Functional Theory. J. Chem. Phys. 2013, 138, 124305. (35) Celli, M.; Zoppi, M.; Zaghloul, M. A. S.; Ulivi, L. High Pressure Optical Cell for Synthesis and In Situ Raman Spectroscopy of Hydrogen Clathrate Hydrates. Rev. Sci. Instrum. 2012, 83, 113101. (36) del Rosso, L.; Celli, M.; Ulivi, L. Raman Measurements of Pure Hydrogen Clathrate Formation from a Supercooled Hydrogen-Water Solution. J. Phys. Chem. Lett. 2015, 6, 4309−4313. (37) Dion, M.; Rydberg, H.; Schroder, E.; Langreth, D. C.; Lundqvist, B. I. Van der Waals Density Functional for General Geometries. Phys. Rev. Lett. 2004, 92, 246401. (38) Roman-Perez, G.; Soler, J. M. Efficient Implementation of a van der Waals Density Functional: Application to Double-Wall Carbon Nanotubes. Phys. Rev. Lett. 2009, 103, 096102. (39) Wang, J.; Roman-Perez, G.; Soler, J. M.; Artacho, E.; FernandezSerra, M. V. Density, Structure, and Dynamics of Water: The Effect of van der Waals Interactions. J. Chem. Phys. 2011, 134, 024516. (40) Zhang, C.; Wu, J.; Galli, G.; Gygi, F. Structural and Vibrational Properties of Liquid Water from van der Waals Density Functionals. J. Chem. Theory Comput. 2011, 7, 3054−3061. (41) Corsetti, F.; Artacho, E.; Soler, J. M.; Alexandre, S. S.; Fernandez-Serra, M. V. Room Temperature Compressibility and Diffusivity of Liquid Water from First Principles. J. Chem. Phys. 2013, 139, 194502. (42) Bankura, A.; Karmakar, A.; Carnevale, V.; Chandra, A.; Klein, M. L. Structure, Dynamics, and Spectral Diffusion of Water from FirstPrinciples Molecular Dynamics. J. Phys. Chem. C 2014, 118, 29401− 29411. (43) English, N. J. Structural Properties of Liquid Water and Ice Ih from Ab-Initio Molecular Dynamics with a Non-Local Correlation Functional. Energies 2015, 8, 9383−9391. (44) Gillan, M. J.; Alfe, D.; Michaelides, A. Perspective: How Good is DFT for Water? J. Chem. Phys. 2016, 144, 130901. (45) Nosé, S. A. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519. (46) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (47) Soler, J. M.; Artacho, E.; Gale, J. D.; Garcia, A.; Junquera, J.; Ordejon, P.; Sanchez-Portal, D. The SIESTA Method for Ab Initio Order-N Materials Simulations. J. Phys.: Condens. Matter 2002, 14, 2745−2779. (48) Borysow, J.; del Rosso, L.; Celli, M.; Moraldi, M.; Ulivi, L. Spectroscopic and Thermodynamic Properties of Molecular Hydrogen Dissolved in Water at Pressures up to 200 MPa. J. Chem. Phys. 2014, 140, 164312. (49) English, N. J.; Carroll, D. G. Prediction of Henry’s Law Constants by a Quantitative Structure-Property Relationship and Neural Networks. J. Chem. Inf. Comput. Sci. 2001, 41, 1150−1161. (50) Khintchine, A. Korrelationstheorie der Stationaren Stochastischen Prozesse. Math. Ann. 1934, 109, 604−615. (51) Gaussian 09, Revision E.01: Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian, Inc.: Wallingford, CT, 2009. (52) Kazuo, N. Infrared and Raman Spectra of Inorganic and Coordination Compounds. Part A, Theory and Applications in Inorganic Chemistry; Wiley: New York, 2009. (53) Becke, A. D. Density Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648. (54) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula Into a Functional of the Electron Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785.

(55) Vosko, S. H.; Wilk, L.; Nusair, M. Accurate Spin-Dependent Electron Liquid Correlation Energies for Local Spin Density Calculations: A Critical Analysis. Can. J. Phys. 1980, 58, 1200−1211. (56) Klamt, A.; Schuurmann, G. COSMO: A New Approach to Dielectric Screening in Solvents with Explicit Expression for the Screening Energy and its Gradient. J. Chem. Soc., Perkin Trans. 2 1993, 2, 799−805.

G

DOI: 10.1021/acs.jpcc.6b11029 J. Phys. Chem. C XXXX, XXX, XXX−XXX