Molecular dynamics study on fast diffusion of hydrogen molecules in filled ice II

Journal of Molecular Liquids 292 (2019) 111316

Contents lists available at ScienceDirect

Journal of Molecular Liquids journal homepage: www.elsevier.com/locate/molliq

Molecular dynamics study on fast diffusion of hydrogen molecules in filled ice II Asuka Harada a,⁎, Yudha Arman b, Shinichi Miura c a b c

Division of Mathematics and Physics, Graduate School of Natural Science and Technology, Kanazawa University, Kakuma, Kanazawa 920-1192, Japan Faculty of Mathematics and Natural Sciences, Universitas Tanjungpura, Pontianak, Kalimantan Barat 78115, Indonesia Faculty of Mathematics and Physics, Kanazawa University, Kakuma, Kanazawa 920-1192, Japan

a r t i c l e

i n f o

Article history: Received 26 April 2019 Received in revised form 1 July 2019 Accepted 4 July 2019 Available online 05 July 2019 Keywords: Hydrogen hydrates Filled ice Ice II Molecular dynamics method

a b s t r a c t Molecular dynamics calculations have been performed for the hydrogen hydrate C1 that is the ice II filled by the hydrogen molecules. Fast diffusion of the hydrogen molecules in the C1 crystal was observed; the calculated selfdiffusion coefficient is in good agreement with the experimental value. The hydrogen molecules were observed to migrate only inside the tube-like voids in the C1 and the resulting diffusion processes were highly anisotropic in space. To diffuse the hydrogen molecules, defect or vacancy in the tube-like voids was found to be essential; no diffusion was observed for the system without vacancies. © 2019 Elsevier B.V. All rights reserved.

1. Introduction At low temperature and high pressure, water molecules together with other gas molecules are known to be transformed into solid states called hydrates [1]. Molecules incorporated into the framework of hydrogen bonded water molecules are called guest molecules; typical examples are given by methane, hydrogen, and other gas molecules [2–4]. The hydrates can be classified into two types, clathrate hydrate and filled ice. In the clathrate hydrate [1], water molecules form cage-like structures by hydrogen bonding; guest molecules are confined in the cages. There is no chemical bond between the guest and water molecules; the structure of the clathrate hydrate is maintained by only weak van der Waals interaction. The clathrate hydrates can store gas molecules at high density; for example, the methane hydrate can encapsulate gas molecules up to 160 times the volume of the clathrate itself [2]. On the other hand, in the case of the clathrate hydrate of the hydrogen, very severe conditions are required for synthesis; however, the conditions can be relaxed by incorporating other molecules such as THF [5] and methane [6]. By applying pressure to the clathrate hydrate, the cage structure is transformed into the filled ice structure where guests occupy the space of the ice framework [7] if the guest molecule is very small such as hydrogen or a light noble gas [8–10]. This structure allows even higher density storage of guest molecules. Among various hydrates, the hydrogen hydrate has been attracting attention as next generation energy source, since the hydrogen ⁎ Corresponding author. E-mail address: [email protected] (A. Harada).

https://doi.org/10.1016/j.molliq.2019.111316 0167-7322/© 2019 Elsevier B.V. All rights reserved.

molecules produce only water as the byproduct of fuel reaction. At low pressure condition, the hydrogen molecules can be filled in the ice Ih. A previous study has revealed the ability to store a weight fraction of hydrogen molecules of 3.8% [11]. However, the experimental study showed that 0.1% is the limit of filling the ice Ih with the hydrogen [12]. It is known that other ice structures, ice II and ice Ic can stably store hydrogen molecules at higher pressures [13–15]. Substances packed with hydrogen molecules in the structures of the ice II and Ic are called C1 and C2, respectively; it means “Compound 1” and “Compound 2”. When the hydrogen molecules are filled in all of the void sites, those can be stored at high density such that the molar ratio of water to hydrogen is 6:1 in C1 and is 1:1 in C2 [16]. Recent experiments have revealed that the hydrogen molecules can diffuse rapidly in the filled ices [17,18]; this phenomenon could be important to efficient hydrogen-storage [19]. However, the detailed diffusion mechanism in the dense crystals has not been well understood. In this paper, we report molecular dynamics simulation results to elucidate the diffusion mechanism of the hydrogen molecules in the hydrogen hydrate C1. 2. Computational details In the present study, the hydrogen molecule was modeled as a rigid diatomic molecule consisting of two hydrogen atoms; the bond length was fixed to be 0.7414 Å [20]. To describe the intermolecular interaction, three interaction sites were introduced, two hydrogen atoms and center-of-mass of the molecule. Two positive charges were located at the atomic sites and a negative charge was put at the center-of-mass of the molecule together with the Lennard-Jones potential. The TIP4P/

2

A. Harada et al. / Journal of Molecular Liquids 292 (2019) 111316

Table 1 Parameters of models used in this study. TIP4P/Ice [21] is a rigid 4-site model of water molecule having dummy site on a bisecting line of bond angle, which is a massless dummy site carrying negative charge. Hydrogen molecule is modeled as a 3-site rigid model consisting of two hydrogen atoms and a massless center-of-mass site carrying negative charge [20]. Model

Site

q (e)

TIP4P/Ice

H O M H Center

+0.5897

H2

−1.1794 +0.4932 −0.9864

σ (Å)

ε (kcal/mol)

3.1668

0.21084

3.038

0.06816

1 V systems using LAMMPS Molecular Dynamics Simulator [25]. Molecular dynamics calculations have been carried out 50 ns with time step 1 fs. Rigid body condition was handled by the SHAKE method [26]. The pressure and temperature of each system were controlled to be 1 GPa and 291 K, corresponding to the thermodynamic condition of the experiment where the fast diffusion of the hydrogen molecules was found [17]. NPT condition was maintained by the standard LAMMPS procedure [27–32].

3. Results

Fig. 1. Initial arrangement of water (red and white) and hydrogen molecules (blue) in the hydrogen hydrate C1 drawn by VMD molecular graphic software [35]. The c-axis is directed perpendicular to the paper surface. There exist hexagonal voids along the c-axis.

Ice model [21] was adopted as the water molecule. All the potential parameters are listed in Table 1. The Lorentz-Berthelot combination rule [22] was applied to the Lennard-Jones parameters between water and hydrogen molecules. The cut-off distance of van der Waals interaction was 10.0 Å and long-range Coulomb potential was calculated by particle-particle particle-mesh method [23]. In the initial configuration, water molecules were arranged to satisfy the Bernal-Fowler rule [24] on 3 × 3 × 3 unit cell, containing 324 water molecules. The molar ratio of water molecules to hydrogen molecules was 6:1; there were 54 hydrogen molecules in the simulation box. For comparison, we examined another system where one hydrogen molecule was removed from the filled ice; the system contained one vacancy site of the hydrogen molecule. The former system is called fully occupied (FO) system, and the latter 1 vacancy (1 V) containing system. We have performed molecular dynamics calculations for the FO and

We first show the snapshot of the hydrogen hydrate C1 in Fig. 1. This crystal has the rhombohedral structure. The hexagonal tube-like voids containing hydrogen molecules are clearly seen; the tube-like voids are extended along the c-axis. In Fig. 2, the orientational distribution of the hydrogen molecules in the C1 for the fully occupied (FO) system is presented. The orientation is described by the bond vector of each hydrogen molecule, and is represented by the azimuthal and polar angles of the bond vector. The polar angle is measured from the z-axis and the azimuthal angle is introduced by the projected bond vector on the x-y plane. The figures show that there is almost no bias on either of the azimuthal and polar angle distribution, which indicates the hydrogen molecules rotate freely in the C1 crystal. We also analyzed the orientational distribution for the 1 vacancy (1 V) containing system; no difference was found for the FO and 1 V systems. In the present study, therefore, we focus on the center-of-mass behavior of hydrogen molecules in the C1 crystal. Here, we summarize observations found in the trajectory both for the FO and 1 V systems. We found no diffusion of hydrogen molecules for the FO system during the 50 ns simulation; all the hydrogen molecules were fluctuated around the sites where the molecules were initially located. On the other hand, active diffusion processes were observed for the 1 V system; the diffusing molecules were found to move in a highly anisotropic manner. The molecules move only inside the tube parallel to the c-axis. The diffusion was observed only for the tube where the vacancy was generated; no diffusion was occurred inside the fully occupied tubes. The snapshot of the tube containing hydrogen molecules is presented in Fig. 3. The snapshots are presented for the initial and final configurations in a 5 ns duration. A vacancy site is found at the top of the tube in the initial configuration. The vacancy was found to migrate downward in the tube and reach at the bottom of the tube as seen in the final configuration. To see more detailed behavior, the diffusion or hopping processes of the hydrogen molecules are diagramed in Fig. 4. The hydrogen molecule next to the vacancy is found to hop to the vacancy. Some molecules hop one time and then, the vacancy is occupied another molecule. Other molecules repeat hopping forward and backward during several nanoseconds. Two or more hydrogen molecules were not observed to occupy the single site simultaneously.

Fig. 2. Orientational distribution of hydrogen molecules in the hydrogen hydrate C1: (a) azimuthal ϕ and (b) polar angle θ distributions.

A. Harada et al. / Journal of Molecular Liquids 292 (2019) 111316

3

7

6

MSD (Å2)

5

4

3

2

1

0

0.5

1.0

1.5

t (ns)

2.0

2.5

3.0

Fig. 5. Mean square displacement (MSD) of a tagged hydrogen molecule in fully occupied (FO) and 1 vacancy (1 V) containing hydrogen hydrate C1. The red curve is for the FO system and the blue one for the 1 V system.

Fig. 3. Snapshots of five hydrogen molecules in the tube-like vacancy along the c-axis in the hydrogen hydrate C1 drawn by VMD molecular graphic software [35]. Each sphere represents the center-of-mass position of each hydrogen molecules with different colors. Left one is a configuration at an instant and right one is the configuration after 5 ns calculation from the left configuration.

In order to quantitatively discuss the diffusion processes, we have evaluated the self-diffusion coefficient D for the hydrogen molecules. We show the mean square displacement (MSD) of a tagged molecule 〈|Δr(t)|2〉 where t is the time interval of the measurement. While the average was only taken over the hydrogen molecules in the tube where the vacancy exists for the 1 V system, all the hydrogen molecules is taken into account for the FO system. The calculated MSDs are presented both for the FO and 1 V systems in Fig. 5. For the 1 V system, the MSD increases linearly in time, which is the characteristic behavior of the diffusion processes. On the other hand, for the FO system, the MSD does not grow in time and converges to a value less than 1 Å2for

6

Lattice site

5

4

3

the long time regime, reflecting all the hydrogen molecules oscillate around the sites where the molecules were initially located. This MSD function is related with the self-diffusion coefficient D by the following relation: D

E jΔrðt Þj2 ¼ 2nDt

ð1Þ

for the long time regime; here, n is the dimensionality of the diffusion processes. The calculated self-diffusion coefficient of the molecule is listed in Table 2 together with the experimental value. As mentioned above, however, the diffusion is highly anisotropic in space and could be better to be regarded to occur in a one-dimensional tube; the selfdiffusion coefficient with n = 1 is also presented. Both calculated values are found to be in the same order of the experimental value, but somewhat larger. It would be difficult to identify the source of the discrepancy, since the system size dependence is expected to be large for the vacancy assisted diffusion. Larger system calculations would be needed for quantitative comparison between calculated and experimental values. To confirm this anisotropic diffusion along the c-axis direction inside the tube, the displacement Δr(t) is decomposed into the component parallel and perpendicular to the c-axis, respectively denoted by (Δr(t))∥and (Δr(t))⊥. The MSDs for both components are presented in Fig. 6. While the MSD regarding the parallel component is found to increase linearly in time, the perpendicular MSD converges to a value as short as the square of the tube width. The figure clearly shows the diffusion purely occurs along the c-axis. The MSD regarding the parallel component is essentially the same as that for the 1 V system presented in Fig. 5, since the perpendicular component is not contributed to the full MSD at all. Indeed, the full MSD can be written by D

E D D 2 E 2 E D 2 E jΔrðt Þj2 ¼ ðΔrðt ÞÞ∥  þ 2 ðΔrðt ÞÞ⊥  ≅ ðΔrðt ÞÞ∥ 

ð2Þ

2

1 0.0

1.0

2.0

t(ns)

3.0

4.0

5.0

Fig. 4. Time evolution of the center-of-mass positions of the five hydrogen molecules presented in Fig. 3 where the left configuration is given at t = 0 and the right one is given at t = 5 ns. The same colour scheme is applied. The vertical axis gives the number of the lattice sites in Fig. 3 labeled from the bottom.

Table 2 Self-diffusion coefficient of hydrogen molecules in the hydrogen hydrate C1 in units of cm2/s by our molecular dynamics calculations; n is the dimensionality of the diffusion processes in Eq. (1). Experimental value is taken from the literature [17]. Present work

Experiment

n=1

n=3

8.6 × 10−8

2.9 × 10−8

1.2 × 10−8

4

A. Harada et al. / Journal of Molecular Liquids 292 (2019) 111316 7

6

MSD (Å2)

5

4

3

2

1

0

0.5

1.0

1.5

t (ns)

2.0

2.5

3.0

Fig. 6. Mean square displacement (MSD) of a tagged hydrogen molecule in 1 vacancy containing hydrogen hydrate C1. The red curve is for the MSD parallel to the c-axis and the blue and green curves are for the MSDs perpendicular to the c-axis.

The MSD and associated self-diffusion coefficient are completely determined by the c-axis parallel component of the displacement. To shed light on roles of water molecules in the hydrogenic diffusion, we have also performed molecular dynamics calculations for a system where one water molecule is removed. At the early stage of the simulation, the water vacancy is occupied by a neighboring hydrogen molecule; then, hydrogen molecules start to move in the tube to fill the generated vacancy as found in the 1 V system in a couple of nanoseconds, indicating the vacancy of the water molecules gives a minor effect in this time scale. 4. Concluding remarks In the present study, molecular dynamics calculations have been performed on the hydrogen hydrate C1 that is the ice II filled by hydrogen molecules. Attention was paid on the fast diffusion mechanics of the hydrogen molecules in the crystal. We found the importance of the vacancy site of the hydrogen molecules; the vacancy assists the hydrogenic diffusion, while no diffusion was found for fully occupied no vacancy containing crystal. The self-diffusion coefficient calculated for the vacancy containing crystal is in good agreement with the experimental value. We also found that the diffusion is highly anisotropic in space and occurs along the tube-like voids directed to the c-axis. It is known that the fast diffusion of the hydrogen molecules is observed in the hydrogen hydrate C2 [18]. Due to the lightness of the hydrogen, nuclear quantum effect could play important role in describing the properties of the hydrogen hydrates at low temperature [33,34]. These issues will be addressed in the near future. References [1] W.L. Mao, C.A. Koh, E.D. Sloan, Clathrate hydrates under pressure, Phys. Today 60 (10) (2007) 42–47. [2] A. Lenz, L. Ojamäe, Structures of the I-, II- and H-methane clathrates and the icemethane clathrate phase transition from quantum-chemical modeling with forcefield thermal corrections, J. Phys. Chem. A 115 (2011) 6169–6176. [3] W.L. Mao, H. Mao, Hydrogen storage in molecular compounds, Proc. Natl. Acad. Sci. 101 (3) (2004) 708–710.

[4] T. Nakajima, T. Kudo, R. Ohmura, S. Takeya, Y.H. Mori, Molecular storage of ozone in a clathrate hydrate: an attempt at preserving ozone at high concentrations, PLoS One 7 (11) (2012). [5] L.J. Florusse, et al., Stable low-pressure hydrogen clusters stored in a binary clathrate hydrate, Science 306 (2004) 469–471. [6] R.V. Belosludov, R.K. Zhdanov, O.S. Subbotin, H. Mizuseki, Y. Kawazoe, V.R. Belosludov, Theoretical study of hydrogen storage in binary hydrogen-methane clathrate hydrates, J. Renew. Sustain. Energy 6 (2014), 053132. [7] W.L. Vos, L.W. Finger, R.J. Hemley, H. Mao, Pressure dependence of hydrogen bonding in a novel H2O-H2 clathrate, Chem. Phys. Lett. 257 (1996) 524–530. [8] D. Londono, W.F. Kuhs, J.L. Finney, Enclathration of helium in ice II: the first helium hydrate, Nature 332 (10) (1988) 141–142. [9] Y.A. Dyadin, et al., Clathrate formation in water-Noble gas (hydrogen) systems at high pressures, J. Struct. Chem. 40 (5) (1999) 790–795. [10] L. Hakim, K. Koga, H. Tanaka, Novel neon-hydrate of cubic ice structure, Physica A 389 (2010) 1834–1838. [11] T.A. Pascal, C. Boxe, W.A. Goddard, An inexpensive, widely available material for 4 wt % reversible hydrogen storage near room temperature, J. Phys. Chem. Lett. 2 (2011) 1417–1420. [12] V.S. Efimchenko, V.E. Antonov, O.I. Barkalov, A.I. Beskrovnyy, V.K. Fedotov, S.N. Klyamkin, Phase transitions and equilibrium hydrogen content of phases in the water–hydrogen system at pressures to 1.8 kBar, High Press. Res 26 (4) (2006) 439–443. [13] S. Machida, H. Hirai, T. Kawamura, Y. Yamamoto, T. Yagi, A new high-pressure structure of methane hydrate surviving to 86 GPa and its implications for the interiors of giant icy planets, Phys. Earth Planet. Inter 155 (2006) 170–176. [14] S. Machida, H. Hirai, T. Kawamura, Y. Yamamoto, T. Yagi, Raman spectra for hydrogen hydrate under high pressure: intermolecular interactions in filled ice Ic structure, J. Phys. Chem. Solids 71 (2010) 1324–1328. [15] N.J. English, P.D. Gorman, J.M.D. MacElroy, Mechanisms for thermal conduction in hydrogen hydrate, J. Chem. Phys. 136 (2012), 044501. [16] W.L. Vos, L.W. Finger, R.J. Hemley, H. Mao, Novel H2-H2O clathrates at high pressures, Phys. Rev. Lett. 71 (19) (1993) 3150–3153. [17] T. Okuchi, M. Takigawa, J. Shu, H. Mao, R.J. Hemley, T. Yagi, Fast molecular transport in hydrogen hydrates by high-pressure diamond anvil cell NMR, Phys. Rev. B 75 (144104) (2007). [18] T. Okuchi, Collision and diffusion dynamics of dense molecular hydrogen by diamond anvil cell nuclear magnetic resonance, J. Phys. Chem. C 116 (2012) 2179–2182. [19] H.P. Veluswamy, R. Kumar, P. Linga, Hydrogen storage in clathrate hydrates: current state of the art and future directions, Appl. Energy 122 (2014) 112–132. [20] L. Hakim, K. Koga, H. Tanaka, Thermodynamic stability of hydrogen hydrates of ice Ic and II structures, Phys. Rev. B 82 (144105) (2010). [21] J.L.F. Abascal, E. Sanz, R.G. Fernández, C. Vega, A potential model for the study of ices and amorphous water: TIP4P/ice, J. Chem. Phys. 122 (234511) (2005). [22] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press, 1989. [23] R.W. Hockney, J.W. Eastwood, Computer Simulation Using Particles, special student ed. Adam Hilger, Bristol [England], 1988. [24] J.D. Bernal, R.H. Fowler, A theory of water and ionic solution, with particular reference to hydrogen and hydroxyl ions, J. Chem. Phys. 1 (1933) 515–548. [25] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117 (1995) 1–19. [26] J.P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, Numerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of nalkanes, J. Comput. Phys. 23 (1977) 327–341. [27] S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys. 81 (1984) 511–519. [28] W.G. Hoover, Canonical dynamics: equilibrium phase-space distributions, Phys. Rev. A 31 (3) (1985) 1695–1697. [29] W. Shinoda, M. Shiga, M. Mikami, Rapid estimation of elastic constants by molecular dynamics simulation under constant stress, Phys. Rev. B 69 (134103) (2004). [30] G.J. Martyna, D.J. Tobias, M.L. Klein, Constant pressure molecular dynamics algorithms, J. Chem. Phys. 101 (5) (1994) 4177–4189. [31] M. Parrinello, A. Rahman, Polymorphic transitions in single crystals: a new molecular dynamics method, J. Appl. Phys. 52 (1981) 7182–7190. [32] M.E. Tuckerman, J. Alejandre, R. López-Rendón, A.L. Jochim, G.J. Martyna, A Liouvilleoperator derived measure-preserving integrator for molecular dynamics simulations in the isothermal-isobaric ensemble, J. Phys. A Math. Gen. 39 (2006) 5629–5651. [33] M.M. Conde, C. Vega, C. McBride, E.G. Noya, R. Ramírez, L.M. Seś, Can gas hydrate structures be described using classical simulations? J. Chem. Phys. 132 (114503) (2010). [34] C.J. Burnham, Z. Futera, N.J. English, Quantum and classical inter-cage hopping of hydrogen molecules in clathrate hydrate: temperature and cage-occupation effects, Phys. Chem. Chem. Phys. 19 (2016) 717–728. [35] W. Humphrey, A. Dalke, K. Schulten, VMD: Visual Molecular Dynamics, J. Mol. Graph. 14 (1996) 33–38.