cij: A Python code for quasiharmonic thermoelasticity
Chenxing Luo, Xin Deng, Wenzhong Wang, Gaurav Shukla, Zhongqing Wu, Renata M. Wentzcovitch
11 cij : A Python code for quasiharmonic thermoelasticity Chenxing Luo a , Xin Deng b , Wenzhong Wang b,c , Gaurav Shukla d , Zhongqing Wu b,e,f , Renata M. Wentzcovitch a,g,h,* a Department of Applied Physics and Applied Mathematics, Columbia University, New York, NY 10027, USA b Laboratory of Seismology and Physics of Earthβs Interior, School of Earth and Space, University of Science and Technology of China, Hefei, Anhui 230026, China c Department of Earth Sciences, University College London, London WC1E 6BT, United Kingdom d Department of Earth Sciences, Indian Institute of Science Education and Research Kolkata, Mohanpur, West Bengal, India e CAS Center for Excellence in Comparative Planetology, USTC, China f National Geophysical Observatory at Mengcheng, USTC, China g Department of Earth and Environmental Sciences, Columbia University, New York, NY 10027, USA h LamontβDoherty Earth Observatory, Columbia University, Palisades, NY 10964, USA * Corresponding author at: Lamont-Doherty Earth Observatory, Columbia University in the City of New York, 61 Route 9W, Palisades, NY 10964, USA. E-mail address: [email protected] (R.M. Wentzcovitch).
Abstract
The Wu-Wenzcovitch semi-analytical method (SAM) to compute the thermoelastic tensor (Cij) (SAM-Cij)is a concise and predictive formalism to calculate the high-pressure and high-temperature (high- PT ) Cij of crystalline materials. This method has been successfully applied to materials across different crystal systems in conjunction with ab initio calculations to compute static elastic coefficients and phonon frequencies. Such results have offered first-hand insights into the composition and structure of the Earthβs mantle.Here we introduce the cij package, a Python implementation of the SAM-Cij formalism. It enables a thermoelasticity calculation to be initiated from a single command and fully configurable from a calculation settings file to work with solids within any crystalline system. These features allow SAM-Cij calculations to work on a personal computer and to be easily integrated as a part of high-throughput workflows. Here we show the codeβs performance for three minerals from different crystal systems at their relevant PT s: diopside (monoclinic), akimotoite (trigonal), and bridgmanite (orthorhombic). Keywords: thermoelasticity; acoustic velocity; diopside; akimotoite; bridgmanite; high pressure and high temperature
PROGRAM SUMMARY
Program title: cij
Developer's repository link: https://mineralscloud.github.com/cij
Licensing provisions:
GNU General Public License 3 (GPL)
Programming language:
Python 3
Nature of problem (approx. 50-250 words):
Experimental measurements of full elastic tensor coefficients under high-pressure and high-temperature conditions are challengingand susceptible to uncertainties. Computations of thermoelastic coefficients with the conventional density functional theory (DFT) plus quasiharmonic approximation (QHA) or molecular dynamics (MD) methods are computationally extremely demanding,especially for materials with low symmetries because of the revaluation of free energy for strained configurations.
Solution method (approx. 50-250 words):
Based on a semi-analytical method proposed by Wu and Wentzcovitch [1], the present code only needs static-state elastic coefficients and phonon vibrational density of states for equilibrium configurations as input to calculate the thermal elasticity. This method avoids the reevaluation of free energy for strained configurations and can be applied to all crystal systems.
Reference: [1] Z. Wu, R.M. Wentzcovitch, Phys. Rev. B 83 (2011) 184115.
1. Introduction
Elasticity is a fundamental property of solids that characterizes their mechanical response to external stress. Determination of elastic coefficients, especially at high pressures and temperatures ( PT ), has dramatically advanced geophysics. However, despite recent methodological developments, measurement of the full elastic tensor at high PT has remained a challenging undertaking and susceptible to uncertainties [1,2]. Although computational methods are regularly resorted to as alternative, fully numerical ab initio approaches based on the density functional theory (DFT), such as DFT + QHA (quasiharmonic approximation) (e.g., Ref. [3,4]) or DFT + MD (molecular dynamics) (e.g., Ref. [5]), are computationally demanding, considering the numerous strained configurations involved, especially for crystals with lower symmetries [6].To overcome such a computational challenge, a semi-analytical method (SAM) was proposed to compute the thermoelastic tensor (Cij) (hereafter SAM-Cij) [7]. Compared to the traditional QHA approach, SAM-Cij adopts an analytical expression for the thermal part of Cij to circumvent the reevaluation of vibrational density of states (VDoS) for slightly strained configurations, which drastically reduces the calculation cost by at least one order of magnitude. This formalism offers an overall improved agreement with experimental measurements for high- ππ elasticity compared to the fully numerical approach [7]. Such improvement is possibly a benefit from the imposed isotropic thermal pressure.Other formalisms and codes are also proposed to computationally resolve thermoelasticity, among which some are under active development [8β11]. However, those formulations either do not attempt to avoid tedious phonon or MD calculations for strained configurations [9], which becomes nearly unfeasible for low symmetry structures, or oversimplify calculations to address extreme- ππ conditions of geophysical interests (such as Ref. [8,10,11]). In contrast, the SAM-Cij formalism remains predictive up to the ππ boundary of validity of the QHA and has already been extensively tested forlower mantle minerals at their relevant conditions [12β15]. It has also been recently extended and validated for low crystal symmetries such as monoclinic and trigonal[16,17].Here we introduce the cij package, a Python implementation of SAM-Cij. Unlike some other thermoelasticity calculation methods (such as Ref. [9,18]), this package is decoupled from a particular DFT software suite. As a standalone package, cij requiresonly total static energies, VDoS, and static Cij at a series of volume points as input, obtainable with most DFT software suites. One can initiate a SAM-Cij calculation from a single command and configure it within a single settings file to work with materials across different crystal systems. Since most cij calculations only need a few minutes to complete on a desktop-level computer, high-performance computing (HPC) setup is not imperative. Therefore, this package is easy to use on a personal computer and is ready for integration into high-throughput workflows.This paper is organized as follows: the next section briefly reviews the SAM-Cij method; Secs. 3 and 4 describe the structure and usage of cij ; Sec. 5 shows its application to systems of different symmetries: diopside (monoclinic), bridgmanite (orthorhombic), and akimotoite (trigonal); Sec. 6 summarizes the paper.
2. The SAM-Cij formalism
The isothermal elastic tensor elements, or elastic coefficients, π πππππ , are second-order strain derivatives of the Helmholtz free energy πΉπ πππππ = 1π ( π πΉππ ππ ππ ππ ) + 12 π(2πΏ ππ πΏ ππ β πΏ ππ πΏ ππ β πΏ ππ πΏ ππ ), (1) the QHA Helmholtz free energy πΉ(π, π, π) under strain state π is [19β22] πΉ(π, π, π) = π st (π, π) + β 12 ππ βπ ππ (π, π)+βπ π΅ π β ln ππ [1 β exp [β βπ ππ (π ππ , π)π B π ]] (2) where π ππ are phonon frequencies of the π -th mode at the π -th wave-number. The first, second, and third terms on the r.h.s. of Eq. (2) are respectively the static total energy, π st (π, π) , the zero-point energy, πΈ zpm (π, π) , and thermal excitation energy. The vibrational or phonon energy is πΈ ph (π, π, π) = πΈ zpm (π, π) + πΈ th (π, π, π). (2 β² ) The adiabatic elastic moduli π πππππ can be converted from the isothermal one, π πππππ , as [23] π πππππ = π πππππ + πππΆ π ππ ph ππ ππ ππ ph ππ ππ πΏ ππ πΏ ππ (3) where π ππ ( π, π = 1,2,3 ) are infinitesimal strains and π πβ (π, π, π) is the phonon entropy at relevant strain state. Computing the strain derivatives in Eq. (1) requires knowledge of the variation of mode frequencies π ππ w.r.t. strains π ππ , namely strain GrΓΌneisen parameters, πΎ ππππ ππ ππ π ππ = βπΎ ππππ π ππ . (4) SAM-Cij avoids the expensive frequency calculations for strained configurations by deriving analytical relationships between the mode average of πΎ ππππ (π) and the mode average of volume-GrΓΌneisen parameters, πΎ ππ (π) , which can be readily obtained from phonon calculations under hydrostatic compression. The derivation makes QHA thermal stresses hydrostatic, which is only an approximation for anisotropic materials but has the beneficial effect of producing the thermal component of Cij under this desirable stress condition. The mode-averaged strain-GrΓΌneisen parameters necessary for longitudinal and off-diagonal Cij calculations are given by πΎ ππ = π + π + π ππ πΎ , (5) where the averages are over all ππ vibrational modes; π ππ ( π = 1,2,3 ) are longitudinal lattice strains produced under static hydrostatic compression. They express the crystal anisotropy ignored in the isotropic approximation. Similarly, their products and derivatives are given by πΎ ππ πΎ ππ = { 15 (π + π + π ) π ππ π ππ (πΎ) if π = π115 (π + π + π )π ππ π ππ (πΎ) if π β π (6)ππΎ ππ ππ ππ = { 15 (π + π + π ) π ππ π ππ π ππΎππ if π = π115 (π + π + π )π ππ π ππ 2 π ππΎππ if π β π (7)
The thermoelastic coefficients can be analytically expressed using mode-averaged strain-GrΓΌneisen parameters, their products, and derivatives. The isothermal elastic coefficients π πππππ (π, π) = π ππππst (π) + π ππππph (π, π) are the sum of static and phonon contributions (Eqs. 1-2). π ππππst (π) is obtained by straightforward static DFT calculations of stress vs. strain relations. In the next two sub-sections, we show how SAM-Cij evaluates π ππππph (π, π) for longitudinal ( π ππππph , π = 1,2,3 ), off-diagonal ( π ππππph , π, π = 1,2,3 , π β π ) or shear ( π ππππph , π β π or π β π ) elastic coefficients. For longitudinal and off-diagonal terms of the elastic tensor, Eq. (1) reduces to π πππππ = 1π ( π πΉππ ππ ππ ππ ) + π(π, π)(1 β πΏ ππ ). (8) The phonon contribution terms are therefore π ππππph (π, π) = π ππππzpm (π) + π ππππth (π, π) + (1 β πΏ ππ )π ph (π, π), (9) where the phonon contribution to the pressure is π ph (π, π) = π(π, π) β π st (π) . Combining Eqs. (2) and (4), the zero-point term becomes π ππππzpm = β2π β π π ππ (π)ππ ππ ππ ππππ = β2π β (πΎ ππππ πΎ ππππ β ππΎ ππππ ππ ππ + πΏ ππ πΎ ππππ ) ππ π ππ . (10) The thermal contribution π ππππth is π ππππth (π, π) = π π΅ ππ β π [ ln (1 β π βπ ππ )]ππ ππ ππ ππππ = π π΅ ππ β [ ππ β π ππ2 π π ππ (π π ππ β 1) πΎ ππππ πΎ ππππ + π ππ (π π ππ β 1) (πΎ ππππ πΎ ππππ β ππΎ ππππ ππ ππ + πΎ ππππ πΏ ππ )] (11) where π ππ = βπ ππ /ππ . These two expressions are then simplified by using the mode-averaged values given in Eqs. (5-7) (see Ref. [7] for details). The ππ/ππ ππ term necessary to compute Eq. (3) is given by ππππ ππ = π B β π ππ2ππ π π ππ (π π ππ β 1) πΎ ππππ . (12) For the elastic tensor components that have not been addressed so far, SAM-Cij employsaxis rotations to convert shear strains back to longitudinal or off-diagonal ones discussed above. To solve for π ππππ , this SAM-Cij implementation applies a symmetric strain, π ππππ , in a rotated crystal system with components given by π πΌπ½ππππ = [1 β (1 β πΏ ππΌ πΏ ππ½ )(1 β πΏ ππΌ πΏ ππ½ )(1 β πΏ ππΌ πΏ ππ½ )(1 β πΏ ππΌ πΏ ππ½ )]π. (13) For example, π = [0 0 00 0 π0 π 0]βπ = [0 π 0π 0 π0 π 0]βπ = [π 0 π0 0 0π 0 0]. It is always possible to find a rotation matrix, π , that diagonalizes the symmetric tensor π ππ , i.e., a real orthonormal matrix π that gives π β1 π ππππ π = π πβ²πβ²πβ²πβ² , where π πβ²πβ²πβ²πβ² is diagonal. Under this rotation, the invariance of strain energy gives: β π πΌπ½πΎπΏπΌπ½πΎπΏ π πΌπ½ππππ π πΎπΏππππ = β π πΌ β² πΌ β² π½ β² π½ β² πΌ β² π½ β² π πΌ β² πΌ β² π β² π β² π β² π β² π π½ β² π½ β² π β² π β² π β² π β² , (14) where πΌ, π½, πΎ, πΏ, πΌβ², π½β² = 1,2,3 .We note that the r.h.s. of the Eq. (14) contains longitudinal or off-diagonal terms only, which can be analytically resolved as is discussed in subsection 2.3.1 with rotated strains π πβ²π π ππ π ππβ² = π πβ²πβ² (π, π, πβ², πβ² = 1,2,3) containing negligible off-diagonal terms when π : π : π β 1: 1: 1 .The l.h.s. of Eq. (14), being a little more complicated, will fall into one of two scenarios:1. For π ππππ -like terms ( π β π ), the only term we have on the l.h.s of Eq. (14) is the unknown term, so the equation is solvable.2. For other π ππππ , the l.h.s. is a combination of π ππππ -like terms ( π β π ) (solved in situation 1), π ππππ -like terms (solved analytically in subsection 2.3.1), and the unknown term π ππππ . So, again, this is solvable.A recursive algorithm is currently implemented to solve these shear terms.
3. The cij distribution
The cij package is written in Python 3. After decompressing the cij.zip zip file, onesees the Python source code in the cij sub-folder, input for three examples in the examples sub-folder, documentation in the docs sub-folder, and the installation script setup.py . The cij package runs on all major platforms supported by the qha package .The Python code is organized into several modules. A description of essential modules and scripts are shown in Table 1.
The package can be installed with the pip package manager: at the command prompt, one should navigate to the directory that contains the setup.py script and execute " pip install . ". Then, the package should be ready for use.
After preparing the input files (to be discussed in Sec. 4), one can navigate to the directory that contains the YAML settings file (hereafter settings.yaml ) and execute 0" cij run settings.yaml " to perform the calculation. This command invokes the main.py script under the cij/cli directory. The flowchart in Fig. 1 will help understandthe procedure of a SAM-Cij calculation.
A typical thermoelasticity calculation for an orthorhombic crystal with the cij command-line program finishes in less than one minute on a desktop computer. Each output variable specified in the output section of settings.yaml will be saved to a separate file with the same tabular (π, π) - or (π, π) -grid format as in the qha code [21] .The available output variables are listed in Table 2.The cij package provides three utilities to inspect calculation results right out-of-the-box: cij plot converts a data table to a PNG plot; cij extract extracts data from the original (π, π) table files and prepares data table with multiple variables at specified π or π for further analysis, e.g., table with π ππ βs, πΎ , and πΊ vs. π at 300 K; cij extract-geotherm extracts data and creates a data table along ππ of a geotherm. Detailed documentation of this program will be available online at https://mineralscloud.github.io/cij/. The source of this documentation is located in the docs sub-folder and can be built locally with Sphinx.
4. Input files
At the beginning of a calculation, the cij run program reads the settings.yaml file and two input data files that contain phonon data (hereafter input01 ) and static elasticity data (hereafter elast.dat ). Instructions on how to prepare these files are given below(see Secs. 4.1 to 4.3). settings.yaml ) The settings.yaml file is home to all calculation settings. One needs to specify calculation parameters, such as thermal EoS fitting parameters, phonon interpolation 1settings, input data location, and output variables to store. The available parameters and their detailed descriptions are listed in Table 3. input01 ) The QHA input data file contains the static energy and phonon frequencies at various volume points. The general structure of this file is identical to the one used by the qha program as described in Ref. [21], but the number of formula units ( nm ) and atoms ( na ) need to be additionally appended to the end of the fourth line, after the number of volumes ( nv ), π -points ( nq ), and modes ( np ). The ordering of phonon mode frequencies should be matched between different volume points according to the mode symmetry to ensure proper interpolation. The package also includes a cij modes utility that plots the interpolated frequencies vs. volume at a given π -point. elast.dat ) The static elasticity input data file tabulates the static elastic coefficients ( π ππst , π, π = 1,2,3 ) and axial length along three axes in Cartesian coordinates ( π ππ , π = 1,2,3 ) at a series of volume points. This file format is specified in Table 4. The required elastic tensor components for each crystal system can be found in Ref. [24,25]. To compute Voigt-Reuss-Hill (VRH) averaged elastic moduli and acoustic velocities, unless all non-zeroterms are listed one needs to either specify the crystal system in settings.yaml ormanually preprocess the elast.dat with the cij fill utility to generate a new elast.dat which contains all non-zero terms. The column names in this table also determine the output thermoelastic coefficient terms.
5. Examples
Here we show the high- ππ elasticity of three important minerals in geophysics: diopside, akimotoite, and bridgmanite. These materials' thermoelastic properties have been well-studied with SAM-Cij in Ref. [13,14,16,17]. Here, we revisit these minerals using thenew cij package to demonstrate its reliability.2 Diopside, the primary host of Ca in the upper mantle, is a rock-forming pyroxene mineral with a chemical composition of MgCaSi O . Its structure belongs to the monoclinic crystal system, with a πΆ2/π space group. The elastic tensor of diopside contains 13 independent terms. Results shown here use the local-density approxiamations (LDA)exchange-correlation functional [26]. Details of these DFT calculations are given inRef. [17] and Appendix A.Isobaric and isothermal equations of state (EoS) for diopside are shown in Figs. 2 (a, b). LDA + QHA reproduces the EoS measured at high- π [27], high- π [28], and high- ππ [29]with sufficient accuracy. The calculated thermal expansivity πΌ in Fig. 2 (c) shows no obvious superlinear π dependence up to 1500 K at high pressurs. However, at 0 GPa, an inflection point develops gently at ~1000 K. At higher pressures inflection points develop at approximately π = 1170 + 37 π ( π in K, π in GPa). At pressures higher than these, the validity of the QHA might be questionable. Fig. 3 shows the ππ dependence of individual elastic coefficients of diopside. The nine elastic coefficients of orthorhombic systems increase with π and decrease with π , the other four, π , π , π , and π decrease with π and increase with π . In terms of π -dependence, our results in Fig. 3 (aβc) are consistent with 300 K measurements fromRef. [27]. There are somewhat larger discrepancies for off-diagonal and shear πΆ ππ terms, but after contrasting them with the 2β3% experimental uncertainties and the more significant 24% overestimation seen in previous MD simulations [30], the discrepancies are relatively insignificant. The comparison between the calculated and measure [37] π -dependence of πΆ ππ terms made in Fig. 3 (dβf) shows good agreement up to 1000 K, except that π is systematically underestimated by ~10 GPa, and π has contradictory π -dependence against Ref. [28]. At π > 1000 K, measurements [28] continue to change linearly with π , while our results start deviating. This is likely caused by anharmonic effects and this behavior is consistent with that of the inflection point in the QHA thermal expansion coefficient shown in Fig. 2 (c).3Fig. 4 shows positive π dependence, and negative π dependence of the VRH averagedadiabatic bulk and shear moduli ( πΎ π and πΊ ) and compressional and shear velocities ( π£ p and π£ s ). Our results agree well with high- ππ ultrasonic measurements on polycrystalline samples from Ref. [29]. Compared to Ref. [27,28], the similar deviating behavior seen in Fig. 3, which is likely caused by anharmonicity, is also observed here at π beyond that of the inflection points in the QHA thermal expansion coefficient shown in Fig. 2 (c). MgSiO akimotoite is a high- π polymorph of pyroxene. It is stable at transition zone and uppermost lower mantle conditions in the Earth [31]. Its strong elastic anisotropypredicted by static calculations [32] makes akimotoite an outstanding candidate for the source of large acoustic-wave anisotropy observed at the bottom of the transition zone[33,34]. A recent study [16] used SAM-Cij to investigate the ππ dependence of its anisotropy. It has an ilmenite-like structure and has an π 3βΎ space group with trigonal symmetry. Its elastic properties are fully described by 7 independent elastic tensor components. DFT details are described in Ref. [16] and Appendix A.Fig. 5 (a, b) shows the
πππ
EoS of akimotoite. Our LDA + QHA EoS has compatible ππ dependence compared to measurements at high- π [35] and at high- ππ [36]. The systematic 1% overestimation in volume compared to measurements [35,36] is common for LDA results and can be easily reconciled with EoS corrections, if necessary [37,38].The present LDA results are also compared to high- ππ generalized gradient approximation (GGA) [39] + MD EoS results from Ref. [40] after their proposed EoS correction. The comparatively superior agreement of the LDA + QHA EoS with measurements (less than 2 GPa difference) over their uncorrected GGA + MD EoS (a β π necessary to match experimental data in Ref. [41]) justifies thechoice of LDA to study thermoelasticity. This has been the case since the first elasticity calculations [3,4]. Fig. 5 (c) shows the π dependence of πΌ at various π βs. The inflection points in πΌ vs. π at ~ 1500β2000 K (approximated by π = 22.5 π + 1400 , π in K, π in GPa) and the πΌ βs superlinear dependence of π beyond this boundary suggest that the 4QHA may be unreliable and anharmonicity might start impacting these results. Beyondthis boundary, results should be treated with caution.Fig. 6 shows the akimotoiteβs elastic coefficients π ππ as functions of π and π . Here, π , π , π , π and π , increase with π and decrease with π ; π and π decrease with π and increase with π . The only major conflict here is the inverted sign of π in Ref. [35],which, according to Ref. [35] is caused by differences in crystal setting; once inverted, their results and ours are consistent. Other than that, our results agree well with 300 KBrillouin spectroscopy measurements in Ref. [35,41]. The smaller π and π , and slightly larger π have discrepancies comparable with experimental uncertainties. Theaforementioned EoS correction [37,38], if applied, would increase π and π , which mitigates the discrepancies further.Fig. 7 shows the ππ dependence of πΎ π , π£ p , πΊ and π£ s . Compared to Ref. [35], our 300 K πΎ π and π£ p values agree soundly, but πΊ and π£ s depart downwards from measurements by ~2 % when compressed. A similar feature in the π -dependence exists in a previous LDA calculation [32], so this is a consistent LDA prediction. The earlier reports of stiffer πΎ π and πΊ , and faster π£ p and π£ s measured at high- ππ with ultrasound [36] compared to both our results and Ref. [35], it is likely caused by a partially transformed sample, according to Ref. [35]. These measurements in Ref. [36] are, thus, deemed unreliable. Nevertheless, the π -gradient of their πΎ π and π£ p are mostly aligned with ours but the that of their πΊ and π£ s are slightly larger than to ours. MgSiO -perovskite (Mg-Pv) is the Mg-endmember of bridgmanite, the most abundant mineral in the Earth's lower mantle. The thermal πΆ ππ tensor of Mg-Pv calculated with SAM-Cij was used as a reference for comparison with its iron-bearing counterparts to understand the effect of alloying and iron spin-crossover [13,14,42]. Experimentaldetermination of the thermal πΆ ππ tensor, especially at high- ππ , is involved with substantial certainties [2]. High- ππ experimental data for pure Mg-Pv has not beenpublished yet.5Mg-Pv has a Pbmn space group with orthorhombic symmetry. 9 individual elastic coefficients are required to describe its elastic properties. Calculations reported here werecarried out with LDA. Calculation details are described in Ref. [13] and Appendix A.Figs. 8 (a, b) show the
πππ
EoS of Mg-Pv. LDA + QHA reproduces the EoS obtained with XRD-DAC measurements at high- π [43β45] and high- ππ [44,46,47] faithfully. Compared to our LDA + QHA EoS, GGA + MD simulations [5] report a roughly less than 5% overestimated π of π, π , due to GGAβs under-binding. The reliable compression curves here allow us to proceed to calculate πΆ ππ at high- ππ . Fig. 8 (c) plots πΌ as a function of π at several π βs. Emperically defined by the non-superlinear dependence of πΌ on π , the region outlined by π πΌ/ππ | π = 0 , the black line, corresponds to the ππ -range of QHA validity [4,48]. The ππ validity range generally resembles that of Ref. [48] and small variations on this boundary are caused by numerical errors in high-order π -derivatives, interpolation, and different choices of ππ -grids. Fig. 9 shows elastic coefficients π ππ of Mg-Pv as functions of π and π . These π ππ βs increase almost linearly with π and decrease linearly with π . Our results agree well with those in Ref. [49] determined at ambient conditions. GGA + MD calculations fromRef. [5] disagree more with measurements. The scarcity of high- π or high- π measurements of elastic coefficients on pure Mg-Pv single-crystal does not allow further comparisons here, which shows why SAM-Cij results are crucial.Fig. 10 shows the ππ dependence of VRH-averaged πΎ π , πΊ , as well as π£ s , and π£ p . Derived from πΆ ππ , these properties show uniformly positive π -dependence and negative π -dependence. Our results can be verified against 300 K and 2700 K ultrasonic measurements from Refs. [45,50] up to 100 GPa. Although Refs. [44,46] suggest a slightly larger ππ -gradient within a narrower ππ range, i.e., 0β20 GPa, up to 1200 K, the inconsistencies among these measurements are either comparable or more significant than their deviation from our results. Compared to the MD simulations in Ref. [5], our SAM-Cij calculation is not only less time-consuming, but also offer much-improved consistency with experimental measurements.6
6. Conclusion
In summary, this paper presented cij , an easy-to-use Python package that calculates thermal Cij, elastic moduli, and acoustic velocities for crystalline materials at high- ππ based on the SAM-Cij formalism. The code presented here is tested on three minerals with different crystal symmetry. Consistency between our high- ππ results with measurements highlights the performances of the code.
7. Acknowledgments
This work was supported by DOE DESC0019759; National Natural Science Foundation of China 41925017; and IISER-K Start-up Research Grant. 7
Appendix A. DFT details
All DFT calculations were performed using the Quantum ESPRESSO [51] and the LDAexchange-correlation functional [26]. Detailed calculation parameters for these three minerals are described in the next three subsections.
Diopside
Calculations on diopside were performed using norm-conserving pseudopotentials. For Mg, the pseudopotential was generated using the von BarthβCarβs method [52,53], for Siand O the Troullier-Martins method [54] was used, for Ca, an ultrasoft pseudopotential[55] was used. The plane-wave kinetic energy cutoff was 70 Ry. Structural optimizations at 7 different π s were performed using variable cell-shape damped molecular dynamics(VCS-MD) [56,57] with a 2 Γ 2 Γ 2 π -point mesh. Dynamical matrices for optimized structures were obtained using density functional perturbation theory (DFPT) [58] on a 2Γ 2 Γ 2 π -point mesh grid. Force constants obtained from these dynamical matrices were later interpolated to an 8 Γ 8 Γ 8 mesh grid to obtain the VDoS. Akimotoite
For akimotoite, the pseudopotential of Mg was generated using the BarthβCarβs method,and the pseudopotentials of O and Si were generated using Troullier-Martinsβ method. The plane wave cutoff energy was 70 Ry. The structure of akimotoite were optimized using the VCS-MD. with a 4 Γ 4 Γ 4 π -point mesh at 8 different π s. The dynamical matrices for akimotoite were calculated using DFPT [58] on a 2 Γ 2 Γ 2 π -point mesh and then extrapolated to a denser 4 Γ 4 Γ 4 π -mesh to obtain VDoS. Bridgmanite
Calculations on Mg-Pv were performed on a 40-atom supercell. Ultrasoft [61] pseudopotentials were used for Al, Fe, Si, and O. A norm-conserving pseudopotential generated with von Barth-Carβs method was used for Mg. The electronic states weresampled on a 2 Γ 2 Γ 2 π -point grid with a plane-wave kinetic energy cutoff of 40 Ry, respectively. The Mg-Pv structure was optimized at 10β12 relevant π points with VCS-8MD. Dynamical matrices were calculated using DFPT on a 2 Γ 2 Γ 2 q -point grid for these structures and then interpolated to an 8 Γ 8 Γ 8 π -point grid to obtain VDoS.9 Figures
Figure 1. The flowchart for thermoelasticity calculations with SAM-Cij.0Figure 2. (a, b) Unit-cell volume (2 f.u.) of diopside (a) vs. π at various π s, and (b) vs. π at various π s. Teal diamonds in (a) correspond to measurements reported in Ref. [27] at 300 K, solid purple circles in (b) correspond to measurements reported in Ref. [28] at 0 GPa. Colored squares in (a, b) correspond to measurements reported in Ref. [29]. In (a) π and (b) π are represented by colors in the color-bars. (c) Thermal expansivity of diopside vs. π at various π s.1Figure 3. Elastic tensor components ( π ππ ) of diopside vs. (aβc) π at various π , and (dβf) vs. π s at various π s. Teal diamonds in (aβc) correspond to measurements reported in Ref. [27] at 300 K, solid purple circles in (dβf) correspond to measurements reported in Ref. [28] at 0 GPa.2Figure 4. (a, c) Elastic moduli and (c, d) acoustic velocities of diopside vs. (a, b) π and (c, d) π . Teal diamonds in (a, b) correspond to measurements reported in Ref. [27] at 300 K, solid purple circles in (c, d) correspond to measurements in Ref. [28] at 0 GPa. Colored squares correspond to calculations in Ref. [29] at (a, b) π and (c, d) π represented by colors in the color-bars. 3Figure 5. (a, b) Unit-cell volume (6 f.u.) of akimotoite (a) vs. π at various π s, and (b) vs. π at various π s. Dark blue triangles in (a) correspond to measurements reported in Ref.[35] at 300 K. Colored symbols in (a, b) correspond to calculations reported in Ref. [40]and measurements in Ref. [36] at (a) π s and (b) π s represented by colors in the color-bars. (c) Thermal expansivity of akimotoite vs. π at various π s.4Figure 6. The π ππ βs of akimotoite vs. (a, b) π at various π s and (c, d) vs. π at various π s. Blue squares and triangles in (a, b) correspond to measurements in Ref. [41] and Ref. [35] at 300 K, purple squares in (c, d) correspond to measurements in Ref. [41] at 0 GPa.5Figure 7. (a, c) Elastic moduli and (c, d) acoustic velocities of akimotoite vs. (a, b) π and (c, d) π . Dark blue triangles in (a) correspond to measurements reported in Ref. [35] at 300 K. Colored symbols correspond to measurements reported in Ref. [36] at (a) π s and (b) π s represented by colors in the color-bars.6Figure 8. (a, b) Unit-cell volume (16 f.u.) of Mg-Pv (a) vs. π at various π s, and (b) vs. π at various π s. Dark blue pentagons and triangles in (a, b) correspond to measurements at 300 K reported in Refs. [43,45]. Colored squares, diamonds, pluses, and crosses correspond to high ππ GGA + MD calculations reported in Ref. [5] and measurements in Ref. [44,46]. Their π s (a) and π s (b) are represented by colors in color-bars. 7Figure 9. The π ππ βs of Mg-Pv vs. (aβc) π at various π s, and (dβf) vs. π at various π s. Dark blue (aβc) and purple (dβf) circles correspond to 300 K and 0 GPa measurements reported in Ref. [49]. Colored squares correspond to GGA + MD results reported in Ref. [5] . Their π s (aβc) and π s (dβf) are represented by colors in colorbars. 8Figure 10. (a, c) Elastic moduli and (c, d) acoustic velocities of Mg-Pv vs. (a, b) π and (c, d) π . Dark blue (aβc) and purple (dβf) circles in (a, b) correspond to 0 GPa and 300 K measurements reported in Ref. [49]; dark blue pentagons in (a) correspond to 300 K measurements [45]. Colored squares, triangles, crosses, and diamonds correspond to high ππ GGA + MD results reported in Ref. [5] and measurements reported by [44,46,50].Their π s (aβc) and π s (dβf) are represented by colors in color-bars. 9 Tables
Table 1. List of modules and command-line utilities in the cij distribution
Module Description cij.core
Core functionalities β’ calculator β The calculator that controls the workflow. β’ mode_gamma β Interpolate phonon frequencies and calculate mode GrΓΌneisen parameters. β’ phonon_contribution β Calculate π ππph . β’ full_modulus β Interpolate π ππst vs. π , and calculate π πππ and π πππ . β’ tasks β Handles the ordering of π ππph calculation. cij.util Methods used in the main program β’ voigt β Convert between Voigt ( π ππ ) and regular ( π ππππ ) notations of elastic coefficients. β’ units β Handle unit conversion. cij.io Input output functions and classes. cij.plot
Plotting functionalities. cij.cli
Command-line programs β’ cij run ( main.py ) β Perform a SAM-Cij calculation. β’ cij run-static ( static.py ) β Calculate static elastic properties. β’ cij extract ( extract.py ) β Extract calculation results for a specific π or π to a table. β’ cij extract-geotherm ( geotherm.py ) β Extract calculation results along geotherm ππ (given as input) to a table. β’ cij plot ( plot.py ) β Convert output data table to PNG plot. β’ cij modes ( modes.py ) β Plot phonon frequency interpolation results. β’ cij fill ( fill.py ) β Fill all the non-zero terms for elastic coefficients given the constraint of a crystal system. cij.data Data distributed with the program, e.g., the relationship between π ππ βs for different crystal systems, the naming scheme for output files, etc. cij.misc Miscellaneous functionalities that are not used in the main program, e.g., methods that facilitate the preparation of input files.0Table 2. Available keywords and their options in settings.yaml . Under " qha.settings " DELTA_P_SAMPLE number
Pressure-sampling interval, used for output, the default value is 1 GPa.
DELTA_P number
The interval between two nearest pressures on the grid, in GPa.
P_MIN number
The minimum pressure in GPa.
NTV integer
Number of volumes (or equivalently, pressures) on the grid.
T_MIN number
The minimum temperature in Kelvin.
DT number
The interval between two nearest temperatures on the grid in Kelvin.
NT integer
The number of temperatures on the grid.
Under " elast.settings " system string The crystal system used, one of: triclinic , monoclinic , hexagonal, trigonal7 , trigonal8 , orthorhombic , tetragonal6 , tetragonal7 , cubic . mode_gamma.interpolator string The method to interpolate phonon frequencies vs. volume, one of: spline , lsq_poly , lagrange , krogh , pchip , hermite , akima . mode_gamma.order integer The order of phonon frequencies spline interpolation.
Under " qha " and " elast " input string The location of the input files.Table 3. The output variable, their keyword in settings output section, default output unit,and their output file naming conventions.Property Keyword Unit Output file naming convention (i,j = 1β6; base = tp or t v ) Adiabatic elastic modulus cij_s GPa c{ij}s_{base}_gpa.txt
Isothermal elastic modulus cij_t
GPa c{ij}t_{base}_gpa.txt bm_V
GPa bm_V_{base}_gpa.txt
Reuss average of bulk modulus bm_R
GPa bm_R_{base}_gpa.txt
Voigt-Reuss-Hill average of bulk modulus bm_VRH
GPa bm_VRH_{base}_gpa.txt
Voigt average of shear modulus
G_V
GPa
G_V_{base}_gpa.txt
Reuss average of shear modulus
G_R
GPa
G_R_{base}_gpa.txt
Voigt-Reuss-Hill average of shear modulus
G_VRH
GPa
G_VRH_{base}_gpa.txt
Shear acoustic wave velocities v_s km/s v_s_{base}_km_s.txt
Compressive acoustic wave velocities v_p km/s v_p_{base}_km_s.txt
Pressure vs. volume p GPa v_tp_gpa.txt
Volume vs. pressure and temperature V Γ p_tv_ang3.txt Table 4. The structure of static elastic coefficients input data file ( elast.dat ) Structure of input data Description V N m cell
The calibration volume π for static elastic moduli interpolation; number of volumes included in this data file π ; total cell mass, π cell , in Dalton for calculation of acoustic wave velocities calculationsβ v c11 c22 c33 β¦ β Column names of the input data.2The output elastic moduli are named after this list.V c [V ] c [V ] β¦ The first volume and elastic moduli at this volume; the order corresponds to the column names specified above.V c [V ] c [V ] β¦ Similarly organized data for subsequent volumes.β¦V N c [V N ] c [V N ] β¦ β lattice_a lattice_b lattice_c β Column names for axial lengths tablea [V ] a [V ] a [V ] The axial length along three axes of Cartesian coordinates for the first volume.a [V ] a [V ] a [V ] The axial length for subsequent volumes.β¦a [V N ] a [V N ] a [V N ]3 References [1] D. Fan, Z. Mao, J. Yang, J.-F. Lin, Am. Mineral. 100 (2015) 2590β2601.[2] J.-F. Lin, Z. Mao, J. Yang, S. Fu, Nature 564 (2018) E18βE26.[3] B.B. Karki, Science 286 (1999) 1705β1707.[4] R.M. Wentzcovitch, B.B. Karki, M. Cococcioni, S. de Gironcoli, Phys. Rev. Lett. 92 (2004) 018501.[5] A.R. Oganov, J.P. Brodholt, G.D. Price, Nature 411 (2001) 934β937.[6] J. Zhao, J.M. Winey, Y.M. Gupta, Phys. Rev. B 75 (2007) 094105.[7] Z. Wu, R.M. Wentzcovitch, Phys. Rev. B 83 (2011) 184115.[8] K. KΓ‘das, L. Vitos, R. Ahuja, B. Johansson, J. KollΓ‘r, Phys. Rev. B 76 (2007) 235109.[9] C. Malica, A.D. Corso, J. Phys. Condens. Matter 32 (2020) 315902.[10] S.-L. Shang, H. Zhang, Y. Wang, Z.-K. Liu, J. Phys. Condens. Matter 22 (2010) 375403.[11] Y. Wang, J.J. Wang, H. Zhang, V.R. Manga, S.L. Shang, L.-Q. Chen, Z.-K. Liu, J. Phys. Condens. Matter 22 (2010) 225404.[12] W. Wang, Y. Xu, D. Sun, S. Ni, R. Wentzcovitch, Z. Wu, Nat. Commun. 11 (2020) 64.[13] G. Shukla, M. Cococcioni, R.M. Wentzcovitch, Geophys. Res. Lett. 43 (2016) 5661β5670.[14] G. Shukla, Z. Wu, H. Hsu, A. Floris, M. Cococcioni, R.M. Wentzcovitch, Geophys. Res. Lett. 42 (2015) 1741β1749.[15] G. Shukla, K. Sarkar, R.M. Wentzcovitch, J. Geophys. Res. Solid Earth 124 (2019) 2417β2427.[16] S. Hao, W. Wang, W. Qian, Z. Wu, Earth Planet. Sci. Lett. 528 (2019) 115830.[17] F. Zou, Z. Wu, W. Wang, R.M. Wentzcovitch, J. Geophys. Res. Solid Earth 123 (2018) 7629β7643.[18] M. Destefanis, C. Ravoux, A. Cossard, A. Erba, Minerals 9 (2019) 16.[19] M. Born, Math. Proc. Camb. Philos. Soc. 36 (1940) 160β172.[20] D.C. Wallace, Thermodynamics of Crystals, Wiley, New York, 1972.[21] T. Qin, Q. Zhang, R.M. Wentzcovitch, K. Umemoto, Comput. Phys. Commun. 237 (2019) 199β207.[22] R.M. Wentzcovitch, Z. Wu, P. Carrier, Rev. Mineral. Geochem. 71 (2010) 99β128.[23] G.F. Davies, J. Phys. Chem. Solids 35 (1974) 1513β1520.[24] P.R.C. da Silveira, L. Gunathilake, A. Holiday, D.A. Yuen, M.N. Valdez, R.M. Wentzcovitch, in: Proc. Conf. Extreme Sci. Eng. Discov. Environ. Gatew. Discov. - XSEDE 13, ACM Press, San Diego, California, 2013, p. 1.[25] M.J.P. Musgrave, Crystal Acoustics: Introduction to the Study of Elastic Waves and Vibrations in Crystals, Holden-day, 1970.[26] J.P. Perdew, A. Zunger, Phys. Rev. B 23 (1981) 5048β5079.[27] L. Sang, J.D. Bass, Phys. Earth Planet. Inter. 228 (2014) 75β79.[28] D.G. Isaak, I. Ohno, P.C. Lee, Phys. Chem. Miner. 32 (2006) 691β699.[29] B. Li, D.R. Neuville, Phys. Earth Planet. Inter. 183 (2010) 398β403. [30] A.M. Walker, Phys. Earth Planet. Inter. 192β193 (2012) 81β89.[31] H. Sawamoto, High-Press. Res. Miner. Phys. 209 (1987) 219.[32] C.R.S. Da Silva, B.B. Karki, L. Stixrude, R.M. Wentzcovitch, Geophys. Res. Lett. 26 (1999) 943β946.[33] J. Wookey, J.-M. Kendall, G. Barruol, Nature 415 (2002) 777β780.[34] A. Nowacki, J.-M. Kendall, J. Wookey, A. Pemberton, Geochem. Geophys. Geosystems 16 (2015) 764β784.[35] N.C. Siersch, D. Frost, The Effect of Fe and al on the Elasticity of Akimotoite, UniversitΓ€t Bayreuth, 2019.[36] C. Zhou, Phys. Earth Planet. Inter. (2014) 9.[37] M.L. Marcondes, R.M. Wentzcovitch, J. Appl. Phys. 117 (2015) 215902.[38] Z. Wu, R.M. Wentzcovitch, K. Umemoto, B. Li, K. Hirose, J.-C. Zheng, J. Geophys. Res. 113 (2008) B06204.[39] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865β3868.[40] L. Li, D.J. Weidner, J. Brodholt, D. AlfΓ¨, G.D. Price, Phys. Earth Planet. Inter. 173 (2009) 115β120.[41] D.J. Weidner, E. Ito, Phys. Earth Planet. Inter. 40 (1985) 65β70.[42] Z. Wu, R.M. Wentzcovitch, Earth Planet. Sci. Lett. 457 (2017) 359β365.[43] D. Andrault, N. Bolfan-Casanova, N. Guignot, Earth Planet. Sci. Lett. 193 (2001) 501β508.[44] B. Li, J. Zhang, Phys. Earth Planet. Inter. 151 (2005) 143β154.[45] M. Murakami, S. Sinogeikin, H. Hellwig, J. Bass, J. Li, Earth Planet. Sci. Lett. 256 (2007) 47β54.[46] J. Chantel, D.J. Frost, C.A. McCammon, Z. Jing, Y. Wang, Geophys. Res. Lett. 39 (2012) n/a-n/a.[47] G. Fiquet, A. Dewaele, D. Andrault, M. Kunz, T. Le Bihan, Geophys. Res. Lett. 27 (2000) 21β24.[48] P. Carrier, R. Wentzcovitch, J. Tsuchiya, Phys. Rev. B 76 (2007) 064116.[49] A. Yeganeh-Haeri, Phys. Earth Planet. Inter. 87 (1994) 111β121.[50] M. Murakami, Y. Ohishi, N. Hirao, K. Hirose, Nature 485 (2012) 90β94.[51] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, Davide Ceresoli, G.L. Chiarotti, M. Cococcioni, I. Dabo, A.D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, Anton Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, Stefano Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A.P. Seitsonen, A. Smogunov, P. Umari, R.M. Wentzcovitch, J. Phys. Condens. Matter 21 (2009) 395502.[52] B.B. Karki, R.M. Wentzcovitch, S. de Gironcoli, S. Baroni, Phys. Rev. B 61 (2000) 8793β8800.[53] K. Umemoto, R.M. Wentzcovitch, Y.G. Yu, R. Requist, Earth Planet. Sci. Lett. 276 (2008) 198β206.[54] N. Troullier, J.L. Martins, Phys. Rev. B 43 (1991) 1993β2006.[55] D. Vanderbilt, Phys. Rev. B 41 (1990) 7892β7895.5