DFT-continuum characterization of third-order elasticity of sI methane hydrates under pressure

Methane gas hydrates (GHs) are polyhedral crystalline guest-host materials found under high pressure and low-temperature conditions, which can serve as an energy source. Previous work on methane GH material physics was limited to simple linear models, which only involves second-order elasticity. However, this is not fully suited to high-stress load conditions in technological applications and fundamental material physics. For other material systems, it has been demonstrated that third-order elasticity and pressure derivatives of second-order elasticity have a strong and hence significant correlation. To narrow a critical theory-simulation gap in gas hydrates materials research, in this work we expand prior work from second-order elastic constants (SOECs) to third-order elastic constants (TOECs). By using the open-source Python tool Elastic3rd and the DFT calculation software Vienna Ab initio Simulation Package (VASP), we found that the non-linear fitting involving TOECs gave a better overall prediction and a smaller root-mean-square deviation on pressure-strain evaluation when compared with linear fitting. In addition, the non-linear fitting provides robust results on the piezo-effect on the shear constant C44 and the ductile-to-brittle transition (P = −0.5 GPa). These results are not achievable from previous work based on a linear model and these findings prove that non-linear models, including TOECs, are needed under high pressures. In addition, this research includes a detailed analysis of the calculation of TOECs and mechanical properties to study pressure stability limits and ductile-brittle transitions. Together the results, findings, and analyses from this work are a novel and significant contribution to the material physics knowledge of gas hydrates and hydrogen-bonded crystalline materials.


Introduction
Extreme environments are always an important area of material science and technology investigation. These severe conditions (e.g., high pressures and low temperatures) preclude routine laboratory research. As a result, theoretical and computational modelling are indispensable to characterize material properties under high pressure-low temperature conditions. For routine material characterization, linear elastic theory has been widely applied to evaluate elastic constants/ mechanical properties. However, a wide range of research shows that the second-order elastic constants (SOECs) have intrinsic limitations due to linearity. On the other hand, higher-order elastic constants can overcome the SOECs' limitations and describe non-linear responses, which can provide more accurate high-fidelity results. Third-order elasticity has been proven that they have a strong relationship with the pressure derivative of second-order elastic constants (SOECs) [1]. It is noteworthy that TOEC calculations have already been performed on metallic materials [2,3]. Several research groups have also found that third-order elastic constants can provide reliable predictions on piezo-elastic effects on the elastic constants [1,4]. De Jong M et al [5] showed that the SOECs and TOECs can provide reliable indications on the ductile-to-brittle transitions and the failure mode under uniaxial loading. The type of mechanical instability can The organization of this paper is as follows. In section 2, all the methodologies and computational software used in this work will be discussed. Section 2.1 discusses the algorithm of the energy-strain method. Section 2.2 presents the detailed information on the DFT simulations, and sI methane hydrates' structure section 2.3 introduces three material characterization methods that were applied in this study. Sections 2.4 and 2.5 shows the equations for Born Stability criteria and Pugh's ratio and Pettifor criterion respectively. In section 3, all the results obtained for SOECs and TOECs under zero pressure will be presented. Section 3.1 demonstrates the values of SOECs and TOECs under zero pressure for the sI methane hydrates. Section 3.2 shows the relationship between hydrostatic pressures and strains and the relationship between SOECs under pressures and strains. Section 3.3 demonstrates the values of SOECs under pressure by three different characterization approaches. Section 3.4 presents the values of mechanical properties under pressures estimated by non-linear fitting involving TOECs. Section 3.5 shows the hydrate's stability evaluated by Born Stability Criteria. Section 3.6 shows the prediction on the ductile-to-brittle transition by using classical Pugh ratio and Pettifor criteria. Section 4 presents the conclusions, novelty, and significance of this work. All the technical details are given in the Supporting Information (SI), including specific examples.
We wish to emphasize that we use DFT at zero Kelvin and that the material symmetries of the gas hydrates remain cubic at all pressures, and that the studied thermodynamic phase known sI methane gas hydrate does not transform into other gas hydrates such sII or other structures through complex nucleation and growth and/or spinodal decomposition and/or complex polyhedral transformations. The methane sI hydrate has full gas occupancy and periodic boundary conditions are always applied in the simulations, representing an infinitely large material volume. Conditions, symmetries, chemical composition and processes beyond those described above are outside the scope of the present work.

Energy-strain Methodology
This section briefly introduces the basics of the energy-strain method used for the calculations of elastic constants in this paper [27][28][29]. While stress-strain analysis can be applied, the energy-strain method is preferable based on the results from previous work [14,30,31]. Within the stability region, the lattice energy E is given by an expansion of the applied strains ij h [32], as shown in equation (1). The equation (1) presented below, is expanded to the third-order, which is the maximum order of applied strain : h where 0 is the ground-state energy of the unstrained system, V 0 is the unstrained lattice volume, the strain subscripts { } ij kl mn , , represent spatial directions, C ijkl are the second order elastic constants and C ijklmn the third order ones. The symmetry properties of these latter tensors are dictated by the Voigt symmetries and lattice symmetry [33,34], which in our case is always cubic. sI methane hydrate has a cubic structure, and thus it only has three SOECs and six TOECs. The detailed definition of elastic constants and the simplification process for the cubic material are shown in the Supporting Information 1 (SI1).
To obtain the values of elastic constants ( ) C C , , ijkl ijklmn imposing cubic symmetry, using Voigt notation [33,34] and known Lagrangian strains, equation (1) can be expressed in scalar form as equation (2) [4]: E E 0 and V 0 are found from our DFT simulations, the applied strains h are known independent variables and the material property functions (A 2 , A 3 ) are the unknowns, resulting in a classical inverse problem [2,4]. In this work, the Lagrangian strain h is the applied strain for energy-strain calculation. A 2 and A 3 are linear functions of second-order elastic constants and the third-order elastic constants, respectively; see table 1 for the correspondence between A 2 and SOECs, and the correspondence between A 3 and TOECs. Table 1 is valid if equation (2) terminates at cubic terms. If equation (2) ends at quadratic terms, A , 2 then different strain sets must be applied. Details are shown in the Supporting Information (SI2). Table 1, shows the definitions A 2 and A 3 for the six selected strain modes (M1-M6) for gas hydrates, again using the Voigt notation. Six strain modes were applied because there are six third-order elastic constants for cubic structures. In this work, the applied strains h vary from −6% to 6% with increments of 2%. This range and step increase captures the non-linear response of the system and ensure the results' accuracy.
In partial summary, the energy-strain method results in an inverse problem where the unknowns are the second-and third-order elastic constants and the known values are the applied strain (h) inputs and DFTcomputed energy outputs. A representative example on calculating the SOECs and TOECs is shown in Supporting Information (SI2).
2.2. DFT simulations and sI methane hydrate structure DFT computations are performed using the VASP software [35] which can relax the lattice to its lowest-energy state. To start a DFT simulation, a unit lattice of sI methane hydrate must be constructed. Periodic boundary conditions are applied. The water molecules build the hydrate lattice, where two small cages (5 12 ) and six large cages (5 12 6 2 ) are involved in one unit lattice. The abbreviation 5 12 stands for polyhedrons that contain 12 pentagonal faces. The abbreviation 5 12 6 2 stands for polyhedrons that contain 12 pentagonal faces and two hexagonal faces. The initial atomic positions of oxygens and hydrogens were proposed by Takeuchi F et al [36]. The oxygen atoms' positions were obtained from the x-ray diffraction and hydrogen atom's positions were followed the Bernal-Fowler ice rules [37]. Methane molecules were selected to be the only guest molecules in the simulation box. As above-mentioned we assume 100% occupancy, such that each cage has a methane molecule at the centre in the initial structure. The details of structural information are summarized in Supporting Information 3 (SI3) and our previous work [14].
Besides the structural configuration, constraint parameters are fixed to specify the relaxation mode. A conjugate gradient algorithm was used, and the cut-off energy, 520 eV was selected, which was at 30% above the maximum element energy. The electronic energy tolerance used was 10 −4 eV, and the ionic optimization force tolerance used was 0.5 meV Å, in accordance with our previous work [12,23,38,39]. Projector augmented wave potential [40,41] were employed for all calculations. The revised Perdew-Burke-Ernzerhof (revPBE) exchange correlation functional with DFT-D2 dispersion correction were applied based on their great performance in our previous work [14,15,31]. We used a gamma-centered 4 4 4´mesh to ensure proper convergence with respect to the total energy for all systems [14,15]. A summary of these parameters is shown in Supporting Information 3 (SI3).

Three material characterization methods
In this paper we use, integrate and compare three calculation methods to evaluate gas hydrates' performance under arbitrary pressure: The direct numerical simulation (DNS) method uses the energy-strain method to obtain SOECs at zero pressure and under load. In other words, under zero pressure or load, the lattice system will first completely relax to its ground state. Then the strains would be applied to calculate the strained system energy. The predictions of this methods serve as a reference frame. The details of executing the energy-strain equation, shown as equation (3), and the highest order is quadratic. A 2 is the coefficient, which involves the SOECs values.
The linear model (L) uses the same energy-strain equation as DNS, equation (3), to obtain SOECs at zero pressures. Then equations (4)-(7) [42] were applied for the predictions of SOECs under pressure. The deformed strain , e is defined as the uniformed deformation in x, y, z directions due to hydrostatic pressure P: Table 1. Lagrangian strain modes ( ) Mi and the corresponding coefficients for cubic symmetry [4] (see equation (2)).

No.
Note: In M2 and M5, the elastic constant C 112 equals to C 113.
However, for the non-linear method (NL), a different energy-strain equation was applied, equation (8), for the calculation of SOECs and TOECs under zero pressure. Compared with equation (3), an extra term is involved, A , 3 which includes the TOECs: Once the SOECs and TOECs under zero pressures are obtained, the predictions of values under pressure are found using equations (9)-(12) [42]. The definition of deformed strain e is the same as in the equations shown above. Summarizing, we find: The main goal to use three different characterization methods to evaluate the SOECs values under pressures is to assess how well the linear and non-linear models perform when compared with DNS. DNS approach is considered the most reliable one because there are no fitting process applied. However, the time and computational cost for DNS are significant. Thus, if the linear and non-linear models generate accurate results of SOECs under pressure, then more cost-effective high-fidelity methods can be deployed with confidence.

Born stability criteria
Born stability criteria [43,44] are widely used for estimating a material's stability under different pressures using SOECs. Since SOECs are known to change with pressure [14], the pressure term should be included while evaluating stability. If one of these criteria is violated, it indicates material failure The simplified Born stability conditions for cubic material under load are given by [45,46]: where B T is the bulk stiffness modulus, G ¢ the tetragonal shear stiffness modulus, G and is the rhombohedral shear stiffness modulus. The entire resistance of the material to tensile and compressive pressures is determined by equation (13). Equations (14)- (15) incorporates the shear stress applied on the system.

Pugh's ratio and pettifor criterion
Two common metrics for determining a material's ductility and brittleness are the Pugh's ratio and the Pettifor criterion. The two criteria have undergone numerous revisions, but the core principle of evaluating the brittleductile transition has remained unchanged. The bulk-shear modulus ratio is determined by Pugh's ratio: B G.
/ If the ratio is more than 1.75, the material is said to be ductile, and otherwise it is considered brittle [47]. In the Pettifor criterion C 12 is compared to C 44 : C C , 12 44 where this difference is the Cauchy stress coefficients. A non-metallic system is referred to as ductile when the C 12 -C 44 > 0, and otherwise it is considered brittle [48].The goal of this paper is to integrate the two techniques in order to determine the pressure threshold for the ductileto-brittle transition in sI methane hydrate.

Results and discussion
This section will present all the data acquired using the methods described above, as well as discuss their physical significance and contribution. The predictions of various methods are identified by superscripts as defined above. Figure 1 shows the organization of the presentation, the flow of information and the connections between results given in the six subsections.

SOECs and TOECs at zero pressure
The energy-strain approach is used to determine the SOECs and TOECs of sI methane hydrates with maximum strain changes of 6%. The coefficients for the quadratic and cubic components are presented in table 2. These coefficients were obtained using six different strain modes (M1-M6). The details on the calculation methods and examples are shown in the SI2. In table 2, all quadratic components were found to be positive (A 2 > 0), while all cubic components were found to be negative (A 3 < 0). Furthermore, A 3 appears to be one order of magnitude larger than A 2 . The difference in sign between A 2 and A 3 is explained as follows. The physical meaning of ground-state and energy-strain lead to A 2 > 0. All three degrees of freedom must be relaxed to reach the ground state. Since the lattice volume, shape and atomic positions would be relaxed to their lowest-energy state first. Then the strains are applied to the relaxed system. To evaluate the energy of the strained system, the lattice volume and shape will remain unchanged. Otherwise, the strained system will return to its initial state. As a result, the energy-strain relationship is always concave-up, with minima at zero strain and as a result, A 2 > 0. On the other hand, A 3 < 0. The negative cubic components would change the sign of the product of the third-order coefficients and the cubic term . As a result, for positive strains, the estimated energy involving the cubic term The values of SOECs and TOECs for sI methane hydrates, as well as the impacts on SOECs brought on by considering the third order, will be presented in section 3.1. The established relationships between pressure and strains and between SOECs and strains will be shown in section 3.2. Then using strains as the link, we establish the relationship between pressure and SOECs with three different approaches in section 3.3. Once the pressure-constants relationship is established, we use the SOECs under pressures for the further mechanical property's calculations and stability determination in section 3.4 and section 3.5 respectively. Also, by using the calculated mechanical properties, the ductile-to-brittle transition can be determined as shown in section 3.6. will be lower than the model with only quadratic terms, and vice versa. Taking this feature into account, the lattice will attain a local maximum at a particular positive strain and subsequently decline to the lowest-energy state or even lower value, indicating system failure. When considering the negative strain side, when more stresses are applied with a negative third-order coefficient, the energy increases. This negative sign can also be rationalized by the ductile-to-brittle transition of hydrates. Previous DNS results [14] suggest that the sI methane hydrate displays a brittle behaviour when the tensile forces are applied and becomes ductile when compressive forces are applied. Thus, when the compressive pressures are applied, there would be a plastic region beyond the elastic region. Hence, the energy would accumulate. On the other side, when the tensile pressures are applied, the system would fail after reaching the threshold, which is analogous to the curve with negative coefficients for the cubic terms.
The fact that TOECs have a higher order of magnitude than SOECs can be explained by the additional strain term of the cubic terms. In the selected strain interval (−0.06 to + 0.06) the large magnitude values of A 3 can compensate the decreasing effect of strain. As shown in table 2, the A 3 values for M2 and M5 are significantly greater than for other A . 3 M2 and M5 are two strain sets that applied shear strains, which are in diagonal directions. The large values of A 3 imply that the nonlinear terms are significant in order to simulate shear deformation and emphasize the necessity of TOEC.
The exact expressions of A 2 and A 3 are given in the section 2. By using these expressions, the values of three SOECs and six TOECs can be obtained, as presented in table 3. The values of TOECs are negative except for C . 123 This can be contributed to the negativity of A . 3 Since TOECs have positive correlations with A , 3 they share the same sign. C 123 is the only TOEC that has a positive value. This exception has been noticed for other materials as well [2,4,49].
For metals and some carbonates, such as silicon, cadmium carbonate and calcite, C 144 > 0 [2,4], and for germanium and carbon, C 123 > 0 [4,49]. The alkali halides, such as LiF and NaI, have three positive TOECs 166 under 25°C [50]. The semiconductor GaAs has a positive value of C 144 under 25°C [51]. For other semiconductor, like Ge, under 25°C, the results obtained from two research groups are different. Bateman T et al [49] showed that C 123 is a positive value; however, McSkimin, H. J. and P. Andreatch [52] found all TOECs for Ge under 25°C were negative. The sign difference has been discovered for more than half a decade but does not have a physical explanation. Thus, further investigations on the sign different would be our future work.
To estimate the effects of TOEC on the prediction of SOEC NL under pressure, we compared the results with our previous data [14] of SOECs DNS . Under 0 K, the sI methane hydrate has the following values for C , 17.802 GPa, 7.231 GPa and 6.086 GPa, respectively [14]. The superscript 'oDNS' indicates results obtained at zero pressure by DNS approach. The SOECs NL are higher than SOECs DNS , and the differences range from 10% to 18.55%. This difference can be explained by the presence of TOECs. The negative TOECs will have an opposite effect on the value. Thus, greater SOECs values are required to compensate this effect.

Hydrostatic pressures as a function of strains
If the TOECs at zero pressure are known, it is possible to calculate the hydrostatic pressures P surrounding the sI methane hydrates, by imposing a deformed strain .
e Since the elastic constants may characterize the material reaction under a given stress/pressure, the surrounding pressures can be deduced from the TOECs and strains. In this section, three (DNS, L, NL) approaches described in section 2.3 are applied to evaluate the hydrostatic pressure with a given deformed strain. Firstly, the direct numerical simulation (DNS) relaxes the entire system under certain pressure. The strains can be calculated based on the deformed lattice volume and the initial lattice volume. The L and NL methods were shown in equations (4)- (9). The linear fitting used equation (4) but stopped at the first-order of e where only SOECs were involved. The nonlinear fitting used the whole equation (9), which contains firstand secondorder calculations ( ) , 2 e e and involves both SOECs and TOECs. Figure 2 presents the (DNS, L, NL) hydrostatic pressures at a given strain. DNS provides the raw simulation data, which can serve as reference. This is because, through DNS approach the lattice was fully relaxed to its lowest-energy state under a given hydrostatic pressure and the relaxed volume was calculated by VASP directly. for the DNS approach. Specifically, the hydrostatic pressure is the input and the deformed strain is the output. However, for the linear model and non-linear model, the deformed strain e is an input parameter, and the hydrostatic pressure P becomes output, which is opposite from the DNS approach. Figure 2 shows that the linear fitting underestimates the pressures at a specific strain and the nonlinear fitting always overestimates the pressures at a given strain compared with DNS. When the absolute strain values increase, the difference between linear and nonlinear approaches becomes more significant. From the curvature aspect, the nonlinear fitting has a closer nonlinear curvature with the DNS, compared with linear model. The maximum relative difference between the nonlinear and DNS under compressive pressures is 24.55%; however, the maximum relative discrepancy increases to 46.85% between linear and DNS. For the tensile pressures, the maximum relative discrepancy for both linear and nonlinear are approximately 40%. As partial conclusion, if the uniform strains e are less than 4%, there is no advantage to use nonlinear fitting for tensile pressure prediction. In other words, linear fitting is preferable under tensile pressure because it requires less computational effort. In the case where the strain is less than 3% under compressive pressure, linear fitting can produce a similar or even better forecast than nonlinear fitting for compressive loading. If the strain exceeds 3%, the nonlinear fitting should be explored for a more accurate pressure forecast. The difference between the linear and nonlinear approaches from the DNS is quantified using root-mean-square deviation estimates. The linear fitting's root-mean-square deviation is 1.25 over the whole strain range, while the nonlinear fitting's is 0.676. These findings strongly suggest that nonlinear fitting can produce a considerably more accurate estimate of hydrostatic pressures, given strain changes, than the linear fitting, especially when the strains are greater than 3%.

SOECs under pressure by three characterization methods
Elastic properties are frequently used to describe how a material behaves under different load conditions. SOECs are widely acknowledged as critical data for estimating other mechanical properties. The results of section 3.2 reveal that TOECs have a substantial advantage when estimating pressures with specified strains. Thus, this section will illustrate how well TOECs evaluate piezo-elastic effects on SOECs.
For the SOECs under pressure, C , C and C 11 12 44 can be computed using equations (5)-(7) from linear model or equations (10)- (12) from non-linear model. Engineering sensitivity analysis would be required for the examination of piezo-elastic effects, according to our earlier work [14]. The standard normalized sensitivity factor of response Q to changes in parameter P, using method m is shown below as equation (16)   For S , C 44 the DNS approach are based on the raw data C 44 with fluctuations around 7.2957 GPa under compressive pressures; the particular computational issues and existing data for gas hydrates for C 44 was previously discussed [14]. Thus, it was assumed that C 44 is not subject to compressive pressures and thus the S C DNS 44 is zero. However, it might not be the fact. S C NL 44 increase from 0 to 0.058 under pressures ranging from 0 GPa to 5 GPa. This would be more reasonable and acceptable than S C DNS 44 from the physical aspect. This is because, firstly from the curvature aspect, the positive S C showed straight lines along with pressures. From either curvature or magnitude aspects, linear fitting could not provide a good prediction of the normalized pressure sensitivity factor. Firstly, the non-linear part could not be captured by the linear fitting. Secondly, the values of S C NL DNS 11 12 44 / / / are always positive; however, S C L 11 12 44 / / has negative values along with pressures. This indicates that linear fitting is a poor tool on estimating the piezo-effects on SOECs. On the other hand, the non-linear fitting gives a very good estimation in terms of magnitude and curvature aspect.

Mechanical properties under pressure
SOECs are always served as the fundamental parameters for mechanical properties' calculations. Thus, the SOECs values under pressure become rather important on the evaluations on the mechanical properties under load. In this section, four parameters are selected to characterize for the comparisons among three approaches. All mechanical properties in terms of pressures are summarized in terms of normalized sensitivity factors. The Voigt-Reuss approximation equations used to calculate the mechanical properties are shown in the SI4. Figure 4 shows the normalized sensitivity factors values of bulk modulus (B), Poisson ratio (v), transverse wave speed (V s ), and longitudinal wave speed (V p ). The blue/orange/yellow curves stand for non-linear, linear and DNS respectively. calculated through linear fitting approach by using equations (5)-(7) and through non-linear fitting approach by using equations (10)- (12). Then the obtained results were applied into equation (10)   however; showed an increasing trend initially and converged to a constant value after 1 GPa. From these results it can be concluded that TOECs has an advantage in predicting the longitudinal wave speed. For the transverse wave speed, more studies need to be performed to rationalize the presented differences.  (12) to obtain the values of C C and C , , 11 12 44 and then applied in the corresponding equations shown in SI4. Then the parameter values under pressure equation. will are applied through equation (10) in section 2 to get their sensitivity factors. The linear fitting approach used the same methodology but used equations (5)- (7) to obtain the values of C C and C , . 11 12 44 The results show that the NL model predictions are consistent with DNS except for the transverse wave speed.

Born stability criteria
The most common material criterion for evaluating a system's stability under pressure are Born stability criteria. To access different features of pressures, such as shear stress and tensile stress, three criteria are used. The criteria have evolved into many versions over time, but the essential notion remains the same. The details are shown in equations (13)- (15).
Based on equations (13)-(15), figure 5 depicts the values of assessed terms as well as pressures. Equation (13) relates to the blue curve, which shows the bulk stiffness resistance ability. It evaluates the material's total resistance to external loads in all directions. Equation (14), shown as a red curve, is mostly used to analyze the term C 11 -C 12 , which is a different shear stress than C 44 . The shear stress term C 44 is calculated using equation (15), shown as a yellow curve. Because equation (15) is the only one that cannot hold a value greater than zero under compressive forces, the C 44 plays a substantial role under compressive pressures, as seen in figure 5. It crosses the dashed grey line (y = 0) at 7.8 GPa. This indicates that the compressive stability limits of the system have been reached. On the tensile side, there was a competition between the C 11 and the C 11 -C 12 . Furthermore, the pressures at the places where they cross the barrier are around −0.5 GPa. These phenomena, on the other hand, differs from previous observations. We previously established that the terms C 11 -C 12 had a considerable impact on both compressive and tensile stresses [14]. In the previous study, a decreasing trend of C 11 was predicted under high compressive pressures and the values of C 12 were larger than C 11 leading the term C 11 -C 12 to cross the grey dashed line on the compressive side.
To conclude, the Born stability criteria with SOECs evaluated from the non-linear method show that the compressive and tensile limits are 7.5 GPa and −0.5 GPa, respectively. Compared with the data from DNS, the non-linear fitting results show the same prediction on the compressive limits, both being 7.5 GPa. For the tensile limit, there is a small discrepancy between DNS and non-linear model. The tensile pressure limit from the DNS approach is −1 GPa [14], while the value from the non-linear model is −0.5 GPa. Nevertheless, the two tensile limits are of the same magnitude.

Ductile-to-brittle transition
To examine the fracture mechanisms of a material, it is necessary to understand the ductile or brittle properties. Two commonly used criteria are presented in this section: the Pugh ratio and Pettifor criterion. The SOECs obtained from the non-linear approach was used to determine the necessary quantities including the bulk modulus, shear modulus. To investigate the significance of non-linear fitting, we compared the results with the results from DNS approach from previous work.
The Pugh ratio and Pettifor criterion are used to assess the ductile-to-brittle transition in figure 6. The ductile-to-brittle transition threshold of Pugh ratio and Pettifor criterion are 1.75 and 0, respectively, which are shown as grey dashed lines in figure 6. If the assessed terms are greater than the threshold value, the material is considered as ductile, and vice versa. Specifically, if the B/G value is greater than 1.75 or the C 12 -C 44 value is greater than 0, the material can be recognized as ductile material. The insert of figure 6 (left) shows that the B/G value has an increasing tendency under compressive pressures and is always above the grey dashed lines. Even though the tendency and assessed term values differ between the two criteria, sI methane hydrate exhibits ductile-to-brittle transition at similar pressures. Pugh's and Pettifor's ductile-to-brittle transition threshold pressures are −0.4941GPa and −0.4806GPa, respectively. Although there is a minor difference between the two criteria, it is small enough to be considered acceptable. In comparison to prior findings [14], which showed a ductile feature under compressive pressures and a brittle feature under tensile stresses, The SOECs values under pressure from the non-linear model reflect not only general pattern and but also provide more quantitative results.
We can conclude that sI methane hydrate would be ductile under compressive forces based on this solid quantitative data. The ductile-to-brittle transition occurs under tensile pressures in the low-pressure zone, and brittle behaviour occurs when the pressure reaches −0.5 GPa.

Conclusions
This paper's objective is to characterize the elasticity of monocrystal ideal methane sI gas hydrate under pressure using linear and nonlinear models; the scope, assumptions and restrictions are specified in the introduction and emphasized throughout the paper. First, we calculated the SOECs and TOECs values at zero pressure, at zero Kelvin for sI methane hydrate through three approaches. The results show that the differences of SOECs are between 10% and 18.55%. Secondly, we investigate the relationship between the TOECs and the SOECs under pressure. The results show that the non-linear model involving the TOECs gives a closer prediction of the DNS results between the pressures and strains than the linear model. In addition, SOECs under pressure from nonlinear models in terms of normalized sensitivity factors have revealed smoother shifting patterns. The mechanical properties, e.g., bulk modulus, shear modulus, Young's modulus, and Poisson ratio, in terms of pressure, were calculated based on the SOECs from the non-linear model. The prediction of the sensitivity factor from the non-linear model can be only 0.0078 from the DNS approach. Born Stability criteria are utilised to assess the hydrate's structure stability by using SOECs under pressures from the non-linear model. It shows that the compressive and tensile limits are 7.5 GPa and −0.5 GPa from the non-linear model, which is consistent with previous work. The results are satisfactory and can provide reliable information about the pressure stability limitations. The ductile-to-brittle transition was determined using Pugh's ratio and the Pettiford criterion. The resultant transitions, at approximately −0.5 GPa, are remarkably similar across the two techniques.
The significance and contribution of TOECs can be highlighted with these excellent outcomes. Instead of simulating each pressure to determine how it affects the attributes, TOECs at zero pressure could provide all the information needed for the prediction of SOECs under pressure. Since the mechanical properties, stability determinations and ductile-to-brittle transitions are all based on the SOECs values under corresponding conditions, the ability of TOECs on the evaluation of piezo-effect on SOECs would significantly save simulation and time costs in piezo-effects research. It is noteworthy that all the findings presented here are transferable and can be applied to other hydrogen-bonded lattice structure systems.