Performance of SCAN Meta-GGA Functionals on Nonlinear Mechanics of Graphene-Like g-SiC

Although meta-generalized-gradient approximations (meta-GGAs) are believed potentially the most accurate among the efficient first-principles calculations, the performance has not been accessed on the nonlinear mechanical properties of two-dimensional nanomaterials. Graphene, like two-dimensional silicon carbide g-SiC, has a wide direct band-gap with applications in high-power electronics and solar energy. Taken g-SiC as a paradigm, we have investigated the performance of meta-GGA functionals on the nonlinear mechanical properties under large strains, both compressive and tensile, along three deformation modes using Strongly Constrained and Appropriately Normed Semilocal Density Functional (SCAN) as an example. A close comparison suggests that the nonlinear mechanics predicted from SCAN are very similar to that of Perdew-Burke-Ernzerhof (PBE) formulated functional, a standard Density Functional Theory (DFT) functional. The improvement from SCAN calculation over PBE calculation is minor, despite the considerable increase of computing demand. This study could be helpful in selection of density functionals in simulations and modeling of mechanics of materials.

Strain engineering is a common approach to tune a material's properties, especially two-dimensional nanomaterials due to their outstanding flexibility [19]. In general, a nanomaterial is vulnerable to strain [20][21][22][23][24][25], which could be induced by thermal vibration, surface adhesion, epitaxial grown, substrate deformation, thermal expansion mismatch, prestretched substrate, pressurized blisters, substrate topography modification, tip indentation, and many other process. It was predicted that the indirect-direct band gap transition of the SiC nanosheet can be tuned through changing the lattice constant [3] which is equivalent to biaxial straining. Control of the electronic and optical properties of multi-layer g-SiC sheets by strain engineering up to 9% strain was also reported recently [26]. The knowledge the mechanical properties is crucial to the applications of these nanomaterials [27][28][29].
We have investigated the non-linear mechanical behaviors by means of first-principles calculations within the framework of Density Functional Theory (DFT) [30]. DFT is a feasible first-principles method with less computationally demanding than other computational methods with a similar accuracy. Therefore, DFT is the most widely-used method of electronic structure calculations. The biggest challenge of DFT is the description of the exchange-correlation term, which represents the correction to the kinetic energy and the electron-electron repulsion energy [31]. There are many approximations to the exchangecorrelation energy, with different level of accuracy, which are described by various rungs of "Jacob's Ladder" [32]. The first three rungs are the local density approximation (LDA), generalized gradient approximation (GGA), and meta-GGA approximation (MGGA) [33]. Among them, Perdew-Burke-Ernzerhof (PBE) formulated GGA [34] is widely used as a "standard" of DFT calculations for its combined good performance of accuracy and computation demand. Recently, Strongly Constrained and Appropriately Normed Semilocal Density Functional (SCAN) rises up with promising improvement of accuracy whereas not-much increasing of computing demands [35,36]. There is a concern about the precision of first-principles modeling of the nonlinear elastics at large strain. In general, it is believed that the selection of exchange-correlation functionals affect the prediction precisions in DFT calculations. To which extent, however, is still obscure.
In this work, we aim to assess the performance of the rungs of DFT exchangecorrelation approximations on the outcome precision of nonlinear mechanics of 2D materials. We focus on the GGA (second rung) and the MGGA (the third rung) for g-SiC , taking PBE-GGA and SCAN-MGGA as an example. We deform the g-SiC mono-layers to various large strains and compute the high order (up to 5th order) elastic constants. These elastic constants are important for an accurate nonlinear elastic description. Furthermore, the pressure dependence properties are predicted. These properties include sound wave speed (both p-wave and s-wave) and the second order elastic constants. Our method is introduced in the Section 2. Our results and discussion are presented in Section 3, followed by Conclusions in Section 4.

Computation Method
We adopted the model of a 6-atom conventional unitcell, which has 3 carbon atoms and 3 silicon atoms, as shown in Figure 1. This large unitcell is relatively larger than the primitive unitcell which only contains two atoms (one carbon and one silicon). Our model is capable to capture the soft mode in hexagonal rings [37]. We performed the DFT calculations using the program of Vienna Ab-initio Simulation Package (VASP) [38]. We adopted the PBE-GGA [34] and SCAN-MGGA [35] for exchange-correlation functions in DFT calculations in this study. We denoted them as "PBE" and "SCAN", respectively, in the rest of text for convenience of description. We have employed projector augmented wave (PAW) and pseudo-potential approach [39] for modeling electrons. For carbon atoms, the 2s 2 2p 2 electrons are explicitly considered. For silicon atoms, the 3s 2 3p 2 electrons are included. The van der Waals interactions are treated with DFT+q method [40]. The periodic boundary conditions are applied along all three directions to avoid the surface effect. A cutoff of energy of 600 eV for the basis of plane waves is used in all the DFT calculations. The criterion of electronic relaxation is the convergence of total energy within a fluctuation less than 1.0 × 10 −5 eV. The ionic relaxation is achieved when the atomic forces are less than 0.001 eV/Å. The irreducible Brillouin Zone of the unitcell is sampled with a Gammacentered 24 × 24 × 1 k-mesh. A 15 Å thick vacuum region is included to model the single layer structure. The second Piola-Kirchhoff (K-P) stress [27] is used for the 2D forces per length with units of N/m. The mechanical deformations are applied in thee modes: uniaxial along the armchair direction (i.e., nearest neighbor direction), uniaxial along the zigzag direction (i.e., second nearest neighbor direction), and biaxial. Both tensile (namely positive strain) and compressive (i.e., negative) strains are examined in all the three modes. The fourteen non-zero elastic constants of g-SiC are obtained by a least-squares fit to stress-strain results [41]. This method has been widely used to study the mechanical properties of various graphene-like materials, including g-BN, g-GeC, g-GaN, g-TlN, g-ZnO, g-Boron, g-BNC, g-ZnS, g-MoS2, graphene-oxide, and g-InN [23][24][25]29,[42][43][44][45][46][47][48][49][50][51][52][53][54].

Atomic Structure
We have fully relaxed the system with all degree of freedom including both the cell geometry and individual atoms. All atoms of a g-SiC mono-layer are coplanar, same as that of graphene. The lattice constant of a = 3.103 Å and 3.096 Å for SCAN and PBE, respectively. These values are in an excellent agreement with previous DFT studies [1][2][3]. The bond length of C-Si bond is 1.792 (1.787) Å from SCAN (PBE) calculations, which agrees well with previous DFT calculations of 1.77 Å [2] and 1.79 Å [1]. The bond angles of C-Si-C and Si-C-Si are 120 • . Our result of atomic structure of g-SiC agrees well with literature [1][2][3]. The lattice constant of g-SiC mono-layers is compared to graphene and silicene as summarized in Table 1.

Strain Energy Profile
The energy profile is the stored potential energy under applied strain. The energy profile illustrates the change of energy with respect to strain, which can be used to describe the mechanical stability. When a strain is mechanically loaded, all the atoms in the supercell are allowed to relax. We have studied both the compression (negative strains) and tension (positive strains). The Lagrangian strains are used due to describe the applied large strains. For PBE calculations, the range is applied strain is −0.1 to 0.4. The tensile strain is up to 0.3 for SCAN calculations because the critical strains are less than 0.3. The step size of applied strain is 0.01 all three deformation modes. The reduced strain energy is defined as E s = (E tot − E 0 )/n for convenience of comparison between different systems including graphene and silicene [25], where E tot is the total energy of the strained system, E 0 is the total energy of the strain-free system, and n = 6 is the number of atoms in the unitcell. The strain energy profile (the strain energy as a function of the strain) of g-SiC are displayed in Figure 2. We conclude from the comparison in Figure 2 that SCAN predicts similar behaviors of strain energy profiles for uniaxial armchair, uniaxial zigzag and equibiaxial deformations to that of PBE. The critical strains (denoted by dashed lines) in strain energy profile are the ultimate tensile strains, where the second derivative of the strain energy to the strain is zero. The critical strains of the three deformation modes are summarized in Table 1 for comparison. The armchair critical strain from SCAN calculation is identical to that of PBE calculation. However, the zigzag and biaxial critical strain of SCAN is 0.01 larger than that of PBE. The difference is trivial because the error bar in strain is 0.01, same as the minimum increment strain step in this study.

Stress-Strain Curves
The second P-K stress is used due to the two-dimension nature of the system. The stress-strain relationships of g-SiC are shown in Figure 3 for three deformation modes, which are uniaxial strains along the armchair and zigzag directions and biaxial strains. The ultimate tensile strength is the maximum stress of the stress-strain curve, which means the upper limit that a material can withstand in a tensile test. The corresponding strain is named the ultimate tensile strain. Therefore, the system under strains beyond the ultimate strains is unstable. The ultimate strain serves as an upper limit of elastic stability. As a result, the ultimate tensile strength is a key parameter for strain engineering and applications. As listed in Table 1, the ultimate tensile stress at the three deformation modes are compared to that of graphene [45] and silicene [49].
The ultimate strains of SCAN are similar to that of PBE. However, the ultimate tensile stress of SCAN are larger than that of PBE, which is 6.7%, 9.0%, and 9.9% for armchair, zigzag, and biaxial deformation, respectively.

High Order Elastic Constants
The stress-strain relationship is the central part of the elastic theory of continuum mechanics, where the elastic constants are the key parameters. These elastic constants are critical to predict the mechanical behaviors at macro-scale. Among them, the high order elastic constants are essential in nonlinear elastic mechanics. We have predicted the 14 non-zero elastic constants from SCAN calculations, which includes normal and high order elastic constants. There are 2 SOECs (second order elastic constants), 3 TOECs (third order elastic constants), 4 FOECs (fourth order elastic constants), and 5 FFOECs (fifth order elastic constants) as listed in Table 1. The SOECs is for linear elastic mechanics and the high order elastic constants (TOECs, FOECs, and FFOECs) are for nonlinear elastic mechanics. Our results of elastic constants and the modulus from SCAN calculations are compared to that of PBE calculations as summarized in Table 1, along with that of graphene and silicene. It is worth noting that the in-plane Young's modulus Y s are defined as Y s = (C 2 11 − C 2 12 )/C 11 and Poisson ratio is computed as ν = C 12 /C 11 .

Pressure Effect on Elastic Constants
One of the applications of the high order (>order 2) elastic constants is to describe the pressure effects on mechanical properties. The second-order elastic moduli of g-SiC are not constants when the in-plane pressure is introduced. Instead, second-order elastic moduli are a function of in-plane pressure. The pressure-dependent second-order elastic moduli (C 11 ,C 12 ,C 22 ) can be obtained from second-order elastic moduli at zero pressure (C 11 , C 12 , C 22 ,Y s , and ν), and third order elastic constants (C 111 , C 112 , and C 222 ) as [30,41]: The results of the pressure effects are illustrated in Figure 4 for both SCAN and PBE calculations. The predicted second-order elastic moduli, in-plane Young's modulus, and Poisson ratio as a function of the in-plane pressure from SCAN are similar to that from PBE calculations.

Pressure Effect on the Soundwave Speed
One of the applications of the high order elastic constants is to predict the pressure effects on the speed of soundwaves propagating in the g-SiC. We have investigated the speed of anharmonic sound waves as a function of the in-plane pressure. The pressure-dependent soundwaves are interesting in many applications. Depending on the deformation mode during propagation, there are two main classes of soundwaves: compressive waves or p-waves, and shear waves or s-waves. With the pressure-dependent elastic constants, the speed of these two types of soundwave can be computed as The dependences of both u p and u s on in-plane pressures, are plotted in Figure 5. It is well demonstrated that the speed of both p-wave and s-wave predicted from SCAN calculation is very similar to that of PBE.
Besides the speed along, the velocity ratio of the two, u p /u s , is a key parameter in the mechanics of materials [41]. The velocity ratio has been used for many purposes, such as a lithology indicator, identifying pore fluid, determining degree of consolidation, and predicting velocities [55]. The velocity ratio relates to the Poisson's ratio as It is well demonstrated in Figure 5 that the velocity ratio predicted from SCAN is similar to that of PBE calculations.

Elastic Stability
We have also examined the elastic stability according to the Born elastic stability criteria [56]. The necessary and sufficient elastic stability conditions are illustrated in the following equation reflecting that the strain energy density is positive for a strain [57]: Our results of SCAN calculation of strain energy profiles shown in Figure 2 implies the elastic stabilities of g-SiC mono-layers. The same conclusion can be drawn from PBE calculations.
In terms of elastic stiffness tensor, furthermore, we examine elastic stability from the aspect of elastic stiffness tensor. The criteria are expressed in the following equations for a 2D hexagonal lattice: Both the SCAN and PBE results of the second order elastic constants summarized in Table 1 suggest the satisfactory of the stiffness tensor criteria (Equation (8)) for g-SiC mono-layers.

Conclusions
We have assessed the performance of the SCAN calculations on the nonlinear mechanical properties and mechanical stabilities of g-SiC under large strains with point-by-point comparison to that of PBE calculations, including ultimate strength, ultimate strain, 14 non-trivial elastic constants up to 5th order, strain energy profile, stress-strain behaviors, and pressure effect on normal modulus, shear modulus, Poisson ratio, p-wave, s-wave, and velocity ratio, under large deformations in three modes. We demonstrate that results from SCAN are similar to that of PBE calculations. The armchair ultimate strain from SCAN is identical to that of PBE. The difference in zigzag and biaxial ultimate strain is only 0.01, within the precision of strain we applied in this study. The differences in ultimate stress are all less than 10%, which fall in the precision of a DFT calculations.
The agreement between SCAN and PBE calculations is strikingly good in prediction of the pressure on the normal modulus, shear modulus, Poisson ratio, p-wave, s-wave, and velocity ratio. The results of the elastic stability from SCAN and PBE are also identical. From these comparisons, we could conclude that the SCAN calculations have very similar performance in predicting the nonlinear mechanics. In other words, the improvement in precision in predicting nonlinear mechanics from SCAN calculation over PBE calculation is minor. It is worth noting that the computing demand of SCAN calculation is about double that of PBE calculation. Considering the considerable increase of computing demand, the improvement of the precision seems not worthwhile. The PBE results are still reliable in predicting mechanical properties as the standard of DFT calculations. Data Availability Statement: The datasets generated during and/or analyzed during the current study are available from the corresponding author on request.

Conflicts of Interest:
The authors declare no conflict of interest.