{'fraction_non_alphanumeric': 0.02730852279944113, 'fraction_numerical': 0.017147211990346755, 'mean_word_length': 4.846534653465347, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Thermal decomposition and diffusion of methane in clathrate hydrates from quantum mechanics simulations

Clathrate hydrates are ice-like crystalline substances in which small gas molecules are trapped inside the polyhedral cavities of water molecules. They are of great importance in both scientific research and the petroleum industry because of their applications in modern energy and environmental technologies. To achieve an atomistic-level understanding of the diffusion and decomposition of trapped molecules in clathrate hydrate, we used methane hydrates (MHs) as the prototype system and examined the methane diffusion and decomposition mechanism by employing quantum mechanics (QM) and quantum mechanics molecular dynamics (QMD) simulations. Our QMD simulations illustrated that the initial decomposition reaction in MHs initiates from hydrogen transfer among water molecules and attacks by fragments of O and OH on CH4 molecules are responsible for the destruction of the methane molecules. Next, our QM simulations revealed that the methane molecule prefers to escape from the ice cage through the hexagonal face at low temperature. To suppress the methane diffusion, we demonstrated that the diffusion barrier is significantly enhanced by adding electron or hole carriers. This is because the extra electrons and holes enhance the electrostatic interaction between methane and water molecules, leading to an increased diffusion barrier. Thus, the clathrate hydrates could be stabilized by adding extra free electron or hole carriers.


Introduction
Clathrate hydrates are ice-like nonstoichiometric crystalline compounds in which small non-polar molecules or polar molecules with large hydrophobic moieties are encapsulated in hydrogen-bonded three-dimensional network consisting of polyhedral cages of water molecules. 1,2 Clathrate hydrates are of great interest in both scientic research and industry because of their important energy applications. 3 Huge deposits of natural gas in the form of clathrate hydrates are present in petrochemical production lines and hostile environments such as marine sediments and permafrost regions. Clathrate hydrates exhibit high stability under these complex nonequilibrium environments, which is a self-preservation effect. This selfpreservation effect arises from the coating ice layer that prevents the further dissociation of the hydrate. 4 Clathrate hydrates provide new energy sources of natural gas in the future as the current accessible fossil fuel resources have become depleted. Another important application is the storage and transportation of natural gas in the form of frozen hydrates, which is considered as an alternative to conventional storage approaches such as liqueed or compressed natural gas. 1 Moreover, the CO 2 hydrates have potential applications in the CO 2 sequestration to stabilize the atmospheric concentrations of greenhouse gases. 5 Among the gas hydrates, methane hydrates (MHs) are especially attractive for many countries with limited fossil energy resources as they are abundantly detected in seaoor and permafrost zones and the estimated total amounts of carbon are twice than those in the current proved fossil fuels on earth. 6 Therefore, MHs are one of the most promising energy alternatives to the fossil energies in the future for a long period. 7,8 However, the exploitation techniques of MHs still confront many problems. Particularly, the rapid extract of MHs from the sea oor inevitably increases the local temperature, leading to fast chemical reactions and causing serious safety issues. Moreover, the increasing temperature may destabilize the structures of MHs so that large-scale methane molecules would enter atmosphere, aggravating global warming and resulting in dramatic geology problems. 6 The safety issues are also of serious concerns for the engineering applications of MHs in energy storage and transportation, and petrochemical production. In particular, the accidental re should be prevented as the natural gas hydrates are stored and transported using huge tankers. Thus, it is essential to understand the diffusion and decomposition mechanism of methane in MHs under external stimuli (e.g. heat, shock, or impact), which is crucial to prevent the spontaneous ignition and explosion of natural gas under extraction, storage and transportation.
For exploiting and utilizing the enormous energy resources contained in clathrate hydrates, it is of great importance to control the stability and combustion of clathrate hydrates, which are closely related to the diffusion and decomposition processes of trapped molecules. Many experiments have been conducted to investigate the decomposition and combustion behavior of the clathrate hydrates. [9][10][11] For example, the combustion experiments on tetrahydrofuran (THF) hydrate indicated that the ame propagates in the direction against the air ow at a constant velocity of 0.2 mm s À1 . 12 Another combustion experiment was performed to investigate the propagation of a ame above a methane hydrate cake in an air laminar boundary layer, and the ame propagation velocity was determined to be $2 mm s À1 . 13 Recently, the combustion of methane hydrates has attracted intensive attention because complete combustion of methane hydrate is difficult to achieve by direct combustion. 14 The reason is that the methane hydrate ame contains a high percentage of water vapor, leading to a low ame temperature. For example, the maximum ame temperature achieved experimentally is 1615 K with the water to methane molar ratio (H 2 O/CH 4 ¼ 1.4). 15 In addition, Hwang et al. 16 found that methane diffusion ame diluted with an appropriate amount of H 2 O may cause the enhancement of overall reaction rate. This might be because the radicals produced by water molecules, such as H, O, and OH, participate in the methane reactions and modify the methane reaction pathways from the reactants to the products. The involvement of radicals decelerates the methane oxidation reaction rate and leads to the unstable ame and the ame extinction. However, the fundamental mechanisms associated with decomposition processes of CH 4 remain not fully understood. This is due to the complex environment and coupling of thermal, chemical, and mechanical degree of freedom, which is extremely difficult to unravel experimentally.
In theoretical studies, the atomistic mechanisms of nucleation, structure stabilities, and chemical and physical properties of MHs are extensively examined by employing molecular dynamics (MD) and quantum mechanics (QM) simulations. [17][18][19][20][21][22][23][24] Previous MD simulations revealed that both water and methane molecules play an important role in the nucleation process of MHs. In particular, Guo et al. demonstrated that the 5 12 cage is able to absorb the methane molecules nearby and these methane molecules would in turn signicantly prolong the lifetime of the labile clusters. 17,18 Jacobson et al. claimed that clathrate cages only form and aggregate on the guest methane clusters, which causes the initiation of the nucleation. 19,20 Khan showed that the CH 4 molecule, stabilized in the 5 12 cages, has a stabilization energy of around 7 kcal mol À1 . 21, 22 Ida et al. showed that the proton of methane in the cavities could form the weaken hydrogen bond with four water molecules based on the electronic states of enclathrated 5 12 and 5 12 6 2 cages. 23 Besides MHs, Willow et al. showed that the presence of methane could block the released hydrogen gases in hydrogen hydrates. 24 However, many aspects remain unclear of the diffusion mechanism and decomposition reaction mechanism of methane in MHs, which is important for their energy applications.
In this work, we employed the quantum mechanics (QM) simulations to illustrate the diffusion and decomposition mechanism of methane in clathrate hydrates. Firstly, in order to examine the plausible decomposition mechanisms between CH 4 and H 2 O molecules, we carried out the quantum molecular dynamic (QMD) simulations by heating up the methane clathrate system. Then we examined the diffusion mechanism of methane in clathrate using the QM simulations. The CH 4 diffusion is expected to occur at much lower temperature than its decomposition, even before the melting of clathrate. Finally, to improve the stability of the MH, we added extra electrons or holes to the MH and examined how these free carriers inuence the diffusion mechanism and barriers. Our study provides important information to understand the thermal decomposition and diffusion mechanisms of clathrates.

Quantum mechanics simulations
All the simulations were performed using the periodic DFT code Vienna ab initio simulation package (VASP). [25][26][27] The Perdew-Burke-Ernzerhof functional and the projector augmented-wave (PAW) method were applied for the exchange-correlation interaction and the core-valence interaction, respectively. 28,29 The 2s 2 2p 2 electrons, 1s 1 electrons and 2s 2 2p 4 electrons were constructed as the valence states for the C, H, and O elements, respectively in the pseudopotentials. The electron partial occupancies were determined using the tetrahedron method with Blöchl corrections. 30 The kinetic energy cutoff was set to 500 eV for the plane-wave expansions. The convergence criteria were set to 1 Â 10 À5 eV and 1 Â 10 À2 eV$Å À1 for the electronic self-consistent eld (SCF) procedure and ionic relaxation loop, respectively. The Brillouin zone integration was performed on G-centered symmetry-reduced 1 Â 1 Â 1 Monkhorst-Pack kpoint grid meshes because of the large unit cell of the MH. We considered the van der Waals interactions using the DFT-D3 method with the Becke-Johnson damping approach. 31 We applied the Born-Oppenheimer molecular dynamics (BO-MD) approach implemented in the VASP package to examine the decomposition reactions during the cook-off simulations. Firstly, the equilibrium crystal structure was obtained at ambient conditions by heating the optimized system from 20 K to 300 K within 2 ps. Then, the system is equilibrated at 300 K for 2 ps using the NVT ensemble (constant volume, constant temperature, and constant number of atoms) with the Nose-Hoover thermostat. Finally, to accelerate the dissociation process of methane molecules, the system was heated from 300 K to 4350 K over a period of 30 ps with a constant heating rate of 135 K ps À1 . A timestep of 0.5 fs was applied for integrating the equations of motion in the QMD simulations.
To analyze the fragments in the QMD simulations, we applied the molecular fragment recognition analysis algorithm using the connectivity matrix and bond orders at 0.1 ps intervals. 32 Independent molecules are identied if their bond orders are smaller than 0.3, and then were assigned with specic ID numbers to track the reaction paths. The cutoff values of 0.3 are tested to make sure that the fragments during the reactions are appropriately described. The breaking or formation of bonds arising from thermal uctuations was avoided by a time window of 1.0 ps in which every bond must exist for at least 1.0 ps. The bond orders were computed using the reactive force eld (ReaxFF) 33 for C-H-O systems by converting the QMD trajectory to the ReaxFF trajectory.

Methane diffusion
There are three common crystal structures of clathrate hydrates, such as type sI, type sII and type sH. Among these three types of structures, the MH crystal prefers to adopt the type I structure with two small 5 12 cages and six large 5 12 6 2 cages per unit cell because of the size of guest molecules and specic temperaturepressure formation conditions. 2 Thus, we selected the sI cubic structure MH, in which one methane molecule is trapped in each cage. The unit cell of methane hydrate consists of 46 water molecules, forming two types of cages: two small (5 12 ) cages which have the shape of a pentagonal dodecahedron, and six large (5 12 6 2 ) ones which have a hexagonal truncated trapezohedron. In order to examine the methane diffusion mechanism, as well as the effects of free carriers (electrons and holes), we rst constructed a single 5 12 6 2 cage with one methane molecule inside as the initial state. Then we constructed the nal state with the methane molecule outside the cage. Finally, we carried out the climbing image nudged elastic band (CI-NEB) method to determine the transition state by nding the saddle points and minimum energy paths between the initial and nal states. 34,35 Five intermediate images were optimized along the reaction path, with the constraints implemented by adding spring forces along the band between images and projecting out the component of the force because of the potential perpendicular to the band. The spring force constant was set up to À5.0 eV A À2 . The transition state was conrmed by the frequency calculations using the density functional perturbation theory approach. 36 The migration mechanism was examined at zero temperature.
To examine the effects of electrons (or holes), we repeated the CI-NEB calculations by adding (or removing) one addition electron in the system. The extra electrons (or holes) are simulated by adding a homogeneous background-charge offset the unbalanced charge. It is possible to add or subtract a certain number of electrons even fractional charge. To illustrate the mechanism of electron (or hole) effects, we performed the analyses of charge density difference between charged and neutral systems. The VESTA soware 37 is used to visualize crystal structures and plot charge densities. In our simulations, we added electrons or holes into the neutral system, which may lead to the interaction of periodic charged images. Thus, we create a large supercell with the lattice parameter of 25 A along three perpendicular axes. Then we add one clathrate molecule into this supercell. To check the energy convergence of the simulation box, we expanded the supercell twice along one direction and added two electrons or holes to the system. The calculated electronic energies for electron and hole states in the large supercell are almost identical to the small supercell with the difference by less than 0.1%. Thus, our current model is accurate to predict the charged systems.

MH crystal structure
The unit cell of MHs was optimized using a conjugate gradients method in the VASP package with a plane-wave basis set, as shown in Fig. 1. [25][26][27] The optimized cubic cell parameters are a ¼ b ¼ 11.44 A, and c ¼ 11.43 A (at 0 K), which lead to a density of 0.92 g cm À3 , agreeing well with the X-ray experimental values of a ¼ b ¼ c ¼ 12.00 A, leading to a density of 0.81 g cm À3 at 273 K and 0.1 MPa. 2 Thus, the DFT-D3 methods can accurately describe the crystal structure of MHs.

Reaction mechanisms at high temperatures from QMD simulations
To uncover the chemical reactions of the methane hydrate, we rst examined the important fragments including HO, O and CO 2 as representatives during the QMD cook-off simulation, as shown in Fig. 2. During the initial heating-up process, from 300 to 2000 K, no chemical reactions are observed. At $2000 K, we observed the rst reaction in which the hydrogen atom is transferred between two water molecules, leading to the formation of OH and (H 3 O) + . This hydrogen transfer phenomenon was commonly observed in the hydrated materials in QMD simulations. 38 As the temperature continuously increases to $3400 K, hydrogen transfer between water molecules occurs more frequently and rapidly, producing a large number of isolated H, HO, H 3 O and O fragments. However, all CH 4 molecules remain inert during this period. This is because of the strong C-H bond with the bond energy of $98.7 kcal mol À1 . 39 To illustrate the details of how the methane molecules are decomposed, we plot key snapshots at several temperatures during the cook-off simulation, as shown in Fig. 3. The initial reaction of CH 4 molecule, occurring at 3480 K, is a bimolecular reaction in which one O 2À ion attacks CH 4 molecule, leading to the formation of OH À and CH 3 (Fig. 3a). This mutilated CH 3 becomes more chemically active than CH 4 : it rst combines with one OH fragment and forms methanol (CH 3 OH) at 3520 K (Fig. 3b); then it is attacked by another OH À ion and loses its second H atom (Fig. 3c); aer the third attacking by OH À , the carbon atom loses all of its hydrogen atoms at 3570 K and forms the second C-O bond (Fig. 3f); the intermediate CH 2 O 2 product nally turns into CO 2 at 4000 K (Fig. 3g). We only observed one CH 4 molecule converting into the nal product CO 2 . While other seven CH 4 molecules remain intact because of the strong C-H bond in CH 4 . Our QMD simulation indicated that the interactions between CH 4 and fragments of O and OH are responsible for the thermal decomposition pf methane molecules through the continuous assaults of O or OH to the carbon atom and its affiliated hydrogen atoms. In addition, in order to obtain detailed chemical reactions information, we heat up the system with a gradually linearly increased temperature from 300 K to above 4000 K. In fact, the chemical reactions would happen at a lower temperature in real experiment than our simulations here due to the long timescale, as discussed by Ye et al. 40 The initial decomposition temperature of CH 4 in my QMD simulations is 3480 K which is higher than the experimental measured ame temperature ($1500 K). 14,15 This is due to the extremely high heating rate in our QMD simulations.

Methane migration mechanism and electron excitation impaction
Methane diffusion plays an important role in the stability of MHs, as well as CH 4 exploitation from MHs. To uncover how the methane molecule escapes from the hydrate cage, we calculated the energy barrier of the escaping process from a 5 12 6 2 cage through QM simulations. Two plausible diffusion paths were considered here, as shown in Fig. 4: (I) the migration via the hexagonal face (Fig. 4b), and (II) the diffusion through pentagonal ring (Fig. 4c) by breaking and reforming a hydrogen bond of that pentagonal ring. 24 The relationship between energy and reaction coordinate (the summation of all atom movements) of these two possible pathways are calculated from CI-NEB and shown in Fig. 4a. For the hexagonal path I, the diffusion barrier is only 0.25 eV, which is much lower than that of the pentagonal path II (0.58 eV). Thus, the lowest barriers for methane diffusion is through hexagonal face, which is consistent with previous periodic calculations. 41,42 The transition state structures at both paths, shown in Fig. 4, indicate that two hydrogen bonds are broken for path I while in path II three hydrogen bonds are broken to have enough space for CH 4 diffusion. This explains the lower diffusion barrier for hexagonal path I.
The drilling in offshore regions to access the natural gas may lead to charge and discharge of clathrate hydrates. Thus, it is important to examine how electrons and holes affect the diffusion process. To illustrate how electron excitation affects the CH 4 migration behavior of MHs, we rst computed the diffusion path I at extra-electron and electron-decient (hole) conditions, respectively. Fig. 5a displays the energies for charged (and discharged) systems and the comparison with the neutral state. The e À represents the state with one extra electron, while h + stands for the state with one extra hole (or one electron decient). The diffusion barrier for one extra electron is increased from 0.25 to 0.35 eV, suggesting that adding electron will make the cage more stable. More promising, in the hole state the diffusion barrier is signicantly increased to 0.61 eV, much higher than the neutral state and e À state. Our results suggest that the storage of the CH 4 in the clathrate can be stabilized by adding electron (charging) or hole (discharging) to the MHs system. The discharging is more effective than charging in stabilizing MHs.
In order to explain the increased energy barrier by electrons and holes, we analyzed the charge density difference for the initial, transition and nal states by subtracting the charge density of the neutral system to the charged (and discharged) systems. As shown in Fig. 5b1-b3, the extra electron is distributed on both water molecule and CH 4 molecule. As the CH 4 migrates to the transition state (Fig. 5b2), the repulsion is increased by the electrostatic interaction between CH 4 and nearby H 2 O molecular, leading to an increased energy barrier. Under the hole condition, the positive charge is also distributed on both CH 4 and water molecules (Fig. 5c1-c3), leading to an increased barrier due to the electrostatic repulsion. For the e À case, the negative charges are more delocalized than the h + case, leading to a weaker electrostatic repulsion while for the h + case the positive charges are mainly concentrated in the long pair nonbonding orbitals of oxygen atoms, making a stronger electrostatic repulsion force between water and methane molecules, leading to a higher energy barrier.

Conclusions
In summary, we have performed QM and QMD simulations to provide a detailed atomic mechanism revealing the decomposition process, methane migration mechanism, as well as the inuence of the electron and hole carriers. This mechanism will be useful for understand the diffusion during the rapid extract and exploit processes and decomposition during combustion of MHs. Our major ndings are: (1) Chemical reactions in MHs initiate from the hydrogen transfer among water molecules and the reaction rates are increased dramatically with the accelerated temperature. The fragments from the water dissociation attack the methane molecules, leading to the initial decomposition reaction of CH 4 + O 2À / CH 3 + + HO À and the fast destruction of methane molecule leads to the nal product of CO 2 in our QMD simulations.
(2) The diffusion mechanism of methane molecule prefers to move through hexagonal face rather than through pentagonal side.
Charge characteristics are helpful to increase the stability of clathrate for the storage of methane molecules because of the enhanced electrostatic repulsion force between water and methane molecules. The energy barrier of extra electron (or hole) is higher than that of neutral system due to the localized partial charge distribution. The hole has a more signicant effect than the extra electron.

Conflicts of interest
The authors declare no competing nancial interests. Fig. 4 The diffusion mechanism of methane: (a) minimum energy path of hexagonal (blue) and pentagonal (black) migration path as a function of the summation of displacements of all atoms including four hydrogen ones and one carbon one. The total displacements were obtained from CI-NEB calculations; (b1-b3) the initial, transition state, and final configurations of hexagonal escape path; and (c1-c3) the initial, transition state, and final configurations of pentagonal escape path. Fig. 5 (a) Energy of electron excited (red), hole excited (green) and the ground state (blue) through hexagonal path as a function of the summation of displacements of all atoms. Here n, e À and h + represent the results for the ground state, one-electron excited state and one-hole excited state, respectively. (b1-b3) The partial charge densities for electron excited (yellow isosurfaces) migration state; (c1-c3) the partial charge densities for hole excited (blue isosurfaces) migration state. The isosurface value for plotted hole and electron densities is set to 0.004 e A À3 .