ARTICLES

Sci. Rep. Kanazawa Univ.
Vol. 62, pp. 87–96, 2018

Molecular dynamics study of guest hydrogen moleculesin filled ices II and Ic
Yudha A RMAN∗1 and Shinichi M IURA†2
1 Graduate School of Natural Science and Technology, Kanazawa University, Kakuma, Kanazawa

920-1192, Japan
2 Faculty of Mathematics and Physics, Kanazawa University, Kakuma, Kanazawa 920-1192, Japan(Received July 25, 2018 and accepted in revised form August 21, 2018)

Abstract We have performed molecular dynamics calculations on fully occupied and one
vacancy systems of hydrogen-filled ices II (C1 ) and Ic (C2 ) with temperature controlled at
291 K. The pressure of 1 GPa and 4 GPa for C1 and C2 , respectively, was adopted, where thefast diffusion of guest hydrogen molecules was experimentally observed. In our calculations,
hopping processes between lattice sites were observed only in one vacancy system of C1
and C2 , indicating the importance of occupation defect for the fast diffusion. The calculated
guest diffusion coefficients for one vacancy systems were found to be in good agreementwith experimental results. The guest molecules in C1 show highly anisotropic diffusion,
moving only along the optical c axis. On the other hand, the hydrogen molecules in C2 arefound to diffuse in space isotropically.
Keywords.

1

Filled Ice Ic and II, C1 , C2 , hydrogen guest.

Introduction

Extensive research of hydrogen hydrates have reported the ability of this compound to store a large
amount of hydrogen molecules. The compound can be in the form of classical clathrate hydrates
and ice-like structures. In the classical clathrate hydrates, the hydrogen molecules are encapsulated in polyhedral cavities formed by hydrogen-bonded water molecules. The encapsulation only
involves physical interaction between cage-like structure and the accommodated gas moleculeswith no chemical mixing. Theoretical and experimental studies have revealed the potential encapsulation of hydrogen molecules within the classical clathrate hydrates with the dissolvation
through single and multiple accommodations [1–7]. The hydrogen molecules encasing have beenstudied within sII structure in the range of atmospheric pressure up to 0.36 GPa at the temperatureof 150 K to 240 K [8]. The high-pressure condition to perform the synthesize can be reduced byintroducing other molecules such as THF [9] and methane [10].
Similar to the classical clathrate hydrates, the molecule’s encapsulation in ice-like structures isperformed by filling lattice sites of ice with hydrogen molecules to obtain a stable stoichiometric∗ Corresponding author

Current address: Faculty of Mathematics and Natural Sciences, Universitas Tanjungpura,Jl. Prof. Dr. H. Hadari Nawawi, Pontianak, Kalimantan Barat, Indonesia. Email: yudhaarman@gmail.com† smiura@mail.kanazawa-u.ac.jp

88

Molecular dynamics study of guest hydrogen molecules in filled ices II and Ic

solid-compound. Rather than named as the clathrate hydrate, these hydrogen inclusions are known
to be filled ice. A high-pressure condition is needed to compensate the repulsive interaction between the host and guest molecules to keep the compound stable. Not only hydrogen, other lightnoble gaseuos molecules [11–13] can also be occupied. Compared to the respective unfilled ice
structures, filled ices are more stable in an extended range of pressure and temperature [14–18].
At low-pressure condition, hydrogen molecules can be stored in hexagonal ice. A theoretical
study has revealed the ability to store a weight fraction of hydrogen molecules of 3.8 % [19].
However, the experimental study showed the maximum value only up to 0.1 % [20–22]. It has beenreported that under a high-pressure condition, ice II and Ic can accommodate hydrogen molecules
to give stable hydrate compounds [23]. These materials are also known as hydrogen hydrate C1
and C2 , respectively. The lattice site diameter of these ices is acceptable to enclose hydrogenmolecules. The total of H2 :H2 O stoichiometry for C1 is 1:6 while 1:1 is for C2 . This guest-hostratio of C2 is equivalent to 10 % wt of H2 capacity which is potential to be used as hydrogenstorage. Under 2 GPa and 295 K, the hydrate C1 has been reported to have a hexagonal unit
cell containing 36 water molecules [23]. The full arrangement consists of the hexagonal arrayof hydrogen-bonded water molecules which lies on the surface known as the basal plane. The
planes are piled along the optical c axis to form hexagonal tube-like columns in a similar way as
the ice II. On the other hand, under 3 GPa of pressure, hydrogen hydrate C2 has been reported tohave a cubic crystal structure [23]. The arrangement of water molecules gives eight faces to form
diamond-like crystal structures. Hexagonal-shaped column of hydrogen-bonded water molecules
can be viewed on (11̄0) and (110) plane due to the stacking disorder of (111) layers. The surfacesof (111) are similar to the unoccupied cubic ice structure and perfectly coherent to the basal planeof the hexagonal ice [24].
Compared to the diffusion coefficient of the hydrogen molecule in the solid hydrogen, fastdiffusion of hydrogen molecules has experimentally been observed in the dense C1 and C2 athighly compressed systems [25, 26]. Another interesting feature is provided by the mode of theH2 molecules to diffuse concerning the crystallographic planes. It has been suggested that thehydrogen molecules diffuse anisotropically in the hydrate C1 while isotropically in the C2 [25].
However, the molecular mechanism of the guest diffusion has not been studied theoretically. Thisknowledge could be valuable in delivering useful features to design more efficient hydrogenstorage systems [27].
In this paper, we report the molecular dynamics (MD) simulation results of the hydrogen hydrate C1 and C2 . We aim to provide microscopic information on inter-site diffusion mechanism ofhydrogen molecule regarding the experimental observations in C1 and C2 structures [25, 26]. We
have performed long MD calculations to investigate the hopping events of accommodated hydrogen molecules. We prepared two types of setups on each of the filled ices. Full occupation systems
were set to adopt experimental H2 :H2 O stoichiometry, in which 1:6 for C1 and 1:1 for C2 [23]. Ina similar way to the previous study [28–30], we also simulated the one vacancy system preparedby taking out a hydrogen molecule from a randomly chosen lattice site of the full occupied system;
here, one lattice site left empty is considered as an occupation defect that plays an essential role
for the hydrogen diffusion in the filled ices.
The rest of this paper is organized as follows. In Sec. 2, we describe the MD methods, systems
setup complemented with potential parameters and pair-potential energy calculation. In Sec. 3, we
present the results and discussion. Finally, concluding remarks are given in Sec. 4.

89

Yudha A RMAN and Shinichi M IURA

2

Computational Details

We carried out molecular dynamics (MD) calculations of hydrogen hydrate C1 and C2 systems
by using LAMMPS molecular dynamics package [33]. The system’s pressure was 1 GPa for C1
and 4 GPa for the C2 . Meanwhile, the temperature for all simulated system was controlled to be
291 K. These thermodynamic values were derived from experimental conditions where the fastdiffusion of H2 molecules was found [25, 26]. The NPT condition was applied by the LAMMPS
standard procedure [34–39]. Subsequently, the desired temperature and pressure were controlled
by implementing a chain of 4 thermostats coupled to the system particles, while the other four tothe barostat variables.
To model the intermolecular interaction, Lennard-Jones (LJ) and the Coulomb potential()"#
 σ 12  σ 6
qi q j e2ijij+
φ = ∑ 4εi j
−
(2.1)
ri j
ri j
ri j
i, j
were used; long-range Coulomb interaction was evaluated by using a particle-particle particlemesh method [40]. To describe the guest molecules, one site H2 model [31] was chosen with
TIP4P/Ice [32] to specify the water-water interaction. SHAKE algorithm [41] was used to handle
rigid water molecules. All of the potential parameters are listed in Table 1. The parametersof intermolecular interaction between water and hydrogen guest molecules were given by using
Lorentz-Berthelot combination rules [42, 43].
We performed MD calculations in several stages. Firstly, guest molecules were inserted into
the lattice site of ice structure to subsequently entering energy minimization step. All of thechosen ice configurations satisfy Bernal-Fowler [44] rules on 3 × 3 × 3 unit cells, containing
324 water molecules for C1 while 216 for C2 [45]. We simulated two different setups for every
filled ice system, namely full occupation and one vacancy implemented with periodic boundary
conditions. For the maximumly occupied system, the ratio between water and guest molecules
for C1 was 6:1 while 1:1 for C2 . For comparison, one vacancy was prepared by letting a single
lattice site unoccupied. Secondly, all minimum energy configurations were equilibrated for 50 ps.
During the equilibration process, energy conservation and controlled thermodynamic conditions
were carefully observed. Finally, long MD calculation of 100 ns was performed for full occupation
setup, while 50 ns for the one vacancy; longtime simulation run was needed to study the diffusion
process in C1 and C2 systems. In the equilibrium and production run, thermostat and barostat were
implemented with 0.1 ps and 1.0 ps of relaxation time, respectively. We used 1 fs time step for the
production run of all simulated systems. A van der Waals cut off radius of 10.0 Å was used in all
simulated systems and stages. All of the system’s structure presented here were drawn by using
Visual Molecular Dynamics (VMD) [46].
Table 1: The potential parameters for water and hydrogen molecules are collected. For H2 molecules, LJ
site is located at the center of mass [31]. For water molecules, TIP4P/Ice water model is adopted [32]. A
negative charge is placed at an M site, located between the center of H-H and oxygen.

Model

Site

q(|e|)σ (Å)ε(kcal/mol)

H2 [31]
TIP4P/Ice [32]CenterH
OM

0.0
0.58972.96

0.067962

3.16680.21084-1.1794

90

Molecular dynamics study of guest hydrogen molecules in filled ices II and Ic

(a)

(b)

Figure 1: The arrangement of oxygen atoms (red spheres) and H2 molecules (green spheres) in (a) hydrate
C1 and (b) C2 . Hydrate C1 is viewed parallel to the basal plane and C2 on (110) plane. Bond connecting
the adjacent oxygen atoms are drawn to clarify the structural configuration with the distance criteria 3.0 Å
for hydrate C1 and 2.8 Å for C2 .

3

Results and Discussion

We first mention the molecular dynamics results in full occupation setups. In these arrangements,
each lattice site accommodates one hydrogen molecule. The structure can be considered to beideally synthesized filled ices in which no site left unoccupied. Initially, the guest molecules havebeen placed in the center of each site in the ice structure. During 100 ns of simulation time, all theguest molecules are found to oscillate only inside its site. Eventually, we observed no inter-site
movement during the simulation time on neither hydrate C1 nor C2 .
We then present diffusion process occurred in one vacancy setups. As can be seen in the fulloccupation, the lattice site is supposed to be engaged by one hydrogen molecule. In contrast withthe full encapsulation, during 50 ns of simulation, many inter-site jumping events were observed
for the hydrogen molecules. The inter-site hopping schemes are noticeable in hydrate C1 andC2 . The movement involving only the hydrogen molecules located next to the empty site wasobserved. This clearly indicates the need for the vacancy to facilitate the inter-site translation of theencapsulated guest molecules. To quantitatively investigate the process, we calculated the meansquare displacement of a tagged hydrogen guest molecule in one vacancy systems as a function oftime t, h|∆r(t)|2 i. The mean square displacement is linked with the diffusion coefficient D byh|∆r(t)|2 i = 6Dt(3.1)

in a long time regime. In the case of the C1 structure, the hydrogen molecules are anisotropicallycoordinated by other hydrogen molecules as shown in Fig. 1 (a). Since the lattice site distancealong the c axis is about three times shorter compared to the adjacent neighboring distance on
the basal plane, the molecules are found to move only along the optical c axis. We observed thatthe inter-site diffusion of the hydrogen molecules could be seen only in a specific channel wherethe vacancy is initially set up. To take account of this observation to evaluate the diffusion coefficient, the mean square displacement was calculated for the hydrogen molecules in the specificchannel. On the other hand, as shown in Fig. 1 (b) for the hydrate C2 , the hydrogen molecules
are tetrahedrally coordinated by other hydrogen molecules. The molecules are supposed to movealong the direction of the coordinated molecules. Under this circumstances, hydrogen molecules
would travel along the tetrahedral direction on the inter-site translation. Therefore, all of hydrogenmolecules movements are considered to contribute to the diffusion process. The calculated meansquare displacements are presented in Fig. 2. For 2 ns, all plots show linearly increasing MSD.

Yudha A RMAN and Shinichi M IURA

91

Figure 2: Mean square displacement of the hydrogen molecule as a function of time t in one vacancy
configurations.

The evaluated diffusion coefficients are collected in Table 2. The diffusion coefficients for thehydrate C1 and C2 are found to be in good agreement with the experimental results [25, 26]; here,
the fast diffusion of hydrogen molecules is successfully reproduced by our simulations.
The guest jumping among sites occur sequentially. In more details, these inter-site movements
are presented in Fig. 3. Time evolution of one tagged hydrogen molecule on each filled ice system
is selected to demonstrate details of the inter-site movement. In the case of the hydrate C1 , theinter-site translation only occurred along the optical c axis; an example of this anisotropic hydrogen diffusion is depicted in Fig. 3 (a). The figure shows the time evolution of a tagged hydrogenmolecule from 110.8 ps to 111.0 ps pictured along the basal plane of ice II structure. Shown byyellow spheres, the movement of a tagged hydrogen molecule in the channel took 200 fs to complete. During the first 100 fs, the hydrogen molecule already on its halfway towards the new latticesite position. Distance covered from this early stage is 2.55 Å while the second is 1.49 Å. With
only one vacancy in the whole system, most of the hydrogen molecules oscillate within its latticesite. Inter-site movement is observed only when the adjacent lattice site is unoccupied.
In the case of the hydrate C2 system, the evolution of a tagged hydrogen molecule from 24.8748
ns to 42.3493 ns is depicted in Fig. 3 (b). All of the molecules are viewed on the hexagonal
puckered-plane of the cubic ice structure. The trajectory drawn in the figure corresponds to thepathway of one hydrogen molecule when performing inter-site jumping during the time period.
To clarify the jumping process, three different sphere colors are used to describe the trajectorytaken by the guest molecule on performing inter-site movement at a consecutive time. The yellow and blue spheres show the jumping parallel to the (011) plane observed at 24.8748 ns and42.3491 ns, respectively. The green particles are used to demonstrate the movement along the(01̄1) plane at 36.2858 ns. All of the inter-sites translation took 200 fs to complete with distinctive
covered distance throughout stages in the period; here, the guests are hopping between the adjacent hexagonal vertex formed by tetrahedrally coordinated hydrogen molecules as shown in Fig. 1(b). The inter-site diffusion is terminated since no available vacancy to prohibit the occupation. In
the meantime, the hydrogen molecule oscillates inside its new lattice site during the time interval
Table 2: Diffusion coefficient of hydrogen molecules (D) in units of cm2 s−1 . C1 and C2 are for filled ice
II and Ic, respectively.

System

D

Exp.

C1
C2

3.852(3)×10−8
5.1798(9)×10−9

1.3×10−8 [25]' 10−9 [26]92

Molecular dynamics study of guest hydrogen molecules in filled ices II and Ic

24.8748ns

111.0ps

24.875ns

24.8749ns

110.9ps36.2858ns
36.2859ns
42.3491ns

42.3492ns

110.8ps36.286ns

42.3493ns

(a)

(b)

Figure 3: Diffusion mechanism of a tagged hydrogen molecule in (a) hydrate C1 viewed along ice II basal
plane and (b) hydrate C2 drawn parallel to (110) plane of ice Ic. Red spheres are drawn to describe oxygen
atoms of water molecules, connected by hydrogen bond shown by the red stick. Blue arrows indicate thedirection of movement along the trajectory. Colored spheres (green, blue and yellow) are for hydrogenmolecules. In both figures, similar color is used to represent the same hopping scheme in successive time.
In the case of hydrate C1 , yellow spheres indicate the hydrogen movement along the optical c axis. In
hydrate C2 , the inter-site diffusion mechanism parallel to the (01̄1) plane is shown by the yellow and blue
spheres, while for the scheme parallel to the (011) plane, it is represented by the green spheres.

between the observed jumping process.
To provide more information, we present a potential field map of pair interaction between atagged H2 molecule and its surrounding water molecules. The calculation was performed by usinga system snapshot at a consecutive time before an inter-site jumping process is observed. We
adopted a similar method developed in Ref. [47]. We calculated the map as follows. Firstly, asnapshot of the filled ice system is selected. Then, the projection surface is calculated concerning
the hexagonal-puckered plane of ice structure. We then divided the surface into small rectangularspaces of 0.1 Å. Subsequently, we move the tagged molecule to the grid on the projection surface
and calculate the total potential energy. Finally, the energy values are contoured with the interval
defined by a × 4n , where a is the starting value of potential energy, and n is the number of thecontour line. The results are presented in Figs. 4 and 5. We calculate the potential maps to show
the availability of a pathway through sites concerning the potential energy. This information willlead to the chosen jumping trajectories performed by hydrogen molecules, and the track indeed
comparable to the unit cells of hydrogen hydrate. Within these figures, the local maxima areshown when the trial tagged hydrogen molecule position coincides with the occupied lattice site
location and every water molecules coordinate along the calculated crystallographic plane. On
the other hand, the local minima area can be found when the trial position of a tagged hydrogenmolecule post in an unoccupied calculation grid of lattice point and its actual position. It also canbe seen in the figure that there are energy barriers found between local minima’s points and theposition of a tagged hydrogen molecule. This feature inhibits the tagged molecule to move acrossthe confinement to accommodate other low potential energy areas. The similar pattern is observedon all calculated planes. However, we find relatively wide-flat-area that covers tagged hydrogenmolecule location. This low potential field area encompasses the surface where the vacancy islocated. The area can be found in Fig. 4 (b), Fig. 5 (a) and Fig. 5 (b) where hydrogen molecules
perform inter-site diffusion along the depicted crystallographic planes. To connect the localized
area of potential minima feature with the inter-site translation, we observed the path followed by
the tagged molecule in the next consecutive time-frame. It is found that the diffusing molecule

93

Yudha A RMAN and Shinichi M IURA

20
20

15
15

10
10

5
5

0

5

10

5.0

(a)

15

20

0

5

10

5.0

15

20

(b)

Figure 4: Potential energy contour of one vacancy C1 system, snapshots are drawn at 0.36619 ns on thesurface parallel to (a) basal and (b) optical plane. On the basal plane, ten contour lines are used with thelowest represent energy value of 1 kcal/mol and the highest is 2.6 × 105 kcal/mol. Meanwhile, there are
14 contour lines used on the optical plane; here, the lowest energy line represents 6.25 ×10−3 kcal/mol.
Projected hydrogen molecules and oxygens of water molecules are represented by gray rectangles and black
spheres, respectively. The tagged hydrogen molecule to perform inter-site translation is represented by ablack triangle, and the vacancy location is described by a red sphere. The dashed blue line shows the pathtaken by the tagged hydrogen molecule when performing inter-site jumping in the next consecutive time.

moves across this low potential field area. In constrast, a slightly different pattern occurred in thebasal surface of the hydrate C1 . As can be seen in Fig. 4 (a), energy barriers are built between
the tagged hydrogen molecule and the vacancy location, thus prohibiting the inter-site movement
along the plane; the guest hopping is only noticeable along the optical c axis.

4

Concluding Remarks

In the present study, we have performed molecular dynamics calculations for the hydrogen hydrate filled ice II and Ic; the former is called C1 , and the latter is C2 . Our calculated diffusion
coefficients are found to be in good agreement with experimental results; the fast diffusion of thehydrogen molecules has successfully been reproduced. We observed the anisotropic diffusion of
the hydrogen molecules for C1 and the isotropic for C2 . We found that the occupation defect of
guest molecules or vacancy is essential to facilitate the inter-site translation; no energy barrier on
the movement is seen when the adjacent nearest neighbor lattice site is unoccupied.
Our next work would be the hydrogen hydrates at a low temperature. The recent experimental study has revealed the ability to synthesize filled ice system at low temperature and pressure
with partial occupation [48]. Under the low-temperature condition, quantum effect of hydrogenmolecule dynamics is expected to be significant. It has been reported in hexagonal ice that quantum effect at 220 K gives noticeable contrast in lattice vibrations and potential energy in comparison to the classically approximated systems [49]. This issue will be addressed in the near future.

94

Molecular dynamics study of guest hydrogen molecules in filled ices II and Ic

20

20

15

15

10

10

5

5

0

0

-10

-5

0

(a)

5.0

5

10

5

10

15

5.0

20

25

(b)

Figure 5: Contour plot of one vacancy hydrogen hydrate C2 potential energy drawn on the (a) (110) surfaceand (b) (1̄10) plane of cubic ice. The snapshot is depicted at 0.40106 ns for the (110) surface while on the(1̄10) is at 0.60572 ns. There are eight isolines to describe potential energy map for the (110) surface, withthe lowest value is 64 kcal/mol. Similarly, eight increasing contour lines are used to describe the potential
field on the (1̄10) face. The line represents energy value started at 16 kcal/mol. In both figures, the taggedhydrogen molecule is represented by a black triangle and the vacancy position is denoted by a red sphere.
Black spheres and gray rectangles indicate the projected oxygen position of water molecules and otherhydrogen molecules on the plane, respectively. The dashed blue line indicates the inter-site diffusion pathfollowed by the tagged hydrogen molecule at the next consecutive time.

Acknowledgments
YA would like to show his gratitude to Directorate General of Higher Education, Ministry of
National Education, Indonesia for the support.

References[1] V. V. Struzhkin, B. Militzer, W. L. Mao, H.-k. Mao, and R. J. Hemley, “Hydrogen Storage in Molecular
Clathrates,” Chemical Reviews, vol. 107, pp. 4133–4151, Oct. 2007.
[2] W. L. Mao and H.-k. Mao, “Hydrogen Storage in Molecular Compounds,” Proceedings of the National Academyof Sciences, vol. 101, pp. 708–710, Jan. 2004.
[3] H. Lee, J.-w. Lee, D. Y. Kim, J. Park, Y.-T. Seo, H. Zeng, I. L. Moudrakovski, C. I. Ratcliffe, and J. A. Ripmeester,“Tuning Clathrate Hydrates for Hydrogen Storage,” Nature, vol. 434, pp. 743–746, Apr. 2005.
[4] Y. H. Hu and E. Ruckenstein, “Clathrate Hydrogen Hydrate - A Promising Material for Hydrogen Storage,”
Angewandte Chemie International Edition, vol. 45, pp. 2011–2013, Mar. 2006.
[5] N. I. Papadimitriou, I. N. Tsimpanogiannis, C. J. Peters, A. T. Papaioannou, and A. K. Stubos, “Hydrogen Storagein sH Hydrates: A Monte Carlo Study,” The Journal of Physical Chemistry B, vol. 112, pp. 14206–14211, Nov.
2008.
[6] K. Katsumasa, K. Koga, and H. Tanaka, “On the Thermodynamic Stability of Hydrogen Clathrate Hydrates,” The
Journal of Chemical Physics, vol. 127, p. 044509, July 2007.
[7] T. A. Strobel, C. A. Koh, and E. D. Sloan, “Water Cavities of sH Clathrate Hydrate Stabilized by Molecular
Hydrogen,” The Journal of Physical Chemistry B, vol. 112, pp. 1885–1887, Feb. 2008.
[8] R. V. Belosludov, R. K. Zhdanov, O. S. Subbotin, H. Mizuseki, M. Souissi, Y. Kawazoe, and V. R. Belosludov,“Theoretical modelling of the phase diagrams of clathrate hydrates for hydrogen storage applications,” Molecular
Simulation, vol. 38, pp. 773–780, Sept. 2012.
[9] L. J. Florusse, “Stable Low-Pressure Hydrogen Clusters Stored in a Binary Clathrate Hydrate,” Science, vol. 306,
pp. 469–471, Oct. 2004.

Yudha A RMAN and Shinichi M IURA

95[10] R. V. Belosludov, R. K. Zhdanov, O. S. Subbotin, H. Mizuseki, Y. Kawazoe, and V. R. Belosludov, “Theoretical
study of hydrogen storage in binary hydrogen-methane clathrate hydrates,” Journal of Renewable and Sustainable
Energy, vol. 6, p. 053132, Sept. 2014.
[11] D. Londono, W. F. Kuhs, and J. L. Finney, “Enclathration of Helium in Ice II: the First Helium Hydrate,” Nature,
vol. 332, pp. 141–142, Mar. 1988.
[12] Y. A. Dyadin, E. G. Larionov, E. Y. Aladko, A. Y. Manakov, F. V. Zhurko, T. V. Mikina, V. Y. Komarov, and E. V.
Grachev, “Clathrate Formation in Water-Noble Gas (Hydrogen) Systems at High Pressures,” Journal of Structural
Chemistry, vol. 40, pp. 790–795, Sept. 1999.
[13] L. Hakim, K. Koga, and H. Tanaka, “Novel Neon-Hydrate of Cubic Ice Structure,” Physica A: Statistical Mechanics and its Applications, vol. 389, pp. 1834–1838, May 2010.
[14] W. L. Vos, L. W. Finger, R. J. Hemley, and H.-k. Mao, “Pressure Dependence of Hydrogen Bonding in a Novel
H2 O-H2 clathrate,” Chemical Physics Letters, vol. 257, pp. 524–530, Aug. 1996.
[15] V. F. Petrenko and R. W. Whitworth, Physics of Ice. Oxford: Oxford University Press, 2006. OCLC: 455956454.
[16] A. D. Fortes, I. G. Wood, J. P. Brodholt, and L. Vočadlo, “Ab Initio Simulation of the Ice II Structure,” The Journalof Chemical Physics, vol. 119, pp. 4567–4572, Aug. 2003.
[17] L. Hakim, K. Koga, and H. Tanaka, “Phase Behavior of Different Forms of Ice Filled with Hydrogen Molecules,”Physical Review Letters, vol. 104, Mar. 2010.
[18] L. Hakim, K. Koga, and H. Tanaka, “Thermodynamic Stability of Hydrogen Hydrates of Ice Ic and II Structures,”Physical Review B, vol. 82, Oct. 2010.
[19] T. A. Pascal, C. Boxe, and W. A. Goddard, “An Inexpensive, Widely Available Material for 4 wt % Reversible
Hydrogen Storage Near Room Temperature,” The Journal of Physical Chemistry Letters, vol. 2, pp. 1417–1420,
June 2011.
[20] V. S. Efimchenko, V. E. Antonov, O. I. Barkalov, A. I. Beskrovnyy, V. K. Fedotov, and S. N. Klyamkin, “Phasetransitions and equilibrium hydrogen content of phases in the water–hydrogen system at pressures to 1.8 kbar,”
High Pressure Research, vol. 26, pp. 439–443, Dec. 2006.
[21] N. J. English, P. D. Gorman, and J. M. D. MacElroy, “Mechanisms for thermal conduction in hydrogen hydrate,”
The Journal of Chemical Physics, vol. 136, p. 044501, Jan. 2012.
[22] J.-Y. Wu, L.-J. Chen, Y.-P. Chen, and S.-T. Lin, “Molecular Dynamics Study on the Equilibrium and Kinetic
Properties of Tetrahydrofuran Clathrate Hydrates,” The Journal of Physical Chemistry C, vol. 119, pp. 1400–
1409, Jan. 2015.
[23] W. L. Vos, L. W. Finger, R. J. Hemley, and H.-k. Mao, “Novel H2 -H2 O Clathrates at High Pressures,” PhysicalReview Letters, vol. 71, pp. 3150–3153, Nov. 1993.
[24] M. A. Carignano, “Formation of Stacking Faults during Ice Growth on Hexagonal and Cubic Substrates,” The
Journal of Physical Chemistry C, vol. 111, pp. 501–504, Jan. 2007.
[25] T. Okuchi, M. Takigawa, J. Shu, H.-k. Mao, R. J. Hemley, and T. Yagi, “Fast Molecular Transport in Hydrogen
Hydrates by High-Pressure Diamond Anvil Cell NMR,” Physical Review B, vol. 75, p. 144104, Apr. 2007.
[26] T. Okuchi, “Collision and Diffusion Dynamics of Dense Molecular Hydrogen by Diamond Anvil Cell Nuclear
Magnetic Resonance,” The Journal of Physical Chemistry C, vol. 116, pp. 2179–2182, Jan. 2012.
[27] H. P. Veluswamy, R. Kumar, and P. Linga, “Hydrogen Storage in Clathrate Hydrates: Current State of the Art andFuture Directions,” Applied Energy, vol. 122, pp. 112–132, June 2014.
[28] A. Demurov, R. Radhakrishnan, and B. L. Trout, “Computations of Diffusivities in Ice and CO2 Clathrate Hydrates via Molecular Dynamics and Monte Carlo Simulations,” The Journal of Chemical Physics, vol. 116,
pp. 702–709, Jan. 2002.
[29] S. Alavi and J. Ripmeester, “Hydrogen-Gas Migration through Clathrate Hydrate Cages,” Angewandte ChemieInternational Edition, vol. 46, pp. 6102–6105, Aug. 2007.
[30] H. Lo, M.-T. Lee, and S.-T. Lin, “Water Vacancy Driven Diffusion in Clathrate Hydrates: Molecular Dynamics
Simulation Study,” The Journal of Physical Chemistry C, vol. 121, pp. 8280–8289, Apr. 2017.
[31] M. B. Sevryuk, J. P. Toennies, and D. M. Ceperley, “Why are Para-Hydrogen Clusters Superfluid? A Quantum
Theorem of Corresponding States Study,” The Journal of Chemical Physics, vol. 133, p. 064505, Aug. 2010.
[32] J. L. F. Abascal, E. Sanz, R. Garcı́a Fernández, and C. Vega, “A Potential Model for the Study of Ices andAmorphous Water: TIP4p/Ice,” The Journal of Chemical Physics, vol. 122, p. 234511, June 2005.
[33] S. Plimpton, “Fast Parallel Algorithms for Short-Range Molecular Dynamics,” Journal of Computational Physics,
vol. 117, pp. 1–19, Mar. 1995.

96

Molecular dynamics study of guest hydrogen molecules in filled ices II and Ic[34] M. E. Tuckerman, J. Alejandre, R. López-Rendón, A. L. Jochim, and G. J. Martyna, “A Liouville-Operator Derived Measure-Preserving Integrator for Molecular Dynamics Simulations in the Isothermal–Isobaric Ensemble,”
Journal of Physics A: Mathematical and General, vol. 39, no. 19, p. 5629, 2006.
[35] S. Nosé, “A Unified Formulation of the Constant Temperature Molecular Dynamics Methods,” The Journal ofChemical Physics, vol. 81, pp. 511–519, July 1984.
[36] W. G. Hoover, “Canonical Dynamics: Equilibrium Phase-Space Distributions,” Physical review A, vol. 31, no. 3,
p. 1695, 1985.
[37] W. Shinoda, M. Shiga, and M. Mikami, “Rapid Estimation of Elastic Constants by Molecular Dynamics Simulation Under Constant Stress,” Physical Review B, vol. 69, p. 134103, Apr. 2004.
[38] G. J. Martyna, D. J. Tobias, and M. L. Klein, “Constant Pressure Molecular Dynamics Algorithms,” The Journalof Chemical Physics, vol. 101, pp. 4177–4189, Sept. 1994.
[39] M. Parrinello and A. Rahman, “Polymorphic Transitions in Single Crystals: A New Molecular Dynamics
Method,” Journal of Applied Physics, vol. 52, pp. 7182–7190, Dec. 1981.
[40] R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles. Bristol [England] Adam Hilger,
special student ed ed., 1988.
[41] J.-P. Ryckaert, G. Ciccotti, and H. J. C. Berendsen, “Numerical Integration of the Cartesian Equations of Motionof a System with Constraints: Molecular Dynamics of n-Alkanes,” Journal of Computational Physics, vol. 23,
pp. 327–341, Mar. 1977.
[42] H. A. Lorentz, “Ueber die Anwendung des Satzes vom Virial in der kinetischen Theorie der Gase,” Annalen der
Physik, vol. 248, pp. 127–136, Jan. 1881.
[43] D. Berthelot, “Sur la Détermination Rigoureuse des Poids Moléculaires des Gaz en Partant de leurs Densités et del’écart que celles-ci Présentent par Rapport á la Loi de Mariotte, Compt,” Rendus, vol. 126, pp. 954–956, 1898.
[44] J. D. Bernal and R. H. Fowler, “A Theory of Water and Ionic Solution, with Particular Reference to Hydrogen andHydroxyl Ions,” The Journal of Chemical Physics, vol. 1, pp. 515–548, Aug. 1933.
[45] V. Buch, R. Martoňák, and M. Parrinello, “A New Molecular-Dynamics Based Approach for Molecular Crystal
Structure Search,” The Journal of Chemical Physics, vol. 123, p. 051108, Aug. 2005.
[46] W. Humphrey, A. Dalke, and K. Schulten, “VMD: Visual Molecular Dynamics,” Journal of Molecular Graphics,
vol. 14, pp. 33–38, Feb. 1996.
[47] T. Ikeda-Fukazawa, S. Horikawa, T. Hondoh, and K. Kawamura, “Molecular Dynamics Studies of Molecular
Diffusion in Ice Ih,” The Journal of Chemical Physics, vol. 117, pp. 3886–3896, Aug. 2002.
[48] R. Kumar, D. D. Klug, C. I. Ratcliffe, C. A. Tulk, and J. A. Ripmeester, “Low-Pressure Synthesis and Characterization of Hydrogen-Filled Ice Ic,” Angewandte Chemie International Edition, vol. 52, pp. 1531–1534, Jan.
2013.
[49] L. Hernández de la Peña, M. S. Gulam Razul, and P. G. Kusalik, “Quantum Effects in Ice Ih,” The Journal ofChemical Physics, vol. 123, p. 144506, Oct. 2005.