Capturing Anharmonicity in a Lattice Thermal Conductivity Model for

Article pubs.acs.org/cm

Capturing Anharmonicity in a Lattice Thermal Conductivity Model for High-Throughput Predictions Samuel A. Miller,† Prashun Gorai,‡,§ Brenden R. Ortiz,§ Anuj Goyal,‡,§ Duanfeng Gao,∥ Scott A. Barnett,† Thomas O. Mason,† G. Jeffrey Snyder,† Qin Lv,∥ Vladan Stevanović,‡,§ and Eric S. Toberer*,‡,§ †

Northwestern University, Evanston, Illinois 60208, United States National Renewable Energy Laboratory, Golden, Colorado 80401, United States § Colorado School of Mines, Golden, Colorado 80401, United States ∥ University of Colorado Boulder, Boulder, Colorado 80309, United States ‡

S Supporting Information *

ABSTRACT: High-throughput, low-cost, and accurate predictions of thermal properties of new materials would be beneficial in fields ranging from thermal barrier coatings and thermoelectrics to integrated circuits. To date, computational efforts for predicting lattice thermal conductivity (κL) have been hampered by the complexity associated with computing multiple phonon interactions. In this work, we develop and validate a semiempirical model for κL by fitting density functional theory calculations to experimental data. Experimental values for κL come from new measurements on SrIn2O4, Ba2SnO4, Cu2ZnSiTe4, MoTe2, Ba3In2O6, Cu3TaTe4, SnO, and InI as well as 55 compounds from across the published literature. To capture the anharmonicity in phonon interactions, we incorporate a structural parameter that allows the model to predict κL within a factor of 1.5 of the experimental value across 4 orders of magnitude in κL values and over a diverse chemical and structural phase space, with accuracy similar to or better than that of computationally more expensive models.

■

als.11 However, to make high-throughput predictions, it is necessary to make various assumptions. Heat is transported through multiple-phonon scattering processes, which are commonly determined by solving the Boltzmann transport equation and calculating group velocities, phonon frequencies, and both harmonic and anharmonic interatomic force constants. This is especially expensive when considering three phonon scattering processes, which requires calculation of third-order anharmonic interatomic force constants.12,13 Several studies have attempted to predict κL for large sets of materials. Seko et al.14 used expensive anharmonic lattice-dynamics calculations to identify low-lattice thermal conductivity materials; Carrete et al.15 used ab initio calculations and machine-learning algorithms to predict κL values for halfHeusler compounds, while Madsen et al.16 and Toher et al.17 employed the quasiharmonic Debye approximation to compute thermal properties, including κL. We have previously developed a semiempirical model for determining κL, which combines first-principles calculations and experimental measurements and

INTRODUCTION Thermal properties of materials, namely, thermal conductivity and heat capacity, are important for applications ranging from thermal barrier coatings to integrated circuits and light-emitting diodes.1−3 This is especially true in thermoelectric materials where the thermoelectric efficiency is determined by the thermal conductivity and thus a low lattice thermal conductivity (κL) is desirable.4,5 However, finding low-κL materials remains a challenge from an experimental standpoint as synthesis procedures are costly and time-consuming.6 To drive these technologies forward, efficient and accelerated searches with accurate predictions of κL of new functional materials are necessary.7 The ability to accurately predict thermal properties across a wide range of materials with simple calculations would be beneficial for decreasing the time to realization of novel materials. Intuitively, one would expect that κL will decrease with an increasing atomic mass (heavier atoms) and atoms in the primitive unit cell (complex structure).8,9 Figure 1 shows that this is generally the case, though there still exists an order of magnitude variation in κL for compounds with a similar atomic mass and a similar number of atoms. This demonstrates that intuition is limited and robust calculations require more nuanced theory. Recent advances in high-throughput computations have accelerated experimental searches of novel functional materi© 2017 American Chemical Society

Special Issue: Computational Design of Functional Materials Received: September 30, 2016 Revised: January 6, 2017 Published: January 6, 2017 2494

DOI: 10.1021/acs.chemmater.6b04179 Chem. Mater. 2017, 29, 2494−2501

Article

Chemistry of Materials

first-principles calculations such as density functional theory (DFT). In this semiempirical treatment of κL, the Grüneisen parameter (γ) is treated as a constant and incorporated into A1. By assuming it is the same for all materials, the computational costs associated with explicit calculations of γ are eliminated. Computational Details. The DFT calculations were performed as previously described in refs 18 and 20 on 26 additional compounds and can be found in the Supporting Information. Accurate computational assessment of γ is expensive, making it intractable for high-throughput predictions of complex materials.17 We propose to use local coordination as an estimate of γ. The local coordination of an atomic site in a crystal is established by estimating the number of neighboring atoms directly bonded to it. Neighbors directly bonded to an atomic site are determined on the basis of the shortest neighbor distance. For each atomic site, all neighbors within 10% of the nearest neighbor length are counted toward the site coordination. An average coordination number is then computed for a given crystal structure by summing the coordination number for all the atomic sites and dividing by the total number of atoms. Synthesis and Characterization. To augment the data set used for fitting the semiempirical model, we synthesized Ba3In2O6, Ba2SnO4, Cu3TaTe4, Cu2ZnSiTe4, InI, MoTe2, SnO, and SrIn2O4. These compounds were chosen to represent a diverse chemical and structural space and because of their predicted κL values of 4 are overestimated. With the assumption that bonding and local coordination have the largest effect on the Grüneisen parameter, γ can be approximated on the basis of the average coordination number. A functional form for this relationship is chosen such that the estimated γ is lower at small coordination and higher for large coordination as is observed in experimental data: γmodeled = γ0[1 − e−a(CN − CN0)]

κL,tot = A1

y Mv 3k b ⎛ π ⎞1/3 vs ⎛ ̅ s 1 ⎞ ⎜ ⎟ + z ⎜1 − 2/3 ⎟ 2 z x ⎝ ⎠ ⎝ 2 6 V Tγ V n n ⎠

(4)

2

To avoid overfitting of data, the adjusted R was used. A single term in the model was allowed to vary, each one at a time; it was found that the free parameter leading to the lowest AFD was A1. Radj.2 was determined, and a second free parameter was added to the model, allowing both A1 and the second parameter to vary. This process was repeated until the fifth free parameter

(3)

This equation was fit to all experimental γ measurements in this study using the average coordination as previously defined. The 2497


Article


exponents on each property, neither considers the κL,opt contribution to the total κL. Because of the use of three different models each using a distinct way of computing the Grüneisen parameter, there is variation in the modeled γ for each experimental value, as seen in Figure 8. Though this work appears to underestimate γ, the

led to a decrease in Radj.2, indicating that four fitted variables, A1, x, y, and z, lead to the best model without overfitting the data (see the Supporting Information). The resulting predictions using this model in Figure 7 show good agreement between experimental and predicted κL values. Note that the predicted κL values for this figure are generated using the “leave one out” method. The new model with four fitting variables was fit to 62 data points and then used to predict the κL of the 63rd compound (i.e., the model is fit to all compounds in the study except GaN, and then the κL of GaN is predicted using this model) and is iterated over the entire data set. Furthermore, 4-fold cross validation was performed to determine the validity of the semiempirical model in not only describing but also predicting κL. The results are summarized in Table 1. Within this data set, both the leave one out method and cross validation indicate that the model is predictive rather than just descriptive. Table 1. Summary of the AFD and Each of the Fitting Variables in the Cross Validation Study for Each Test Set as Well as the Average and the Results of Fitting the Entire Data Set test set 1 test set 2 test set 3 test set 4 average reported model

AFD

A1

x

y

z

1.63 1.56 1.39 1.41 1.49 1.48

0.00217 0.00202 0.00196 0.00230 0.00225 0.00230

1.04 1.06 1.07 0.994 1.04 1.03

4.59 4.64 4.61 4.36 4.55 4.51

0.356 0.357 0.329 0.340 0.345 0.334

Figure 8. Comparison of γ approximations between the semiempirical model and previous high-throughput efforts. Though simple, the estimation of γ using the coordination number is no less accurate than other high-throughput methods.

fit used here produces the minimum possible average factor difference due to numerous γexp measurements at 0.75 and 1.5 for tetrahedral and octahedral compounds, respectively. The mode-averaged γ predictions of both this work and ref 16 have strong correlations with experiment (Table 3). For both the Pearson and Spearman rank correlation coefficients, a value of 1 shows ideal positive correlation with the experimental values. However, these correlation coefficients are only a measure of linear or monotonic relationships between two values and do not account for the offset between experimental and predicted properties. As such, the average factor difference is needed to fully grasp the accuracy of the various models. The inexpensive Mie−Grüneisen method employed by ref 17 performs poorly on all metrics, systematically overestimating γ and having a weak correlation with experimental data. The high correlation coefficients and low AFD demonstrate that the simple method presented here is at least as accurate in estimating γ as more computationally expensive high-throughput methods. Though approximations for γ are needed, whether using the quasiharmonic Debye or semiempirical model, the property of interest is κL. To further verify the predictive power of this semiempirical model and compare it to other models, the same four randomized test sets used for cross validation were compared with other models. Because the training set for the semiempirical model does not include the test set, this allows for comparison of the AFD between the semiempirical model and those of refs 16 and 17 within each of the test sets. The results for each of the four test sets are summarized in Table 2, which shows that our model performs similar to the model of ref 16 and better than the model of ref 17 within each test set and on average. Additionally, there is not significant variation in the ability of the semiempirical model to predict κL in any given randomized set taken from the complete data set. Also note that the κL prediction accuracy is similar for both simple and complex compounds. As summarized in Table 2, the AFD

The model presented here accurately predicts lattice thermal conductivity across a diverse range of compounds over 4 orders of magnitude in κL. The resulting model has a factor difference of 1.48, a significant improvement over the previous work. While the expansion and diversification of the data set provide justification for increased fitting parameters according to the adjusted R2 analysis, the majority of the improvement is driven by the treatment of γ as seen in Figure 4. Here we approximate γ using the average coordination number (eq 3) and note that the semiempirical model simplifies the speed of sound as vs ≃ (B/d)1/2. Comparison of a Refined Model with Previous Efforts. Accurate, low-cost, and high-throughput predictions of κL are useful for a variety of applications. To highlight the advantages of this semiempirical model, it is compared to previous highthroughput methods used in the literature. Both Madsen et al.16 and Toher et al.17 employ the quasiharmonic Debye approximation and calculate κL according to the equation proposed by Slack.8 In ref 16, the authors perform full Grüneisen parameter calculations followed by mode averaging, resulting in fairly accurate but expensive calculations that limit the data set to only rocksalt, zincblende, and diamond structures (two atoms per unit cell). In ref 17, an alternative approach is taken, screening a much larger set of compounds to rank different classes by order of magnitude. Reference 17 uses the Mie−Grüneisen equation to calculate γ, requiring just the volume, pressure, number of atoms, and Debye temperature. The result is a much less computationally expensive method that allows a much larger number of calculations to be performed on more complex cells but produces less accurate predictions. Though refs 16 and 17 both use the classical 2498


Article

Chemistry of Materials Table 2. Comparison of the Average Factor Difference Using 4-fold Cross Validation Test Sets and the Leave One Out (LOO) Methoda

test set 1 test set 2 test set 3 test set 4 average LOO method

data set 1, 24 compounds



ref 16

SE

ref 17

SE

SE

2.03 1.57 1.81 1.53 1.73 −

1.68 1.52 1.18 1.50 1.47 1.49

3.19 3.07 4.08 3.69 3.51 −

1.55 1.50 1.42 1.42 1.43 1.45

1.71 1.61 1.23 1.23 1.51 1.52

a Data set 1 includes the compounds in ref 16 and the present semiempirical (SE) model, data set 2 those in ref 17 and the SE model, and data set 3 those compounds used in this work but neither ref 16 nor 17.

values for more complex compounds not considered in other models are 1.51 and 1.52 for the average of the 4-fold cross validation and the leave one out method, respectively. These are both similar to the AFD of 1.48 for the data set as a whole, demonstrating the versatility of this semiempirical model. We have demonstrated through 4-fold cross validation that the fitting parameters do not significantly change between the test sets (Table 1) and that the predictions for these tests sets compare favorably with those of other models (Table 2). Therefore, κL predictions are made for a single compound using the leave one out method and compared on a case by case basis with those of the other models. The resulting predictions are shown in Figure 9 by taking the ratio of predicted to experimental lattice thermal conductivity. Thus, predictions of κL equal to the experimental value have a value of 1.0, whereas under- and overestimations by a factor of 2 are shown as 0.5 and 2.0, respectively. Figure 9a shows the compounds considered in all three models while Figure 9b shows those more complex structures only in ref 17 and this work. The compounds treated only in this semiempirical model are shown in Figure S5. In all cases, compounds are shown in order of increasing experimental κL. As seen in Figure 9a and Table 3, the semiempirical model developed herein performs equally well at predicting κL. On average, the factor difference between predicted and experimental κL values for the model presented here and that used by ref 16 is half that of ref 17. Similarly, the Pearson and Spearman rank correlation coefficients for these two models are comparable and closer to the ideal value of 1 than that for the predictions from ref 17. The κL predictions from ref 16 are accurate for the simple rocksalt, zincblende, and diamond structures, but their method is much more computationally expensive (a few orders of magnitude for simple structures and even more costly for complex structures). The method of ref 17 is useful for screening large data sets but overestimates γ and underestimates the lattice thermal conductivity for compounds with experimental κL values of ≥10 W m−1 K−1. Limitations. Given the semiempirical nature of the model, the predictive power is limited by the data to which the model is fit. We have included covalent and ionic materials with chemistries ranging from simple elements to complex quaternaries and structural diversity spanning simple crystals (d-C) to complex Zintl phases, cage-like structures to highly anisotropic quasi-two-dimensional structures. Our model can be applied to bulk materials but not micro- or nanostructured

Figure 9. Comparison of predicted κL values with experimental values finds the method herein combines accuracy with low computational cost. (a) Values of rocksalt, zincblende, and diamond structured materials calculated in all models and (b) values of more complicated structures handled by only the model presented here and ref 17. Horizontal dashed lines show over- and underestimation by a factor of 2 from experimental κL.

Table 3. Summary of the Comparison between This Work and Other Previously Described High-Throughput Models, for γ and κL γ factor difference γ Pearson correlation γ Spearman correlation κL factor difference κL Pearson correlation κL Spearman correlation

this work

ref 16

ref 17

1.21 0.815 0.799 1.48 0.983 0.953

1.54 0.834 0.654 1.71 0.983 0.957

2.17 0.214 0.0917 3.36 0.969 0.807

materials. The model is meant to handle stoichiometric materials, so while intermetallic compounds such as halfHeuslers can be handled by the model, the κL values of solid solution alloys such as Si1−xGex cannot be accurately predicted. Small amounts of substituted elements, whether as an alloy or dopant, are not expected to drastically change the κL but could affect descriptors used here such as the number of atoms per unit cell. As such, application of our model should be limited to 2499



■

well-defined crystal structures found in databases such as the ICSD, allowing prediction of κL of the “parent” material. From these predictions of the stoichiometric compound, alloy scattering models such as those proposed by Klemens25 or Yang et al.26 can be used for substituted materials.

CONCLUSIONS Here we “close the loop” between theory and experiment, yielding a predictive model of thermal conductivity that has an average factor difference of

Capturing Anharmonicity in a Lattice Thermal Conductivity Model for

Recommend Documents