Proton Disorder and the Dielectric Constant of Type II Clathrate Hydrates

Computational studies are presented examining the degree of proton disorder in argon and molecular hydrogen sII clathrate hydrates. Results are presented using a variety of model potentials for the dielectric constant, the proton order parameter, and the molecular volume for the clathrate systems. The dielectric constant for the clathrate systems is found to be lower than the dielectric constant of ice in all models. The ratio of the clathrate to ice dielectric constant correlates well with the ratio of the densities, which is not the case for comparisons to the liquid, so that differences in the dielectric constants between ice and the clathrates are most likely due to differences in densities. Although the computed dielectric constant is a strong function of the model potential used, the ratio of the dielectric constant of ice to that of the clathrates is insensitive to the model potential. For the nonpolar guest molecules used in the current study, the degree proton of disorder is found to depend weakly on the identity of the guest but the dielectric constant does not appear to be sensitive to pressure or the type of guest.


I. INTRODUCTION
Water is not only one of the most important compounds in nature, but among the most intriguing as well. The solid phase diagram of water is particularly interesting with at least 12 distinct crystal structures. In addition to pure water, the solid phase of water has a particularly rich phenomenology and structure in the presence of gas phase solutes. Under suitable thermodynamic conditions, dissolved gases in water can transform into solid inclusion compounds where the gas solute molecules occupy sites in aqueous cage structures formed by the hydrogen bonded water molecules. Such inclusion compounds are known as clathrate hydrates, 1 and the clathrate cages principally take on one of three lattices often called sI, sII, and sH with more complex structures observed 2,3 especially for mixed systems. 4 Methane, [5][6][7][8][9][10][11][12] nitrogen, [13][14][15][16][17] argon, 18 hydrogen, [19][20][21][22][23][24][25][26] and hydrocarbons 27 are among the important gas phase species that can form stable clathrate hydrates.
Of particular interest to the current work is the sII clathrate hydrate structure that consists of cages having two distinct structures. The smaller of the two cages has dodecahedral symmetry and is often denoted 5 12 where the "5" represents the pentagonal faces of the dodecahedron and the "12" indicates that there are 12 such faces for a total of 20 water molecules. The larger of the two cages has hexakaidecahedral symmetry and is denoted 5 12 6 4 where, again, the "5 12 " denotes 12 pentagonal faces and the "6 4 " represents four hexagonal faces for a total of 28 water molecules. In the sII structure there are 136 water molecules in a unit cell of fcc symmetry. The number of solute molecules per unit cell varies depending on the details of the lattice. In the case of argon sII clathrate hydrates, every small and large cage contains one argon atom whereas in the case of molecular hy-drogen solutes, it is believed that the large hexakaidecahedral cages contain four hydrogen molecules and the small dodecahedral cages contain either one or two hydrogen molecules. 19,21 It has been known for many decades that the most stable structure of ice under ambient conditions has a residual entropy at low temperatures of approximately S = Nk B ln 3 / 2. The residual entropy in the Ih ice structure was first rationalized in a mean-field sense by Pauling,28 who showed that the protons in Ih ice are disordered because there are ͑3 / 2͒ N ways for N water molecules to satisfy the Bernal-Fowler ice rules. Briefly, the Bernal-Fowler ice rules state that a water molecule remains neutral as H 2 O and forms four hydrogen bonds, two as a donor, and two as an acceptor. 28,29 This disorder adds stability to the crystal by increasing the entropy and the disorder gives ice a high dielectric constant. Other crystal structures of ice have similar residual entropies and high dielectric constants. While to our knowledge there is no experimental information about the possible residual entropy in the clathrate hydrates, the dielectric constant of the clathrate hydrates is high and proton disorder in the clathrate systems is likely. A purpose of the current work is to study computationally the proton disorder in some important clathrate hydrates having the sII structure.
The residual entropy and associated proton disorder is complicated by the fact that not all hydrogen bond arrangements are energetically equivalent. Rather, the differing hydrogen bond arrangements depend on the orientation of the hydrogen bonded dimer, 30 as well as the interactions from other neighbors. [31][32][33] For the kind of dimers that make up a type sII clathrate, there are two possible dimer orientations, which can be labeled inverse mirror ͑IM͒ or oblique mirror ͑OM͒ ͑see Fig. 1͒. 30 The repulsive interactions between the hydrogens make the IM lower in energy; 30 however the next nearest neighbor interactions tend to favor the OM, 31  tice is a complicated balance of interactions involving distant neighbors resulting in only a small difference in energy. 32,33 Several mathematical approaches have been developed to enumerate the number of configurations for various clathrate cages. [34][35][36][37] For example, there are 30 0226 for the 20 molecules making up the 5 12 ͑Ref. 34͒ and 61 753 344 for the 28 molecules in 5 12 6 4 cage. 35 These studies also indicate the IM mirror structure should be lower in energy, when only considering interactions from water molecules in the same cage. 36,37 The mechanism of proton rearrangement involves extremely rare ͑less than 1 per 1 ϫ 10 6 water molecules at 273 K͒ defects ͑Bjerrum D and L defects͒. 38 Simulations which do not contain enough water molecules to have defects at the normal concentrations do not undergo any changes in the underlying proton structure. A simulation would remain in a single proton-disordered structure and would not sample over all the significant structures. Experimentally, the concentration of defects effectively goes to zero at a temperature below 50 K, so the dielectric constant drops to near 1 as the system becomes trapped in a single proton configuration. 2 Computer simulations can sample, in principle, all relevant proton configurations using a Monte Carlo method in which many particle moves are attempted between different structures that all obey the ice rules. 33,[39][40][41][42] This method has been applied to calculate the dielectric constant and proton order parameters for ice Ih. 33,42 In this study, this method will be used to examine type II clathrates with different guests. The goal of this work is to find answers to the following questions: • Do the clathrates show the same low-temperature disorder that is observed in Ih ice?
• Can we account for the differences in the measured dielectric constants between Ih ice and the clathrate hydrates?
• Given the known sensitivity of the calculated dielectric constant of Ih ice to the model potential used, 33 what is the sensitivity of the ratio of the ice to clathrate dielectric constants to model potential?
• What is the sensitivity of the computed dielectric constant of the clathrates to the specific identity of nonpolar gas solutes?
The contents of the remainder of this paper are as follows. In Sec. II we review, briefly, the model potentials used in this work along with the computational methods. In Sec. III we present our numerical results, and in Sec. IV we summarize our conclusions.

II. METHOD
The simulations of the argon clathrates use the Anderson, et al. potential, 9 with an exponential-6 potential between argon and the water oxygen atoms and an r −12 repulsive potential between argon and water hydrogen atoms. For the argon-argon interactions we use a Lennard-Jones potential with = 0.2375 kcal/ mol and = 3.405 Å. 43 Three different water potentials are used: SPC/E, 44 TIP4P/Ice, 45 and TIP4P-FQ/Ice. 46 The SPC/E potential is a commonly used potential, which has been applied in many studies of clathrate hydrates. The TIP4P/Ice model is a reparametrization of the TIP4P model, 47 which accurately reproduces the densities and phase coexistence properties of many of the ice phases. 45 The TIP4P-FQ/Ice model is a polarizable model, a reparametrization of the TIP4P-FQ model, 48 which gives an accurate density and melting point for ice Ih. 46 The simulations use a single unit cell with 136 water molecules, containing eight large ͑5 12 6 4 ͒ cavities and 16 small ͑5 12 ͒ cavities. For the argon clathrates, 24 argon atoms are used, one for each cavity. For the H 2 clathrates, we include 48 molecules ͑one hydrogen molecule in the 5 12 and four hydrogen molecules in the 5 12 6 4 cages͒. Three different temperatures ͑150, 200, and 250 K͒ and two different pressures ͑1 and 2 kbar͒ are used for the argon clathrates. The H 2 clathrates use the same range of temperature but only one pressure ͑2 kbar͒.
All simulations are done in the isothermal-isobaric ͑constant T-P-N͒ ensemble by coupling to a pressure bath and a Nosé-Hoover temperature bath. [50][51][52][53][54] An orthorhombic simulation box is used, with each side of the box treated as an independent variable for the constant pressure dynamics. 55 A time step of 1 fs, SHAKE to enforce bond constraints, and Ewald sums for long-ranged interactions are used. 56 Each simulation is simulated for 1 ns or more.
Proton disorder is simulated using the method described previously. 33,42 This method runs conventional molecular dynamics for a number of steps ͑here, 10 or 0.01 ps͒ then attempts a many-particle Monte Carlo move which generates new proton configurations. The generation of new proton positions involves two steps, a "walk" step, in which a random walk is made on the lattice until a closed hydrogen bonded loop is found, and a "roll" step, in which an alternate hydrogen bonded loop is made by rotating each molecule in the loop. This method is both ergodic, so that it can, in principle, generate all proton disorder configurations, and satisfy detailed balance. The dielectric constant ⑀ is calculated from the fluctuations in the total dipole moment of the system oblique mirror inverse mirror where M is the total dipole moment of the system and ⑀ 0 is the optical dielectric constant. Details associated with Eq. ͑1͒ including the choice of ⑀ 0 are given in Ref. 42. From Eq. ͑1͒ it is evident that for a given model potential there are two important sources of systematic discrepancies with experiment. The fluctuations in the total dipole moment depend both on the magnitude of the dipole moments of individual water molecules as well as the fluctuations arising from the interactions between the water molecules in the system. The variations in computed dielectric constants with different model potential can be significant, 42 owing to both contributions. In this work, our concerns relate to the nature of the proton disorder in the clathrates and comparisons with the observed disorder in ice Ih. Given that goal, we can define an order parameter that can provide a useful picture of the hydrogen bond configuration in the lattice. Recognizing that each water molecule in the sII clathrate lattice has four near neighbors, we define the proton order parameter X OM as the number of hydrogen bond near neighbors that are in the OM orientation ͑see Fig. 1͒. Because there are four nearest neighbors, the number that is in the IM orientation X IM is 4−X OM . There are twice as many ways to be in an OM orientation than the number of ways to be in an IM orientation. Consequently, a completely random lattice has 2/3 of the four hydrogen bonds as OM and 1/3 as IM, so that X OM =8/ 3.

III. RESULTS
The dielectric constant and the proton order parameter are given in Table I. The dielectric constant is fairly independent both of pressure and the identity of the guest, but depends strongly on the water model. The SPC/E and TIP4P/ Ice models are both nonpolarizable and both underestimate the dielectric constant ͑see Fig. 2͒. The polarizable model, TIP4P-FQ/Ice, appears to overestimate the dielectric constant. The results in Table I show that the hydrogen positions are not completely random, with a preference for the lower energy configuration IM ͑because X OM Ͻ 8 / 3͒. The different water models give different order parameters, but the order parameters do not appear to be pressure dependent, for the argon clathrates. The order parameters do seem to be somewhat dependent on the type of guest.
As the temperature decreases, the preference for IM increases. The ratio of X IM to X OM can be understood from the expression where ͗E ␣ ͘ is the average energy of hydrogen bonds of type ␣ at a particular temperature and density and the factor of 1/2 is required because there are twice as many OM as IM hydrogen bonds. If ͗E ␣ ͘ is independent of temperature then a plot of ln͑X IM / X OM ͒ versus 1/T should be give straight line ͑Fig. 3͒. From this analysis, the energy difference between the two types of dimers is small. For the argon clathrates at 2 kbar, the energy difference, ͗E IM ͘ − ͗E OM ͘ is 0.08, 0.03, and 0.05 kcal/mol for SPC/E, TIP4P/Ice, and TIP4P-FQ/Ice, respectively. The values at 1 kbar are the same. For the H 2 clathrates, ͗E IM ͘ − ͗E OM ͘ is 0.05 kcal/mol. These energy differences are much smaller than the estimate of 1.08 kcal/mol from Bjerrum, 30 based just on the dimer, or the value of 0.24  kcal/mol from Pitzer and Polissar, 31 which include the nearest neighbors. It is also different than the value for the ice Ih lattice, which, for the TIP4P-FQ/Ice model is 0.08 kcal/mol. These results indicate that the relative energies of the two orientations are strongly influenced by the lattice. To a smaller extent, the type of guest influences the relative energies. Table I reports the molecule volume ͑the volume per water molecule͒ for all the clathrates studied. There is a difference of about 4% in the molar volumes among the models ͑from the 2 kbar, 250 K Ar clathrate data͒. The variations among the models are consistent with the results for ice Ih. The SPC/E model underestimates the molar volume by 3%-4%, 33,45 the TIP4P/Ice model overestimates the molar volume by 1%, 45 and the TIP4P-FQ/Ice model, by construction, gives the correct molar density. Experimentally, the molar volume is 35.57 Å 3 / molecule with an argon guest ͑at 100 K, prepared at 0.1 kbar and determined at 1 atm͒. 57 Since the experimental conditions are different, it is difficult to asses which model is most accurate, but consistent with the ice results the molar volume is largest for the TIP4P/Ice model and smallest for SPC/E. For clathrates with H 2 guests, the experimental value is 36.42 Å 3 / molecule ͑at 234 K and 2.2 kbar͒. 19 The SPC/E is within 1% of this value ͑using the nearest state point at 250 K and 2 kbar͒.

IV. CONCLUSIONS
The dielectric constant for type II clathrates appears not to depend strongly on pressure or the type of solute, whether the guest is argon or H 2 . The dielectric response arises only from the water molecules because the guest molecules do not have a dipole moment. The guests could affect the dielectric constant by influencing the fluctuations of the water molecules, but such an influence is not evident. Perhaps with guests that interact more strongly with the water molecules, the contribution from the guest would be more pronounced.
We can relate the higher dielectric constant of Ih ice compared to the clathrates from their differences in density. The ratio of the experimental dielectric constants for ice and type II clathrates at 273K ͑⑀ ice / ⑀ clathrate =94/ 58= 1.6͒ is close to the ratio of the water number molar volumes ͑V ice / V clathrate = 65.1 Å 3 molecule −1 / 38.2 Å 3 molecule −1 = 1.7͒. 2 From Fig. 2 it is apparent that the nonpolarizable models SPC/E and TIP4P/Ice underestimate the dielectric constant, while the TIP4P-FQ/Ice overestimates. Simulations of ice also reveal that nonpolarizable models underestimate the dielectric constant, 33,42 even for models such as TIP5P ͑Ref. 58͒ and TIP5P-E, 59 which have accurate dielectric constants for the liquid. For example, SPC/E gives a dielectric constant of 71 for the liquid, 60 close to the experimental value of 79, but gives a value for ice equal to 50Ϯ 4 at 200 K. 42 The nonpolarizable models predict a value for ice that is appreciably less than that of the liquid, but both polarizable and nonpolarizable models correctly predict a smaller value for clathrates than for ice Ih. Using the values at 200 K, the SPC/E model gives a ratio of the ice to clathrate dielectric constants equal to 50Ϯ 4 / 40Ϯ 4 = 1.3Ϯ .2 and the TIP4P-FQ/Ice gives 130Ϯ 16/ 90Ϯ 8 = 1.4Ϯ .2 This enhancement is less than the experimental value of 1.6, and this difference may reflect differences in the molar volumes of the various models.
The clathrate results, along with the ice results, indicate that the dielectric response may be influenced by intermolecular interactions differently for the liquid than for the solid phases. The value of the dielectric constant for the clathrates increases as the dipole moment of the model increases ͑the dipole moments of the individual water molecules in these models are 2.35, 2.426 and 3.0 D for SPC/E, TIP4P/Ice, and TIP4P-FQ/Ice, respectively͒. Most nonpolarizable models predict a dipole moment of a water molecule in ice that is too small, which is estimated to be around 3.0 D ͑Refs. 39 and 61-64͒ and likely to be greater than that of the liquid. 65 The polarizable TIP4P-FQ/Ice model gives a dipole moment of about 3.0 D and has a dielectric constant close to the experimental value for ice. 42 The analysis of the hydrogen bond order parameters reveals that the protons are disordered but show a preference for the low energy dimer arrangement ͑the OM arrangement from Fig. 1͒. The energy difference between the two is small, less than 0.10 kcal/mol, which is much smaller than the estimate of 1.08 kcal/mol from Bjerrum, 30 based just on the dimer, or the value of 0.24 kcal/mol from Pitzer and Polissar, 31 which includes the nearest neighbors. The energy difference is also different than the value for the ice Ih lattice, which, for the TIP4P-FQ/Ice model is 0.08 kcal/mol. Consequently, the relative energies of the two orientations are strongly influenced by the lattice. The clathrates are ͑slightly͒ more disordered than ice, since the energy differences are smaller and there would be a smaller energetic driving force to form the IM dimer structure. To a lesser extent, the type of guest influences the relative energies. The H 2 guests, which have a quadrupole moment, appear to increase the amount of proton disorder ͑X OM is closer to the purely random value of 8/3 for H 2 than for argon guests͒.
In addition to the results reported in this work, we have attempted to examine the dielectric constant and order parameter for the H 2 clathrate hydrates with double occupancy of the small cages. As with the single occupancy results that we have reported, we have used the same potential models used in Ref. 22. Over the time scales of our simulations, we find the lattice to be unstable with respect to dissociation at 2 kbar pressure and temperatures as low as 150 K. Whether the instability is a consequence of the model or the physics of the real system is unknown.