Ab initio calculations of pressure-dependence of high-order elastic constants using finite deformations approach

Accepted Manuscript Ab initio calculations of pressure-dependence of high-order elastic constants using finite deformations approach I. Mosyagin, A.V. Lugovskoy, O.M. Krasilnikov, Yu. Kh. Vekilov, S.I. Simak, I.A. Abrikosov

PII: DOI: Reference:

S0010-4655(17)30190-X http://dx.doi.org/10.1016/j.cpc.2017.06.008 COMPHY 6242

To appear in:

Computer Physics Communications

Received date : 27 December 2016 Revised date : 26 May 2017 Accepted date : 7 June 2017 Please cite this article as: I. Mosyagin, A.V. Lugovskoy, O.M. Krasilnikov, Y.K. Vekilov, S.I. Simak, I.A. Abrikosov, Ab initio calculations of pressure-dependence of high-order elastic constants using finite deformations approach, Computer Physics Communications (2017), http://dx.doi.org/10.1016/j.cpc.2017.06.008 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Revised manuscript w/o linenumbers

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Ab initio calculations of pressure-dependence of high-order elastic constants using finite deformations approach. Mosyagin I. Yu.a,b , Lugovskoy A. V.b , Krasilnikov O. M.c , Vekilov Yu. Kh.b,c , Simak S. I.a , Abrikosov I. A.a a Department

of Physics (IFM), Link¨ oping University, SE-58183 Link¨ oping, Sweden Modelling and Development Laboratoy, NUST “MISIS”, RU-119991 Moscow, Russia c Department of Theoretical Physics and Quantum Technologies, NUST “MISIS”, RU-119991 Moscow, Russia

b Materials

Abstract We present a description of a technique for ab initio calculations of the pressure dependence of second- and third-order elastic constants. The technique is based on an evaluation of the corresponding Lagrangian stress tensor derivative of the total energy assuming finite size of the deformations. Important details and parameters of the calculations are highlighted. Considering body-centered cubic Mo as a model system, we demonstrate that the technique is highly customizable and can be used to investigate non-linear elastic properties under high-pressure conditions. Keywords: ab initio calculations, elastic moduli, pressure effects in solids and liquids PACS: 31.15.A-, 62.20.de, 62.50.-p 1. Introduction There is a notable increasing interest in higher order elastic constants (HOEC) of solids [1, 2, 3, 4, 5, 6]. The study of nonlinear elasticity helps to reveal the details of complex behavior of materials. For example, investigating the trends of HOECs one can understand mechanisms of structural instabilities [3] or incorporate proper description of nonlinear elasticity in case where there is a discrepancy between more general theory and experiment [4]. Usually the elasticity theory is considered in an approximation of infinitesimal deformations [7]. Such approach is the most logical choice, when the applied deformation is small compared to inter-atomic distances of undeformed material. But inter-atomic potentials in real materials are anharmonic. The explicit accounting of anharmonicity in study of solids becomes more and more prominent. This is quite evident when theoretical investigation of materials is performed for realistic conditions including the extreme ones, which are, for instance, of interest for cutting edge studies [8, 9]. Standard theoretical methods Email address: [email protected] (Mosyagin I. Yu.)

Preprint submitted to Elsevier

May 26, 2017

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

of defining components of elastic constants tensor typically operate in a linear limit of elasticity theory. Although there are comprehensive works present in the literature, which are dedicated to finite deformation and calculation of higher order elastic constants [10, 11] the techniques, applied in these works do not take into account the external conditions, e.g. pressure and temperature. We also would like to note, that there is a lack of systemic study of numerical aspects of higher order elastic constants calculation. In these work we expand and clarify most important numeric aspects of methodology of non-linear elastic properties investigation for solids under pressure, using finite strain theory to calculate second- and third- order elastic constants under pressure. Doing so we concentrate on the algorithm and the strategy for finding optimal calculation parameters, both regular convergence tests in regards to simulation parameters, as well as internal algorithm parameters, such as deformation range. The importance of careful adjustment of the parameters naturally comes from the fact, that even the second order elastic constants require highly converged mesh of the k-points and plane wave cut-off energy [12]. Obviously the importance of this aspect increases, if the HOEC are calculated. Although discussion of convergence of third and even fourth order elastic constants calculation parameters is present in [10, 11], we additionaly focus on the impact of external loading and its effect on the calculation procedure. To demonstrate the application of this technique we choose molybdenum, a 6th group transitional metal with a body-centered cubic structure. We investigate the second and third order elastic constants in the pressure range of up to 300 GPa. The investigation of elastic properties at pressures this high is motivated by the recent breakthrough in experimental techniques of achieving static pressures [13, 14, 15]. The measurement of elastic constants at extreme static compression is a challenging task that contains some uncertainties [16]. This makes theoretically obtained elastic constants trends under pressure a valuable addition to the knowledge about studied materials. We would also like to mention, that transition metals, including molybdenum, are often used in high pressure experiments [17, 15, 13, 14, 18]. 2. Method We start with an outline of basic ideas of elasticity theory behind the proposed method. In this work we focus on calculation-specific technical issues and then give a detailed outline of the computational algorithm. A short description of a convetional infinitesimal strains method can also be found in this section for completeness. A detailed theoretical description, including the full set of deformation schemes, can be found in [19, 20, 21, 22]. 2.1. Basic concepts of elasticity theory behind the method 2.1.1. Elastic constants tensor, infinitesimal strain approach In this subsection we outline the basic methods of elastic constants calculation, with an emphasis on differences between a standard method that uses infinitesimal distortions and our method. Consider a crystal with elementary cell of volume V , that is subjected to deformation described by strain tensor

2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

eij . We can write down the following expression for the energy of the crystal under deformation applied [23]: 1 X U (e) = U (0) + V0 cijkl eij ekl + . . . , 2

(1)

ijkl

where U (0) is the energy of undeformed state, and i, j, k, l are coordinate indices with values range 1–3. Usually infinitesimal strain tensor is applied which assumes the applied deformation to be small, compared to inter-atomic distances of undeformed material. If such assumption is made, higher order terms in (1) are considered negligible. cijkl are usually expressed using Voigt’s notation, which maps them to the form of cij , where the new i and j are integers in the range from 1 to 6: 11 → 1, 22 → 2, 33 → 3, 23 → 4, 13 → 5, 12 → 6. Total amount of non-zero elements in cijkl tensor depends on the system’s symmetry, but is limited to maximum 21 for triclinic lattices [7]. In the case of tetragonal lattices, Bravais lattice symmetry allows one to drop that number down to 6, while hexagonal lattices have 5 independent non-zero elements and there’s only 3 of them for cubic crystals. To obtain the desired number of elastic constants one should take as many different deformation tensors. This paper focuses on the case of cubic crystals, for which the only non-zero constants are c11 , c12 and c44 (all other constants being equal to one of them or zero), and bulk modulus is defined as B = 31 (c11 + 2c12 ). In practice, when studying cubic crystals, it is common to introduce additional elastic constant c0 = (c11 − c12 )/2. This is one of the two so-called shear moduli, the second one is c44 . The other common combination of these elastic constants is Cauchy pressure c00 = (c12 − c44 ), which is sometimes used to estimate ductility [24]. Depending on application, there might be other additional combinations when investigating cubic system. Nevertheless, all of them are based on a combination of secondorder elastic constants c11 , c12 and c44 . If deformation is not volume-conserving, than the corresponding volume change, in general, produces significantly higher energy changes than the deformation corresponding to the change in shape. In order to separate two effects, one usually considers volume-conserving deformations. In that case, energy change with deformation can be expressed using (1) in the form of ∆E = A · V · c · δ 2 , where A is some numeric constant, depending on the stress tensor σij applied, c is a corresponding elastic constant (or their combination) and δ is a parameter of deformation (see Ref. [23] for the explicit relation between δ and ε in Eq. (1)). The algorithm in this case is the following: for the volume of interest, calculate the deformed state energies E(δ) for a number of δ values (most often, these 6 points are recommended: 0.00, 0.01, 0.02, 0.03, 0.04, and 0.05 [23]). The E(δ) is assumed to be quadratic and thus E(δ 2 ) is approximated with a linear function, the slope of E(δ 2 ) provides the second-order elastic constant corresponding to this deformation. In this paper, the approach, based on infinitesimal deformations was used to calculate second-order elastic constant c44 and shear modulus c0 , in order to provide the benchmark for comparison with the method where the infinitesimal deformations assumption is lifted and the second, as well as higher order elastic constants are calculated. To calculate c0 , the following orthorhombic deformation was used: 3

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65



1+δ Do + I =  0 0

0 1−δ 0

 0 0 .

1 1−δ 2

(2)

By substituting the corresponding deformation tensor components in (1), we obtain the following expression for deformation energy: ∆E(δ) = 2V c0 δ 2 + O(δ 4 ),

(3)

To calculate c44 , the following monoclinic deformation was used:   1 δ 0 0  Dm + I = δ 1 1 0 0 1−δ 2

(4)

Similarly, using 4 in (1), we obtain:

∆E(δ) = 2V c44 δ 2 + O(δ 4 )

(5)

Note, that on the right hand side of both (2) and (4) the third diagonal 1 element is 1−δ 2 . This is to ensure that the determinants of both tensors is equal to 1 and thus such deformation preserves the cell’s volume. In this work, we used distortions (2) and (4) with corresponding energy expansions (3) and (5) to calculate second-order elastic constants c0 and c44 of bcc molybdenum at the set of volumes V in the range of V /V0 = [0.65 : 1.00], where V0 is the equilibrium volume corresponding to 0 GPa, while 0.65V0 corresponds to roughly 300 GPa (see Fig. 1). Both distortions were calculated using four values of the deformation parameter δ: {0.00, 0.01, 0.02, 0.04}. We would like to note that expressions (2) and (4) are symmetrical with respect to the sign of δ. In practice, the negative part is completely optional, but can be used for a verification. In most cases it’s enough to carry out calculations at positive values of δ. 2.1.2. Finite deformations approach The expansion on the right-hand side of Eq. (1) is intentionally limited to second order. To correctly take into account higher order expansion terms of the energy of deformed crystal, we need to use finite strain tensor [25] 1X ηij = (αik αkj − δij ), (6) 2 k

where αij = ∂ri /∂Rj is a tensor of transformation coefficients of the system under investigation; δij is a Kronecker delta and R = (Rx , Ry , Rz )T and r = (rx , ry , rz )T being spacial configurations of, respectively, undeformed and deformed state. According to the nonlinear theory of elasticity [25], the internal energy U and a free energy F of deformed crystal can be expanded in a power series of η: U (R, ηij , S) = U (R, S) + V

1 X S 1 X S Cijkl ηij ηkl + V Cijklmn ηij ηkl ηmn 2! 3! ijkl

1 +V 4!

ijklmn

X

ijklmnpq

S Cijklmnpq ηij ηkl ηmn ηpq + · · · , and

(7) 4

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

F (R, ηij , T ) = F (R, T ) + V

1 X T 1 X T Cijkl ηij ηkl + V Cijklmn ηij ηkl ηmn 2! 3! ijkl

1 +V 4!

ijklmn

X

ijklmnpq

T Cijklmnpq ηij ηkl ηmn ηpq + · · · ,

(8)

T S where Cijkl... and Cijkl... represent isothermal and adiabatic elastic constants, and we temporary switch back from Voigt’s notation. Using these expressions, one can introduce the thermodynamic definition of elastic constants of n−th order (n > 2):   ∂nF 1 T Cijkl... = V0 ∂ηij ∂ηkl . . . η=0   (9) 1 ∂nU S Cijkl... = V0 ∂ηij ∂ηkl . . . η=0

Here we once again highlight that the energy derivatives are taken with respect to the undeformed configuration corresponding to ηij = 0, in line with Eq. (6). We would like to also note that it is possible to transform (7) and (8) back to infinitesimal formalism, if one takes ηij = eij + 12 e2ij and suppresses all terms higher than the second order of eij .

2.1.3. Pressure effects Trends of elastic constants under pressure and temperature are always of interest, since observing this dependencies allows to evaluate mechanical properties, e.g. brittleness/ductility of material, as well as provide valuable information about the mechanical stability of material at given conditions. For the case of non-zero temperature T 6= 0 and applied pressure P 6= 0, the Gibbs free energy G = F + P V and the enthalpy H = U + P V are used for crystal description. Using (9), we can now obtain the expression for the case of external load:  ∂nG ∂ηij ∂ηkl . . . η=0   1 ∂nH = V0 ∂ηij ∂ηkl . . . η=0

1 T eijkl... C = V0 S eijkl... C



(10)

Again, as in Eq. (9), the derivatives are taken with reference state corresponding to ηij = 0. Elastic constants defined by (10) take into account the work that system does against applied hydrostatic pressure. When there is no external load (P = 0), we revert back to (9). For the case of T = 0 K, the free energy F = U − T S equals T S to internal energy U , and, therefore we can define Cijkl... = Cijkl... = Cijkl... . In order to simplify notation, we drop the super-index and in order to take into account that we limit the calculation in this paper to zero temperature, we will generally use H instead of G as a thermodynamic potential. 5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

The change in enthalpy of the deformed crystal ∆H = H(η) − H(η = 0), therefore, can be expressed as: ∆H 1 Xe 1 X e = Cijkl ηij ηkl + Cijklmn ηij ηkl ηmn V 2! 3! ijkl

ijklmn

1 + 4!

X

ijklmnpq

eijklmnpq ηij ηkl ηmn ηpq + · · · C

(11)

Denoting ∆H = H(P, η) − H(P, 0), ∆U = U (P, η) − U (P, 0), and ∆V = V − V0 , where V0 is an equilibrium volume for given P and T , we have: ∆H ∆U ∆V ∆U = +P = + P (J − 1), V0 V0 V0 V0

(12)

where J = det kαij k is a Jacobian of transformation, that happens in the system under applied deformation ηij . Note that, since the deformation schemes applied are not conserving the volume of the crystal, in actual calculations we always obtain ∆U and hence Cijkl . eijkl provide the actual practical value since they can be comElastic constants C pared with experimentally obtained elastic constants. In many cases, secondeijkl and second-order elastic constants calculated using infinitesimal order C strain tensor (1) can be also compared directly, however when non-linear elastic effects become significant at given conditions, they can differ. We show later that in case of body-centered cubic molybdenum this discrepancy in secondorder elastic constants is indeed observable at higher pressures. eijkl and Cijkl can easily be obtained from the subThe relations between C stitution of (7) into (12), since corresponding form of αij is known: these relations are provided in AppendixC. From now on we denote the elastic constants as C˜ijkl , since this provides unified description of elastic moduli obtained by all covered methods. 2.1.4. Defining deformations for finite strain tensor Using (6) with (12), (11), (8), and thus expressing αij as a function of ηij we can explicitly take into account displacement terms up to desired order (in this case the terms up to the third order are taken into account): 1 1 αij = δij + ηij − ηki ηkj + ηrk ηri ηkj + o(η 4 ) 2 2

(13)

If all components of ηij are expressed using one parameter η, for deformed state we obtain ∆U (ηij ) = ν1 η + ν2 η 2 + ν3 η 3 + o(η 4 ) (14) where ν1 is either null or proportional to the pressure in the system, and coefficients νi with i > 2 would contain linear combinations of the elastic constants of corresponding order. For the case of cubic crystals, we have three 2nd -order EC, six 3rd -order EC and eleven 4th -order EC. Explicit energy expansions for up to fourth order elastic constants are provided in AppendixA, while deformations for up to fourth order in η can be found in AppendixB. If one is interested in obtaining all 4th-order elastic constants for cubic crystal, we need to have 11 distinct deformations applied to our system to solve the resulting system of 6

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

equations, while for a full set of 3rd-order elastic constant 6 distinct deformations are sufficient. Deformation schemes, necessary to obtain all second and third order elastic constants are given in Table 1. Additional deformations for fourth order elastic constants can be found in Table A.2 of AppendixA, while expressions for fourth order terms ν4 can be found in Table A.3 of the same appendix. A few extra schemes were taken to verify the consistency of the technique applied. The methodology description and parameters convergence in the article is performed using these exact nine schemes. Deformation scheme N (η1 , η2 , η3 , η4 , η5 , η6 ) S1 (η, 0, 0, 0, 0, 0) S2 (η, η, 0, 0, 0, 0) S3 (η, −η, 0, 0, 0, 0) (η, η, η, 0, 0, 0) S4 S5 (0, 0, 0, η, 0, 0) (η, 0, 0, η, 0, 0) S6 S7 (η, 0, 0, 0, 0, η) S8 (0, 0, 0, η, η, η) S9 (η, 0, 0, η, η, η)

Deformation energy ν1 η + ν2 η 2 + ν3 η 3 + o(η 4 ) ν1 ν2 ν3 1 e 1 e −P 2 C11 6 C111 1 e e11 + C e12 e112 −2P C C +C 111 3 e11 − C e12 0 C 0 3 e 1 e e e e −3P 2 C11 + 3C12 2 C111 + 3C112 + C123 e 0 2C44 0 1 e 1 e e e155 −P C + 2 C C + 2C 44 2 11 6 111 1 e 1 e e e −P 2 C11 + 2C44 6 C111 + 2C144 e e 0 6C44 8C456 1 e e44 1 C e111 + 2C e144 + 4C e155 + 8C e456 −P C11 + 6C 2

6

Table 1: Example of possible deformations and corresponding energy expansions in the series of η.

2.1.5. Calculation details To calculate elastic constants we need to obtain the energy of deformed and undeformed states of material. For this purpose we used density functional theory as implemented in the VASP code [26, 27, 28, 29]. The implementation is, of course not limited to this exact code, however, our technical discussion is heavily focused on this particular package. In this work the exchangecorrelation part was treated with Generalized gradient approximation (GGA) with Perdew-Burke-Ernzerhof parametrization [30]. Projector augmented waves method (PAW), as implemented in VASP [31], was used. In general, the amount of electrons in valence band is also a convergence parameter, and should be increased if one wants to perform simulations under high pressure. In case of bcc Mo, we used semi-core PAW-potential which includes 12 electrons in the valence band. That number was sufficient in our case. The plane wave cut-off energy and k-points mesh were obtained from careful convergence test, described below. The final set of this parameters for molybdenum is 700 eV cut-off and 30 × 30 × 30 k-points mesh. 2.1.6. Verifying the description of compressibility: equation of state Checking the correctness of equation of state parameters is an essential step to correctly model the effects of pressure on the material. In our case the chosen PAW-potential, as well as all computational parameters should provide realistic description of material compressibility, which is required for the relevance of high pressure elastic properties predictions we make. We use Birch–Murnaghan equation of state, which gives: 7

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

3B0 P = 2

"

V0 V

 73

−



V0 V

 35 #

"

3 × 1 + (B00 − 4) 4

"

V0 V

 23

−1

##

,

(15)

where V0 is the equilibrium volume of the system, B0 is its bulk modulus, and B00 is it’s derivative. Bulk modulus obtained from (15) is of the most important elastic constants: it represents the relative change in pressure when its volume changes isotropically: ∂P (16) B0 = − V 0 ∂V V =V0

In this P (V ) form, Birch-Murnaghan equation of state is often used by experimentalists to fit their P (V ) and obtain B. When doing theoretical calculations, one obtains E(V ) expression. To go from P (V ) in eq. (15) to E(V ), one needs to integrate it: 9 E(V ) = E0 − B0 V0 16



V0 V

 23

!2

−1

B00

+3



V0 V

 32

+6

!

(17)

where E0 is the equilibrium energy of the system, E0 = E(V0 ). To ensure, that our calculations provide the realistic description of bcc molybdenum compressibility we provide the comparison of Birch-Murnaghan equation of state, obtained from (17), with coefficients (B0 , B00 , V0 ) then used in (15) to obtain P (V /V0 ) dependence (see Figure 1). We compare these results with the P (V /V0 ) dependence, obtained using the proposed finite strain tensor methodology and computational scheme S1 (see Table. 1). We, of course, compare our results with the set of experimentally obtained data available from the literature. Agreement between data obtained by different methods in the framework of this study is good enough to conclude, that our method provides results consistent with the standard technique (see detailed discussion in the “Results and discussion” section). 2.2. Algorithm outline The algorithm is shown graphically in Fig. 2. Overall, our algorithm can be expressed as follows: 1. Define the set of volumes that will give us the pressure in the system in the desired range This is done by calculating undeformed system energies on a grid of volumes fitted with EOS. An example of such calculation is shown in Fig. 1. 2. Depending on the desired set of elastic constants, decide upon deformation matrices to use. For example, if one’s interest consist in finding only the second-order elastic constants as a function of pressure for the system with cubic symmetry, three distinct schemes is sufficient. Possible deformations can be taken in the form expressed in Table 1, e.g. 1st, 3rd and 5th deformation scheme.

8

300

Mo EOS this work Birch-Murnaghan fit Exp. Duffy[32] Exp. Liu[33] Exp. Nguyen[34] Exp. Wang[35] Exp. Hixson[36]

250

Pressure, GPa

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

200 150 100 50 0 0.60

0.65

0.70

0.75

0.80 V /V0

0.85

0.90

0.95

1.00

Figure 1: bcc Mo equation of state. Black solid triangles denote the results obtained using the proposed technique of finite strain tensor and computational scheme S1 (Table 1). Conventional calculation of pressure using Birch-Murnaghan equation of state fit (15) is shown with red line. All other points represent available experimental data. Birch-Murnaghan equation of state fitting was carried over the same set of points that was used to get equation of state by applying our method and extracting pressure from the regression coefficients ν1 .

3. Decide on deformation η range, for example [−0.045 . . . 0.045], and its step value, e.g. 0.003. This is a matter of balance between smooth enough E(η) dependence and minimization of the total amount of calculations. The effect of η range on fit coefficients in equation (14) is discussed in subsection 2.4. If TOECs are of interest, usually the broader range of η is required to reach the desired precision. 4. For each of the volumes from step 1, we calculate E(η) dependence for each deformation scheme from step 2 with η in a range defined in step 3. For each of E(η) dependencies, we obtain the corresponding coefficients ν1,2,3 using least-square regression. 5. First-order coefficients ν1 are used to get the equation of state P (V ). Knowing those, it is possible to switch from νi (V ) to νi (P ). For example, if using scheme S1 from Table 1, one obtains set of fitting coefficients as a functions of volume: ν1 (V ), ν2 (V ) and ν3 (V ), that are used to get e11 (V ) = 2ν2 (V ) pressure P (V ) = −ν1 (V ), second order elastic constant C e and third order elastic constant C111 (V ) = 6ν3 (V ). Each scheme in Table 1 would have similar corresponding relations. Combining obtained P (V ) with regression coefficients ν2,3 (η) for multiple schemes, we get eieijkl as a function of ther direct dependencies or linear combinations of C pressure. 9

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

A set/range of volumes to be used

pressure range of interest

Calculation parameters, including: 1. η range, 2. plane-wave cut-off energy, 3. k-point mesh. Deformation schemes selection Calculate E(η) for each volume and for each deformation scheme used. Obtain corresponding regression coefficients ν1,2,3,4

eijkl Choice of C

Regression coefficients ν1,2,3,4 as a functions of volume P (V ) dependence

eijkl (V ) A set of C

Elastic constants as a funceijkl (P ) tions of pressure C

Figure 2: Algorithm for calculation of elastic constants as a function of pressure using finite strain tensor.

2.3. Accuracy of calculations using finite strain tensors 2.3.1. Convergence in terms of cut-off energy and k-points mesh density It is important not only to converge in terms of cut-off energy and k-points mesh, but also to take into account that at higher pressures those parameters need to be increased. For the undeformed cell, the energy convergence is typically achieved at lower energy cut-off, than for the deformed case, even at high pressureses. When carrying out calculations with undeformed cell it is safe to use the default setting for cut-off energy (i.e. 225 eV for the bcc Mo PAW-potential used). For the case of pressures that are comparable in their magnitude with the bulk modulus of the material under investigation, it is important to increase the calculation precision. In particular, to avoid potential frozen core overlapping, it might be necessary to use higher cut-off energy. Effectively this lowers the real space cut-off radius and allows for ions to come closer to each other without producing computational artifacts. e11 and Figure 3 demonstrates the convergence of two elastic constants C e111 , both calculated using deformation S1 from Table 1 with increasing size C e11 at two of automatically generated k-points mesh. The left frames show C pressures: 0 GPa and 300 GPa. The points are connected as a guide for an eye. All values calculated with default cut-off energy for the PAW potential used. 10

e111 at 0 and 300 GPa. While The right frames show third-order elastic constant C e e111 , even differences in absolute values exceed 10 GPa for C11 and 100 GPa for C at very dense k-points meshes, the relative differences are, in fact quite small. Indeed, we would like to note that the relative dispersion is around 3% for 26 × 26 × 26 mesh size and higher k-mesh densities for both second- and thirdorder elastic moduli. Figure 4 depicts the similar convergence test carried out for e11 , and the difference between the cut-off energy. The left pair of frames show C e C11 at 400 eV and 500 eV is below 1 GPa, which is a reasonable precision. e111 one would require cut-off energy However to obtain the same accuracy for C of 700 eV or higher. For any considered type of the deformation scheme, the values of k-mesh size and cut-off energy parameters should at least provide the convergence of third order elastic constants value, calculated at highest compression ratio. The final choice of parameters, which equals to 30 × 30 × 30 for k-mesh size and 700 eV e111 at the highest pressure for cut-off energy is dictated by the convergence of C studied. Convergence of third order elastic constant automatically guarantees that lower order elastic constant convergence should be reached as well. 480 470 465 460 455 1725

−4250 −4275 −4300 −4325 −4350 −4375 −4400 −4425

0 GPa

−10500

300 GPa

e111 , GPa C

0 GPa

e11 , GPa C

475

300 GPa

−10750

e111 , GPa C

1700

e11 , GPa C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

−11000

1675

−11250

1650

−11500

1625

−11750

12 14 16 18 20 22 24 26 28 30 kx x ky x kz

12 14 16 18 20 22 24 26 28 30 kx x ky x kz

e11 (green squares, left panels) and C e111 (blue triangles, right Figure 3: Convergence of C panels), obtained using the scheme S1 from Table 1 in terms of k-mesh density. The tests were performed at two compressions: 0 GPa (top) and 300 GPa (bottom). Cutoff energy is set to 224 eV.

11

466.5

-4350

0 GPa

0 GPa

-4400

e11 , GPa C

e111 , GPa C

466.0

-4450

465.5

-4500

465.0

-4550

464.5

-4600 300 GPa

300 GPa

-11140 e111 , GPa C

1698

e11 , GPa C

-11160

1696

-11180

1694

Cutoff energy, eV

900

800

600

700

500

400

300

200

900

700

800

600

500

300

400

-11200 200

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Cutoff energy, eV

e11 (red squares, left panels) and C e111 (purple triangles, right Figure 4: Convergence of C panels) with respect to cut-off energy parameter for 0 GPa (top) and 300 GPa (bottom). e11 converges at lower cut-off energies than C e111 , especially at higher pressure. Note that C These tests were carried out using k-point mesh of 30 × 30 × 30.

2.4. Effect of η range on fit coefficients Another important convergence parameter is the maximum value of distortion parameter η. Each deformation in Table 1 depends on η in a similar way as expressions (2) and (4) depend on δ. Figure 5 shows how fitting coefficients ν1 , ν2 and ν3 vary with η. The figure uses the same notation as in Table 1, and the points are connected as a guide for an eye. While pressures obtained from first-order fit coefficients ν1 , as well as second order elastic constants, obtained from ν2 do not depend on the range of η (except for values higher than 0.06 for schemes S4, S8 and S9), this is not the case for third-order parameter ν3 . There is a significant effect of ηmax on the curvature of these dependencies, especially for schemes S8 and S9, which feature non-zero non-diagonal elements. Due to the specific nature of these schemes (see AppendixB), they are plotted as a function of 2η instead of η, which stretches the dependence along the x-axis. Finding a proper range is a question of balance, and, unlike other calculation parameters, increasing this one would not necessary lead to a better description. The necessity to check for this parameter even for the case of calculation of second-order elastic constants was highlighted before by Golesorkhtabar et al. [37]. One should carefully choose proper η range in order to keep the strain finite, but small. For bcc molybdenum our tests show that, in order to obtain reliable second- and third-order elastic constants, it is necessary to stay in range of ηmax = 0.06 and not go over this parameter, when using schemes S1–S7 from the Table 1. For schemes akin to S8 and S9 with non-zero non-diagonal elements, ηmax should be taken twice as small (see AppendixB for details).

12

0 GPA

3

300 GPA

S1 S2 S3 S4 S5 S6 S7 S8∗ S9∗

2

1

0 −100

ν1 ∝ Pressure, GPa

ν1 (∝ Pressure), GPa

100

−200 −300 −400 −500 −600 −700 −800

0

5500

1200 1100 1000 900 800 700 600 500 400 300 200

5000 4500 4000 ν2 , GPa

ν2 , GPa

−900

3500 3000 2500 2000 1500 1000

5000

0

0 −5000

−2000

−10000

ν3 , GPa

ν3 , GPa

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

−15000

−4000

−20000

−6000

−25000

−8000

−35000

−10000

−30000 −40000 0.01

0.02

0.03

0.04

0.05 η

0.06

0.07

0.08

0.09

−45000

0.01

0.02

0.03

0.04

0.05 η

0.06

0.07

0.08

0.09

Figure 5: Behavior of the fit coefficients in (14) for different schemes given in Table 1. At first, the wider the range, the smoother the resulting curve becomes. However, with further increase of range width, one can see that the functions ill-behaves. Note that schemes S8 and S9 are extremely sensitive to η range, and plotted in stretched mode as a function of 2η because of the specific features of these two schemes (see text and AppendixB for more discussion).

13

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

3. Results and discussion Using the described method, second- and third- order elastic constants were calculated in the range of volumes corresponding to pressures from 0 GPa up to 300 GPa. e 0 and C e44 are shown in FigSecond-order constants of bcc molybdenum, C ure 6. These constants were obtained using distortions from the Table 1, namely S3: (η, −η, 0, 0, 0, 0) and S5: (0, 0, 0, η, 0, 0). Exact expressions for distortion vectors can be found in AppendixB. In the same range of volumes, second-order elastic constants c0 and c44 were also calculated using standard infinitesimal deformations method using volume-conserving deformations (eqs. (2) and (4)). The pressure for the case of infinitesimal deformations method was obtained using equation of state for undeformed bcc Mo. Moreover, the use of volumeconserving deformations allows to neglect the pressure contribution P ∆V to the energy of the deformed system. The perfect agreement between the pressure dependencies of c44 obtained e44 obtained by finite strains method is observed by standard method and C in the bottom panel of Fig. 6. However on top panel of Fig. 6 we note the e 0 from difference between the values of c0 provided with standard method and C finite deformations technique. At 100 GPa the difference of 3 GPa between those two values is observed. It goes up with pressure up to 22 GPa at highest presure of 300 GPa. While the discrepancies are small and do not exceed the expected accuracy, we once again note the perfect agreement of c44 (P ) and e44 (P ), obtained by two different methods. Worth mentioning are the difference C e44 (P ), which is nearly linear, and more shallow C e 0 (P ), which was between C e 0 and c0 from observed previously [21]. We assume, that the difference between C the standard method may be connected with the higher-order contributions by eij (see Equation (1)), which in case of bcc molybdenum at high pressure are implicitly included in values obtained using infinitesimal strains. However the justification of such an assumption requires careful check and comparison, and lies beyond the scope of this methodological paper. 3.1. Comparison with experiment Discussing the absolute values of the SOECs at 0 GPa, our results are quite consistent with the literature data, from both experimental [38, 39, 40] and theoretical [41] sources. The difference between the calculated and experimentally obtained values does not exceed 1% for c11 : our calculations yield 465.9 GPa, while [38] gives 465±0.6 and [39] provides 461.7±7.9 GPa. The obtained value of c12 at zero pressure also shows good agreement with experiment: 163 in our work with 164.7±9.9 in [39] and 163±2 in [38]. The c44 is slightly underestimated in comparison with experiment, which is well known issue [41] for PAW calculations: our calculations provide 100 GPa, 108.4±3.1 is given in [39] and 109.09 ± 0.9 in [38]. Taking into account the experimental error, these results are rather reasonable. Obtained bulk modulus value (260.5 GPa), is in excellent agreement with the experiment (261.8 GPa and 269±9 GPa in [39] and [40] respectively). 3.2. Comparison with other calculations The ab initio (PAW+GGA) calculations results, provided in work [41] are nearly identical to our results, with exception of small disagreement in the value 14

e 0 and c0 , GPa C

350 300 250 200 150 100 50 0

400 350 300 250 200 150 100 50 0

0

e0 C c0 50

100

150 Pressure, GPa

e44 and c44 , GPa C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0

200

e44 C c44 50

100

150 Pressure, GPa

Mo

200

250

300

250

300

Mo

e 0 (top panel, red triangles) and C e44 (bottom panel, Figure 6: Effective elastic constants C orange triangles), calculated using finite deformations method. In addition, results obtained using standard infinitesimal deformations method are depicted as blue circles for c0 and as blue triangles on for c44 . The points are connected with lines as a guide for an eye.

of V0 and bulk modulus, which still can be considered important for high pressure calculations. The pressure dependence of SOECs presented in that paper e44 and 200 GPa for C e 0 . At higher agrees with our data up to 100 GPa for C pressures one observes that there are some deviations between two sets of calculations. There’s a difference in the values of equilibrium volume V0 : 15.78 ˚ A3 3 ˚ in our work compared to 16.01 A (both differ from experimental value [42] of 15.58 ˚ A3 ). The values for bulk modulus B0 are also different. The reference paper is skimping on details, therefore it could also be a completely different equation of state. These two factors lead to a difference in EOS in two calculations and thus will provide different pressure value at a given volume. We 0 note that first pressure derivatives of bulk modulus B0 are very close for two calculations. Another reason is that Koˇci et al. [41] use smaller value of the cutoff energy, which demonstrated to be very important for calculations at high compression ratios. The method used for elastic constants calculation may be a source of mismatch as well. Although it is not directly mentioned in [41], the elastic properties were most probably calculated using infinitesimal strains method, which doesn’t determine higher order elastic constants explicitly. Thus, nonlinear terms can affect the value of SOECs calculated with this technique. The mismatch of such nature is quite possible at conditions where elastic properties peculiarities are present and structural instability is expected and more subtle effects start to play important role. Third-order elastic constants, calculated with finite strains method are shown in Figure 7. In order to calculate them, all 9 deformations from the Table 1 were employed. We refer again to AppendixB for explicit expressions of deformation

15

matrices. There are 7 deformations in the table that have non-zero third-order term, but there are only 6 third-order elastic constants for bcc lattice. Theree456 was obtained by averaging from the values of ν3 (V ) obtained from fore C (0, 0, 0, η, η, η) and (η, 0, 0, η, η, η) distortions. This was done to compensate for bad convergence of ν3 (V ) in terms of η range. e111 , that is roughly 8 times bigger in absolute Left panel of Figure 7 depicts C magnitude than any other third order constant of bcc molybdenum. The right e155 (green contoured tripanel depicts all other third-order elastic constants: C e112 (blue squares), C e123 (pink diamonds), C e144 (orange filled triangles) angles), C e and C456 (blue circles). The points are connected as a guide for an eye. Note that the second and third order elastic constants of Mo presented in Figs. 6 and 7 differ from those published in [22]. The present values are obtained with much higher precision using more advanced methodology, described above. e111 , C e123 , and C e112 , that keep It is known (see for example [22]), that C the growing trend with pressure, are connected with instability, associated with e0 . softening of C 0

0

e111 C

−2000

Mo

−4000

e155 C e112 C e123 C e144 C e456 C

−500

−1000

−6000

Mo

−1500

−8000

−2000

−10000 −12000

Third-Order Elastic Constants, GPa

Third-Order Elastic Constants, GPa

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

−2500

0

50

−3000

100 150 200 250 300 Pressure, GPa

0

50

100 150 200 250 300 Pressure, GPa

Figure 7: Third-order effective elastic constants of molybdenum.

Molybdenum exhibits no peculiarities of elastic behavior at pressure up to e 0 at higher pressures (see 300 GPa, with the exception of slight softening in of C Fig. 6). At the same time there is no evidence of peculiarities in the behavior e44 in the studied pressure range. of C 3.3. Conclusions

Higher order elastic constants can be used to describe lattice dynamics, electronic structure peculiarities, acoustic properties of material, or anharmonic effects in the vicinity of structural transformations. However, their calculations under pressure is not a trivial task. This paper presents a detailed description 16

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

and careful verification of an algorithm for calculation of high-order elastic constants under pressure based on finite deformation technique, applied to molybdenum. We calculate elastic constants of second and third order in a range of pressures from 0 GPa to 300 GPa with the aim to cover all possible caveats in the process of calculation of elastic moduli. One must be sure that calculations are converged. Specifically, attention must be paid to the deformation parameter η which acts as a computational input parameter, even when there is no pressure applied [11, 10]. Good agreement obtained for the second order elastic moduli of Mo calculated using the proposed technique and conventional method of infinitely small displacements demonstrate the reliability of the developed algorithm and ensures its predictive power for the calculated higher order elastic moduli. The methodology is supposed to be useful in studies of the elastic behavior of materials in the vicinity of phase transitions, including electronic phase transitions, where the increased importance of non-linear elastic response can be expected [43, 44]. 3.3.1. Acknowledgments We are grateful for the support provided by the Swedish Foundation for Strategic Research (SSF) program SRL Grant No. 10-0026, the Swedish Research Council grants No 2015-04391 and No 2014-4750, the Swedish Government Strategic Research Area in Materials Science on Functional Materials at Link¨oping University (Faculty Grant SFO-Mat-LiU No 2009 00971), as well as by the Ministry of Education and Science of the Russian Federation in the framework of Increase Competitiveness Program of NUST “MISIS” (No K22016-013) implemented by a govermental decree dated 16 March 2013, No 211 . Calculations were performed on supercomputer cluster at NUST “MISIS”. AppendixA. Energy expansions up to 4th order of η It is possible to expand Eq. (13) up to fourth order of η and this would allow to calculate corresponding fourth-order elastic constants. In this case one should perform additional convergence tests similar to the ones provided in this article for third-order elastic constants, including sensitivity of particular schemes to ηmax value (see subsection 2.4) Since there are 11 fourth-order elastic constants, one requires 11 equations to find all of them. Table 1 has only 9 schemes, therefore two more should be added. Suggested schemes S10 and S11 are shown in Table A.2. Table A.3 shows the corresponidng ν4 fourth-order terms assuming the same deformations as in Table 1 and Table A.2. Deformation scheme N (η1 , η2 , η3 , η4 , η5 , η6 ) S10 (η, −η, 0, η, 0, 0, 0) S11 (η, −η, 0, 0, η, 0)

Deformation energy ν1 η + ν2 η 2 + ν3 η 3 + o(η 4 ) ν1 ν2 ν3 e e e 0 C11 − C12 + 2C44 0 e11 − C e12 + 2C e44 e144 − 2C e155 0 C 2C

Table A.2: Example of additonal deformations necessary to define full set of fourth-order elastic constants and their corresponding energy expansions in the series of η. Fourth order coefficients ν4 are provided in Table A.3

17

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

N S1 S2 S3 S4 S5 S6 S7 S8 S9

(η1 , η2 , η3 , η4 , η5 , η6 ) (η, 0, 0, 0, 0, 0) (η, η, 0, 0, 0, 0) (η, −η, 0, 0, 0, 0) (η, η, η, 0, 0, 0) (0, 0, 0, η, 0, 0) (η, 0, 0, η, 0, 0) (η, 0, 0, 0, 0, η) (0, 0, 0, η, η, η) (η, 0, 0, η, η, η)

S10

(η, −η, 0, η, 0, 0)

S11

(η, −η, 0, 0, η, 0)

ν4 1 e 24 C1111 1 e 1 e 1 e C + 1111 12 3 C1112 + 4 C1122 1 e 1 e 1 e 12 C1111 − 3 C1112 + 4 C1122 3 e 3 e 1 e e 8 C1111 + C1112 + 4 C1122 + 2 C1123 2 e 3 C4444 1 e e1155 + 2 C e C + C 24 1111 3 4444 2 e 1 e e 24 C1111 + C1144 + 3 C4444

e4444 + 12C e4455 4C e1144 + 2C e1155 + 8C e1456 + 2C e4444 + 12C e4455 +C  1 e e e e e e4444 C − 4 C + 3 C + 24 C − 24 C + 8 C 1111 1112 1122 1155 1266 12   1 e e e e e e 12 C1111 − 4C1112 + 6C1122 + 12C1155 − 24C1255 + 8C4444 1 e 24 C1111 

Table A.3: Fourth-order coefficients ν4 corresponding to a same possible set of deformations as provided in Table 1.

AppendixB. Deformation matrices We provide explicit deformation matrices for bcc crystals using schemes in Table 1. These expressions were obtained using Eq. (13) and utilize the same notation scheme as used in Table 1 and provided up to fourth order in η. Each deformation provides finite strain tensor in both element-wise and shorthand Voigt’s notation forms (η11 = η1 , η22 = η2 , η33 = η3 , η23 = η32 = η4 /2, η13 = η31 = η5 /2, η12 = η21 = η6 /2). Note that in this notation there is a factor of two for non-diagonal elements, which is the reason why we recommend to take twice as small ηmax for those schemes where these elements are not zero. Eq. (13), rewritten in compact form for up to 4th order of η can be expressed as follows: 1 1 5 αij = δij + ηij − ηki ηkj + ηrk ηri ηkj − ηkj ηmk ηmn ηni 2 2 8

(B.1)

Below are explicit expansions for αij using the same notation for schemes as in Tables 1, A.2, and A.3. S1. η11 = η (η, 0, 0, 0, 0, 0) α11 = 1 + η − 0.5η 2 + 0.5η 3 − 0.625η 4 α22 = α33 = 1

α12 = α13 = α23 = 0 S2. η11 = η22 = η (η, η, 0, 0, 0, 0) α11 = α22 = 1 + η − 0.5η 2 + 0.5η 3 − 0.625η 4

α33 = 1

α12 = α13 = α23 = 0

18

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

S3. η11 = η; η22 = −η (η, −η, 0, 0, 0, 0) α11 = 1 + η − 0.5η 2 + 0.5η 3 − 0.625η 4

α22 = 1 − η − 0.5η 2 − 0.5η 3 − 0.625η 4 α33 = 1

α12 = α13 = α23 = 0 S4. η11 = η22 = η33 = η (η, η, η, 0, 0, 0) α11 = α22 = α33 = 1 + η − 0.5η 2 + 0.5η 3 − 0.625η 4

α12 = α13 = α23 = 0 S5. η12 = η21 = η (0, 0, 0, η, 0, 0)

α11 = 1 − 0.5η 2 − 0.625η 4

α22 = α11 α33 = 1

α12 = η + 0.5η 3 α13 = α23 = 0 S6. η11 = η; η12 = η21 = η (η, 0, 0, 0, 0, η) α11 = 1 + η − η 2 + 1.5η 3 − 3.125η 4

α22 = 1 − 0.5η 2 + 0.5η 3 − 1.25η 4

α33 = 1

α12 = η − 0.5η 2 + η 3 − 1.875η 4

α13 = α23 = 0

S7. η11 = η; η23 = η32 = η (η, 0, 0, η, 0, 0) α11 = 1 + η − 0.5η 2 + 0.5η 3 − 0.625η 4

α22 = α33 = 1 − 0.5η 2 − 0.625η 4

α12 = α13 = 0

α23 = η + 0.5η 3 S8. η23 = η32 = η13 = η31 = η12 = η21 (0, 0, 0, η, η, η) α11 = α22 = α33 = 1 − η 2 + η 3 − 3.75η 4

α12 = α13 = α23 = η − 0.5η 2 + 1.5η 3 − 3.125η 4 S9. η11 = η; η23 = η32 = η13 = η31 = η12 = η21 = η (η, 0, 0, η, η, η) α11 = 1 + η − 1.5η 2 + 3.5η 3 − 10.625η 4

α22 = α33 = 1 − η 2 + 1.5η 3 − 5.625η 4

α12 = α13 = η − η 2 + 2.5η 3 − 7.5η 4

α23 = η − 0.5η 2 + 2η 3 − 5η 4 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

S10. η11 = −η22 = η; η12 = η21 = η (η, η, 0, 0, 0, η) α11 = 1 + η − η 2 − η 3 − 2.5η 4

α22 = 1 − η − η 2 − η 3 − 2.5η 4

α33 = 1

α12 = η + η 3 α13 = α23 = 0 S11. η11 = −η22 = η; η23 = η (η, −η, 0, η, 0, 0) α11 = 1 + η − 0.5η 2 + 0.5η 3 − 0.625η 4 α22 = 1 − η − η 2 − 1.5η 3 − 3.125η 4 α33 = 1 − 0.5η 2 − 0.5η 3 − 1.25η 4

α12 = α13 = 0

α23 = η + 0.5η 2 + η 3 + 1.875η 4 ˜ijkl and Cijkl AppendixC. Relations between C The relations between C˜ijkl and Cijkl is obtained from the substitution of (7) into (12), and are shown in Table C.4. We would like to reiterate that pressure P is defined in this method as a fitting parameter ν1 in (14), for deformations that feature non-zero diagonal elements in the applied strain tensor. ˜αβ.. and Cαβ.. Table C.4: Relations between C

e11 = C11 − P C e111 = C111 + 3P C e144 = C144 − P C e12 = C12 + P C e112 = C112 − P C e155 = C155 + P C e e e456 = C456 + P C44 = C44 − P C123 = C123 + P C

[1] J. M. Lang, Y. M. Gupta, Experimental Determination of Third-Order Elastic Constants of Diamond, Physical Review Letters 106 (12) (2011) 125502. doi:10.1103/PhysRevLett.106.125502. [2] C. S. G. Cousins, Elasticity of carbon allotropes. I. Optimization, and subsequent modification, of an anharmonic Keating model for cubic diamond, Physical Review B 67 (2) (2003) 024107. doi:10.1103/PhysRevB. 67.024107. [3] U. Aschauer, R. Braddell, S. A. Brechb¨ uhl, P. M. Derlet, N. A. Spaldin, Strain-induced structural instability in FeRh, Physical Review B 94 (1) (2016) 014109. doi:10.1103/PhysRevB.94.014109. [4] A. Hmiel, J. M. Winey, Y. M. Gupta, M. P. Desjarlais, Nonlinear elastic response of strong solids: First-principles calculations of the third-order elastic constants of diamond, Physical Review B 93 (17) (2016) 174113. doi:10.1103/PhysRevB.93.174113.

20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[5] D. G. Clerc, H. Ledbetter, Second-order and third-order elastic properties of diamond: An ab initio study, Journal of Physics and Chemistry of Solids 66 (10) (2005) 1589–1597. doi:10.1016/j.jpcs.2005.05.075. [6] Y. Hiki, Higher Order Elastic Constants of Solids, Annual Review of Materials Science 11 (1) (1981) 51–73. doi:10.1146/annurev.ms.11.080181. 000411. [7] L. D. Landau, E. M. Lifshitz, A. M. Kosevich, L. P. Pitaevski, Theory of Elasticity, 3rd Edition, Vol. 7 of Theoretical Physics, ButterworthHeinemann, Oxford, 1986. [8] O. Hellman, I. A. Abrikosov, S. I. Simak, Lattice dynamics of anharmonic solids from first principles, Physical Review B 84 (18) (2011) 180301. doi: 10.1103/PhysRevB.84.180301. [9] O. Hellman, P. Steneteg, I. A. Abrikosov, S. I. Simak, Temperature Dependent Effective Potential Method for Accurate Free Energy Calculations of Solids, Phys. Rev. B 87 (10) (2013) 104111. doi:10.1103/PhysRevB.87. 104111. [10] M. Lopuszy´ nski, J. Majewski, Ab initio calculations of third-order elastic constants and related properties for selected semiconductors, Physical Review B 76 (4) (2007) 045202. doi:10.1103/PhysRevB.76.045202. [11] H. Wang, M. Li, Ab initio calculations of second-, third-, and fourthorder elastic constants for single crystals, Physical Review B 79 (22) (2009) 224102. doi:10.1103/PhysRevB.79.224102. [12] E. M. Ann, a. S. Peter, P. D. Michael, R. M. Thomas, L. Kevin, Designing meaningful density functional theory calculations in materials science—a primer, Modelling and Simulation in Materials Science and Engineering 13 (1) (2005) R1. doi:10.1088/0965-0393/13/1/R01. [13] L. Dubrovinsky, N. Dubrovinskaia, V. B. Prakapenka, A. M. Abakumov, Implementation of micro-ball nanodiamond anvils for high-pressure studies above 6 Mbar., Nature communications 3 (2012) 1163. [14] L. Dubrovinsky, N. Dubrovinskaia, E. Bykova, M. Bykov, V. Prakapenka, C. Prescher, K. Glazyrin, H.-P. Liermann, M. Hanfland, M. Ekholm, Q. Feng, L. V. Pourovskii, M. I. Katsnelson, J. M. Wills, I. A. Abrikosov, The most incompressible metal osmium at static pressures above 750 gigapascals., Nature 525 (7568) (2015) 226–9. doi:10.1038/nature14681. [15] T. Lukinov, S. I. Simak, A. B. Belonoshko, Sound velocity in shock compressed molybdenum obtained by ab initio molecular dynamics, Physical Review B 92 (6) (2015) 060101. doi:10.1103/PhysRevB.92.060101. [16] A. Bosak, M. Krisch, A. Chumakov, I. Abrikosov, L. Dubrovinsky, Possible artifacts in inferring seismic properties from X-ray data, Physics of the Earth and Planetary Interiors 260 (2016) 14–19. doi:10.1016/j.pepi. 2016.08.007.

21

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[17] A. Dewaele, M. Torrent, P. Loubeyre, M. Mezouar, Compression curves of transition metals in the Mbar range: Experiments and projector augmented-wave calculations, Physical Review B - Condensed Matter and Materials Physics 78 (10) (2008) 1–13. doi:10.1103/PhysRevB.78. 104102. [18] X. Huang, F. Li, Q. Zhou, Y. Meng, K. D. Litasov, X. Wang, B. Liu, T. Cui, Thermal equation of state of Molybdenum determined from in situ synchrotron X-ray diffraction with laser-heated diamond anvil cells., Scientific reports 6 (August 2015) (2016) 19923. doi:10.1038/srep19923. [19] O. M. Krasilnikov, Yu. Kh. Vekilov, I. Yu. Mosyagin, Elastic constants of solids at high pressures, Journal of Experimental and Theoretical Physics 115 (2) (2012) 237–241. doi:10.1134/S1063776112070096. [20] O. M. Krasilnikov, Yu. Kh. Vekilov, I. Yu. Mosyagin, E. I. Isaev, N. G. Bondarenko, Elastic phase transitions in metals at high pressures., Journal of Physics: Condensed Matter 24 (19) (2012) 195402. doi: 10.1088/0953-8984/24/19/195402. [21] O. M. Krasilnikov, M. P. Belov, A. V. Lugovskoy, I, Yu. Mosyagin, Yu. Kh. Vekilov, Elastic properties, lattice dynamics and structural transitions in molybdenum at high pressures, Computational Materials Science 81 (2014) 313–318. doi:10.1016/j.commatsci.2013.08.038. [22] O. M. Krasilnikov, Yu. Kh. Vekilov, A. V. Lugovskoy, I. Yu. Mosyagin, M. P. Belov, N. G. Bondarenko, Structural transformations at high pressure in the refractory metals (Ta, Mo, V), Journal of Alloys and Compounds 586 (2014) S242–S245. doi:10.1016/j.jallcom.2013.05.151. [23] L. Vitos, Computational Quantum Mechanics for Materials Engineers: The EMTO Method and Applications, Springer London Ltd, 2007. [24] D. G. Pettifor, Theoretical predictions of structure and related properties of intermetallics, Materials Science and Technology 8 (4) (1992) 345–349. doi:10.1179/mst.1992.8.4.345. [25] D. C. Wallace, Thermoelastic Theory of Stressed Crystals and HigherOrder Elastic Constants, in: F. S. Henry Ehrenreich, D. Turnbull (Eds.), Solid State Physics, Vol. Volume 25, Academic Press, 1970, pp. 301–404. [26] G. Kresse, J. Furthm¨ uller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Computational Materials Science 6 (1) (1996) 15 – 50. doi:http://dx.doi.org/10.1016/ 0927-0256(96)00008-0. [27] G. Kresse, J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B 47 (1993) 558–561. doi:10.1103/PhysRevB.47.558. [28] G. Kresse, J. Hafner, Ab initio molecular-dynamics simulation of the liquidmetal–amorphous-semiconductor transition in germanium, Phys. Rev. B 49 (1994) 14251–14269. doi:10.1103/PhysRevB.49.14251.

22

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[29] G. Kresse, J. Furthm¨ uller, Efficient iterative schemes for ab initio totalenergy calculations using a plane-wave basis set, Phys. Rev. B 54 (1996) 11169–11186. doi:10.1103/PhysRevB.54.11169. [30] J. P. Perdew, K. Burke, M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77 (18) (1996) 3865–3868. doi: 10.1103/PhysRevLett.77.3865. [31] G. Kresse, D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59 (3) (1999) 1758–1775. doi: 10.1103/PhysRevB.59.1758. [32] T. S. Duffy, G. Shen, J. Shu, H.-K. Mao, R. J. Hemley, A. K. Singh, Elasticity, shear strength, and equation of state of molybdenum and gold from x-ray diffraction under nonhydrostatic compression to 24 GPa, Journal of Applied Physics 86 (12) (1999) 6729. doi:10.1063/1.371723. [33] W. Liu, Q. Liu, M. L. Whitaker, Y. Zhao, B. Li, Experimental and theoretical studies on the elasticity of molybdenum to 12 GPa, Journal of Applied Physics 106 (4) (2009) 043506. doi:10.1063/1.3197135. [34] J. H. Nguyen, M. C. Akin, R. Chau, D. E. Fratanduono, W. P. Ambrose, O. V. Fat’yanov, P. D. Asimow, N. C. Holmes, Molybdenum sound velocity and shear modulus softening under shock compression, Physical Review B 89 (17) (2014) 174109. doi:10.1103/PhysRevB.89.174109. [35] J. Wang, F. Coppari, R. F. Smith, J. H. Eggert, A. E. Lazicki, D. E. Fratanduono, J. R. Rygg, T. R. Boehly, G. W. Collins, T. S. Duffy, Xray diffraction of molybdenum under ramp compression to 1 TPa, Physical Review B 94 (10) (2016) 104102. doi:10.1103/PhysRevB.94.104102. [36] R. S. Hixson, J. N. Fritz, Shock compression of tungsten and molybdenum, Journal of Applied Physics 71 (4) (1992) 1721. doi:10.1063/1.351203. [37] R. Golesorkhtabar, P. Pavone, J. Spitaler, P. Puschnig, C. Draxl, Elastic: A tool for calculating second-order elastic constants from first principles, Computer Physics Communications 184 (8) (2013) 1861–1873. doi:10. 1016/j.cpc.2013.03.010. [38] K.-H. Hellwege, O. Madelung (Eds.), Numerical Data and Functional Relationships in Science and Technology, Vol. 18 of Group {III}: Crystal and Solid State Physics, Springer-Verlag, Berlin Heidelberg New York Tokyo, 1984. [39] F. F. Voronov, V. M. Prohorov, E. L. Gromnickaya, G. G. Ilina, Second and third order elastic moduli of molybdenum monocrystal, The Physics of Metals and Metallography 45 (6) (1978) 94105. doi:http://impo.imp. uran.ru/fmm/Electron/vol45\_6/main.html. [40] Y. Zhao, A. C. Lawson, J. Zhang, B. I. Bennett, R. B. Von Dreele, Thermoelastic equation of state of molybdenum, Phys. Rev. B 62 (13) (2000) 8766–8776. doi:10.1103/PhysRevB.62.8766.

23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[41] L. Koˇci, Y. Ma, A. Oganov, P. Souvatzis, R. Ahuja, Elasticity of the superconducting metals V, Nb, Ta, Mo, and W at high pressure, Physical Review B 77 (21) (2008) 214101. doi:10.1103/PhysRevB.77.214101. [42] K. W. Katahara, M. H. Manghnani, E. S. Fisher, Pressure derivatives of the elastic moduli of BCC Ti-V-Cr, Nb-Mo and Ta-W alloys, Journal of Physics F: Metal Physics 9 (5) (1979) 773–790. doi:10.1088/0305-4608/9/5/006. [43] M. I. Katsnelson, I. I. Naumov, A. V. Trefllov, Singularities of the electronic structure and pre-martensitic anomalies of lattice properties in β-phases of metals and alloys, Phase transitions 49 (1994) 143–191. doi:10.1080/ 01411599408201172. [44] Ya. M. Blanter, M. I. Kaganov, A. V. Pantsulaya, A. A. Varlamov, The theory of electronic topological transitions, Physics Reports 245 (4) (1994) 159–257. doi:10.1016/0370-1573(94)90103-1.

24