Screening of Complex Layered Chalcogenide Structures as High-Performance Thermoelectrics by High-Throughput Calculations

Abstract

Thermoelectric materials have drawn much attention over the last two decades due to the increase in global energy demand. However, designing efficient thermoelectrics reveals itself as a tough task for their properties (Seebeck coefficient, electrical conductivity, thermal conductivity) are mutually opposed. Hence, most recently, new design approaches have appeared, among which high-throughput methods have been implemented either experimentally or computationally. In this work, a high-throughput computer program has been designed to generate over 4000 structures based on a small set of complex layered chalcogenide compounds taken from the mA^IVB^VI nA₂^VB₃^VI homologous series, where A^IV is Ge, A^V is Sb and B^VI is Te. The computer-generated structures have been investigated using density-functional theory methods, and the electronic and transport properties have been calculated. It has been found, using the quantum theory of atoms in molecules and crystals, that a wide variety of bond types constitutes the bonding network of the structures. All the structures are found to have negative formation energies. Among the obtained final structures, 43 are found with a wide band gap energy (>0.25 eV), 358 with semi-conductor/metal characteristics, and 731 with metallic characteristics. The transport properties calculations, using the Boltzmann equation, reveal that two p-type and 86 n-type structures are potentially promising compounds for thermoelectric applications.

1. Introduction

The limited efficiency of current thermoelectric devices and of their constitutive materials hinders their widespread commercial use. The efficiency of materials to convert heat into electricity is measured by the figure of merit ZT, where T is the operating temperature. Z is proportional to the square of the Seebeck coefficient S and to the electrical conductivity σ and to the reciprocal of the thermal conductivity κ (both electronic and phononic). For a ZT of 0.5, the efficiency reaches only 10%. Most of the state-of-the-art materials do not exceed ZT = 1.2 despite many efforts engaged in their improvement. By adjusting the carrier concentration [1,2,3,4,5], carrier mobility [5,6,7], and effective mass of the carriers through band convergence [8,9,10], it is possible to enhance the electrical conductivity or the Seebeck coefficient of the material, thereby optimizing its electrical transport performance. Alternatively, increasing ZT can be done by decreasing the lattice thermal conductivity by playing on, e.g., the atomic mass fluctuation [7,11,12], the phonon group velocity [13,14] and the phonon scattering mechanisms [15,16,17,18]. However, to date, most of the material improvements have been based on trial-and-error approaches.

As a recent method of finding new materials with enhanced properties, high-throughput computational material design has emerged as a new field in materials science, enabling the screening of large sets of compounds. It combines advanced thermodynamic and electronic-structure methods with intelligent data mining and supercomputer capabilities [19]. To avoid the significant cost of experimental synthesis and testing, high-throughput computational screening provides an effective way to identify high-efficiency thermoelectric candidates. This approach enables scientists to create, manage, and analyze extensive data repositories, revolutionizing materials research and unlocking new paths for material development. Several studies have contributed to identifying promising thermoelectric materials [20,21,22]. Yang et al. [23] employed the constant relaxation time approximation (CRTA) (with relaxation time τ of 10 fs) to explore the maximum power factors of 36 half-Heusler compounds, suggesting Co-/Rh-/Fe-based half-Heuslers and LaPdBi as promising candidates for p-type and n-type thermoelectric materials, respectively. Carrete et al. [24] combined a constant mean-free-path with the BoltzTraP code to predict the ZT values of 75 nanograined half-Heusler compounds, with approximately 15% showing potential to surpass ZT ≈ 2 at high temperatures. Gan et al. [25] investigated promising TE chalcogenides using a machine-learning-based approach with high-throughput ab initio calculations, achieving a peak ZT of 1.21 for n-type Pb₂Sb₂S₅. Jin et al. [26] conducted high-throughput calculations on 11,993 materials within the MatHub-3d database, identifying 9957 compounds with converged electrical transport properties. Among these, 156 compounds exhibited promising characteristics for both n- and p-type thermoelectric transport.

Complex layered chalcogenides are of great interest owing to the possibility of obtaining thermoelectric materials with low lattice heat conductivity [27], for example, the ternary layered compounds in the quasibinary A^IVB^VI–A₂^VB₃^VI systems (A^IV = Ge, Sn, Pb; A^V = Bi, Sb; B^VI = Te, Se). In these systems, homologous series of layered compounds of the mA^IVB^VI nA₂^VB₃^VI type with composite crystal lattices are formed. They are well known for their outstanding structural and electronic properties and include a wide variety of mixed-layer materials with more complex crystal structures than their parent A^IVB^VI and A₂^VB₃^VI compounds. Their structures are derived from tetradymite (Bi₂Te₂S) [28]; however, most of them are characterized by distorted rocksalt-type slabs of varying thickness. By using various combinations of layer stacking sequences, the system behaves like a compound-generating device. It offers large possibilities for optimizing the thermoelectric properties. Hence, in this work, a home-made python-based code has been developed to generate 4307 layered structures and, among these, 1132 structures have been successfully optimized to their ground state. Electronic properties, transport properties, and topological properties have been investigated. Two p-type and 86 n-type TE materials are identified as promising candidates for TE applications.

2. Computational Methods

The electronic structure calculations have been performed within the frame of the density functional theory using projector-augmented waves (PAW) [29] and planewaves techniques, as implemented in the VASP program [30,31,32]. The PBEsol [33] exchange-correlation functional has been used. As 4307 structures have been generated with a number of atoms in each unit cell varying from 5 to 63, full structural optimizations (atomic positions and cell parameters) have been performed along 3 steps from coarse level to accurate ones. For the coarse level, the kinetic energy cutoff has been set to the maximum value among the three atoms as recommended in the POTCAR file of VASP, i.e., 174.982 eV. For the more accurate levels, the kinetic energy cutoff was set to 300 eV and 400 eV. The energy and force thresholds employed for the structure optimizations have been set to 10⁻⁵ eV and 10⁻² eV/Å, and the energy threshold for the subsequent electronic properties calculations has been set to 10⁻⁶ eV. The Brillouin zone (BZ) has been sampled with a k-point grid determined from the KSPACING parameter of VASP that is used to calculate the number of k-points from the formula

$$N_i=\max \left(1, \mathrm{ceiling}\left(\frac{2\pi \left|b_i\right|}{\mathrm{KSPACING}}\right)\right),$$

(1)

where b_i is the lattice constant of the optimized structures. The values of KSPACING have been set to 0.8, 0.8, and 0.4 for each step. For the subsequent self-consistent (SCF) calculations, the k-spacing and planewave kinetic energy cutoff have been set to 0.1 and 450 eV. The electronic structure and charge density calculated from the SCF step have been used to analyze the topological properties within the quantum theory of atoms in molecules and crystals (QTAIMAC) [34], as implemented in Critic2 code [35]. The thermoelectric properties have been obtained from the semi-classical Boltzmann transport theory within the constant relaxation time approximation (CRTA) and the rigid band structure approximation as implemented in BoltzTraP2 code [36]. For the resolution of the Boltzmann transport equation, full BZ band structures are rebuilt for the k-points sampling and eigenvalues, for which a very dense k-point mesh is required, and the k-spacing has been set to 0.05.

3. Structure Generation

The workflow of the high-throughput calculations performed in this work is illustrated in Figure 1. First, the structures are generated by the python script and fully optimized with VASP. The electronic properties are subsequently calculated, and the band gap energies are extracted. If the gap is found to be opened, the transport and topological properties are calculated; otherwise, the structure is discarded. At this point, one should underline that GGA functionals underestimate band gap energies; hence, compounds with an opened bandgap may be erroneously discarded from the list of candidates. This is, however, a weaker drawback than if the process had led to keeping metals in the list.

The slab candidates used to build the layered structures stem from naturally existing compounds from the n(GeTe)-m(Sb₂Te₃) system: Sb₂Te₃, GeSb₂Te₄, and Ge₂Sb₂Te₅ with two stacking sequences. The atomic slabs are illustrated in Figure 2b–e. The bulk Sb has also been taken into account as it also possesses a layered structure (Figure 2a). The bulk Sb, Sb₂Te₃, and GeSb₂Te₄ crystallize in a rhombohedral lattice system (space group R$\overline{3}$m) with 2, 5, and 7 atoms in the primitive cell stacked along the c-axis, and they can also be modeled in hexagonal unit cells. Inside the bulk Sb, Sb₂Te₃, and GeSb₂Te₄, three different atomic slabs with 2, 5, and 7 atomic layers are found, respectively. The bulk Ge₂Sb₂Te₅ crystallizes in the P3m1 space group with two possible stacking sequences (stacking S1 and stacking S2, see Figure 2d,e), and consists of a nine-atomic-layer slab. The two sequences differ by the complete interchange of Sb and Te atoms in the layers. To provide a better illustration, a 1 × 1 × 2 supercell is depicted in Figure 2d,e. Hence, five atomic slabs consisting of 2, 5, 7, 9, and 9 atomic layers have been chosen as candidates to build the complex layered structures, the slabs being labeled as slab1, slab2, slab3, slab4, and slab5 in Figure 2. Each atom inside the slab is bonded with six adjacent atoms (Figure 2f), while the outermost Te atom in each slab is bonded with two atoms and presents a van der Waals interaction with the adjacent slab.

From the five slab candidates, many-layered structures can be built by combinations and perturbations. However, there are also many choices for the slab stacking sequences in the ab plane, which means that the coordinates of the outermost Te atoms from adjacent slabs should be constrained. To constrain the structures and maintain the space group of the complex structures, the coordinates in the ab plane have been constrained according to the following sequence: →(1/3, 2/3)→(2/3, 1/3)→(0,0)→. Three examples are shown in Figure 3 and the sequence of slabs with slab4—slab4—slab5—slab5, slab5—slab4—slab4—slab3—slab2, and slab3—slab1—slab3—slab4—slab5—slab2 in the unit cells are labeled as 4455, 54432, and 313452, respectively. Even though the atomic layers in a slab can also exist in many different possibilities, in this work, the slabs have been treated as a whole and the combinations have only been applied for the slab stacking sequences. The number of generated structures increases sharply with the increase in the number of slabs in the unit cells. Therefore, the maximum number of slabs has been set to seven. An amount of 4307 complex layered structures have been generated for the optimization step.

All the structures have been generated simply from geometric stacking, which cannot ensure the rationality of the structures. After generating the structures, we have optimized the so-generated structures from coarse to accurate levels by improving the kinetic cutoff and the k-mesh sampling of BZ. In the end, 1132 structures have been optimized with success.

4. Results

This section deals with the electron density topological analysis of the structures, the stability of the calculated structures, and their electronic and transport properties.

4.1. Electron Density Topology

Electron charge density values and distribution play an important role in the bonding trend characteristics and affect the thermoelectric properties. Therefore, the analysis of electron charge density and chemical bonding as developed by Bader and subsequent QTAIMAC methods can help us to explore the electronic character of the selected compounds. The key fields are the electron density ρ and its Laplacian ∇²ρ, from which the total energy density H, the kinetic energy density G, and potential one V can be derived. Besides, closed shell interactions (ionic, H-bonds, and van der Waals) have a large positive value of ∇²ρ over the entire interaction region, G/ρ > 1, |V|/G < 1, and a small ρ. Conversely, ∇²ρ < 0, G/ρ < 1, |V|/G > 2, and a large ρ are expected for shared interactions (covalent or polar bonds). Within the generated compounds, four types of bonds are found: Ge-Te, Sb-Sb, Sb-Te, and Te-Te. The corresponding QTAIMAC parameters are plotted in Figure 4. In agreement with the structure of the compounds, the QTAIMAC results evidence one type of bond for Ge-Te, Te-Te, and Sb-Sb and two types of bonds for Sb-Te: one corresponds to the bonds inside the slab and the other corresponds to a bond between slab1 and any other one. Small values of ρ and positive ∇²ρ with 1 < |V|/G < 2 have been evidenced at the BCPs of the Ge-Te and Sb-Te bonds. These characteristics agree with those defining the transit region according to Espinosa’s work [37]. In this region, bonds are neither pure covalent nor pure closed-shell interactions. The Te-Te bonds exhibit a positive bond degree (H/ρ) with |V|/G < 1 at their bond critical point. According to Espinoza’s bonding classification, these weak Te-Te bonds are associated with closed shell interactions and are expected to be of a van der Waals type. By contrast, most of the Sb-Sb bonds (those located in the slab) show a negative ∇²ρ, indicating an electron concentration at the bond critical point. The |V|/G ratio of the Sb-Sb bonds is also larger than two, suggesting that the Sb-Sb bonds involve covalent interactions.

4.2. Structure Stability

Stability is an essential property for any structure. In DFT, the current standard way to verify structure stability is to perform phonon calculations by means of DFPT. However, these calculations would require formidable computer resources for our many structures. Therefore, to investigate the relative stability of the layered structures, we start from the cohesive energy. Cohesive energy is the energy gained by arranging the atoms in a crystalline state, as compared with the gas state, and is given by the following:

$$E_{coh}^{{\mathrm{Ge}}_x{\mathrm{Sb}}_y{\mathrm{Te}}_z}=\frac{E_{\mathrm{tot}}^{{\mathrm{Ge}}_x{\mathrm{Sb}}_y{\mathrm{Te}}_z}-x{E}_{\mathrm{atom}}^{\mathrm{Ge}}-y{E}_{\mathrm{atom}}^{\mathrm{Sb}}-z{E}_{\mathrm{atom}}^{\mathrm{Te}}}{x+y+z}$$

(2)

where $E_{\mathrm{tot}}^{{\mathrm{Ge}}_x{\mathrm{Sb}}_y{\mathrm{Te}}_z}$ is the total energy and $E_{\mathrm{tot}}$ is the energy of the isolated atoms. Cohesive energy can provide information on the relative stability of the compounds of interest. Indeed, insulators and semiconductors, which have large cohesive energies, are strongly bonded and have good stability [38]. The calculated cohesive energy is shown in Figure 5. The value of the cohesive energies is projected on the left and right planes. The composition of the generated structures is projected to the bottom plane and shown by the ternary diagram. As discussed above, the layered structures are based on only five slabs: Sb, Sb₂Te₃, GeSb₂Te₄, Ge₂Sb₂Te₅-S1, and Ge₂Sb₂Te₅-S2; hence, the compositional range of each element is 22.2–100% for Sb, 0–22.2% for Ge, and 0–60% for Te. The projection of composition only occupies part of the ternary diagram. As shown in Figure 5, the range of the cohesive energy varies from −2.585 eV/at. to −2.922 eV/at.

As shown in the projection on the right plane (Figure 5), the value of cohesive energy decreases with the increase in the Ge content. According to the previous electron density topological analysis, all the bonding interactions inside the slabs are located at the transit region with neither pure covalent nor pure ionic interactions. Comparing the pure elements Sb, Ge, and Te, Ge exhibits stronger interactions, which is consistent with the calculated cohesive energy of −4.16 eV/at., −2.87 eV/at., and −3.12 eV/at. for Ge, Te, and Sb, respectively. As for Sb and Te, they have a similar impact on the cohesive energy.

Enthalpies of formation correspond to the difference between the total energy of a target compound A_xB_y and that of its constitutive elements A, B, and so on, taken in their standard reference state. The formation energy of a complex structure of Ge_xSb_yTe_z is calculated according to Equation (3):

$$E_{\mathrm{f}}^{{\mathrm{Ge}}_x{\mathrm{Sb}}_y{\mathrm{Te}}_z}=\frac{E_{\mathrm{tot}}^{{\mathrm{Ge}}_x{\mathrm{Sb}}_y{\mathrm{Te}}_z}-x{\mu }^{\mathrm{Ge}}-y{\mu }^{\mathrm{Sb}}-z{\mu }^{\mathrm{Te}}}{x+y+z}$$

(3)

where μ is the chemical potential of the pure elements. For most cases, chemical potentials are equal to the DFT total energies of their ground states.

Stevanović et al. [39] pointed out that when all the compounds and elements of interest pertain to similar classes of materials (e.g., metals, semiconductors), the calculation of the studied systems can be performed within one of the standard approximations of DFT, namely, the LDA or GGA, which benefit from the cancellation of errors associated with the similarly imperfect description of bonding in A_xB_y and its constitutive solid elements A and B. In our case, Ge, Sb, and Te elements are semiconductors located in the p block of the periodic table. All the generated structures are oxygen atom-free, which suppresses the problem of the inadequate calculation of the ground state of this atom by DFT at 0K. The calculated formation energies are shown in Figure 6. All formation energies present negative values varying from −0.53 eV/at. to −0.27 eV/at. Combining the projections on the right and left planes, we can observe that the formation enthalpy increases with the increase in Sb content. Based on the configuration of the structural generation, all layered structures are built from different slabs maintained together by van der Waals forces. Two new Sb-Te bonds are created when inserting slab1 (Sb-Sb) between two slabs terminated with Te (slabs two to five), which strengthen the new structures.

Table 1. Lattice parameters (in Å) and formation and cohesive energies (in eV/atom) of selected structures.

Compounds	31	222	333	4455	54432	313452
a	4.34	4.32	4.30	4.28	4.29	4.30
c	17.40	31.94	42.84	75.27	81.47	81.94
Formation energy	−0.32	−0.27	−0.28	−0.31	−0.30	−0.30
Cohesive energy	−2.80	−2.70	−2.84	−2.91	−2.87	−2.84

shows the calculated formation and cohesive energies of six selected structures that have been chosen based on their narrow band gap and high power factor (see electronic and transport properties sections).

4.3. Electronic Properties

Based on the optimized structures, the electronic band structures have been calculated with the PBEsol functional. Although the spin-orbit coupling (SOC) may have a sizable effect on the band structure for compounds with heavy elements, the calculations have been performed without SOC. Indeed, for our complex structures, SOC would have consumed too many calculation resources. The results are presented in Figure 7. The negative values are a consequence of the valence and conduction band overlapping and reflect the metallic character of the compound. As can be seen, the band gap varies from −0.535 to 0.655 eV. The distribution of band gap energies can be divided into three parts: the wide-gap region (gap > 0.25 eV), the semi-conductor/metal region (0 eV < gap < 0.25 eV), and the metallic region (gap < 0 eV).

An amount of 43 layered compounds are located in the wide-gap region. They are reported in

Table 2. Band gap values and types of compounds in the wide band gap region.

Compounds	CBM-VBM (eV)	CBM-VBM Type	Compounds	CBM-VBM (eV)	CBM-VBM Type
4	0.6553	Indirect	32	0.4727	Direct
34334	0.6457	Indirect	34232	0.4548	Direct
43334	0.6424	Indirect	432234	0.4441	Direct
3433	0.6242	Indirect	23342	0.4372	Direct
2434	0.5899	Indirect	42242	0.4364	Indirect
44423	0.5899	Indirect	23432	0.4246	Indirect
3443443	0.5833	Indirect	232233	0.4243	Direct
342	0.5755	Indirect	422424	0.4235	Indirect
432	0.5749	Indirect	3223	0.4166	Direct
333432	0.5748	Direct	4222	0.3872	Direct
44243	0.5747	Indirect	233223	0.3855	Direct
424233	0.5707	Indirect	242232	0.3752	Direct
4423	0.5574	Indirect	224232	0.3732	Direct
333423	0.5566	Indirect	422232	0.3717	Direct
333	0.5526	Indirect	222432	0.3411	Direct
434232	0.5429	Indirect	22244	0.3406	Direct
32324	0.5359	Direct	22232	0.3199	Direct
24233	0.5294	Indirect	333222	0.3156	Direct
33233	0.5251	Direct	2323222	0.2964	Direct
433323	0.5083	Direct	222	0.2922	Direct
2442442	0.4940	Direct	342222	0.2882	Indirect
4332334	0.4841	Direct

with their band gap value and gap type. Among these compounds, none of them contain slab1 and slab5. The compounds generated from only one type of slab (4, 333, and 222) belong to this region. The gap is decreasing from 4 to 333 and to 222. The decreasing band gap from the highest value is obtained starting from compound four. Generally, the compounds with high band gap values are located at the low-Sb part. The slab5 can seriously decrease the band gap for the models. Following the band theory, the Seebeck coefficient can be given as the following:

$$S=-\frac{k_B^2}{e}\frac{1}{n{\mu }_n+p{\mu }_p}\left\{\left(2-\frac{E_f}{k_{B}T}\right)n{\mu }_n-\left(2-\frac{E_f+E_g}{k_{B}T}\right)p{\mu }_p\right\}$$

(4)

where k_B is the Boltzmann constant, e is the elementary charge, E_f and E_g are the Fermi level and the band gap energy, μ_n and μ_p are the mobility of electrons and holes, and n and p are the number of electrons and holes, respectively. There is a positive relationship between the Seebeck coefficient and band gap values. However, the electrical conductivity decreases sharply with the increase in the band gap and deteriorates the power factor (PF = S²σ).

An amount of 358 layered compounds belong to the semi-conductor/metal region, among which, the more they are metallic, the more slab1 and slab5 are represented. The metallic character of Ge₂Sb₂Te₅-S2, from which slab5 is extracted, has been assessed (see Ref. [10]). The compounds in the semi-conductor/metal region are expected to present promising TE properties with moderate Seebeck coefficients and electrical conductivity. The remaining 731 layered compounds have been evidenced to be metallic ones.

For the following, we have chosen six compounds located in the wide-gap and semi-conductor/metal region and calculated their electronic structures. The band structures computed with the PBEsol functional are shown in Figure 8 and the corresponding DOS is illustrated in Figure 9. The band structures range from −2 eV to 2 eV. The space group numbers of the generated compounds are 156, 164, or 166. To keep the k-path consistent, the conventional unit-cells (Figure 10b) are used to calculate the band structures for the compound with space group 166 (Figure 10a). The shapes of the Brillouin zone of their primitive and conventional cells are represented, as well as the path between high symmetry points chosen for the band structure calculations, which is Γ-K-M-Γ-A-H-L-A. Due to the excessive delocalization of the occupied states by the GGA functional, a serious underestimation of the band gap is usually observed, which is a common problem in GGA functionals. However, since there are too many compounds to investigate, we only have used the PBEsol functional to calculate the transport properties.

The layered compound 31 exhibits characteristics of an indirect band gap semiconductor with an energy band gap of 0.006 eV, with CBM located on the A-H line, and VBM located on the Γ-K line. There are four second-highest VBMs and four second-lowest CBMs within an energy range of 0.1 eV relative to VBMs and CBMs, respectively. However, the band near the Fermi level is very soft, and the slope of the DOS is moderate in both the valence and conduction parts.

Sb₂Te₃ (222) belongs to the rhombohedral system with a space group R$\overline{3}$m. When the structure is converted to the space group P3m1, the conventional cell contains three five-atom-layer slabs. After optimization, the lattice constants are a,b = 4.32 Å and c = 31.94 Å. As shown in Figure 8b, Sb₂Te₃ is a multi-valley semiconductor with a direct band gap of 0.292 eV. Both the CBM and the VBM are situated at the Γ point (0, 0, 0). Within the conduction bands, three second-highest minima can be found, two of them along the Γ-M line located at (0.058, 0.058, 0) and (0.153, 0.153, 0) and one along the L-A line located at (0.104, 0.104, 0.500). Multi-valley band structures can boost the TE properties: a multi-gap between the valley providing several electron transition pathways, and a band alignment boosting the Seebeck coefficient. From the density of states shown in Figure 9, we observe that the slope of the total DOS near the VBM is much more cliffy than that of the CBM, indicating a higher Seebeck coefficient for p-type doping. The DOS of 222 has a strong Te-p state character near the Fermi level, while the conduction band orbitals near the CBM are equally contributed by the hybridization between Te-p and Sb-p states.

In the case of 333 (Figure 8c and Figure 9c), the band structure evidences the features of a semiconductor with a direct band gap. The lowest energy value of the conduction band and the highest one of the valence band are both located at point A, and the band gap energy value is 0.553 eV. Moreover, at the bottom of the conduction band, there are multiple valleys with similar energy minima, resulting in high degeneracy. Therefore, the Seebeck coefficient of the n-type doped compound should theoretically be better than those of the p-type doped one. The valence and conduction bands are mainly contributed by Te-p and Te-p+Sb-p+Ge-p orbitals, respectively.

The band gap of 4455 (Figure 8d and Figure 9d) is an indirect one, with the lowest energy value of the conduction band at point A and the highest energy value of the valence band between L and A. The energy gap value is 0.011 eV. Near Γ, the valence band is very flat and has small parabolic characteristics (Te-p type orbitals); hence, holes should have a large effective mass resulting in a high Seebeck coefficient for p-type doped 4455.

For 54432 (Figure 8e and Figure 9e), the CBM is located at the Γ point, the VBM is located between L and A, and the energy band gap is 0.021 eV. Notably, near the top of the valence band, several states bear a similar energy, resulting in converging energy bands. So, the p-type doped compound may have a higher Seebeck coefficient than the n-type doped one. However, near Γ, the conduction band is very flat; hence, electrons should have a large effective mass resulting in a larger Seebeck coefficient for the n-type doped 54432 than for the p-type one. As we can see, these two effects on the Seebeck coefficient contradict each other, so it is difficult to conclude about the thermopower on the basis of the description of the band structure only.

For 313452, the lowest energy values of the conduction band are all located between Γ and M, the highest energy values of the valence band are all located between A and H, and the band gap energy value is 0.003 eV. The valence band is contributed by a combination of the p orbitals of Te, Sb, and Ge.

It can be seen from the density of states diagrams that the DOS is zero in a small range near the Fermi level. In the valence band or conduction band, the density of states increases, and it has a steep peak near the Fermi level, especially in 4455 and 54432, and to a lesser extent in 313452. This behavior should result in a large Seebeck coefficient, which is beneficial to achieving a high power factor.

4.4. Transport Properties

The electronic conductivity (σ) and the electronic thermal conductivity (κ_e) calculated by BolzTrap2 are given as the ratios σ/τ and κe/τ, respectively. The relaxation time τ depends on the charge carrier concentration, the temperature, and the electron energy, which can be evaluated from deformation potential (DP) theory. However, due to the high computational burden of DP calculations, we keep τ undetermined in the following. The best TE properties are obtained for chalcogenides with carrier concentrations around 10¹⁹ carriers/cm³, so the doping dependence has been investigated in the domain of carrier concentrations ranging from 10¹⁷ cm⁻³ to 10²² cm⁻³. The calculations have been performed for both electron-doping and hole-doping. The temperature and doping dependence of transport properties S, σ/τ, and PF (S²σ) are shown in Figure 11, Figure 12 and Figure 13 for all the layered compounds having a GGA band gap higher than 0. eV (around 401 compounds). In these figures, the diameter and the color of the balls are representative of the magnitude of the power factor. The graphs show an uneven spreading of points along the x-axis. Along the y-axis, opposite signs of the Seebeck coefficient can be found on the two halves of the graph with respect to the zero due to the two types of doping. The color gradient shows a steady increase in PF toward maximization when either the Seebeck coefficient or the electrical conductivity increases.

By contrast to what is observed for intermediate and high doping levels (Figure 12 and Figure 13), a more pronounced asymmetrical distribution of the Seebeck coefficient values between the positive and negative ones is observed at low doping levels (Figure 11). This could be related to an insufficient doping level to compensate for the pristine carrier concentration. With the increase in temperature, the electrical conductivity increases, while the Seebeck coefficient decreases. It is difficult to reach both high Seebeck and high conductivity at the same time, as shown by the absence of points in that region. The highest Seebeck coefficients can be reached at moderate temperatures and low doping levels. The distribution of points according to their color suggests that the Seebeck coefficient has stronger effects on the PF than the electrical conductivity, as the points with the same color parallel the x-axis.

At low doping levels, irrespective of the temperature, (see Figure 11a), the n-type compounds present a higher power factor than the p-type ones. Nonetheless, with the increasing temperature, more p-type compounds show high PF.

At moderate- and high-carrier doping levels (Figure 12 and Figure 13), irrespective of the charge carrier types, the Seebeck coefficient values spread over a slightly larger range as the temperature increases, and the compounds bearing a low Seebeck coefficient value (around zero) tend to migrate toward the extremes of the domain. The largest number of compounds bearing high PF values has been found for the intermediate doping level.

The maximum PF/τ is found in the layered compound 3541153, which is 2.79 × 10¹² Wm⁻¹K⁻²s⁻¹ for an electron doping level of 3.16 × 10²¹ cm⁻³ and at 850 K. Two p-type and 86 n-type layered compounds have been identified as promising candidates (PF/τ > 10¹² Wm⁻¹K⁻²s⁻¹) for TE applications.

5. Conclusions

In this work, 4307 layered structures have been generated by a high-throughput computer program, among which 1132 structures have been successfully optimized to their ground states. Electronic properties, transport properties, and topological properties have been investigated. All Ge-Te and Sb-Te bonds are located at the transit region with small values of electron density ρ and positive Laplacian ∇²ρ with 1 < |V|/G < 2, while Te-Te BCPs exhibit a positive bond degree (H/ρ) with |V|/G < 1. The maximum PF/τ is found in the layered compound 3541153, which is 2.79 × 10¹² Wm⁻¹K⁻²s⁻¹ for an electron doping level of 3.16 × 10²¹ cm⁻³ and at 850 K. Two p-type and 86 n-type layered compounds have been identified as promising candidates for TE applications.