Fifth-degree elastic energy for predictive continuum stress–strain relations and elastic instabilities under large strain and complex loading in silicon

Materials under complex loading develop large strains and often phase transformation via an elastic instability, as observed in both simple and complex systems. Here, we represent a material (exemplified for Si I) under large Lagrangian strains within a continuum description by a 5th-order elastic energy found by minimizing error relative to density functional theory (DFT) results. The Cauchy stress—Lagrangian strain curves for arbitrary complex loadings are in excellent correspondence with DFT results, including the elastic instability driving the Si I → II phase transformation (PT) and the shear instabilities. PT conditions for Si I → II under action of cubic axial stresses are linear in Cauchy stresses in agreement with DFT predictions. Such continuum elastic energy permits study of elastic instabilities and orientational dependence leading to different PTs, slip, twinning, or fracture, providing a fundamental basis for continuum physics simulations of crystal behavior under extreme loading.

Notably, third-order [24][25][26] and seldom fourth-order elastic constants 27,28 are known for different crystals, as determined at small strains (e.g., 0.02-0.03). As such, fourth-order elastic constants "should be treated as an estimation only," e.g., for Si 28 . Extrapolation to large strain is unreliable to describe the lattice instability (e.g., at 0.2 for Si 10 or 0.3-0.4 for B 4 C 29,30 ). Thus, to describe correctly elasticity, including any lattice instability, higher-order elastic energies are required, and must be calibrated for a range of strain including lattice instability. For some loadings, stress-strain curves at finite strains are obtained 4,5,10,18,19,[29][30][31] , yet this is insufficient for simulation of material behavior or describing lattice instabilities under arbitrary complex loadings.
Here, fifth-degree elastic energy for Si I (cubic-diamond phase, space group Fd3m) under large strain was determined in terms of Lagrangian strains (all 6 independent components) by minimizing error with respect to density functional theory (DFT) results for large strain ranges that includes instability points. The Cauchy stress-Lagrangian strain curves for multiple complex loadings are in excellent agreement with DFT results, including elastic instability that drives the PT to Si II (β-tin structure, space group I4 1 /amd) and shear instabilities. Conditions for Si I → Si II PT under action of cubic axial stresses are found to be linear in Cauchy stresses, as predicted by DFT. Importantly, lower-order energy cannot yield similar precision in the description of stress-strain curves and elastic instabilities. Obtained elastic energy opens the possibility to study all elastic instabilities leading to different PTs, fracture, slip, and twinning, and represents a fundamental basis for continuum simulations of crystal behavior under extreme static and dynamic loading, including all the above processes and their orientational dependence.
Silicon I, as a semiconducting material, has broad applications in electronics, solar cells, and nano/micro electromechanical systems. Knowledge of rules of plasticity and transformational behavior is an important part of ensuring mechanical response and lifetime reliability in electronic devices under contact loading. The bulk hardness of semiconducting Si I phase is 12 GPa; such loads cause plastic flow and various high-pressure PTs. Machining of a brittle semiconducting Si I is accompanied by microcrack propagation inside the bulk. Utilizing strain-induced PTs to ductile metallic Si II and amorphous Si during machining allows one to realize and optimize ductile regimes of machining 32 . This also may eliminate the necessity of using chemical additives during machining, which will bring definite environmental benefits by reducing pollution.
The cubic-diamond structure of Si I can be transformed to the tetragonal β-tin structure of Si II by the transformation deformation gradient with just diagonal components, compressive −0.486 and two equal tensile 1.243, and small shuffles 10 . Thus, this transformation can occur martensitically but with large misfit strains that requires relaxation, which determines the kinetics of transformation. Under quasi-hydrostatic conditions and indentation at room temperature, the phase transformation Si I → Si II occurs in the range of 9-16 GPa 33 . Temperature T weakly affects the phase equilibrium pressure between Si I and Si II; it slightly increases or decreases the equilibrium pressure depending on different literature sources 34 . Also, increase in temperature reduces difference between the phase equilibrium pressure and pressure for Si I → Si II phase transformation. Note that there are experimentally determined data for Si I on the temperature dependence of the second-order elastic constants from 0 K up to the melting temperature 35,36 , as well as results of reactive molecular dynamics simulations 36 . The effect of temperature is relatively weak. A weak T-dependence is also expected due to the high frequencies (≈12-16 THz) of optical phonons 37 . Temperaturedependence of all elastic constants up to fifth order can be determined by utilizing ab initio molecular dynamics, see e.g., refs. [38][39][40] .
Si I → Si II transformation is irreversible. During slow unloading, Si II transforms to Si XII (rhombohedral structure r8, space group R3) and then to mixture of Si XII and Si III (bcc b8 structure, space group Ia3) under quasi-hydrostatic conditions and indentation. At rapid decompression, Si II transforms to Si IX (tetragonal st12 structure, space group P6 3 /mmc) under hydrostatic conditions and to amorphous Si after indentation. Under plastic shear in rotational diamond anvils, phase transformation Si I → Si III occurs at 3-4 GPa and Si III → Si II at 5.4 GPa 34 . Thus, large shear stresses and plastic deformation during friction, cutting, and polishing, may change the desired semiconducting Si I phase to other phases, under pressure or normal stresses, much lower than traditionally accepted 10-12 GPa. Also, conditions for semiconductor-metal transitions under complex triaxial loading were determined by DFT simulations in ref. 10 .
Due to the technological importance, the deformation and PT properties of silicon have been studied intensively. The third-order elastic constants were found with DFT 24,25 and experiments 41,42 ; however, higher-order elastic constants were not reported. The lattice instability under two-parametric loadings was studied in refs. 4,5,18,19 . Lattice instability conditions driving the Si I → II PT under action of the Cauchy stress tensor (6 independent stresses) and corresponding transformation paths were obtained in refs. 8,10 , utilizing predictions from the phase field approach 43 .

RESULTS AND DISCUSSION
Nonlinear elastic energy under large strains Motion of an elastic body is described by vector function x i (X j , t), where t is time and x i (deformed) and X j (undeformed reference state) are the Cartesian coordinates of the position vector in a natural cubic coordinate system. The deformation gradient and Lagrangian strain are then F ij = ∂x i /∂X j and η ij ¼ 1 2 ðF ki F kj À δ ij Þ, respectively, where δ ij is the Kronecker delta and Einstein summation notation is assumed. Using Voigt notation, i.e., η ii → η i (for i = 1, 2, 3), and η 23 → η 4 /2, η 31 → η 5 /2 and η 12 → η 6 /2, the specific internal energy per unit undeformed volume is, as a power-series expansion: The second Piola-Kirchhoff (PK2) stress and the true Cauchy stress are defined as We performed DFT simulations supplementing our simulations in ref. 10 , especially for shear strains and complex combined compression-shear loadings. Parameter identification procedure is carried out and results are presented in the natural cubic coordinate system. Fitting procedure Rather than determining certain set of elastic moduli from the distinct deformations 27,45 , we find all elastic moduli by the leastsquares regression using all of the DFT data we have. The error Z is a weighted sum of two terms related to the energy and PK2 stresses: Here parameters without superscript 0 designate results from approximate Eqs. (1) and (7) and those with superscript 0 are from DFT; summation over the k means summation over all sets of deformation gradients and corresponding energy u k and components of the PK2 stresses S k i , for which DFT simulations were performed, M is the number of sets of results of DFT simulations, and w is the weight factor.
Elastic moduli Fitted elastic moduli for strains η i < 0.35 are listed in Tables 1 and  2. The second-and third-order elastic moduli are also calculated for small strains η i < 0.05 for comparison with other DFT results 24 and experiments 41,42 at room temperature. Our results (while performed at 0 K) are in excellent agreement with experiments, within an experimental error and better than in ref. 24 , which validates our DFT simulations. The fourth-and fifth-order elastic moduli have no corresponding parameters from experiments and calculations to compare with. As our primary focus is on the large strain and an elastic instability, we tolerate small discrepancies between the second-and third-order constants, which we use for large strain η i < 0.35, and the corresponding values at small strains. Otherwise, the difference between stress-strain curves from the elastic energy and DFT will be worse for large strain, and the sixthorder energy will be required. Our continuum model accounts for elasticity and acoustic phonons, but ignores optical phonon modes, which become important at T > 575 K for ≈12 THz vibrations near L and X, and at T > 750 K for the 16 THz optical phonon excitations at Γ 37 .
Validation for energy Comparison of the energy contours from the elastic energy and DFT results in the plane of strains η 0 1 ¼ η 0 2 and η 3 is given in Fig. 1a (η 0 1 ¼ η 0 2 are rotated by 45°around axis 3 coordinate system, as in DFT unit cell). The stress-free Si I from elastic approximation has lattice parameters a 1 = 3.89 Å, c 1 = 5.47 Å, within 1% of DFT results (a 1 = 3.8653 Å, c 1 = 5.4665 Å), and close to the recommended value of 5.431 020 511(89) Å 46 . The saddle point (SP: η 0 1 = 0.1777 and η 3 = −0.2584) has energy 3.2976 J/mm 3 vs. 3.2893 J/mm 3 from DFT. The ability to yield the SP is crucial in capturing the elastic instabilities driving the PT. Furthermore, in Fig. 1b, the gradients of elastic energy in η 0 1 À η 3 plane (with components equal to the PK2 stresses S 0 1 ¼ S 0 2 and S 3 ) from nonlinear elastic approximation correspond well to those from DFT. Deviations between the analytical results and DFT are quite small. Note that we did not aim to fit points far from the SP toward Si II as they should be fitted to the elastic energy for Si II.
Stress-strain curves for triaxial loading We compare the Cauchy stress σ 3 -η 3 curves for different fixed lateral stresses σ 1 = σ 2 along the path toward Si I → Si II PT (Fig. 2). Corresponding transformation paths in the (η 1 = η 2 , η 3 ) plane are found iteratively using Newton method both for elastic energy and DFT simulations and are presented in Fig. 3.
It is clear from Figs. 2 and 3 that the fifth-order elastic energy captures the loading paths and stress-strain curves from DFT calculations correctly for 0 ≤ −η 3 ≤ 0.3, including peak points of the stress-strain curves, corresponding to elastic instabilities. We use the same definition as in ref. 10 : elastic lattice instability at prescribed true stress σ occurs at stresses above which the crystal cannot be at equilibrium. To determine elastic instability of Si I, nonlinear elastic properties are sufficient. However, to find the final stable or metastable state, to where the system evolves, one needs to determine positions of other local energy minima and the elastic properties of corresponding phases or states, and to model the transition process using, e.g., ab initio molecular dynamic (MD) simulations. Thus, in ref. 10 we found the PT  conditions and transformation paths for Si I → Si II PT, which will be used here. All stress-strain curves in Fig. 2 are smooth, except the one for hydrostatic loading. For nonhydrostatic loading, after instability point, elastically distorted tetragonal lattice of Si I continues transformation to tetragonal Si II. However, for hydrostatic loading, a primary isotropic deformation of cubic Si I is getting unstable with respect to a secondary tetragonal perturbation leading to Si II. Such a bifurcation of the deformation path causes discontinuity of the first derivative at the instability point. This bifurcation and jump in slope are captured correctly in Fig. 2.
Combining lattice instability points from DFT and elastic energy, we present the lattice instability criterion in the form of the critical value A of the modified transformation work: Here ε t1 = ε t2 = 0.243 and ε t3 = −0.514 are transformation strains mapping stress-free crystal lattice of Si I into stress-free lattice of Si II, and b 1 and b 3 are modifying constants. This criterion was derived in refs. 8,9,43 via phase field and was verified and quantified by both molecular dynamics simulation using Tersoff potential 8  Shear stress-strain curves and instabilities under complex loading Shear stress-strain curves for simple shears (without normal strains) and for complex loading (shear plus normal strains) are shown in Fig. 5. The elastic shear instability starts at 12.84 GPa (DFT: 12.97) for single shear, reduces to 10.7 GPa (DFT: 11) for double shear (η 4 = η 5 ), and then to 8.71 GPa (DFT: 8.56) for triple shear, below 3% error with DFT for strains beyond the instability points. Due to symmetry with respect to sign change, there are fewer nonzero elastic constants for shear than for normal strains; for single shear η 4 , c 444 = c 44444 = 0, and third and fifth degrees of η 4 , η 5 and η 6 are absent. Expectedly, deviation of elastic approximation from DFT grows for strains beyond the shear instability points much faster than for normal strains in Fig. 2. This is not critical, as for unstable branch a phase transformation occurs, which is better described by the order parameter 47,48 . In a molecular dynamics simulation 15 with a Stillinger-Weber potential 49 , the instability for simple shear along <112> in the 111 ð Þ plane and along the <111> in the 110 À Á plane lead to amorphization.
Note that double and triple shears along the <100> in the 001 ð Þ plane in Fig. 5a represent single shear in <110> in the 001 ð Þplane with η 0 4 ¼ ffiffi ffi 2 p η 4 and triaxial normal-strain loading in <111> and in the 111 ð Þ plane with η 0 1 ¼ 2η 4 and η 0 2 ¼ η 0 3 ¼ Àη 4 , respectively. Then curves in Fig. 5a can be analyzed in terms of the effect of crystallographic anisotropy. Generally, by rotating coordinate system and transforming elastic energy accordingly, one can study the effect of the anisotropy for an arbitrary complex loading.
For the shearing in combination with compressive normal strains (Fig. 5b), the DFT results are described by our elastic energy even better than just for shearing, i.e., with smaller deviation for  10 , such instabilities and paths lead to Si II. DFT (circles) and elastic energy (triangles) designate results with maximum σ 3 . The excellent agreement between elastic potential and DFT is evident. Fig. 3 Comparison of the loading paths obtained using DFT simulations and elastic energy in the (η 1 = η 2 , η 3 ) plane for compression/tension along η 3 for different fixed stresses σ 1 = σ 2 corresponding to the results in Fig. 2. The "hydrostatic" line designates loading path for σ 1 = σ 2 = σ 3 . Fig. 4 Comparison of the elastic instability condition σ 3 versus 0.5 (σ 1 + σ 2 ) for Si I → Si II PT from the fifth-order elastic energy (circles) and DFT results (*). The inset shows the same instability points from the fifth-order elastic energy (circles) in 3D space σ 1 , σ 2 , and σ 3 , which lie very close to the instability plane calibrated with DFT in ref. 10 .
larger strains even exceeding 0.35. Interestingly, superposing uniaxial compression η 1 = −0.5η 4 orthogonal to shear plane in Fig. 5b slightly increases ultimate (theoretical) shear strength but slightly reduces corresponding shear strain in comparison with Fig. 5a. At the same time, superposing uniaxial compression η 2 = −0.5η 4 in the shear η 4 direction reduces ultimate shear strength by~2 GPa, but increases corresponding shear strain. Superposing biaxial compression η 1 = η 2 = −0.5η 4 further reduces ultimate shear strength down to 6.77 GPa (6.79 from DFT) with corresponding shear strain between two previous cases. Shape of shear stress-strain curves changes also with superposition of different compressive strains. Also, superposing isotropic compression η 1 = η 2 = η 3 = −0.5η 4 = −0.5η 5 = −0.5η 6 on the triple shearing reduces ultimate shear strength from 8.71 GPa (8.56 from DFT) in Fig. 5a to 4.4 GPa (4.31 from DFT) in Fig. 5b and also strongly reduces corresponding shear strain. The tendency in reducing shear stability under hydrostatic loading in combination with presence of the dislocations with local stress concentrators may lead to pressure-induced amorphization observed experimentally 50 . The observed coupling between shear and normal stresses is very nontrivial and well captured. Typically, shear instabilities do not lead to Si II but rather to possible amorphization, hexagonal diamond Si IV, slip, or twinning.
Note that presence of the plateau-like portion in the stress-strain curves for diamond was coined in ref. 51 as "atomic plasticity" and was considered as an indicator of desired combination of high strength with sufficient ductility. Such an atomic ductility is observed for Si under compression (Fig. 2a) but not for shears (Fig. 5a). However, superposition of certain normal strains (e.g., η 2 = −0.5η 4 and especially η 1 = η 2 = −0.5η 4 ) significantly increases plateau.
Some examples and applications DFT simulations can be applied to the systems with up to 10 4 atoms with size up to 10 nm. Classical interatomic potentials used in molecular dynamics simulations treat the systems with up to 10 8 atoms with size up to 1 micron. They require calibration based on DFT simulations or/and experiments. Continuum theories and simulations are applicable to the samples with size from several nanometers (e.g., for the phase field approach to phase transformations in nanoparticles and surface-induced phenomena 52,53 ) and without upper bound. Nonlinear elasticity based on the second-and third-degree elastic constants is determined even for monolayer graphene 54,55 . Continuum theories are based on the constitutive equations describing each separate phenomena and their interactions, which are calibrated using DFT or MD simulations, experiments, and multiscale continuum or atomisticcontinuum modeling. Finite element method (FEM), finite difference and spectral methods are mostly used for solution of the physical and engineering large-scale problems with heterogeneous and complex fields.
Elastic energy in terms of components of an elastic strain tensor and corresponding stress-strain curves is the main part of any continuum theory involving mechanics. There are many cases when elastic strains are finite, i.e., exceed 0.1. For them, nonlinear, high-order elasticity should be used; in many cases, nonlinear elasticity should be used even at much smaller strains, 0.01-0.03 [24][25][26][27][28] . Under extreme loading, e.g., in shock waves 22,23 and under high static pressure 21 , volumetric strain is large and since the yield strength in shear significantly increases with pressure, shear or deviatoric elastic strains are also finite. For dislocation-free or almost free crystals, the shear stresses and strains are limited not by a macroscopic yield strength (0.1-1 GPa) but by the theoretical strength, which is one-two orders of magnitude larger. For example, diamond (the hardest material) of few micron size in a diamond anvils under pressure of 300 GPa has normal strain of 0.15 56,57 . The fourth-order elastic energy was used to model behavior of diamond anvils 21 and the third-order elastic energy for tungsten was used as a part of the total system of partial differential equations of continuum mechanics for large elastoplastic deformations under pressure up to 400 GPa. It was shown that the higher-order elastic constants significantly affect results of FEM solutions, and their values were refined by fitting to the experimental fields. The characteristic size of the diamond and tungsten sample is several mm, i.e., MD simulations cannot treat this problem. Still, equivalent stresses in diamond anvil reached 0.43 of the theoretical strength under the same biaxial lateral stresses, which means that the diamond was free of nanocracks and other strong stress concentrators.
Similarly, the third-order elasticity is used to model shock wave propagation in macroscopic samples 22,23 . Large elastic strains can also be reached in other (almost) defect-free crystals and nanoregions of a macroscopic samples that are free of defects, at coherent matrix-precipitate interfaces with large misfit strain, etc.
In the phase-field models, nonlinear elasticity is described by elastic energy in terms of elastic strains, and inelastic processes (like nucleation and motion of dislocations, twins, and cracks, and phase transformations) are described by internal variables, called order parameters 43,47,48 . As material instability and softening are involved in the models (for the second or higher order elasticity), leading to an ill-posed problem formulation and strain localization, gradient energy is introduced to regularize the problem and specify a characteristic size for the localized strain regions. For multiphase materials, elastic energy (and all other material properties) is defined for each phase independently and then interpolated between phases using order parameters. For dislocations, twins, and cracks, described by the corresponding order parameters 43,[58][59][60] , defect nucleation at the nanoscale occurs at reaching theoretical strength in shear and tension, i.e., at large elastic strains. While theory is developed for higher-order energies, the second-order elasticity is used currently 47,48,[59][60][61] , due to lack of reliable data. Utilizing higher-order elastic energy developed in the current paper will significantly improve phase fields models for phase transformations, dislocations, and fracture, and their interaction. As it is shown in the paper, our energy reproduces well DFT results and, consequently, is much better than any MD results based on a classical interatomic potential. At the same time, utilizing analytical expression for elastic energy and stresses for complex large-scale and long-time continuum simulations is orders of magnitude faster than any atomistic simulations.
For Si I, there are hundreds of publications on experiments and continuum modeling of nano-and microindentation, deformation of nano-and microspheres and micropillars, use of Si substrate for epitaxial crystal growth with large misfit strain, surface processing (polishing, scratching, cutting, etc.), and compression and torsion in diamond anvils and rotational diamond anvils. The pressure/ stress level in these processes exceeds 10 GPa, strains are finite, and stresses reach theoretical strength for nucleation of dislocations or cracks. Thus, utilizing nonlinear elastic energy found in the paper should significantly increase the accuracy of the continuum modeling. Atomistic simulations cannot treat such large samples and actual process time.
Besides, our recent results for studying lattice instability and phase transformation Si I → Si II for heterogeneous perturbations and fields 61 are entirely related to nonlinear elasticity and continuum simulations. Indeed, while for small samples that can be treated by atomistic simulations, results are usually sizedependent.
Similar elastic energy can be determined for other phases of Si, e.g., Si II, III, IX, XII, and amorphous. Due to lower symmetry, e.g., of Si II, larger number of elastic constants need to be determined, but it is doable. Elastic instability of Si II in the normal stress space along the way toward Si I was determined in ref. 10 . However, in experiment, Si II transforms much earlier to Si XII under slow unloading or to Si IX or amorphous Si under rapid decompression 32 . Thus, other types of elastic and phonon instabilities should be studied for the description of reality of transformation of Si II.
In summary, the fifth-degree elastic energy for Si I under large strain including instability points was obtained in terms of Lagrangian strains by minimizing error relative to DFT results. Elastic energy and true stress-strain curves for arbitrary complex loadings (including elastic instability) reproduce DFT results very well. Phase transition conditions for Si I → Si II under three normal cubic stresses are found to be linear in true stresses, in perfect agreement with DFT. Any lower-order energies (less than fifthdegree) cannot derive a similar precision in description of elastic instabilities and stress-strain curves, whereas, in contrast, they are currently found mostly using third-order elastic constants determined at small strains. Our results also show the potential of controlling the stress-strain curves and phase transitions by applying optimized, multidimensional loading to control desirable properties and to drastically reduce phase transition pressures (1-2 orders of magnitude) 10,17,31,62 .
Besides being generally applicable, the elastic energy contains in convenient analytical form a plethora of information and now permits a direct continuum study of all elastic instabilities under complex loading driving different phase transitions (allotropic and amorphization), fracture, slip, and twinning. Using higher-order energies and large strains that include instabilities yields qualitatively and quantitatively better predictive capability, improving the entire continuum model-based simulations, which are much faster than DFT and molecular dynamics. Notably, our approach represents a fundamentally new basis for continuum simulations of crystal behavior under extreme static and dynamic loadings involving multiple the above mentioned orientationaldependent mechanisms. In particular, higher-order elasticity is required for determination of the stress-strain states and optimization of the diamond anvil cell for reaching maximum possible pressures 21 . This approach is general and will significantly improve phase-field models for phase transformations, in contrast to the second-order elasticity used currently 47,48,61 . It also provides a basis for the description of the competition between different instabilities at different loadings.

METHODS
We used DFT as implemented in VASP [38][39][40] with the projector augmented waves (PAW) basis 63,64 and PBE exchange-correlation functional 65 . The PAW-PBE pseudo-potential of Si had 4 valence electrons (s 2 p 2 ) and 1.9 Å cutoff radius. The plane-wave energy cutoff (ENCUT) was 306.7 eV, while the cut-off energy of the plane wave representation of the augmentation charges (ENAUG) was 322.1 eV. We used a Davidson block iteration scheme (IALGO = 38) for the electronic energy minimization. Electronic structure was calculated with a fixed number of bands (NBANDS = 16) in a tetragonal 4-atom unit cell (a supercell of a 2-atom primitive cell). Brillouin zone integrations were done in k-space (LREAL = FALSE) using a Γcentered Monkhorst-Pack mesh 66 containing 55-110 k-points per Å −1 (fewer during atomic relaxation, more for the final energy calculation). Accelerated convergence of the self-consistent calculations was achieved using a modified Broyden's method 67 .
Atomic relaxation and energy minimization in a unit cell (ISIF = 2) fixed by the different prescribed values of the components of the deformation gradient tensor F ij (for more than 10 4 different combinations of F ij and corresponding E ij ) were performed using the conjugate gradient algorithm (IBRION = 2), allowing symmetry breaking (ISYM = 0). The transformation path was confirmed by the nudged-elastic band (NEB) calculations, performed using the C2NEB code 68 . We used DFT forces in ab initio molecular dynamics (MD) to verify stability of the relaxed atomic structures. Si atoms were assumed to have mass POMASS = 28.085 atomic mass units (amu). The time step for the atomic motion was set to POTIM = 0.5 fs.
Convergence vs. plane-wave energy cutoff DFT codes, including VASP, typically return energy with a fairly high precision, but have larger errors in stress components calculated using the Hellmann-Feynman theorem 69 within a less precise linear-response method. We avoided this issue by calculating the energy versus deformation on a grid and used exclusively finite differences to find derivatives of energy with respect to deformation and stress components.
As an example of convergence versus plane-wave energy cutoff (ENCUT) for the structure with a = b = 4.1279 Å and c = 4.4638 Å, identified as the stress barrier under uniaxial loading at σ 1 = σ 2 = 0, we fixed a = b, varied c by ±1%, obtained the finite-difference derivative of energy dE/dc ≈ [E(c + δ) − E(c − δ)]/[2δ], and plotted σ 3 versus ENCUT (E cut ) in Fig. 6. The chosen plane-wave energy cutoff of 306.7 eV is sufficient to achieve convergence within ±0.1 GPa (1 kBar) for the finite-difference method, which we use. The converged DFT data contained over 10 4 entries and was processed in the format, outlined in Table 3 in ref. 70 H. Chen et al.
The transformation formulas of the deformation gradient and Cauchy and second Piola-Kirchhoff (PK2) stresses to the natural cubic coordinate system are Parameter identification procedure is carried out and results are presented in the natural cubic coordinate system.

DATA AVAILABILITY
Data sets for all figures and all raw data from DFT simulations are available at https:// doi.org/10.25380/iastate.12668843.
Received: 5 April 2020; Accepted: 2 July 2020; 1000 E cut (eV) 11 11.2 11.4 -3 (GPa) Finite difference Linear response Fig. 6 Convergence of stress component σ 3 obtained from energy derivatives calculated using finite-difference and Hellmann-Feynman linear-response methods. We use exclusively the finite differences. The chosen energy cutoff of 306.7 eV (vertical dashed green line) is sufficient to achieve convergence within ±0.1 GPa (horizontal dashed lines) within the finite-difference calculations.