Design Principles Guided by DFT Calculations and High-Throughput Frameworks for the Discovery of New Diamond-like Chalcogenide Thermoelectric Materials

Design Principles Guided by DFT Calculations and High-Throughput Frameworks for the Discovery of New Diamond-like Chalcogenide Thermoelectric Materials Adolfo E. Rosado-Miranda, § Victor Posligua, § Javier Fdez. Sanz, Antonio M. Márquez, Pinku Nath,* and Jose J. Plata* Cite This: ACS Appl. Mater. Interfaces 2024, 16, 28590−28598 Read Online ACCESS Metrics & More Article Recommendations * sı Supporting Information ABSTRACT: Rational design principles are one pathway to discovering new materials. However, technological breakthroughs rarely occur in this way because these design principles are usually based on incremental advances that seldom lead to disruptive applications. The emergence of machine-learning (ML) and high-throughput (HT) techniques has changed the paradigm, opening up new possibilities for efficiently screening large chemical spaces and creating on-the-fly design principles for the discovery of novel materials with desired properties. In this work, the approach is used to discover novel thermoelectric (TE) materials based on quaternary diamond-like chalcogenides. A HT framework that integrates density functional theory calculations, ML, and the solution of the Boltzmann transport equation is used to efficiently rationalize the transport properties of these compounds and identify those with potential as TE materials, achieving ZT values above 2. KEYWORDS: thermoelectrics, chalcogenides, design principles, ZT, transport properties 1. INTRODUCTION A large percentage of technological breakthroughs is sustained on materials discovery. Traditionally, there have been three main pathways to the synthesis of new materials. 1 Some compounds are fortuitously discovered through serendipity, while pursuing different objectives. Others have been known for many years until their potential applications have been recognized and explored. The third pathway involves the use of design principles to predict and synthesize materials with welldefined properties. However, the emergence of highthroughput (HT), first-principles simulations and artificial intelligence (AI) has revolutionized the process of materials discovery and design. 2 These methods facilitate the efficient screening of large compositional spaces, leading to the identification of novel compounds with targeted properties. Although these methods are often viewed as substitutes for chemical and physical intuition, they actually accelerate and enrich the development of new design principles. While the extraction of materials design principles from the combination of HT calculations and machine learning (ML) has led to significant advances in many technological areas, this strategy faces significant challenges in the development of thermoelectric (TE) materials. Establishing design rules for materials applications with a single objective is trivial; however, the scenario is completely different when multiple and competing objectives are targeted. 3 TE materials are particularly challenging due to the complex interplay of electrical conductivity (σ), thermal conductivity (κ), and Seebeck coefficient (S). Achieving high TE performance Received: March 12, 2024 Revised: May 2, 2024 Accepted: May 13, 2024 Published: May 21, 2024 Research Articlewww.acsami.org © 2024 The Authors. Published by American Chemical Society 28590 https://doi.org/10.1021/acsami.4c04120 ACS Appl. Mater. Interfaces 2024, 16, 28590−28598 This article is licensed under CC-BY 4.0 Downloaded via UNIV DE SEVILLA on October 14, 2024 at 14:03:54 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles. requires a delicate balance between these properties, which is often difficult to achieve through conventional materials design strategies. Additionally, the HT calculation of TE material properties as the main source for the application of AI/ML techniques presents severe limitations. The accurate prediction of electronic and thermal transport properties is several orders of magnitude more expensive than the other material properties linked to other technological applications. 4 Additionally, transport properties strongly depend on many variables that can be modified during the synthesis and processing of the material but are not included in the simulation of bulk, stoichiometric, and defect-free materials. During the past decade, the development of new methods and frameworks for the accurate prediction of thermal and electronic transport properties has accelerated the search for new and more efficient TE materials. 5,6 Approaches based on the ML-aided extraction of high-order force constants have reduced the computational cost of simulating thermal transport properties without compromising accuracy. 7,8 Similarly, a new whole family of packages has been developed for the accurate prediction of electronic transport properties, 9,10 without the need for extensive computational resources required by the combination of Wannier functions and density functional pertubation theory (DFPT). 11 For instance, Pal et al. have conducted an extensive screening of thousands of quaternary chalcogenides (AMM’X3) predicting their stability and computing the thermal conductivity of a reduced set of them. 12,13 Chalcogenides are one of the largest family of compounds in the TE field due to their thermal and electronic features. 14 Despite the large variety of compositions and structural prototypes, diamond-like (DL) compounds play a preponderant role due to their usually low cost and simple synthesis. DL chalcogenides are a rich set of materials in which all cations present tetrahedral coordination. Although binary DL chalcogenides have limited applications as TE materials, several high-performance alternatives can be found among ternary DL chalcogenides, such as CuFeS2, 15 AgInSe2, 16 and AgInTe2. 17 However, quaternary DL chalcogenides are the best example of how these compounds can reach a figure of merit, ZT, of above 1.5. There are some examples in which the TE properties of some subset of quaternary DL chalcogenides have been explored guided by density functional theory (DFT) calculations. 18,19 However, there is a need for design principles to effectively explore such an extensive chemical space. In this work, we integrate HT frameworks and ML to accurately and systematically predict the lattice thermal conductivity of quaternary DL chalcogenides. These findings serve as a foundation for establishing a set of design rules to guide the search for new and more efficient TE materials. Through the application of these principles, we have not only explored but successfully optimized the TE performance of novel quaternary DL chalcogenides, surpassing ZT values of 2. 2. RESULTS AND DISCUSSION 2.1. Validation. While most of the methodology applied here has been validated in previous studies using DL ternary chalcogenides 20 and pnictides, 21 it is essential to verify the accuracy and reliability of our computational framework. Accurate determination of interatomic distances is crucial for obtaining precise κlvalues. The calculated lattice parameters are compared with experimental data for a selected group of materials to assess the efficacy of the DFT methodology in describing their structural properties. The calculated lattice parameters show a strong agreement with experimental values, with errors below 3% in all cases (Figure S1). For instance, the local density approximation (LDA) and the Perdew−Burke− Ernzerhof (PBE) functionals tend to underestimate and overestimate the experimental lattice parameter of Cu2CdSnTe4, respectively. In contrast, the PBE-D3 functional provides a much closer prediction to the experiment. Optimizing the calculation of κlfrom first principles involves considering factors beyond just the exchange−correlation functional. Interaction cutoff radii are used to reduce the computational costs. While harmonic (second-order) force constant cutoffs were selected based on the largest sphere that can be fit inside the supercell, third and fourth-order interatomic force constants (IFCs) were chosen according to previous study on ternary chalcogenides. 20 A good agreement was achieved between phonon dispersion curves obtained using the finite-differences approach and the ML-aided methodology for Cu2CdSnTe4(Figure 1a). The number of distorted supercells used in training the ML model can also modify the accuracy of the IFCs prediction, especially for thirdand fourth-order force constants. It was found that using 30 supercells was sufficient to achieve a root mean squared error (RMSE) in forces below 22 meV/Å for the validation step as well as converge κlvalues (Figure S2). The calculated values of κlshow a slight overestimation when compared to the experimental values reported for polycrystalline samples by Nolas et al. 22,23 (Figure 1b). However, note that these calculated values do not take into account scattering processes caused by grain boundaries. Experimental reports have described the grains as large, with sizes exceeding 100 nm. 23 Figure 1. (a) Phonon dispersion curves for stannite Cu2CdSnTe4 using finite difference approach, full DFT (solid blue line) and ML regression, ML (dashed green line). (b) Lattice thermal conductivity, κl, convergence with the density of the q-point mesh at 300 K. κl values are depicted for the Cu2CdSnTe4single-crystal (blue solid line) and polycrystalline sample (dashed blue line). Green dashed line and green area represent reported experimental values. ACS Applied Materials & Interfaces www.acsami.org Research Article https://doi.org/10.1021/acsami.4c04120 ACS Appl. Mater. Interfaces 2024, 16, 28590−28598 28591 If grain boundary scattering is included in the model, taking an average grain size of 500 nm, our predictions align much better with experimental results, showing less than a 5% error (Figure 1b). 2.2. Monovalent Atom Role. Extensive research has been focused on Cu-based chalcogenides due to their promising TE properties. 20,24 However, it has been reported that Ag-based ternary chalcogenides present a significantly lower thermal conductivity than their Cu counterparts. 25 Nevertheless, it remains unclear whether this trend can be generalized to other DL chalcogenides, and the mechanisms responsible for the reduced thermal conductivity in Ag-based chalcogenides are still not fully understood. To understand the role of the monovalent cation, the lattice thermal conductivity of Cu2CdSnTe4and Ag2CdSnTe4in their stannite structure was calculated. Phonon density of states (pDOS) already shows a clear difference between both solids due to the substitution of Cu by Ag (Figure 2a,b). The larger atomic mass of Ag compared to that of Cu significantly influences the contribution of the monovalent atom to the acoustic modes and low-energy optic modes. Ag-based compound presents an important contribution of the monovalent atom to the acoustic modes, while its contribution in the Cu counterpart is only noticeable for the low-energy optical modes. The substitution of Cu by Ag also introduces severe modifications in the anharmonicity of the materials through the phonon scattering rates, especially at low energies (Figure 2c). Ag2CdSnTe4 presents phonon scattering rates around 1 order of magnitude larger than those of Cu2CdSnTe4in the low-energy range, anticipating lower lattice thermal conductivity values. Cumulative lattice thermal conductivity corroborates that the four times lower κlof the Ag-based compound is due to the higher scattering rates of both acoustic and low-energy optic modes (up to 1.75 THz). The origin of this difference in scattering rate values can be explained by comparing the phase space and the third-order IFCs of both compounds (Figure S3). Interestingly, Cu2CdSnTe4presents larger values for the phase space than Ag2CdSnTe4, which would lead to higher scattering rates. This could be due to a more effective overlapping between vibrational modes in Cu2CdSnTe4as observed in pDOS (Figure 2a,b). However, third-order IFCs are generally larger for Ag2CdSnTe4, explaining the larger scattering rates compared to the Cu-based compound. 2.3. Polymorphs. One notable contrast between ternary and quaternary DL chalcogenides is the potential of several patterns or ordering motifs within the cation sublattice, which leads to the presence of different structural prototypes. In this work, all previous materials have been explored using the stannite polymorph (I42m); however, it is worth noting that kesterite (I4) and primitive mixed CuAu-PMCA (P42m) polymorphs have been also experimentally reported. Indeed, stannite and kesterite 0 K energies rarely differ by more than 10 meV/atom, making it challenging to differentiate them experimentally without advanced structural characterization techniques. To explore the role of different structural prototypes in the thermal properties of quaternary DL chalcogenides, the κlof Cu2ZnSnSe4in both stannite and kesterite structures (Figure 3a) was calculated. The comparison of the pDOS for both polymorphs exhibited negligible differences (Figure S4). Both materials exhibit similar acoustic bands below 3 THz and two well-defined optic bands. The Cu2ZnSnSe4stannite polymorph shows slightly higher lattice thermal conductivity compared to the kesterite structure in the whole temperature range (Figure 3b). This slight variation is attributed to the higher scattering rates of the acoustic band in the kesterite structure (Figure 3c). Approximately, 80% of the total thermal conductivity, κl, in both compounds comes from vibrational modes with frequencies below 1.5 THz. The difference in thermal conductivity between the kesterite and stannite compounds can be attributed to their crystal symmetries. Kesterite space group is lower than stannite. This symmetry break promotes the appearance of vibrational modes with different frequencies, enhancing the number of scattering processes and reducing the κlvalue. Besides the difference between both polymorphs, it is worth noting the anisotropy of both compounds in terms of thermal conductivity. Both structural prototypes are constructed by stacking tetrahedra layers in the z-direction, creating a more homogeneous environment in the xy-plane. Consequently, κl(xx) is consistently higher than κl(zz) for both compounds. 2.4. Combining Design Principles. Based on the results presented in this study, optimizing κlto enhance ZT in quaternary DL chalcogenides is associated with two factors: (1) including Ag atoms and (2) synthesizing these compounds primarily in their kesterite form. However, a key challenge in maximizing ZT lies in the interconnection between the thermal Figure 2. Phonon density of states (pDOS) for (a) Cu2CdSnTe4and (b) Ag2CdSnTe4. Total pDOS is depicted using solid blue lines and the monovalent atom projection using dashed green lines. Cumulative lattice thermal conductivity as a function of phonon frequencies is included in dash-dotted grey lines. (c) Scattering rates vs mode frequency for Cu2CdSnTe4(blue) and Ag2CdSnTe4(green). ACS Applied Materials & Interfaces www.acsami.org Research Article https://doi.org/10.1021/acsami.4c04120 ACS Appl. Mater. Interfaces 2024, 16, 28590−28598 28592 and electronic transport properties of the material. For instance, the complete substitution of Cu by Ag in chalcogenides drastically reduces the carrier concentration in the compound, resulting in poor power factors (PF). Consequently, the presence of both Cu and Ag simultaneously has been reported in quaternary DL compounds, exhibiting outstanding performance. 26 Stability and cost are two additional factors that should be considered beyond efficiency. Te-based compounds have been characterized in detail because of their high performance; 27,28 however, the volatility of Te reduces their stability. Besides offering higher stability, the use of selenides instead of tellurides opens the door to a reduction in production costs. Selenium is approximately 50 times more abundant than Te in the Earth’s crust. Combining these design principles based on optimizing performance, stability, and production cost, the list of potential quaternary DL chalcogenides is drastically reduced. In this work, kesterite In2AgCuSe4is proposed as a candidate as it fulfills all the aforementioned requirements. In2AgCuSe4dispersion curves exhibit characteristics consistent with low thermal conductivity (Figure 4a). These features include the presence of low-lying optical modes and flat bands. Optical modes overlap with acoustic modes in the low-frequency region (around 1 THz), resulting in increased phonon scattering and decreased thermal conductivity. At 300 K, In2AgCuSe4exhibits a lattice thermal conductivity of around 2.6 W/mK, with approximately 70% attributed to acoustic modes below 1 THz (Figure 4b). Projected pDOS shows that Ag, In, and Se are the main contributors to these low energy modes (Figure 4b). Flat dispersion bands indicate low group velocities. Only a few modes in the 0−1 and 2−3 THz intervals have group velocities greater than 1 nmTHz (Figure 4c). Large scattering rates, similar to the ones obtained for Ag2CdSnTe4, are a clear sign of the anharmonicity of the material (Figure 4d). The lattice thermal conductivity tensor as a function of the temperature is depicted in Figure 4e. The anisotropy of κl, measured as κl(zz)/κl(xx), for this compound at 300 K is around 0.83. This value is significantly larger than the values reported for ternary chalcopyrites 20 and closer to those of ternary pnictides. 21 The lattice thermal conductivity is mainly influenced by phonon−phonon Umklapp scattering, resulting in a T−1dependence (−1.075 when data is fitted). This reflects the growing number of phonons available for scattering as temperature increases. Lattice thermal conductivity can also be drastically reduced by nanostructuring. To quantify the effect of nanostructuring on κlreduction at a specific temperature, we calculated the particle size (or more precisely, the mean-freepath threshold) that results in a 50% decrease in the bulk value of lattice thermal conductivity, L0.5 (Figure 4f). Although In2AgCuSe4presents low L0.5 values (below 100 nm for both 300 and 700 K), it is possible to reduce the κlof this compound around 20−40% by nanostructuring it to a particle size between 100 and 300 nm at both 300 and 700 K. Low lattice thermal conductivity alone is not sufficient to ensure high TE performance; electronic transport properties also play a critical role. The electronic band structure presents two main features that govern the electronic transport properties: the band gap and effective masses. The PBE functional predicts a band gap of less than 0.1 eV (Figure 5a), while the HSE functional obtained a gap of approximately 0.6 eV, in close agreement with the predicted band gaps for CuInSe2and AgInSe2, which ranged between 0.5 and 1 eV. 29,30 Due to significant differences in the topology of bands at the edge of the valence and conduction bands, variations should be observed between pand n-type behavior of the material. The conduction band presents a Γ-centered main pocket with a relatively low effective mass, indicating high electron mobility, μ. However, the valence band presents different pockets close in energy at the Γpoint. These pockets have different effective masses, indicating a more complex behavior for hole conductivity based on their occupancy. The qualitative predictions from the band structure are supported by the behavior of Sand σ(Figure 5b,c). For n-type In2AgCuSe4, high mobility, moderate S, high σ, and moderate PF are observed (Figure 5d). However, p-type In2AgCuSe4shows a more complex behavior, as anticipated by the band structure. At low hole concentrations and temperatures, this material exhibits a Figure 3. (a) Crystal structures for Cu2ZnSnSe4kesterite (left) and stannite (right). Color code: black, Cu; gray, Sn; red, Zn; yellow, Se. (b) Lattice thermal conductivity as a function of temperature for Cu2ZnSnSe4kesterite (blue) and stannite (green). Solid and dashed lines are used for in-plane (xx) and perpendicular (zz)κl, respectively. (c) Scattering rates vs mode frequency for Cu2ZnSnSe4kesterite (blue) and stannite (green). ACS Applied Materials & Interfaces www.acsami.org Research Article https://doi.org/10.1021/acsami.4c04120 ACS Appl. Mater. Interfaces 2024, 16, 28590−28598 28593 large Seebeck coefficient in agreement with the large effective masses and low mobility of states at the edge of the valence band (Figure 5b). However, bands with higher mobility near energy become populated at higher temperatures, leading to a drastic drop in Sat higher temperatures and low carrier concentrations. This trend changes at higher carrier concentrations where more bands become populated at lower temperatures. The presence of quasi-degenerate states also enhances electrical conductivity, resulting in remarkable values of PF, especially at high carrier concentrations (Figure 5d). Once the primary trends have been examined based on temperature and carrier concentration, it is essential to investigate the significance of each type of scattering mechanism incorporated in the theoretical model. This step is crucial for rationalizing how the electronic transport properties can be altered in In2AgCuSe4. An effective method for this purpose is to evaluate the contributions of each scattering mechanism to the carrier mobility (Figure S5). It looks like polar optical phonon (POP) scattering prevails as the dominant mechanism affecting carrier mobility in In2AgCuSe4across all temperatures and carrier concentrations (Figure S5). The electrons are scattered through the interaction of the Coulomb field of the lattice polarization waves due to optical vibrational modes. Only at very high carrier concentrations, ionized impurity (IMP), and scattering processes contribute to total mobility for p-type In2AgCuSe4. This finding once again emphasizes a strong correlation between phonon vibrational structure and transport properties. While acoustic modes are crucial for thermal transport properties, optic modes also significantly influence electron mobility within this compound. The TE figure of merit, ZT, for p-type and n-type In2AgCuSe4has been calculated by combining κand PF over a wide range of temperatures, T, and carrier concentrations, n (Figure 6). Similar to other DL chalcogenides, p-type In2AgCuSe4exhibits higher ZT values compared to those of its n-type counterpart, particularly at high carrier concentrations. While the maximum ZT for bulk n-type In2AgCuSe4is around 0.6−0.7 at moderate carrier concentrations and high temperatures, bulk p-type In2AgCuSe4shows a maximum ZT of around 1.8 at high carrier concentration and temperature (Figure 6). This represents not only quantitative but also qualitative differences between both semiconductors especially in the TE optimization process. These differences stem from the distinct topology described at the edges of the conduction and valence band. The parabolic single band at the bottom of the conduction band makes temperature the main variable for maximizing ZT by reducing κl. Nevertheless, the complex Figure 4. (a) Phonon dispersion curves for kesterite In2AgCuSe4. (b) Phonon density of states (pDOS) for kesterite In2AgCuSe4. Total pDOS is depicted using solid blue lines. In, Cu, Ag, and Se atom projections are depicted using yellow, red, purple, and green solid lines, respectively. Cumulative lattice thermal conductivity as a function of phonon frequencies is included in dashed grey lines. (c,d) In2AgCuSe4group velocities and scattering rates as a function of phonon frequencies. (e) Lattice thermal conductivity as a function of temperature for In2AgCuSe4kesterite. Solid, dashed, and dotted lines are used for in-plane (xx), perpendicular (zz), and average κl, respectively. (f) Cumulative lattice thermal conductivity from mean-free-path contributions up to distance Lfor In2AgCuSe4, illustrating the effect that nanostructuring would have on thermal conductivity. Solid and dashed lines denote T= 300 and 700 K, respectively. ACS Applied Materials & Interfaces www.acsami.org Research Article https://doi.org/10.1021/acsami.4c04120 ACS Appl. Mater. Interfaces 2024, 16, 28590−28598 28594 topology of the top of the valence band, featuring multiple pockets with similar energies, yet significantly different effective masses, opens up opportunities for ZT maximization. This can be achieved not only by reducing κlas the temperature increases but also by increasing hole concentration. To the best of our knowledge, there are no previous studies on the synthesis and carrier concentration of In2AgCuSe4. However, there are many reports about the dopability of CuInSe2and AgInSe2. As expected, CuInSe2 Figure 5. (a) Electronic band structure for In2AgCuSe4. (b−d) Seebeck coefficient (S), electrical conductivity (σ), and power factor (PF) for p- (green) and n-type (blue) In2AgCuSe4as a function of temperature and carrier concentration, respectively. Figure 6. TE figure of merit, ZT, dependence on temperature, T, carrier concentration, n, and grain size, L, for p- (green) n-type (blue) In2AgCuSe4. ACS Applied Materials & Interfaces www.acsami.org Research Article https://doi.org/10.1021/acsami.4c04120 ACS Appl. Mater. Interfaces 2024, 16, 28590−28598 28595 presents higher carrier concentration levels (1019 cm−3as ptype) than AgInSe2(2.5 ×1018 cm−3as n-type). 31 Optimizing ZT requires taking into account factors beyond the temperature and carrier concentration. The microstructure plays a critical role in determining the TE properties of DL chalcogenides. Many strategies have been reported for the synthesis and processing of DL chalcogenides, resulting in substantial improvements in their TE performance. The methodology used in this study can consider the influence of grain boundaries on the electronic and thermal transport properties of materials. The primary aim is to optimize ZT through nanostructuring by identifying a grain size window where κlcan be reduced without significantly reducing the PF. For instance, contour lines aligned parallel to the x-axis suggest that reducing the grain size, L, does not have a significant impact on ZT. This is the case for n-type In2AgCuSe4(Figure 6). However, for p-type In2AgCuSe4, the contour lines show an improvement of the ZT values (above 2) with average grain sizes below 100 nm (Figure 6). This indicates that nanostructuring can effectively enhance the TE performance of ptype In2AgCuSe4by reducing the thermal conductivity without compromising the PF. 3. CONCLUSIONS The application of design principles in materials discovery has become increasingly powerful with the emergence of HT techniques and ML. These methods open doors to on-the-fly design principles while exploring large chemical and physical spaces. This is particularly interesting in fields where multiple and competing objectives are targeted, such as thermoelectricity. In this work, DFT calculations, ML techniques, and the use of the Boltzmann transport equation for phonons were combined into an HT framework to explore the TE performance of DL quaternary chalcogenides. This framework has been used to investigate the lattice thermal conductivity of these compounds and simultaneously establish design principles to accelerate the discovery of new TE materials. First, we discussed the critical role of substituting Cu for Ag to reduce the lattice thermal conductivity. Next, we explored the differences in thermal transport between various polymorphs and identified the structural prototype kesterite as the best candidate. By combining these design principles with other factors such as stability and cost-effectiveness, In2AgCuSe4was proposed as a potential TE candidate. This material exhibits low lattice thermal conductivity (around 1 W/m K at 700 K) due to the scattering between acoustic and low-energy optical modes. The presence of different pockets with similar energies close to the top edge of the valence band makes In2AgCuSe4 an interesting p-type semiconductor with a large PF when the carrier concentration is around 1020 cm−3. While n-type In2AgCuSe4presents a ZT maximum around 0.7, p-type In2AgCuSe4can reach ZT values higher than 1.8 at 900 K. ZT can also be optimized by including the average grain size as a variable. For grain sizes around 100 nm, lattice thermal conductivity is drastically reduced, and PF is only slightly modified, resulting in ZT values larger than 2 for p-type In2AgCuSe4. 4. METHODOLOGY AND COMPUTATIONAL DETAILS 4.1. Thermal Transport Properties. 4.1.1. Geometry Optimization. All ground state structures at 0 K were fully optimized (atoms and lattice parameters) using the VASP package, 32,33 with projector-augmented wave potentials. 34 The PBE exchange−correlation functional, 35 along with Grimme D3 van der Waals correction, 36 were combined to obtain the energies. Core and valence electrons were selected following standards proposed by Calderon et al. 37 A high-energy cutoff of 500 eV and a dense k-point mesh of 2300 k-points per reciprocal atom were used. Convergence of the wave function was determined when the energy difference between consecutive electronic steps fell below 10−9eV. The geometry and lattice vectors were relaxed until forces on all atoms reached values below 10−7eV/Å, using a conventional cell containing 16 atoms. To minimize the noise in the forces, an additional support grid was included for the evaluation of the augmentation charges. 4.1.2. Supercell Single-point Calculations and Force Constants. Interatomic force constants (IFCs) were calculated using the hiPhive package, which combines the forces calculated for distorted atoms in supercells with ML regression. 38 A 3 ×3×2 supercell containing 288 atoms was used for the calculation of forces, following the same setup as the geometry optimizations. To ensure the accurate calculation of the IFCs, a two-step approach was implemented. 20 First, small random distortions were applied to all atoms in three supercells, and then secondand third-order IFCs were extracted using the hiPhive package. Subsequently, an additional set of 18 distorted supercells was generated by superimposing normal modes with random phase factors and amplitudes corresponding to a temperature of 300 K by employing the second-order IFCs obtained from the previous step. To maximize the use of data and minimize training bias, we incorporated 5-fold cross-validation into our workflow. This involved training the model on 80% of the data while rotating through different subsets for validation at each fold. The IFCs were determined using the recursive feature elimination (RFE) algorithm through multilinear regression to the DFT forces. Despite being more computationally expensive compared to ordinary least-squares regression, RFE has shown its effectiveness in achieving convergence with fewer structures. 8 Moreover, employing an RFE reduces the number of parameters and simplifies the model by retaining only the most significant interaction terms. IFCs were calculated including cutoffs for second-, third-, and fourth-order terms. To ensure transferability across different compounds, these cutoff values were determined based on coordination shells. 4.1.3. BTE Solver. The lattice thermal conductivity, κl, is calculated using the ShengBTE code, through the iterative solution of the BTE for phonons. 39 This method provides more accurate results compared with the relaxation time approximation. Scattering rates are computed by considering isotopic and three-phonon scattering processes. To balance memory demand and ensure convergence of κlwith respect to q-points, a Gaussian smearing of 0.1 is applied along with a dense mesh of 16 ×16 ×8q-points. 4.2. Electronic Transport Properties. The AMSET package, 9 consistent with the methodology outlined in our previous study, 17 was employed to calculate the electrical conductivity, the Seebeck coefficient, and the electronic contribution to thermal conductivity. This tool solves the Boltzmann transport equation to predict electronic transport properties using Onsager coefficients, with the wave function derived from a DFT calculation serving as the primary input. Scattering rates for each temperature, doping concentrations, bands, and k-points were determined by considering scattering ACS Applied Materials & Interfaces www.acsami.org Research Article https://doi.org/10.1021/acsami.4c04120 ACS Appl. Mater. Interfaces 2024, 16, 28590−28598 28596 due to deformation potentials, POPs, and ionized impurities. Wave function coefficients were calculated using the HSE06 functional 40,41 employing the primitive cell (8 atoms) and a 12 ×12 ×12 k-point mesh. Elastic constants and deformation potential, crucial to calculating different scattering contributions, were obtained following the same configuration used for geometry optimization and force constant calculations. ■ASSOCIATED CONTENT * sı Supporting Information The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acsami.4c04120. Lattice parameters comparison, validation RMSE in the forces, phonon density of states, and electron mobility. Further data sets are available at the ZENODO repository (doi: 10.5281/zenodo.10687721) (PDF) ■AUTHOR INFORMATION Corresponding Authors Pinku Nath −Institute for Chemical Reaction Design and Discovery (WPI-ICReDD), Hokkaido University, Sapporo 060-0808, Japan; Email: [email protected] 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; Email: [email protected] Authors Adolfo E. Rosado-Miranda −Departamento de Química Física, Facultad de Química, Universidad de Sevilla, Seville 41012, Spain Victor Posligua −Departamento de Química Física, Facultad de Química, Universidad de Sevilla, Seville 41012, Spain; orcid.org/0000-0003-3375-3706 Javier Fdez. Sanz −Departamento de Química Física, Facultad de Química, Universidad de Sevilla, Seville 41012, Spain; orcid.org/0000-0003-2064-7007 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 Complete contact information is available at: https://pubs.acs.org/10.1021/acsami.4c04120 Author Contributions § A.E.R.M. and V.P. contributed equally to this paper. Notes The authors declare no competing financial interest. ■ACKNOWLEDGMENTS This work was funded by Spanish MICIU/AEI/10.13039/ 501100011033 and by “European Union Next Generation EU/ PRTR” (grants PID2019-106871GB-I00, TED2021-130874BI00 and PID2022-138063OB-I00). We thankfully acknowledge the computer resources at Lusitania and the technical support provided by Cénits-COMPUTAEX and Red Espanola de Supercomputación, RES (QHS-2022-3-0032, QHS-2022-10031). This work was partially supported by the Institute for Chemical Reaction Design and Discovery (ICReDD), which was established by the World Premier International Research Initiative (WPI), MEXT, Japan. ■REFERENCES (1) Cheetham, A. K.; Seshadri, R.; Wudl, F. Chemical synthesis and materials discovery. Nat. Syn. 2022,1, 514−520. (2) Tabor, D. P.; Roch, L. M.; Saikin, S. K.; Kreisbeck, C.; Sheberla, D.; Montoya, J. H.; Dwaraknath, S.; Aykol, M.; Ortiz, C.; Tribukait, H.; Amador-Bedolla, C.; Brabec, C. J.; Maruyama, B.; Persson, K. A.; Aspuru-Guzik, A. Accelerating the discovery of materials for clean energy in the era of smart automation. Nat. Rev. Mater. 2018,3, 5−20. (3) Jablonka, K. M.; Jothiappan, G. M.; Wang, S.; Smit, B.; Yoo, B. Bias free multiobjective active learning for materials design and discovery. Nat. Commun. 2021,12, 2312. (4) Butler, K. T.; Frost, J. M.; Skelton, J. M.; Svane, K. L.; Walsh, A. Computational materials design of crystalline solids. Chem. Soc. Rev. 2016,45, 6138−6146. (5) Blancas, E. J.; Plata, J. J.; Santana, J.; Lemus-Prieto, F.; Márquez, A. M.; Sanz, J. F. Unraveling the role of chemical composition in the lattice thermal conductivity of oxychalcogenides as thermoelectric materials. J. Mater. Chem. A 2022,10, 19941−19952. (6) Santana-Andreo, J.; Márquez, A. M.; Plata, J. J.; Blancas, E.; González-Sánchez, J. L.; Sanz, J. F.; Nath, P. High-Throughput Prediction of the Thermal and Electronic Transport Properties of Large Physical and Chemical Spaces Accelerated by Machine Learning: Charting the ZT of Binary Skutterudites. ACS Appl. Mater. Interfaces 2024,16, 4606−4617. (7) Zhou, F.; Nielson, W.; Xia, Y.; Ozolins, V. Lattice Anharmonicity and Thermal Conductivity From Compressive Sensing of FirstPrinciples Calculations. Phys. Rev. Lett. 2014,113, 185501. (8) Fransson, E.; Eriksson, F.; Erhart, P. Efficient construction of linear models in materials modeling and applications to force constant expansions. npj Comput. Mater. 2020,6, 135. (9) Ganose, A.; Park, J.; Faghaninia, A.; Woods-Robinson, R.; Persson, K.; Jain, A. Efficient calculation of carrier scattering rates from first principles. Nat. Commun. 2021,12, 2222. (10) Graziosi, P.; Li, Z.; Neophytou, N. ElecTra code: Full-band electronic transport properties of materials. Comput. Phys. Commun. 2023,287, 108670. (11) Poncé, S.; Margine, E.; Verdi, C.; Giustino, F. EPW: Electron− phonon coupling, transport and superconducting properties using maximally localized Wannier functions. Comput. Phys. Commun. 2016, 209, 116−133. (12) Pal, K.; Xia, Y.; Shen, J.; He, J.; Luo, Y.; Kanatzidis, M. G.; Wolverton, C. Accelerated discovery of a large family of quaternary chalcogenides with very low lattice thermal conductivity. npj Comput. Mater. 2021,7, 82. (13) Pal, K.; Park, C. W.; Xia, Y.; Shen, J.; Wolverton, C. Scaleinvariant machine-learning model accelerates the discovery of quaternary chalcogenides with ultralow lattice thermal conductivity. npj Comput. Mater. 2022,8, 48. (14) Xi, L.; Pan, S.; Li, X.; Xu, Y.; Ni, J.; Sun, X.; Yang, J.; Luo, J.; Xi, J.; Zhu, W.; Li, X.; Jiang, D.; Dronskowski, R.; Shi, X.; Snyder, G. J.; Zhang, W. Discovery of High-Performance Thermoelectric Chalcogenides through Reliable High-Throughput Material Screening. J. Am. Chem. Soc. 2018,140, 10785−10793. (15) Ge, B.; Lee, H.; Zhou, C.; Lu, W.; Hu, J.; Yang, J.; Cho, S.-P.; Qiao, G.; Shi, Z.; Chung, I. Exceptionally low thermal conductivity realized in the chalcopyrite CuFeS2via atomic-level lattice engineering. Nano Energy 2022,94, 106941. (16) 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. Adv. Sci. 2018,5, 1700727. (17) Plata, J.; Blancas, E.; Márquez, A. M.; Posligua, V.; Fdez Sanz, J.; Grau-Crespo, R. Harnessing the unusually strong improvement of thermoelectric performance of AgInTe2with nanostructuring. J. Mater. Chem. A 2023,11, 16734−16742. (18) Ortiz, B. R.; Peng, W.; Gomes, L. C.; Gorai, P.; Zhu, T.; Smiadak, D. M.; Snyder, G. J.; Stevanovic, V.; Ertekin, E.; Zevalkink, A.; Toberer, E. S. Ultralow Thermal Conductivity in Diamond-Like ACS Applied Materials & Interfaces www.acsami.org Research Article https://doi.org/10.1021/acsami.4c04120 ACS Appl. Mater. Interfaces 2024, 16, 28590−28598 28597 Semiconductors: Selective Scattering of Phonons from Antisite Defects. Chem. Mater. 2018,30, 3395−3409. (19) 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. (20) Plata, J. J.; Posligua, V.; Márquez, A. M.; Fernandez Sanz, J.; Grau-Crespo, R. Charting the Lattice Thermal Conductivities of I− III−VI2Chalcopyrite Semiconductors. Chem. Mater. 2022,34, 2833− 2841. (21) Posligua, V.; Plata, J. J.; Márquez, A. M.; Sanz, J. F.; GrauCrespo, R. Theoretical investigation of the lattice thermal conductivities of II-IV-V2pnictide semiconductors. ACS Appl. Electron. Mater. 2023. (22) Dong, Y.; Khabibullin, A. R.; Wei, K.; Ge, Z.-H.; Martin, J.; Salvador, J. R.; Woods, L. M.; Nolas, G. S. Synthesis, transport properties, and electronic structure of Cu2CdSnTe4.Appl. Phys. Lett. 2014,104, 252107. (23) Souna, A. J.; Wei, K.; Nolas, G. S. The effect on the optical modes of quaternary chalcogenides upon metal and chalcogen substitution. Appl. Phys. Lett. 2020,116, 082103. (24) Deng, T.; Xing, T.; Brod, M. K.; Sheng, Y.; Qiu, P.; Veremchuk, I.; Song, Q.; Wei, T.-R.; Yang, J.; Snyder, G. J.; Grin, Y.; Chen, L.; Shi, X. Discovery of high-performance thermoelectric copper chalcogenide using modified diffusion-couple high-throughput synthesis and automated histogram analysis technique. Energy Environ. Sci. 2020, 13, 3041−3053. (25) Coughlan, C.; Ibanez, M.; Dobrozhan, O.; Singh, A.; Cabot, A.; Ryan, K. M. Compound copper chalcogenide nanocrystals. Chem. Rev. 2017,117, 5865−6109. (26) Xie, H.; Li, Z.; Liu, Y.; Zhang, Y.; Uher, C.; Dravid, V. P.; Wolverton, C.; Kanatzidis, M. G. Silver Atom Off-Centering in Diamondoid Solid Solutions Causes Crystallographic Distortion and Suppresses Lattice Thermal Conductivity. J. Am. Chem. Soc. 2023, 145, 3211−3220. (27) Xie, H.; Bozin, E. S.; Li, Z.; Abeykoon, M.; Banerjee, S.; Male, J. P.; Snyder, G. J.; Wolverton, C.; Billinge, S. J.; Kanatzidis, M. G. Hidden Local Symmetry Breaking in Silver Diamondoid Compounds is Root Cause of Ultralow Thermal Conductivity. Adv. Mater. 2022, 34, 2202255. (28) Xie, H.; Liu, Y.; Zhang, Y.; Hao, S.; Li, Z.; Cheng, M.; Cai, S.; Snyder, G. J.; Wolverton, C.; Uher, C.; Dravid, V. P.; Kanatzidis, M. G. High Thermoelectric Performance in Chalcopyrite Cu1−xAgxGaTe2-ZnTe: Nontrivial Band Structure and Dynamic Doping Effect. J. Am. Chem. Soc. 2022,144, 9113−9125. (29) Rincón, C. Order-disorder transition in ternary chalcopyrite compounds and pseudobinary alloys. Phys. Rev. B: Condens. Matter Mater. Phys. 1992,45, 12716−12719. (30) Gutay, L.; Regesch, D.; Larsen, J. K.; Aida, Y.; Depredurand, V.; Redinger, A.; Caneva, S.; Schorr, S.; Stephan, C.; Vidal, J.; Botti, S.; Siebentritt, S. Feedback mechanism for the stability of the band gap of CuInSe2.Phys. Rev. B: Condens. Matter Mater. Phys. 2012,86, 045216. (31) Miller, S. A.; Dylla, M.; Anand, S.; Gordiz, K.; Snyder, G. J.; Toberer, E. S. Empirical Modeling of Dopability in Diamond-Like Semiconductors. npj Comput. Mater. 2018,4, 71. (32) Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1993,47, 558−561. (33) Kresse, G.; Furthmuller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996,54, 11169−11186. (34) Blöchl, P. E. Projector Augmented-Wave Method. Phys. Rev. B: Condens. Matter Mater. Phys. 1994,50, 17953−17979. (35) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996,77, 3865−3868. (36) 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. (37) Calderon, C. E.; Plata, J. J.; Toher, C.; Oses, C.; Levy, O.; Fornari, M.; Natan, A.; Mehl, M. J.; Hart, G. L. W.; Buongiorno Nardelli, M.; Curtarolo, S. The AFLOW Standard for HighThroughput Materials Science Calculations. Comput. Mater. Sci. 2015,108 (Part A), 233−238. (38) Eriksson, F.; Fransson, E.; Erhart, P. The Hiphive Package for the Extraction of High-Order Force Constants by Machine Learning. Adv. Theory Simul. 2019,2, 1800184. (39) 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. (40) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid Functionals Based on a Screened Coulomb Potential. J. Chem. Phys. 2003,118, 8207−8215. (41) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Erratum: “Hybrid Functionals Based on a Screened Coulomb Potential” [J. Chem. Phys.118, 8207 (2003)]. J. Chem. Phys. 2006,124, 219906. ACS Applied Materials & Interfaces www.acsami.org Research Article https://doi.org/10.1021/acsami.4c04120 ACS Appl. Mater. Interfaces 2024, 16, 28590−28598 28598