The performance of F4 order parameter

There was not any research using F4 order parameter to describe the population distribution of every each oxygen atom of whole system.[97, 98, 106, 107] In order to well distinguish the CMI position of the system, we have to carefully select the optimize value of F4, the result of F4 value of oxygen atoms distribution are shown in Figure 4.2.1(a)~(d). The average value of the F4 distribution that did not do any block-average analysis (Figure 4.2.1-(a)) for Ih, water and sI hydrate is -0.45, 0, and 0.86, respectively.

These values correspond to data published, indicating that the distribution profile really reflects to the system.[97, 98, 106, 107] Figure 4.2.1 shows that if the block size below 20 (100 ps), the overlap area for Ih-water and sI hydrate-water is very large and the degree of the overlapping of Ih-water is much large than sI hydrate-water (see Figure 4.2.1-(b)).

The reason we do not want to see the overlapping area is because it is hard to use F4 order parameter to differentiate the crystal and fluid-like H2O at this condition. There will not be a critical F4 value can discriminate the difference of material. Finally, we decide to use the 40 of block size (200 ps) as our baseline of resolution. However, there is a main problem exist – Does the fluctuation of wave be averaged if using block-average technique or will the block-average generate the artificial CMI? We will discuss this problem in the last part, but we can spoil the answer at this time, the answer is – Block-average only affects the Rayleigh wave of the CMI, the capillary wave will not be affected

Figure 4.2-1 (a) F4 probability distribution of the system with no block-average. (b) block-average size is set as 4. (c) block-average size is set as 20. (d) block-average size is set as 40. The line colored by azure, orange and silver denote Ih, water, and sI hydrate, respectively.

4.2.2 Thermodynamics properties of CMI system of Ih

We have already shown that using F4 to find out the location of CMI is adequate.

Then we input the h(vn,t) of each orientation system into the Fourier-transform formula,

not smooth interface, by representing ln[ <|hq|²>A/(kBT) ] versus ln(q), according to Equation.(37), we should obtain a straight line with -2 slope and the intercept – ln(̃).

The data profiles are shown in Figure 4.2.2 for every orientation also be focused in E.Sanz et al.[36] The blue solid-circle in Figure 4.2.2 correspond to our simulation data and red solid-line denotes a straight fitting line with -2 slope at low-q region.

Obviously, our simulation data points in low-q region are perfectly described by equation (37). These data points allow us to calculate the reasonable ̃ from intercept, and the result of calculation are shown in Table 4.2-1. We only use first 4 points, determined by observing discrepancy level, in fitting for every orientation.

Once the stiffness is known for a set of different orientations, we can get the interfacial free energy by solving Equation (38). In order to solve equation (38), the relationship between  and normal vector direction 𝑛̂ should first determined. Since the point group of Ih is 6/mmm (D6h), the  and 𝑣̂ can be written as an expansion in terms of Spherical Harmonics [108, 109] :

(,) = ₀[1 +₂₀𝑦₂₀(,) +₄₀𝑦₄₀(,) +₆₀𝑦₆₀(,) +₆₆𝑦₆₆(,)] (27) where ₀ is the interfacial free energy averaged over all orientations, , are spherical coordinate defined in Figure 4.2-3 and lm are anisotropy parameters.

The spherical ylm(,) represent normalized real harmonics functions taking into the form shown in Table 4.2-2, and the results of applying the parameters , into the Table 4.2-2 are listed in Table 4.2-3. Using Table 4.2-2 and Table 4.2-3, the equation (24) could be solved for interfacial free energy. The expression of expansion of spherical harmonics for the orientation in Table 4.2-1 are listed in Table 4.2-4. The results are 3.963, 0.4297, -0.3575, 0, and infinite for 20, 40, 60, 66, and 0, respectively.

Unfortunately, the equations of Table 4.2-4 are not totally linear independent, in other

words, the four anisotropy parameter and 0 cannot be solved. There are same problems indicated in ref. [109] and E.Sanz’s work. Therefore, in both of their work, they treated

40 and60 as zero because these two factors are too small comparing to 20, 66. Finally, we can accurately fit the remaining parameters with this strategy by using the equations listed in Table 4.2-4, and 20 = -0.027,66 = -0.0027 and 0 = 29.24(0.91) mN/m. The 0

value of TIP4P/Ice calculated in this method is in good agreement with both other simulation method and experiment data.[47, 65, 110] Since we also validate this

simulation results by other simulation method, all simulation works including this work were sorted out in section 4.2.5. The interfacial energy of different orientations have also been calculated shown in Table 4.2-5, and the basal plane is the lowest energy plane, same as the experimental result.[111]

Figure 4.2-2 Plots of ln[ <|hq|²>A/kBT ] versus ln(q) for all Ih/water with different orientations studies. The unit of <|hq|²>A/kBT is used in m³*N^-1 and q is given in m^-1. The red solid-lines shown in (a) ~ (e) are linear fits with slope -2 at low-q region. The intercept of the fitting is –ln(̃), and the unit of ̃ calculating by intercept is N/m.

Table 4.2-1 the interfacial stiffness (unit: mN/m) of each orientation for Ih/Water.

Orientation

This work E.Sanz (Tip4p/2005)

[pI](basal)

28.10 28.04(0.43)

[basal](pI)

26.57 25.30(2.45)

[pII](pI)

30.19 29.01(0.15)

[basal](pII)

33.37 33.70(1.55)

[pI](pII)

27.13 26.83(0.42)

Figure 4.2-3 Coordinate system of Ih, , are spherical coordinate angles used to define the orientation of normal vector direction (𝑛̂) to CMI.[109]

Table 4.2-2 The expression of normalized spherical harmonics function [108]

y₂₀(,) = √ 5

Table 4.2-3 Interfacial free energy in terms of normalized spherical harmonics for different orientations

Table 4.2-4 Stiffness in terms of normalized spherical harmonics for different orientations

Table 4.2-5 Interfacial free energy of different orientation of Ih

This work (mJ/m²)

Basal 28.73(1.00)

pI 29.55(0.88)

pII 29.44(0.86)

4.2.3 Thermodynamics properties of CMI system of CH4 hydrate