Models

We constructed 22 different models, A1 ~ G3 are labeled, using commercial software.

Material Studio 5 was used to create the systems and MD were performed by Gromacs 4.5.5.[88, 89]

Table 3.1-1 The composition and simulation targets of model A1~G3.

Target of

There are two things should be notableFirst, the initial ice structure is zero-dipole and generated based on the rule proposed by Hayward et al. [90] Second, the zero-dipole of sI methane hydrate structure created based on Takeuchi. F et al.[91]

3.1.1 Models for melting point of sI CH4 hydrate & Ih

The two- phase model A1, as shown in Figure 3.1-1, consisting of a Ih solid phase and a pure H2O phase were used in this study to observe the melting point of Ih. The solid phase contained 6×3×3 unit cells of perfect Ih hydrate and was created following the both Hayward rules. [90] and Bernal-Fowler ice rules [92] , and the pure H2O zone totally composed of 864 molecules, which stayed in liquid phase under our simulation conditions.

The size of this model was 27.26 Å ×23.62 Å ×65.00 Å .

The four- phase model A2, as shown in Figure 3.1-2, consisting of a CH4 solid phase, a pure CH4 phase, a pure H2O phase and a H2O/CH4 mixture phase were used in this study to observe the melting point of CH4 hydrate. The solid phase contained 2×2×1 unit cells of perfect sI methane hydrate and the cavities filled with 32 CH4 molecules and created following the both Bernal-Fowler ice rules [92] , the pure H2O zone composed of 738 molecules, the pure CH4 zone composed of 90 molecules, and the H2O/CH4 mixture phase was consisted of 736 molecules and 128 molecules, respectively. The size of this model was 23.74 Å ×23.74 Å ×150.00 Å .

Figure 3.1-2 Models for melting point of CH4 hydrate crystal.

3.1.2 Models for diffusivity of H2O at different melting condition After finding out the melting point of both Ih and CH4 hydrate phase at specific pressure, the diffusivity of H2O molecules at the each of two melting point conditions is also needed. Calculating the diffusivity is not only can check the accuracy of force field, but also used in the calculation of dynamics property for capillary fluctuation theory.

The single- phase model B1&B2, as shown in Figure 3.1-3, consisting of a pure H2O phase, composed of 2112 molecules, was used in this study to observe the diffusivity of H2O at the melting point condition. The size of this model was 40.00 Å ×40.00Å ×40.00 Å .

Figure 3.1-3 Models for diffusivity of fluid-like H2O at melting condition of Ih and CH4

hydrate

3.1.3 Models for heat of fusion of Ih & sI methane hydrate

There is another property, heat of fusion of solid phase, can be also used in confirming the accuracy of force field. Moreover, heat of fusion of crystal is also another parameter used in the calculation of dynamics property for capillary fluctuation theory.

There are two sets of single- phase model C1/C2 & C3/C4, as shown in Figure 3.1-4 and Figure 3.1-5. C1 and C3 consisting of a pure solid Ih and pure CH4 hydrate crystal phase, respectively. C2 and C4 are the totally melting configuration of C1 and C3, respectively. In other words, the heat of fusion of crystal of Ih and CH4 hydrate can be obtained by comparing the difference of enthalpy energy of both C1/C2 and C3/C4 system, respectively. The size of C1/C2 model was 27.05 Å ×31.23Å ×29.44 Å , and that of C3/C4 model was 23.74 Å ×23.74Å ×23.74 Å .

Figure 3.1-4 Models for heat of fusion of Ih hydrate.

Figure 3.1-5 Models for heat of fusion of CH4 hydrate.

3.1.4 Models for mold integration method of H2O system

The one-phase model D1 ~ D3, as shown in Figure 3.1-6, consisting of a pure H2O phase was used in this study to observe the interfacial free energy of different orientation Ih/water interface. The pure H2O zone stayed in melting point condition during whole

simulation. The parameters of size and composition about this set of models are listed in Table 3.1-2.

Table 3.1-2 System size and composition for the models of mold integration Interface

orientation

Lx (Å ) Ly (Å ) Lz (Å )

H2O numbers

Mold numbers

pI 36.37 50.43 29.63 1792 128

pII 25.44 31.22 29.44 768 128

basal 36.37 31.35 47.66 1792 128

Figure 3.1-6 Models for mold integration of pI (D1), pII (D2) and basal (D3) plane of Ih/water. The snapshots of all Figures are viewed from side of normal direction of interface.

3.1.5 Models for capillary fluctuation method of Ih/water

The two-phase model E1 ~ E5, as shown in Figure 3.1-7, consisting of a pure H2O phase and pure Ih phase were used in this study to observe the fluctuation of interface between fluid/solid phase, then interfacial free energy of different orientation Ih/water interface

condition during whole simulation. The parameters of size and composition about this set of models are listed in Table 3.1-3.

Figure 3.1-7 Models for capillary fluctuation of different orientation interface of Ih/water.

Table 3.1-3 The system size and composition for the models of capillary fluctuation for Ih/water system. The Miller index in the square bracket and parenthesis indicate the direction parallel to the direction of propagation of surface wave, and the direction parallel to the normal of the interface, respectively.

Interface orientation

Ls (Å ) Lv (Å ) Ln (Å )

H2O numbers

[pI](basal) 18.151 165.242 85.805 8064

[basal](pII) 25.44 177.755 90.932 9216

[pII](pI) 29.587 181.548 90.97 15360

[basal](pII) 23.608 177.532 69.781 9216

[pI](pII) 22.190 196.688 78.776 10800

3.1.6 Models for capillary fluctuation of sI methane hydrate/water The three-phase model F1 ~ F3, as shown in Figure 3.1-8, consisting of a pure CH4 phase, H2O/CH4 mixture phase and pure sI methane phase were used in this study to observe the fluctuation of interface between water/hydrate-crystal phase, then interfacial free energy of different orientation methane hydrate/water interface can be calculated from the fluctuation data. The methane hydrate/water systems stayed in melting point condition during whole simulation. The parameters of size and composition about this set of models are listed in Table 3.1-4.

Figure 3.1-8 Models for capillary fluctuation of different orientation interface of methane hydrate/water

Table 3.1-4 The system size and composition for the models of capillary fluctuation for methane hydrate/water system.

Interface orientation

Ls (Å ) Lv (Å ) Ln (Å )

H2O numbers

[110](1̅10) 24.212 171.363 117.836 11170

[010](100) 24.091 191.105 94.186 9464

3.1.7 Models for capillary fluctuation of SDS/methane hydrate/water The three-phase model G1~G3, as shown in Figure 3.1-9, consisting of a pure CH4 phase, H2O/CH4 mixture phase and pure sI methane phase which covered by around 30% of full coverage of SDS on surface were used in this study to observe the fluctuation of interface between water/hydrate-crystal phase, then interfacial free energy of different orientation methane hydrate/water interface can be calculated from the fluctuation data. The SDS/methane hydrate/water systems stayed in melting point condition during whole simulation. The parameters of size and composition about this set of models are listed in Table 3.1-5.

Figure 3.1-9 Models for capillary fluctuation of different orientation interface of SDS/methane hydrate/water

Table 3.1-5 The system size and composition for the models of capillary fluctuation for

3.2 The setting of temperature and pressure condition

Before starting the main simulation, doing the following steps can let the system shown as Figure 3.2-1. The initial structure was first energy minimized. A short, 200 ps, NVT simulation was then conducted at 200 K to relax extra-stresses in the system. The temperature was then increased to the desired value at a rate of 0.5 K ps^-1 using NPT simulation at desired pressure. This was followed by a long, up to 10 ns NPT simulation to approach equilibrium state. The Nose–Hoover thermostat is used to maintain the temperature of the system throughout the simulation. The Berendsen barostat was applied to the third step in order to efficiently close to target pressure, and Parrinello-Rahman barostat was applied in the fourth and fifth steps to maintain the system at equilibrium. It is notable that the leap frog algorism was used for all ensemble, integrator time fixed at 1fs. PME and Lennard-Jones potential energy were used to calculate the intermolecular interaction and the long-range Coulomb energy, respectively. The cut-off of PME and LJ is 0.95nm.

For different purpose, different setting of target temperature and pressure was given

Figure 3.2-1 The flowchart of simulation.

Table 3.2-1. The composition and simulation targets of model A1~A2, B1~B2, C1~C4, D1~D3, E1~E5, F1~F3 and G1~G3