Charting the Lattice Thermal Conductivities of I-III-VI2Chalcopyrite Semiconductors

Charting the Lattice Thermal Conductivities of I−III−VI2Chalcopyrite Semiconductors Jose J. Plata, Victor Posligua, Antonio M. Márquez, Javier Fernandez Sanz, and Ricardo Grau-Crespo* Cite This: Chem. Mater. 2022, 34, 2833−2841 Read Online ACCESS Metrics & More Article Recommendations * sıSupporting Information ABSTRACT: Chalcopyrite-structured semiconductors have promising potential as low-cost thermoelectric materials, but their thermoelectric figures of merit must be increased for practical applications. Understanding their thermal properties is important for engineering their thermal conductivities and achieving better thermoelectric behavior. We present here a theoretical investigation of the lattice thermal conductivities of 20 chalcopyrite semiconductors with an ABX2composition (I−III−VI2) (A = Cu or Ag; B = Al, Ga, In, or Tl; X = S, Se, or Te). To afford accurate predictions across this large family of compounds, we solve the Boltzmann transport equation with force constants derived from density functional theory calculations and machine learning-based regression algorithms, reducing by between 1 and 2 orders of magnitude the computational cost with respect to conventional approaches of the same accuracy. The results are in good agreement with available experimental data and allow us to rationalize the role of chemical composition, temperature, and nanostructuring in the thermal conductivities across this important family of semiconductors. 1. INTRODUCTION The chalcopyrite structure, with the formula unit ABX2and tetrahedrally coordinated A and B cations and X anions (Figure 1), is a versatile structure for semiconducting materials. It can be obtained by inexpensive synthesis in a wide variety of chemical compositions, because of the ability of this structure to accommodate a range of cations and anions in either a I− III−VI2(A = Cu or Ag; B = Al, Ga, In, or Tl; X = S, Se, or Te) or a II−IV−V2(A = Zn or Cd; B = Si, Ge, or Sn; X = P, As, or Sb) stoichiometry. The former can be regarded as the ternary analogues of the II−VI binary semiconductors (like ZnS), and the latter as the ternary analogues of the III−V binary semiconductors (like GaP). The ternary chalcopyrite-struc- tured semiconductors first attracted interest in the 1960s in the hope of thermoelectric applications, because they are cheap, tend to have low lattice thermal conductivities, and can be obtained as both p- and n-type semiconductors. Despite some of the materials exhibiting good carrier mobility, sufficiently high thermoelectric figures of merit could not be found for any of the pure compounds. 1 In recent years, interest in this family of materials has returned, primarily driven by their photovoltaic applications. Chalcopyrite-based thin films are important in the solar cell sector, with copper indium gallium selenide (CIGS) compounds, with a CuInGaSe2composition, showing some of the highest photovoltaic efficiencies among all thin-film technologies. 2−4 Interest in thermoelectric applications has also come back, as researchers have learned to improve the dopability, synthesize more complex solid solutions, and decrease the thermal conductivity with nanostructuring strategies. 5−10 For example, high thermoelectric figures of merit (zT) have been found for p-type Ag0.95GaTe2(zT = 0.8 at 850 K), 10 n-type Ag0.9Cd0.1InSe2(zT = 1.1 at 900 K), 9 and p-type CuGaTe2(zT = 1.4 at 950 K). 8 Computer simulations, based on density functional theory (DFT), are improving our understanding of the thermoelectric behavior of chalcopyrite compounds 11−13 and are even used to predict new chalcopyrite compositions of interest for thermoelectric applications. 14 The thermal conductivity (κ) is important for thermoelectric applications because the thermoelectric figure of merit zT is inversely proportional to κ. Therefore, a systematic understanding of phonon transport in chalcopyrite semiconductors, across compositions and temperatures, and of the effect of nanostructuring on the thermal conductivity, is essential for the design of better device components based on these materials. We present here a computational investigation of the thermal conductivity of chalcopyrite semiconductors. We consider 20 experimentally reported compositions of I−III−VI2chalcopyr- Received: February 1, 2022 Revised: February 17, 2022 Published: March 1, 2022 Articlepubs.acs.org/cm © 2022 American Chemical Society 2833 https://doi.org/10.1021/acs.chemmater.2c00336 Chem. Mater. 2022, 34, 2833−2841 Downloaded via UNIV DE SEVILLA on May 26, 2023 at 12:53:31 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles. ite-structured nonmagnetic semiconductors, as listed in ref 1. They correspond to all ABX2compositions with A = Cu or Ag; B = Al, Ga, or In; and X = S, Se, or Te, and two compositions with B = Tl (CuTlS2and CuTlSe2). The accurate prediction of the lattice thermal conductivity, which makes the dominant contribution to κin semiconductors and insulators, is computationally challenging. One of the most accurate methods available is based on Boltzmann’s transport equation (BTE) for phonons, 15 which requires the calculation of the second-order and, at least, thirdorder interatomic force constants (IFCs). Traditionally, these IFCs are predicted by calculating atomic forces in supercells in which one or two atoms are systematically displaced from their equilibrium positions, using DFT. 16,17 The number of DFT calculations required to obtain third-order IFCs is usually very high, and thus, this step represents the bottleneck in the prediction of κfrom first principles. 16−20 Alternatively, empirical expressions can be applied to efficiently estimate κ, without using IFCs, 21−26 but these methods rely on experimental data that are not always available and suffer from accuracy issues that make them unsuitable when quantitative information is required, such as for the calculation of the thermoelectric figure of merit. New algorithms have recently become available to accelerate the calculation of IFCs, using techniques such as compressive sensing and machine learning (ML) to extract the necessary information from a much smaller number of DFT calculations. 27−29 These new methods open the door to the accurate calculation of κ, not only for more complex materials but also for large families of materials, as described below for I−III−VI2chalcopyrites. Because of the sensitivity of κto synthetic conditions, which control grain size and defect chemistry, any trends in thermal conductivities extracted from experiments under different conditions should be taken very cautiously. In contrast, computer modeling allows a direct comparison of intrinsic thermal conductivity behavior across compositions and temperatures. 2. METHODOLOGY DFT Calculations. Tetragonal conventional cells (16 atoms) were fully relaxed using the VASP code 30 with the projector-augmented wave method, 31 where the number of valence electrons for each atom was selected following standards proposed by Calderon et al. 32 Energies and forces were obtained using the generalized gradient approximation (GGA) functional proposed by Perdew, Burke, and Ernzerhof (PBE), 33 combined with Grimme’s D3 van der Waals corrections 34 and a plane-wave basis set with a cutoffof 500 eV. Structures were relaxed until the forces acting over each atom were <10−7eV/Å. An additional support grid for the evaluation of the augmentation charges was included to reduce noise in the forces. The forces for obtaining the IFCs were calculated using a 4 ×4×2 supercell, as recently used by Park et al. for CuFeS2. 35 Reciprocal space integrations were performed only at the Γpoint of the supercell, but increasing the grid density to a Γ-centered 2 ×2×2 mesh did not have any significant effect on the results. Force Constant Prediction and Machine Learning. We used a machine learning-based approach implemented in the HiPhive package 28 to extract second-, third-, and fourth-order force constants within optimized cutoffdistances. Although the fourth-order force constants do not enter into our BTE model, their inclusion leads to a better regression for the force constant potential (FCP) model. The force constants are determined from multilinear regression to the DFT forces using the recursive feature elimination (RFE) algorithm. 36 A more detailed description of the training of the model can be found in the Supporting Information. The convergence of the parameters (number of distorted structures and cutoffdistances) involved in the FCP model was tested by assessing the variation of the force errors, the phonon frequencies, and the lattice thermal conductivity itself. We used CuGaTe2as a representative case for these tests because of the availability of experimental data for that composition. The comparison with κvalues obtained from the “full-DFT”approach, where the IFCs are determined directly from the displacement of one or two specific atoms in the supercell, thus requiring many more DFT calculations, is also presented in the case of CuGeTe2. For this composition, cutoff distances of 11, 6.2, and 4 Å for the second-, third-, and fourth-order force constants were found to be sufficient for convergence of the ML-based method, and these cutoffs were extrapolated to other compositions based on coordination shells (rather than distances) for the sake of consistency. We developed a wrapper code for the hiPhive program that automates the distorted supercell creation, force calculation using VASP, and the construction of the ML FCPs, which is available for public use. 37 Boltzmann’s Transport Equation Solution. Once the FCP model is built, lattice thermal conductivities were obtained by solving the Boltzmann transport equation using the ShengBTE code. 17 We used the full iterative procedure to go beyond the relaxation time approximation. Scattering times were computed, including isotopic and three-phonon scattering. A Gaussian smearing of 0.1 eV and a dense mesh of 20 ×20 ×10 qpoints were used in all of the calculations, balancing the memory demand and the convergence of κ with the number of qpoints. The effects of including non-analytical contributions (NACs) on κwere tested in two compounds (CuGaS2 and AgGaS2), and only small changes (<2.5%) were found; therefore, the results reported below do not include NACs, to avoid the extra cost of computing Born effective charges. The scalar values reported correspond to one-third of the trace of the thermal conductivity tensor; small anisotropic effects will be discussed below. In what follows, we refer to the calculated lattice thermal conductivity as κ, which we compare with the experimental total (lattice + electronic) values, because the electronic contribution can be expected to be very small. 3. RESULTS AND DISCUSSION Obtaining accurate cell parameters is an important condition for computing accurate phonon properties. The lattice parameters of all compounds are well reproduced by our PBE-D3 calculations, as shown in Figure 2. The mean absolute Figure 1. Tetragonal unit cell of the chalcopyrite structure with an ABX2composition. Chemistry of Materials pubs.acs.org/cm Article https://doi.org/10.1021/acs.chemmater.2c00336 Chem. Mater. 2022, 34, 2833−2841 2834 errors are 0.025 Å (average 0.42% deviation) for aand 0.170 Å (1.5%) for c, which are well below the differences in cell parameters across compositions, leading to the good correlations observed in Figure 2. The inclusion of dispersion corrections improves the prediction of cell parameters in comparison with the uncorrected PBE. For example, for CuGaTe2, the PBE functional overestimates the aand clattice parameters by 1.51% and 2.00%, respectively, with respect to the experiment, whereas for PBE-D3, the discrepancies are −0.03% and 0.65%, respectively. For comparison, using the more expensive meta-GGA functional SCAN for this compound leads to discrepancies of −0.15% and 0.23%, respectively. The overestimation of cell volume by PBE leads to underestimation of phonon frequencies. By improving the description of the dispersion interactions and correcting the cell volume overestimation, the PBE-D3 provides an excellent balance between the computational cost and accuracy for the calculation of phonons (even when it does not correct for the limitations of the GGA in the prediction of the electronic structure). This is consistent with the conclusions of the investigation of phonon properties in semiconductors by Skelton et al., 38 who showed that standard GGA functionals overestimate cell volumes and predict overly soft phonons, but GGA-based functionals that correct the cell volumes (PBEsol, or D2/D3-corrected PBE) tend to give a much better description of the phonons and of the finite-temperature behavior. Within this family of compositions, the nature of the anion has the strongest effect on the cell parameters; tellurides have the largest and sulfides the smallest cell parameters. In agreement with the experiment, the nature of the A cation affects the aand cparameters differently, with the c/aratio being significantly lower for Ag-based (average c/a= 1.86) than for Cu-based (average c/a= 1.98) chalcopyrites. This has been observed elsewhere and related to the fact that the Ag− X−Ag bonds are softer than the Cu−X−Cu bonds. 3 The deviation from c/a= 2 characterizes the structural anisotropy in the system. Therefore, we can expect that Ag-based chalcopyrites will have slightly more anisotropic properties than Cu-based ones, which will be confirmed in terms of thermal conductivities below. Figure 3 illustrates that the force constant potentials obtained from the ML algorithms in HiPhive give results very similar to those obtained with the more computationally expensive “full-DFT”method, where force constants are determined from DFT calculations for all symmetrically distinct displacements of individual atoms. Using CuGaTe2 as an example, Figure 3a gives the dispersion curves obtained by the two methods, which are practically indistinguishable Figure 2. Comparison of DFT-calculated vs experimental cell parameters. 1 The solid line in each plot represents perfect agreement, and the green shaded area represents deviations of ≤2% from experiment in either direction. The horizontal error bars represent the range of values reported experimentally. Figure 3. (a) Comparison of CuGaTe2dispersion curves obtained using the machine-learned force constant potential (ML) via HiPhive with those from the full-DFT method via Phonopy. (b) Convergence of κwith the number of distorted structures used to obtain force constants (with a low third-order cutoff), in comparison with the full- DFT limit. The shaded area represents deviations of ≤10% in either direction. (c) Convergence of κwith the number of structures but for a higher third-order cutoff. In panels b and c, the points in red represent the evolution of the root-mean-square error (RMSE) in fitting the forces. Chemistry of Materials pubs.acs.org/cm Article https://doi.org/10.1021/acs.chemmater.2c00336 Chem. Mater. 2022, 34, 2833−2841 2835 from each other, demonstrating the equivalence between the two sets of second-order force constants. The third-order force constants are also very close. The calculated lattice thermal conductivity becomes closer to the full-DFT result when the number of distorted structures used in the ML fitting of the forces increases, and they get within 10% of each other when 16−20 structures are used (Figure 3b). This comparison was made using a relatively short cutoff rthird cut of 4.7 Å for the third-order force constants, in such a way that the full-DFT result could also be obtained. The full-DFT calculation with this cutoffrequired the evaluation of 348 distorted structures, i.e., ∼20 times more DFT calculations than when using the ML-based method. Because the DFT step represents the bottleneck in the calculation of the lattice thermal conductivities, the ML-based approach is almost 20 times faster (or more, depending on the cutoff) than the DFT- based approach. This huge saving of computing effort allows the computation of lattice thermal conductivities for the whole family of compounds in a high-throughput fashion. Furthermore, using the ML-based method, we can now afford to increase the third-order cutoff, which leads to a much better agreement between the fitted forces and the DFT forces (provided that enough data for training are available; if a very small number of distorted structures is used, a shorter cutoff leads to a better result). For example, for CuGaTe2, the rootmean-square error (RMSE) of the fitted forces went below 0.01 eV/Å with an rthird cut of 6.2 Å when four or more distorted structures were used (Figure 3c). Increasing the third-order cutofffrom 4.7 to 6.2 Å also affects the calculated lattice thermal conductivity, decreasing its value from 12.3 to 11.8 W m−1K−1when using forces from 18 structures for fitting, bringing the result slightly closer to the experimental value (10.7 W m−1K−1) measured by Bodnar et al. for CuGaTe2 single crystals. 39 The possibility of using larger cutoffs constitutes another advantage of the ML-based method over the full-DFT approach, which is necessarily constrained to small cutoffs due to the huge computational cost involved (e.g., >600 distorted structures would be needed for CuGaTe2when rthird cut = 6.2 Å). We have given here the cutoffdistances used for CuGaTe2, but it should be noted that the final cutoffs were kept consistent across different compositions by defining them by the coordination shells rather than by distance. The calculated average lattice thermal conductivities for all of the chalcopyrite compositions are summarized in Table 1 at 300 and 700 K, and the full dependence with temperature is plotted in Figure S1. We will first focus on discussing the room-temperature bulk values before examining the effect of temperature and/or nanostructuring. The most obvious trend is that the lattice thermal conductivity of Ag-based chalcopyrites is much lower than that of Cu-based ones. The average of the nine compositions corresponding to B = Al, Ga, or In, and X = S, Se, or Te, is ∼6.5 times higher for Cu-based (9.8 W m−1K−1) than for Ag-based (1.5 W m−1K−1) compositions. To examine the origin of this distinction, in panels a and b of Figure 4, we have plotted the phonon density of states corresponding to CuGaTe2and AgGaTe2, including the contributions from the A metal atoms. Clearly, Ag contributes modes at low frequencies, <1 THz, which does not happen in the Cu compounds. While the group velocity distribution is similar for both compounds (Figure 4c), the low-lying optical modes introduced by Ag lead to higher scattering rates (Wanh) at those frequencies, which dominate the behavior of the thermal conductivity (Figure 4d). On the contrary, the trend with the nature of the chalcogen X atom is not monotonous down the group. The thermal conductivity of selenides is lower than those of the corresponding sulfides and tellurides. This is illustrated in Figure 5a, by showing that the lattice thermal conductivities of CuBX2and AgBX2(averaged over B = Al, Ga, or In) reach minimum values when X = Se. This behavior is reminiscent of the nonmonotonous variation of the bandgap of lead chalcogenides (the bandgap of PbSe is lower than those of both PbS and PbTe), which results from the delicate balance of several electronic factors. 40 The three Ag-based selenide chalcopyrites (AgAlSe2, AgGaSe2, and AgInSe2) are then predicted to have remarkably low lattice thermal conductivities, all <1 W m−1K−1. The very low thermal conductivity of AgInSe2has been recently observed experimentally and rationalized in terms of an “avoided crossing”feature and low-lying optical modes in the phonon dispersion. 41 The effect of the B3+ cation is less regular. For example, making B = Ga leads to the lowest κamong the Ag-based selenides AgBSe2(0.77 W m−1K−1at 300 K, which is also the lowest value obtained in this study) but also leads to the composition with the highest κamong the Ag-based sulfides, AgBS2. In the Cu-based compounds, B = Ga always leads to the highest κacross the series of sulfides, selenides, and tellurides. This does not mean that κis insensitive to the nature of the B3+ cation, just that trends cannot be generalized as in the case of the ions in the A and X sites. In fact, we have also calculated two Cu-based compounds with B = Tl (CuTlS2and CuTlSe2; other combinations including Tl were not considered as they do not seem to be stable, or there is very little experimental information about them), and the thermal conductivities are much lower than those of the CuBS2and Table 1. Calculated Lattice Thermal Conductivities (κ) for 20 Chalcopyrite Compositions at 300 and 700 K a 300 K 700 K κ(W m−1K−1)L0.5 (nm) κ(W m−1K−1)L0.5 (nm) CuAlS28.27 335 3.39 159 CuAlSe27.25 159 3.06 63 CuAlTe210.2 192 4.34 76 CuGaS215.6 278 6.33 110 CuGaSe29.16 192 3.86 76 CuGaTe211.8 192 5.03 91 CuInS212.3 175 5.05 76 CuInSe26.85 132 2.88 63 CuInTe26.70 101 2.85 52 CuTlS25.12 76 2.12 52 CuTlSe22.52 43 1.07 17 AgAlS22.03 43 0.86 17 AgAlSe20.94 12 0.41 5 AgAlTe21.46 43 0.63 17 AgGaS22.23 63 0.95 25 AgGaSe20.77 6 0.33 3 AgGaTe21.43 43 0.61 17 AgInS22.05 30 0.88 12 AgInSe20.89 10 0.38 4 AgInTe21.69 110 0.72 43 a L0.5 is the mean-free-path cutoff(and therefore approximate nanostructure size) that leads to halving of κwith respect to the bulk value. Chemistry of Materials pubs.acs.org/cm Article https://doi.org/10.1021/acs.chemmater.2c00336 Chem. Mater. 2022, 34, 2833−2841 2836 CuBSe2compounds with B = Al, Ga, or In. In particular, CuTlSe2has a κvalue of 2.52 W m−1K−1, which is the lowest among the Cu-based chalcopyrites, and comparable with the Ag-based ones. The analysis presented above is based on average values along all directions (from the tensor trace). However, the calculation gives the full κtensor, so the slightly anisotropic behavior, due to the tetragonal symmetry of the structure, can be discussed. We pointed out above that Ag compounds present a larger structural anisotropy than Cu compounds, with a more pronounced deviation from a c/aratio of 2. This is reflected in the anisotropy found in the thermal conductivity. The average κz/κx,yratio for Cu-based compounds (0.91) is closer to 1 than that for the Ag-based compounds (0.84). Panels b and c of Figure 5 summarize the comparison of our results, and of previous theoretical results, with roomtemperature experimental measurements. Clearly, the calculations reported in our work offer the closest agreement so far with experimental data across the family of chalcopyrite semiconductors. This is not surprising as we have employed a more sophisticated model for the evaluation of κ, compared to previous work. For example, the results reported by Rincon et al. 42 have the highest mean absolute error (MAE), because their calculations were based on a very simple (and computationally inexpensive) analytical model, using a modified version of Leibfried and Schlomann’s equation 52 for κas a function of the Debye temperature and the Gruneisen parameter, which were obtained either directly from available experimental data or via extrapolation. Their model widely overestimates the thermal conductivities of the Ag-based chalcopyrites. More recently, Toher et al. 26 used an approach not requiring any experimental data. Thermal conductivity was calculated by combining the Slack equation 53 with the Debye temperature and Gruneisen parameter, which were both obtained from DFT calculations by using a quasi-harmonic Debye model. The lack of experimental parameters and the low computational cost make this method suitable for application to a large set of materials; however, it severely underestimates the thermal conductivity of the Cu-based chalcopyrites. Yan et al. 25 developed an approach based on the Debye− Callaway model, 54 which is predictive within 1 order of magnitude across a large range of experimental data. Although this approach improved results over those obtained from Slack’s equation, it required the fitting of parameters from experimental data. Still, the method tends to overestimate the thermal conductivities of the Ag-based chalcopyrites, and underestimate it for the Cu-based chalcopyrites, as shown in Figure 5b. The approach presented in this work improves the accuracy of the calculation (our MAE is less than half of that from the best previous theoretical work) and does so without relying on experimental data, and at a very low computational cost compared with the traditional DFT-based approach to obtain the force constants needed to solve Boltzmann’s transport equation. We can also compare the predicted evolution of κwith temperature with available experimental data. Figure 6a shows typical variations of κwith temperature T, using CuInSe2and AgInSe2as examples. The lattice thermal conductivity is dominated by phonon−phonon Umklapp scattering and therefore exhibits a T−1dependence, which reflects the increasing number of phonons available for scattering when Figure 4. Phonon density of states for (a) CuGaTe2and (b) AgGaTe2, showing the projections on the A metal atoms. (c) Group velocities vs mode frequency. (d) Scattering rates vs mode frequency for the same two compounds. Chemistry of Materials pubs.acs.org/cm Article https://doi.org/10.1021/acs.chemmater.2c00336 Chem. Mater. 2022, 34, 2833−2841 2837 the temperature increases. The theoretical and experimental curves are fitted using the expression κ κ=+ m T 0 where the parameter mcharacterizes the rate of decrease of the thermal conductivity with temperature. Figure 6b shows that our calculation gives a reasonable prediction of mvalues in comparison with experiment. Part of the discrepancy between theoretical and experimental values stems from the presence of other scattering mechanisms such as grain boundaries or fourphonon scattering in the actual compounds, which can modify the T−1relationship. 55,56 For example, in the case of CuInSe2, shown in Figure 6a, the authors of the experimental work concluded that the scattering on the crystal lattice began to dominate at temperatures above 400 K, giving rise to a deviation from the T−1law. 57 The thermal conductivities of Ag-based chalcopyrites, which are already relatively low, Figure 5. (a) Variation of the room-temperature lattice thermal conductivities of CuBX2and AgBX2with the nature of the chalcogen atom X = S, Se, or Te (ZXterms are the atomic numbers; empty symbols represent individual compounds, and solid symbols represent averages over B = Al, Ga, or In). (b) Comparison of the roomtemperature κvalues calculated in this work and in previous theoretical work with available experimental data. Black symbols (squares for Ag-based and circles for Cu-based systems) represent the results from this work, whereas symbols of other colors represent previous theoretical determinations. (c) Mean absolute error vs experiment in this work compared to those in previous theoretical determinations. The colors of the bars match the color of the symbols in panel b for each previous theoretical work reference (Rincon 1995 is ref 42; Toher 2014 is ref 26; Yan 2015 is ref 25). Experimental data from refs 9,12, and 43−51. Figure 6. (a) Calculated temperature variation of κfor representative Cu-based and Ag-based chalcopyrite semiconductors (solid symbols), in comparison with experiment (empty symbols; ref 57 for CuInSe2 and ref 9for AgInSe2). (b) Calculated temperature variation coefficients in comparison with available experimental data (ref 58 for AgGaSe2, ref 9for AgInSe2, refs 5,48, and 59 for AgGaTe2, refs 9 and 48 for AgInTe2, ref 49 for CuInTe2, ref 57 for CuInS2and CuInSe2, ref 49 for CuAlTe2, and ref 39 for CuGaTe2). Chemistry of Materials pubs.acs.org/cm Article https://doi.org/10.1021/acs.chemmater.2c00336 Chem. Mater. 2022, 34, 2833−2841 2838 decrease more slowly with temperature than those of the Cubased chalcopyrites. Finally, we discuss the effect of nanostructuring on the thermal conductivities, based on the decomposition of the contributions to κby the phonon mean free path. 60 In this approach, which has been widely used in the theoretical investigation of nanostructuring effects on thermal transport in thermoelectric materials, 61−65 the value of κcorresponding to a particular particle size Lis approximated as the cumulative contributions for all mean free paths up to L; i.e., the contributions from mean free paths longer than the particle size are subtracted. We find that the long mean-free-path contributions are large for the Cu-based chalcopyrites. For example, for CuGaSe2, particle sizes on the order of micrometers already have a significant effect on the thermal conductivity (Figure 7a). However, the size effect is much weaker in the low-κAg-based chalcopyrites, which therefore do not benefit as much from nanostructuring. To quantify the behavior of κreduction with nanostructuring at a given temperature, we have calculated the particle size (or, more exactly, the mean-free-path threshold) that leads to halving of the bulk value (Figure 7a). The values, which we call L0.5, are listed in Table 1. We observe that the lower the bulk value of κ, the smaller the L0.5 (Figure 7b). A rule of thumb emerges from that plot. The particle size that is needed to reduce the thermal conductivity by half is roughly 20 nm/W m−1K−1of the bulk value of κ,andthe proportionality is not significantly affected by temperature. This “universal”behavior for the chalcopyrite semiconductors is useful for designing nanostructuring strategies, and it is not too surprising. Scaling laws and universal behaviors in the cumulative lattice thermal conductivity have been explored and rationalized before for semiconductors. 61,66 In our case, there are some outliers, notably CuAlS2and AgInTe2, which are well above the regression line. These are, in fact, interesting cases, because they represent compositions at which nanostructuring leads to a much faster reduction of the thermal conductivity than for the average chalcopyrite, which is useful for thermoelectric applications and deserves further investigation. 4. CONCLUSIONS A combination of density functional theory simulations and machine learning regression algorithms have allowed us to calculate force constant potentials efficiently and accurately and obtain lattice thermal conductivities for a large set of semiconductors with the chalcopyrite structure. The computational cost linked to the calculations of the force constants was >1 order of magnitude lower than that of traditional approaches, which are based on systematic atomic displacements. These calculations lead to several important insights about the behavior of the thermal conductivity of I−III−VI2 chalcopyrite as a function of chemical composition, temperature, and microstructure. First, we have demonstrated that Ag-based chalcopyrites present considerably lower values of thermal conductivity than Cu-based chalcopyrites, and that this is mainly due to the lower frequencies of the vibrational modes in which Ag atoms participate. These vibrational modes overlap more effectively with acoustic modes, increasing the scattering processes and decreasing the scattering times. While no clear trends are found when the B cation is substituted, the trend with the nature of the chalcogen X atom is not monotonous down the group. The thermal conductivity of selenides is lower than those of the corresponding sulfides and tellurides. Moderate anisotropy is found for κ, being larger in Ag-based chalcopyrites than in Cu-based chalcopyrites, thus mirroring the structural anisotropy given by the c/aratio. Roomtemperature lattice thermal conductivities are accurately predicted, and the dependence of κon temperature is also in good agreement with experimental data. Finally, the effects of grain size on κhave been explored by calculating the cumulative κvalue up to a certain phonon mean free path. We showed that the particle size needed to reduce the thermal conductivity by half is roughly 20 nm/W m−1K−1of the bulk value of κ, at all temperatures of interest. This relationship provides a useful rule of thumb to facilitate the design of nanostructuring strategies to reduce thermal conductivity within this important family of semiconducting materials. Figure 7. (a) Cumulative lattice thermal conductivity from mean-free- path contributions up to distance Lfor a Cu-based chalcopyrite, indicating the effect that nanostructuring would have on their thermal conductivity. Black lines at T= 300 K; red lines at T= 700 K. (b) Correlation between L0.5 and κ, including points at 300 and 700 K. The solid line corresponds to a proportionality constant of 20 nm/W m−1K−1between L0.5 and κ. The CuAlS2and AgInTe2outliers are labeled. Chemistry of Materials pubs.acs.org/cm Article https://doi.org/10.1021/acs.chemmater.2c00336 Chem. Mater. 2022, 34, 2833−2841 2839 ■ASSOCIATED CONTENT * sıSupporting Information The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.chemmater.2c00336. Details of methods (multilinear regression and recursive feature elimination, generation of distorted structures, and selection of cutoffs) and figures showing the variation of thermal conductivities for all compounds with temperature (PDF) ■AUTHOR INFORMATION Corresponding Author Ricardo Grau-Crespo −Department of Chemistry, University of Reading, Reading RG6 6DX, U.K.; orcid.org/0000- 0001-8845-1719; Email: [email protected] Authors Jose J. Plata −Departamento de Química Física, Facultad de Química, Universidad de Sevilla, Seville 41012, Spain; orcid.org/0000-0002-0859-0450 Victor Posligua −Department of Chemistry, University of Reading, Reading RG6 6DX, U.K.; Department of Chemistry, Imperial College London, London W12 0BZ, U.K.; orcid.org/0000-0003-3375-3706 Antonio M. Márquez −Departamento de Química Física, Facultad de Química, Universidad de Sevilla, Seville 41012, Spain; orcid.org/0000-0001-6699-064X Javier Fernandez Sanz −Departamento de Química Física, Facultad de Química, Universidad de Sevilla, Seville 41012, Spain; orcid.org/0000-0003-2064-7007 Complete contact information is available at: https://pubs.acs.org/10.1021/acs.chemmater.2c00336 Notes The authors declare no competing financial interest. Input and output files from our simulations are available at the online repository 10.5281/zenodo.5510656. Our wrapper code for HiPhive can be downloaded from https://github.com/ NewMaterialsLab/hiPhive_wrapper_NML. ■ACKNOWLEDGMENTS This work was funded by the Ministerio de Ciencia e Innovación (PID2019-106871GB-I00). J.J.P., A.M.M., and J.F.S. thankfully acknowledge the computer resources at Lusitania and the technical support provided by Cénits- COMPUTAEX and Red Espanola de Supercomputación, RES (QS-2021-1-0027 and QHS-2021-2-0022). V.P. and R.G.-C. are grateful to the UK Materials and Molecular Modelling Hub for computational resources (Young facility), which is partially funded by EPSRC (EP/P020194/1 and EP/T022213/1). ■REFERENCES (1) Shay, J. L.; Wernick, J. H. Ternary chalcopyrite semiconductors: growth, electronic properties, and applications; Pergamon Press: Oxford, U.K., 1975. (2) Nakamura, M.; Yamaguchi, K.; Kimoto, Y.; Yasaki, Y.; Kato, T.; Sugimoto, H. Cd-free Cu(In,Ga)(Se,S)2thin-film solar cell with record efficiency of 23.35%. IEEE Journal of Photovoltaics 2019,9, 1863−1867. (3) Wada, T. CuInSe2and related I-III-VI2chalcopyrite compounds for photovoltaic application. Jpn. J. Appl. Phys. 2021,60, 080101. (4) Thomas, S. R.; Chen, C.-W.; Date, M.; Wang, Y.-C.; Tsai, H.-W.; Wang, Z. M.; Chueh, Y.-L. Recent developments in the synthesis of nanostructured chalcopyrite materials and their applications: a review. RSC Adv. 2016,6, 60643−60656. (5) Cao, Y.; Su, X. L.; Meng, F. C.; Bailey, T. P.; Zhao, J. G.; Xie, H. Y.; He, J.; Uher, C.; Tang, X. F. Origin of the Distinct Thermoelectric Transport Properties of Chalcopyrite ABTe2(A = Cu, Ag; B = Ga, In). Adv. Funct. Mater. 2020,30, 2005861. (6) Lu, Y.; Chen, S.; Wu, W.; Du, Z.; Chao, Y.; Cui, J. Enhanced thermoelectric performance of a chalcopyrite compound CuIn3Se5−xTex(x=0∼0.5) through crystal structure engineering. Sci. Rep. 2017,7, 40224. (7) Nomura, A.; Choi, S.; Ishimaru, M.; Kosuga, A.; Chasapis, T.; Ohno, S.; Snyder, G. J.; Ohishi, Y.; Muta, H.; Yamanaka, S.; Kurosaki, K. Chalcopyrite ZnSnSb2: A Promising Thermoelectric Material. ACS Appl. Mater. Interfaces 2018,10, 43682−43690. (8) Plirdpring, T.; Kurosaki, K.; Kosuga, A.; Day, T.; Firdosy, S.; Ravi, V.; Snyder, G. J.; Harnwunggmoung, A.; Sugahara, T.; Ohishi, Y.; et al. Chalcopyrite CuGaTe2: a high-efficiency bulk thermoelectric material. Adv. Mater. 2012,24, 3622−3626. (9) Qiu, P.; Qin, Y.; Zhang, Q.; Li, R.; Yang, J.; Song, Q.; Tang, Y.; Bai, S.; Shi, X.; Chen, L. Intrinsically High Thermoelectric Performance in AgInSe2n-Type Diamond-Like Compounds. Advanced Science 2018,5, 1700727. (10) Yusufu, A.; Kurosaki, K.; Kosuga, A.; Sugahara, T.; Ohishi, Y.; Muta, H.; Yamanaka, S. Thermoelectric properties of Ag1−xGaTe2 with chalcopyrite structure. Appl. Phys. Lett. 2011,99, 061902. (11) Wang, B.; Xiang, H.; Nakayama, T.; Zhou, J.; Li, B. Theoretical investigation on thermoelectric properties of Cu-based chalcopyrite compounds. Phys. Rev. B 2017,95, 035201. (12) Parker, D.; Singh, D. J. Thermoelectric properties of AgGaTe2 and related chalcopyrite structure materials. Phys. Rev. B 2012,85, 125209. (13) Kumar Gudelli, V.; Kanchana, V.; Vaitheeswaran, G.; Svane, A.; Christensen, N. E. Thermoelectric properties of chalcopyrite type CuGaTe2and chalcostibite CuSbS2.J. Appl. Phys. 2013,114, 223707. (14) Li, R.; Li, X.; Xi, L.; Yang, J.; Singh, D. J.; Zhang, W. Highthroughput screening for advanced thermoelectric materials: diamond-like ABX2compounds. ACS Appl. Mater. Interfaces 2019, 11, 24859−24866. (15) Ziman, J. M. Electrons and Phonons: The Theory of Transport Phenomena in Solids; Oxford University Press: Oxford, U.K., 2001; p 568. (16) Togo, A.; Chaput, L.; Tanaka, I. Distributions of phonon lifetimes in Brillouin zones. Phys. Rev. B 2015,91, 094306. (17) Li, W.; Carrete, J.; A. Katcho, N.; Mingo, N. ShengBTE: A solver of the Boltzmann transport equation for phonons. Comput. Phys. Commun. 2014,185, 1747−1758. (18) Plata, J. J.; Nath, P.; Usanmaz, D.; Carrete, J.; Toher, C.; de Jong, M.; Asta, M.; Fornari, M.; Nardelli, M. B.; Curtarolo, S. An efficient and accurate framework for calculating lattice thermal conductivity of solids: AFLOWAAPL Automatic Anharmonic Phonon Library. npj Comput. Mater. 2017,3, 45. (19) Chernatynskiy, A.; Phillpot, S. R. Phonon Transport Simulator (PhonTS). Comput. Phys. Commun. 2015,192, 196−204. (20) Tadano, T.; Gohda, Y.; Tsuneyuki, S. Anharmonic force constants extracted from first-principles molecular dynamics: applications to heat transfer simulations. J. Phys.: Condens. Matter 2014,26, 225402. (21) Toher, C.; Oses, C.; Plata, J. J.; Hicks, D.; Rose, F.; Levy, O.; de Jong, M.; Asta, M.; Fornari, M.; Buongiorno Nardelli, M.; Curtarolo, S. Combining the AFLOW GIBBS and elastic libraries to efficiently and robustly screen thermomechanical properties of solids. Physical Review Materials 2017,1, 015401. (22) Nath, P.; Plata, J. J.; Usanmaz, D.; Toher, C.; Fornari, M.; Buongiorno Nardelli, M.; Curtarolo, S. High throughput combinatorial method for fast and robust prediction of lattice thermal conductivity. Scripta Materialia 2017,129,88−93. (23) Bjerg, L.; Iversen, B. B.; Madsen, G. K. H. Modeling the thermal conductivities of the zinc antimonides ZnSb and Zn4Sb3. Phys. Rev. B 2014,89, 024304. Chemistry of Materials pubs.acs.org/cm Article https://doi.org/10.1021/acs.chemmater.2c00336 Chem. Mater. 2022, 34, 2833−2841 2840 (24) Toberer, E. S.; Zevalkink, A.; Snyder, G. J. Phonon engineering through crystal chemistry. J. Mater. Chem. 2011,21, 15843−15852. (25) Yan, J.; Gorai, P.; Ortiz, B.; Miller, S.; Barnett, S. A.; Mason, T.; Stevanovic, V.; Toberer, E. S. Material descriptors for predicting thermoelectric performance. Energy Environ. Sci. 2015,8, 983−994. (26) Toher, C.; Plata, J. J.; Levy, O.; De Jong, M.; Asta, M.; Nardelli, M. B.; Curtarolo, S. High-throughput computational screening of thermal conductivity, Debye temperature, and Gruneisen parameter using a quasiharmonic Debye model. Phys. Rev. B 2014,90, 174107. (27) Zhou, F.; Nielson, W.; Xia, Y.; Ozolins,V.Lattice Anharmonicity and Thermal Conductivity from Compressive Sensing of First-Principles Calculations. Phys. Rev. Lett. 2014,113, 185501. (28) Eriksson, F.; Fransson, E.; Erhart, P. The Hiphive Package for the Extraction of High-Order Force Constants by Machine Learning. Advanced Theory and Simulations 2019,2, 1800184. (29) Xia, Y.; Hegde, V. I.; Pal, K.; Hua, X.; Gaines, D.; Patel, S.; He, J.; Aykol, M.; Wolverton, C. High-Throughput Study of Lattice Thermal Conductivity in Binary Rocksalt and Zinc Blende Compounds Including Higher-Order Anharmonicity. Physical Review X2020,10, 041029. (30) Kresse, G.; Furthmuller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B1996,54, 11169−11186. (31) Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B 1999,59, 1758− 1775. (32) Calderon, C. E.; Plata, J. J.; Toher, C.; Oses, C.; Levy, O.; Fornari, M.; Natan, A.; Mehl, M. J.; Hart, G.; Buongiorno Nardelli, M.; Curtarolo, S. The AFLOW standard for high-throughput materials science calculations. Comput. Mater. Sci. 2015,108, 233− 238. (33) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996,77, 3865−3868. (34) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (35) Park, J.; Xia, Y.; Ozolins, V. First-principles assessment of thermoelectric properties of CuFeS2.J. Appl. Phys. 2019,125, 125102. (36) Guyon, I.; Weston, J.; Barnhill, S.; Vapnik, V. Gene selection for cancer classification using support vector machines. Machine learning 2002,46, 389−422. (37) Plata, J. J.; Posligua, V.; Marquez, A. M.; Fernandez Sanz, J.; Grau-Crespo, R. hiPhive Wrapper by NML. https://github.com/ NewMaterialsLab/hiPhive_wrapper_NML (accessed 2022-02-15). (38) Skelton, J. M.; Tiana, D.; Parker, S. C.; Togo, A.; Tanaka, I.; Walsh, A. Influence of the exchange-correlation functional on the quasi-harmonic lattice dynamics of II-VI semiconductors. J. Chem. Phys. 2015,143, 064710. (39) Bodnar, I. V.; Victorov, I. A.; Dabranski, V. M. Crystal growth and investigation of CuAlxGa1−xTe2solid solutions. J. Cryst. Growth 2004,265, 214−219. (40) Wei, S.-H.; Zunger, A. Electronic and structural anomalies in lead chalcogenides. Phys. Rev. B 1997,55, 13605−13610. (41) Zhu, Y.; Wei, B.; Liu, J.; Koocher, N. Z.; Li, Y.; Hu, L.; He, W.; Deng, G.; Xu, W.; Wang, X.; Rondinelli, J. M.; Zhao, L.-D.; Snyder, G. J.; Hong, J. Physical insights on the low lattice thermal conductivity of AgInSe2.Materials Today Physics 2021,19, 100428. (42) Rincon, C.; Valerigil, M. L.; Wasim, S. M. Room-Temperature Thermal-Conductivity and Gruneisen-Parameter of the I-III-VI2 Chalcopyrite Compounds. Physica Status Solidi a-Applications and Materials Science 1995,147, 409−415. (43) Bodnar, I. V. Growth and properties of CuAlSe2single crystals. Inorg. Mater. 2005,41, 560. (44) Carr, W. D. Electronic and thermal transport in copper-based chalcopyrite semiconductors for thermoelectric applications. PhD Thesis, Michigan State University, East Lansing, MI, 2016. (45) Berger, L. I.; Prochukhan, V. D. Ternary Diamond-Like Semiconductors; Springer, 1995. (46) Wasim, S. M.; Noguera, A. Thermal conductivity of CuGaTe2. Solid State Commun. 1987,64, 439−442. (47) Beasley, J. D. Thermal conductivity measurements in nonlinear optical materials. Proceedings Volume 1869, Intracavity and Extracavity Control of Laser Beam Properties; 1993; Vol. 1869. (48) Charoenphakdee, A.; Kurosaki, K.; Muta, H.; Uno, M.; Yamanaka, S. Thermal Conductivity of the Ternary Compounds: AgMTe2and AgM5Te8(M= Ga or In). Materials Transactions 2009, 50, 1603−1606. (49) Bodnar, I. V. Optical and thermal properties of CuAlxIn1−xTe2 solid solutions. Semiconductors 2003,37, 1247−1251. (50) Wasim, S. M. Thermal conductivity of ternary Compounds. Phys. Status Solidi A 1979,51, K35−K40. (51) Gherouel, D.; Gaied, I.; Amlouk, M. Effect of heat treatment in air on physical properties of AgInS2sprayed thin films. J. Alloys Compd. 2013,566, 147−155. (52) Leibfried, G.; Schlomann, E. Warmeleitung in elektrisch isolierenden kristallen. Nachr. Ges. Wiss. Goett. Math. -Phys. Kl 1954, 2, 71. (53) Slack, G. A. The Thermal Conductivity of Nonmetallic Crystals. In Solid State Physics; Ehrenreich, H., Seitz, F., Turnbull, D., Eds.; Academic Press, 1979; Vol. 34,pp1−71. (54) Callaway, J. Model for Lattice Thermal Conductivity at Low Temperatures. Phys. Rev. 1959,113, 1046−1051. (55) Zhu, Y.; Xia, Y.; Wang, Y.; Sheng, Y.; Yang, J.; Fu, C.; Li, A.; Zhu, T.; Luo, J.; Wolverton, C.; Snyder, G. J.; Liu, J.; Zhang, W. Violation of the T−1Relationship in the Lattice Thermal Conductivity of Mg3Sb2with Locally Asymmetric Vibrations. Research 2020,2020, 4589786. (56) Xia, Y.; Pal, K.; He, J.; Ozolins, V.; Wolverton, C. Particlelike Phonon Propagation Dominates Ultralow Lattice Thermal Conductivity in Crystalline Tl3VSe4.Phys. Rev. Lett. 2020,124, 065901. (57) Makovetskaya, L. A.; Bodnar, I. V.; Korzun, B. V.; Yaroshevich, G. P. Thermal Conductivity, Thermoelectric Power, and Thermal Expansion of CuInS2xSe2(1−x).Phys. Status Solidi A 1982,74, K59− K62. (58) Aggarwal, R. L.; Fan, T. Y. Thermal diffusivity, specific heat, thermal conductivity, coefficient of thermal expansion, and refractiveindex change with temperature in AgGaSe2.Appl. Opt. 2005,44, 2673−2677. (59) Yusufu, A.; Kurosaki, K.; Ohishi, Y.; Muta, H.; Yamanaka, S. Effect of Cu Doping into the Ga Site on the Thermoelectric Properties of AgGaTe2with Chalcopyrite Structure. J. Jpn. Soc. Powder Powder Metallurgy 2012,59, 206−209. (60) Dames, C.; Chen, G. Thermal conductivity of nanostructured thermoelectric materials. In Thermoelectrics Handbook: Macro to Nano; Rowe, D., Ed.; CRC Taylor & Francis: Boca Raton, FL, 2006; p 421. (61) Aketo, D.; Shiga, T.; Shiomi, J. Scaling laws of cumulative thermal conductivity for short and long phonon mean free paths. Appl. Phys. Lett. 2014,105, 131901. (62) Aziz, A.; Mangelis, P.; Vaqueiro, P.; Powell, A. V.; Grau-Crespo, R. Electron and phonon transport in shandite-structured Ni3Sn2S2. Phys. Rev. B 2016,94, 165131. (63) Shafique, A.; Samad, A.; Shin, Y.-H. Ultra low lattice thermal conductivity and high carrier mobility of monolayer SnS2and SnSe2:a first principles study. Phys. Chem. Chem. Phys. 2017,19, 20677− 20683. (64) Romano, G.; Esfarjani, K.; Strubbe, D. A.; Broido, D.; Kolpak, A. M. Temperature-dependent thermal conductivity in silicon nanostructured materials studied by the Boltzmann transport equation. Phys. Rev. B 2016,93, 035408. (65) Ma, J.; Li, W.; Luo, X. Intrinsic thermal conductivity and its anisotropy of wurtzite InN. Appl. Phys. Lett. 2014,105, 082103. (66) Freedman, J. P.; Leach, J. H.; Preble, E. A.; Sitar, Z.; Davis, R. F.; Malen, J. A. Universal phonon mean free path spectra in crystalline semiconductors at high temperature. Sci. Rep. 2013,3, 2963. Chemistry of Materials pubs.acs.org/cm Article https://doi.org/10.1021/acs.chemmater.2c00336 Chem. Mater. 2022, 34, 2833−2841 2841