Fifth

Jun 13, 2024

npj Computational Materials volume 6, Article number: 115 (2020) Cite this article

2310 Accesses

13 Citations

1 Altmetric

Metrics details

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.

Nonlinear, anisotropic elastic properties of single crystals determine material response to extreme loading, e.g., in shock waves, under high static pressure, and in defect-free crystals and nanoregions. Elastic nonlinearity ultimately results in elastic lattice instabilities1,2,3,4,5,6. Such instabilities dictate various phenomena, including phase transitions (PT, i.e., crystal-crystal7,8,9,10, amorphization11,12,13,14,15, and melting16,17), slip, twinning, and fracture, in particular, theoretical strength in tension, compression, or shear3,4,5,18,19,20. Besides, nonlinear elastic properties are necessary for continuum simulations of material behavior under extreme static21 or dynamic22,23 loadings and near interfaces with significant lattice mismatch.

Notably, third-order24,25,26 and seldom fourth-order elastic constants27,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 Si28. Extrapolation to large strain is unreliable to describe the lattice instability (e.g., at 0.2 for Si10 or 0.3–0.4 for B4C29,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 obtained4,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 I41/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 machining32. 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 shuffles10. 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 GPa33. 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 sources34. 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 temperature35,36, as well as results of reactive molecular dynamics simulations36. The effect of temperature is relatively weak. A weak T-dependence is also expected due to the high frequencies (≈12–16 THz) of optical phonons37. Temperature-dependence 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 \(R\bar{3}\)) and then to mixture of Si XII and Si III (bcc b8 structure, space group \(Ia\bar{3}\)) under quasi-hydrostatic conditions and indentation. At rapid decompression, Si II transforms to Si IX (tetragonal st12 structure, space group P63/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 GPa34. 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 DFT24,25 and experiments41,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 approach43.

Motion of an elastic body is described by vector function xi(Xj, t), where t is time and xi (deformed) and Xj (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 Fij = ∂xi/∂Xj and \({\eta }_{ij}=\frac{1}{2}({F}_{ki}{F}_{kj}-{\delta }_{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:

where the c’s are elastic moduli of second, third, fourth, fifth and higher order. For crystals with cubic symmetry, Eq. (1) is specified in cubic axes, with 3 second-, 6 third-, 11 fourth-, and 18 fifth-order moduli44, found here using DFT. Thus, with assuming the zero energy for strain free case, i.e., u0 = 0, and

in which

All nontrivial terms in Eqs. (3)–(6) are determined in the same way: general polynomial expressions are taken as the starting point and then all known symmetry operations for the cubic lattice are applied.

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.

Rather than determining certain set of elastic moduli from the distinct deformations27,45, we find all elastic moduli by the least-squares 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 uk and components of the PK2 stresses \({S}_{i}^{k}\), for which DFT simulations were performed, M is the number of sets of results of DFT simulations, and w is the weight factor.

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 results24 and experiments41,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 sixth-order 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.

Comparison of the energy contours from the elastic energy and DFT results in the plane of strains \({\eta }_{1}^{\prime}={\eta }_{2}^{\prime}\) and η3 is given in Fig. 1a (\({\eta }_{1}^{\prime}={\eta }_{2}^{\prime}\) 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 a1 = 3.89 Å, c1 = 5.47 Å, within 1% of DFT results (a1 = 3.8653 Å, c1 = 5.4665 Å), and close to the recommended value of 5.431 020 511(89) Å46. The saddle point (SP: \({\eta }_{1}^{\prime}\) = 0.1777 and η3 = −0.2584) has energy 3.2976 J/mm3 vs. 3.2893 J/mm3 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 \({\eta }_{1}^{\prime}-{\eta }_{3}\) plane (with components equal to the PK2 stresses \({S}_{1}^{\prime}={S}_{2}^{\prime}\) and S3) 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.

a Fifth-order elastic energy and b energy gradients in \({\eta }_{1}^{\prime}={\eta }_{2}^{\prime}\) – η3 plane. Components of gradients are PK2 stresses \({S}_{1}^{\prime}={S}_{2}^{\prime}\) and S3.

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.

From ref. 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.

The “hydrostatic” line designates loading path for σ1 = σ2 = σ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 b1 and b3 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 potential8 and DFT simulations10. Instability lines can be approximated by σ3 = 0.4144(σ1 + σ2) − 10.9121 for nonlinear elasticity and by σ3 = 0.4066(σ1 + σ2) − 11.4493 for DFT results, see Fig. 4. Thus, our fifth-order elastic energy developed here successfully reproduces the lattice instability found in DFT over a range 0.5(σ1 + σ2) ⊂ [−73.8; 16]. The strong effect of the nonhydrostatic stresses on the lattice instability is evident: the transformation pressure under hydrostatic loading is ~75 GPa and transformation stress σ3 under uniaxial loading is ~11 GPa (or mean stress of 3.7 GPa).

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.

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, c444 = c44444 = 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 parameter47,48. In a molecular dynamics simulation15 with a Stillinger–Weber potential49, the instability for simple shear along \(<\bar{11}2>\) in the \(\left(111\right)\) plane and along the \(<111>\) in the \(\left(1\bar{1}0\right)\) plane lead to amorphization.

a For single, double, and triple simple shear strains (η1 = η2 = η3 = 0); b combination of normal and shear strains (all non-mentioned strains are zero). Fifth-order energy describes DFT results well, including shear instabilities.

Note that double and triple shears along the \(<100>\) in the \(\left(001\right)\) plane in Fig. 5a represent single shear in \(<110>\) in the \(\left(001\right)\) plane with \({\eta }_{4}^{\prime}=\sqrt{2}{\eta }_{4}\) and triaxial normal-strain loading in \(<111>\) and in the \(\left(111\right)\) plane with \({\eta }_{1}^{\prime}=2{\eta }_{4}\) and \({\eta }_{2}^{\prime}={\eta }_{3}^{\prime}=-{\eta }_{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 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 experimentally50. 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.

DFT simulations can be applied to the systems with up to 104 atoms with size up to 10 nm. Classical interatomic potentials used in molecular dynamics simulations treat the systems with up to 108 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 phenomena52,53) and without upper bound. Nonlinear elasticity based on the second- and third-degree elastic constants is determined even for monolayer graphene54,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 atomistic-continuum 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.0324,25,26,27,28. Under extreme loading, e.g., in shock waves22,23 and under high static pressure21, 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.1556,57. The fourth-order elastic energy was used to model behavior of diamond anvils21 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 samples22,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 parameters43,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 parameters43,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 currently47,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 fields61 are entirely related to nonlinear elasticity and continuum simulations. Indeed, while for small samples that can be treated by atomistic simulations, results are usually size-dependent.

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 decompression32. 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 fifth-degree) 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 orientational-dependent 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 pressures21. This approach is general and will significantly improve phase-field models for phase transformations, in contrast to the second-order elasticity used currently47,48,61. It also provides a basis for the description of the competition between different instabilities at different loadings.

We used DFT as implemented in VASP38,39,40 with the projector augmented waves (PAW) basis63,64 and PBE exchange-correlation functional65. The PAW-PBE pseudo-potential of Si had 4 valence electrons (s2p2) 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 mesh66 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 method67.

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 Fij (for more than 104 different combinations of Fij and corresponding Eij) 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 code68. 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.

DFT codes, including VASP, typically return energy with a fairly high precision, but have larger errors in stress components calculated using the Hellmann–Feynman theorem69 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 (Ecut) 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 104 entries and was processed in the format, outlined in Table 3 in ref. 70.

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.

The unit cells in the DFT simulations contained 4 atoms and were oriented along \({A}^{* }=\left\langle 110\right\rangle\), \({B}^{* }=\left\langle 1\bar{1}0\right\rangle\), and \({C}^{* }=\left\langle 001\right\rangle\). The transformation matrix for transforming the current simulation coordinate system to the natural 8-atom cubic cell, oriented along \(A=\left\langle 100\right\rangle\), \(B=\left\langle 010\right\rangle\), and \(C=\left\langle 001\right\rangle\), is:

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 sets for all figures and all raw data from DFT simulations are available at https://doi.org/10.25380/iastate.12668843.

Hill, R. & Milstein, F. Principles of stability analysis of ideal crystals. Phys. Rev. B 15, 3087 (1977).

CAS Google Scholar

Grimvall, G., Magyari-Köpe, B., Ozoliņš, V. & Persson, K. A. Lattice instabilities in metallic elements. Rev. Mod. Phys. 84, 945 (2012).

CAS Google Scholar

De Jong, M., Winter, I., Chrzan, D. & Asta, M. Ideal strength and ductility in metals from second-and third-order elastic constants. Phys. Rev. B 96, 014105 (2017).

Google Scholar

Pokluda, J., Černy`, M., Šob, M. & Umeno, Y. Ab initio calculations of mechanical properties: methods and applications. Prog. Mater. Sci. 73, 127–158 (2015).

CAS Google Scholar

Telyatnik, R., Osipov, A. & Kukushkin, S. Ab initio modelling of nonlinear elastoplastic properties of diamond-like C, SiC, Si, Ge crystals upon large strains. Mater. Phys. Mech. 29, 1–16 (2016).

CAS Google Scholar

Wang, J., Yip, S., Phillpot, S. & Wolf, D. Crystal instabilities at finite strain. Phys. Rev. Lett. 71, 4182 (1993).

CAS Google Scholar

Mizushima, K., Yip, S. & Kaxiras, E. Ideal crystal stability and pressure-induced phase transition in silicon. Phys. Rev. B 50, 14952 (1994).

CAS Google Scholar

Levitas, V. I., Chen, H. & Xiong, L. Lattice instability during phase transformations under multiaxial stress: modified transformation work criterion. Phys. Rev. B 96, 054118 (2017).

Google Scholar

Levitas, V. I., Chen, H. & Xiong, L. Triaxial-stress-induced homogeneous hysteresis-free first-order phase transformations with stable intermediate phases. Phys. Rev. Lett. 118, 025701 (2017).

Google Scholar

Zarkevich, N. A., Chen, H., Levitas, V. I. & Johnson, D. D. Lattice instability during solid-solid structural transformations under a general applied stress tensor: Example of Si I to Si II with metallization. Phy. Rev. Lett. 121, 165701 (2018).

CAS Google Scholar

Binggeli, N. & Chelikowsky, J. R. Elastic instability in α-quartz under pressure. Phys. Rev. Lett. 69, 2220 (1992).

CAS Google Scholar

Kingma, K. J., Meade, C., Hemley, R. J., Mao, H.-K. & Veblen, D. R. Microstructural observations of α-quartz amorphization. Science 259, 666–669 (1993).

CAS Google Scholar

Brazhkin, V. & Lyapin, A. Lattice instability approach to the problem of high-pressure solid-state amorphization. High Press. Res. 15, 9–30 (1996).

Google Scholar

Zhao, S. et al. Shock-induced amorphization in silicon carbide. Acta Mater. 158, 206–213 (2018).

CAS Google Scholar

Chen, H., Levitas, V. I. & Xiong, L. Amorphization induced by 60° shuffle dislocation pileup against different grain boundaries in silicon bicrystal under shear. Acta Mater. 179, 287–295 (2019).

CAS Google Scholar

Tallon, J. L. A hierarchy of catastrophes as a succession of stability limits for the crystalline state. Nature 342, 658 (1989).

CAS Google Scholar

Levitas, V. I. & Ravelo, R. Virtual melting as a new mechanism of stress relaxation under high strain rate loading. Proc. Natl. Acad. Sci. USA 109, 13204–13207 (2012).

CAS Google Scholar

Umeno, Y. & Černy`, M. Effect of normal stress on the ideal shear strength in covalent crystals. Phys. Rev. B 77, 100101 (2008).

Google Scholar

Černý, M., Řehák, P., Umeno, Y. & Pokluda, J. Stability and strength of covalent crystals under uniaxial and triaxial loading from first principles. J. Phys. Condens. Matter 25, 035401 (2012).

Google Scholar

Tang, M. & Yip, S. Lattice instability in β-SiC and simulation of brittle fracture. J. Appl. Phys. 76, 2719–2725 (1994).

CAS Google Scholar

Levitas, V. I., Kamrani, M. & Feng, B. Tensorial stress-strain fields and large elastoplasticity as well as friction in diamond anvil cell up to 400 GPa. Npj Comput. Mater. 5, 1–11 (2019).

CAS Google Scholar

Clayton, J. D. Analysis of shock compression of strong single crystals with logarithmic thermoelastic-plastic theory. Int. J. Eng. Sci. 79, 1–20 (2014).

Google Scholar

Clayton, J. D. Crystal thermoelasticity at extreme loading rates and pressures: analysis of higher-order energy potentials. Extreme Mech. Lett. 3, 113–122 (2015).

Google Scholar

Zhao, J., Winey, J. M. & Gupta, Y. M. First-principles calculations of second- and third-order elastic constants for single crystals of arbitrary symmetry. Phys. Rev. B 75, 094105 (2007).

Google Scholar

Łopuszyński, M. & Majewski, J. A. Ab initio calculations of third-order elastic constants and related properties for selected semiconductors. Phys. Rev. B 76, 045202 (2007).

Google Scholar

Cao, T., Cuffari, D. & Bongiorno, A. First-principles calculation of third-order elastic constants via numerical differentiation of the second Piola-Kirchhoff stress tensor. Phys. Rev. Lett. 121, 216001 (2018).

CAS Google Scholar

Wang, H. & Li, M. Ab initio calculations of second-, third-, and fourth-order elastic constants for single crystals. Phys. Rev. B 79, 224102 (2009).

Google Scholar

Telichko, A. V. et al. Diamond’s third-order elastic constants: Ab initio calculations and experimental investigation. J. Mater. Sci. 52, 3447–3456 (2017).

CAS Google Scholar

Guo, D. & An, Q. Transgranular amorphous shear band formation in polycrystalline boron carbide. Int. J. Plast. 121, 218–226 (2019).

CAS Google Scholar

An, Q., Goddard III, W. A. & Cheng, T. Atomistic explanation of shear-induced amorphous band formation in boron carbide. Phys. Rev. Lett. 113, 095501 (2014).

Google Scholar

Gao, Y. et al. Shear driven formation of nano-diamonds at sub-gigapascals and 300 k. Carbon 146, 364–368 (2019).

Google Scholar

Patten, J., Cherukuri, H. & Yan, J. In High Pressure Surface Science and Engineering (eds Gogotsi, Y. & Domnich, V.) Section 6, 521–540 (Institute of Physics, Bristol, 2004).

Gogotsi, Y. & Domnich, V. In High Pressure Surface Science and Engineering (eds Gogotsi, Y. & Domnich, V.) Section 5, 381–442 (Institute of Physics, Bristol, 2004).

Blank, V. D. & Estrin, E. I. Phase Transitions in Solids Under High Pressure (CRC Press, 2013).

Vanhellemont, J., Swarnakar, A. K. & der Biest, O. V. Temperature dependent Young’s modulus of Si and Ge. ECS Trans. 64, 283–292 (2014).

CAS Google Scholar

Schall, J. D., Gao, G. & Harrison, J. A. Elastic constants of silicon materials calculated as a function of temperature using a parametrization of the second-generation reactive empirical bond-order potential. Phys. Rev. B 77, 115209 (2008).

Google Scholar

Kulda, J., Strauch, D., Pavone, P. & Ishii, Y. Inelastic-neutron-scattering study of phonon eigenvectors and frequencies in Si. Phys. Rev. B 50, 13347–13354 (1994).

CAS Google Scholar

Kresse, G. & Hafner, J. Ab initio molecular dynamics for liquid metals. Phys. Rev. B 47, 558 (1993).

CAS Google Scholar

Kresse, G. & Hafner, J. Ab initio molecular-dynamics simulation of the liquid-metal–amorphous-semiconductor transition in germanium. Phys. Rev. B 49, 14251 (1994).

CAS Google Scholar

Kresse, G. & Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci. 6, 15–50 (1996).

CAS Google Scholar

Hall, J. J. Electronic effects in the elastic constants of n-type silicon. Phys. Rev. 161, 756 (1967).

CAS Google Scholar

McSkimin, H. & Andreatch Jr, P. Measurement of third-order moduli of silicon and germanium. J. App. Phys. 35, 3312–3319 (1964).

CAS Google Scholar

Levitas, V. I. Phase-field theory for martensitic phase transformations at large strains. Inter. J. Plast. 49, 85–118 (2013).

CAS Google Scholar

Teodosiu, C. Elastic models of crystal defects (Springer Science & Business Media, 2013).

Mosyagin, I. et al. Ab initio calculations of pressure-dependence of high-order elastic constants using finite deformations approach. Comput. Phys. Commun. 220, 20–30 (2017).

CAS Google Scholar

Tiesinga, E., Mohr, P. J., Newell, D. B. & Taylor, B. N. The 2018 CODATA Recommended Values of the Fundamental Physical Constants, NIST database developed by Baker, J., Douma, M. & Kotochigova, S. (National Institute of Standards and Technology, Gaithersburg, MD 20899, USA, 2020).

Levitas, V. I. Phase field approach for stress-and temperature-induced phase transformations that satisfies lattice instability conditions. Part I. General theory. Inter. J. Plast. 106, 164–185 (2018).

CAS Google Scholar

Babaei, H. & Levitas, V. I. Phase-field approach for stress-and temperature-induced phase transformations that satisfies lattice instability conditions. Part 2. simulations of phase transformations Si I - Si II. Inter. J. Plast. 107, 223–245 (2018).

CAS Google Scholar

Stillinger, F. H. & Weber, T. A. Computer simulation of local order in condensed phases of silicon. Phys. Rev. B 31, 5262 (1985).

CAS Google Scholar

Deb, S. K., Wilding, M., Somayazulu, M. & McMillan, P. F. Pressure-induced amorphization and an amorphous–amorphous transition in densified porous silicon. Nature 414, 528 (2001).

CAS Google Scholar

Liu, C., Song, X., Li, Q., Ma, Y. & Chen, C. Smooth flow in diamond: atomistic ductility and electronic conductivity. Phys. Rev. Lett. 123, 195504 (2019).

CAS Google Scholar

Levitas, V. I. & Samani, K. Size and mechanics effects in surface-induced melting of nanoparticles. Nat. Commun. 2, 1–6 (2011).

Google Scholar

Levitas, V. I. & Javanbakht, M. Surface-induced phase transformations: multiple scale and mechanics effects and morphological transitions. Phys. Rev. Lett. 107, 175701 (2011).

Google Scholar

Lee, C., Wei, X., Kysar, J. W. & Hone, J. Measurement of the elastic properties and intrinsic strength of monolayer graphene. Science 321, 385–388 (2008).

CAS Google Scholar

Cadelano, E., Palla, P. L., Giordano, S. & Colombo, L. Nonlinear elasticity of monolayer graphene. Phys. Rev. Lett. 102, 235502 (2009).

Google Scholar

Hemley, R. J. et al. X-ray imaging of stress and strain of diamond, iron, and tungsten at megabar pressures. Science 276, 1242–1245 (1997).

CAS Google Scholar

Feng, B., Levitas, V. I. & Hemley, R. J. Large elastoplasticity under static megabar pressures: Formulation and application to compression of samples in diamond anvil cells. Int. J. Plast. 84, 33–57 (2016).

Google Scholar

Levitas, V. I. & Javanbakht, M. Thermodynamically consistent phase field approach to dislocation evolution at small and large strains. J. Mech. Phys. Solids 82, 345–366 (2015).

Google Scholar

Javanbakht, M. & Levitas, V. I. Phase field approach to dislocation evolution at large strains: Computational aspects. Int. J. Solids Struct. 82, 95–110 (2016).

Google Scholar

Levitas, V. I., Jafarzadeh, H., Farrahi, G. H. & Javanbakht, M. Thermodynamically consistent and scale-dependent phase field approach for crack propagation allowing for surface stresses. Int. J. Plast. 111, 1–35 (2018).

Google Scholar

Babaei, H. & Levitas, V. I. Stress-measure dependence of phase transformation criterion under finite strains: Hierarchy of crystal lattice instabilities for homogeneous and heterogeneous transformations. Phys. Rev. Lett. 124, 075701 (2020).

CAS Google Scholar

Ji, C. et al. Shear-induced phase transition of nanocrystalline hexagonal boron nitride to wurtzitic structure at room temperature and lower pressure. Proc. Natl Acad Sci USA 109, 19108–19112 (2012).

CAS Google Scholar

Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B 50, 17953 (1994).

Google Scholar

Kresse, G. & Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B 59, 1758 (1999).

CAS Google Scholar

Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 77, 3865 (1996).

CAS Google Scholar

Monkhorst, H. J. & Pack, J. D. Special points for Brillouin-zone integrations. Phys. Rev. B 13, 5188 (1976).

Google Scholar

Johnson, D. D. Modified Broyden’s method for accelerating convergence in self-consistent calculations. Phys. Rev. B 38, 12807 (1988).

CAS Google Scholar

Zarkevich, N. A. & Johnson, D. D. Nudged-elastic band method with two climbing images: finding transition states in complex energy landscapes. J. Chem. Phys. 142, 024106 (2015).

Google Scholar

Güttinger, P. Das verhalten von atomen im magnetischen drehfeld. Z. Physik 73, 169–184 (1932).

Google Scholar

Zarkevich, N. A. Structural database for reducing cost in materials design and complexity of multiscale computations. Complexity 11, 36–42 (2006).

Google Scholar

Download references

V.I.L. and H.C. are supported by NSF (CMMI-1943710 & MMN-1904830), ONR (N00014-16-1-2079), & XSEDE (MSS170015). N.A.Z. and D.D.J. are supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences, Materials Science & Engineering Division. Ames Laboratory is operated for DOE by Iowa State University under contract DE-AC02-07CH11358. X.C.Z. and H.C. are also sponsored by the National Key Research and Development Program of China (2018YFC1902404), the National Natural Science Foundation of China (51725503, 51975211), and Innovation Program of Shanghai Municipal Education Commission (2019-01-07-00-02-E00068).

Key Laboratory of Pressure Systems and Safety, Ministry of Education, School of Mechanical and Power Engineering, East China University of Science and Technology, Shanghai, 200237, China

Hao Chen & Xiancheng Zhang

Ames Laboratory, U.S. Department of Energy, Iowa State University, Ames, IA, 50011-3020, USA

Nikolai A. Zarkevich, Valery I. Levitas & Duane D. Johnson

Department of Aerospace Engineering, Iowa State University, Ames, IA, 50011, USA

Valery I. Levitas

Department of Mechanical Engineering, Iowa State University, Ames, IA, 50011, USA

Valery I. Levitas

Department of Materials Science & Engineering, Iowa State University, Ames, IA, 50011, USA

Duane D. Johnson

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

H.C. developed the fifth order elastic potential and performed all post processing of DFT simulations to obtain all elastic constants, stress–strain curves, and instability conditions. N.A.Z. performed all DFT simulations. V.I.L. designed research and supervised H.C. work. D.D.J. supervised N.A.Z. work. All authors participated in analysis and writing and editing the paper.

Correspondence to Hao Chen, Valery I. Levitas or Duane D. Johnson.

The authors declare no competing interests.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and Permissions

Chen, H., Zarkevich, N.A., Levitas, V.I. et al. Fifth-degree elastic energy for predictive continuum stress–strain relations and elastic instabilities under large strain and complex loading in silicon. npj Comput Mater 6, 115 (2020). https://doi.org/10.1038/s41524-020-00382-8

Download citation

Received: 05 April 2020

Accepted: 02 July 2020

Published: 04 August 2020

DOI: https://doi.org/10.1038/s41524-020-00382-8

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative

Communications Materials (2022)

Nonlinear Dynamics (2022)